Other effectors of SAR11 biogeography

Beyond temperature and latitude per se, many other factors are likely to influence the abundance of members of the SAR11 clade at a local level. Particularly in temperate and coastal zones the physical interactions of convergent water bodies with different histories, including rivers, may have a role in determining the final SAR11 community composition. Both Chesapeake (GS012) and Delaware (GS011) bays have very complex local hydrographic regimes with relatively extreme local variations in temperature and salinity, as well as large seasonal fluctuations (e.g., Shiah and Ducklow, 1994). Previous work in Chesapeake Bay has highlighted the presence of P1a.1, P1a.2 and P3.2 phylotypes (Kan et al, 2008). We identified these marine phylotypes, as well as P2.2, in GS012, GS011 and GS013.

Several time series analyses have also shown that the SAR11 clade as a whole (Carlson et al, 2009; Eiler et al, 2009; Gilbert et al, 2012), and the composition of subgroups (Brown et al, 2005; Carlson et al, 2009), varies in relative abundance in association with physical parameters (e.g., convective overturning—Morris et al, 2005; Carlson et al, 2009), environmental parameters (Brown et al, 2005; Carlson et al, 2009; Eiler et al, 2009; Gilbert et al, 2012), and with basin scale climatic events (Eiler et al, 2009). Despite methodological differences for studying SAR11 community composition, our data are concordant with many findings from these studies. Of note, at the Bermuda Atlantic Time Series (BATS) in the Sargasso Sea, pronounced vertical and seasonal patterns were observed in the distributions of three SAR11 subgroups: S1a, S1b and S2 (Carlson et al, 2009). S1a became more abundant than S1b in the surface waters (integrated 40m depth) 2 months after convective overturning events, and remained more abundant during the highly oligotrophic and high UV summer period. Our analysis shows P1a.3 is the dominant phylotype within the S1a subgroup in samples taken across a wide range of warm oligotrophic waters (Figure 1) indicating P1a.3 is likely the organism observed by Carlson et al (2009) and other studies (Wilhelm et al, 2007). Of further relevance, during periods of winter mixing at BATS (early January and late February), the S1b subgroup was most abundant (Carlson et al, 2009). In our analysis, S1b was most abundant in Sargasso Sea samples taken in February (Venter et al, 2004), and at an upwelling site (GS031) near Fernandina Island (Rusch et al, 2007). Deep mixing appears to have a significant effect on the relative abundance of S1b in a number of oceanic regions.

However, given the distributions of wind-driven upwelling/downwelling (i.e., Ekman pumping) and of mixed layer depths across the world’s ocean, upwelling or mixing effects alone do not provide the best explanation for the observed SAR11 phylotype distribution and abundance across the 128 samples analysed. The large-scale pattern of wind-driven vertical velocity consists of upward motion in cold high latitude regions and in the warm tropics, and downward motion in the subtropical gyres. If the phylotype distributions were strongly linked to the pattern of Ekman pumping, then we would expect a weak relationship with temperature (as both the warmest and coldest waters upwell, and waters of intermediate temperature downwell). Other factors, like eddies and convergent/divergent surface currents, can drive upwelling and downwelling, but there is no data with which to assess their global distribution. Similarly, the depth of the surface mixed layer and its seasonal variation is not a strong function of latitude or surface temperature. Deepest winter mixed layers are observed on the equatorward side of strong mid-latitude currents (the ‘mode’ waters), while relatively shallow mixed layers are observed in the tropics and polar latitudes. There is also a strong contrast in winter mixed layer depth between the North Pacific and North Atlantic, despite similar surface temperatures (e.g., Tomczak and Godfrey, 1994). Summer mixed layer depths are relatively constant from the tropics to high latitude while summer sea surface temperature is a strong function of latitude. The observed strong association between phylotypes and temperature/latitude argues against a strong dependence on the depth of the surface mixed layer.

