
==== Front
FEMS Microbiol Lett
FEMS Microbiol Lett
femsle
FEMS Microbiology Letters
0378-1097
1574-6968
Oxford University Press

39165128
10.1093/femsle/fnae068
fnae068
Research Article
Taxonomy, Systematics & Evolutionary Microbiology
AcademicSubjects/SCI01150
A universal and constant rate of gene content change traces pangenome flux to LUCA
Trost Katharina Faculty of Mathematics and Natural Sciences, Institute of Molecular Evolution, Heinrich Heine University Düsseldorf, 40225 Düsseldorf, Germany

Knopp Michael R Faculty of Mathematics and Natural Sciences, Institute of Molecular Evolution, Heinrich Heine University Düsseldorf, 40225 Düsseldorf, Germany

Wimmer Jessica L E Faculty of Mathematics and Natural Sciences, Institute of Molecular Evolution, Heinrich Heine University Düsseldorf, 40225 Düsseldorf, Germany

https://orcid.org/0000-0003-4291-2115
Tria Fernando D K Faculty of Mathematics and Natural Sciences, Institute of Molecular Evolution, Heinrich Heine University Düsseldorf, 40225 Düsseldorf, Germany
Los Alamos National Laboratory, Los Alamos, NM, United States

Martin William F Faculty of Mathematics and Natural Sciences, Institute of Molecular Evolution, Heinrich Heine University Düsseldorf, 40225 Düsseldorf, Germany

Corresponding author. Faculty of Mathematics and Natural Sciences, Institute of Molecular Evolution, Heinrich Heine University Düsseldorf, 40225 Düsseldorf, Germany. E-mail: katharina.trost@hhu.de
2024
20 8 2024
20 8 2024
371 fnae06822 3 2024
15 5 2024
19 8 2024
12 9 2024
© The Author(s) 2024. Published by Oxford University Press on behalf of FEMS.
2024
https://creativecommons.org/licenses/by-nc-nd/4.0/ This is an Open Access article distributed under the terms of the Creative Commons Attribution-NonCommercial-NoDerivs License (https://creativecommons.org/licenses/by-nc-nd/4.0/), which permits non-commercial reproduction and distribution of the work, in any medium, provided the original work is not altered or transformed in any way, and that the work is properly cited. For commercial re-use, please contact journals.permissions@oup.com

Abstract

Prokaryotic genomes constantly undergo gene flux via lateral gene transfer, generating a pangenome structure consisting of a conserved core genome surrounded by a more variable accessory genome shell. Over time, flux generates change in genome content. Here, we measure and compare the rate of genome flux for 5655 prokaryotic genomes as a function of amino acid sequence divergence in 36 universally distributed proteins of the informational core (IC). We find a clock of gene content change. The long-term average rate of gene content flux is remarkably constant across all higher prokaryotic taxa sampled, whereby the size of the accessory genome—the proportion of the genome harboring gene content difference for genome pairs—varies across taxa. The proportion of species-level accessory genes per genome, varies from 0% (Chlamydia) to 30%–33% (Alphaproteobacteria, Gammaproteobacteria, and Clostridia). A clock-like rate of gene content change across all prokaryotic taxa sampled suggest that pangenome structure is a general feature of prokaryotic genomes and that it has been in existence since the divergence of bacteria and archaea.

Analysis of 5655 prokaryotic genomes and 2872 metagenomic assemblies reveals a constant long-term gene flux rate across all prokaryotic groups that traces pangenome flux back to the last universal common ancestor (LUCA).

gene flux
pangenomes
prokaryotes
accessory genome
core genome
metagenomes
European Research Council 10.13039/501100000781 101018894
==== Body
pmcIntroduction

The evolution of genome diversification in eukaryotes is mostly driven by gene duplication and differential loss (Albalat and Cañestro 2016, Stull et al. 2021). In contrast, prokaryotic genome evolution is driven mostly by gene loss and gene acquisition via lateral (or horizontal) transfer (LGT), while gene duplication is rare (Mira et al. 2001, Treangen and Rocha 2011, Tria and Martin 2021). Genetic recombination and gene transfer in prokaryotes involves unidirectional transfer of genes from donors to recipients via transformation, transduction, conjugation, gene transfer agents, or membrane vesicles (Arnold et al. 2021). Over time, these mechanisms generate a process of DNA flux through constant gene loss and gain in prokaryotic chromosomes. How much gene flux they generate and whether such flux has been in operation throughout evolutionary history are questions of interest. In early work using codon bias as a proxy for laterally acquired genes, Lawrence and Ochman (1998) estimated that about 18% of the genes in Escherichia coli MG1655 genome correspond to recent lateral acquisitions. Today, it is recognized that the gene content of prokaryotic species is typically organized as pangenomes (Tettelin et al. 2005), with a conserved core genome consisting of genes present in all genomes of a given taxonomic sample, such as strains, and an accessory genome consisting of genes that are differentially present in all sample members that is, in evolutionary terms, in a state of continuous flux (Medini et al. 2005, Tettelin et al. 2005, 2008, Vernikos et al. 2015, Brockhurst et al. 2019).

The rate of gene flux between genomes of prokaryotes has been extensively studied at the species and genus level and several studies uncovered a clear relationship between phylogenetic distance and the frequency of gene differences that may arise by gain or loss (gain/loss). For example, Hao and Golding (2006) examined the association between the rate of gene-flux and point mutations in the core genome of seven Bacillus cereus strains; they estimated a rate of 4.4 gene gain or loss per nucleotide substitution per site in the strains’ core genome. Applications of the same method yielded similar inferences of 1.17 and 1.18 gene gain/loss events per point mutation per site for 12 Streptococcus genomes (Marri et al. 2006) and five Corynebacterium genomes (Marri et al. 2007), respectively. Higher rates of gene gain/loss were inferred for 27 Pseudomonas syringae strains, where a comparison of gene gain/loss among closely related strains showed that up to 5000 gain/loss events may have occurred before 1% amino acid sequence divergence in the core genomes was reached (Nowell et al. 2014). Other studies reported a positive association between sequence divergence of core genes and divergence in gene content (Wolf et al. 2016). For example, the more distant the genomes of two E. coli strains are, the fewer genes they share, although the relationship was weak (Touchon et al. 2009, Rocha 2018, Haudiquet et al. 2022). Similarly, a study of 22 Myxococcus xanthus strains showed that the number of gene differences is increased with phylogenetic distance as inferred from amino acid differences in the core genome (Wielgoss et al. 2016). Since the estimation of diverged gene content may be biased by differences in genome size, the “genome fluidity” metric was proposed as an unbiased estimate (Kislyuk et al. 2011), which is positively associated with the core genome sequence divergence at synonymous sites (Andreani et al. 2017).

While these studies arrived at similar conclusions regarding the rate of gene flux and phylogenetic relatedness, a direct comparison of the estimated rates across studies and taxa is challenging for several reasons. First, the size and the composition of the core genome varies strongly depending on the pangenome taxonomic composition (Vernikos et al. 2015), such that estimates of core genome divergence vary across studies. Second, the criteria for quantifying the number of shared genes between genomes (or differences in shared genes) can be very similar or even identical across studies, but the gene sets used to estimate sequence divergence are not. As a consequence, most estimates for rates of gene flux relative to sequence divergence are not comparable across species, samples, taxonomic levels, or studies.

Here, we ask whether prokaryotic genomes harbor evidence for a general correlation between the rate of gene flux into and out of genomes as a function of sequence divergence. For this purpose, we compared sequence divergence in a universal set of 36 proteins that are present in almost all genomes across higher taxa and are sufficiently conserved to be useful for comparisons at the deepest taxonomic levels (Hansmann and Martin 2000, Charlebois and Doolittle 2004, Dagan and Martin 2006). We then tested for an association between sequence divergence in these core genes and gene content divergence in the complete genomes of 5655 taxonomically diverse isolates as well as in 2872 metagenomic assemblies (MAGs).

Methods

Prokaryotic dataset

