Abiotic selection of microbial genome size in the global ocean

Non-prokaryotic metagenomic sequences confound average genome size estimations

In this work, we employed MicrobeCensus²² for de novo estimation of the average genome size (AGS) of microorganisms captured in shotgun metagenome sequences (Fig. 1a; Supplementary Data 1). Briefly, MicrobeCensus optimally aligns metagenomic reads to a set of 30 conserved single-copy gene (CSCG) families found in prokaryotes ²². Based on these mappings, the relative abundance of each CSCG is then computed and used to estimate AGS based on the proportionality principle—that is, the AGS of the community is inversely proportional to the relative abundance of each marker genes²². Finally, a weighted average AGS is calculated that excludes outliers to obtain a robust AGS estimate for a given metagenomic sample²².

Of note, the AGS of complete prokaryotic genomes increases with the cumulative number of associated phages and other mobile genetic elements³⁷. Similarly, AGS estimates derived from metagenomic sequences of uncultured “free-living” microbes (captured in 0.1–3 µm-size filters) may also be affected by putative phage and eukaryotic microbiomes sequenced concurrently in fractionated seawater samples (see,^8,22). To evaluate this possibility in our AGS predictions, we compared AGS estimates obtained directly from quality-controlled metagenomes with estimates from the same metagenomes iteratively subjected to three (de novo) decontamination procedures to filter out potential eukaryotic and viral sequence reads (Fig. 1a; see details in the “Methods” section). Overall, putatively ‘contaminating’ viral and eukaryotic reads accounted for 1% to 20% (average 7.5%) of the high-quality trimmed sequences in the four microbial metagenome collections (Fig. 1b; Supplementary Data 1). As expected, the average proportion of contaminating sequences in metagenomes from large (0.2–3.0 µm) and small (0.1–1.2 µm) size fraction filters were the highest (~11%) and lowest (~5%), respectively (Fig. 1b). In addition, the proportion of contaminating reads was significantly dependent on the depth layer of the ocean (Kruskal-Wallis χ² = 32.40, df = 2, p < 0.001); the lowest values were in the bathypelagic (Fig. 1c).

Crucially, significant differences were observed in AGS estimates in the presence and absence of contaminating reads (repeated measures ANOVA; p < 0.0001; Supplementary Data 2), regardless of the metagenome sample, its geography, or the range of filters used (0.1–3 µm; Fig. 1d). The vast majority of metagenomes (85% of 259 samples) showed AGS estimates that were ~1–19% lower (average 5%) after the putative eukaryotic and viral metagenomic sequences were removed (Fig. 1e). These results suggest that environmental eukaryotic and viral genomic sequences affect AGS predictions for prokaryotes in marine metagenomes. Therefore, the AGS estimates reported and discussed below are based on high-quality metagenomes lacking putative viral and eukaryotic sequences (i.e., AGS4; Fig. 1), which we refer to as ‘free-living’ prokaryotic communities unless otherwise indicated.

Prokaryote genome sizes increase with ocean depth and lifestyle

In three independent microbial metagenomic collections (Supplementary Data 1) comprising temporal (ALOHA, n = 83)⁸, regional (Red Sea, n = 45)³³, and global (Malaspina Vertical Profiles, MProfile, n = 81)³⁶ ocean microbiome surveys, we found distributions of AGS estimates that were consistently unimodal (Hartigan’s dip test for unimodality, p = 0.567–0.918) in the ‘free-living’ prokaryotic communities (0.1–3 µm) sampled from the surface (3 m) to the bathypelagic ocean (~4000 m; Fig. 2a). Notably, the Malaspina Vertical Profile metagenomes (Supplementary Data 1) covered a much wider range of ocean depths (3 to 4000 m) and latitudes than previous studies estimating AGS from oceanic metagenomes, but with relatively few samples^8,13,16, providing a unique resource for our objectives. Combined analysis of all three metagenome datasets (Red Sea, ALOHA, and MProfile) revealed significant differences (Kruskal-Wallis χ² = 72.762, df = 2, p < 2.2 × 10^–16) in AGS estimates of marine prokaryotes inhabiting the three canonical ocean depth layers (Fig. 2b). Epipelagic prokaryotes (~3–200 m) had mean (± SD) and median AGS estimates of 2.21 ± 0.33 and 2.15 Mbp, respectively (n = 110). Mean (± SD) and median AGS estimates continued to increase for mesopelagic prokaryotes (> 200–1000 m) with values of 2.45 ± 0.29 and 2.43 Mbp, respectively (n = 73), and were highest for bathypelagic prokaryotes (> 1000–4000 m) with values of 2.92 ± 0.27 and 2.95 Mbp, respectively (n = 26).