It has been suggested that the differences in productivity in coastal versus open-ocean environments may be a primary driver in the distribution of organisms associated with phylotypes P1a.1 and P1a.3 (Schwalbach et al, 2010). This study (Schwalbach et al, 2010) examined the presence and abundance of genes associated with the glycolysis operon (which is present in the genome of HTCC1062, P1a.1 but not in HTCC7211, P1a.3) in 46 surface metagenomes, reporting the overrepresentation of these genes in coastal samples. In our analysis, on a global scale P1a.1 is predominant in polar regions and P1a.3 in tropical regions (Figure 1). Coastal versus open-ocean environments had a small but significant impact in temperate, but not in tropical or polar biomes, and this was caused primarily by the relative abundances of P3.2 and P1a.2. The overabundance of coastal metagenomes from temperate regions is therefore likely to bias analyses of SAR11 representation.

Chlorophyll did not have a strong role in defining the biogeography of SAR11 phylotypes (Supplementary Table S4). Although cold waters tended to contain higher chlorophyll levels than tropical waters, many of the Antarctic samples, typical of High Nutrient Low Chlorophyll regions, had low chlorophyll allowing us to de-convolute these parameters. Irrespective of chlorophyll levels, in cold water the polar phylotypes remained abundant, and in tropical waters the tropical phylotypes remained abundant.

Phylotype definition

The term ‘SAR11 clade’ is used to describe the entire collection of SAR11 type organisms, the term ‘subgroup’ (S) to describe a branch of the 16S rRNA gene tree (Table I; Figure 1; Supplementary Figure S1), and the term ‘phylotype’ (P) to describe the branches within subgroups defined using ITS sequences (Table I; Figure 1; Supplementary Figure S1). Subgroup and phylotype designations are based on previous work: Morris et al (2005) and Carlson et al (2009) for 16S rRNA gene-based SAR11 subgroups, and Garcia-Martinez and Rodriguez-Valera (2000) and Brown and Fuhrman (2005) for ITS-based phylotypes. The 16S rRNA gene tree in Supplementary Figure S1 was generated using the neighbour-joining algorithm in the software package ARB (Ludwig et al, 2004). Near full-length (>1300, bp) 16S rRNA gene sequences representing SAR11 subgroups (Morris et al, 2005) were obtained from Greengenes version 513274 (12 October 2010) alignment (DeSantis et al, 2006) and used to make the scaffold tree, and boostrap values were calculated for 1000 replicate trees using the programs Seqboot, DNAdist, Neighbour and Consense from the PHYLIP package (Felsenstein, 1993). Shorter sequences that also contained ITS regions (that were included in the ITS tree) were manually aligned and added using the maximum parsimony algorithm. The ITS trees in Supplementary Figure S1 were generated using neighbour-joining analysis of a 657 base-pair alignment containing 862 ITS sequences derived from SAR11 16-ITS clones and metagenomic data sets (Supplementary Table S1; Supplementary Figure S1). Individual SAR11 ITS sequences were 379–454bp in length. Bootstrap values were calculated as above. The annotated database and ITS alignment files are available from the authors upon request.

SAR11 subgroup 1a (S1a) consists of a number of phylotypes (Garcia-Martinez and Rodriguez-Valera, 2000; Brown and Fuhrman, 2005) we now define as P1a.1, P1a.2 and P1a.3. The 16S rRNA gene sequences associated with these three phylotypes all branch closely (~98–100% sequence similarity) with the first described species Candidatus Pelagibacter ubique, isolated off the coast of Oregon (Rappe et al, 2002). The separation of P1a.1 sequences from other ITS sequences was first identified using sequences from Antarctic waters (Garcia-Martinez and Rodriguez-Valera, 2000) and named SAR11-A. Ca. P. ubique, along with the isolates HTCC1002 (Rappe et al, 2002), HTCC8022, HTCC8010 and HTCC8048 (Stingl et al, 2007), all isolated from cool waters off the coast of Oregon, belong to P1a.1, as do ITS sequences originating predominantly from surface waters in the Arctic (Supplementary Table S1). A mosaic genome from metagenome data from Ace Lake, Antarctica, ACE_P1a.1 also belongs to this phylotype. P1a.2 was previously defined from ITS sequences originating primarily from the temperate waters at the San Pedro Ocean Time Series, off California and from the Mediterranean Sea. Given the clustering of sequences from different environments, along with moderate bootstrap support (44%) this group of sequences was previously designated as a separate phylotype named S1a (Brown and Fuhrman, 2005). P1a.2 also contains sequences from other Arctic, sub-Arctic and temperate surface waters, along with the strains HTCC8038, HTCC8041, HTCC8045, HTCC8046 and HTCC8049 which were isolated off the coast of Oregon. The phylotypes P1a.1 and P1a.2 are distinguished by the conserved region starting at position 422 in the alignment:

