Parallel evolution of genome streamlining in Marinimicrobia.

We found that the habitat preferences observed throughout the marinimicrobial tree were strongly correlated with patterns of genomic organization. Genomes of epipelagic Marinimicrobia exhibited signatures of streamlining such as lower percent GC content and shorter intergenic regions (16, 36) (Fig. 1a and ​andcc and ​and2b).2b). Moreover, epipelagic Marinimicrobia also exhibited fewer nitrogen atoms per residue side chain (N-ARSC) in their encoded proteins, consistent with the hypothesis that this represents an adaptation to reduce nitrogen demand in oligotrophic surface waters (32, 37). In contrast, mesopelagic Marinimicrobia contained lower carbon content in their encoded proteins (C-ARSC), consistent with higher nitrogen but lower carbon availabilities in deeper waters (32) (Fig. 2b). Many of these features were correlated, suggesting the presence of distinct genomic modalities in epipelagic versus mesopelagic Marinimicrobia (see Fig. S3 at figshare.com/projects/Marinimicrobia_Pangenomics/30881), which is consistent with observations of a genomic transition zone between these two regions (32).

Overall, we found GC content, mean intergenic spacer length, N-ARSC, and C-ARSC to be significantly different between all Marinimicrobia in the two habitat categories (Mann-Whitney U test, P <0.01, “All” category in Fig. 2a), and our intraclade comparisons demonstrate that these disparities evolved independently in several different clades (Fig. 2b). Percent GC content was the most prominent feature that shifted with habitat preference, with clades 2, 5, 6, and 7 all showing significantly lower values in epipelagic versus mesopelagic genomes. C-ARSC was the next most prevalent feature distinguishing between habitat groups, with 3 clades showing significant differences. Interestingly, clade 1 did not show genome features that were significantly different between groups despite the epipelagic genomes in this group displaying marked indications of streamlining (Fig. 1a and ​andb).b). This is likely because only 5 genomes in this clade are more abundant in mesopelagic waters, which limits the statistical power of comparisons. Moreover, the genome with the highest GC content, second highest N-ARSC, and lowest C-ARSC in clade 1 belongs to Marinimicrobia NORP180, which is the genome in this clade found to be most abundant in the mesopelagic habitat (Fig. 1), suggesting that genomic transitions have begun to evolve in this lineage. Other lineages that may have recently switched between habitats may have not had enough time to acquire the genomic features typical of their new environment, indicating that these traits require long periods of time to evolve.

Although C-ARSC and N-ARSC are not typically considered indicators of genome streamlining, our findings indicating that these metrics vary consistently with other aspects of streamlined genomes in Marinimicrobia suggest that they represent salient features that future studies should consider when assessing bacterioplankton genome evolution. Recent analysis of whole-community genomic differences between epipelagic and mesopelagic microbes has also shown that C-ARSC and N-ARSC are strongly connected to other features associated with streamlining, such as percent GC content and mean intergenic region length (32). The evolution of efficient nutrient utilization strategies is an important aspect of streamlining theory (16), and because proteins comprise a large pool of cellular carbon and nitrogen, it is likely that shifts in N-ARSC or C-ARSC that aid in the efficient allocation of macronutrients are highly advantageous in oligotrophic environments. Recent modeling of genome evolution in marine bacteria has indicated that changes in nutrient allocation can exert a strong influence on other genomic features such as percent GC content (38).

Interestingly, we did not observe a consistent reduction in genome size in epipelagic versus mesopelagic Marinimicrobia (Fig. 2b), which is perhaps paradoxical, considering that the genomes in the former group appear more streamlined in most other aspects. We would expect that the genomes of epipelagic Marinimicrobia would experience at least a modest decrease in genome size due to their shorter intergenic regions, but this reduction may be minor and not statistically significant given that coding regions comprise the vast majority of total DNA. The lack of large differences in genome size between epipelagic and mesopelagic Marinimicrobia is not entirely surprising given that streamlining has been hypothesized to take place over a range of genome sizes, since adaptation to a given environment requires particular coding potential that would in turn dictate genome size (16). Our results are consistent with this hypothesis and suggest that in the genome streamlining that we observed in the Marinimicrobia, which appears to be driven largely by differential selection in distinct environments and nutrient regimes, we would not necessarily expect to observe large differences in genome size but rather changes in features such as intergenic spacer length, percent GC content, and N-ARSC and C-ARSC, for which we saw strong and repeated shifts.

Convergence of functional repertoires in Marinimicrobia.