The median AGS estimate range of 2.2 to ~3.0 Mbp in the sampled free-living (0.1–3 µm in size) marine prokaryotic communities (n = 209 metagenomes) is consistent with other large-scale metagenome sequence-based estimates and the sizes of metagenome-assembled prokaryotic genomes (MAGs; in 0.22–3 µm filters) from the photic ocean (surface to mesopelagic) based on the Tara Oceans Expedition (1.5–2.3 Mbp)^15,16. Overall, our metagenome sequence-based AGS estimates support the unimodal distribution of prokaryotic genome sizes recently demonstrated in environmental genomes in several biomes³⁸ and on cultured isolates (including marine bacterioplankton)^14,39. However, estimates from isolates are likely biased since current cultivation approaches tend to favor copiotrophs (see, ref. ³).

We next tested whether the derived AGS estimates depended on microbial cell size by analyzing 25 paired bathypelagic metagenomes (MDeep; Supplementary Data 1) sampled during the global Malaspina Expedition⁴⁰ in which both prokaryotic life strategies, free-living (0.2–0.8 µm size) and particle-associated (0.8–20 µm size), were sampled simultaneously³⁵. The analyzed metagenomes (MDeep) were from the Atlantic, Pacific, and Indian Ocean provinces and cover a spatial distance of 9437 km with an average depth (± SD) of 3688 ± 526 m at the tropical and subtropical latitudes (–33.55°N to 32.0788°N). These microbial metagenomes were also screened for contaminating eukaryotic and viral sequences as indicated in Fig. 1a (see details in the “Methods” section and Supplementary Data 1). The genomes of bathypelagic prokaryotes associated with marine particles (5.6 ± 0.97 Mbp) were twice as large (paired two-sided Wilcoxon test, p < 0.0001) as those of their free-living counterparts (2.65 ± 0.3 Mbp; Fig. 2c). Crucially, these AGS estimates are also consistent with read-based predictions from bathypelagic waters (4000 m) in the Pacific Ocean (Station ALOHA)¹³ and estimates based on MAGs (~3.6 Mbp) recently compiled from both cell size fractions of the same bathypelagic metagenomes³⁵.

The significant increase in AGS estimates with depth (Fig. 2b) and the twofold larger AGS for particle-associated compared to free-living bathypelagic prokaryotes (Fig. 2c) suggest larger genome size patterns in the hadal biosphere. Extending our analysis to metagenomes spanning hadal to abyssal depths (4000–10,500 m) based on seven recent Pacific Ocean metagenomes sampled from the Challenger Deep of the Mariana Trench^41,42 yielded AGS estimates in the range of 3.46–4.19 Mbp and 3.88–4.92 Mbp for free-living (0.2–3 µm) and particle-associated (> 3 µm) prokaryotes, respectively (Supplementary Data 3). These estimates are also consistent with those of MAGs reconstructed from the same metagenomes in the Challenger Deep (Mariana Trench)⁴³. Overall, this reinforces the patterns of larger AGS in particle-associated compared to free-living bathypelagic prokaryotes, and larger microbial genomes in the deep ocean compared to the upper ocean.

AGS patterns are not geographically constrained