P1a.1 --ATCATTTATATCAATATCTATATCCGAACATT-AAAA-GT-TATTATTTA

P1a.2 TAATAATTAATATCAATATCTATATCCGAACATTAGTAATGT-TATTAGCTA

P1a.3 was originally defined from ITS sequences from temperate and tropical waters including the San Pedro Ocean Time Series, and Sargasso Sea and was previously named S2 (Brown and Fuhrman, 2005). P1a.3 also contains sequences from strains HTCC8051 and HTCC8047 isolated off the coast of Oregon, and HTCC7211, HTCC7215, HTCC7216 and HTCC7217, isolated from warm waters off Bermuda (Stingl et al, 2007). The genome sequence for HTCC7211 is available but unpublished (GenBank Accession: EF616619). SAR11 subgroup S1b contains the 16S rRNA gene sequence from the original SAR11 clone from the Sargasso Sea, and includes a single phylotype P1b. No cultures or genome sequences appear to be available for this phylotype. SAR11 subgroup 2 contains three phylotypes. P2.1 contains surface water sequences from temperate and tropical waters and was originally designated as S3 (Brown and Fuhrman, 2005). P2.2 was initially defined based on ITS sequences from the Antarctic and named SAR11-A21 (Garcia-Martinez and Rodriguez-Valera, 2000). P2.3 was defined based on ITS sequences from the deep waters in the Mediterranean Sea (Garcia-Martinez and Rodriguez-Valera, 2000) and originally named SAR11-D (deep). It contains sequences from waters of 70–3000, m depth (the majority from deep water) from tropical, temperate and polar biomes. No cultures or genome sequences appear to be available for these phylotypes. SAR11 subgroup S3 contains two marine phylotypes along with freshwater SAR11 LD group that are not discussed in our study. P3.1 includes the strain/genome HIMB114, isolated from warm coastal waters of Kanaohe Bay, Hawaii. P3.2 includes the Arctic strain/genome, IMCC9063, isolated from seawater off the coast of Svalbard, Norway (Oh et al, 2011) and the mosaic genome from Ace Lake, Antarctica, ACE_P3.2.

Other potential phylotypes were determined from our ITS phylogeny using our database. Several sequences that branch deeply near the SAR11 subgroup S1b (P1b.a, NA1, NA2 in Figure 1) are defined as SAR11 ITS sequences in GENBANK but generally lack a clear connection to the 16S rRNA gene phylogeny (i.e., a lack of 16S-ITS sequences). Only P1b.a contains a 16S-ITS linked clone (SPOTS_May03_160m_16) which is associated with the S1b subgroup (Supplementary Figure S1, upper panel). Hence the group of sequences within this collapsed node were designated as P1b.a (Figure 1). All sequences in collapsed nodes NA1 and NA2 are derived from deep Arctic waters (3000, m in the Greenland Sea and 500m in the Sub-Arctic North Pacific). The sequences for P1b.a, NA1 and NA2 were included in our BLAST analysis but received no hits from the metagenomes we analysed and were thus excluded from discussion.

Determination of phylotype abundance in metagenomic samples and biogeography

The 862 SAR11 ITS sequences were used to generate an annotated sequence database (each ITS sequence name was annotated with its phylotype designation) against which all metagenomic data sets were compared by blastn (e=0.0001). The phylotype affiliation of the high-scoring pairs for metagenomic sequences that matched the database with >95% sequence identity over at least a continuous 200 base-pair region was recorded, resulting in 2983 records. This multivariate data set (relative abundance of SAR11 phylotypes at each location) was used in further statistical analysis.

