The sponge microbiome project

1Sponges (phylum Porifera) are an ancient metazoan clade, with more than 8500 formally described species. Sponges are benthic organisms that have important ecological functions in aquatic habitats. Marine sponges are often found in symbiotic association with microorganisms, and these microbial communities can be very diverse and complex. Sponge symbionts perform a wide range of functional roles, including vitamin synthesis, production of bioactive compounds, and biochemical transformations of nutrients or waste products. The diversity of microorganisms associated with sponges has been the subject of intense study (the search of “sponge microbial diversity” returned 348 publications in the Scopus database). Most of these studies were performed on individual species from restricted geographic regions [e.g., ]. A comparative assessment of these studies is often hindered by differences in sample processing and 16S rRNA gene sequencing. However, 2 recent studies incorporating a large number of sponge microbiomes (>30) revealed the potential of large-scale, standardised, high-throughput sequencing for gaining insights into the diversity and structure of sponge-associated microbial communities. The purpose of this global dataset is to provide a comprehensive 16S rRNA gene-based resource for investigating and comparing microbiomes more generally across the phylum Porifera.

2Sample collection and processing, species identification, and DNA extractions were conducted as previously described. A total of 3569 sponge specimens were collected, representing at least 268 species, including several yet unidentified taxa (hereafter collectively referred to as species) (Supplementary Table S1). Of all species, 213 were represented by at least 3 specimens. Carteriospongia foliascens had the highest replication, comprising 150 individuals. Seawater (n = 370), sediment (n = 65), algae (n = 1), and echinoderm (n = 1) samples as well as biofilm swabs (n = 21) of rock surfaces were collected in close proximity to the sponges for comparative community analysis. Six negative control samples (sterile water) were processed to identify any potential contaminations. Of the samples included in this current dataset, 973 samples had been analysed previously. Samples were collected from a wide range of geographical locations (Fig. 1; Supplementary Table S1). Total DNA was extracted as previously described and used as templates to amplify and sequence the V4 region of the 16S rRNA gene using the standard procedures of the Earth Microbiome Project (EMP).

3Raw sequences were demultiplexed and quality controlled following the recommendations of Bokulich et al.. Quality-filtered, demultiplexed fastq files were processed using the default closed-reference pipeline from QIIME v. 1.9.1 (QIIME, RRID:SCR_008249). Briefly, sequences were matched against the GreenGenes reference database (v. 13_8 clustered at 97% similarity). Sequences that failed to align (e.g., chimeras) were discarded, which resulted in a final number of 300 140 110 sequences. Taxonomy assignments and the phylogenetic tree information were taken from the centroids of the reference sequence clusters contained in the GreenGenes reference database (Greengenes, RRID:SCR_002830). This closed-reference analysis allows for cross-dataset comparisons and direct comparison with the tens of thousands of other samples processed in the EMP and available via the Qiita database.

4Quality-filtered, demultiplexed fastq files were also processed using Mothur v. 1.37.6 (Mothur, RRID:SCR_011947) and Python v. 2.7 (Python Programming Language, RRID:SCR_008394) custom scripts with modifications from previously established protocols. Detailed descriptions and command outputs are available at the project notebook (see Availability of supporting data). Briefly, sequences were quality-trimmed to a maximum length of 100 bp. To minimize computational effort, the dataset was reduced to unique sequences, retaining total sequence counts. Sequences were aligned to the V4 region of the 16S rRNA gene sequences from the SILVA v. 123 database (SILVA, RRID:SCR_006423). Sequences that aligned at the expected positions were kept, and this dataset was again reduced to unique sequences. Further, singletons were removed from the dataset, and the remaining sequences were preclustered if they differed by 1 nucleotide position. Sequences classified as eukaryote, chloroplast, mitochondria, or unknown according to the Greengenes (v. 13_8 clustered at 99% similarity) and SILVA taxonomies were removed. Chimeras were identified with UCHIME (UCHIME, RRID:SCR_008057) and removed. Finally, sequences were de novo clustered into operational taxonomic units (OTUs) using the furthest neighbour method at 97% similarity. Representative sequences of OTUs were retrieved based on the mean distance among the clustered sequences. Consensus taxonomies based on the SILVA, Greengenes, and RDP (v. 14_03 2015; Ribosomal Database Project, RRID:SCR_006633) databases were obtained based on the classification of sequences clustered within each OTU. The inclusion of these taxonomies is helpful considering that they have substantial differences, as recently discussed. For example, Greengenes and RDP have the taxon Poribacteria, a prominent sponge-enriched phylum, which did not exist in the SILVA version used.

5Recently, sub-OTU methods that allow views of the data at single-nucleotide resolution have become available. One such method is Deblur, which is a de-noising algorithm for identification of the actual bacterial sequences present in a sample. Using an upper bound on the polymerase chain reaction and read-error rates, Deblur processes each sample independently and outputs the list of sequences and their frequencies in each sample, enabling single nucleotide resolution. For creating the deblurred biom table, quality-filtered, demultiplexed fasta files were used as input to Deblur using a trim length of 100 and min-reads of 25 (removing sOTUs with <25 reads total in all samples combined). Taxonomy was added to the resulting biom table using QIIME, RDP classifier, and Greengenes v. 13.8.

