Protein-Level Selection and Microbial Genome Size: Model Prediction and Genomic Data.

It has been shown previously that, in prokaryotes, genome size is negatively and significantly correlated with dN/dS, suggestive of a stronger protein-level selection in larger genomes (13). For the purpose of this work, we reproduced these findings on a substantially expanded set of groups of closely related bacterial and archaeal genomes compiled in the latest version of the Alignable Tight Genomic Cluster (ATGC) database (14) (Fig. 1). Specifically, genome size (measured by the number of genes) and dN/dS values are correlated with Spearman’s rank correlation coefficient ρ=−0.397. These findings imply that, in addition to local influences of different environments and lifestyles, there are underlying universal components of gene gain and loss probabilities with respect to genome size that are common across all microbes.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 1.

Genome size and selection in prokaryotes. Mean observed genome sizes (number of genes) are plotted against (A) the selection strength and (B) the estimated effective population size for each ATGC. Error bars correspond to 1 SD. ATGCs with more than 20 species are indicated in orange. The mean genome size, calculated using the mathematical model, is indicated in B by a solid red line. The model was applied with a constant selection coefficient, the gain and loss rates of Eq. 11, and optimized parameters shown in Table 1.

To quantitatively analyze the relations between selection strength and genome size, a mathematical model was implemented using the Moran process framework (15). In the model, a genome is represented as a collection of x genes, and the variation in genome size is attributed to stochastic gain and loss of genes (Fig. 2). New genes are assumed to be acquired at rate α, which in principle, depends on x and accounts for the probability that the acquired genetic material is inserted in a nondeleterious locus and expressed. Genes can be acquired by either horizontal gene transfer (in the model, assumed to draw genes from an infinite pool) or duplication, and α represents a weighted average over all possible evolutionary processes that lead to gene gain. Reconstructions of genome evolution performed on the ATGCs indicate that the contribution of horizontal gene transfer by far exceeds that of duplication (14). To model gene loss, the deleted gene is picked from the genome at random with the deletion rate β, which analogous to α, in principle, depends on x. Acquisitions and deletions of genetic material are either fixed or eliminated stochastically, with a fixation probability that depends on the selection coefficient s, which is a function of genome size as well s=s(x). Specifically, s(x) denotes the mean selective advantage of an individual with x+1 genes over an individual with x genes. Conveniently, selection coefficients associated with gene acquisition and gene deletion are related as (Methods)sdeletion=−sacquisition.[1]The fixation probability F can be approximated byF(x)≈s(x)1−e−Nes(x),[2]where Ne is the effective population size (16). The effective population size is estimated for each prokaryotic group from the dN/dS ratio (Methods). Fixed acquisition and deletion events, hereinafter “gene gain” and “gene loss,” respectively, occur with probabilities given by the multiplication of the event rate and the fixation probability:P+(x)=α(x)⋅s(x)1−e−Nes(x)[3]andP−(x)=β(x)⋅s(x)⋅e−Nes(x)1−e−Nes(x),[4]where Eq. 2 is used together with the relation of Eq. 1.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 2.

The model of genome evolution. (Upper) Illustration of the mathematical model of genome size evolution. The number of genes changes stochastically via gene gains and losses, which occur with probabilities P+ and P−, respectively. (Lower) Gene gain (solid purple curve) and loss (dashed purple curve) probabilities. Gain and loss probabilities are equal at x0, indicated by a vertical solid line. For values of x smaller than x0, the gain probability is larger than the loss probability, and therefore, the extremum of the steady-state genome size distribution (orange curve) at x0 is a maximum. The distribution is moderately skewed, and the mean value of x, indicated by a vertical dashed line, is close to the value of x0 (indicated by the solid vertical line).

Gain and loss probabilities govern the stochastic genome size dynamics, with the intuitive relation (derivation is given in Methods)x˙=P+(x)−P−(x),[5]meaning that, at each time step, there is a probability P+ that one gene is gained and a probability P− that one gene is lost, where time is measured in fixation time units (a unit is the time that it takes, on average, for an acquired gene to be fixed in the population) rather than generations. Because of the stochastic fluctuations, genome sizes form a distribution. If, for a certain genome size value x0, the gain and loss probabilities are equal, there is a steady-state genome size distribution with an extremum point at x0 as illustrated in Fig. 2 and SI Appendix, Fig. S1. The extremum point depends on the deletion–acquisition rates ratio r(x),r(x)=β(x)α(x),[6]and the equality of gain and loss probabilities impliess(x0)=lnr(x0)Ne.[7]As reflected in the above equation (and intuitively), steady state is possible only when the more frequent event, gene acquisition or gene deletion, is counterselected at genome size x0 (Fig. 3). For example, positive s(x0) implies selective advantage for larger genomes. In this case, steady state is possible only when deletion rates are higher than acquisition rates [i.e., r(x0)>1]. Formally, the sign of s(x0) in Eq. 7 is determined by the sign of lnr(x0). In the special case of acquisition and deletion rates being equal at x0, s(x0)=0, which implies that either there is no selection with respect to genome size or the fitness function has an extremum at x0 (Eq. 16). The genome size distribution extremum point at x0 is a maximum only if an additional condition is satisfied: for x<x0, gain and loss probabilities must satisfy P+>P− (SI Appendix, Fig. S1). This condition is met whenP+′(x0)<P−′(x0).[8]The case where x0 is a minimum of the steady-state genome size distribution is biologically irrelevant, corresponding to genome sizes tending toward either zero or infinity (SI Appendix, Fig. S1).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 3.