The prokaryotic clustering set was used from Brueckner and Martin (2020) including 5655 prokaryotic genomes, 19 050 992 protein sequences from the Reference Sequence database (RefSeq), September 2016 from the National Center for Biotechnology Information (NCBI; O’Leary et al. 2016). The clustering was created using the Markov Cluster Algorithm (MCL; van Dongen 2008) as previously described (Brueckner and Martin 2020, Nagies et al. 2020). In total, 450 283 protein families were detected and for protein families with at least four protein sequences multiple alignments were made with Mafft L-INS-I version 7.130 (Katoh 2002).

Due to large sample numbers, Proteobacteria and Firmicutes were divided into classes. Archaea were divided into orders to allow multiple groups to be assessed. The resulting 59 prokaryotic taxa comprise 41 bacterial and 18 archaeal groups that are called higher taxa in the following (Table S1).

The dataset including MAGs was obtained from Garg et al. (2021) including 103 assemblies from NCBI BioProject PRJN270657 and 2546 assemblies from BioProject PRJNA288027 downloaded in 2018. Additionally, 223 assemblies form the Microbial dark matter project were added (Rinke et al. 2013). The dataset was clustered using the following pipeline: to search for local alignments, an all versus all blastp was conducted using diamond version 2.0.11 (Buchfink et al. 2015). Reciprocal best blast hits (Wolf and Koonin 2012) with an e-value ≤1E-10 were aligned globally with the Needleman–Wunsch Algorithm (Emboss Needle version 6.6.0.0; Rice et al. 2000). All global alignments with an Identity ≥25% were used for clustering into protein families, using MCL (van Dongen 2008) version 14–137 with pruning parameters -P 180000, -S 19800, and -R 25200 (P = Pruning, S = Selection, and R = Recovery). In total, 285 787 protein families were detected. The completeness and contamination of all MAGs was measured by using CheckM v1.2.1 (Parks et al. 2015). Two out of the 2872 MAGs were not included in the clustering.

Protein family annotation via KEGG

For protein family annotation, all clustered sequences from the RefSeq dataset were blasted against the KEGG database using diamond 2.0.1 (Buchfink et al. 2015, Kanehisa et al. 2017). All best hits with at least 25% identity and a maximum e-value of 1E-10 were used for annotation. Based on these hits, a KO and a name were assigned to the clusters based on majority rule. Protein families, which contained equal to or more than 75% of unknown sequences, were not annotated to preserve the strict nature of the annotation method.

Verticality values

Verticality is a measure for how often members of a given protein family tend to recover monophyly of prokaryotic phyla in their respective gene trees (Nagies et al. 2020). Verticality values used in this analysis were obtained from the study by Nagies et al. (2020), based on the same dataset used in this study.

Determination of IC gene set

For the 260 972 prokaryotic protein families of the RefSeq dataset with a determined verticality value, a weight was calculated for the number of genomes, the number of higher taxa as well as the verticality. To obtain only one value per protein family, the average weight was calculated. To avoid taking genes that are particularly universal but not sufficiently vertical, or vice versa, the 50 protein families with the best average weights were compared to the 100 most vertical, the 100 most universal protein families based on number of genomes and to the 100 protein families that are most widely distributed across all higher taxa. Protein families present in all three lists and among the lowest 50 weights were assigned as informational core (IC) genes (Table S2).

The corresponding metagenomic IC clusters were obtained by using the reciprocal best cluster approach described in Ku et al. (2015): if 50% of all sequences of a prokaryotic cluster have their best hit in another cluster and if in this cluster also 50% of all sequences have their best hit in the prokaryotic cluster, it can be defined as a reciprocal best cluster. In the dataset, 36 IC genes were found. For every IC gene a multiple sequence alignment was calculated using Mafft L-INS-I version 7.505 (Katoh 2002).

Calculation of sequence divergence in the IC gene set (ICD)

For all 15 986 685 prokaryotic genome comparisons in the RefSeq dataset and for 4 100 275 genome pairs in the metagenomic dataset, the average sequence divergence in the IC gene set was calculated. In every genome comparison, for each IC gene (x), the proportion of different sites (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $IC{{D}_x}$\end{document}) in the multiple sequence alignment (a, lengtha = number of sites in the multiple sequence alignment) was calculated with the following formula:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*} IC{{D}_x} = \ \frac{{\textit{diff}\_site{{s}_a}}}{{\textit{lengt}{{h}_a}}}. \end{eqnarray*}\end{document}

The average proportion of different sites of all IC genes (n) is then determined as IC gene divergence ICD per genome comparison.

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*} ICD = \ \frac{{\mathop \sum \nolimits_{x = 1}^n IC{{D}_x}}}{n}. \end{eqnarray*}\end{document}

If a genome was not present in all 36 IC genes, the average divergence (ICD) was calculated for all IC genes in which the genome is represented by at least one protein.

Calculation of gene sharing distance (dsg)

For all 15 986 685 prokaryotic genome comparisons in the RefSeq dataset and for 4 100 756 genome pairs in the metagenomic dataset, the difference in gene content was calculated by comparing the protein families corresponding to the compared genomes (Fig. 1B). Two groups were defined, the paired protein families and the unpaired protein families. Paired protein families are protein families in which both genomes were represented by at least one protein. If only one genome is present in the protein family, the protein family corresponds to the group of unpaired protein families. The gene sharing distance dsg is defined as one minus the number of paired protein families divided by the average number of all protein families in the compared genomes.

Figure 1. Comparison of number of higher taxa, number of genomes, and verticality for the IC genes (A) and calculation of gene sharing distance dsg (B). (A) The number of higher taxa (x-axis) is plotted in relation to the number of genomes (z-axis) and verticality (y-axis) for the 36 defined IC genes (red) and for all other protein families from Nagies et al. (2020) (black). Higher taxa include all phyla present in the dataset except of Proteobacteria and Firmicutes, which were divided into classes due to large sample sizes and Archaea that were divided into orders to allow multiple groups to be assessed. Gene names for the 36 IC genes are shown on the right part of (A) and are sorted by verticality (B). The gene sharing distance (dsg) was calculated by assigning protein families from the compared genomes into two groups. Protein families present in both genomes are defined as paired protein families, the others as unpaired protein families. The gene sharing distance (dsg), which represents differences in gene content, was then calculated by subtracting the proportion of paired protein families from one. We expect that genomes belonging to the same species (intraspecies comparisons) have a higher proportion of paired protein families than genomes belonging to different higher-levels and thus show a lower gene sharing distance (PG = protein families of genome).

Filtering for genome size and strains

To avoid overplotting, the original prokaryotic RefSeq dataset was filtered for genomes with genome sizes within average higher taxa genome size plus/minus one standard deviation. Additionally, to reduce phylogenetic bias, only one random strain per species was retained in the dataset. The resulting dataset comprises 1630 strains (Table S1).

Calculation of average nucleotide identity values

The average nucleotide identities (ANI) were calculated using FastANI v1.34 (Jain et al. 2018).

Calculation of rRNA sequence divergence

For every species represented by at least 10 strains in the prokaryotic RefSeq dataset, rRNA sequences were downloaded from the RefSeq database of NCBI. A multiple alignment for all sequences was made with Mafft L-INS-I version 7.505 (Katoh 2002). Equal sites were then counted and the proportion of equal sites from all sites was subtracted from one to obtain the rRNA sequence divergence per genome pair. Values for dsg calculated from the prokaryotic dataset were assigned by using a random strain per species.

Statistical tests

To analyse the relationship between the ICD and dsg across prokaryotic taxa, Spearman rank correlations and linear regression lines were calculated. The regression lines were evaluated doing an analysis of residuals, including residual plots, mean square residuals, QQ plots of residuals, and the Kolmogorov–Smirnov test to test whether residuals are normally distributed. The Spearman rank correlation and the linear regression were also applied for the relationship between species-level accessory genome proportion and y-axis intercept as well as for the relationship between ICD and verticality of paired and unpaired protein families. To compare the distributions of paired and unpaired protein families, a paired t-test was applied. The cloud gene analysis was evaluated using a one-way ANOVA. All statistical tests were performed using python.

Calculation of average species level accessory genome proportion per higher taxon

The average species level accessory genome proportion per higher taxon are calculated based on intraspecies genome comparisons. For every species, the average gene sharing distance (dsg) of intraspecies genome comparisons was calculated. The mean of all species level accessory genome proportions corresponding to a specific higher taxa represents the average species level accessory proportion per higher taxon. The analysis was made for all higher taxa present in the filtered dataset, that are represented by at least 10 strains, 2 species, and with a P-value lower or equal to .05 from the Spearman correlations between ICD and dsg.

