Insights into global diatom distribution and diversity in the world’s ocean

At a cutoff level of 85% identity to sequences in our reference database (39), a total of 63,371 V9 rDNA ribotypes (represented by ∼12 million sequence reads) from 293 communities could be assigned to diatoms. Rarefaction analysis indicated that these ribotypes approached saturation at a global scale (Fig. 2A) although individual oceanic regions, such as the NAO and RS, were far from saturation. Preston log-normal distribution extrapolated the true diatom ribotype richness to 96,710 ribotypes (fitted red curve in Fig. 2B), suggesting that our survey retrieved ∼66% of diatom ribosomal diversity in the photic zone of the global ocean (shaded region in Fig. 2B). All of the ribotypes were clustered (36, 40) into biologically meaningful operational taxonomic units (OTUs), which yielded 3,875 distinct OTUs. Each OTU was represented by the most abundant ribotype in the OTU cluster. For these OTUs, Preston’s veil revealed the completion in sampling to be 81.6%, with an extrapolated number of OTUs to be 4,748 (SI Appendix, Fig. S3).

Based on ribotype abundance, diatoms were found to be one of the most represented eukaryotic lineages [number two in eukaryotic phototrophic lineages (after the Dinophyceae, although note that there are many taxa of Dinophyceae that are not photosynthetic at all or perform photosynthesis only facultatively) and number five with respect to all marine eukaryotic lineages] (36). Overall, diatom reads accounted for about 2.86% of total eukaryotic reads and 4.86% of protist ribotypes in our set of samples but represented more than 25% of the total eukaryotes at some locations: e.g., in the SO (Fig. 2C). Diatoms contributed ∼75% to the total photosynthetic community at station 11 (MS), more than 78% and 65% at polar stations 84 and 85, respectively (SO), 44% at subpolar station 82 (SO), and more than 38% and 44% at stations 122 and 123, respectively (Marquesas Islands; tropical SPO), and globally represented 27.7% of the total eukaryotic photosynthetic planktonic community. The mean percentages of diatom reads across 46 stations were 2.6%, 5%, and 19.9% with respect to the total eukaryotic reads, protistan reads, and photosynthetic reads, respectively (Fig. 2C). Many tropical and subtropical stations in the MS (stations 18, 20, and 30), inner RS (stations 31, 32, and 33), IO (stations 41, 45, and 48), subtropical SAO (stations 72, 76, and 78), and in the SPO subtropical gyre (station 98) were found to be very scarce in diatom sequences in comparison with other photosynthetic groups, such as dinoflagellates and haptophytes (Fig. 2C and ref. 36).

Nearly 58% of the reads (corresponding to 33,314 ribotypes) could be assigned at least down to genus level, and the large majority (>90%) of these assigned sequences belonged to known planktonic genera (SI Appendix, Fig. S2). Of the 79 genera found, Chaetoceros was the most abundant genus, representing 23.1% of total assigned sequences. Fragilariopsis accounted for 15.5% of total assigned sequences, followed by Thalassiosira (13.7%) Corethron (11%), Leptocylindrus (10.1%), Actinocyclus (8.7%), Pseudo-nitzschia (4.4%), and Proboscia (3.9%) (Fig. 3, column a and Dataset S1). Only a few sequences were assigned to genera known from freshwater or benthic environments, and in most cases only with low similarity (e.g., Fragilariforma and Epithemia) (SI Appendix, Fig. S2), likely because of the lack of reference sequences for a number of marine planktonic genera (see Unassigned Sequences and Comparison Between Light Microscopy and V9 Ribotype Counts).

The Marine Ecosystem Biomass Data (MAREDAT) project previously provided global abundance and biomass data for all major planktonic diatoms of the global ocean ecosystem (41). Our dataset showed an overlap of 45 diatom genera with MAREDAT (SI Appendix, Fig. S4 A–C) whereas 34 genera from our study are not present in MAREDAT. A total of 23 genera present in both MAREDAT and the reference database were not found in our dataset. Most of the unmapped genera were either freshwater (e.g., Tabellaria, Ulnaria, Urosolenia) or benthic and marine littoral species (e.g., Amphiprora, Caloneis, Ardissonea, Hyalodiscus, Pseudostriatella, Entomoneis, Phaeodactylum), except for only a few pelagic marine genera (e.g., Bacterosira) (7). Some of these unmapped genera have been reported only in northern latitudes, which may explain their absence in our dataset, which is principally from the Southern Hemisphere (Fig. 1A). A comparison of Bacillariophyta distributions in the OBIS database (20) similarly revealed little overlap because of the lack of previous data from the locations sampled during the Tara Oceans expedition (SI Appendix, Fig. S4D).