We also identified clear differences in genomic repertoires between epipelagic and mesopelagic Marinimicrobia, with our pangenomic analyses revealing 758 orthologous groups that were enriched in either of the habitat groups (Fisher’s exact test; corrected P values of <0.01). Many epipelagic Marinimicrobia in clades 1 to 7 have acquired proteorhodopsin proton pumps, photolyases associated with UV stress, peroxide stress genes, and phosphate starvation genes, consistent with convergence toward similar mechanisms for life in oligotrophic surface waters in which UV radiation, peroxides, and low nutrient levels are prevalent stressors (39) (Fig. 3). Our phylogenetic analysis of marinimicrobial photolyases and proteorhodopsins indicated that these genes were acquired independently in different clades, supporting the hypothesis of parallel evolution of distantly related Marinimicrobia species toward similar ecological niches (Fig. 4a and ​andb).b). In contrast, several mesopelagic Marinimicrobia genomes have independently acquired the cellular machinery for nitrate respiration (Fig. 1d and ​and4c4c and ​andd),d), mirroring the independent gene acquisitions observed in epipelagic groups and consistent with findings that many Marinimicrobia are poised to exploit alternative electron acceptors under conditions of low oxygen concentrations (21, 41).

Open in a separate window

FIG 3

Presence of selected marker genes and bioenergetic complexes across the Marinimicrobia. The dendrogram on top shows the habitat-based genome clustering, and the color strips below it show the habitat groups and clades of the genomes. The colors used to denote habitat groups and clades are identical to those in Fig. 1 and ​and2.2. Unique identifiers for the protein clusters are indicated on the right.

Open in a separate window

FIG 4

Evolutionary history of metabolic marker genes in the Marinimicrobia. Phylogeny of marinimicrobial photolyases (COG0415; panel a), proteorhodopsins (ENOG4111G9N; panel b), NarG (COG5013; panel c), NarH (COG1140; panel d), NuoH (COG1005; panel e), and NqrA (COG1726; panel f). Each phylogeny also contains reference sequences, which were obtained from the EggNOG website for all phylogenies except those of proteorhodopsins, which were obtained from the MicRhoDE database (40) (see Materials and Methods). Solid circles denote nodes with support values of >0.95. Interactive phylogenies are available online at http://itol.embl.de/shared/faylward.

We also identified several cytochrome-associated proteins and pyruvate:ferredoxin/flavodoxin oxidoreductases (PFORs) that were differentially enriched in epipelagic versus mesopelagic Marinimicrobia (Fig. 3; see also Table S1 at figshare.com/projects/Marinimicrobia_Pangenomics/30881), with all PFOR subunits and most cytochrome subunits more prevalent in mesopelagic groups. Recent work has shown that cytochrome c oxidases are coexpressed with anaerobic respiratory genes in some Marinimicrobia under conditions of low levels of dissolved oxygen, suggesting that these cytochromes are either involved in the coreduction of electron acceptors other than oxygen or involved in the rapid switching between aerobic and anaerobic metabolism (21). Another study of microbial communities in a subseafloor aquifer that included Marinimicrobia also identified genomic signatures of anaerobic respiration despite oxic conditions (41), further suggesting that switching between electron acceptors may be a dynamic process in deep marine environments that is dictated by prevailing environmental conditions. Overall, the presence of a wide array of cytochromes in mesopelagic Marinimicrobia is consistent with their use of a variety of terminal electron acceptors, which is similar to what has been observed in other well-studied microbes such as Shewanella oneidensis (42). The prevalence of PFORs in mesopelagic Marinimicrobia is potentially consistent with the metabolic versatility of these bacteria, since the ability to shuttle electrons through alternative carriers such as ferredoxin or flavodoxin may allow a broader range of respiratory complexes to be used. Our phylogenetic analyses of PFORs are consistent with multiple independent acquisitions by mesopelagic Marinimicrobia (see Fig. S4 at figshare.com/projects/Marinimicrobia_Pangenomics/30881), further indicating convergence toward similar bioenergetic modalities among habitat groups.