Calculation of absolute gene flux rates

Absolute gene flux rates were calculated by use of the relative gene flux rates, which show the percentage of gene differences in the genome per one % ICD. The absolute number of substitutions in IC genes then correspond to one % of the average number of sites present in the IC gene alignments. Gene differences at one % ICD were then calculated by multiplying the relative gene flux rate with the average number of genes present in the higher taxa. Then the number of gene differences at one substitution in the IC genes could be calculated by dividing the number of gene differences at one % ICD by the number of substitutions at one % ICD. The analysis was made for all higher taxa present in the filtered dataset, that are represented by at least 10 strains and with a P-value lower or equal to .05 from the Spearman correlations between ICD and dsg.

Eukaryotic organelle dataset

From the RefSeq Database (NCBI) 207 plastid and 68 mitochondrial proteins were downloaded in October 2022 for the species Porphyra umbilicalis and Reclinomonas americana as well as 189 089 nuclear protein sequences for the higher eukaryotic group Discoba. Additionally, for 150 eukaryotic genomes, 3 420 731 protein sequences were downloaded from the RefSeq Database, January 2018 (NCBI; O’Leary et al. 2016). The nuclear genome of P. umbilicalis was downloaded from the universal protein knowledgebase (UniProt Consortium 2021) in October 2022, including 13 559 protein sequences. To obtain mitochondrial homologous for prokaryotic IC genes, alphaproteobacterial IC gene sequences were searched in mitochondrial protein sequences of model organism R. americana by using diamond 2.0.1 (Buchfink et al. 2015). After filtering for identities higher or equal to 25%, a maximum e-value of 1E-10 and a search for best hits, 11 mitochondrial homologous for IC genes were defined based on majority rule. The remaining genes were searched in the Discoba proteome (24 IC genes) and the eukaryotic genomes (25 IC genes). The same search was made for cyanobacterial IC genes in plastid proteins of P. umbilicalis (14 IC genes) and in nuclear P. umbilicalis proteins (19 IC genes) as well as in eukaryotic genomes (22 IC genes). This resulted in four groups that represent the IC genes in mitochondria or plastids combined from organelle and nuclear protein sequences.

Informational core gene divergence between organelle genes and corresponding alphaproteobacterial or cyanobacterial genes

For the four datasets of mitochondrial or plastid universal genes, multiple alignments were made for every IC gene, using Mafft L-INS-I v7.505 (Katoh 2002) including the organelle (mitochondrial or plastid) protein sequence and the alphaproteobacterial or cyanobacterial protein sequences corresponding to the IC gene family. Average IC gene divergence (ICD) were then calculated between the organelle sequences and the corresponding prokaryotic sequences as described above (see the section “Calculation of sequence divergence in the IC gene set (ICD)”).

Results and discussion

The IC is vertically inherited

For the estimation of sequence divergence, we selected genes that have a nearly universal distribution and a low degree of LGT , that is, a high level of verticality (Nagies et al. 2020). These proportions are plotted for all protein families including at least four genomes from two or more higher taxa, shown in Fig. 1(A). Higher taxa include all phyla present in the dataset except for Proteobacteria and Firmicutes, which were divided into classes due to large sample sizes and Archaea that were divided into orders to allow multiple groups to be assessed. The 36 genes (red) that are widely distributed across higher prokaryotic taxa and genomes and exhibit high levels of verticality are mostly involved in information processing (Rivera et al. 1998), hence we call this set the informational core (IC) (Table S2). On average, IC genes display a verticality of 17.09 and are present in over 5585 prokaryotic genomes and 58 higher taxa. Verticality is a measure for how often members of a given protein family tend to recover monophyly of prokaryotic phyla in their respective gene trees (Nagies et al. 2020). In the data set of Nagies et al. (2020), verticality values can vary between 0 and 42, however, the highest calculated verticality value was 24.0 for 30S ribosomal protein S10 (RP-S10). By contrast, most of the other non-IC genes (black) are distributed across few genomes (mean number of genomes = 66, mean number of higher taxa = 2) and are inherited much less vertically (mean = 0.15). Functional annotations of the IC gene set using KEGG (Kanehisa et al. 2017) showed that the IC, selected here on the basis of universal distribution and verticality, comprises mainly genes encoding for ribosomal proteins, in addition to amino acid tRNA synthetases, translation elongation factors, RNA modifications enzymes, and several genes involved in metabolism (enolase, 3-phophoglycerate kinase, and pyrimidine synthesis enzymes). From the multiple amino acid sequence alignments of the 36 IC genes, we used the proportion of amino acid differences in pairwise comparisons as a robust measure for evolutionary divergence between genome pairs, termed here informational core gene divergence (ICD). The functional distribution of IC genes as well as their tendency to be vertically inherited is in line with previous studies about the properties of universally conserved (core) genes (Tettelin et al. 2005).

To quantify gene content differences for each genome pair, we determined the proportion of paired protein families (Fig. 1B, light red), which is the number of protein families in our clustered sample that is present in both of the compared genomes divided by the average number of protein families present in the two genomes. The proportion of paired genes, subtracted from one, yields an estimate for gene sharing differences (dsg) between prokaryotic genome pairs.

A clock-like rate of gene content change

Plots of ICD versus gene sharing distance (dsg) for pairwise genome comparisons reveal a highly significant positive correlation between evolutionary divergence and gene content differences (Fig. 2, left-hand panels; Fig. S1). Genome comparisons with similar IC gene sequences show little difference in gene content, whereas genome pairs with a high proportion of substitutions in amino acid sequences in pairwise IC gene comparisons show large differences in gene content. The plots for the raw data containing all genomes per higher taxon are shown in the left-hand panels of Fig. 2.

Figure 2. Correlation between ICD (x-axis) and gene sharing distance (dsg; y-axis) of pairwise genome comparisons in prokaryotic higher taxa. Every point represents one genome comparison and are colored based on the density of points. Spearman rank correlations and linear regressions between ICD and dsg are calculated for every higher taxon. The resulting r2 of the linear regression as well as the slope a are shown in the lower right corner (* = Spearman rank P-value ≤ .01). The upper part of the figure shows six highly sampled bacterial higher taxa and the lower part four archaeal ones. For each prokaryotic higher taxon, two data sets were examined. Every left-hand panel, per taxon, shows the correlation in the original dataset (u). The right-hand panel contains the correlation in the filtered data set (f). Thereby, only genomes were used whose genome size was within the average higher taxon genome size ± one standard deviation. Additionally, only one random genome per species was selected. Statistics of all prokaryotic higher taxa sampled are listed in Fig. S1.

In Alphaproteobacteria and Gammaproteobacteria, some genome comparisons have low gene content differences at high ICD, which is caused by the presence of many genomes of endosymbionts having highly reduced genomes (Moran and Bennett 2014, Wernegreen 2015) in these higher taxa (Fig. 2B and F). To mitigate the influence of reduced endosymbiont genomes, we generated data sets containing only genomes whose genome size was within the average higher taxon genome size ± 1 SD. Furthermore, to avoid the effect of oversampling and phylogenetic bias, only one genome per species was retained in the dataset. These plots for the data filtered for size and strains are shown in the right-hand panels of Fig. 2. For almost all higher taxa, the correlation remains almost unchanged relative to the unfiltered dataset (Fig. 2, right-hand panels; Fig. S1), as only outliers are removed.

We used the slopes of the linear regressions of plots of ICD versus dsg as an estimate for the relative gene flux rates per higher taxon (Table S3). The calculated gene flux rates represent the proportion of gene changes per 1% amino acid differences in the IC gene set per higher taxon. When we compared taxa with sufficient sample sizes, the rate of gene flux is similar across higher taxa, varying within Bacteria from 2.04 for Chlamydiae to 3.47 for Bacilli. In archaeal higher taxa, the lowest rate is 1.82 for Thermococcales and the highest is 3.27 for Halobacteriales. On average, the bacterial gene flux rate is 2.90 and the archaeal gene flux rate 2.57, yielding an average rate of 2.83% gene content change per 1% ICD. The similar gene flux rates across higher taxa show that a regular, and averaged over long timescales, almost clock-like behavior of gene content change in genomes relative to amino acid sequence divergence in the most universal and vertically inherited genes exists in prokaryotes. The archaeal dataset used here is limited to the largest groups of Euryarchaea and does not include genomes from Asgard archaea or fast-evolving groups like DPANN and CPR. The analysis is based on referenced cultured genomes, whereby Asgard sequences from enrichment cultures comprise a very small sample. The gene flux rates within archaea are similar across groups sampled here, future studies will reveal whether genomes from other archaeal taxa show similar or anomalous rates of change for ICD versus dsg.

