Microbiome differences between river-dwelling and cave-adapted populations of the fish Astyanax mexicanus (De Filippi, 1853)

Study area

We included one cavefish population (Pachón) from the Sierra de El Abra, a mountain range running north to south along the eastern margin of the Sierra Madre Oriental, a karstic region in northeastern Mexico (Elliott, 2016). El Abra is considered as one of the highest flowing karst springs worldwide (Gary & Sharp, 2006) and the presence of A. mexicanus cavefish populations have been reported in at least 30 caves in the region (Elliott, 2016; Espinasa et al., 2018). Additionally, we included two surface populations located at the vicinity of El Abra, between the Sierras Las Crucitas and Sierra de la Colmena (San Luis Potosí, Fig. 1). All locations form part of the Pánuco river basin, which drains to the Gulf of Mexico (Hudson, 2003).

Sample collection and DNA extraction

Sampling was made from August 5th to 8th, 2016, during the wet season. A total of eight samples were collected in cave and surface habitats (Fig. 1): Pachón cave (n = 3, two females and one male), Rascón river (n = 2, both males), and Micos river (n = 3, two females and one male). All the samples were collected under the auspices of the permit SGPA/DGVS/02438/16, delivered by SEMARNAT. Fish specimens were captured using hand nets, kept alive and rapidly transported in their environment water to a place with semi-sterile conditions to avoid environmental and cross-contamination. Fish were euthanized following conditions stated in the collecting permit SGPA/DGVS/02438/16. The complete digestive tract was dissected and rinsed with sterile saline solution to remove any adherent material and preserved in 2 ml cryovials with Invitrogen™ RNAlater™ Stabilization Solution. The samples were kept in liquid nitrogen during sampling (5 days), and once in the laboratory were stored at −72 °C. All dissected specimens were preserved as vouchers in 95% ethanol and deposited in the Colección Nacional de Peces, IBUNAM (Mexico). Under sterile conditions and in a laminar flow hood, the stomach sack was excised from the rest of the intestine (approximately 1 cm²) and immediately transferred into the lysing matrix provided by the FastDNA Isolation Kit for Soil (MP Bio, Santa Ana, CA, USA). The genomic DNA was extracted and purified following manufacturer’s instructions. DNA quality and quantity were assessed by both, agarose gel electrophoresis and spectrophotometrically with a NanoDrop 2000 (Thermo Fisher Scientific, Waltham, MA, USA).

At each sampling point, seven physicochemical parameters were measured: pH, water temperature, dissolved oxygen (DO), electrical conductivity, total dissolved solids, salinity, and oxido-reduction potential (ORP), using a multiparameter Hanna HI9828 (HANNA Instruments).

16S rRNA library preparation and sequencing

The stomach microbiome for the eight A. mexicanus samples were assayed using the Ion 16S Metagenomic kit (LifeTech A2616), which uses seven hypervariable regions: V2, V3, V4, V6-7, V8 and V9. Amplification occurred in two multiplex pools (one for the V2, V4, and V8, and one for the V3, V6-7 and V9), with 25 PCR cycles, using the Ion Xpress Barcoded adapters. Amplicon sizes of the 6 fragments averaged 254 bp, within the range of 215–295 bp in E. coli (2014; Ion Toffent, Gilford, NH, USA). Libraries were quantified using Real-Time PCR, normalized and pooled. Emulsion PCR used Ion PGM Template OT2 400 kit. Sequencing was performed on an Ion Torrent PGM, using Ion 316 chip kit v2 BC with 4 and 8 chips.

Sequences analysis

Bacterial taxonomic assignments were performed using Ion Reporter Metagenomics 16S software v5.2 (Life Technologies, Carlsbad, CA, USA) with the following analysis pipeline: (a) a minimum Phred score of Q20; (b) only reads with no mismatches in the barcode, up to three mismatches in the primer sequence, and a minimum length of 150 bp were included, and (c) reads were trimmed by primers at both ends. After primer trimming and length checking, reads were assigned to a hash table to get unique sequences and their abundance. Taxonomic assignment for unique sequences with a minimum of 10 copies was performed using MicroSEQ ID 16S Reference Library v2013.1 (Thermo Fisher Science) and Greengenes v13.5 (McDonald et al., 2012) databases. Reads were grouped into two bins: 90–97% match (approximating the family level) and >97–99% identity (genus level). This tool also provided information on relative taxonomic abundance by consensus and by primer within the dataset.

Rarefaction curves and alpha diversity indices were calculated using the alpha diversity script in the Quantitative Insights into Microbial Ecology software (QIIME v. 1.9.0; Caporaso et al., 2010). For the diversity estimates, sequences were rarefied at 88,076 reads, which was the smallest read count across all fish samples (sample Micos female 1). The species richness was estimated with the Chao index, and diversity was determined by Shannon index (that quantifies the uncertainty in predicting the species identity in a random sampling) and Simpson index (that quantifies evenness and dominance of the species in the sample) (Hill et al., 2003).

