Structure, dynamics and predicted functional role of the gut microbiota of the blue (Haliotis fulgens) and yellow (H. corrugata) abalone from Baja California Sur, Mexico

Sample collection and genetic analyses

Wild abalone (n = 31 HF, n = 35 HC) were collected from the commercial harvest along the Pacific coast Baja California (Mexico) from Isla Natividad (27°53′33″N 115° 13′19″W) to Punta Abreojos (26°42′52″N 113°34′20″W), during two field expeditions organized in April and November 2012 (Table S1). In the field, approximately 30 mg of post esophageal tissue were dissected from visually healthy animals bearing no signs of the withering syndrome (Friedman, 2012), and immediately transferred to sterile 1.5 ml microcentrifuge tubes containing molecular grade ethanol, until further analysis. Total DNA was extracted and purified from preserved tissues using DNeasy blood & tissue kit (Qiagen, Valencia, CA, USA) following manufacturer’s protocols.

A fragment of the bacterial ribosomal 16S rRNA spanning V1–V3 regions was PCR amplified using universal eubacterial primers 28F: 5′-GAGTTTGATCNTGGCTCAG- 3′ (Ludwig, Mittenhuber & Friedrich, 1993) and 519R: 5′-GTNTTACNGCGGCKGCTG - 3′ (Ruff-Roberts, Kuenen & Ward, 1994). PCR reactions (20 µl) contained: 1X PCR Buffer (Kapa Biosystems, Woburn, MA, USA), 1.5 mM magnesium chloride (Kapa Biosystems, Woburn, MA, USA), 0.2 mM dNTPs (New England Biolabs, Beverly, MA, USA), 0.5 µM each primer, 0.4 mM bovine serum albumin (New England Biolabs, Beverly, MA, USA), 1U Taq polymerase (Kapa Biosystems, Woburn, MA, USA), and 100 ng purified DNA. Thermal cycling consisted of an initial incubation at 94 °C for 4 min, followed by 40 cycles of: 94 °C for 1 min; 62 °C for 30 s and 72 °C for 30 s, and a final incubation of 8 min at 72 °C. Confirmation of amplification was carried out by 1.5% agarose gel electrophoresis. Amplicons were subsequently 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. (2008). Following normalization, Roche 454 pyrosequencing was carried out in a GS FLX Titanium platform by Research and Testing Laboratory (Lubbock, TX, USA).

Finally, given that different DNA extraction and/or purification methods may introduce bias against either gram positive or negative bacteria in purified DNA (Yuan et al., 2012; Gill et al., 2016), we evaluated the introduction of such bias in our samples. We tested in a random subset of eight DNA samples (four from each species) their suitability for amplifying gram negative and gram positive DNA using gram-specific forward primers (16S_68d for gram negative and 16S_143 for gram positive) (Klausegger et al., 1999). PCR reactions (15 µl) contained: 1×  PCR Buffer (Kapa Biosystems, Woburn, MA, USA), 1 mM magnesium chloride (Kapa Biosystems, Woburn, MA, USA), 0.4 mM dNTPs (New England Biolabs), 0.4 µM each primer, 1U Taq polymerase (Kapa Biosystems, Woburn, MA, USA), and 90 ng purified DNA. The same thermal cycling conditions were used as described for universal eubacterial primers, except that the annealing temperatures were 55 and 60 °C for gram positive and negative specific primers, respectively. Confirmation of amplification was carried out by 1.5% agarose gel electrophoresis.

Bioinformatic analyses

The 16S rRNA reads were analyzed using Quantitative Insights Into Microbial Ecology (QIIME) software version 1.9.1 (Caporaso et al., 2010b). The first step consisted in demultiplexing; subsequently, reads were filtered according to Phred quality scores obtained from the 454 pyrosequencing. Acceptance quality criteria consisted of: (i) minimum and maximum lengths of 250 and 550 bp, respectively; (ii) default minimum quality Phred score of 25; (iii) maximum homopolymer length of 6 bp.