Intragenus diversity was found to vary from as low as one ribotype per genus (e.g., Nanofrustulum, Asteroplanus, Bellerochea) to as high as 6,094 ribotypes (Chaetoceros) (Fig. 3, columns a and b and Dataset S1). Chaetoceros was found to be the most abundant and diverse genus, with 73.3% of the ribotypes (and 59.6% of the sequences) belonging to the subgenus Phaeoceros and the remainders to Hyalochaetae (Dataset S1). Chaetoceros (both subgenera), Thalassiosira, Corethron, and Pseudo-nitzschia accounted for the highest number of OTUs (Fig. 3, column c and Dataset S1). As expected, the 5- to 20-µm-size and 20- to 180-µm-size fractions contained the highest numbers of diatom ribotypes although an unexpectedly high number were also found in the smaller size fractions, belonging to smaller species (e.g., Nanofrustulum, Cyclotella, and Minutocellus) but also to larger species (e.g., Attheya, Ditylum, and Bellerochea) (7), perhaps derived from broken cells, broken fecal pellets, or from gametes. The 180- to 2,000-µm-size fraction contained the lowest number of ribotypes, including from chain-forming diatoms (e.g., Hyalosira, Fragilaria) and epizoic species (e.g., Pseudohimantidium), but also from small cells (e.g., Nanofrustulum), possibly having been ingested by larger organisms or otherwise associated with them or with microplastics, or retained in this fraction because of net clogging. A clear distinction was seen in the distribution among different size fractions: e.g., small and mainly solitary Minidiscus, Attheya, and Minutocellus were found highly restricted to the smallest size fractions whereas larger, chain-forming Asterionellopsis, Lauderia, and Odontella were found principally in the 20- to 180-µm-size fractions (Fig. 3, column d).

Different genera were also found to prefer different depths, such as Actinoptychus, Corethron, Coscinodiscus, Fragilariopsis, Leptocylindrus, and Rhizosolenia in subsurface (SRF) samples, whereas Asterionellopsis, Bellerochea, Helicotheca, Nanofrustulum, and Lithodesmium were seen mostly in deep chlorophyll maximum (DCM) samples (Fig. 3, column e). The level of percentage identity to the reference sequence also varied across genera (Fig. 3, column f). Pseudo-nitzschia, Actinocyclus, Attheya, Chaetoceros, Eucampia, Fragilariopsis, Minutocellus, and Thalassiosira were among the most cosmopolitan genera whereas many others (mainly benthic and freshwater genera) were restricted to only a few stations (Fig. 3, column g and Dataset S1).

We performed manual annotation on the top unassigned sequences (representing ∼87% of the unassigned reads) and compared GenBank annotations with those in the PR2 database, which resulted in our being able to assign an additional 13 ribotypes (representing ∼8% of the unassigned reads) from the 113 most abundant sequences to genus or species level. The best assignments and percent identity of these sequences to those present in the reference databases are shown in Dataset S2. Overall the ribotypes that could not be unambiguously assigned to any diatom genus but could be classified only as araphid or raphid pennate, polar, or radial centric, or unassigned diatom on the basis of V9 rDNA annotation (Fig. 3) represented between 31% and 81% of the total number of unique diatom ribotypes at different sampling stations (SI Appendix, Fig. S5). The best assignments and percent identity of these sequences to those present in the reference database are shown in Dataset S2. In general, unassigned ribotypes were particularly common in the SPO, where most of the stations are in the high nutrient low chlorophyll (HNLC) region downstream of the equatorial and Peruvian upwellings, in the IO, and in the warm and salty RS, with almost similar percentages at both depths. The diatoms in the smallest size fraction contributed most to the unknown sequences, with depth having no significant impact (Fig. 4A). On the other hand, the larger size fractions (20–180 µm and 180–2,000 µm) contained the lowest percentage of unassigned ribotypes, consistent with microplanktonic diatoms being the most common and the best studied. The number of unassigned sequences also varied among sampling sites, with the MS, the Benguela upwelling (station 67) (SI Appendix, Fig. S5), and the SO containing the best characterized diatom communities (Fig. 4B).