Based on the average number of genes per higher taxon and the relative gene flux rates (Table S3) we calculated the number of different genes per substitution in the IC genes, to estimate absolute gene flux rates (Table S4). The average number of gene differences per substitution in the IC gene sequences is five for bacteria and four for archaea. However, in contrast to the relative gene flux rates (Table S3) the absolute rates vary substantially from 1.29 (Chlamydiae) to 9.03 (Betaproteobacteria). For archaea the highest rate of gene flux is found for Halobacteriales (6.27) and the lowest for Thermococcales (2.41; Table S4). The absolute gene flux rates are higher for free-living prokaryotes than for prokaryotic endosymbionts due to different genome sizes and the isolated nature of endosymbiont genome evolution (Clark et al. 1999, Shigenobu et al. 2000, van Ham et al. 2003, Moya et al. 2008, Moran and Bennett 2014, Martínez-Cano et al. 2015).

The question arises whether the similar gene flux rates are also present across different taxonomic levels. Therefore, we analysed the relationship between ICD and dsg at lower taxonomic levels. The correlation remained nearly constant except of taxa on the species level (Fig. S2a and Table S5) because ICD between genomes corresponding to the same species is too close (Fig. 3D). In order to better quantify gene flux at or near the species level, we used a genome-based measure for species-level divergence, ANI (Konstantinidies and Tiedje 2005, Goris et al. 2007, Wright and Baum 2018). We calculated the ANI for all genome pairs, using FastANI (Jain et al. 2018). For direct comparison to the ICD measure, we calculated the average nucleotide differences by subtracting the ANI from 1. This measure provides a higher resolution for the intraspecies genome comparisons compared to ICD, showing a range of 0%–12% 1-ANI at an ICD range of 0%–1% (Fig. 3). The relationship between 1-ANI and gene sharing distance from species- to order-level shows that the positive correlation holds for all taxonomic levels (Fig. S2b and Table S5). As noted in previous studies at the species and genera level (Wright and Baum 2018, Hassler et al. 2022), a correlation between dsg and 16S rRNA sequence divergence is not observed, as the latter tend to saturate at pairs distributed near 0.3 16S rRNA sequence divergence (Fig. S3).

Figure 3. Relationship between 1-ANI and dsg (A–C) as well as ICD and dsg (D–F) for E. coli, Escherichia, and Enterobacteriaceae. Every point represents one genome comparison and are colored based on the density of points. Spearman rank correlations and linear regressions between ICD and dsg are calculated for every taxon. The resulting r2 of the linear regression as well as the slope a are shown in the lower right corner (* = Spearman rank P-value ≤ .01). For each prokaryotic taxon, two data sets were examined. Every left-hand panel, per taxon, shows the correlation in the original dataset (u). The right-hand panel contains the correlation in the filtered data set (f). Thereby, only genomes were used whose genome size was within the average higher taxon genome size ± 1 SD. Additionally, only one random genome per species was selected.

By comparing the estimated relative rates of gene flux based on the slopes of the regression lines between ICD and dsg, we see that the rate of gene content change is not constant across taxonomic levels. The relative gene flux rates decrease with increasing phylogenetic divergence even if the differences are small except for the species level where the rates are much higher than on all other taxonomic levels. (Fig. S4a and Table S6). This could indicate that the mechanism of gene content change within species and between species are different, as suggested by Baumdicker and Kupczok (2023). From the genus level onwards, the slopes decrease slowly with increasing phylogenetic depth, which can result from sequence saturation in distance comparisons (or potentially from, genes that are gained, lost, and regained). The highly disparate gene flux rates within species and between species are also shown by comparing the relative gene flux rates based on the slopes of the linear regression between 1-ANI and dsg (Fig. S4b and Table S6). However, at taxonomic levels of genus, family and order the gene flux rates are similar.

Since the estimation of gene flux rates is based on linear regressions, which assume a normal distribution of residuals and low variance of y-values at all x-values, we performed residual analysis for the regression between ICD or 1-ANI and dsg for all higher taxa and on lower taxonomic levels (Figs S5–S8). The correlations of dsg with 1-ANI and ICD, though highly significant, were not strictly linear. High values of 1-ANI and ICD as well as low values of ICD show deviation from linearity (Figs S5 and S6). The deviation of high values of 1-ANI can be explained by the limitations of the method: ANI values can only be estimated up to about 20% 1-ANI because the values saturate (Fig. 4). ICD seems to be a robust measure at higher divergence. However, at values of ICD close to zero, no correlation is visible between ICD and dsg because no difference in ICD values is detectable (Fig. 4B). To investigate improved fit of the regression line between ICD and dsg, we performed logarithmic fits on genome comparisons of higher taxa as performed in previous studies (Touchon et al. 2009, Wielgoss et al. 2016, Andreani et al. 2017, Rocha 2018) including log–log transformations, log transformation of y-values (gene content differences, dsg) and log transformation of x-values (ICD). The log–log transformations of ICD and dsg values improved the r2 values form the regression line in most higher taxa sampled and thus the fit (Fig. S7), especially in the unfiltered dataset (Fig. S8). However, the residual analysis from linear regressions on log–log transformations were often worse than (or equal to) that observed when analysing the linear regressions between ICD and dsg without transformations (Figs S7 and S8). Linear regression on transformations of either ICD values or dsg values resulted in worse fit than that found for linear regression between ICD and dsg (Figs S7 and S8). Even though a power function (log–log transformation) gave a better fit of the regression line, we show the linear regressions of raw values without transformations, because neither are strictly linear (though obviously close to linear), and the residual analysis of raw values is better than (or not worse than) that observed with transformed data. Saturation of dsg estimates at higher ICD values might stem from genes that are gained, lost, and regained at appreciable frequency in the current sample.

Figure 4. Relationship between 1-ANI and ICD for all genome pairs with calculated ANI from FastANI. (A) Relationship between 1-ANI (x-axis) and ICD (y-axis) for all genome pairs with calculate ANI (n = 397 578). Every point represents one genome comparison and is colored based on the density of points. (B) Relationship between ICD or 1-ANI (x-axis) and dsg (y-axis) for prokaryotic taxon E. coli. Every point represents one genome comparison and is colored based on the density of points.

At low values of ICD in the raw dataset, vertical “icicle”-like structures, equidistant point groups (EPGs) can be observed (Fig. 2, left-hand panels). Genome pairs that generate EPGs have little or no ICD, but vary with respect to gene content differences. We colored the genome pairs according to their taxonomic affiliation to highlight the dependence of EPGs upon phylogenetic structure of the data (Fig. 5). The colors indicate the lowest common taxonomic level of genome pairs. Comparisons between genomes from the same species are colored pink and comparisons between genomes from the same phylum are colored gray. EPGs are only observed when genome pairs belong to the same species or genus (Fig. 5A, upper left enlarged section). When we compared one genome of a species with all other genomes corresponding to the same species, the differences in their gene content become greater from genome to genome, forming the EPG (Fig. 5B). This indicates that every EPG reflects the growth of the accessory genome as more genomes are added to the species-level pangenome. At higher evolutionary divergence, the EPGs are lost (Fig. 5A, lower right enlarged section), because as taxonomic divergence approaches the phylum level, the shared component of accessory genomes in pairwise comparisons also approaches zero. This leads to a narrower distribution of gene sharing values for taxa with shared genes in the accessory genome in comparison to other species with increasing sequence divergence. One might argue that EPGs can also be formed by low-quality genomes. However, all genomes used in this portion of our analysis are complete genomes. The presence or absence of the EPGs represents the phylogenetic structure of the data, even though no trees were used for the calculation of the ICD. In the filtered dataset (Fig. 2, right-hand panels), the EPGs are no longer visible due to lacking intraspecies comparisons, as only one strain per species was retained in the dataset.