Examination of the geographic patterns of AGS estimates showed that AGS distribution was independent of geographic distance in both the regional (Red Sea, Mantel statistic r = 0.01824, p = 0.2971) and global (MProfile, r = –0.01413, p = 0.7924) ocean metagenomes. Furthermore, AGS estimates in the vertically profiled global Malaspina metagenomes (MProfile, n = 81) were significantly independent of the Longhurst biogeochemical province sampled (n = 9; Kruskal-Wallis χ² = 1.0006, df = 8, p = 0.9982; Supplementary Data 1). The lack of covariance between the patterns of AGS estimates and geographic distance or Longhurst province sampled may reflect the high connectivity of microbial communities throughout the global ocean, particularly the redistributive effects of circulation by ocean currents and other transport processes, as well as the enormous population sizes of plankton that allow dispersal constraints to be overcome^44,45. This is consistent with the relatively small differences in microbial assemblages recently found in different ocean basins^23,46. Another possible explanation is the effect of seasonality, which can cause selection of different taxa, resulting in the succession of microbial communities and affecting their distribution (see, ref. ⁴⁷), and thus influence AGS patterns.

An assessment of the relationship between AGS and measured environmental variables (Supplementary Fig. S1; Data 1)—separately for the Red Sea metagenomes (regional scale) and Malaspina Vertical Profiles metagenomes (global scale), showed that the cumulative effect of temperature, salinity, dissolved oxygen, and depth on AGS patterns was significant at both the regional scale (n = 45; Mantel statistic r = 0.1944, p = 0.0057) and the global scale (n = 81; Mantel statistic r = 0.1779, p = 1 × 10^–4). This result suggests that environmental conditions are a driving force behind predicted AGS patterns in the marine microbiome. While no significant interaction effect was evident between many environmental variables (i.e., salinity, depth, oxygen, nitrate, and phosphate) in controlling AGS patterns (one-way ANOVA, p < 0.05; Supplementary Data 4), we found that depth and temperature covaried significantly, as expected (Spearman’s r = –0.72; p = 1.3 × 10^–14).

Further statistical tests of the relative importance of environmental factors in linear regression analyses based on variance decomposition⁴⁸ showed that temperature and depth (and nitrate in the case of the Red Sea) had similar but higher relative importance as predictors of AGS patterns than salinity or dissolved oxygen (Supplementary Fig. S2). The strong relationship between temperature and depth is consistent with evidence for a strong temperature dependence of pelagic bacterioplankton diversity and composition from single-cell genomics and metagenomics^3,23. However, our results extend this assessment from the photic ocean studied previously to the deep ocean interior, where the maximum depth of the studied samples (4000 m; Supplementary Fig. S1) reflects the average global ocean depth²⁰. This in turn captures the extensive microbiome of uncultured marine microorganisms in the dark and cold deep ocean, the largest habitat by volume in the ocean²⁰. Overall, these results suggest strong environmental selection on microbial genome sizes, with comparatively little constraint on the dispersal of the prokaryotic community of average genome size in the global ocean.

Scaling laws predict microbial genome sizes in the thermally stratified ocean

Based on the above findings, we used the gathered data at global (MProfile), regional (Red Sea), and temporal (ALOHA) scales and performed curve-fitting analysis for each data set individually and for all three together to find the best regression models that accurately explained the relationship between the AGS estimates and depth or temperature. We found that power laws best represented the relationship between AGS estimates and these two variables (Fig. 3; Supplementary Data 5). AGS estimates increased with a power of 0.072 of depth (Fig. 3a) and decreased with a power of –0.165 of temperature (Fig. 3b) on the global scale (Malaspina). Examination of these correlations for a local Red Sea data set exhibiting extreme surface temperatures (up to 32 °C; Supplementary Fig. S1; Data 1) and remarkably uniform but high temperature (~22 °C) from 200 m to the seafloor (2000–3000 m)⁴⁹, effectively reproduced the patterns observed at the global scale (Fig. 3c, d; Supplementary Data 5). Although the samples spanned a one-month expedition³³, which was not possible during the global Malaspina expedition (from December 2010 to July 2011)^34,35,36, this consistency of results increases confidence in the observed global trends and suggests that the inferred relationships are robust predictors of the underlying ecological factors. However, a steeper rate of AGS change with temperature was observed in the regional dataset (z = –0.698; n = 45; Fig. 3d) than in the global dataset (z = –0.165; n = 81; Fig. 3b), suggesting that AGS changes more rapidly with temperature extremes.

