Decision letter: Single-amino acid variants reveal evolutionary processes that shape the biogeography of a global SAR11 subclade
Article Figures and data Abstract Introduction Results and discussion Materials and methods Data availability References Decision letter Author response Article and author information Metrics Abstract Members of the SAR11 order Pelagibacterales dominate the surface oceans. Their extensive diversity challenges emerging operational boundaries defined for microbial 'species' and complicates efforts of population genetics to study their evolution. Here, we employed single-amino acid variants (SAAVs) to investigate ecological and evolutionary forces that maintain the genomic heterogeneity within ubiquitous SAR11 populations we accessed through metagenomic read recruitment using a single isolate genome. Integrating amino acid and protein biochemistry with metagenomics revealed that systematic purifying selection against deleterious variants governs non-synonymous variation among very closely related populations of SAR11. SAAVs partitioned metagenomes into two main groups matching large-scale oceanic current temperatures, and six finer proteotypes that connect distant oceanic regions. These findings suggest that environmentally-mediated selection plays a critical role in the journey of cosmopolitan surface ocean microbial populations, and the idea ‘everything is everywhere but the environment selects’ has credence even at the finest resolutions. https://doi.org/10.7554/eLife.46497.001 Introduction The SAR11 order Pelagibacterales (Thrash et al., 2011; Ferla et al., 2013) is one of the most ubiquitous free-living lineages of heterotrophic bacteria in the world’s oceans (Giovannoni et al., 1990; Morris et al., 2002; Carlson et al., 2009; Eiler et al., 2009; Schattenhofer et al., 2009; Treusch et al., 2009). Successful cultivation efforts and single amplified genomes from the environment have led to studies revealing their critical role in marine carbon cycling (Rappé et al., 2002; Giovannoni et al., 2005; Stingl et al., 2007; Oh et al., 2011; Tsementzi et al., 2016; White et al., 2019), and environmental sequencing surveys have offered detailed insights into the ecology of this ancient branch of life in aquatic environments across the globe (Zinger et al., 2011; Brown et al., 2012). The evolution of SAR11 is an active area of research (Giovannoni, 2017) that is critically important to understanding the determinants of its remarkable ability to maintain abundant populations in the global ocean. The evolutionary origins of SAR11 and thus its precise placement in the Tree of Life is debated (Thrash et al., 2011; Rodríguez-Ezpeleta and Embley, 2012; Ferla et al., 2013; Viklund et al., 2013), and our understanding of the evolutionary processes that define the biogeography of SAR11 cells is not complete. At the level of major SAR11 clades, previous studies have attributed markedly distinct patterns of distribution in the global ocean to both niche-based (Brown et al., 2012; Eren et al., 2013a) and neutral processes (Manrique and Jones, 2017). At the level of individual populations, a key simulation by Hellweger et al. (2014) showed that the intra-population sequence divergence that reflects the geographic patterns of distribution for SAR11 cells could emerge solely as a function of ocean currents, without selection (Hellweger et al., 2014). Between the extremes of inter-clade and intra-population diversity lies a wealth of variation that potentially can yield insights into the ecological and genetic forces that determine genomic diversity and fitness between closely-related, naturally occurring SAR11 populations. High-throughput sequencing of metagenomes provides access to genome-wide heterogeneity within environmental populations (Simmons et al., 2008), and current computational strategies can reveal associations between ecological parameters and microdiversity patterns at various levels of resolution (Eren et al., 2015; Scholz et al., 2016; Nayfach et al., 2016; Costea et al., 2017; Truong et al., 2017). However, SAR11 poses multiple challenges for such investigations, including their remarkable intra-population genomic diversity and the limited success of reconstructing SAR11 genomes from metagenomic data. Comprehensive investigations of the genetic contents of naturally occurring microbial populations (see Denef, 2018) for a review) often rely on population genomes directly reconstructed from metagenomes (Simmons et al., 2008; Bendall et al., 2016; Anderson et al., 2017; Garcia et al., 2018). While advances in genome-resolved metagenomics have made microbial clades more accessible without cultivation (Spang et al., 2015; Brown et al., 2015; Anantharaman et al., 2016), reconstructing SAR11 genomes from the surface ocean remains a difficult endeavor, as evident in recent comprehensive surveys of metagenome-assembled genomes (MAGs) from seawater samples from around the globe (Tully et al., 2018; Delmont et al., 2018). In the absence of population genomes recovered directly from the environment, genomes from isolates can also offer insights into environmental populations through genome-wide recruitment analyses in which short metagenomic reads are aligned to a reference (Denef, 2018). Using metagenomic read recruitment to investigate the structure of environmental populations is confounded by the challenge of defining the boundaries of microbial populations. Without an established species concept in microbiology, defining units of microbial diversity and their boundaries is a significant challenge (see Shapiro, 2018 and Cohan, 2019 for discussions). Nevertheless, from analyses of isolated microbial strains with formal taxonomic descriptions, a genome-wide average nucleotide identity (gANI) cutoff of 95% emerged as an operational delineation of species (Konstantinidis and Tiedje, 2005; Varghese et al., 2015) and was confirmed in a recent analysis of eight billion pairwise comparisons of whole genomes (Jain et al., 2018). Both gANI calculations using complete genomes, as well as the average nucleotide identity of metagenomic short reads (ANIr) recruited from environmental metagenomes using reference genomes, show an interesting discontinuity among sequence‐discrete populations at sequence identity levels between 80% and 90–95% (Konstantinidis and DeLong, 2008; Caro-Quintero and Konstantinidis, 2012; Jain et al., 2018). Regardless of their theoretical significance, these cutoffs are essential for multiple practical purposes, such as the identification and subsequent exclusion of metagenomic reads that originate from non-target environmental populations, to avoid inflating variants arising from contaminating non-specific reads in microbial population genetics studies. Interestingly, the boundaries of environmental SAR11 populations appear to not comply with the 95% ANIr cutoff. For instance, Tsementzi et al. (2016) observed substantial sequence diversity within sequence-discrete SAR11 subclades in the environment, and suggested that an ANIr as low as 92% would be required to adequately define the boundaries of the SAR11 populations recovered in their study (Tsementzi et al., 2016). These findings are consistent with a comprehensive study of isolate genomes and marine metagenomes by Nayfach et al. (2016), which suggested that SAR11 is one of the most genetically heterogeneous marine microbial clades (Nayfach et al., 2016). The substantial sequence diversity within environmental SAR11 populations not only explains the absence of SAR11 population genomes in genome-resolved metagenomics studies, but also challenges conventional approaches to the study of population genetics in microorganisms. For instance, the multiple occurrence of single-nucleotide variants in individual codon positions would render commonly used computational strategies that classify synonymous and non-synonymous variations based on independent nucleotide sites (such as in Schloissnig et al., 2013; Bendall et al., 2016) unfeasible. Despite these challenges, SAR11, with its ubiquity in surface seawater samples, extensive diversity in sequence space, and unique evolutionary history, remains one of the exciting puzzles of contemporary microbiology. Here we investigated the evolutionary processes that maintain genetic diversity within a natural SAR11 lineage accessible through a single isolate genome that recruited more than 1% of surface ocean metagenomic reads from a global dataset. Using single-amino acid variants, we were able to (1) delineate multiple proteotypes whose distributions were more closely linked to large-scale oceanic current temperatures than they were to geographic proximity, and (2) resolve positive and negative selection mediated by temperature and its co-variables. Our findings suggest that environmentally mediated selection, rather than neutral processes, dominate the biogeographic partitioning of SAR11 at fine scales of taxonomic resolution. Our study also offers new computational approaches to characterize variation within complex microbial populations, including additional means to integrate amino acid and protein biochemistry into microbial population genetics. Results and discussion To find the most appropriate SAR11 isolate genome to study the population genetics of naturally occurring SAR11, we used the complete genomes of 21 SAR11 isolates in a competitive recruitment of short reads from 103 metagenomes. Most of these metagenomes were from the TARA Oceans Project (Sunagawa et al., 2015), and correspond to 93 stations across four oceans and two seas. We also included an additional 10 metagenomes from the Ocean Sampling Day Project (Kopf et al., 2015) to cover high-latitude areas of the Northern hemisphere. All metagenomes correspond to small planktonic cells (0.2–3 μm in size) from the surface (0–15 meters depth; n = 71) and deep chlorophyll maximum (17–95 meters depth; n = 32) layers of the water column (Supplementary file 1a). The isolates we used belonged to SAR11 subclades Ia.1 (n = 6), Ia.3 (n = 11), II (n = 1), IIIa (n = 2) and the related alphaproteobacterium Va (n = 1) (Supplementary file 1b), which collectively recruited 1,029,716,339 reads from all metagenomes, or 3.3% of the dataset (Supplementary file 1c). The metapangenome of SAR11 To investigate associations between ecology and gene content of SAR11 lineages, we first performed a pangenomic analysis in conjunction with read recruitment from the metagenomic data. The pangenome of SAR11 genomes consisted of all 29,719 genes grouped into 6175 gene clusters (Supplementary file 1d). The clustering of genomes based on shared gene clusters (Supplementary file 1e) matched that of the previously described phylogenetic clades (Grote et al., 2012) (Figure 1A; an interactive version of which is available at http://anvi-server.org/p/4Q2TNo). The SAR11 pangenome across metagenomes (i.e., the SAR11 metapangenome) revealed distinct distribution patterns for each clade within SAR11 (Figure 1A). Clade Ia recruited the most reads compared to other clades (Supplementary file 1b), consistent with previous studies that found this clade to be highly abundant in surface seawater (Field et al., 1997; Brown et al., 2012; Eren et al., 2013a; Manrique and Jones, 2017). Gene clusters divided clade Ia into two main clusters corresponding to the high-latitude subclade Ia.1 and the low-latitude subclade Ia.3 (Figure 1A). While all high-latitude genomes displayed a bi-polar geographic distribution in the metagenomic dataset, gene clusters in low-latitude genomes revealed multiple sub-groups that also showed different patterns of geographic distribution (Figure 1A). This emphasized the need to further refine subclade 1a.3, in which each genome pair had over 98.6% sequence identity at the 16S rRNA gene level (Supplementary file 1f). Our consideration of geographical co-occurrence patterns, phylogenomic characteristics, and pangenomic properties in this metapangenome revealed six subclades within 1a.3 with cultured representatives (Figure 1A, also see Supplementary file 1g for gANI estimates between SAR11 genomes). We tentatively name them SAR11 subclade 1a.3.I (HTCC7211, HTCC7214 and HTCC7217; gANI of >93% and 16S rRNA gene identity of >99.4%), 1a.3.II (HIMB5), 1a.3.III (HIMB4 and HIMB1321; gANI of 94.8% and 16S rRNA gene identity of 100%), 1a.3.IV (HTCC8051 and HTCC9022; gANI of 86.9% and 16S rRNA gene identity of 100%), 1a.3.V (HIMB83) and 1a.3.VI (HIMB122 and HIMB140; gANI of 94.6% and 16S rRNA gene identity of 99.7%). Overall, the refinement of SAR11 subclades reveals a striking agreement between phylogeny, pangenome, and the ecology of the members of the SAR11 clade Ia. Figure 1 with 1 supplement see all Download asset Open asset The SAR11 metapangenome. Panel A describes the pangenome of 21 SAR11 isolate genomes based on the occurrence of 6175 gene clusters, in conjunction with their phylogeny (clade level) and relative distribution of recruited reads in 103 metagenomes ordered by latitude from the North Pole to the South Pole (top right heat map). The relative distributions were displayed for a minimum value of 0.1% and a maximum value of 1%. The layer named ‘Core 1a.3.V genes’ displays the occurrence of the 799 core 1a.3.V genes (in green) and those found in HIMB83 but not in the 1a.3.V lineage (in purple). Panel B describes the relative distribution of reads the 799 core 1a.3.V genes recruited across surface metagenomes from TARA Oceans. https://doi.org/10.7554/eLife.46497.002 A remarkably abundant and widespread SAR11 lineage at low latitudes While Ia.3 was the most abundant SAR11 subclade in our dataset, the new subclades we defined in this group differed remarkably in their competitive recruitment of short reads from metagenomes (Figure 1A, Supplementary file 1b). For example, while the least abundant subclade (1a.3.II; represented by HIMB5) recruited 22.6 million reads, the most abundant one (1a.3.V; represented by HIMB83), recruited 390.9 million reads, or 1.18% of the entire metagenomic dataset (Supplementary file 1b). For perspective, this is roughly two times more reads than the most abundant Prochlorococcus isolate genome recruited from the same dataset (Delmont and Eren, 2018) (Supplementary file 1h). Strain HIMB83 contains a 1.4 Mbp genome with 1470 genes, and was isolated from coastal seawaters off Hawai’i, USA. But it also recruited large numbers of reads from locations that were distant to the source of isolation (Supplementary file 1c). The gANI between HIMB83 and the most similar genome in our dataset, HIMB122 (1a.3.VI) was 82.6%, and the remarkable abundance of HIMB83 has also been recognized by others (Brucks, 2014; Nayfach et al., 2016). To the best of our knowledge, 1a.3.V is the most abundant and widespread SAR11 subclade in the euphotic zone of low-latitude oceans and seas. Although it is a member of the subclade 1a.3.V, the genomic context HIMB83 provides does not exhaustively describe the gene content of all members of 1a.3.V. Nevertheless, it gives access to the core 1a.3.V genes through read recruitment. To identify core 1a.3.V genes, we used a conservative two-step filtering approach. First, we defined a subset of the 103 metagenomes within the main ecological niche of 1a.3.V using genomic mean coverage values (Supplementary file 1c). Our selection of 74 metagenomes in which the mean coverage of HIMB83 was >50X encompassed three oceans and two seas between −35.2° and +43.7° latitude, and water temperatures at the time of sampling between 14.1°C and 30.5°C (Figure 1—figure supplement 1, Supplementary file 1i). We then defined a subset of HIMB83 genes as the core 1a.3.V genes if they occurred in all 74 metagenomes and their mean coverage in each metagenome remained within a factor of 5 of the mean coverage of all HIMB83 genes in the same metagenome. This criterion accounted for biological characteristics influencing coverage values in metagenomic surveys of the surface ocean such as cell division rates and variations in coverage as a function of changes in GC-content throughout the genomic context. Figure 1—figure supplement 1 displays the coverage of all HIMB83 genes across all metagenomes, and Supplementary file 1j reports underlying coverage statistics. While the 799 genes that met these criteria systematically occurred within the niche boundaries of 1a.3.V, 40% of the remaining 671 HIMB83 genes that were filtered out were present in five or fewer metagenomes and coincided with hypervariable genomic loci (Figure 1—figure supplement 1). Hypervariable genome regions are common features of surface ocean microbes (Coleman et al., 2006; Zaremba-Niedzwiedzka et al., 2013; Kashtan et al., 2014; Delmont and Eren, 2018) that are not readily addressed through metagenomic read recruitment but do influence pangenomic trends. Here, less than 10% of gene clusters unique to HIMB83 were among core 1a.3.V genes (Figure 1A), indicating HIMB83’s unique genes are mostly accessory to the members of 1a.3.V. In contrast, more than 80% of gene clusters that were core to the 21 SAR11 genomes matched to the core 1a.3.V genes. The overlap between environmental core genes of 1a.3.V revealed by the metagenomic read recruitment and the genomic core of SAR11 revealed by the pangenomic analysis of isolate genomes suggests that these genes represent a large fraction of the 1a.3.V genomic backbone (Figure 1A). Core 1a.3.V genes recruited on average 1.25% of reads in the 74 metagenomes (Figure 1B, Supplementary file 1j). The broad geographic prevalence of core 1a.3.V genes represents a unique opportunity to study the population genetics of an abundant marine microbial subclade across distant geographies. SAR11 subclade 1a.3.V maintains a substantial amount of genomic heterogeneity To investigate the amount of genomic heterogeneity within 1a.3.V, we first studied individual short reads that the HIMB83 genome recruited from metagenomes. The percent identity of reads that matched to the 799 core 1a.3.V genes ranged from 88% to 100% (Figure 2), which is considerably more diverse than those observed in similar reference-based metagenomic studies (Konstantinidis and DeLong, 2008; Tsementzi et al., 2016; Meziti et al., 2019). Notably, we also observed similar trends for the other SAR11 genomes included in this study (Figure 2—figure supplement 1), suggesting that the relatively high sequence diversity observed among core 1a.3.V genes may be a characteristic shared with other SAR11 lineages in the surface ocean. Figure 2 with 2 supplements see all Download asset Open asset Statistics of recruited reads. Left panel shows percent identity distributions in each of the 74 metagenomes. Curves are colored based on height. Metagenomes are ordered according to how the percent identity distributions hierarchically cluster based on Euclidean distance (dendrogram). Right panels display a summary of distribution statistics for each percent identity distribution compared against in situ temperature in a linear regression (correlations to all other available parameters are summarized in Figure 2—figure supplement 2). Each point is a metagenome and black lines are lines of best fit. For visual clarity, the data in left panel considers only the median read length and interpolates between data points, whereas the data in right panels consider all read lengths with no interpolation. https://doi.org/10.7554/eLife.46497.004 Overall, our data confirm that ANIr values of >95% used previously to delineate sequence-discrete populations does not apply to SAR11. One immediate implication of this substantial amount of sequence diversity that defies previous empirical observations is our inability to explicitly define what we are accessing in the environment. This challenge is partially because a precise and exhaustive description of what constitutes a ‘population’ remains elusive (Cohan and Perry, 2007; Shapiro and Polz, 2015; Cohan, 2019), which creates significant practical challenges (Rocha, 2018), such as the accurate determination of the boundaries of naturally occurring microbial populations especially in metagenomic read recruitment results. Nevertheless, the term ‘population’ is frequently used in literature (Simmons et al., 2008; Kashtan et al., 2014; Bendall et al., 2016), which implies that Charles Darwin’s observation in his historical work ‘On the Origin of Species’ continues to summarize our struggle in life sciences to describe theoretical boundaries of fundamental units of life even though contemporary enviornmental microbiology has gone beyond the term species in this pursuit: 'no one definition of species has yet satisfied all naturalists; yet every naturalist knows vaguely what [they mean] when [they speak] of a species’ (Darwin, 1859). Our study is not well-positioned to offer a precise theoretical definition for the term 'population', either. Instead, similar to previous studies, we resort to an operational definition that suggests a population is 'an agglomerate of naturally occurring microbial cells, genomes of which are similar enough to align to the same genomic reference with high sequence identity’ (Delmont and Eren, 2018; also see Denef, 2018 and references therein for a comprehensive discussion of what constitutes a population from a metagenomic perspective). By outsourcing the hypothetical radius of a population in sequence space to the minimum sequence identity of short reads recruited from metagenomes, this approach offers a practical means to study very closely related environmental sequences without invoking theoretical The broad heterogeneity that no sequence-discrete we observed within the sequence defined this the metagenomic reads that to HIMB83 genes (Figure 2), the that this within a population (Figure 1—figure supplement 1). However, to the theoretical and with the of short metagenomic reads, in we more that our reads originate from multiple closely related yet SAR11 populations within subclade 1a.3.V. Both high rates between cells low gANI values and of genes between distinct clades could the of between SAR11 populations in the surface ocean et al., The high of closely related 1a.3.V cells in the surface ocean suggests the of these two forces could be high within populations as At least two extensive SAR11 sequence diversity and in understanding its One is that the members of 1a.3.V we access are in the of into multiple sequence-discrete populations and we are an emerging in the evolutionary journey of SAR11. the observed diversity may represent a of sequence variants to a et al., 2012). To these we the between properties of these (i.e., and and environmental parameters linear regression (Figure 2—figure supplement Supplementary file This analysis revealed a significant between in situ temperature and distribution which suggests a influence of temperature and its on the sequence heterogeneity within 1a.3.V (Figure 2) and is with the of sequence of non-synonymous variation identity distributions are to statistics of short reads to a they do not information their significance, or with To this we a to characterize amino acid in metagenomic data and to study genomic variation that amino acid sequences (see Materials and our approach only metagenomic short reads that cover all three in a codon to determine the of single-amino acid variants (SAAVs) in protein While is a codon in it is often from a single-nucleotide with the that the two remaining are However, populations with extensive nucleotide variation can this in the of the core 1a.3.V genes, on average of metagenome with other in the same of codon sequences as in the is a to the 799 core 1a.3.V genes and 74 metagenomes, we SAAVs in which of amino from the (i.e., the most amino acid for a codon and The of codon positions that a of core 1a.3.V genes and with on across the 74 metagenomes Figure 1—figure supplement 1 and Supplementary file and not in metagenomes to the source of isolation (Supplementary file and suggesting that the of isolation for HIMB83 does not the biogeography and population genetics of 1a.3.V. To we codon positions if their coverage in of the 74 metagenomes was which in a of SAAVs occurred in codon positions that a in at least one metagenome among the of codon positions within the core 1a.3.V genes (Supplementary file We a protein to be (i.e., absence of variation to purifying or positive in a metagenome if it were in our in we of all across the 74 that encompassed only genes (Supplementary file In all genes, one nucleotide at least one in at least one metagenome (Supplementary file revealing a of amino acid sequence among core 1a.3.V the of purifying selection on amino To how commonly each amino acid was found in we compared the amino acid of SAAVs to the amino acid of the core 1a.3.V genes (see Materials and In a in which amino are as common in SAAVs as they are across all 799 core genes, the that an amino acid occurred in SAAVs would with its within the core genes While these were we observed large from this occurrence of amino in SAAVs relative to their occurrence in core genes (Figure Figure supplement 1, Supplementary file All and amino were in SAAVs compared to the core 1a.3.V genes (Figure For instance, while made only of all amino in the core genes, on average of SAAVs across the 74 metagenomes (Supplementary file Interestingly, amino amino not substantial between core 1a.3.V genes and amino were or in SAAVs with to their within core genes. In contrast, all amino with the very of and were in SAAVs (Figure Figure supplement 1, Supplementary file within the core of are to be critical for the required for and which a purifying selection on occurring in positions et al., and 2005; et al., 2009). amino the of they are on average purifying selection, which is the for the of amino within the other in positions on the surface of are as they are less to protein Overall, our analysis revealed that the occurrence of amino in SAAVs is roughly with the occurrence of amino within the core 1a.3.V genes, and that from this are in by levels of purifying selection that the of an amino for a environment (Figure supplement 1). acid rates reveal of and
Read more