Figure 5. Correlation of ICD and dsg of Chloroflexi. (A) Enlarged view of the correlation between ICD (x-axis) and dsg (y-axis) of the bacterial taxon Chloroflexi. Points represent genome comparisons and are colored by the lowest equal taxon of the genomes. Equidistant point groups (EPGs) are shown in intraspecies comparisons in the upper left enlarged section. In lower taxonomic groups, like intraspecies and intragenus comparisons, EPGs are formed by genome comparisons that have a nearly identical average frequency of substitutions in IC genes but differ in their gene content. In intraphylum comparisons, the EPGs are lost (lower right enlarged section). Panel (B) shows possible tree structures that could represent EPGs of intraspecies comparisons.

The y-axis intercept of the regression lines in Fig. 2 rarely goes through the origin. In most higher taxa, y-axis intercept assumes a value between 0.2 and 0.4. This means that in pairwise genome comparisons with zero ICD, the differences in pairwise gene content are 20%–40%. That is, at zero ICD, roughly 60%–80% of genes in a given comparison are shared. These 60%–80% paired genes are similar to the 70% value of DNA–DNA sequence hybridization (DDH), used in the 1960s and 1970s to delineate prokaryotic species (Wayne et al. 1987). The value of DDH correlates positively with the proportion of conserved DNA sequences (>90% sequence identity) in genomes pairs (Goris et al. 2007). In the present study, the y-axis intercept reflects the gene content divergence in genome pairs with identical ICD. In DDH, nonconserved sequences (<90% sequence identity) correspond to the proportion of nonshared genes in pangenome analysis, which belong to the accessory genome (Tettelin et al. 2005). In previous species level pangenome analyses, the proportion of accessory genes per genome typically varied between 0.2 and 0.5, calculated as 1 minus the proportion of species-level core genes in prokaryotic genomes (Tettelin et al. 2005, Rasko et al. 2008, Schoen et al. 2008, Scaria et al. 2010, van Schaik et al. 2010, Budroni et al. 2011, Park et al. 2019). This value is close to the average proportion of gene content differences between genomes having zero ICD across higher taxa found here.

Thus, the y-axis intercept provides a rough estimate for the proportion of accessory genes in genomes of the same species within the higher taxon. To test this, we calculated average species-specific accessory genome proportions per higher taxon based on the proportion of unpaired genes in pairwise genome comparisons corresponding to the same species. To avoid bias from small samples, only higher taxa represented by at least 10 strains, 2 species and a P-value ≤ .05 from the correlations calculated in Fig. 2 and Fig. S1 were used. The average species accessory genome proportions as well as the y-axis intercept were taken from the filtered dataset to avoid phylogenetic bias. The average species-level accessory genome proportions and the y-axis intercept values inferred from the regression lines are quite similar (Table S7) and correlate positively at P ≤ .01 (Fig. S9 and Table S7). That is, the y-axis intercept inferred from the entire history of the higher taxon delivers an estimate for the average accessory genome proportion in current genomes. This is interesting per se, but it is also evidence for the existence of pangenomes throughout the evolutionary history of prokaryotes.

To analyse the influence of cloud genes, genes that are present in only a few genomes or in only one genome, we excluded protein families containing proteins from less than 5%, 10%, 15%, or 20% of all genomes from the gene sharing distance calculation. We then compared the statistical results of the correlation and regression between the ICD and each filtered gene sharing distance dataset and the dataset including all protein families. For the Spearman rank correlation coefficient (rs), the slope a, and r2 of the regression line, no differences are shown (Fig. S10 and Table S8). The y-axis intercept decreases as more genes are excluded. This makes sense because the y-axis intercept represents the proportion of accessory genes per genome, including genes that are only present in a few genomes or only in one. Therefore, we cannot see any effect of cloud genes on the analysis.

The accessory genome is always more affected by LGT

Accessory genes are generally thought to be more often transferred between genomes and to contribute to functions important for the adaptation to a specific niche, whereas core genes are more conserved and tend to include housekeeping functions (Tettelin et al. 2005, 2008, Kung et al. 2010, Vernikos et al. 2015, Brockhurst et al. 2019). To see whether this aspect is captured by our approach, we plotted the average verticality of genes that are present in both genomes (paired) versus the average verticality of those that are missing in one genome of the pair (unpaired) for all genome pairs (Fig. 6). The verticality distributions between paired (red) and unpaired genes (gray) are significantly different for all genome comparisons corresponding to the same higher taxon (P = .0; Table S9) and are furthermore nonoverlapping sets. With evolutionary divergence, the verticality of paired genes as well as for unpaired genes increases. However, the slope for paired genes is much higher (a = 22.57) than the slope for unpaired genes (a = 2.89). The average verticality of paired genes increases because the more different two genomes are, the fewer genes they share. These genes are highly conserved and furthermore exhibit high verticality. The higher y-axis intercept for paired genes (paired genes y-axis intercept = 1.84, unpaired genes y-axis intercept = 0.34) shows that at zero ICD, accessory (unpaired) genes are always less vertical than core (paired) genes. This also holds for higher ICD, because all 124 447 compared genome pairs have a higher verticality for paired protein families than for unpaired, except for one comparison between two Mollicutes genomes (Spiroplasma mirum and Spiroplasma atrichopogonis). The analysis of the unfiltered dataset gave the same result (Fig. S11 and Table S9), whereby 88 of 2 229 958 comparisons have a higher average verticality for unpaired protein families than for paired ones, yet always belonging to the same species or genus.

Figure 6. Comparison between average verticality (y-axis) of paired and unpaired protein families and divergence in IC genes (x-axis) in the filtered dataset. Paired protein families are defined as protein families in which both genomes of a genome comparison were represented by at least one protein. If only one genome is present in the protein family, the protein family corresponds to the group of unpaired protein families. For every genome pair two points are plotted. The average verticality of paired genes and the average verticality of unpaired genes in the compared genome pair relative to the sequence divergence in the IC genes from the compared genomes. The color of the points represents the density in the region near the point. Spearman rank correlations and linear regressions between ICD and dsg are calculated for paired protein families and unpaired protein families (number of genome pairs = 124 447). The resulting r2 of the linear regression as well as the slope a are shown next to the regression lines (* = Spearman rank P-value ≤ .01).

The nonoverlapping verticality distribution for paired and unpaired genes in every genome comparison (Fig. 6) might suggest that the two gene sets are drawn from different samples, but because one and the same gene can be paired (core) in one comparison but unpaired in another, the two sets delineated in Fig. 6 simply indicate that the accessory genome is vastly more prone to LGT than the core genome in every genome comparison. This does not mean that core genes cannot be affected by LGT, however, it shows that accessory genes are far more readily transferred (Nesbø et al. 2001). An analysis of the frequency of functional categories (KO) from KEGG (Fig. S12; Kanehisa et al. 2017) revealed that, as expected genes belonging to the category translation were more frequent in the paired set than in the unpaired. The most common KEGG B functional annotation for unpaired genes is virus information processing (Table S10). It is known that prokaryotic genomes have core and accessory components (Medini et al. 2005, Tettelin et al. 2005, 2008, Vernikos et al. 2015, Brockhurst et al. 2019). What Fig. 6 shows, is that in all comparisons, the accessory (unpaired) component is always comprised of genes that are transferred more frequently (they are less vertical). This is consistent with the observation that the tendency for a gene to undergo LGT is restricted by the presence of a preexisting copy (Nagies et al. 2020). It is also consistent with the early observation by Roger Milkman (1996) that “The structure of genetic variation in a bacterial species is the result of recombination superimposed upon the repeated formation and spread of clones”, whereby recombination in prokaryotes is never reciprocal and need not require donor and acceptor to belong to the same species, genus, phylum, or domain.

High-quality MAGs show similar rates of gene flux