DistLM and dbRDA methods were used to analyse and model the relationship between the multivariate SAR11 ITS data set and available predictor variables. Not all variables were available for all samples. To make use of the most available data, missing data for salinity and chlorophyll were estimated using the expectation-maximization algorithm. However, all patterns reported were checked for consistency against the subset of 79 metagenomes for which all variables were present (Supplementary Tables S3 and S4). Variables were normalized (the values for each variable had their mean subtracted and were then divided by their standard deviation) to reduce the effect of different measurement scales. Bray-Curtis similarity matrices were generated from normalized, square root transformed data composed of the relative abundance of each SAR11 phylotype in each sample. MDS, ANOSIM, DistLM, dbRDA and expectation-maximization of variables analysis were performed using the Primer V6+PERMANOVA+ software (Clarke and Gorley, 2006).

Genome assembly and annotation

Initial shotgun metagenomic hybrid assemblies were generated as previously described (Ng et al, 2010; DeMaere et al, 2011; Lauro et al, 2011) from the 0.1μm fraction of the 5m depth sample from Ace lake, Antarctica (Lauro et al, 2011) using a combination of pyrosequencing and Sanger paired-end reads. Assemblies were manually inspected and validated using AMOS 3.1.0 (Phillippy et al, 2008) and Hawkeye 2.0 (Schatz et al, 2007). The mosaic draft genome of ACE_P1a.1 (~1.26Mbp) contained 79 contigs (N50=44.3kbp) that were assembled into 12 scaffolds (N50=150.2kbp) with 87.8% of the total base-pairs in contigs larger than 10kbp. The mosaic draft genome of ACE_P3.2 (~1.27Mbp) contained 129 contigs (N50=15.2kbp) that were assembled into 8 scaffolds (N50=329.6kbp) with 73.1% of total base-pairs in contigs larger than 10kbp. The scaffolds containing a 16S rRNA marker gene were used as phylogenetic anchors for manually binning the resulting scaffolds using a combination of paired-end reads, GC%, tetranucleotide frequencies and phylogenetic congruence using custom perl scripts. This approach yielded 12 scaffolds that could be confidently assigned to ACE_P1a.1 and 8 scaffolds that could be confidently assigned to ACE_P3.2. The scaffolds were annotated with the automatic annotation pipeline SHAP (DeMaere et al, 2011), as previously described (Lauro et al, 2011).

Genome alignment and recruitment

The genomic scaffolds were oriented and joined in alignments to known reference genomes with the 6-frame stop-codon spacer ‘NNNNCACACACTTAATTAATTAAGTGTGTGNNNN’ using the custom perl script scaffolding.pl to create a contiguous pseudomolecule. Whole genome alignments were performed by comparing the pseudomolecule using blastn (run with standard parameters and the -m8 flag and only matches of at least 30 nucleotides or longer used) and ACT (Carver et al, 2005) against reference SAR11 genomes. After removing all N’s from the query sequence, ANI and PCD between representative genomes of the SAR11 clade were computed using 1020 nucleotide fragments (Goris et al, 2007) by employing custom perl scripts. blastn was run using settings as previously described (Goris et al, 2007): parameters were used at the default settings, except X=150 (where X is the drop-off value for gapped alignment), q=21 (where q is the penalty for nucleotide mismatch) and F=F (where F is the filter for repeated sequences). Genome recruitment plots were obtained by comparing marine metagenomes using blastn (e-value <1e−10, sequence identity >90%) and visualized with the custom perl and R script plot_coverage.pl using a sliding window size of 1000, bp and a natural logarithmic scale for recruitment depth. We reported the recruitment of four representative metagenomes from the different temperature biomes (polar, temperate and tropical). The metagenome data sets were chosen because they were the most similar to each other in terms of sampling strategy, sequencing technology and sample size: Ace Lake polar (GS232, 539536 reads), Southern Ocean polar (GS362, 508628 reads), temperate (GS368, 512395 reads) and tropical (S_35155+S_35163, 576392 reads). Recruitment of metagenome data from other sampling sites gave similar results. The recruitment depth averaged over the length of each reference genome was used as a proxy for genome abundance in each data set. For statistical analyses an additional polar metagenome was included (GS358, 491340 reads) so that the polar data could be compared with tropical plus temperate data with equivalent data set sizes, and compared using a two-tailed t-test.

