Evolutionary jumps in bacterial GC content

Abstract

Genomic GC (Guanine-Cytosine) content is a fundamental molecular trait linked with many key genomic features such as codon and amino acid use. Across bacteria, GC content is surprisingly diverse and has been studied for many decades; yet its evolution remains incompletely understood. Since it is difficult to observe GC content evolve on laboratory time scales, phylogenetic comparative approaches are instrumental; but this dimension is rarely studied systematically in the case of bacterial GC content. We applied phylogenetic comparative models to analyze GC content evolution in multiple bacterial groups across 2 major bacterial phyla. We find that GC content diversifies via a combination of gradual evolution and evolutionary "jumps." Surprisingly, unlike prior reports that solely focused on reductions in GC, we found a comparable number of jumps with both increased and decreased GC content. Overall, many of the identified jumps occur in lineages beyond the well-studied peculiar examples of endosymbiotic and AT-rich marine bacteria and do not support the predicted role of oxygen dependence. Our analysis of rapid and large shifts in GC content thus identifies new clades and novel contexts to further understand the ecological and evolutionary drivers of this important genomic trait.

Keywords: Lévy jumps; ecological drivers; phylogenetic models.

Fig. 1.

Summary of methods used in this study. Each major step in the analysis is numbered in the order in which it was performed. The analysis was performed independently for each of 10 order-level clades belonging to 2 bacterial phyla. Step 1: We derived the phylogenies of major bacterial clades and GC content of taxa from the genome taxonomy database (GTDB). Step 2: We obtained the ML and parameter estimates for different phylogenetic models using the phylogenies and the GC content distributions as the input. For the Brownian motion (BM) and Ornstein-Uhlenbeck (OU) models, these were obtained by exact analytical solutions implemented in the geiger package in R; while for the Lévy jumps model, these were obtained by expectation-maximization (EM) + markov chain monte carlo (MCMC) sampling method implemented in the levolution software. Step 3: For each branch in a phylogeny, we obtained the posterior probability of a jump in GC content using the phylogeny, GC contents, and the best-fit estimates of the parameters of the Lévy jumps model (obtained from step 2). These probabilities were obtained using an empirical Bayes approach implemented in the levolution software. Step 4: To calculate posterior probability thresholds to decide the presence or absence of jumps, we first simulated data with GC content jumps. The simulations were performed on the original phylogenies using the best-fit estimates of the parameters of the Lévy jumps model. We then attempted to detect the known jumps in simulated data using various posterior probability thresholds. We chose posterior probability thresholds that led to an optimal choice between precision and recall of the simulated jumps. Step 5: We deemed branches whose posterior probability of experiencing a jump (calculated from actual data in step 3) was greater than the optimal probability threshold (calculated from simulated data in step 4) as those having experienced a jump in GC content.

Fig. 2.

GC content map and location of inferred jumps in Rhizobiales. GC content was mapped onto a phylogeny of Rhizobiales using the contMap function from R package phytools. This mapping itself is only indicative of trends since it assumes a Brownian model of evolution. Branches with inferred jumps i.e. where the posterior probability of observing jump(s) is greater than the chosen threshold are indexed in filled circles. Two interesting examples of jumps in Rhizobiales are highlighted in dashed boxes, which occur in the stem branches of Liberibacter (jump index 14), an obligate plant pathogen and Bartonella (jump index 19), an obligate animal pathogen, respectively. Within the genus Bartonella, the lineage leading to B. australis experienced an upward jump (index 20). Mapping for other clades is shown in Fig. 3 and Supplementary Figs. 5–12.

Fig. 3.

GC content map and location of inferred jumps in Acetobacterales and related orders. GC content was mapped onto a phylogeny of Acetobacterales and related orders as noted in Fig 2. Branches with inferred jumps i.e. where the posterior probability of observing jump(s) is greater than the chosen threshold are indexed in filled circles. An example of a downward jump in an endosymbiont (Endolissoclinum, jump index 6) and an upward jump in an uncultivated bacterial lineage (jump index 4) are highlighted in dashed boxes. Mapping for other clades is shown in Fig. 2 and Supplementary Figs. 5–12.

Fig. 4.

Direction and magnitude of GC jumps. We estimated the magnitude of each GC content jump as the difference in median GC content of descendant taxa of a branch affected by a jump (“affected”) and descendant taxa of a sister branch not affected by a jump (“unaffected”). Jumps involving increased GC content in affected taxa were designated as upward jumps, and those involving decreased GC content of affected taxa were designated as downward jumps. a) The relation between the fraction of total jumps that were upwards and the ancestral GC content of each order-level clade. Ancestral GC content was estimated as a parameter of the Lévy jumps model using the procedure implemented in levolution. b) Distribution of jump magnitudes. Black arrows denote the median magnitudes of upward and downward jumps. c) Relation between jump magnitude and the estimated ancestral GC content, with the best-fit regression line (excluding nonendosymbiont clades). Ancestral GC content was estimated as the median GC content of unaffected taxa of the sister clade.

Fig. 5.

Proportions of upward and downward jumps across habitat and lifestyle categories. The number of upward and downward jumps are shown across 4 categories of datasets: (1) all datasets that could be analyzed for habitat or lifestyle related changes (n = 91) (2) a subset of datasets where the affected taxa (where a GC jump occurred) were host-associated, but related unaffected taxa were free-living (n = 19) (3) a subset of datasets where the affected taxa were not associated with hosts, but related unaffected taxa were host-associated (n = 12), (4) a subset of datasets where the taxa affected by the jump were isolated from marine habitats (n = 45). Numbers in the bars denote the percentage of upward and downward jumps within each category.