The number of unique and shared OTUs (Operational Taxonomic Units) at genus level between sexes and populations was visualized with Venn diagrams using the Euleer package in R software (R Core Team, 2015). We evaluated the capability of each sequenced region to describe the stomach microbiota composition through the survey of the relative proportion of bacterial taxonomic groups retrieved by each hypervariable region of the 16S rRNA gene. Posteriorly, a non-metric multidimensional scaling (NMDS) ordination plot using Bray–Curtis dissimilarity matrix was constructed with the vegan package in R to depict community structure patterns in two dimensions. Permutation tests (n = 1,000) were used to fit significant water physico-chemical properties (P < 0.1) onto the ordination plot as vectors. Analysis of similarity (ANOSIM) was used to test for significant differences in microbiome structure between fish samples (999 permutations, P < 0.001). To analyze the distribution of OTUs at family level with relative abundances >1% of the total reads across all the samples, a heatmap was constructed using the Manhattan dissimilarity matrix and Ward’s hierarchical clustering algorithms with the plots package in R. The 16S rRNA sequences were deposited in the NCBI Sequence Read Archive (SRA) under the Bioproject PRJNA487659 (SAMN09907737 –SAMN09907742).

Diversity and composition of bacterial community in wild populations obtained by 16S rRNA variable regions

A total of 1,969,556 reads were recovered, with an average length of ∼200 bp. We observed that all rarefaction curves from fish samples attain an asymptote (Fig. S1), indicating that the sequencing effort was sufficient to cover the bacterial community. Alpha diversity indices varied among the fish samples, but didn’t show a clear pattern among populations (Table 1). In terms of species richness, samples from Pachón cave showed the highest Chao index values, followed by Micos and finally Rascón. Diversity indices across samples were similar, with a range for the Shannon diversity index between 2.26 and 3.65, while for the Simpson diversity index the samples varied between 0.614 and 0.884, with the lowest values in both indices for the Rascón male 1. Finally, sex did not seem to impact biodiversity indexes (Fig. S2).

In total, we recovered 16 phyla, 120 families, and 178 genera, where the four following lineages accounted for 93.8% of total sequences across all the samples: Gammaproteobacteria (31.7%), Firmicutes (29.9%), Bacteroidetes (16.2%) and Betaproteobacteria (15%). Differences in composition and distribution of phylogenetic groups across populations were recovered (Fig. 2). The stomach microbiome of surface Rascón samples was dominated by Bacteroidetes (61.4% male 1 and 29.5% male 2) and Firmicutes (29% male 1 and 63.5% male 2), with a lesser proportion of Actinobacteria, Verrucomicrobia and Fusofacteria. In contrast, a different microbiota was recovered from Pachón cave samples, which were dominated by Gammaproteobacteria (37.6% male 1 and 47.9% female 1), Betaproteobacteria (17% male 1, and 34.5% female 1) and Firmicutes (26.3% male 1 and 11.9% female 2). The microbiome from Micos surface population was more similar to the Pachón cave population than to the Rascón surface population, being Gammaproteobacteria (68.5% female 1 and 32.2% female 2), Betaproteobacteria (19.4% female 1 and 13.6% female 2) and Firmicutes (7.8% female 1 and 40.5% female 2) the most abundant groups, with a lesser proportion of other groups such as Cyanobacteria and Fusobacteria.

Comparison of 16S rRNA hypervariable regions for assignment of microbial diversity

The percentage of mapped reads was not homogenous among the hypervariable regions (Fig. S3). The region with the highest percentage of mapped reads was V3 (from 31.8% in female 1 to 44.9% in female 2, both from Micos river), followed by V8 (from 9.5% in Micos female 2 to 29.6% in Micos female 1), V4 (from 11.7% in Rascón male 2 to 18% in Rascón male 1 and Pachón female 1), and V6–V7 (from 6.6% in Pachón female 1 to 20.9% in Rascón male 1). The regions with the lowest percentage of mapped reads were V2 (from 1.4% in Rascón male 2 to 10.9% in Pachón female 1) and V9 (from 0.8% in Rascón male 1 to 12% in Micos female 1). In addition, not all the regions could distinguish the different taxonomic groups in the same way. For instance, V3 was the region with the greatest coverage of taxonomic groups (from six to nine groups), allowing the detection of the four most abundant groups (Gammaproteobacteria, Firmicutes, Betaproteobacteria, and Bacteroidetes) across all the fish samples. On the other hand, V9 covered the least number of taxonomic groups (from 1 to 2 groups) and could only detect Gammaproteobacteria and Betaproteobacteria groups.

Shared and dominant stomach microbiome between A. mexicanus populations