Detection of homology, COG identification and analysis

The orthologous predicted proteins were identified from each pair of SAR11 genomes using the reciprocal smallest distance (RSD) algorithm (Wall and Deluca, 2007) with a blastp e-value threshold of 1e−15 and a sequence divergence for the alignment of 0.5. The protein sequences that were not identified as orthologues during the first RSD run were then subjected to subsequent rounds of reciprocal comparisons and added to the list of paralogues until no more homologues could be identified. The remaining list of protein-encoding genes was considered genome specific. The list of genome-specific proteins was searched against version 2 of the database of all metagenomic 454 reads in CAMERA (Sun et al, 2011) (excluding Ace Lake metagenome data) with tblastn (e<1e−50 and >100 amino acids alignment). The hits were subdivided based on water temperature and statistically significant differences were analysed by 1000 resamplings with replacement of the blast hits at a confidence level of 99%. Briefly, 1000 hits were randomly sampled from all the tblastn hits for ACE_P1a.1 to the metagenomes in CAMERA. This subsample was compared with a similar subsample obtained from the hits of HTCC7211 and the difference in number of hits (binned according to water temperature) was calculated. This process was repeated 1000 times and the median difference in number of hits was calculated for each water temperature bin. To verify that this difference was unlikely to happen by random chance, this process was repeated except that each of the 1000 tblastn hits were sampled from the combined hits of HTCC7211 and ACE_P1a.1. This process was once again repeated 1000 times and the results were ordered from the lowest to the largest difference in hits for each temperature bin. From the combined data sets, the low confidence cutoff for 99% significance was the 10th median difference and the high cutoff was the 991st. If the results for the test comparison fell outside this interval, then the difference was considered statistically significant. This process was also followed for comparisons of ACE_P3.2 to HIMB114. Each protein was assigned to a COG as described previously (Lauro et al, 2011). Overrepresentation of specific COG categories compared with the genome background was verified by a resampling method as previously described (Allen et al, 2009; Lauro et al, 2011) over a confidence interval of 97% using the custom perl script COG_scrambler.pl.

Detection of proteins under positive selection

Each protein on the list of orthologues was aligned with clustalW (Thompson et al, 1994) with protein pairs that differed by >120 amino acids in length discarded from the analysis. The alignments were then back-translated to the corresponding nucleotide coding sequence using RevTans ( http://www.cbs.dtu.dk/services/RevTrans/download.php) and the dN/dS ratios were computed by the method of Nei-Gojobori (Nei and Gojobori, 1986) using yn00 from the PAML package (Yang, 2007) over the whole alignment and over a 60-codon sliding window. The highest dN/dS ratio was recorded. Orthologous protein pairs were considered under positive selection only if dN/dS ratio was >1 over at least one window of 60 amino acids. For the analysis of predicted membrane proteins, the protein topology was inferred by computing antigenicity (Hopp and Woods, 1983) and hydrophobicity (Eisenberg et al, 1984) indexes over a sliding window of 60 amino acids.

We acknowledge technical support for computing infrastructure and software development from Intersect, and in particular assistance from Joachim Mai, and acknowledge Matthew Lewis from the JCVI for his assistance with DNA sequencing. This work was supported by the Australian Research Council and the Australian Antarctic Division. Funding for sequencing was provided by the Gordon and Betty Moore Foundation to the JCVI.

Author contributions: RC initiated the research program. RC, TT, FML, MVB and JMH designed and performed the field expeditions. AA, JMH, CAP and ML undertook activities leading to the generation of the DNA sequence data. MVB, FML, SSR, LM, DW and RC conceived, designed or performed the data analysis. MVB, RC and FML wrote, and SSR, TT, MZD, MJR, LM and JAF contributed to the manuscript.