|Home | About | Journals | Submit | Contact Us | Français|
In this work, we study the diversity of bathypelagic microbial eukaryotes (0.8–20μm) in the global ocean. Seawater samples from 3000 to 4000m depth from 27 stations in the Atlantic, Pacific and Indian Oceans were analyzed by pyrosequencing the V4 region of the 18S ribosomal DNA. The relative abundance of the most abundant operational taxonomic units agreed with the results of a parallel metagenomic analysis, suggesting limited PCR biases in the tag approach. Although rarefaction curves for single stations were seldom saturated, the global analysis of all sequences together suggested an adequate recovery of bathypelagic diversity. Community composition presented a large variability among samples, which was poorly explained by linear geographic distance. In fact, the similarity between communities was better explained by water mass composition (26% of the variability) and the ratio in cell abundance between prokaryotes and microbial eukaryotes (21%). Deep diversity appeared dominated by four taxonomic groups (Collodaria, Chrysophytes, Basidiomycota and MALV-II) appearing in different proportions in each sample. Novel diversity amounted to 1% of the pyrotags and was lower than expected. Our study represents an essential step in the investigation of bathypelagic microbial eukaryotes, indicating dominating taxonomic groups and suggesting idiosyncratic assemblages in distinct oceanic regions.
The bathypelagic region of the deep ocean, defined as the water column between 1000 and 4000m, is the largest habitat on Earth, but remains relatively unexplored (Arístegui et al., 2009). Despite its huge volume, it is a relatively homogeneous environment in terms of environmental conditions (Angel, 1993). Thus, pressure for a given depth is constant, as well as temperature (range of −1 to 3°C), salinity (34.3 to 35.1) and dissolved oxygen concentration (2.4 to 5.7mgl−1). The deep ocean holds many and varied life forms, mostly linked to microbial food webs (Azam et al., 1983; Michaels and Silver, 1988; Tanaka et al., 2005), dependent on resources like particulate and dissolved organic matter and inorganic nutrients. These resources vary at broad regional scales across the bathypelagic zone depending on the flux of organic compounds from the surface, the age of the water mass and the occasional presence of benthic inputs, such as those from hydrothermal vents (Hansell et al., 2009; Hamme and Emerson, 2013), conforming the niche structure supporting microbial life (Nagata et al., 2010).
Heterotrophic microbial eukaryotes, generally considered as bacterial grazers, are the most important component of bathypelagic microbial communities after prokaryotes. In a recent study, we estimated the average cell abundance (14±1 cells ml−1) and biomass (50±14pg C ml−1) of microbial eukaryotes in bathypelagic samples (Pernice et al., 2015). These deep ocean assemblages were numerically dominated by cells between 2 and 5μm in size. It is worth to remember that these cells counts may be underestimated given that a percentage of cells may not survive the pressure changes during recovery (Edgcomb et al., 2011a). Yet, the diversity of these bathypelagic assemblages remains still poorly resolved.
Relatively few papers have analyzed the diversity of bathypelagic microbial eukaryotes using clone libraries of 18S ribosomal RNA (rRNA) genes and Sanger sequencing, an approach that is inherently limited in the amount of sequences generated. Most of these studies focused on the water column, always in a regional scale, either by using universal eukaryotic primers (López-García et al., 2001; Stoeck et al., 2003; Countway et al., 2007; Not et al., 2007; Edgcomb et al., 2011b) or group-specific primers (Bass et al., 2007; Lara et al., 2009). Other diversity studies were done in benthic environments, such as deep-sea subsurface and surface sediments (Edgcomb et al., 2011a; Salani et al., 2012) or in hydrothermal vents (Edgcomb et al., 2002; Sauvadet et al., 2010). The bathypelagic microbial eukaryotic diversity has also been studied by fluorescence in situ hybridization staining (Morgan-Smith et al., 2011; Morgan-Smith et al., 2013). Although fluorescence in situ hybridization provides useful quantitative information of given taxa, it remains restricted to the taxa targeted by the available probes, which depend on previous sequencing information, and therefore never targets the complete assemblage. High-throughput sequencing now allows a much more exhaustive assessment of microbial diversity. The Roche 454-pyrosequencing technique (Margulies et al., 2005) has been used in studies targeting the 18S rRNA gene of microbial eukaryotes inhabiting marine epipelagic waters (Edgcomb et al., 2011b; Logares et al., 2012; Kilias et al., 2013; Massana et al., 2015) and freshwaters (Charvet et al., 2012). This approach includes a PCR step, known to be prone to biases like DNA polymerase errors or primer selectivity, that may affect final amplicon ratios (Acinas et al., 2005). Yet, metagenomics now allows a PCR-free approach to investigate microbial diversity, which is based on extracting 18S ribosomal DNA (rDNA) reads (miTags when using Illumina) from metagenomes (Logares et al., 2014). Combining tag-sequencing and miTags provides a more robust view of microbial diversity, as tag-sequencing can recover low abundance taxa, whereas miTags provides a less biased representation of most abundant taxa.
Here we use pyrosequencing to investigate deep-sea microbial eukaryotic assemblages from the 0.8 to 20μm size fraction originating from 27 stations located in the Atlantic, Indian and Pacific oceans. Samples were taken during the Malaspina-2010 Circumnavigation Expedition (Pernice et al., 2015). Our study is the first attempt to describe exhaustively the global diversity of bathypelagic heterotrophic microbial eukaryotes by high-throughput sequencing.
Samples of bathypelagic heterotrophic microbial eukaryotes were collected in 27 stations visited during the Malaspina-2010 Circumnavigation Expedition (Figure 1), conducted from December 2010 to July 2011 on board R/V Hespérides (Duarte, 2015). Bathypelagic seawater was collected with Niskin bottles attached to a rosette sampling system fitted with a Seabird 0911Plus CTD probe, including temperature, salinity and oxygen sensors. Fourteen samples were taken at 4000m, 11 at depths between 3000 and 4000m and 2 were shallower: 2400m in station 62 and 2150m at station 82 because of bottom depth constrains. Seawater was prefiltered through a 200μm mesh placed at the end of the hose of the Niskin, and through a 20μm mesh in a funnel. Between 100 and 120l of the 20μm filtrate was then filtered with a peristaltic pump on 142mm Millipore (Billerica, MA, USA) polycarbonate filters of 0.8μm pore-size. Filters were flash frozen in liquid nitrogen and stored at −80°C until processed in the lab.
We sampled three principal bathypelagic water masses along the track of the Malaspina-2010 Expedition: North Atlantic deep water (NADW), Weddel Sea deep water (WSDW) and Circumpolar deep water (CDW). The proportion of these three water masses in each sample was inferred from its temperature, salinity and oxygen concentration (Catalá et al., 2015). Samples were then clustered in a dendrogram using the inferred deep-water proportions with Pvclust (Suzuki and Shimodaira, 2006), to identify the main water mass types sampled during the cruise. Six types were identified, NADW pure, NADW enriched, CDW pure, CDW enriched, CDW–WSD and WSDW enriched.
Filters were cut into small pieces and soaked into 3ml of lysis buffer (40mM EDTA, 50mM Tris-HCl, 0.75M sucrose) at the beginning of the DNA extraction. Afterward, enzymatic digestions were performed. Filters were incubated with Lysozyme (1mgml−1 final concentration) at 37°C for 45min while slightly shaken, and then Proteinase K (0.2mgml−1 final concentration) and sodium dodecyl sulfate (1% final concentration) were added and tubes were incubated at 55°C for 60min while slightly shaken. The resulting lysate underwent two steps of standard phenol–chloroform extraction to remove lipids and proteins. After the last centrifugation step, the aqueous phase was collected, concentrated in an Amicon Ultra unit (Millipore) and washed three times with 2ml sterile deionized water. After the third wash, between 100 and 250μl of purified total genomic DNA extract was recovered and quantified using a Nanodrop ND-1000 spectrophotometer (NanoDrop Technologies Inc., Wilmington, DE, USA). Typical yields were between 0.08 and 0.58μg DNA, corresponding to ~1.8ngl–1 of seawater.
A two-step PCR was required to obtain enough amplicons for pyrosequencing considering the low density of cells in the bathypelagic ocean. We added 2ng of genomic DNA to triplicate PCR tubes containing dNTPs (0.2mM), the eukaryotic primers TAReukFWD1 (5′-CCAGCASCYGCGGTAATTCC-3′) and TAReukREV3 (5′-ACTTTCGTTCTTGATYRA-3′) (Stoeck et al., 2010) at 0.5mM, as well as PCR buffer (1x) in a final volume of 20μl. Amplicon size was around 380bp (and variable among taxonomic groups). The first PCR consisted of an initial denaturation step of 98°C for 2.5min and then 10 cycles of 45s at 98°C, 35s at 53°C, 35s at 72°C and a final elongation of 10min at 72°C. Triplicate amplicons were pooled and purified using the NucleoSpin Gel and PCR Clean-up kit (Macherey-Nagel, Hoerdt, France) to a final volume of 18μl. Then, we did a second PCR step of 20 cycles (45s at 98°C, 35s at 48°C, 30s at 72°C; with a final step of 10min at 72°C) using 2μl of the previous PCR product as template. This second PCR included primers that were the same as before except that they had the 454-specific adaptors, and the forward primer also carried the barcode. We tested with an agarose gel that this second PCR worked correctly, and after this test we run the second PCR step for five replicates, each with 2μl of template, therefore using most of the first PCR product. Final PCR products were purified, eluted in 30μl of sterile deionized water and quantified with the Qubit 1.0 fluorometer (Invitrogen, Merelbeke, Belgium). About 200ng of PCR product were sequenced on a 454 GS FLX+ system (Lifesequencing S.L., Valencia, Spain; http://www.lifesequencing.com). Amplicon sequences produced in this study are publicly available at the European Nucleotide Archive (http://www.ebi.ac.uk/ena/data/view/PRJEB9943).
Quality control and demultiplexing of 454 data was processed with QIIME (Caporaso et al., 2010). Owing to the version of the 454-sequencer (GS FLX+), it was not possible to run a denoising algorithm (that is, DeNoiser or AmplicoNoise), so we increased the stringency in the cleaning process. Clean sequences had between 150 and 600bp in size, no more than 2 mismatches in the primer and no homopolymers longer than 8bp. Phred quality was analyzed by running windows of 50bp, and pyrotags were truncated to the last window that had an average Phred >25. Truncated pyrotags shorter than 150bp were removed. An operational taxonomic unit (OTU) table was constructed by clustering high-quality pyrotags at a 97% similarity threshold using UCLUST v1.22q (Edgar, 2010). The most abundant sequence within each OTU was selected as the OTU representative sequence. The table was manually curated by removing OTUs represented by short sequences (<250bp) and OTUs that occurred in only one sample, regardless of their abundance (Supplementary Figure S1). OTUs were then taxonomically assigned using three reference databases: SILVA 108 (Quast et al., 2013), PR2 (Guillou et al., 2013) and MAS9013, an in-house database (Pernice et al., 2013; available on request). SILVA 108 was used to exclude prokaryotic 16S rDNA OTUs, whereas PR2 and MAS9013 were used to classify eukaryotic OTUs mostly to the class level. OTUs were assigned to a given group when its representative sequence had a BLAST hit with an e-value below 10−100 with a reference sequence. Above this e-value, OTUs were classified as uncertain. Chimera check was done with UCHIME (Edgar et al., 2011) with default parameters using the MAS9013 database as reference, and the results were carefully evaluated to avoid removing novel OTUs. Thus, OTUs identified as chimeras were kept if their representative sequence had an e-value of 0 against the GenBank database, the similarity between parents was above 90% and the OTU occurred in at least 14 stations or was represented by at least 100 pyrotags.
Community analyses were performed with the R packages Vegan (Okasanen et al., 2013; analyses: non-metric multidimensional scaling (NMDS), Permanova, Shannon index, star plot) and GMT (Wessel et al., 2013; Geodistance test). The number of pyrotags per sample varied from 6625 to 29926, therefore comparative analyses were done using a second OTU table, subsampled at 6625 pyrotags per sample with rrarefy in Vegan. A Mantel test, based on Bray–Curtis distances and using a Pearson correlation coefficient, indicated that this subsampling set represented very well the original pool of sequences (R=0.98, P<0.0001, Supplementary Figure S2). Repeated NMDS plots based on Bray–Curtis distances were done to keep the plot displaying the minimal stress (Clarke and Ainsworth, 1993).
The same genomic DNA extracts were used to obtain metagenomes (except for stations 74, 88 and 109 because of technical reasons) by shotgun sequencing using a Illumina HiSeq2000 (San Diego, CA, USA) platform (2 × 100-bp sequencing) installed at the Department of Energy Joint Genome Institute (JGI) (Walnut Creek, CA, USA). Metagenomic reads including the 18S rDNA (miTags) were extracted as indicated in Logares et al. (2014). Only miTags longer than 100bp were kept (a total of 19434 miTags). These miTags were used as queries for a BLAST search against a data set composed by 454 sequences from the same samples. Following this taxonomic identification, miTags were clustered in supergroups (Alveolata, Rhizaria, Stramenopiles, Fungi, Excavata, Amoebozoa) plus Uncertain and Metazoa, and a cluster of other groups (mainly Apusozoa); and then compared with a similar analysis in pyrotags. In a second and more sensitive step, we blasted the 20 most abundant 454 OTUs against the miTag pool to assess their relative abundance. Comparisons were done by averaging the relative contributions in the same set of 24 samples. Metagenomic sequences are available on the Integrated Microbial Genome website (Gs0053074, http://img.jgi.doe.gov/) and a script for the extraction of miTags is available on request.
The diversity of bathypelagic microbial eukaryotes (size fraction of 0.8–20μm) analyzed by pyrosequencing 18S rRNA genes gave an initial OTU table of 405270 pyrotags clustered at 97% similarity in 7202 OTUs (Supplementary Figure S1). This value was reduced to 362910 pyrotags after eliminating OTUs represented by short sequences (<250bp) and OTUs occurring in only one sample. A posterior chimera check further reduced the number of OTUs by about 10% and that of pyrotags by about 1%. The final OTU table (2481 OTUs and 359140 pyrotags) was then used to analyze the community diversity in—and between—stations (alpha and beta diversity) and the taxonomic affiliation of bathypelagic microbial eukaryotes. The complete set of 27 samples was considered for the taxonomic analysis, whereas the shallower samples collected at stations 62 and 82 were excluded from the comparative community analysis, thereby based on 25 samples.
The richness observed oscillated between 400 and 1000 OTUs per sample, and rarefaction curves relating the OTUs detected depending on the sequencing effort showed that most samples did not show a sign of saturation (excepting perhaps samples from stations 23 and 131, Figure 2a). The Shannon index, providing an estimate of alpha diversity in each sample, ranged from 0.90 to 4.84 with a mean (±s.e.) value of 3.63±0.17 (Supplementary Figure S3). Shannon indices did not differ across basins, with 3.62±0.45 in the Atlantic, 3.56±0.17 in the Indian and 3.68±0.18 in the Pacific. Samples with the highest Shannon values (above 4.5) were sampled in the subtropical Atlantic (stations 10, 17, 19 and 43), a region that also included two of the samples with lowest diversity (stations 23 and 131). The different water masses also did not show distinct Shannon indices (details not shown). We also performed a rarefaction curve for pooled reads from all samples to produce a view of global diversity, considering the deep ocean as single environment (Figure 2b). The curve showed a clear sign of saturation, with a maximal richness estimated around 2500 OTUs, thereby providing a first-order conservative estimate of the global abundance of bathypelagic microbial eukaryotes OTUs.
To explore microbial eukaryotic community composition among the 25 deepest samples, we performed a NMDS analysis (Figure 3a), based on Bray–Curtis distances calculated from the subsampled OTU table of 6625 pyrotags per sample. Several NMDS plots were done and the one displaying the minimal stress (0.17) was kept. In an attempt to interpret the NMDS clustering, we performed a Permanova test with one continuous variable (distance of the sampling depth to the seafloor) plus three categorized variables (abundance ratio of prokaryotes to microbial eukaryotes, water mass and oceanic region) to identify the contribution of these factors to community structuring. The order of parameters matters so we followed a subordination criterion (ratio <water mass <bottom distance <ocean). We excluded station 10 from this test because we did not have microbial abundance data for this sample. The abundance ratio explained 21% of the variability with high significance (P=0.0001), whereas water mass explained 26% with lower significance (P=0.0018). The distance to bottom and the ocean region explained 6% each, with a significance of about 0.03. Thus, the variance explained by this set of variables was 59%.
The NMDS analysis showed a tendency of geographically closer stations to be more clustered in the NMDS plot. To better analyze the relationship between geographic distance and community composition, we performed a Mantel test between the Bray–Curtis dissimilarity matrix and the matrix of geographic distances (Figure 3b). The results showed a weak (R=0.32) but significant (P=0.001) correlation between both variables. In fact, geographically distant samples were almost always different, whereas closer samples were only seldom similar, displaying in some cases high differences (see also Figures 1 and and44).
Taxonomic analysis was performed at two geographical scales (local and global). To better display the differences in taxonomic composition among samples, we selected the 10 most abundant groups and displayed their relative pyrotag contribution in samples grouped by water mass types (Figure 4). The most evident feature is that fungal tags appeared to be abundant in Pacific and Indian bathypelagic waters (see also Figure 1). Taking together, Ascomycota and Basidiomycota represented on average 25% of pyrotags in CDW enriched samples and 26% in CDW–WSDW samples (with a maximal value of 70% in station 88), whereas they were less important in CDW pure samples (9%) and scarce in the pool of Atlantic samples (<5% in each water mass). Polycystinea (Collodaria and Spumellaria together) were dominant in stations belonging to different water masses, being better represented in NADW enriched (50%), WSDW enriched (46%) and in CDW pure (37%). Chrysophytes preferred NADW pure (47% of pyrotags, with a peak of 82% in station 131) and NADW enriched samples (11%). They were also contributors in the northern part of CDW samples (30% in stations 109-121). MALV-II were widespread in all water types, exhibiting a more constrained relative abundance than the previous groups (between 10% and 18%). The remaining groups (Amoebozoa, Dinophyceae, RAD-B and MALV-I) showed sporadic peaks in few selected stations.
Once analyzed the regional pattern, we decided to average these results in order to give a global view of the bathypelagic ocean. The analysis of pyrotags at the supergroup level indicated that the bathypelagic ocean was dominated by Rhizaria (35%) and Alveolata (25%), followed by Fungi (15%) and Stramenopiles (13%) (Figure 5a). Pyrotags included in the uncertain category represented 1% of the total. The analysis of OTU (clustered at 97% of similarity) was quite different (Figure 5b): Alveolata represented 53% of the OTUs, tripling the number of Rhizaria OTUs (18%), whereas Stramenopiles and Fungi had 5% and 3%, respectively. Interestingly, Excavata, with only 1% of pyrotags, accounted for 8% of the OTUs.
Analyzing our data at deeper taxonomic level, we identified 46 class-level groups, also including MAST and MALV ribogroups and the uncertain category (Table 1). At the global level, Collodaria (Rhizaria) was the group with the highest representation (96688 pyrotags in total contributing to 23.1% pyrotags per sample on average), followed by Chrysophytes (Stramenopiles; 12.2%), MALV-II (Alveolata; 13.1%), Basidiomycota (Fungi; 10.9%) and Dinoflagellates (Alveolata; 7.9%). Interestingly four of the five most abundant groups belonged to different Supergroups. In an attempt to show the degree of novelty within the classes, we also calculated the average of the similarity values of the OTUs against reference sequences: the lowest this average the highest the novelty. A parallel analysis was carried out to calculate the same averaged similarity but weighted by tag abundance. In both cases, Collodaria turned out to contain a high level of intragroup novelty, with an averaged similarity of reads toward reference sequences of 91% (Table 1). For the other dominant groups, this value was generally above 97%. In fact, 59% of taxonomically assigned OTUs had a similarity below 97% against reference sequences, highlighting a considerable level of novelty even within these known groups.
We analyzed the 20 most abundant OTUs in detail (Table 2). Surprisingly, seven of them had high similarity (>97%) with cultured organisms and were present in almost all stations. The most abundant was a colonial Collodaria (45,261 pyrotags) that was 90.3% similar to Collophidium ellipsoides, followed by a Basidiomycota (29099 pyrotags) and a Chrysophyceae (23512 pyrotags). The two fungal OTUs in the list were very similar to terrestrial strains (98.5% to Tilletiopsis minor and 98.7% to Engyodontium album). The first OTU of MALV-II appears in the thirteenth position. Whereas MALV-II was the most represented group in most samples, it did not contain any OTU with a high dominance, indicating a large diversity and high evenness of the group. Some OTUs detected are very similar to cultured species known to have a cyst stage in their life cycle, such as Pedospumella encystans and Platyamoeba contorta, indicating that the cells sampled may have been dormant or non active. There are evidences that this life stage could be also present in marine fungi (Richards et al., 2012).
A total of 128 OTUs were classified as uncertain because they presented e-values above 10−100 in hits against reference data sets (Table 1). We visually checked whether or not these pyrotags aligned with reference sequences (therefore having the V4-18S rDNA signature), which the majority did. For each uncertain OTU, we identified the closest environmental match and the closest cultured match in a BLAST search against nr-NCBI (Supplementary Table S1). Many of these uncertain OTUs could be assigned to established taxonomic categories (most of them within Rhizaria and Alveolate) in base of the closest cultured match hit, forming probably novel lineages within them. The most abundant group was Collodaria with 15 OTUs, followed by Cercozoa with 14 and Dinophyceae with 6. Some OTUs (11) were metazoan and are not discussed further here. A group of 19 OTUs formed a second level of novelty, as the sequence coverage in the BLAST search was below 50% even with the closest environmental match. These were assigned tentatively to the taxonomy of the closest cultured match as Novel-Group name. Finally, a third level of novelty was represented by a group of 27 OTUs with no BLAST hit in GenBank. Of these, only 5 appeared >10 times and the most abundant had only 229 pyrotags. These were named just as Novel, being potentially new undescribed high-rank groups.
MiTags from the same samples were used to compare the relative abundance of taxonomic groups inferred by metagenomics (a PCR-independent approach) and pyrosequencing. Of the 19434 miTags retrieved, only 3981 showed a similarity above 97% in an alignment of at least 100bp with the 454-sequencing pool. This percentage (20.5%) is equivalent to the coverage (in percentage) of the amplicon size versus the complete 18S rRNA gene (21.3%). In general, the percentage of supergroups was similar by the two approaches (Figure 6a). Rhizaria and Alveolata had a lower representation in miTags than in pyrotags, whereas Excavata had a much larger representation in miTags (10.7% versus 1.5%). The similarity between the two molecular approaches was generally well supported at local level as well (Supplementary Figure S4 and Supplementary Table S2).
To compare both approaches at a finer phylogenetic scale, we blasted the 20 most abundant OTUs against the pool of miTag sequences using the same similarity criteria (97% in at least 100-bp alignments). Three abundant Chrysophyte OTUs presented a high sequence similarity and many miTags affiliated indifferently to one of the three, so we pooled the data for them. Remarkably, the resulting 18 abundant OTUs recovered a large fraction of miTags. A total of 2581 miTags produced significant hits, representing 64% of the 3981 miTags that covered the V4 region. At the same time, these abundant OTUs comprised 62% of all pyrotags. In particular, OTU_146 belonging to Collodaria (45261 pyrotags) retrieved 829 miTags, OTU_6539 (Basidiomycota, 29099 pyrotags) retrieved 323 miTags and the composite of the three chrysophyte OTUs (23512 pyrotags) retrieved 419 miTags. When expressed as percentage of the respective data sets (normalized per sample), the relative abundance of these 18 OTUs was highly correlated in the two approaches (Figure 6b), with a very high significance (P<0.0001), an R2 of 0.70, a slope of 0.75±0.13 that is not statistically different from 1 (P=0.06), and an intercept of 0.99±0.51 that is not statistically different from 0 (P=0.07).
The results presented comprise, to the best of our knowledge, the first global assessment of the diversity of microbial eukaryotes assemblages in the bathypelagic realm. The geographic coverage sampled during the Malaspina-2010 Circumnavigation Expedition provides a global dimension including all except the polar oceans. Considering the important role of microbial eukaryotes in food webs (Massana, 2011; Richards et al., 2012), our main aim was to characterize global diversity patterns as a first step toward understanding the structure of bathypelagic microbial assemblages and the possible trophic roles of the taxa included. Pyrosequencing has proven to be an adequate tool to assess the composition of microbial communities (Charvet et al., 2012; Logares et al., 2012; Kilias et al., 2013). Our pyrosequencing results are well supported in terms of group relative abundance by miTags. Nevertheless, it is worth mentioning that our study provides a view gathered from screening environmental DNA extracts, whereas a different but complementary view could be obtained by surveying RNA extracts (Stoeck et al., 2007; Not et al., 2009; Massana et al., 2015). The permanence of DNA in the extracellular milieu, different life stages (including dormancy) or differential rDNA copy numbers among species, can generate such different views.
Despite the uniformity of some physical parameters, the deep ocean is not completely a homogeneous environment, as it contains varying amounts of organic matter and other chemical compounds (Hansell et al., 2009; Hamme and Emerson, 2013) that may shape community structure. Considered at a class level, only a few groups dominated the deep microbial eukaryote communities. Thus, Collodaria, Basidiomycota, Chrysophyceae, MALV-II and Dinophyceae accounted together for 69.6% of the pyrotags from the global deep ocean (Table 1) and always explained the majority of reads in each individual sample (Figure 1). The relative proportion of these five groups changed markedly between sites, revealing a high level of variability. At the OTU level, we identified only 42 OTUs that were found in all stations, but together they accounted for 50% of the pyrotags. Although particular classes and OTUs were widespread in the deep ocean, a typical global deep ocean community was not identified, because of the large variation in the relative abundances of these dominant groups.
Our descriptor of alpha diversity, the Shannon index, showed a global average of 3.63±0.17 and did not present a regular pattern along the cruise track nor a relation with environmental parameters. This average was lower than that found in surface marine coastal samples (4.77±0.08, CR Giner, personal communication) at a comparable sequencing effort and using the same sequencing method and data analysis, a difference that had statistical significance (P<0.0001). Rarefaction curves for each individual sample did not saturate, suggesting that our pyrosequencing effort was not exhaustively sampling the microbial eukaryotic diversity in each particular site (Figure 2a). However, this effort appeared satisfactory at a global level (Figure 2b), so we could assume that most abundant taxa present in the deep ocean were encompassed by our samples. In particular, our data suggested a total number of OTUs in the global bathypelagic ocean in the order of 2500. Even assuming that this number could be affected by PCR and sampling biases, we consider this as a relatively low number of taxonomic entities for the largest ecosystem in the planet.
Communities showed a large variability along the cruise, as seen by the differential dominance of taxonomic groups in separate stations (Figures 1 and and4)4) and the distribution of samples in the NMDS (Figure 3a). Geographic distance explained only a small part (10%) of the genetic variability among samples (Figure 3b), and in fact only a few geographically close stations had a very similar community composition (for example, 17 and 19; 112 and 121). Similarity among nearby samples has been found in other studies of marine microbial eukaryotes at spatial scales of thousands of kilometers (Scheckenbach et al., 2010; Salani et al., 2012) and it is likely caused by similar environmental conditions in nearby sites. In our study, the most important factor explaining microbial eukaryotic composition was the water mass type, which explained 26% of the variability observed, and a similar effect has been reported for marine prokaryotes (Agogué et al., 2008; Varela et al., 2008; Galand et al., 2009; Kirchman et al., 2010). If this was due to relevant differences in chemical composition of the water mass, then our data would agree with the general view of microbial biogeography obtained from surface waters (Massana and Logares, 2013), with the environment selecting for specifically adapted microbes. An alternative explanation would be that the particular history of the water mass is determining microbial composition, suggesting geographical restriction to microbial dispersal.
Unexpectedly the abundance of bacteria, which could be seen as the main food source for deep microbial eukaryotes, did not have a direct effect on taxonomic composition (data not shown). However, the situation changed when considering the ratio of prokaryotes to microbial eukaryotes, which explained 21% of the variability in the community structure with high significance. The average bathypelagic prokaryote:micro-eukaryote ratio (5000 prokaryote cells per microbial eukaryotic cell) is very high when compared with surface samples, where it is about 1000–2000 (Sanders et al., 1992; Pernice et al., 2015) suggesting that predator pressure on prokaryotes may be relatively small in the bathypelagic ocean, particularly for the samples with the highest ratio. Despite the abundances of deep prokaryotes and microbial eukaryotes are partially related, it is likely that microbial eukaryotic heterotrophy is also present in the bathypelagic ocean in the form of osmotrophy and parasitism. These alternative lifestyles are generally linked to particular taxonomic classes, such as fungi or MALV clades, yet a more meticulous taxonomic analysis is required to understand the microbial food-web dynamics in the deep ocean.
There are only a handful molecular reports of deep marine microbial eukaryotes produced by Sanger sequencing, and in these studies the assemblages were dominated by Alveolates, mainly MALV-II (López-García et al., 2001; Not et al., 2007) or Rhizarian (Countway et al., 2007; Sauvadet et al., 2010). The most abundant group in our survey was Collodaria, an order of the class Polycystinea (Rhizaria) that includes species with a colonial lifestyle and without silification (Ishitani et al., 2012). In surface waters, Collodaria normally bear photosynthetic endosymbionts and are categorized as organisms with phototrophic behavior (Stoecker et al., 2009), but we do not have information about their lifestyle in the deep, dark ocean where photosynthetic symbionts are not expected. Probably, the presence of colonies could explain the abundance of these tags in our data set. MALV-II were also abundant in our global survey. This ribogroup contains several examples of parasites, such as Amoebophyra spp. infecting dinoflagellates (Coats and Park, 2002), and parasitism probably is the trophic role characteristic of the entire group (Massana, 2011), perhaps including hosts other than dinoflagellates (Small et al., 2012). Interestingly in Atlantic waters, there was a positive correlation (R2=0.96, P=0.0005) between the relative abundances of MALV-II and Metazoan tags, which could indicate metazoans being infected by MALV-II parasites.
Besides the expected finding of rhizaria and MALV-II, our global survey also unveiled Chrysophytes and Fungi as important members of deep microbial eukaryotic assemblages. Chrysophytes, moderately represented in surface marine waters and seldom detected in deep-sea surveys (Massana and Pedrós-Alió, 2008; del Campo and Massana, 2011), was the second most abundant class. These organisms can be facultative or strict phagotrophs, osmotrophs, or phototrophs (Holen and Boraas, 1996). A few cultured isolates (Pedospumella encystans and Spumella danica) explained 79% of chrysophyte sequences. Chrysophytes are good candidates for being important bacterivores in deep waters.
Fungi have been reported to dominate in deep subsurface sediments on the basis of clone libraries (Edgcomb et al., 2011a) and pyrosequencing rRNA surveys (Orsi et al., 2013), and overall deep-sea fungal diversity appears dominated by Basidiomycota and Ascomycota with probable yeast lifestyle (Bass et al., 2007; Richards et al., 2012). As far as we know, our study is the first identifying fungi as a globally important eukaryotic group in deep seawater, with Basidiomycota being particularly abundant. Considering the careful sampling and methodological processing, and the fact that fungal tags followed a comprehensive distribution in the data set, we excluded the possibility of fungal contamination. Finally, our data highlighted that species of Excavata, particularly Diplonemids, can also be important components of bathypelagic microbial eukaryotic assemblages. These include many bacterivorous forms and are typically more abundant in deep waters than at the surface (Lara et al., 2009), representing up to 15% of fluorescence in situ hybridization-counted cells in bathypelagic microbial eukatyotic assemblages (Morgan-Smith et al., 2013). Excavata were poorly represented in the 454 survey, likely due to a putative negative bias against the group by the V4 primers used (Pawlowski et al., 2011), but reached up to 11% of miTags in the parallel metagenomic analysis. It is interesting to highlight that 60% of OTUs having a representative sequence longer than 450bp belonged to Excavata, and this longer sequence could cause negative bias during PCR amplification. Excavata cells, cryptic to our PCR-based approach, could thus be important phagotrophs in the deep ocean.
Compared with the better-studied surface microbial eukaryotic community (Massana, 2011), the first obvious difference in bathypelagic assemblages was the paucity of photosynthetic groups. Although sequences from Bacillariophyta, Bolidophyta, Dictyochophyta, Prasinophyceae, Prymnesiophyceae, Raphidophyta and Eustigmatophyceae were retrieved in our survey, together they only represented 0.14% of the pyrotags and 2% of the OTUs. The phototrophic signal in bathypelagic waters is likely due to sinking particles (Agusti et al., 2015), although we do not totally dismiss the idea of some species being facultative heterotrophs. The absence of photosynthetic groups explains the low contribution of Stramenopiles in the global deep community, whereas Alveolate and Rhizaria were represented with percentages similar as in surface (Moon-van der Staay et al., 2001; Not et al., 2007; Massana, 2011). Surprisingly, in a world dominated by heterotrophy, the MASTs ribogroups that are important phagotrophs in the surface ocean (Massana et al., 2014) were virtually undetected. For instance, one of the most abundant and widely distributed ribogroup, MAST-4 (Rodríguez-Martínez et al., 2009), was represented in our data set by only 18 pyrotags. The difference in the composition of taxonomic groups between surface and deep waters implies that there are different factors in the flux of energy and carbon through the microbial loop.
The amount of pyrotags of uncertain identification was only 1% (Figure 5). A tentative taxonomic assignation based on the closest environmental match showed that many of the uncertain OTUs affiliated with the most abundant taxonomic groups, that is, Rhizaria and Alveolata. The truly novel sequences formed only 27 OTUs (Supplementary Table S1). We also investigated the novelty within each class, and found that the novelty was very high in Collodaria, whereas two of the other dominant groups, Fungi and Chrysophytes, showed lower levels of novelty as seen by the analysis of their dominant OTUs that were rather close to cultured species. Taking into account the peculiarity of the deep ocean environment and the sampling of unstudied oceanic basins, the percentage of novel species discovered was lower than expected, suggesting that the bathypelagic ocean contains a limited set, about 2500 OTUs, of microbial eukaryotic taxa.
This study provides the first insight into general diversity patterns of microbial eukaryotes in the bathypelagic ocean. Thanks to the breadth of the Malaspina-2010 expedition we collected a large number of pyrotags from samples in separate marine regions that were further validated by a parallel metagenomic survey. We identified a few groups (Collodaria, Chrysophytes, MALV-II and Basidiomycota) dominating microbial eukaryotic communities but mixed in different proportions in each sample. Community differences were partially explained by the water mass they belong to and the abundance ratio between prokaryotes and microbial eukaryotes. As far as we know, the fact that Chrysophytes and fungi dominated some bathypelagic water samples was not reported before. Our research describes the global architecture of bathypelagic microbial eukaryotes communities and provides an essential underpinning to assess the role of eukaryotic microbes in the largest ecosystem on Earth.
This project was supported by the Spanish Ministry of Economy and Competitiveness through project Consolider-Ingenio Malaspina-2010 (CSD2008–00077) to CMD and FLAME (CGL2010-16304) to RM. RL was supported by a Juan de la Cierva Fellowship (JCI-2010-06594, Spanish Ministry of Economy and Competitiveness, Spain). Sequencing at the JGI was supported by US Department of Energy (DOE) JGI 2011 Microbes Program grant CSP 602 grant to SGA. The work conducted by the US Department of Energy Joint Genome Institute is supported by the Office of Science of the US Department of Energy under Contract No. DE-AC02-05CH11231. Additional funding was provided by the Spanish Ministry of Economy and Competitivity grant MALASPINOMICS (CTM2011-15461-E). The authors were also supported by the Generalitat de Catalunya Research Group grants 2014SGR/1179. We thank our colleagues Xose A. Álvarez-Salgado, Teresa S. Catalá, Guillem Salazar, Francisco Cornejo-Castillo, Encarna Borrull, Cristina Díez, Elena Lara, Irene Forn, Dolors Vaqué, Eugenio Fraile, Verónica Benitez, the crew and the chief scientists of the various cruises for helpful collaboration.
The authors declare no conflict of interest.
Supplementary Information accompanies this paper on The ISME Journal website (http://www.nature.com/ismej)