Sequences that met quality criteria were initially filtered using USEARCH (Edgar, 2010) to detect and remove chimeras as implemented in QIIME pipeline (Caporaso et al., 2010b) and then clustered in operational taxonomic units (OTUs) at 97% sequence similarity using the UCLUST algorithm (Edgar, 2010). Taxonomic assignment was carried out using the longest sequence from each OTU and SILVA 128 database (Released: 28.09.2016 containing 646,151 sequences; http://www.arb-silva.de/). OTUs with total read number <20 were discarded; remaining OTUs were then aligned using Python Nearest Alignment Space Termination (PyNAST) algorithm (Caporaso et al., 2010a). The alignment was filtered to remove gaps and outliers (e.g., sequences dissimilar to the alignment consensus) using default QIIME parameters.

Additionally, to evaluate the potential bacterial origin of DNA of Unassigned OTUs (with no match in SILVA database), we aligned them with bona fide representative 16S rDNA sequences of Eubacteria, chloroplasts and Archean species from GenBank and reconstructed a Neighbor-joining (NJ) tree using MEGA 6 (Tamura et al., 2013). OTUs that fell outside of the Eubacterial lineage were removed from the PyNAST alignment before further analyses.

The number of reads of each abalone was normalized using the rarefication method (Weiss et al., 2017). Normalization was achieved by randomly subsampling 500 reads per abalone; we selected a sample size as close as possible to the asymptotic plateau of rarefaction curves in most abalone (Fig. S1).

Principal coordinate analysis (PCoA) based on unweighted unique fraction metric (or UniFrac) distance was performed to compare the microbiota among HF and HC abalone species from an ecological-phylogenetic framework. Results were visualized using EMPeror (Vázquez-Baeza et al., 2013). Briefly, unweighted UniFrac is a qualitative β diversity analysis measuring the distance between two or more communities, as the fraction of the branch length in a phylogenetic tree that leads to descendants from either one community, but not both (Lozupone et al., 2007; Lozupone & Knight, 2005). Accordingly, this measure may reflect essential microbial adaptations to one environment (Lozupone & Knight, 2005). Statistical evaluation of UniFrac results was carried out using ANOSIM as implemented in QIIME (Caporaso et al., 2010b).

Ecological analysis

The normalized number of reads was used as abundance proxy to estimate diversity and structure of the abalone gut microbiota. Differences in the bacterial composition between abalone species were evaluated by a linear discriminant analysis (LDA) effect size (LEfSe) (Segata et al., 2011).

To assess how exhaustively bacterial communities of both abalone species were sampled, rarefaction curves of discovered OTUs were generated for increasing numbers of sampled abalone (Fig. 1). Also, OTU abundance was used to compute the non-parametric species richness estimator Chao 1 (Chao, 1984). Rarefaction curves were obtained using the EstimateS V.9.0.1 program (Colwell, 2013).

Microbiome community structure was evaluated using non-parametric multidimensional scaling (MDS) analyses using Bray–Curtis and Sorensen similarity indices based on read abundance and on presence/absence, respectively, as implemented in PRIMER V.6 (Clarke & Warwick, 2001). The statistical comparison of MDS results was performed with ANOSIM as implemented in Past V. 2.17c (Hammer, Harper & Ryan, 2001) using OTUs assigned to bacterial taxa matching SILVA database. To determine which OTUs were primarily responsible for the dissimilarity between HF and HC, a SIMPER analysis (Table S2) with square root transformed data was performed using PRIMER V.6 (Clarke & Warwick, 2001). In order to increase the number of OTUs used in the characterization of microbiomes, we repeated the analyses also including SILVA unassigned OTUs that were phylogenetically identified as bacteria, as described above.

Bacterial interactions in the microbiomes of both abalone species were estimated using Jaccard distance (J_d) as implemented in PRIMER V.6 (Clarke & Warwick, 2001). J_d is a measure of dissimilarity for all pairwise combinations of a data set (Dowd et al., 2008) and was calculated using OTUs presence/absence. J_d values close to 0 (from 0 to 0.33) were interpreted as co-occurrence (or putative mutualistic relationships) and values close to 1 (from 0.68 to 1) as interactions leading to exclusion (or putative competitive) (Rahel, 2000); whereas intermediate values were considered as neutral relationships.

Functional prediction of predominant bacterial species

A phylogenetic investigation of communities by reconstruction of unobserved states (PICRUSt) (Langille et al., 2013) was carried out to predict the functional attributes of metabolic genes from HF and HC GI microbiotas. Briefly, PICRUSt is a bioinformatic approach that uses information from a genetic marker, such as the 16S rRNA gene, to predict the functional content of a bacterial community characterized by metabarcoding (Langille et al., 2013; De Voogd et al., 2015). Functional predictions are obtained by matching 16S rRNA gene sequences against a genomic KEGG database, previous normalization of read numbers (Langille et al., 2013; De Voogd et al., 2015). The central result of PICRUSt consists of a table of the functional gene frequencies known as KEGG Orthologs (or KOs). KOs are hierarchically organized in sets of homologous sequences with known function and assigned to biological pathways. We analyzed the data using the raw KOs counts as well as categorizing them by biological pathway. PICRUSt analyses use Greengenes (last released in 2013 containing 202,421 bacterial and archaeal sequences, http://greengenes.lbl.gov) as taxonomic and functional reference database. To implement quality control, we computed weighted nearest sequenced taxon index (NSTI) values for the metabarcodes of each abalone. NSTI was developed to evaluate the prediction accuracy of PICRUSt, since it reflects the average genetic distance (measured as number of substitutions per site) between an OTU against a reference genome (Langille et al., 2013; De Voogd et al., 2015). Following suggested guidelines (Langille et al., 2013), we eliminated observations with a NSTI higher than 0.17.

In order to assess the contribution of individual OTUs to predicted KO functions, first we focused our analyses on genes involved in metabolic pathways (KEGG IDs from EC:1.1.1 to EC:6.5.1). Next, we categorized the relative importance of KO genes by ranking them according to their counts. These counts were obtained with PICRUSt v.1.1.0 (Langille et al., 2013) and were log transformed to rank their relative abundance (Urbanová & Bárta, 2014). We focused our attention on a chosen subset of KOs with the highest counts (n = 10 selected at random), as a first order analysis to characterize functional differences between these microbiotas. Significance of differences in the contribution of OTUs to these KOs was evaluated non-parametrically using a χ² test.

Finally, to compare the functionality of individual microbiomes in both species, non-parametric MDS analyses based on Bray–Curtis similarity using fourth root transformation were carried out using the count number of identified KEGG orthologs and the count number of predicted metabolic functions in PRIMER V6 (Clarke & Warwick, 2001).

Pyrosequencing and metabarcoding results

Pyrosequencing yielded 451,095 raw 16S rRNA reads of which 239,125 met quality criteria and were non-chimeric, these were assigned to 1,508 OTUs. A number of these OTUs were removed from subsequent analyses, 1,066 because their read number was <20 and an additional 128 because they were genetically dissimilar from bona fide bacterial 16S rRNA sequences. The remaining 314 OTUs (comprising 157,686 quality reads −66%-) included 58 Unassigned-OTUs (Table S3) whose phylogenetic position fell within the Eubacterial linage (Fig. S2), suggesting they originate from unknown or poorly characterized bacterial taxa.

The most diverse microbial phyla in the gut microbiota of both abalone species were (number and % of OTUs): Tenericutes (230, 73.2%), Fusobacteria (16, 5.1%), and Proteobacteria (9, 2.9%). The most abundant OTUs were assigned to class Mollicutes (order Mycoplasmatales). Notably, LEfSe analysis revealed that 130 of 230 Mollicutes OTUs were exclusive or predominant to either in HC or HF (Table S4). Similar results were observed for other predominant bacterial families such Vibrionaceae (Table S4).

The classes Mollicutes, Fusobacteria, Alfa and Gamma-protobacteira comprised 99% of the identifiable reads. Rarefaction curves suggest that the bacterial communities were sufficiently sampled in both abalone species, given their asymptotic shape and the proximity of the observed number of taxa found in each species to CHAO 1 estimates (Fig. 1). Normalization was achieved by standardizing to a subsampling depth of 500 randomly subsampled reads per abalone, which represents the number of reads close to the asymptote for most organisms (Fig. S1). As a result, five HC and two HF abalone with less than 500 reads were excluded from subsequent analyses.

Despite the intra-specific compositional variation (Fig. S3) and similar community structure at the highest taxonomic resolution between abalone species (Fig. 2), MDS analyses performed using only bacterial OTUs, assigned using SILVA database, revealed a clear-cut structural difference between the GI bacterial metabarcodes of both species, both considering read number as a proxy of abundance and solely on the basis of presence/absence data (Fig. 3). Significant interspecific differentiation was corroborated by ANOSIM analyses in both cases (R_abundance = 0.786; p < 0.001; R_{presence∕absence} = 0.788; p < 0.001). Given the large fraction of unassigned OTUs (36%), the inclusion of originally unassigned OTUs but phylogenetically identified as bacteria produced similar results (Fig. S4; R_abundance = 0.792; p < 0.001; R_{presence∕absence} = 0.789; p < 0.001). This suggests that the distinction of both abalone microbiomes can be attributed to known and unknown bacterial taxa. PCoA based on unweighted UniFrac distance also supported the structural separation among HC and HF bacterial communities (Fig. S5; ANOSIM, R_{unweighted UniFrac} = 0.39; p < 0.001).

Jaccard distances revealed that interspecific relationships among OTUs changed by an order magnitude with most involving competition (HC: 23,588 and HF: 30,129), followed by neutral (HC: 1851 and HF: 1879) and a smaller number of mutualistic interactions (HC: 212 and HF: 123). Furthermore, the ${\bar{J}}_d$ of the majority of OTUs decreased significantly with increasing read number in both species (p < 0.001; Fig. 4).

Our use of specific primers corroborated the presence of amplifiable DNA from gram-positive and gram-negative bacteria in the extracted DNA, hence we posit that the observed biased composition is not an artifact of the DNA extraction method (Fig. S6).

Functional profiling

Greengenes database, allowed us to detect 248 OTUs, of which 81% (201) were shared with SILVA. Most abalone possessed mean NSTI (Nearest Sequenced Taxon Index) values from 0.07 to 0.17, except for seven HCs with larger values that suggested unreliable functional assignments; hence they were excluded from further analyses. PICRUSt identified 4,092 KOs genes (Table S5) involved in 261 metabolic functions (Table S6). An order of magnitude drop in log-transformed abundance was observed in the ranking of KOs (Fig. S7); hence, metagenome contributions analysis was carried out on a random set of 10 of the 86 most abundant KOs (i.e., Log (KOs counts) > 4). According to PICRUSt, the metabolic functions in the HF microbiomes were generally enriched by one primary OTU, whereas many more OTUs contributed to the same function in HC (Fig. 5). The contributions of bacterial OTUs to each KO differed significantly between species (χ² >  47.74, p < 0.001; Fig. 5). Given their numerical predominance in the metabarcoding reads involved in functional profiling (76%), Mycoplasma occupied a major role in the taxonomic spectra of the KOs analyzed (72% in HC, 90% in HF).

MDS analyses performed using both metabolic function counts (Table S5) or KOs counts (Table S6) revealed no clear functional distinction of the GI microbiota of both species. Nevertheless, the scatter of individual microbiomes is much larger in HC, which is consistent with its higher microbial diversity (Fig. 6).

Abalone microbiome composition

Mollicutes, mostly represented by Mycoplasma spp., was by far the most abundant class, followed by Fusobacteria, Alphaproteobacteria and Gammaproteobacteria. These bacteria have also been found dominating in the GI microbiota of other abalone species (H. discus hannai, H. diversicolor and H. gigantean (Morales-Bojórquez, Muciño Díaz & Vélez-Barajas, 2008; Tanaka et al., 2003; Macery & Coyne, 2004; Rungrassamee et al., 2014)). Even though the taxonomic composition of HC and HF GI microbiotas bears resemblance at high taxonomic levels, the species level composition showed significant differences. Additionally, unweighted UniFrac suggests that observed differences may reflect microbial evolutionary adaptations to either HC or HF. In this context, several Mycoplasma species predominating in one species of abalone were either absent or at low abundance in the other (Table S4). Furthermore, all Vibrio spp. possessed a higher prevalence in HC. Following our interpretation, bacterial abundance may be controlled by interspecific competition. Indeed ${\bar{J}}_d$ values, suggest that for most bacterial species their abundance increases with decreasing ${\bar{J}}_d$, as observed for Mycoplasma. Furthermore, ${\bar{J}}_d$ for a single bacterium changed according to its host, which suggests that the same bacterial species faces different competitive pressures depending on the microbiota composition. In other words, bacteria that may survive competition in one microbiome, may be out competed in the other.

Some pathogenic bacteria, such as Candidatus Xenohaliotis californiensis were observed in the microbiomes of healthy abalone. This bacterium is a pleomorphic, gram-negative coccobacillus that inhabits abalone gastrointestinal epithelia and is considered an obligate endoparasite, like other Rickettsiales (Friedman, 2012; Friedman et al., 2000; Moore et al., 2002; Crosson et al., 2014). We observed the presence of Candidatus Xenohaliotis californiensis in healthy blue and yellow abalone, which supports that the presence of this pathogen is not sufficient to trigger withering syndrome, as previously suggested (Cáceres-Martínez, Vásquez-Yeomans & Flores-Saaib, 2011). Moreover, its absence and/or low intensity in abalone with morphological and histological signs of withering syndrome has already been reported (Cáceres-Martínez, Vásquez-Yeomans & Flores-Saaib, 2011; Balseiro et al., 2006; Horwitz, Mouton & Coyne, 2016). Accordingly, we posit that further investigations are needed to reveal all the factors involved in withering syndrome outbreaks.

Notably, Vibrio halioticoli was absent in both HF and HC whereas it has been found at a prevalence ranging from 40 to 65% in the GI of other marine invertebrates, including several abalone species (Cox, 1962). The discrepancy may in part be related to the anatomical source of the microbiomes, our samples come from post esophageal tissue only whereas other studies have analyzed the entire GI tract (Iehata et al., 2014; Tanaka et al., 2004). Consequently, the presence of V. halioticoli in the rest of GI tissues of HC and HF should be explored.

Given that our extraction method does not incorporate a bead-beating step, it may introduce a bias against gram-positive bacteria (Yuan et al., 2012; Gill et al., 2016). However, our gram-specific amplifications revealed that both groups are present in the DNA extracted from the post-esophagus of both species of abalone producing bands of comparable strength in agarose gel electrophoresis. This suggests a lack of bias introduced by the DNA extraction method; however, since endpoint PCR results are not quantitative, they are not conclusive on the extent to which the observed gram-negative bias in the pyrosequenced reads relates to the extraction method. This question requires further attention. Nevertheless, our main findings concern the differences between the sympatric abalone species and do not rely as much on the accuracy of the microbiome estimation, important as it is, but rather on the precision required to differentiate them. Our data proved to be robust in this regard.

Predicted functions of abalone microbiomes

According to the MDS of KOs gene counts and ecological functions, the inferred functionality of the GI microbiota of HF appears less variable than that of HC, which may reflect a higher degree of diet specialization. The diet diversity of wild HF in Baja California has been shown to be more limited and dominated by Phyllospadix torreyi (47%) and algae in the order Gelidiales (13%). On the other hand, the diet of sympatric wild HC is more diverse consisting of different species of Phaeophyceae (10–20%), Rhodophyta (20%) and Gelidiales (20%) among others (Guzman del Próo, Serviere-Zaragoza & Siqueiros-Beltrones, 2003). These ecological differences are related to their different bathymetric distribution; indeed, HC are generally found in deeper waters (between 10–20 m), whereas HF generally inhabit shallow subtidal regions of rocky shores (between 3–10 m) (Guzman del Próo, Serviere-Zaragoza & Siqueiros-Beltrones, 2003).

PICRUSt analysis, revealed that Mycoplasma contributed to all predicted KEGG in both HC and HF. Heretofore, Mycoplasma has been considered an obligate parasite and/or commensal due its small genome and limited number of genes (Bano et al., 2007). Nevertheless, according to recent genomic annotation analyses carried out in Mycoplasma strains isolated from deep-sea isopods, a mutualistic relationship has been proposed between Mycoplasma and their hosts (Wang et al., 2016). Indeed, the genome of studied Mycoplasma presented a high number of genes involved in degradation of glycans, proteins and complex oligosaccharides, which supports the hypothesis that Mycoplasma may supply their hosts with amino sugars and simple carbohydrates (Wang et al., 2016). Also, given the presence of sialic acid lyase genes, Mycoplasma likely protects their hosts against microbial pathogen infections, breaking down the sialic acid cell-wall “coat” used by a great variety of bacterial pathogens to avoid the host’s innate immune response (Wang et al., 2016; Severi, Hood & Thomas, 2007).

Our findings support a mutualistic association between Mycoplasma and their abalone hosts; indeed, Mycoplasma was the dominant linage in all explored KOs (Fig. 5) which suggests that these bacterial species may play a pivotal role in several metabolic functions. Moreover, a given KO was generally enriched by a single predominant Mycoplasma OTU in HF, whereas in HC it was generally enriched by two or more Mycoplasma as well as other bacterial-OTUs (Fig. 5). Accordingly, Mycoplasma species may be highly host-specific (Register et al., 2015), and they may bear some degree of specificity to particular metabolic functions and/or to specific steps along a metabolic route. Nevertheless, the functional structures of the GI microbiome of both species of abalone bear significant overlap.