To investigate whether V9-based relative abundance estimates for diatoms are comparable with community composition studies based on classical morphological identification methods using light microscopy (LM), diatom counts were compared between the two methods for 15 sampling stations. A simple comparison was initially disappointing; however, the correlation between the two kinds of data was significantly improved when “unassigned” and “not known” sequences were removed from the V9 dataset and when some specific adjustments were applied (Materials and Methods) (Fig. 5). A few cases of mismatch still persisted: e.g., the surface sample from station 84 was dominated only by Fragilariopsis sp. in LM counts whereas Chaetoceros (Phaeoceros and Hyalochaetae) and Fragilariopsis were equally dominant genera along with unknown centric diatoms in the V9 dataset. However, the overall match between the two datasets was sufficiently close, thus indicating that V9 counts can provide a reliable estimate of diatom relative abundance at the genus level in a given sample.

LM also assisted in samples where we found a high percentage of unknown ribotypes. For instance, station 84 displayed abundant counts of Asteromphalus, a genus for which no sequences are available in the reference database. We also examined samples that contained a large number of V9 sequences that could not be assigned, specifically from stations 122–124 (SI Appendix, Fig. S5). In these samples, we typically observed a large number of pennate diatoms that could not be identified easily, and so we speculated that many of these unassigned sequences could be from pennate diatoms that do not yet have sequence representation in the V9 dataset. Conversely, centric genera identified by LM but not present in the V9 dataset included Asterolampra, Asteromphalus, Climacodium, Dactyliosolen, Hemiaulus, Hemidiscus, and Lauderia.

We next examined intragenus diversity (expressed as exponentiated Shannon Diversity Index) (42) and distribution in different oceanic contexts for the 20 most abundant genera. Of these abundant genera, we found that Pseudo-nitzschia, Chaetoceros (both subgenera), and Thalassiosira were the most diverse genera whereas Corethron, Leptocylindrus, Minidiscus, and Planktoniella were among the least diverse and that this observation also reflected the known differences in species richness for these genera (Fig. 6A). Most diatom genera were seen in all oceanic provinces although their abundance patterns were highly variable: for instance, Chaetoceros (both subgenera), Corethron, and Fragilariopsis were highly abundant in the SO, in accordance with previous data (e.g., ref. 43); Attheya, Planktoniella, and Haslea were seen principally in the SPO; and Leptocylindrus was found to be highly abundant in the MS, especially at station 11, in line with reports from other Mediterranean sites (44). In terms of global biogeography, the diversity of each genus (expressed as the number of ribotypes) was found to be strikingly variable across the oceans (Fig. 6B and SI Appendix, Fig. S6). Three main patterns were found, with some genera having a lower diversity in the tropics (e.g., Fragilariopsis, Proboscia, and Eucampia), others showing lower diversity at high latitudes (e.g., Attheya, Guinardia), and others with a more uniform diversity (e.g., Thalassosira, possibly the most global diatom genus in our dataset). The two Chaetoceros subgenera showed similar distributions, with higher abundance in the SO (SI Appendix, Fig. S6 B and C) and high richness in coastal and open-ocean areas (SI Appendix, Fig. S6D). The subgenus Phaeoceros was more represented in the larger size fractions at almost all sites, including the offshore Atlantic, Pacific, and Mediterranean waters (Dataset S1 and SI Appendix, Fig. S6C).

Among surface samples, diversity (expressed as exponentiated Shannon Diversity Index) and evenness across oceanic provinces varied greatly, attaining the highest values in the RS, whereas, among the DCM samples, the IO showed the highest diversity; the SO was the least diverse at both depths (Fig. 7A). In terms of richness, the SO stations consistently showed the highest values owing to the presence of a majority of very low abundant ribotypes. Considerable variation in terms of overall ribotype diversity in different size fractions was also observed (SI Appendix, Fig. S7A). In contrast with what was observed globally for marine planktonic eukaryotes in the Tara Oceans dataset (36), diatom diversity did not consistently decrease with increasing size (SI Appendix, Fig. S7A). There were also no discernible differences in diatom diversity patterns between SRF and DCM samples.