The core stomach microbiome shared by the three A. mexicanus populations was formed by eight genus OTUs (Fig. 3) belonging to Acinetobacter, Aeromonas, Citrobacter and Klebsiella (class Gammaproteobacteria), Prevotella (Bacteroidetes), Vogesella (class Betaproteobacteria), Clostridium (Firmicutes), and Cetobacterium (Fusobacteria). The abundance distribution of the core microbiome differs among the three populations: Rascón river had the highest abundance of Prevotella (89.9%), while Micos river had more Acinetobacter (60.6%) and Vogesella (11.8%), and Pachón cave had more Vogesella (40.4%) and Acinetobacter (34.3%). Although they come from contrasting environments, Micos river and Pachón cave specimens shared more OTUs between them, than with samples from Rascón river. In addition, the most abundant bacterial groups per locality showed a clear difference across populations (Table 2): the dominant genus OTUs in Pachón cave samples belonged to Proteobacteria, Firmicutes and Actinobacteria, in Micos river samples belonged to Proteobacteria and Firmicutes, and in Rascón river samples belonged to Proteobacteria, Firmicutes and Bacteroidetes.

Environmental determinants for the microbiome composition across populations

As the contrast between environments (cave vs. surface) did not seem to account for the observed differences in stomach microbiota, we assessed the association of water physicochemical variables with respect to β-diversity by means of NMDS analysis (Fig. 4), which grouped the samples in two clusters: bacterial community from the Rascón river samples formed a unique cluster, and those from Micos river and Pachón cave formed another cluster. This strong clustering trend was also confirmed by the ANOSIM (p < 0.001) analysis. Fitting the physicochemical variables onto the ordination plot showed that dissolved oxygen (DO) was the only significant explanatory variable for the bacterial community structure. The DO values were considerably lower in both Pachón cave (2.97 mg/L) and Micos river (4.43 mg/L) relative to that from Rascón river (8.2 mg/L). Interestingly, observed differences in other variables such as conductivity, for which contrasting values were obtained (e.g., 226 µS/cm and 580 µS/cm for Rascón and Micos rivers, respectively), did not emerge as important factor as DO did for stomach microbiota community distribution (Table 3).

The heatmap of family OTUs with relative abundance >1% of the total reads showed the same pattern of grouping as the NMDS analysis (Fig. 5), separating the stomach microbiome in two groups: the first one included the samples from Rascón river, and the second one included a mixture of both surface and cave environments (Micos river and Pachón cave). Moraxellaceae (Gammaproteobacteria) was the most abundant family with 14.5% of the total reads, followed by Enterobacteriaceae (Gammaproteobacteria) with 11.5% of the total reads and Prevotellaceae (Bacteroidetes) with 11.1% of the total reads. In general, fish samples from Micos river and Pachón cave showed greater abundance of bacteria belonging to Moraxellaceae and Enterobacteriaceae, while fish samples from Rascón showed greater abundance of bacteria belonging to Prevotellaceae.

Different hypervariable regions of the 16S rRNA show a differential bacterial identification of the stomach microbiome of A. mexicanus

Next generation sequencing techniques have revolutionized the way of studying microbial diversity. However, inherent biases arise due to insufficient coverage of primers (Hong et al., 2009; Ghyselinck et al., 2013), primer-template mismatches (Sipos et al., 2007), and unequal amplification (Polz & Cavanaugh, 1998), distorting the comprehensive view of complex bacterial communities. The choice of the primer sets and covered regions can deeply affect the profile of microbiome composition (Kumar et al., 2011; Soergel et al., 2012). Therefore, the use of multiple primers targeted to amplify different hypervariable regions of the 16S rRNA gene is an innovative strategy in metagenomic protocols that allows detecting a larger array of bacterial groups and enabling a more accurate microbiome description (Milani et al., 2013; Barb et al., 2016).

Previous studies in fish microbiota only targeted one or two hypervariable regions (typically V1–V3, V4, V6 or V8) (Rawls et al., 2006; Kim, Brunt & Austin, 2007; Sun et al., 2009; Sullam et al., 2012; Franchini et al., 2014; Liu et al., 2016). Here, we reported the first fish microbiota characterization incorporating the sequencing of seven hypervariable regions of the 16S rRNA gene (V2, V3, V4, V6-7, V8, and V9) using the Ion 16S ™ Metagenomics Kit. This methodology has been recently reported for other models (Barb et al., 2016; Cornejo-Granados et al., 2017; Sperling et al., 2017), suggesting that more than one hypervariable region is needed to provide reliable diversity inferences (Barb et al., 2016).

In this study, the V3, V4 and V8 regions were the most accurate both for number of assigned reads per sample and for the OTUs distribution, while the other regions seem to underestimate the bacterial diversity. Among the hypervariable regions of the 16S rRNA gene, V3 and V4 are the less biased and the most used for microbiome studies (Cai et al., 2013; Milani et al., 2013). Additionally, studies using the Ion 16S™ Metagenomics Kit showed differences in the microbiome diversity depending on the region analyzed, recovering similarly to our results that the regions V3, V4 and V8 were the most informative-less biased regions (Cornejo-Granados et al., 2017; Sperling et al., 2017). Our results in the stomach microbiome characterization of A. mexicanus, as previous studies, showed the benefits of using different non-overlapping regions, and the potential biases of single marker identification in the microbiome characterization, which could affect not only the richness but also diversity in microbiome studies.