MAGs have gained much attention in recent years, as cultured genomes are estimated to account for only a very small fraction of all prokaryotic organisms on earth (Lok 2015). Since many prokaryotes cannot be cultivated or only under difficult or expensive laboratory conditions, DNA samples are taken from the environment and are categorized into genomes—MAGs—using methods that assemble genome sized bins (Garza and Dutilh 2015). However, these methods are far from perfect (Garza and Dutilh 2015). As a consequence, programs have been developed to estimate the quality of MAGs. One of the most widely used methods for cross-checking MAGs is CheckM (Parks et al. 2015), which uses (i) marker genes that are specific to the genome’s lineage, inferred from a reference genome tree, and (ii) marker genes that are usually single-copy in genomes of cultivated species, to calculate the completeness and the redundancy (contamination) of the MAG, respectively.

The two parameters estimated by CheckM are related to—but not identical to—the genomic parameters that we have investigated here, even though we do not use a reference tree or score genes as single copy. ICD contains evolutionary information (pairwise distances between genomes), but involves no tree, while dsg gives information about how many genes are shared between two genomes, but not which genes, specifically, are shared.

Because a very similar pattern and correlation between ICD and dsg is observed across a wide spectrum of genomes of cultured strains, it was of interest to see how MAGs appear when viewed from the standpoint of ICD versus dsg. For that, we clustered a dataset of 2872 MAGs that contained a spectrum of assemblies with varying quality. The metagenomic data were clustered using the same methods as for cultured strains, values of ICD and dsg were calculated accordingly. Completeness and contamination values were calculated using CheckM (Park et al. 2015) as the quality measure, defined by Parks et al. (2017), providing a filter for the quality of metagenome assemblies. This quality measure is calculated by subtracting five times the contamination from the completeness (Parks et al. 2017), a procedure that at high stringencies yields metagenome assemblies resembling genomes of cultured strains. We filtered the MAGs for different qualities, ranging from no filter (at least 0% quality) to at least 90% quality and performed correlation and linear regression analysis between ICD and dsg for each dataset (Fig. 7).

Figure 7. Correlation between ICD and dsg using MAGs with a minimum quality ranging from at least 0% to at least 90%. Quality is defined as completeness minus five times the contamination and is shown at the right side over every panel (Parks et al. 2015). The slope (a) and r2 from linear regression between ICD and dsg are shown in every lower-right corner (* = Spearman P ≤ .01). Points are colored based on their density.

At quality thresholds <50%, there is a visible tendency of pairwise comparisons from assemblies with ICD ≈ 0 to span the full range of 0 < dsg < 0.8. This reflects the existence of unusually disparate gene collections in low quality assemblies that are closely related by the measure of ICD divergence (ICD <0.05). Also at quality thresholds <50%, there is a tendency for assemblies to exhibit dsg ≈ 1 across all values of ICD. This reflects a class of assemblies that have gene collections more disparate than that observed for cultured strains, independent of ICD. At the same time, values of ICD often exceed 0.35 for low quality assemblies, something not observed for cultures strains, even in bacterial–archaeal comparisons (Fig. 8), suggesting that in low quality assemblies, IC genes may contain erroneous sequence information. With increasing quality filter stringency, however, MAGs reflect the properties of cultured strains with respect to gene flux rate.

Figure 8. Correlation between ICD and gene sharing distance (dsg) for prokaryotic genome pairs in higher taxonomic ranks. Correlation between ICD (x-axis) and dsg (y-axis) of bacterial and archaeal genome pairs corresponding to the same higher taxon (within higher taxa), bacterial and archaeal genome pairs corresponding to different higher taxa (between higher taxa), and intersuperkingdom genome comparisons. Every point represents one pairwise genome comparison and is colored based on affiliation of the compared genomes to bacteria (blue), archaea (green), or both (purple). Lighter blue and green show pairwise comparisons between genomes corresponding to the same higher taxon. Darker blue and green shows comparisons corresponding to different higher taxa.

In Fig. 7 it is clearly shown that r2 from linear regression between ICD and dsg increases with the quality of metagenomes data. The highest value is calculated for the dataset using MAGs with at least 90% quality (Fig. 7J). Values for 100% quality are not shown because only one MAG was estimated to have a quality of 100%. As r2 in cultured genomes generally varies between 0.60 and 0.95 (Fig. S1) a reliable prediction of statistical values between the correlation and regression of ICD and dsg for dataset including MAGs is only possible by using high-quality MAGs of more than 80%. In lower-quality datasets an overrepresentation of genome pairs that exhibit extremely high ICD values is shown (Fig. 7A–F) as well as genome pairs with low ICD but high gene content differences (dsg; Fig. 7A–D). The similar properties of cultured genomes (Fig. 2 and Table S3) and high-quality MAGs (Fig. 7I–J) including r2, the relative gene flux rates (a) and the y-axis intercept from the regression lines between ICD and dsg, indicates that a constant gene flux rate across higher prokaryotic taxa, with high gene flux within the accessory genome and lower flux in the core, also applies to MAGs, if the quality of the data is high.

Pangenome structure throughout prokaryote evolution

On the basis of pairwise genome comparisons, we can observe a clock-like rate of gene turnover within higher taxa, which can be used to estimate species-level accessory and core genome sizes and which properties reflect the phylogenetic structure of the data. Since LGT events are also detectable between superkingdoms, primarily between archaea and bacteria (Rest and Mindell 2003, Gophna et al. 2004, Boto 2010), the last universal common ancestor (LUCA) (Weiss et al. 2016) might have had the ability to undergo LGT, and therefore a genome containing core genes with mainly vertical inheritance and accessory genes that are in permanent flux (Woese 2002, Nagies et al. 2020).

To test whether LUCA might have had similar rates of gene flux as found in current genomes, we expanded the analysis from Fig. 2 to all pairwise genome comparisons within prokaryotes (Fig. 8), including bacterial and archaeal genome comparisons corresponding to the same higher taxon (bacteria number of genome pairs = 123 940; archaea number of genome pairs = 507), genome comparisons corresponding to different higher taxa (bacteria number of genome pairs = 1 013 846; archaea number of genome pairs = 6753) and intersuperkingdom comparisons between archaea and bacteria (number of genome pairs = 182 589) from the dataset filtered for genome sizes and strains. The rates for bacterial and archaeal comparisons corresponding to only one higher taxon (bacterial slope = 2.51, archaeal slope = 2.07) are close to the average rates calculated in Table S3, with small differences attributable to genome comparisons corresponding to higher taxa with small sample sizes. The correlation between ICD and dsg also holds for deep genome comparisons corresponding to different higher taxa and intersuperkingdom comparisons (P = .0) as well as the verticality distribution for paired (high verticality) and unpaired genes (low verticality; Fig. S13). This indicates that the basic mechanisms of gene content change in prokaryotic genomes—rapid turnover of genes in the accessory genome and slow gene flux in the core genome—has been operational since the divergence of the bacterial and archaeal lineages. However, the slope decreases, yet we have not applied any kind of correction for multiple events on either axis in order to obtain low-parameter estimates. In pairwise comparisons of bacterial genomes corresponding to different higher taxa, the estimated gene flux rate is 1.80% gene content differences at 1% ICD, between superkingdoms 0.78%.

A persistent rate of flux relative to ICD could indicate that a pangenome structure has persisted throughout prokaryotic history. Our analyses provide neither an estimate for the size of LUCA’s pangenome nor an estimate of which genes, specifically, it contained. The data do, however, reveal that the rate of gene content change is roughly constant as far back as our sample probes the bacterial and archaeal lineages, indicating that the average rate of gene content change that we observe today has existed throughout prokaryotic evolution all the way back to the first free-living cells. This, in turn, traces a pangenome structure—a core genome and an accessory genome—back to the common ancestor of bacteria and archaea.

A possible application of our procedure to pangenomes concerns their impact on the origins of organelles (Esser et al. 2007, Ku et al. 2015). We embarked upon estimation of the number of genes that have persisted in any individual bacterial genome that shared common ancestry with the mitochondrial ancestor or with the cyanobacterial ancestor. However, the degree of sequence divergence between the bacterial and nuclear homologues of organelle-encoded or organelle-derived IC genes (0.4–0.5; Table S11–Supplemental Results) exceeds that observed from prokaryotic genomes (0.35; Fig. 8) such that the current method does not permit interpolation.

Conclusion

