Biodiversity is the product of millions of years of evolution and forms the basis of the earth's life support system, but the magnitude and relative diversity of global species richness remain unknown. On earth, there may be over 100 million species1 but fewer than two million have been formally described2. There is also a pronounced bias towards the study of larger organisms, leaving the most speciose communities that are dominated by microscopic organisms understudied. To study diverse environments dominated by small taxa, second-generation sequencing has been used for the quantification of bacteria, archaea3 4 and viruses5 6, but large knowledge gaps still exist regarding the organization of diversity within several eukaryotic kingdoms, including the Metazoa7. The Kingdom Metazoa, also known as Animalia, consists of multicellular heterotrophic organisms ranging from Porifera to Chordata. Contemporary phylogenetic studies routinely recover a monophyletic Bilateria—Triploblasta clade (including deep-branching Ecdysozoa, Lophotrochozoa and Deuterostomia)8 9, but no consensus view exists of the precise relationships between Bilateria, and the basal groups of Porifera, Ctenophora, Placozoa, Cnidaria and Coelenterata8 9. Marine benthic metazoan communities display some of the highest α-diversity on the planet and occupy one of the largest ecosystems on earth, where only 1% of species are estimated to be known10. Benthic meiofauna (small metazoans between 45 and 500 μm in size) comprises members encompassing 60% of animal phyla and represents a major part of marine biodiversity10. Dominated by nematodes and characterized by high abundances (up to 108 individuals per 1 m2) and diversity11, meiofaunal assemblages perform essential roles in marine ecosystem processes, namely, nutrient cycling, secondary production, sediment transport and mineralization12. In this paper, we apply a metagenetic approach (that is, the large-scale analysis of taxon richness through the analysis of homologous genes7) using second-generation sequencing of the 18S nuclear small subunit (nSSU) ribosomal RNA (rRNA) gene to assess simultaneously the relative levels of richness and patterns of diversity of multiple metazoan phyla across an ecological gradient in a temperate benthic ecosystem. The heterogeneous levels of accumulating taxon richness derived from the benchmarked analyses were broadly congruent with those derived from intensive morphological assessments, but MEGABLAST annotation revealed a previously unidentified phylogenetic breadth of microbial metazoan life. Moreover, the finding that the largely predacious turbellarian Platyhelminthes represent a substantial proportion of benthic diversity quantifies their hitherto unrecognized ecological importance in benthic food chains. Annotated metagenetic analyses enable the objective assessment of microbial biodiversity throughout all ecosystems, facilitating understanding of linkages between microbial biodiversity and ecosystem processes. Results Taxon richness estimates nSSU rRNA amplicons were generated from eight benthic samples from the low-tide zone of an estuarine beach near Prestwick on the West coast of Scotland, and from one sample from a beach in Littlehampton in the South of England. The amplicons were processed for sequencing on a Roche 454 FLX platform generating a total of 353,896 sequences that were quality filtered to 305,702 for downstream analysis (Table 1). When performing metagenetic assessments of taxon richness, it is important to cluster taxonomic units in a biologically meaningful manner, as even slight differences in similarity cut-offs and using different algorithms can result in significantly different estimates of richness. The taxon clustering comparisons between ESPRIT13 and OCTUPUS (the Operational Clustering of Taxonomic Units from Parallel UltraSequencing, supplementary software available at http://octupus.sourceforge.net/) (Fig. 1) on the subsampled beach pyrosequencing data show that ESPRIT overestimates Operational Taxonomic Units (OTUs) richness compared with OCTUPUS (between 1.1x–4.4x, over the 90–99% cut-off range). Phylogenetic14 and barcoding15 studies based on the analysis of chain-termination nSSU gene sequences suggest that intraspecific divergence at least in nematode species is low (1–2%). However, the true level of intragenomic and intraspecific variation among rRNA gene repeats is largely unknown, and genome-wide analyses reveal a dominant set of conserved sequences accompanied by rare variant sequences16. Our benchmarking exercise, performed against a reference control data set comprising 41 species17, revealed that the 96% similarity OCTUPUS clustering algorithm with accompanying chimera screening estimated taxon richness that was most closely aligned with species richness. At 96% similarity, OCTUPUS resulted in 37 operational clustered taxonomic units (OCTUs), whereas at 97%, OCTUPUS resulted in 51 OCTUs from the control metagenetic analysis17. Thus, although we may be underestimating richness by at least 10%, we have opted for a more conservative approach of setting a 96% identity OCTU cut-off for all subsequent numerical comparisons. At this cut-off, an OCTU is likely to (at worst) correspond to a group of related species. Following the 96% similarity OCTUPUS clustering strategy, the total number of putatively non-chimeric tag reads and OCTUs was 217,221 and 428, respectively. Before chimera screening, 1013 OCTUs were clustered from the initial quality-screened 305,702 reads. Community richness is closely linked to the environment The peak of standardized OCTU richness for all phyla was within samples 6, 7 and 8 from Prestwick, and cluster analyses indicated clear and fine-scale hierarchical distinctions in OCTU composition within and between the two sites, with clear divergence of samples from Littlehampton (Fig. 2). OCTU richness and sediment grain size were clearly correlated (Spearman's correlation coefficient, n=9, ρ=−0.82, P=0.0108). Dominance of the Nematoda and rise of the Platyhelminthes Plotting phylum richness rank for all nine independent samples (Fig. 3) shows that Nematoda are the most OCTU-rich in all nine samples, with Platyhelminthes and Arthropoda ranking second and third in all, respectively. Annotated metagenetic analyses can yield robust relative richness estimates and here it was possible to assign 374 OCTUs to phylum (Fig. 4a). The PCR primers used are not fully specific to Metazoa, and thus we also sequenced nSSU genes from protist taxa from the Alveolata (2 distinct OCTUs), Cercozoa (3 OCTUs) and Stramenopiles (15 OCTUs). Of the metazoan OCTUs, 182 were from Nematoda, at least three times more than from any other individual meiofaunal taxon (Fig. 4a). Platyhelminthes (61 OCTUs) was the second richest phylum, followed by Arthropoda (29 OCTUs including Copepoda, Ostracoda and Malacostraca), Mollusca (22 OCTUs), Gastrotricha (7 OCTUs), Annelida (6 OCTUs) and five less-rich phyla (for example, Bryozoa, Echinodermata, Cercozoa, Rotifera and Alveolata with between 1 and 3 OCTUs each). Metagenetics reveals phylogenetically distinct taxa According to the comparisons of the OCTU sequences with the NCBI database, the majority (95%) of Nematoda OCTUs have never been sequenced before (Fig. 4a, Table 2). Similarly, only 6.5% of Platyhelminthes and 17.2% of Arthropoda OCTUs had 100% identity to previously sequenced specimens. The Annelida, Mollusca and Stramenopiles, however, had 50, 23 and 26.6%, respectively, of their OCTUs, with a 100% identity to previously sequenced taxa. None of the Gastrotricha OCTUs were identical to previously sequenced taxa. Overall, 54 OCTUs (representing 7,247 individual sequences) had identities below 90% to any reference nSSU sequence (OCTUs outfile) against the downloaded GenBank/EMBL/DDBJ nucleotide database and taxonomic annotation was restricted to matches of 90% and higher. Acknowledging concerns regarding the misinterpretation of levels of richness due to the formation of recombinant DNA molecules in environmental DNA sequencing26 57, we adopted an aggressive putative chimera culling regime embedded within the OCTUPUS pipeline for the primary data set7. The OCTUPUS chimera screening was followed by manual removal of putative false positive OCTUs (that is, those that exhibited more than 10 base pair mismatches with the MEGABLAST reference sequence) and retrieval of clear false-negative chimeras (that is, those that exhibited 100% length matches with already sequenced taxa). To compare independent estimates of taxon richness, the OCTUPUS clustering algorithm was also tested against the HCluster algorithm of ESPRIT13 on a subsample of the larger data set. Finally, OCTU richness generated from the above OCTUPUS algorithm was benchmarked at a range of percentage similarity cutoff against a reference data set comprising the metagenetic analysis of a phylogenetically diverse combination of 41 nematode17 species. Scripts for preprocessing the above data are available from the Natural Environment Research Council's Environmental Bioinformatics Centre Script Repository (http://nebc.nerc.ac.uk/tools/scripts/general-bioinformatics). For the 54 OCTUs exhibiting BLAST hits below 90% identity, manual MEGABLAST searches were conducted and phylogenetic affinities investigated by means of the NCBI taxonomy browser. Sequence data and all associated fusion primer codes have been deposited in the Entrez short read archive under accession number SRA009394. For direct ecological comparisons of between-sample OCTU richness, the original data set was reanalysed using 15,000 randomly picked sequences (over 200 bases in length, n=135,000) from each sample, before OCTUPUS clustering and annotation. The total number of OCTUs generated from the original and standardized data sets was significantly correlated (Spearman's coefficient: ρ=0.783, P=0.0132); nevertheless, interpretations derived from direct comparisons of richness between samples refer to the standardized data set. Sample cluster analyses (UPGMA) were performed using the Multivariate Statistical Package58 using Sorensen's Coefficient on a binary (presence/absence of OCTU) data matrix. Phylum-specific rarefaction curves for the Prestwick transect were generated using EstimateS (Version 8.2.0, R.K. Colwell) using a range of estimators (for example, ACE, Chao1, Jackknife1 and Bootstrap) that yielded very similar results. The ACE abundance-based coverage estimator59 was used because it represents a consensus view of the analyses and has proven to work well for the analysis of metagenetic data sets3. Author contributions S.C., V.G.F., G.R.C. and W.K.T. conceived and designed the experiments. V.G.F., H.F.J., S.P.N., M.P. and S.C. performed the experiments. V.G.F. S.C. and W.S. analysed the data. S.C., G.R.C., M.L.B., M.P., P.J.D.L., W.K.T., S.P.N. and D.M.P. were involved in the critical review of the manuscript. V.G.F., G.R.C., D.M.P., M.L.B. and S.C. wrote the paper. Additional information Accession number: All sequences and associated fusion primer codes have been deposited in the Entrez short read archive under accession number SRA009394. How to cite this article: Fonseca, V.G. et al. Second-generation environmental sequencing unmasks marine metazoan biodiversity. Nat. Commun. 1:98 doi: 10.1038/ncomms1095 (2010). Supplementary Material Supplementary Software OCTUPUS Supplementary Software
