Sample collection, DNA extraction, and PCR amplification
Commercially harvested wild abalone were collected in April and November of 2012 during two field expeditions conducted along the Pacific coast of Baja California Sur (Mexico; Additional file 1: Table S1). In the field, landed abalones were morphologically inspected to identify animals bearing morphological signs of WS including low mobility, color alterations, and pedal muscle and mantle retraction, as proposed by Friedman [34]. After the morphological examination, approximately 30 mg of post-esophageal tissue was collected from each of the 107 abalone included in the study (n = 46 for HC; n = 61 for HF). Specimens were immediately transferred to sterile 1.5-ml microcentrifuge tubes containing molecular grade ethanol until further analysis.
In the laboratory, DNA was extracted and purified from preserved tissues using a DNeasy blood & tissue kit (Qiagen, Valencia, CA, USA) following the protocol of the manufacturer. PCR amplification of an internal fragment of the 16S rRNA gene spanning the V1-V3 regions (~500 bp) was conducted using universal eubacterial primers 28F (5’ – GAGTTTGATCNTGGCTCAG – 3’) [35] and 519R (5’ – GTNTTACNGCGGCKGCTG – 3’) [36]. Amplifications were carried out in 20-μl reactions containing 100 ng of DNA, 1X PCR buffer, 1.5 mM MgCl2 (Kapa Biosystems, Woburn, MA, USA), 0.2 mM dNTPs (New England Biolabs, Beverly, MA, USA), 0.3 μM of each primer, and 1U of Taq polymerase (Kapa Biosystems, Woburn, MA, USA). The thermal cycling conditions were 94 °C for 4 min, followed by 40 cycles of 94 °C for 1 min, 62 °C for 30 s, 72 °C for 30 s, and a final extension of 8 min at 72 °C. Confirmation of amplification was carried out by 1.5% agarose gel electrophoresis.
16S rRNA gene library preparation
PCR amplicons were sent to the Research and Testing Laboratory (Lubbock, TX) for pre-sequencing preparation and Roche 454 pyrosequencing. Briefly, the amplicons were tagged using Roche 454 adaptors and multiplex identifier (MID) tags for each organism, following the bacterial tag-encoded FLX amplicon pyrosequencing (bTEFAP) approach of Dowd et al. [37]. Roche 454 pyrosequencing was carried out in a GS FLX Titanium platform. The raw 16S rRNA gene libraries were processed as follows: 1) the SFF files were formatted into FASTA sequences and quality files; 2) the reads were demultiplexed from SFF files using the Roche sffinfo tool; 3) low quality bases (Phred score < 25), primers, and barcode sequences were removed using a Research and Testing Laboratory internal quality trimming algorithm; and 4) dereplication was carried out using the USEARCH algorithm [38]. The raw reads were deposited in the National Center for Biotechnology Information (NCBI; BioProject: PRJNA494699, accession number: SRR7969002).
The returned reads were processed and analyzed using Quantitative Insights Into Microbial Ecology v. 2019.4 – 2019.7 (QIIME2) [39]. After importing into QIIME2, the reads were de novo clustered into OTUs at 97% identity using the QIIME vsearch cluster-features-de-novo plugin [40]. Chimeras and singletons were detected and removed using UCHIME [41] and the filter-features plug in implement in QIIME2 v.2019-7, respectively. We organized two control steps in order to minimize PCR amplicon noise. Initially, the taxonomic assignments of the OTUs were conducted via a classify-consensus-blast [42] trained on 16S rRNA gene OTUs clustered at 99% similarities within the Silva_132 database [43], and unreliable sequences with no match in the SILVA_132 database (e.g., unassigned) were removed before further analyses. Additionally, we removed reads with frequencies inferior to 0.005% [44]. As a result, OTUs with total read numbers < 4 were removed from both the HF and HC data sets.
Finally, the reads were rarefied by randomly subsampling at 702 and 644 reads for HF and HC, respectively. We selected those sampling depths in order to minimize the loss of both specimens and bacterial diversity. The selected rarefaction values were as close as possible to the asymptotic plateau of the rarefaction curves for both species (Additional file 2: Fig. S1a and Fig. S1b), which also allowed for the retention of most abalone. The number of removed reads and OTUs during each quality control step, as well as the final numbers of reads and OTUs used as input in downstream analyses are reported in Additional file 1: Table S2.
Ecological analyses
To evaluate how exhaustively the bacterial communities were sampled, rarefaction curves of the detected OTUs were generated using the diversity alpha-rarefaction plugin implemented in QIIME2 v. 2019.7 [39] for each abalone. Also, the number of OTUs obtained was compared against the non-parametric species richness estimator Chao 1 (Additional file 2: Fig. S1c) [45].
Microbiome community structure was evaluated by principal coordinate analysis (PCoA) based on Bray Curtis, Jaccard, and weighted and unweighted phylogenetic Unifrac distance metrics obtained from the abundance of each rarefied OTU. The PCoA results were visualized with EMPeror [46]. Statistical differences were evaluated by a permutational multivariate analysis of variance (PERMANOVA) based on 4999 permutations using the beta-group-significance plugin of QIIME2 v.2019.7 [47]. Additionally, the multivariate dispersion based on OTU abundance in morphological categories was estimated by a permutation analysis of multivariate dispersion (PERMDISP) with 4999 permutations implemented in PRIMER+P v.6 [48].
At the OTU taxonomic level and with the aim to include all shared bacterial OTUs and/or exclusive non-rarefied OTUs between morphologically healthy and WS-stressed abalone, the OTUs were visualized using Venn diagrams created with open web-software InteractiVenn [49]. Also, the differential abundance of rarefied OTUs between healthy and WS afflicted abalone was determined by linear discriminant analysis (LDA) effect size (LEfSe) [50]. Finally, to determine which OTUs primarily contributed to the dissimilarity between both proposed categories, a similarity percentage (SIMPER) analysis based on rarefied OTU abundance was conducted using PRIMER+P v.6 [48].
Functional prediction
We used PICRUSt v.2.2.0 beta [51] to predict both the potential functional capabilities and the contributions of distinct OTUs to Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. Briefly, the OTUs were normalized by the predicted 16S rRNA gene copy number of each OTU. To implement quality control, we computed the nearest sequenced taxon index (NSTI). Briefly, the NSTI evaluates the prediction accuracy of PICRUSt since it reflects the average genetic distance (measured as the number of substitutions per site) between each OTU against that of the 16S rRNA gene from a reference genome [51–53]. Following the suggested guidelines [51], we eliminated OTUs with NSTI values higher than 2.
To visualize the functional dispersion of morphologically healthy and WS-stressed abalone, PCoA plots based on fourth root normalized metabolic KEGG pathway counts were carried out using STAMP v. 2.1.3 [54]. Significant differences among enriched pathways (effect size > 2) were evaluated by a Welch test implemented in STAMP v. 2.1.3 [54].