Generally, the western boundary currents of the oceanic basins were the most diverse regions. Furthermore, a sudden drop in diversity was observed in the Agulhas retroflection region between the IO (station 65) and the SAO (stations 66/67/68), and from the SAO (stations 76 and 78) to the SO (stations 82/84/85) (Fig. 7B and SI Appendix, Fig. S7B). Diversity was significantly lower in the samples from the Maldives (station 45, North IO) but increased toward the north and the south (Fig. 7B and SI Appendix, Fig. S7B). Station 11 in the MS displayed the lowest diversity of all, the result of a diatom bloom that was dominated by Leptocylindrus (Figs. 1C, 6B, and 7B). In general, although the standardized abundance of diatoms showed a significant decrease from coastal to open ocean (e.g., from stations 65–67 to stations 68–78) and from surface to DCM, with the exception of the Northern IO and the SPO (Fig. 2C), we found no significant difference in the diversity at open ocean stations versus coastal stations (Fig. 7C). Indeed, diversity showed no correlation with diatom V9 sequence abundance.

We then examined whether diatom diversity follows a latitudinal gradient, as has been observed for other marine organisms (45–49). We indeed found a poleward decrease, although the trend was weak (Fig. 7D), most likely because of the lack of data from 50° to 60° latitudes. Analysis of the complete set of data from Tara Oceans will be required before drawing any concrete conclusions about latitudinal gradients.

Diatom-annotated ribotype distribution patterns were generally consistent across all of the stations, in that only a few ribotypes were abundant and the large majority of the richness was contributed by rare ribotypes (Fig. 8). The number of different ribotypes per station varied from as low as 46 (station 48; IO) to as high as 16,100 (station 85; SO), with a mean richness of 4,927. In general, it was found that the more abundant a ribotype, the more ubiquitous was its distribution (Fig. 8). Several ribotypes with considerable abundance but low occupancy were also seen, possibly indicating endemism or a marked seasonality in their occurrence (blooming species). One of the Leptocylindrus ribotypes was one such example. Only 23 ribotypes were found in ≥90% of the studied sites; however, they represented nearly 24% of the total relative abundance. The majority of these cosmopolitan ribotypes could not be assigned to a known diatom taxon (Fig. 8, Lower). A few selected unassigned ribotypes [marked with an asterisk in Fig. 8, Lower] were identified as Shionodiscus bioculatus (“*4”), Asteromphalus spp. (“*11”), Pseudo-nitzschia delicatissima (“*19”) and Thalassiothrix longissima (“*”) (SI Appendix, SI Materials and Methods). Most ribotypes with intermediate abundance aligned along a line (roughly going from occupancy: 25, evenness: 0 to occupancy: 44, evenness: 0.8), indicating a general tendency toward cosmopolitanism that is directly proportional to a deviation from an opportunistic r-strategy (corresponding to a low evenness) (50–52). Furthermore, the wide set of combinations of evenness and occupancy suggests that diatoms actually occupy all kinds of niches (Discussion).

The total number of ribotypes seen in the MS, RS, IO, SAO, SO, and SPO were 13,119, 4,586, 23,722, 16,269, 26,846, and 29,203, respectively. Most of the ribotypes in the SO (53.3%), SPO (33.7%), and MS (26.9%) were not found elsewhere whereas only a few ribotypes were specific to the RS (2.3%). Similarly, the IO (14.2%) and SAO (10.4%), which are transitional basins between the SPO and NAO, showed only a small number of ribotypes endemic to them (SI Appendix, Fig. S7 C and D). Altogether, nearly 52% (32,850 out of 63,371) of the ribotypes were seen only in one province. Interestingly, a substantial number of ribotypes were shared between two provinces [in particular, the SPO and IO (12,176 ribotypes), where the latter is downstream of the former; the SAO and SPO (9,501 ribotypes), mostly because of the coastal SAO stations; the SAO/IO (8,569 ribotypes); and the SO/IO (7,330 ribotypes)] whereas only 576 ribotypes (out of 63,371; 0.9%) were present in all oceanic provinces (SI Appendix, Fig. S7D). Diatoms thus seem to have a significant association to each oceanic basin or to basins that are physically connected (e.g., the SPO and IO via the Indonesian Passage).