Next, we examined the effect of seasonality on our scaling predictions using a temporally resolved, localised dataset from the subtropical North Pacific ocean (Station ALOHA), where metagenomes were collected over 10 months (August 2010 to December 2011)⁸. Similar scaling relationships of AGS with depth and temperature were also found (Supplementary Fig. S3; Data 5), reinforcing patterns from the regional and global datasets (Fig. 3). However, linear regression analyses applied to individual sampling months (seasons) revealed greater variance in the direction of the relationship between AGS and environmental factors (Supplementary Fig. S4–S7), suggesting that seasonal dynamics may affect AGS patterns in marine prokaryotes. For example, correlation analyses revealed significant (Spearman r = –0.71 to 0.90; p < 0.05) but opposing relationships in the spring, summer, and winter seasons (Supplementary Fig. S4–S7). As expected, the scaling exponent of the rate of change of AGS with temperature (z = –0.07) and the regression strength (p ≤ 0.001) were much lower when all temporal metagenomes (n = 83) were analyzed together (Supplementary Fig. S3) or analyzed together with the previous datasets (Fig. 3e, f; Supplementary Data 6), compared to the separate results from the regional (Red Sea) and global (Malaspina) datasets (Fig. 3a–d; Supplementary Data 5). In contrast, the scaling exponent of the rate of change of AGS with depth remained relatively unchanged for all three comparisons (Fig. 3a, c, e). Overall, these results support the robustness of the scaling laws predicting the rate of change of AGS with depth and temperature in the global ocean microbiome (Fig. 3), despite the inferred seasonal effects.

Rates of change of AGS with respect to temperature and depth are different in photic and aphotic tropical oceans

Seawater column temperature profiles vary at different latitudes; surface water is warmer near the equator and colder at the poles. In low-latitude tropical regions, the sea surface is much warmer, resulting in a highly pronounced thermocline. In addition, the surface temperature changes relatively little seasonally in tropical regions (hence there is little seasonal change in the profiles). At high latitudes (polar regions), there is little difference between surface temperature and temperature at depth, and temperature is fairly constant (and cold) at all depths regardless of season. Thus, if the predicted scaling factor between AGS and temperature or depth is universal, we expect that (1) horizontal transects at a given depth spanning a wide range of temperatures should have the same scaling relationship with AGS in the tropic ocean and that (2) perennially cold polar marine waters should have the same linear relationship between AGS and depth.

Accordingly, we examined the relationship between AGS and depth or temperature in horizontal transects at similar depths in ocean metagenomes sampled worldwide (n = 204; Supplementary Data 1). We found that the relationship between AGS and these two variables (conditioned on the same horizontal depth) clearly separated the photic epipelagic ocean (3 to 200 m; n = 8 depths; Fig. 4a) from the aphotic ocean (mesopelagic and bathypelagic; n = 10 depths; Fig. 4b). In both cases, the dependence of AGS on depth (Fig. 4c) and temperature (Fig. 4d) was significantly linear (p < 0.05; Fig. 4c, e) only in the photic layer, but not in the aphotic layer (Fig. 4d, f). This again suggests that the genome size landscape of uncultured marine prokaryotes in the global ocean is clearly separated along the dimensions of the two abiotic factors. However, the estimated regression slope of AGS with depth (mean ± SD: 0.0038 ± 0.00049; Fig. 4c) was lower than that for temperature (mean ± SD: −0.0591 ± 0.0264; Fig. 4e), suggesting that the rate of change of AGS with temperature is 16-fold higher than with depth. This corresponds to a change of approximately 4 and 59 genes per unit of depth and temperature, respectively, assuming an average-sized microbial protein-coding gene of 1000 base pairs in planktonic prokaryotes found at similar depths up to 200 m. The observed distinct pattern of change in AGS with these environmental variables in photic and aphotic oceans likely reflects the pronounced thermocline conditions in the photic layer at low latitudes.