Perhaps most strikingly, the different evolutionary forces experienced by epipelagic and mesopelagic Marinimicrobia also appear to have altered their cellular bioenergetics, as we observed a prevalence of NQR-(Na⁺) respiratory complexes in epipelagic Marinimicrobia, while mesopelagic groups appear to have largely replaced this with a canonical NDH-(H⁺) respiratory complex (cNDH) (Fig. 1d and ​and3).3). Most genomes in both groups also encoded a noncanonical NDH-(H⁺) respiratory complex (nNDH) for which the NADH reductase subunits were missing, suggesting that alternative electron donors such as flavodoxin may be used, which is similar to what has been observed in other groups (43). These findings indicate that while use of both a sodium motive force and a proton motive force is prevalent across Marinimicrobia, the relative levels of importance of these bioenergetic gradients and how they are used differ between groups. There are a number of possible explanations for these differences between epipelagic and mesopelagic Marinimicrobia. First, the canonical NDH-(H⁺) complex is likely more efficient than the NQR-(Na⁺) pump (44), which is potentially favorable to mesopelagic Marinimicrobia since carbon and energy are less readily available deeper in the water column. For epipelagic Marinimicrobia, an NQR respiratory complex may be sufficient given that nitrogen and phosphorus availabilities and environmental stressors more often limit growth for these bacterioplankton than energy availability. Second, epipelagic waters have slightly higher pH and salinity (45), which may create a more favorable environment for the harnessing of a sodium motive force in surface waters and a proton motive force deeper in the water column. Last, it is possible that the reactive oxygen species (ROS) produced by NDH make the use of this complex disadvantageous in surface waters, where high hydrogen peroxide concentrations already generate substantial quantities of these stressors (46), though it is unclear if the NQR complex of Marinimicrobia produces fewer ROS.

A combination of vertical inheritance and lateral gene transfer (LGT) appears to have shaped the distribution of respiratory complexes throughout the Marinimicrobia. The NQR complex is prevalent throughout the Marinimicrobia phylogeny, including basal branching clades 9 and 10, suggesting that this complex was present in the common ancestor of all Marinimicrobia (Fig. 1d). Phylogenetic analysis of NqrA revealed a topology similar to that of the main Marinimicrobia clades, further suggesting that the NQR-(Na⁺) complex was present in the ancestral Marinimicrobia and has evolved primarily through vertical inheritance (Fig. 4e). The nNDH complex also appears broadly represented in Marinimicrobia, but its absence in basal branching groups 9 and 10 suggests that this complex either was acquired at the last common ancestor of clades 1 to 8 or was present in the last common ancestor of all Marinimicrobia and was then subsequently lost in clades 9 and 10 (Fig. 4f). The evolutionary history of the nNDH complex is broadly consistent with the marinimicrobial phylogeny, consistent with both of these scenarios. The distribution of the cNDH respiratory complex is the most restricted, with only selected mesopelagic Marinimicrobia harboring the gene cluster. Moreover, phylogenetic analysis of the NuoH subunit in cNDH revealed a phylogeny inconsistent with the Marinimicrobia phylogeny, suggesting that LGT is largely responsible for shaping the distribution of this gene cluster across the phylum. Among the members of the clade of cNDH NuoH proteins, clade 2 appears to have the most divergent sequences, suggesting that this gene cluster may have been acquired from the common ancestor of clades 2 to 7 and then transferred between clades afterward (Fig. 4f).

Calculation of genomic characteristics.

