Strong phylogenetic signals in both GC contents and growth temperatures
A significant force shaping prokaryotic evolution is horizontal gene transfer, making the genealogical relationships among bacteria and archaea exhibit a somewhat network-like structure. If bifurcation is not the phylogeny's dominant pattern, most phylogenetic comparative methods will not be necessary for prokaryotic evolutionary studies. We are unsure how much this impression has influenced the researchers in prokaryotic genomic studies, but many papers did not use any phylogenetic comparative methods. Despite the frequent horizontal gene transfers, careful examination of the prokaryotic phylogeny could see a statistical tree (41–43). In principle, the necessity of phylogenetic comparative methods depends on the significance of the phylogenetic signal, a measure of the correlation between the evolution of the analyzed trait and the presumed phylogenetic tree. We first measured the phylogenetic signals of the analyzed traits for the 681 bacteria and 155 archaea obtained from the database TEMPURA (36). As shown in Table 1, all the λ values are close to one, which indicates that simple statistical analysis that does not account for common ancestry's effect would lead to inaccurate results (44, 45).
Table 1
The phylogenetic signals of the variables analyzed in this study.
|
Bacteria
|
|
|
Archaea
|
|
|
Traits
|
n
|
Pagel′s λ
|
P
|
n
|
Pagel′s λ
|
P
|
Tmax
|
681
|
0.957
|
3.5×10−178
|
155
|
1.000
|
1.5×10−72
|
Topt,
|
681
|
0.950
|
1.1×10−196
|
155
|
0.988
|
5.7×10−70
|
Tmin
|
681
|
0.933
|
6.6×10−152
|
155
|
0.966
|
1.9×10−53
|
GCw
|
681
|
1.000
|
2.6×10−294
|
155
|
1.000
|
5.8×10−60
|
GCp
|
681
|
1.000
|
4.7×10−292
|
155
|
1.000
|
2.4×10−59
|
GC4
|
681
|
1.000
|
4.8×10−238
|
155
|
1.000
|
5.1×10−53
|
GCnon
|
681
|
1.000
|
8.6×10−305
|
155
|
1.000
|
2.0×10−65
|
GCtRNA
|
681
|
0.998
|
6.0×10−275
|
155
|
1.000
|
2.1×10−91
|
GC5S
|
646
|
1.000
|
7.0×10−178
|
130
|
1.000
|
9.1×10−51
|
GC16S
|
681
|
0.999
|
8.7×10−250
|
155
|
0.996
|
3.7×10−86
|
GC23S
|
681
|
1.000
|
2.1×10−245
|
155
|
1.000
|
6.8×10−83
|
Tmax, Topt, and Tmin represent maximal, optimal, and minimal growth temperature, respectively; GCw, GCp, GC4, GCtRNA, GC5S, GC16S, GC23S, and GCnon represent the GC contents of the whole genome, the protein-coding sequences, the fourfold degenerate sites, the genes coding tRNAs, the genes coding 5S rRNA, the genes coding 16S rRNA, the genes coding 23S rRNA, and the non-coding DNA (including intergenic sequences and untranslated regions of mRNA), respectively. The phylogenetic signals of the chromosomal genes, the plasmid genes, the core genes, and the accessory genes are also very close to one and deposited in Additional file 1: Tables S1-S4. |
Bacterial but not archaeal genomic GC contents correlated with growth temperatures
We used the phylogenetic generalized least squares (PGLS) regression to examine the relationships between GC contents and growth temperatures. The significant positive and negative slopes of the regressions correspond to significant positive and negative correlations, respectively. The value of the slope represents the phylogenetically corrected rate of change in GC content as growth temperature changes. Four phylogenetic models, the Brownian motion model (BM), the Ornstein-Uhlenbeck model with an ancestral state to be estimated at the root, the Pagel's lambda model, and the early burst model, have been applied in the analysis. Their results are qualitatively identical and quantitatively similar. As the four models lead to the same conclusion, the trivial differences among their results are unrelated to understanding the relationship between GC content and growth temperature. We present the BM model results in the main text and deposit other models' results as Additional file 1.
Consistent with numerous previous studies, we found positive correlations between the GC contents of structural RNA genes and growth temperature in bacteria and archaea (Table 2). We noticed a rank in the slope values, from Tmax, Topt, to Tmin.
Table 2
PGLS regression of GC contents and growth temperatures.
|
Bacteria
|
|
|
Archaea
|
|
|
|
Slope
|
P
|
PBH
|
Slope
|
P
|
PBH
|
GCw-Tmax
|
7.1×10−4
|
7.1×10−4
|
9.4×10−4
|
6.6×10−4
|
0.115
|
0.153
|
GCw-Topt
|
5.7×10−4
|
0.009
|
0.011
|
3.3×10−4
|
0.377
|
0.503
|
GCw-Tmin
|
2.8×10−4
|
0.156
|
0.226
|
5.2×10−4
|
0.126
|
0.168
|
GCp-Tmax
|
6.6×10−4
|
0.002
|
0.002
|
5.6×10−4
|
0.183
|
0.209
|
GCp-Topt
|
5.3×10−4
|
0.015
|
0.016
|
2.4×10−4
|
0.522
|
0.597
|
GCp-Tmin
|
2.5×10−4
|
0.202
|
0.231
|
4.6×10−4
|
0.180
|
0.205
|
GC4-Tmax
|
0.001
|
0.003
|
0.003
|
9.9×10−4
|
0.321
|
0.321
|
GC4-Topt
|
0.001
|
0.016
|
0.016
|
2.2×10−4
|
0.806
|
0.806
|
GC4-Tmin
|
5.5×10−4
|
0.239
|
0.239
|
6.9×10−4
|
0.393
|
0.393
|
GCnon-Tmax
|
8.0×10−4
|
1.7×10−4
|
2.8×10−4
|
9.1×10−4
|
0.025
|
0.041
|
GCnon-Topt
|
6.4×10−4
|
0.004
|
0.006
|
6.4×10−4
|
0.080
|
0.129
|
GCnon-Tmin
|
2.7×10−4
|
0.170
|
0.226
|
6.5×10−4
|
0.048
|
0.077
|
GCtRNA-Tmax
|
4.1×10−4
|
2.2×10−16
|
5.9×10−16
|
7.1×10−4
|
1.8×10−11
|
7.2×10−11
|
GCtRNA-Topt
|
3.9×10−4
|
2.6×10−14
|
6.9×10−14
|
5.0×10−4
|
2.5×10−7
|
6.7×10−7
|
GCtRNA-Tmin
|
1.5×10−4
|
9.1×10−4
|
0.002
|
4.2×10−4
|
1.8×10−6
|
4.7×10−6
|
GC5S-Tmax
|
5.5×10−4
|
1.2×10−6
|
2.4×10−6
|
0.001
|
1.9×10−5
|
3.9×10−5
|
GC5S-Topt
|
4.4×10−4
|
1.4×10−4
|
2.9×10−4
|
8.9×10−4
|
1.6×10−4
|
3.2×10−4
|
GC5S-Tmin
|
3.5×10−4
|
0.001
|
0.002
|
6.1×10−4
|
0.005
|
0.010
|
GC16S-Tmax
|
5.4×10−4
|
2.2×10−16
|
5.9×10−16
|
8.2×10−4
|
3.9×10−11
|
1.0×10−10
|
GC16S-Topt
|
5.2×10−4
|
2.2×10−16
|
8.8×10−16
|
7.2×10−4
|
1.1×10−10
|
4.5×10−10
|
GC16S-Tmin
|
4.6×10−4
|
2.2×10−16
|
8.8×10−16
|
5.5×10−4
|
8.5×10−8
|
3.4×10−7
|
GC23S-Tmax
|
6.6×10−4
|
2.2×10−16
|
5.9×10−16
|
0.001
|
2.2×10−16
|
1.8×10−15
|
GC23S-Topt
|
6.5×10−4
|
2.2×10−16
|
8.8×10−16
|
0.001
|
1.2×10−14
|
9.5×10−14
|
GC23S-Tmin
|
4.9×10−4
|
2.2×10−16
|
8.8×10−16
|
8.3×10−4
|
8.0×10−11
|
6.4×10−10
|
GC contents were the dependent variables, and growth temperatures were the independent variables. The results in this table were obtained using the Brownian motion model. Similar results obtained from three other models are deposited in Additional file 1: Table S5-S7. PBH, Benjamini-Hochberg adjusted P value. Please see Table 1 for the meanings of the other abbreviations. |
Interestingly, we also found positive correlations of Tmax and Topt with various indexes of genomic GC contents, GCw, GCp, GC4, and GCnon, in bacteria (Table 2). Nevertheless, bacterial Tmin is not correlated with three of the four GC content indexes (Table 2). In archaea, none of the three temperature indexes (Tmax, Topt, or Tmin) have any significant correlations with any of the four genomic GC content indexes (Table 2).
If growth temperature could shape GC contents by the stabilities of RNA secondary structures and DNA double helix, a structural RNA or a DNA double helix that is stable at the Tmax or Topt is, of course, stable at the Tmin. This logic makes it reasonable to see that the Tmin has weaker or no significant correlations with GC contents.
The difference in the correlations between bacteria and archaea might be attributed to either unknown intrinsic differences between these two domains or the substantial difference in the sample size, 681 vs. 155.
Sample sizes matter
If the lack of significant correlations between genomic DNA and Tmax and Topt in archaea results from the small sample size, the correlations in bacteria will be lost when the sample size of bacteria is reduced to 155. For this reason, we randomly selected 155 bacteria from the 681 bacterial samples for 1000 rounds. The results of the resampling analysis confirmed the idea, the sample sizes matter (Table 3). In > 950 rounds, the genomic GC content indexes (GCw, GCp, GC4, and GCnon) are not correlated with Tmax or Topt (P > 0.05). This result could also explain the difference between the present study with Hurst and Merchant (2001), which did not find significant correlations between GCw/GC3 and Topt by phylogenetic analysis of about 100 prokaryotes. Meanwhile, a few positive correlation cases happen, indicating that significant positive correlations could also be found by chance when the analyzed sample is small.
Table 3
The appearance of correlations in 1000 rounds of resampling analyses.
|
Significantly Negative
(P < 0.05)
|
Not Significant
(P > 0.05)
|
Significantly Positive
(P < 0.05)
|
GCw-Tmax
|
0
|
974
|
26
|
GCw-Topt
|
0
|
991
|
9
|
GCp-Tmax
|
0
|
976
|
24
|
GCp-Topt
|
0
|
993
|
7
|
GC4-Tmax
|
0
|
962
|
38
|
GC4-Topt
|
0
|
992
|
8
|
GCnon-Tmax
|
0
|
974
|
26
|
GCnon-Topt
|
0
|
992
|
8
|
GCtRNA-Tmax
|
0
|
12
|
988
|
GCtRNA-Topt
|
0
|
21
|
979
|
GC5S-Tmax
|
0
|
308
|
692
|
GC5S-Topt
|
0
|
473
|
527
|
GC16S-Tmax
|
0
|
0
|
1000
|
GC16S-Topt
|
0
|
0
|
1000
|
GC23S-Tmax
|
0
|
0
|
1000
|
GC23S-Topt
|
0
|
0
|
1000
|
In each round of resampling, 155 samples were randomly drawn from the 681 bacteria. PGLS regression analyses were performed for each round. GC contents were the dependent variables, and growth temperatures were the independent variables. The results in this table were obtained using the Brownian motion model. Please see Table 1 for the meanings of the other abbreviations. The datasets for each round of resampling are deposited in Additional file 2: Data S1. |
Besides, the correlations between growth temperature and the GC contents of structural RNA genes might also be lost occasionally when the sample size is severely reduced (Table 3). In the 1000 rounds of resampling, lacking significant correlations happens in 308 (for Tmax) and 473 (for Topt) rounds for 5S rRNA genes, and 12 (for Tmax) and 21 (for Topt) rounds for tRNA genes. However, in the 16S and 23S rRNA genes, positive correlations were consistently observed in all the 1000 rounds of resampling. We suspected that the tens of times more nucleotides in 16S and 23S rRNA than 5S rRNA and tRNAs make the results of 16S and 23S rRNAs less sensitive to small sample sizes.
In statistics, the rule of thumb boundary between small and large samples is n = 30. However, the results in Table 3 indicate that n = 155 is a too-small sample in the phylogenetic comparative analyses of the relationship between growth temperature and genomic GC content. Because of the common ancestor, two closely related lineages with highly similar growth temperatures and GC contents should be regarded as nearly one effective sample rather than two independent samples. The effective sample size in phylogenetic comparative studies should be much lower than the census number of the analyzed lineages.
Positive correlations observed in genes of both chromosomes and plasmids
Previous studies showed that plasmids have significantly lower GC contents than chromosomes (8, 46, 47). Therefore, we examined the correlations between growth temperatures and GC contents separately in chromosomes and plasmids. The separations of plasmids and chromosomes are arbitrary. We strictly followed the classifications of chromosomes and plasmids of the NCBI genome database (48). Among the 681 bacteria and 155 archaea analyzed above, 172 bacteria and 42 archaea have plasmid genomes. The bacterial chromosomes also have GC contents (GCw, GCp, GC4, and GCnon) positively correlated with Tmax and Topt (Table 4). Interestingly, the same pattern was also found in the bacterial plasmids (Table 4) in spite that the correlations of Tmax with GC4 and GCnon are just significant at marginal levels (0.05 < P < 0.1). All these correlations are not significant in archaea.
Table 4
PGLS regression of GC contents and growth temperatures in chromosomes and plasmids.
|
Plasmid
|
|
|
Chromosome
|
|
|
|
Slope
|
P
|
PBH
|
Slope
|
P
|
PBH
|
GCw-Tmax
|
0.001
|
0.009
|
0.043
|
9.6×10−4
|
0.029
|
0.043
|
GCw-Topt
|
0.001
|
0.005
|
0.031
|
9.6×10−4
|
0.023
|
0.031
|
GCp-Tmax
|
0.001
|
0.016
|
0.043
|
9.1×10−4
|
0.038
|
0.046
|
GCp-Topt
|
0.001
|
0.010
|
0.031
|
9.2×10−4
|
0.031
|
0.034
|
GC4-Tmax
|
0.002
|
0.072
|
0.072
|
0.002
|
0.027
|
0.043
|
GC4-Topt
|
0.002
|
0.044
|
0.044
|
0.002
|
0.017
|
0.031
|
GCnon-Tmax
|
8.3×10−4
|
0.055
|
0.060
|
0.001
|
0.021
|
0.043
|
GCnon-Topt
|
9.3×10−4
|
0.025
|
0.031
|
0.001
|
0.021
|
0.031
|
GC contents were the dependent variables, and growth temperatures were the independent variables. The results in this table were obtained using the Brownian motion model. Similar results obtained from three other models are deposited in Additional file 1: Table S9-S11. PBH, Benjamini-Hochberg adjusted P value. Please see Table 1 for the meanings of the other abbreviations. |
In the two previous studies comparing the GC content between plasmids and chromosomes (46, 47), the common ancestor effect was not accounted for. By the way, we performed a phylogenetic paired t-test (49) and confirmed the pattern of lower GC content in plasmids (Table S8).
Positive correlations were observed in both core genes and accessory genes
To correspond to the previous gene-centered studies (35), we examined the correlations in bacterial core genes, i.e., genes present in all the bacteria. With the increase in the number of bacteria, the number of core genes decreases rapidly. With a trade-off between the number of core genes and the number of bacteria, we selected 28 core genes present in 420 genomes, mostly ribosomal protein genes. Significant positive correlations have been found between GC contents (GCp and GC4) and growth temperatures, Tmax, and Topt (Table 5).
Table 5
PGLS analysis of GC contents and growth temperatures in core genes and accessory genes
|
Core Genes
|
Accessory Genes
|
|
Slope
|
P
|
PBH
|
Slope
|
P
|
PBH
|
GCp-Tmax
|
7.6×10−4
|
9.6×10−4
|
0.002
|
9.0×10−4
|
0.001
|
0.002
|
GCp-Topt
|
6.4×10−4
|
0.007
|
0.025
|
6.3×10−4
|
0.026
|
0.030
|
GC4-Tmax
|
0.002
|
6.3×10−4
|
0.002
|
0.002
|
0.003
|
0.003
|
GC4-Topt
|
0.002
|
0.004
|
0.025
|
0.002
|
0.019
|
0.030
|
GC contents were the dependent variables, and growth temperatures were the independent variables. The results in this table were obtained using the Brownian motion model. Similar results obtained from three other models are deposited in Additional file 1: Table S12-S14. PBH, Benjamini-Hochberg adjusted P value. Please see Table 1 for the meanings of the other abbreviations. |
At the opposite side of the core genes, the accessory genes are present in one or a few bacteria. When we define the accessory genes as the genes present in less than 5% of the analyzed bacterial genomes, on average, each bacterium has 152 accessory genes. Positive correlations were observed between GC contents (GCp and GC4) and growth temperatures (Tmax and Topt), although the values of significance are slightly larger than those in core genes (Table 5). Similar patterns were observed when we increased the threshold in defining accessory genes to 10% (P < 0.05 for all cases).
By the way, we compared the GC content between bacterial core genes and accessory genes using a phylogenetic paired t-test (49). Unlike the previous analysis of 36 prokaryotes that did not account for the effect of common ancestors (50), we did not observe significant differences in GC content between the core genes and the accessory genes (Table S15). We also compared the chromosomal accessory genes and plasmid accessory genes. The accessory genes on chromosomes have significantly higher GC contents than those on plasmids (Table S16).
Qualitative data on growth temperature lead to the same conclusion
In the ProTraits database and the IMG database (51, 52), many prokaryotes lack quantitative measures of growth temperature but are qualitatively classified into four categories: psychrophiles/psychrotrophiles, mesophiles, thermophiles, and hyperthermophiles. We constructed a qualitative dataset of prokaryote growth temperature, including data downloaded from these two datasets and the prokaryotes in the TEMPURA database classified into the four categories referring (36). By assigning 1, 2, 3, and 4 to the psychrophiles/psychrotrophiles, mesophiles, thermophiles, and hyperthermophiles, respectively, we transformed the qualitative characters into numerical values. Because only some of the genomes have been completely sequenced, we used their GCw values downloaded directly from the NCBI genome database. PGLS regression revealed a positive correlation between GCw content and growth temperature in bacteria (slope = 0.457, P = 0.001), but not in archaea (slope = −0.582, P = 0.170). Although this dataset (4696 bacteria and 279 archaea) is much larger than analyzed above (681 bacteria and 155 archaea), it lost much information during the qualitative classification. All the differences in growth temperature within each category disappear.
We also examined whether the contrast in the temperature category is correlated with the contrast in the GC content between terminal tips of the phylogenetic tree. Referring (6), 273 pairs of bacteria and 41 pairs of archaea were retrieved from the Genome Taxonomy Database (53). On average, the bacteria with higher ranks in Topt have 1.43% more GC than their paired bacteria with lower ranks. Pairwise comparison showed significantly higher GC contents in the bacteria with higher ranks in growth temperature (Wilcoxon signed rank test, P = 0.019, Fig. 1A). Still, no significant differences were observed between paired archaea with different growth temperature ranks (Wilcoxon signed rank test, P = 0.446, Fig. 1B).
Evolutionary jumps in GC contents are correlated with growth temperature changes
Recently, Mahajan and Agasheand (3) found that the Lévy jumps model (54) could better explain prokaryotic GC content evolution than the Brownian model. The GC content evolves at a constant rate and sometimes experiences discrete changes, i.e., jumps. Following Mahajan and Agasheand (3), we first confirmed that the Lévy jumps model could better explain the GCw and Topt in our dataset than the simple Brownian model.
As the Lévy jumps model has not been integrated into the PGLS packages, we cannot replace the BM model with this one. As an alternate, we retrieved the detected jumps in GCw and examined whether significant changes in Topt accompany them. The phylogenetic locations of jumps were inferred using the levolution software (54). In this procedure, only the posterior probabilities (pp) of the presence of > 0 jumps were estimated, but the exact number or magnitude of jumps on each branch could not be predicted. In practice, the “precision” of jump inference is negatively correlated with the “recall” of actual jumps. By adjusting the threshold of posterior probabilities of the presence of > 0 jumps for a precision > 85%, we obtained the GCw jumps with 88.5% precision and an acceptable recall of 37.0%. As shown in Fig. 2A, the magnitudes of the GCw jumps are positively correlated with the changes in Topt (Spearman′s rank correlation, 2-tailed, n = 108, rho = 0.209, P = 0.030). When the precision of jump inference was increased to 96.9%, the recall decreased to 21.5%, and no significant correlation was observed in the smaller sample (Spearman′s rank correlation, 2-tailed, n = 56, rho = 0.195, P = 0.150).
Meanwhile, we detected the evolutionary jumps in Topt using the same model. By adjusting the threshold of posterior probabilities, we inferred the Topt jumps with 95.3% precision and 21.5% recall. A positive correlation was observed between the magnitudes of the jumps in Topt and the changes of GC contents at the positions of Topt jumps (Spearman's rank correlation, 2-tailed, n = 86, rho = 0.280, P = 0.009, Fig. 2B).