Next, we examined the relationship between AGS and both variables in polar metagenomes (n = 94) extending from the surface to 3800 m depth (Supplementary Data 1) and covering different transects in the Arctic and Antarctic Oceans⁵⁰ and the Tara Ocean Polar Circle Expedition⁵¹. Unlike tropical waters at low latitudes, polar waters do not have a thermocline and temperatures do not exhibit pronounced seasonal variations. This raises the question of whether the same inverse relationship between AGS and temperature exists in permanently cold polar marine waters. However, we found no statistically significant linear correlation (p < 0.05) between AGS and depth or temperature (Supplementary Fig. S8), suggesting that AGS varies invariably in permanently cold waters regardless of seawater depth. To gain further insight, we compared AGS estimates between metagenomes from the tropical ocean (n = 234) and those from high-latitude polar waters (n = 98), both extending from the surface to the bathypelagic ocean. The results show that polar waters contain free-living prokaryotes (0.1–3 µm) with significantly larger genomes (unpaired two-sided Wilcoxon test, p = 0.0008) than their counterparts from low-latitude tropical oceans (mean ± SD: 2.73 ± 0.72 Mbp vs. 2.41 ± 0.39 Mbp; Fig. 5). These results reinforce our earlier finding that prokaryotic genomes from low-latitude tropical marine waters are larger in the cold interior deep ocean than in the photic ocean (see above), and they corroborate reports of large genomes in metagenome-resolved genomes from the polar Arctic⁵¹. However, the findings leave unanswered the question of whether interior ocean prokaryotes harboring large genomes have a greater metabolic capacity due to additional genes in the genome.

Gene lengths and family sizes correlate strongly with microbial genome sizes

If the smaller AGS estimate in prokaryotic communities occupying the photic ocean are due to fewer gene duplications and shorter coding genes in the streamlined bacteria of the surface ocean^3,8, then we expect that the potential for gene family expansion (duplications) and elongation of coding genes predominates in microorganisms with larger genomes colonizing the nutrient-rich interior of the ocean. Therefore, we examined the relationship between AGS and coding gene traits, including the average gene length (AGL), percent GC content, and gene family size, using full-length coding genes previously assembled for global Malaspina Vertical Profiles metagenomes (i.e., MProfile)³⁴. We defined gene families as the number of unique protein-coding genes distinct to each metagenome regardless of water column depth, based on a catalog of 32.7 million full-length non-redundant genes reconstructed from the same metagenomes³⁴. The gene catalog contains microbial coding gene sequences (effectively protein families) clustered at 95% global sequence identity and 80% overlap across the shorter gene. The processed matrix of gene family copies across all MProfile metagenomes (available here, https://doi.org/10.6084/m9.figshare.19673688.v1) was used to compile a binary table of gene presence and absence and the corresponding counts of unique genes per metagenome (Supplementary Data 1).

Both AGL (Fig. 6a) and GC content (Fig. 6b) increased with ocean depth, whereas the number of unique genes per unit coding sequence length (Fig. 6c) or sequenced metagenomic effort (Supplementary Fig. S9) decreased with depth. The median AGL in the three traditional ocean depth layers were 583 bp (epipelagic), 612 bp (mesopelagic), and 637 bp (bathypelagic), indicating a relative length increase of 54 bp from the upper to bathypelagic ocean (Fig. 6a). At the same time, the number of unique genes per coding sequence length (UGPCL) decreased by a factor of two from the epipelagic to bathypelagic ocean (Fig. 6c). These results are consistent with reports of microbial gene length shortening in the low-nitrogen, sunlit ocean^3,8 and gene duplications observed in microbial species with larger genomes, including representative marine prokaryotes^52,53, and in environmental genomes of deep-sea microbes¹³.

Spearman rank correlation analyses also showed that the AGL of coding genes was significantly positively correlated with changes in AGS, GC content, and water column depth (r = 0.48–0.69, p = 10⁻⁴–10⁻¹⁰; Fig. 6d), but significantly negatively correlated with UGPCL and temperature (r = −0.49 to −0.61, p = 10⁻⁶–10⁻⁹; Fig. 6d). In contrast, UGPCL was significantly positively correlated with temperature (r = 0.73, p = 3.1 × 10^–21), but negatively correlated with depth (r = −0.72, p = 4.8 × 10⁻⁸), AGL (r = −0.61, p = 3.3 × 10⁻⁹), GC content (r = −0.39, p = 0.0013), and AGS (r = −0.71, p = 1.2 × 10⁻¹³). Collectively, these results suggest that streamlining selection favors shorter coding genes in the photic ocean and that the size patterns of microbial genomes in the interior ocean reflect the elongation of coding genes and expansion of gene family size, which likely leads to genetic redundancy.

Source: Ecology - nature.com