We predicted GC content, N-ARSC, C-ARSC, and mean intergenic space length data using previously described methods (32). Code for these analyses is available online (https://github.com/faylward/pangenomics/). To estimate genome size (S), we used the following formula:

where α is the number of base pairs in the genome assembly, β is the estimated level of contamination, and γ is the estimated level of completeness. We estimated contamination and completeness for each genome using CheckM v1.0.7 (54).

Protein cluster identification and annotation.

We predicted proteins from all genomes using Prodigal and subsequently identified protein orthologous groups (OGs) using proteinortho v5.16b with default parameters (61). For each OG, we chose the longest member as a representative and compared these proteins to those in the EggNOG release 4.5 (62), Pfam release 31 (63), and TigrFam release 15.0 (64) databases for annotation using HMMER3 (65). For EggNOG, we downloaded all NOG HMMs from the EggNOG website on 1 February 2018 and ran hmmsearch with an E value cutoff of 1e−5. For Pfam and TigrFam annotations, we used the noise cutoffs in each HMM as the lower bounds for annotation.

Respiratory complex annotation.

We annotated the NDH (H⁺) respiratory complex in a manner broadly similar to that previously reported (66). The canonical NDH (H⁺) respiratory complex consists of 14 subunits (nuoA to nuoN [nuoA-N]) that correspond to the COG HMMs COG0838, COG0377, COG0852, COG0649, COG1905, COG1894, COG1034, COG1005, COG1143, COG0839, COG0713, COG1009, COG1008, and COG1007. The subunits are usually syntenic, with the exception of nuoN and nuoM, which can sometimes be found on a distant chromosomal region or adjacent to other respiratory complexes. We identified protein OGs that corresponded to these COGs and considered a genome to contain this complex if at least 6 of the nuoA-L genes were present. We identified a second NDH respiratory complex in many genomes that lacked subunits D to F, and OGs corresponding to these subunits were distinct from those of the canonical NDH complex. We considered this second NDH complex to be present if at least 5 of the OGs corresponding to the nuoABCGHIJK subunits could be identified. For simplicity, we refer to the canonical nuo complex as cNDH and the noncanonical version lacking nuoDEF as nNDH. The canonical NQR-(Na⁺) respiratory complex consists of the 5 nqrA-F subunits, which correspond to COG HMMs COG1726, COG1805, COG2869, COG1347, COG2209, and COG2871. We identified OGs that corresponded to those COGs and considered a genome to encode the NQR complex if at least 3 OGs were present. Detailed information for each protein OG and their annotations and on which genomes encoded members can be found online (figshare.com/projects/Marinimicrobia_Pangenomics/30881).

Marker gene phylogenies.

To assess the evolutionary histories of key marker genes in the Marinimicrobia, we constructed phylogenies with these genes together with available reference sequences. For each marker gene, we identified the NOG to which the OG had been annotated using our EggNOG annotations and then downloaded all proteins belonging to the appropriate NOG on the EggNOG website (62). The one exception to this was the procedure used for the proteorhodopsin phylogeny, for which we used the reference sequences available on the MicRhoDE database (40). Because the reference protein data sets were quite large, we reduced their size by clustering similar proteins using CD-HIT (67) (default parameters). These reference proteins were then combined with the Marinimicrobia proteins into a single FASTA file, and phylogenies were constructed using the ETE Toolkit (56), with the standard_trimmed_fasttree workflow. We refer to these as the “full phylogenies” since they included a large number of reference sequences. For ease of visualization, we manually selected a subset of reference sequences together with all Marinimicrobia sequences from the full phylogenies and then constructed smaller “subset phylogenies.” Subset phylogenies, with Marinimicrobia proteins colored by clade, are provided in Fig. 4 (see also Fig. S4 at figshare.com/projects/Marinimicrobia_Pangenomics/30881). Full phylogenies are available as interactive trees at http://itol.embl.de/shared/faylward.

Marinimicrobia genome distributions and habitat distinctions.

To identify the biogeographic distributions of different Marinimicrobia genomes in global ocean samples, we downloaded 90 metagenome samples from the Tara Oceans expedition and mapped metagenomic reads against the final set of 218 Marinimicrobia genomes. We chose Tara Oceans samples to represent as broad a sampling of environments as possible and to include different depths (surface, deep chlorophyll maximum, and mesopelagic), ocean basins, and Longhurstian provinces. Details for the samples chosen are available online (figshare.com/projects/Marinimicrobia_Pangenomics/30881). We mapped reads using FastViromeExplorer (68), which, although initially intended for identification of viral sequences, includes a rapid and versatile read-mapping utility which contains built-in filters to remove spuriously identified sequences. We report final genome quantifications using the TPM (transcripts per kilobase million) metric (69), which corrects for sample size and reference genome length.

Habitat groups of Marinimicrobia were determined by hierarchical clustering of genome abundances in the Tara Ocean samples. For genome clustering, we loaded a log₁₀-transformed genome TPM abundance matrix into R and calculated pairwise Pearson correlation coefficients for the genomes using the “cor” function. We converted these correlations into distances by subtracting from a value of 1 and then clustered the genomes using the “hclust” command in R using average linkage clustering. We clustered samples using the same method. Heat maps and clustering dendrograms were visualized using the heat map.2 function in the gplots package.

We thank Jessica Bryant, Roderick Jensen, Stephen Melville, and Melanie Spero for helpful discussions and code sharing.

We also thank Benjamin Tully and the participants in a University of Southern California journal club who provided helpful comments on a preprint of this work.

We acknowledge use of the Virginia Tech Advanced Research Computing Center for bioinformatic analyses performed in this study.

This work was supported by an Alfred P. Sloan Research Fellowship in Ocean Sciences to F.O.A.

E.W.G. and F.O.A. performed analyses. S.S.T. and L.Z. contributed software and bioinformatic expertise. E.W.G. and F.O.A. wrote the paper.