Different regimes for the selection and the gain/loss rates ratio (Eq. 7). In the shaded area, the genome size steady state is achieved for s(x0)>0, and accordingly, r(x0)>1. In this regime, gene loss is selected against but occurs at higher rates than gene acquisition and therefore, denoted drift. Genome sizes xmax and x1 (shown by vertical lines) denote values for which s(x0)=0 and r(x0)=1, respectively. These values are not necessarily the same, and the resulting value of x0 depends on the functional forms of s(x) and r(x), which are shown as straight blue lines for convenience.

For the calculation of the steady-state distribution of genome size, it is useful to present the population model of the genome size evolution as a random walk, where the probability of a step up is P+, and the probability of a step down is P−. The equation for the genome size distribution (Methods) has a steady-state solution:f(x)∝[P+(x)+P−(x)]−1e2∫P+(x)−P−(x)P+(x)+P−(x)dx.[9]If α depends on x, f(x) depends on the functional forms of both α(x) and β(x) and unlike the equation for x0 (Eq. 7), cannot be written using r(x) only. The genome size distribution allows one to compare different functional forms for P±(x) in terms of compatibility with observed genome sizes under the assumption that, whatever the gene gain and loss probabilities might be, they are similar in all prokaryotes. Specifically, maximizing the log likelihood of the data given a specific model allows optimization of model parameters and comparison of different cases (details are in Methods). To account for the most general case, where both the selection coefficient and the acquisition–deletion rates ratio vary with genome size, s is taken as linear in x, and r is taken as a power law:s(x)=a+b⋅x[10]andr(x)=r′⋅xλ.[11]These functional forms were chosen to minimize the number of optimized parameters. The linear selection coefficient can be regarded as a first-order expansion, and the power law functional form for r(x) was chosen, because it includes two extreme cases, those with constant and linear rates ratio, as well as all intermediates. The selection coefficient sign is not assumed a priori but is an outcome of the fitting process. Furthermore, it is in principle possible that the selection coefficient sign will be different in different ATGCs because of their different typical genome sizes.

Even for s(x) and r(x), which are approximated by low-order functions, maximizing the log likelihood requires fitting five parameters (Methods), and therefore, in principle, it is not evident that the resulting fit corresponds to a global maximum of the log likelihood rather than a local peak. We, therefore, performed the fitting with different starting points in parameter space chosen as explained in Methods. In brief, Eq. 7 is used to estimate the mean genome size for different effective population sizes, and parameters are optimized, such that goodness of fit R2 with respect to the mean genome sizes of the ATGCs and effective population sizes is maximized. This procedure resulted in five different sets of parameters (SI Appendix, Table S1) that were then used as starting points for additional optimization by log-likelihood maximization (Methods). Three of five starting points converged to similar values (SI Appendix, Table S2), whereas the remaining two converged to local maxima associated with significantly lower likelihood. In all optimized sets of parameters (from both stages), the selection coefficient is positive for all ATGCs, indicating that additional genes are, on average, beneficial as expected given the positive correlation between genome size and effective population size (SI Appendix, Fig. S2).

In all three log-likelihood optimized parameter sets, the selection coefficient is strictly positive for all ATGCs and weakly depends on x. Formally, the variation of the fixation probability with genome size is an order of magnitude smaller than the variation of acquisition and deletion rates with genome size (Methods). Accordingly, additional fittings were performed using constant selection coefficient (namely, independent of the genome size). In this case, the log-likelihood optimization always converged to the same set of parameters, with s≈6×10−12 and r values between 1.0034 and 1.0072, corresponding to smallest and largest genomes, respectively (fitted parameters are summarized in SI Appendix, Table S4). This value of s implies x0⋅s⋅Ne∼1 (i.e., the “resolution of selection” is on the order of a single gene). Notably, the dependence of r on x is weak. The mean steady-state genome size is a function of effective population size for this fit (Fig. 1).

For complementarity, fitting with constant r together with an x-dependent selection coefficient was performed as well (SI Appendix), allowing inference of the selection landscape beyond the first-order expansion for a constant acquisition–deletion rates ratio. In this case, the best fit is achieved for positive and saturating selection landscape, indicating that, on average, additional genes are beneficial but that the benefit decreases with the growth of the genome size.