6For enrichment analysis of metadata terms in a set of sequences, each unique metadata value is tested using both a binomial test and a ranksum test. All analysis is performed on a randomly subsampled (5000 reads/sample) table.

7For a bacterial sequence s and metadata value v, denote N the total number of samples, O(s) the number of samples where s is present, Kv(s) the number of samples with value v where s is present, and T(v) the total number of samples with value v. where PNull(s) = O(s)/N

8The ranksum (frequency aware) P-value is calculated using the Kruskal-Wallis test (implemented in scipy 0.19) as follows.

9For a bacterial sequence s and metadata value v, denote by Fv(s) the vector of relative frequencies of bacteria s in all samples with metadata value v, and denote by the vector of relative frequencies of bacteria s in all samples with metadata other than v. The ranksum P-value is then calculated using the Kruskal-Wallis test for Fv(s) and and shown only if significantly enriched in samples containing v (i.e., rank difference of Fv(s)— > 0).

10We have set up a webserver that performs this enrichment analysis for user-defined sequence submissions. The code for the webserver is also available in Github for a local installation.

11The dataset covers 4033 samples with a total of 1 167 226 701 raw sequence reads. These sequence reads clustered into 39 543 OTUs using QIIME’s closed-reference processing, 518 246 OTUs from de novo clustering using Mothur (not filtered for OTU abundances), and 83 908 sOTUs using Deblur (with a filtering of at least 25 reads total per sOTU). We recommend that data users consider the differences in sequencing depths per sample and abundance filtering for certain downstream analyses, such as when calculating diversity estimates and comparing OTU abundances across samples. In terms of taxonomic diversity, most Mothur OTUs were assigned to the phylum Proteobacteria, although more than 60 different microbial phyla were recovered from the marine sponge samples according to SILVA (n = 63) and Greengenes classifications (n = 72) (Fig. 2).

12This dataset can be utilised to assess a broad range of ecological questions pertaining to host-associated microbial communities generally or to sponge microbiology specifically. These include: (i) the degree of host specificity, (ii) the existence of biogeographic or environmental patterns, (iii) the relation of microbiomes to host phylogeny, (iv) the variability of microbiomes within or between host species, (v) symbiont co-occurrence patterns, and (vi) assessing the existence of a core sponge microbiome. An example of this type of analysis is shown in Fig. 3, where samples were clustered using unweighted UniFrac data with a Principal Coordinates Analysis and visualization in Emperor based on their origins from sponges, seawater, or kelps.

13Project name: The Sponge Microbiome Project

14Project home page: www.spongeemp.com; https://github.com/amnona/SpongeEMP

15Operating system(s): Unix

16Programming language: Python and R

17Other requirements: Python v. 2.7, Biopython v. 1.65, Python 3.5, R v. 3.2.2, Mothur v. 1.37.6, QIIME v. 1.9.1, Deblur

18License: MIT

19Any restrictions to use by non-academics: none

20Raw sequence data were deposited in the European Nucleotide Archive (accession number: ERP020690). Quality-filtered, demultiplexed fastq files, QIIME resulting OTU tables are available at the Qiita database (Study ID: 10 793). The additional datasets that support the results of this article are available in the GigaScience repository, GigaDB and include an OTU abundance matrix (the output “.shared” file from Mothur, which is tab delimited), an OTU taxonomic classification table (tab delimited text file), an OTU representative sequence FASTA file, a table of samples’ metadata, the biom files from QIIME and Deblur analyses, and the QIIME-generated tree file. The project workflow, Mothur commands, and additional scripts are available as HTML in GigaDB.

21The deblurred dataset has also been uploaded to an online server that supplies both html and REST-API access for querying bacterial sequences and obtaining the observed prevalence and enriched metadata categories where the sequence is observed (Figure 4). This allows an interactive view of which sequences are associated with which specific parameters, such as depth or salinity.

22sample.metadata

23EMP: Earth microbiome project; bp: base pairs; OTU: operational taxonomic unit; rRNA: ribosomal RNA.

24T.T. and N.S.W were funded by Australian Research Council Future Fellowships FT140100197 and FT120100480, respectively. T.T. received funds from the Gordon and Betty Moore Foundation. This work was also supported in part by the W.M. Keck Foundation and the John Templeton Foundation. R.K. received funding as a Howard Hughes Medical Institute Early Career Scientist.

25The authors declare that they have no competing interests.

26L.M.-S., N.S.W., and T.T. designed the study. C.A.G., D.S., F.L., G.S., G.K., G.McC., G.-F. F, J.J.B., J.V., J.R.B., J.M.M., J.R., L.S., M.C.P, M.V.M., M.W.T., N.S.W., P.P., P.M.E., P.J.S., R.L.S, R.W.T., R.C., R.T.H., S.L-L., T.D., T.R., U.H., and Z-Y. L. collected samples. C.A.G., D.S., J.V., J.R.B., L.S., M.C.P., M.W.T., N.S.W., P.M.E., R.L.S, R.W.T., S.L-L., and U.H. extracted DNA. G.L.A. and R.K. sequenced DNA. L.M.-S., S.N., A.A., A.G., G.L.A., and T.T. performed data processing and analysis. L.M.-S., N.S.W., and T.T. wrote the manuscript. All authors contributed to the writing of the manuscript.