Using universally applicable proxies for genome divergence, we found very similar rates of gene flux across all higher prokaryotic taxa sampled. Gene content divergence in prokaryotic accessory genomes correlates linearly two measures of sequence divergence, ANI (nucleotides) and ICD (amino acids), which in turn indicates a clock-like rate of gene flux across higher prokaryotic taxa. The analysis shows that a universal rate of gene flux is present at all taxonomic levels, with higher gene flux within species than across higher taxonomic groups. Baumdicker and Kupczok (2023) proposed that within-population gene gain spreads genes within the population, and between-population gene gain leads to the acquisition of novel genes. This is consistent with our observation that within species, gene flux rates differ across taxa, while between species the rates remain similar.

Our use of a measure for lineage divergence that is generally applicable across comparisons of all higher prokaryotic taxa allowed us to look deeper into time than the genus or species level, where ANI performs well. The generally linear correlation between gene content differences and ICD holds as far back as any of these lineages can be traced, including the archaea–bacteria split. This indicates that the rate of genome flux and the pangenome structure of modern prokaryotic genomes is as ancient as prokaryotes themselves. This is an independent line of evidence indicating that LGT is a natural component of genome evolution in prokaryotes and has been in operation throughout the entirety of their evolution.

There is an ongoing debate about the main evolutionary forces shaping gene content in pangenomes (Vos et al. 2015). Here, we have observed similar and constant rates of gene flux across prokaryotic taxa. By analogy to nucleotide substitution rates, this could suggest a neutral or nearly neutral nature of gene flux in prokaryotes (Baumdicker and Kupczok 2023), compatible with Kimura’s neutral theory (Kimura 1968) or the nearly neutral theory proposed by Ohta (1973). However, similar gene flux rates across species can also be obtained using selective models (Barrick et al. 2009). Neither mechanism can be expected to explain all fixation events. Both neutral and selective mechanisms likely influence the rate and frequency of genes fixed by flux between the accessory genomes of prokaryotes.

Supplementary Material

fnae068_Supplemental_Files

Acknowledgements

We thank Cerys Viktoria Wilke for downloading prokaryotic 16S rRNA sequences and Tal Dagan for her helpful comments and suggestions on the manuscript. Computational infrastructure and support were provided by the Centre for Information and Media Technology at Heinrich Heine University Düsseldorf.

Conflict of interest

The authors declare no conflict of interest.

Funding

This work was supported by the European Research Council (ERC) under the European Union’s Horizon 2020 Research and Innovation Program (grant agreement number 101018894 to W.F.M.). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Data Availability

The data underlying this article are available in the article, in its online supplementary material and in https://uni-duesseldorf.sciebo.de/s/YiiyDz57sDala1t.
==== Refs
References