The complex biogeographical patterns become clearer when considering the similarity among surface stations for which all four size fractions were available (37 stations). Stations in the SPO, SO, and MS showed the highest degree of internal similarity (Fig. 9), coherent with their relative homogeneity of conditions (for instance, the actual SPO subset is made up of tropical stations in an HNLC, iron-limited tropical region) and geographical isolation (the SO and MS). The clustering of stations revealed four major groups, including one for the MS (the most isolated case), one for the SPO, and another containing oligotrophic, seasonally stable stations where diatoms were present only at low abundance. Finally, the polar SO stations and the rather coastal, mostly temperate stations around South Africa form a fourth cluster, despite their large distance and, in some cases, huge environmental gradients. This latter observation confirms that the Agulhas region, the choke point of the global circulation, is a region of intense mixing among water masses. With the exception of the low abundance clade, these clades shared a considerable percentage (∼20–37%) of ribotypes within them. The community in the MS, a semienclosed basin, was most distinct from the others whereas the IO, the hub of the global surface circulation, showed the highest similarity with the others (Fig. 9 and SI Appendix, Figs. S7C and S8). The SPO and MS stations were nonetheless each seen to cluster together without any overlapping with each other, and the SO stations showed a very distinct community structure (SI Appendix, Fig. S8). Several specific cases illuminate the limits of this simple geographical approach and need to invoke ecological mechanisms to explain the observed patterns. For instance, station 30 in the Eastern Mediterranean Sea is part of the MS cluster even though it is in a phosphate-limited ultraoligotrophic region, unlike all of the other MS stations. Conversely, the Marquesas Islands (stations 122 and 123; SPO) are clearly under the influence of the far upstream Peruvian Upwelling (stations 100 and 102; SPO) whereas, because of the effect of a natural fertilization (53), they are quite different from the very close-by downstream stations 124 and 125. The equatorial upwelling in turn acts as a barrier, making the station further north (station 109) quite different from the others. This latter station is also upstream of all of the others (except the SO) and in fact is similar to most other stations.

For the present study, 293 global samples encompassing 46 stations from the photic zone [subsurface (SRF) and deep chlorophyll maximum (DCM)] were used that corresponded to four size classes (0.8–5 μm, 5–20 μm, 20–180 μm, and 180–2,000 μm). A total of 63,371 V9 rDNA diatom-assigned ribotypes (represented by ∼12.4 million reads) were retrieved from the 293 communities. Please see de Vargas et al. (36) for details on the sequencing and taxonomic assignation of the V9 sequences used in this study.

Metabarcodes were clustered based on their taxonomic affiliation at the level of genus and were organized under 86 genera. Five additional unassigned classes (unassigned, unassigned polar centric, unassigned radial centric, unassigned raphid pennate, and unassigned araphid pennate) were defined to accommodate those reference sequences (n = 416) for which genus assignment was not available. Genus distribution and diversity were assessed for most represented genera.

Deviations from Preston’s log-normal distribution were used to estimate the completeness of richness sampled. Also, the information from the samples was used to extrapolate the number of ribotypes that might be found if sampling were more intensive. The relation between abundance, occurrence, and evenness of each ribotype was assessed. Pielou's evenness (75) and the exponentiated Shannon–Weiner H′ diversity index (42) were used as estimates of diversity. The percentage of shared ribotypes was calculated for each pair of stations, and a Spearman correlation was used as a distance measure to cluster stations. Compositional similarity between stations was computed based on a Hellinger-transformed abundance matrix and incidence matrix using Bray–Curtis and Jaccard indices, respectively, as a measure of β-diversity. Nonmetric multidimensional scaling was performed to visualize the level of similarity between different stations. For all statistical analyses, a value of P < 0.05 was considered significant. All of the data analyses were performed in R (76).