As an approximation for nonuniversal, ATGC-specific factors that affect the genome size, optimization of the gene acquisition and deletion rates was performed separately for each ATGC. The selection coefficient was taken to be the same for all ATGCs and set to the value obtained in the global fitting at the previous stage. For each ATGC, r′ and λ were optimized, where the fitting of 120 parameters (compared with 5 parameters at the previous stage) is justified by the Akaike Information Criterion (AIC) (Table 1). The fitted values of r' and λ are shown in SI Appendix, Fig. S3, and the resulting genome size distributions for ATGCs with 20 or more species are shown in Fig. 4.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 4.

Comparison of the model predictions with the empirical genome size distributions. The observed genome size distributions are shown by bars for six ATGCs that consist of 20 species or more each. Genome size distributions predicted by the population evolution model are shown by red lines using the selection landscape and deletion–acquisition rate of Eqs. 10 and 11 and optimized r' and λ parameter values for each ATGC separately (optimized values are shown in SI Appendix, Fig. S3). The goodness of fit R2 is indicated for each ATGC. The ATGCs are as follows (the numbers of genomes for each ATGC are indicated in parentheses): (A) ATGC0001 (109), (B) ATGC0003 (22), (C) ATGC0004 (22), (D) ATGC0014 (31), (E) ATGC0021 (45), and (F) ATGC0050 (51). All ATGC genomes are listed in Dataset S1.

Evolution of Distinct Functional Classes of Genes.

In our model, the gene acquisition and deletion rates, α(x) and β(x), respectively, are general characteristics of the organism that do not depend on the content of the acquired or lost genetic material. In contrast, the selection coefficient inferred above represents a local average with respect to the gene content of the organism and the available genetic material in the assumed infinite gene pool. The model can be extended to account for different classes of genes that evolve under distinct selection landscapes. Specifically, the number of class i genes, xi, is determined by the stochastic equation (Methods)x˙i=ki⋅α(x)⋅F+(si(xi))−xix⋅β(x)⋅F−(si(xi)),[12]where x is the total number of genes, ki is the probability of gain of a class i gene, and the selection landscape si is assumed to be a function of xi only. The width of the steady-state distribution is determined primarily by the linear term xi/x in the loss probability. To further test the model consistency with the empirical data, steady-state distributions were calculated for subsets of genes. The subsets were chosen based on the functional classes of genes as classified in the COG (Clusters of Orthologous Genes) database (17, 18). The selection landscape and ki were optimized for the best log-likelihood fit of the distribution predicted by the model to the genomic data (SI Appendix, Table S5). The distributions were calculated using the values of α(x) and β(x) that were obtained by fitting the distributions for complete gene sets (Table 1). The mean value of xi can be approximated by the equilibrium value xi0, for which the gain and loss probabilities are equal (analogous to x0 for the complete genomes as described above):xi0e−Nesi(xi0)=kixr(x).[13]This expression can be regarded as a generalization of previously reported scaling laws for different functional classes of genes with the genome size (19⇓–21), where the ATGC-specific effective population sizes are taken into account (the full implications of this extension of the scaling analysis will be discussed elsewhere). Comparison of the empirical data and the model predictions for the number of genes in most of the functional classes shows a good fit between the model predictions and the genomic data (SI Appendix, Figs. S4 A and B and S5 A and B). For the translation system components and the genes involved in energy transformation, the log-likelihood values were −3,334 and −6,115, respectively, compared with the −6,022 value for complete genomes. Thus, notably, the genes for translation system components, the most conserved, universal functional class (22), are described by the model much better than a random subset of genes.

Finally, analogous to the whole-genome fitting procedure, to account for ATGC-specific effects, model distributions were further optimized by fitting ATGC-specific ki values. The resulting distributions for most of the functional classes showed good fits between the model and the genomic data as illustrated in SI Appendix, Fig. S6 A and B for the translation system components and genes involved in energy transformation of the largest ATGC001 (complete results are in SI Appendix, Figs. S7–S10). However, two classes of genes, namely the components of the “mobilome,” such as prophage genes and transposons, as well as the singletons (genes with no detectable orthologs within the given ATGC), dramatically deviate from model predictions (SI Appendix, Figs. S6 C and D, S9, and S10). For these gene classes, the observed distributions are significantly wider than those predicted by the model, regardless of the model parameters or the selection landscape. Accordingly, the log-likelihood values reflect the disagreement and are −18,870 and −49,384 for mobilome and singletons, respectively. Such a poor fit effectively indicates that these classes of genes evolve under evolutionary regimes qualitatively different from that of the rest of the genomes. Many of the mobilome components, in particular transposons, are prone to active duplication within the genome and therefore, cannot be described with the gene gain rate inferred for the complete gene sets. The singletons dramatically differ from the evolutionarily conserved genes shared by multiple microbes with respect to the tempo and mode of evolution. Most of the singletons encode small proteins and evolve fast, suggesting that they are associated with little (if any) benefit (23, 24). Thus, the key result of this analysis, namely the positive sign of the mean selection coefficient, does not apply to the singletons.