Albalat R , CañestroC. Evolution by gene loss. Nat Rev Genet. 2016;17 :379–91.27087500
Andreani NA , HesseE, VosM. Prokaryote genome fluidity is dependent on effective population size. ISME J. 2017;11 :1719–21.28362722
Arnold BJ , HuangI, HanageWP. Horizontal gene transfer and adaptive evolution in bacteria. Nat Rev Microbiol. 2022;20 :206–18.34773098
Barrick JE , YuDS, YoonSHet al. Genome evolution and adaptation in a long-term experiment with Escherichia coli. Nature. 2009;461 :1243–7.19838166
Baumdicker F , KupczokA. Tackling the pangenome dilemma requires the concerted analysis of multiple population genetic processes. Genome Biol Evol. 2023;15 . 10.1093/gbe/evad067.
Boto L . Horizontal gene transfer in evolution: facts and challenges. Proc R Soc B. 2010;277 :819–27.
Brockhurst MA , HarrisonE, HallJPJet al. The ecology and evolution of pangenomes. Curr Biol. 2019;29 :R1094–103. 10.1016/j.cub.2019.08.012.31639358
Brueckner J , MartinWF. Bacterial genes outnumber archaeal genes in eukaryotic genomes. Genome Biol Evol. 2020;12 :282–92.32142116
Buchfink B , XieC, HusonDH. Fast and sensitive protein alignment using DIAMOND. Nat Methods. 2015;12 :59–60.25402007
Budroni S , SienaE, Dunning HotoppJCet al. Neisseria meningitidis is structured in clades associated with restriction modification systems that modulate homologous recombination. Proc Natl Acad Sci USA. 2011;108 :4494–9.21368196
Charlebois RL , DoolittleWF. Computing prokaryotic gene ubiquity: rescuing the core from extinction. Genome Res. 2004;14 :2469–77.15574825
Clark MA , MoranNA, BaumannP. Sequence evolution in bacterial endosymbionts having extreme base compositions. Mol Biol Evol. 1999;16 :1586–98.10555290
Dagan T , MartinWF. The tree of one percent. Genome Biol. 2006;7 :118. 10.1186/gb-2006-7-10-118.17081279
Esser C , MartinW, DaganT. The origin of mitochondria in light of a fluid prokaryotic chromosome model. Biol Lett. 2007;3 :180–4.17251118
Garg SG , KapustN, LinWet al. Anomalous phylogenetic behavior of ribosomal proteins in metagenome-assembled Asgard archaea. Genome Biol Evol. 2021;13 . 10.1093/gbe/evaa238.
Garza DR , DutilhBA From cultured to uncultured genome sequences: metagenomics and modeling microbial ecosystems. Cell Mol Life Sci. 2015;72 :4287–308.26254872
Gophna U , CharleboisRL, DoolittleWF. Have archaeal genes contributed to bacterial virulence?. Trends Microbiol. 2004;12 :213–9.15120140
Goris J , KonstantinidisKT, KlappenbachJAet al. DNA-DNA hybridization values and their relationship to whole-genome sequence similarities. Int J Syst Evol Microbiol. 2007;57 :81–91.17220447
Hansmann S , MartinW. Phylogeny of 33 ribosomal and six other proteins encoded in an ancient gene cluster that is conserved across prokaryotic genomes: influence of excluding poorly alignable sites from analysis. Int J Syst Evol Microbiol. 2000;50 :1655–63.10939673
Hao W , GoldingGB. The fate of laterally transferred genes: life in the fast lane to adaptation or death. Genome Res. 2006;16 :636–43.16651664
Hassler HB , ProbertB, MooreCet al. Phylogenies of the 16S rRNA gene and its hypervariable regions lack concordance with core genome phylogenies. Microbiome. 2022;10 . 10.1186/s40168-022-01295-y.
Haudiquet M , De SousaJM, TouchonMet al. Selfish, promiscuous and sometimes useful: how mobile genetic elements drive horizontal gene transfer in microbial populations. Phil Trans R Soc B. 2022;377 . 10.1098/rstb.2021.0234.
Jain C , Rodriguez-RLM, PhillippyAMet al. High throughput ANI analysis of 90 K prokaryotic genomes reveals clear species boundaries. Nat Commun. 2018;9 . 10.1038/s41467-018-07641-9.
Kanehisa M , FurumichiM, TanabeMet al. KEGG: new perspectives on genomes, pathways, diseases and drugs. Nucleic Acids Res. 2017;45 :D353–61.27899662
Katoh K . MAFFT: a novel method for rapid multiple sequence alignment based on Fast Fourier transform. Nucleic Acids Res. 2002;30 :3059–66.12136088
Kimura M . Evolutionary rate at the molecular level. Nature. 1968;217 :624–6.5637732
Kislyuk AO , HaegemanB, BergmanNHet al. Genomic fluidity: an integrative view of gene diversity within microbial populations. BMC Genomics. 2011;12 . 10.1186/1471-2164-12-32.
Konstantinidis KT , TiedjeJM. Genomic insight that advance the species definition for prokaryotes. Proc Natl Acad Sci USA. 2005;102 :2567–72.15701695
Ku C , Nelson-SathiS, RoettgerMet al. Endosymbiotic origin and differential loss of eukaryotic genes. Nature. 2015;524 :427–32.26287458
Kung VL , OzerEA, HauserAR. The accessory genome of Pseudomonas aeruginosa. Microbiol Mol Biol Rev. 2010;74 :621–41.21119020
Lawrence JG , OchmanH. Molecular archaeology of the Escherichia coli genome. Proc Natl Acad Sci USA. 1998;95 :9413–7.9689094
Lok C . Mining the microbial dark matter. Nature. 2015;522 :270–3.26085253
Marri PR , HaoW, GoldingGB. Gene gain and gene loss in Streptococcus: is it driven by habitat?. Mol Biol Evol. 2006;23 :2379–91.16966682
Marri PR , HaoW, GoldingGB. The role of laterally transferred genes in adaptive evolution. BMC Evol Biol. 2007;7 :1–14.17214884
Martínez-Cano DJ , Reyes-PrietoM, Martínez-RomeroEet al. Evolution of small prokaryotic genomes. Front Microbiol. 2015;5 . 10.3389/fmicb.2014.00742.
Medini D , DonatiC, TettelinHet al. The microbial pan-genome. Curr Opin Genet Dev. 2005;15 :589–94.16185861
Milkman R . Recombinational exchange between clonal populations. In: NeidhardtFC (ed.), Escherichia coli and Salmonella: Cellular and Molecular Biology. 2nd edn. Washington: American Society for Microbiology Press, 1996.
Mira A , OchmanH, MoranNA. Deletional bias and the evolution of bacterial genomes. Trends Genet. 2001;17 :589–96.11585665
Moran NA , BennettGM. The tiniest tiny genomes. Annu Rev Microbiol. 2014;68 :195–215.24995872
Moya A , PeretóJ, GilRet al. Learning how to live together: genomic insights into prokaryote–animal symbioses. Nat Rev Genet. 2008;9 :218–29.18268509
Nagies FSP , BruecknerJ, TriaFDKet al. A spectrum of verticality across genes. PLoS Genet. 2020;16 :e1009200. 10.1371/journal.pgen.1009200.33137105
Nesbø CL , BoucherY, DoolittleWF. Defining the core of nontransferable prokaryotic genes: the euryarchaeal core. J Mol Evol. 2001;53 :340–50.11675594
Nowell RW , GreenS, LaueBEet al. The extent of genome flux and its role in the differentiation of bacterial lineages. Genome Biol Evol. 2014;6 :1514–29.24923323
O’Leary NA , WrightMW, BristerJRet al. Reference sequence (RefSeq) database at NCBI: current status, taxonomic expansion, and functional annotation. Nucleic Acids Res. 2016;44 :D733–45.26553804
Ohta T . Slightly deleterious mutant substitutions in evolution. Nature. 1973;246 :96–98.4585855
Park SC , LeeK, KimYOet al. Large-scale genomics reveals the genetic characteristics of seven species and importance of phylogenetic distance for estimating pan-genome size. Front Microbiol. 2019;10 . 10.3389/fmicb.2019.00834.
Parks DH , ImelfortM, SkennertonCTet al. CheckM: assessing the quality of microbial genomes recovered from isolates, single cells, and metagenomes. Genome Res. 2015;25 :1043–55.25977477
Parks DH , RinkeC, ChuvochinaMet al. Recovery of nearly 8,000 metagenome-assembled genomes substantially expands the tree of life. Nat Microbiol. 2017;2 :1533–42.28894102
Rasko DA , RosovitzMJ, MyersGSet al. The pangenome structure of Escherichia coli: comparative genomic analysis of E. coli commensal and pathogenic isolates. J Bacteriol. 2008;190 :6881–93.18676672
Rest JS , MindellDP. Retroids in archaea: phylogeny and lateral origins. Mol Biol Evol. 2003;20 :1134–42.12777534
Rice P , LongdenL, BleasbyA. EMBOSS: the European molecular biology open software suite. Trends Genet. 2000;16 :276–7.10827456
Rinke C , SchwientekP, SczyrbaAet al. , Insight into the phylogeny and coding potential of microbial dark matter. Nature. 2013;499 :431–7.23851394
Rivera MC , JainR, MooreJEet al. Genomic evidence for two functionally distinct gene classes. Proc Natl Acad Sci USA. 1998;95 :6239–44.9600949
Rocha EPC . Neutral theory, microbial practice: challenges in bacterial population genetics. Mol Biol Evol. 2018;35 :1338–47.29684183
Scaria J , PonnalaL, JanvilisriTet al. Analysis of ultra low genome conservation in Clostridium difficile. PLoS One. 2010;5 :e15147. 10.1371/journal.pone.0015147.21170335
Schoen C , BlomJ, ClausHet al. Whole-genome comparison of disease and carriage strains provides insights into virulence evolution in Neisseria meningitidis. Proc Natl Acad Sci USA. 2008;105 :3473–8.18305155
Shigenobu S , WatanabeH, HattoriMet al. Genome sequence of the endocellular bacterial symbiont aphids Buchnera sp. APS. Nature. 2000;407 :81–86.10993077
Stull GW , QuX, Parins-FukuchiCet al. Gene duplications and phylogenomic conflict underlie major pulses of phenotypic evolution in gymnosperms. Nat Plants. 2021;7 :1015–25.34282286
Tettelin H , MasignaniV, CieslewiczMJet al. Genome analysis of multiple pathogenic isolates of Streptococcus agalactiae: implications for the microbial “pan-genome.”. Proc Natl Acad Sci USA. 2005;102 :13950–5.16172379
Tettelin H , RileyD, CattutoCet al. Comparative genomics: the bacterial pan-genome. Curr Opin Microbiol. 2008;11 :472–7.19086349
Touchon M , HoedeC, TenaillonOet al. Organised genome dynamics in the Escherichia coli species results in highly diverse adaptive paths. PLoS Genet. 2009;5 :e1000344. 10.1371/journal.pgen.1000344.19165319
Treangen TJ , RochaEPC. Horizontal transfer, not duplication, drives the expansion of protein families in prokaryotes. PLoS Genet. 2011;7 :e1001284. 10.1371/journal.pgen.1001284.21298028
Tria FDK , MartinWF. Gene duplications are at least 50 times less frequent than gene transfers in prokaryotic genomes. Genome Biol Evol. 2021;13 . 10.1093/gbe/evab224.
UniProt Consortium . UniProt: the universal protein knowledgebase in 2021. Nucleic Acids Re. 2021;49 :D480–9.
van Dongen S .Graph clustering via a discrete uncoupling process.Siam Journal on Matrix Analysis and Applications. 2008;30 :121–41.
van Ham RCHJ , KamerbeekJ, PalaciosCet al. Reductive genome evolution in Buchnera aphidicola. Proc Natl Acad Sci USA. 2003;100 :581–6.12522265
van Schaik W , TopJ, RileyDRet al. Pyrosequencing-based comparative genome analysis of the nosocomial pathogen Enterococcus faecium and identification of a large transferable pathogenicity island. BMC Genomics. 2010;11 :1–18.20044946
Vernikos G , MediniD, RileyDRet al. Ten years of pan-genome analyses. Curr Opin Microbiol. 2015;23 :148–54.25483351
Vos M , HesselmanMC, BeekTATet al. Rates of lateral gene transfer in prokaryote: high but why?. Trends Microbiol. 2015;23 :598–605.26433693
Wayne LG , BrennerDJ, ColwellRRet al. Report of the ad hoc committee on reconciliation of approaches to bacterial systematics. Int J Syst Evol Microbiol. 1987;37 :463–4.
Weiss MC , SousaFL, MrnjavacNet al. The physiology and habitat of the last universal common ancestor. Nat Microbiol. 2016;1 :16116. 10.1038/nmicrobiol.2016.116.27562259
Wernegreen JJ . Endosymbiont evolution: predictions from theory and surprises from genomes. Ann NY Acad Sci. 2015;1360 :16–35.25866055
Wielgoss S , DidelotX, ChaudhuriRRet al. A barrier to homologous recombination between sympatric strains of the cooperative soil bacterium Myxococcus xanthus. ISME J. 2016;10 :2468–77.27046334
Woese CR . On the evolution of cells. Proc Natl Acad Sci USA. 2002;99 :8742–7.12077305
Wolf YI , KooninEV. A tight link between orthologs and bidirectional best hits in bacterial and archaeal genomes. Genome Biol Evol. 2012;4 :1286–94.23160176
Wolf YI , MakarovaKS, LobkovskyAEet al. Two fundamentally different classes of microbial genes. Nat Microbiol. 2016;2 :1–6.
Wright ES , BaumDA. Exclusivity offers a sound yet practical species criterion for bacteria despite abundant gene flow. BMC Genomics. 2018;19 :1–12.29291715
