Loss of gut microbial diversity in the cultured, agastric fish, Mexican pike silverside (Chirostoma estor: Atherinopsidae)

Sample collection

Adult Chirostoma estor were collected from Lake Patzcuaro (LP), Mexico, in the spring of 2018 (March) (19°36′14″N, 101°37′56″W, 22 °C, altitude 2,040 masl, N = 13). In May of the same year, cultured adult fish of similar size to those from the lake were also harvested from fiberglass indoor tanks (tank culture = C) and outdoor earth ponds (extensive culture = E) (N = 20 and N = 19, respectively) located in the aquaculture biotechnology laboratory at Universidad Michoacana de San Nicolás de Hidalgo in Morelia, Michoacan, Mexico (19°41′22″N, 101°14′56″W, 23 °C, altitude 1,896 masl). Populations of intensively cultured fish were established twenty years ago using stock from Lake Patzcuaro. Fish from the laboratory were released into the outdoor earth ponds 8 months before sampling, and are designated as E.

At harvest, individual fish from each environment were measured and weighed (Table S1). Immediately after collection, the abdominal cavities of fish were dissected under aseptic conditions, the intestines (guts) extracted and their contents (digesta = D) were gently squeezed out and placed into sterile tubes with 96% ethanol. The guts of C. estor consist of a short intestine with a single loop in the midgut (Fig. S1), and were divided into two sections (anterior = A, and posterior = P gut) and placed in sterile tubes containing 96% ethanol and stored at 4 °C until later DNA extraction. All specimens were treated following EU Directives 2010/63/EU for animal experimentation and 2007/526/EC accommodation and care of animals used for experimental and other scientific purposes.

DNA extraction and sequencing

Metagenomic DNA of fish intestines (D, A, and P) from the 3 different populations was extracted with a modified CTAB extraction protocol (Doyle & Doyle, 1987) by adding lysozyme (100 mg/mL), proteinase k (20 mg/mL), lithium chloride (5M), and sodium acetate (3 M, pH 5.2) to the extraction buffer. After DNA extraction, the V3 variable region of the 16S rRNA gene was PCR-amplified with the primer pair V3-338f and V3-533r (Huse et al., 2008). PCR products were paired-sequenced (300 cycles, 2X150) on an Illumina Miniseq, following Illumina standard protocol at Centro de Investigación en Alimentación y Desarrollo, A.C. (CIAD), Mazatlan, Sinaloa, Mexico.

Bioinformatics

All amplicon sequence data derived from this work was submitted to NCBI Sequence Read Archive (SRA) (SRX11596261 to SRX11596334) under BioProject accession ID: PRJNA750495. After demultiplexing, CIAD sent back data as FastQ files which were preprocessed with PrinSeq (Schmieder & Edwards, 2011). Reads were then assembled and quality filtered with Flash v.1.2.7 software (Magoc & Salzberg, 2011), and VSEARCH (Rognes et al., 2016) was used for further processing, obtaining an abundance matrix of bacterial OTUs (operational taxonomic units clusterized at 97% identity) that was normalized using the metagenomeSeq method (McMurdie & Holmes, 2014). Taxonomic annotation was performed using VSEARCH (Rognes et al., 2016) against the SILVA138 database (Quast et al., 2013). The abundance matrix at the genus taxonomic level was used to calculate Good’s coverage and alpha diversity indexes (i.e., Chao 1 and Shannon indexes), which were estimated using the R Phyloseq library (McMurdie & Holmes, 2013). The genus-level abundance matrix was normalized using the metagenomeSeq method (McMurdie & Holmes, 2014), and the beta diversity distance matrix was calculated using Bray-Curtis dissimilarity. Beta diversity distance matrices were visualized with Non-metric Multidimensional Scaling plots (NMDS) of bacterial data grouped by environments and intestinal components, while ordination was based on between-sample dissimilarities calculated by Bray-Curtis distance and 999 permutations using library vegan 2.5-6 (Oksanen et al., 2019) in R software (R Core Team, 2013).

Although many groups define a core microbiome as those microbes with 100% occupancy across all samples (Turnbaugh et al., 2007), we considered (like other researchers) that a more relaxed threshold of greater than 80% occupancy would allow us to include more physiologically important bacteria in our analysis (Baldo et al., 2015; do Vale-Pereira et al., 2017; Sweet & Bulling, 2017; Rimoldi et al., 2019). Core microbiota was identified at the level of the genus when OTUs were present in at least 80% of samples per defined group (LP, C, E, D, A, and P) and at the level of OTU when present in at least one sample of each group (DLP, ALP, PLP, DC, AC, PC, DE, AE, PE). To compare microbiota between environments or digestive components, data was visualized with Venn diagrams (Chen & Boutros, 2011), in UpSet plots made with the package UpSetR (Conway, Lex & Gehlenborg, 2017).

Statistical analysis

All statistical analyses were performed using R software (R Core Team, 2013) and R library vegan 2.5-6 (Oksanen et al., 2019). Alpha diversity and relative abundance data were tested for normality and homoscedasticity by Shapiro–Wilk’s and Bartlett’s tests, respectively. Data sets grouped by intestinal component (Digesta, Anterior, and Posterior gut) or environment (Lake Patzcuaro, Intensive Culture, and Extensive Culture) were found to not be normally distributed. To compare beta diversity, environmental and intestinal component groups were statistically analyzed using non-metric dimensional scaling (NMDS) based on Bray-Curtis distances (Bray & Curtis, 1957), Analysis of Similarities (ANOSIM), as well as Adonis tests (Permutational Multivariate Analysis of Variance, PERMANOVA) using the Bray-Curtis index at 999 permutations. The number of reads across samples was normalized by sample size, and each taxon’s relative abundance (%) was calculated. Although results were generated for each taxonomic level, only Phylum and Genus levels are shown. Only taxa with the highest relative abundance were considered for statistical analysis, totaling 10 for phylum and 20 for genus. Phylum and genus abundance data were not normally distributed. Non-normal data were analyzed using nonparametric Kruskal–Wallis followed by Mann-Witney-Wilcoxon’s post hoc tests with statistical significance set at p < 0.05. Additional analysis to detect differential abundance was performed using the Linear discriminant analysis Effect Size (LEfSe) method (Segata et al., 2011) integrated within the Galaxy framework (https://huttenhower.sph.harvard.edu/galaxy/). In particular, the non-parametric Kruskal-Wallis sum-rank test was used to detect differentially abundant taxa, while Linear Discriminant Analysis (LDA) was used to estimate the effect size.

Predictive functional analysis by PICRUSt2

The functional metabolic profiles of 16S rDNA data were predicted using the Phylogenetic Investigation of Communities by Reconstruction of Unobserved States 2 (PICRUSt2) v2.3.0 beta software that uses HMMER (http://www.hmmer.org) to place OTUs into a reference phylogeny, followed by the castor R package to predict gene family abundances using hidden-state prediction. Pathway abundances are inferred based on the predicted sample functional profiles linked to pathway reactions using a modified version of MinPath (Douglas et al., 2019). Pathway abundance was then inferred using KEGG (Kyoto Encyclopedia of Genes and Genomes) (pathway_pipeline.py).

The output file (path_abun_unstrat_descrip.tsv) of predicted pathway abundance was loaded to R software (R Core Team, 2013). The 60 most abundant functional pathways (taking into account the total average of the samples) were visualized using a heatmap created by the ComplexHeatmap package in R (Gu, Eils & Schlesner, 2016). To determine specific differences between each of the predicted paths, the output file of PICRUSt2 (path_abun_unstrat_descrip.tsv) was loaded into the MEGAN5 package (Huson et al., 2016) which allowed us to define the metabolic categories predicted in these bacteria. These metabolic predictions were then analyzed using STAMP (statistical analysis of taxonomic and functional profiles) (Parks et al., 2014) software to further interrogate all predicted functional datasets and produce graphical depictions of key functional pathways. To contrast group metabolic predictions by environment (LP-C-E) and intestinal component (D-A-P) samples were compared using an ANOVA and Tukey-Kramer test. Pathways with a p-value < 0.05 were considered to be statistically significant. For paired analysis of all nine metabolic predictions (DLP, ALP, PLP, DC, AC, PC, DE, AE, PE) simultaneously, environmental and intestinal components were compared using a two-sided Welch’s t-test. Pathways with a p-value < 0.05 were considered to be statistically significant.

Bacterial diversity in the gut of C. estor

The bacterial community profile found in C. estor gut is highly variable from individual to individual, especially at the level of genus (Fig. 1A). At the level of phylum (Fig. 1B; Table S2) most samples contained Firmicutes (56.27% ± 4.28) and Proteobacteria (28.68% ± 3.47), of which the two most abundant classes were Gammaproteobacteria (24.35% ± 3.16) and Alphaproteobacteria (4.33% ± 1.25). Other common phyla were Cyanobacteria (4.69% ± 1.57), Actinobacteriota (3.51% ± 1.35), Bacteroidota (1.78% ± 0.37), Desulfobacterota (1.43% ± 0.44), Fusobacteriota (1.39% ± 1.05), Planctomycetota (0.97% ± 0.40), Verrucomicrobia (0.38% ± 0.10) and Acidobacteriota (0.29% ± 0.16).

Of the 1,207 annotated genera (Fig. 1A; Table 1), the most abundant bacterial genera observed within individual fish included Mycoplasma (19.95% ± 4.11), Staphylococcus (10.90% ± 3.30), Spiroplasma (10.23% ± 2.80), Aeromonas (7.51% ± 1.87), Pseudomonas (6.11% ± 1.86), Planomicrobium (4.75% ± 1.90), Clostridium (4.67% ± 1.32), Foliisarcina (2.37% ± 1.36) and Cutibacterium (2.16% ± 1.25). At the genus level, all digesta and Lake Patzcuaro samples contained significantly less Mycoplasma spp. than other intestinal components or environments (Table 1). Other differences were observed in the abundance of Clostridium, Romboutsia, Prochlorococcus, Acinetobacter, Cetobacterium, Anaerobacter, Flavobacterium, Methylobacterium-Methylorubrum, and Dongia within lake samples, while Staphylococcus was most abundant in the intensive culture and Planomicrobium, Erwinia, and Dechloromonas were most abundant in the extensive culture (Table 1).

Alpha-diversity varied considerably but not significantly by both environment and intestinal compartment, except for the Shannon index of fish digesta from Lake Patzcuaro (LP) which was significantly higher than the Shannon index of digesta from cultivated fish (p = 0.021) (Fig. 1C; Table S3). Bacterial diversity and Shannon index values from Lake Patzcuaro samples were significantly higher (p = 0.00477) than those of fish from tank cultures or outdoor ponds (Table S3).

Environmental influence on microbial communities

NMDS of the gut microbiota in fish from different environments (Fig. 2A) revealed a nearly separate cluster of the LP samples concerning the other two environments (C and E), the latter which largely overlapped instead. These observed differences were statistically supported by ANOSIM: R = 0.3119, p = 0.001 (Fig. 2A; Table S4), and PERMANOVA (Table S4). Repeating this analysis with data from the intestinal components showed complete overlapping of the three different groups (Fig. 2B), which surprisingly appeared to be significantly different from each other judging by ANOSIM (R-value of 0.06763, p =0.004 (Fig. 2B; Table S5)). PERMANOVA supported this result, showing a significant difference between digesta and anterior gut with an F value of 2.2341, p =0.001 (Table S5). ANOSIM and PERMANOVA analyses were stratified by intestinal component or environment, showing differences in digesta taken from fish harvested from the three environments (DLP, DC, and DE), differences in anterior guts from fish harvested in ALP vs AC and AE, while in the posterior gut, there was a significant difference only between PLP vs PE (Table S5).

Core microbiota analysis

Analysis to determine the core microbiota suggested that the genera Mycoplasma, Spiroplasma, Aeromonas, Pseudomonas, and Acinetobacter were present in all intestinal components of fish harvested from in all three environments (Fig. 3A). In addition to the aforementioned core genera, digesta samples from all three environments also contained Vibrio and Bacillus (Fig. 3B); all anterior intestine samples contained Mycoplasma, Staphylococcus, Spiroplasma, Aeromonas, Pseudomonas, Acinetobacter, Vibrio, and Weissella genera (Fig. 3C), while posterior intestine samples contained Mycoplasma, Spiroplasma, Aeromonas, Pseudomonas and Acinetobacter genera (Fig. 3D).

An UpSet plot was made to show OTUs (summed across replicates in 9 different groups) shared between the different environments and intestinal components (Fig. 3E). There was a total of 133 OTUs (Fig. 3E), of which only 42 were annotated to the level of genus using the SILVA database (Table S6). The two other UpSet plots show the core microbiota consisting of 4 different genera (Acinetobacter, Aeromonas, Pseudomona s, and Spiroplasma), as well as one OTU which was present in more than 80% of all the samples (Fig. S2).

Differential abundance of taxonomic groups

A total of 11 phyla and 48 genera were present in different samples, with significant taxonomic variation when comparing LP vs C, LP vs E, and C vs E (Fig. 4).

LEfSe analysis showed that fish present in the LP environment had a greater diversity of phyla than fish grown in C (7 vs 2) or E (4 vs 1) and that in general, fish grown in the E environment had a greater diversity of phyla than fish grown in the C environment (2 vs 1) (Figs. 4A, 4B and 4C). Zooming into the level of genus, there were 17 and 15 genera that were more highly abundant in LP than in C (4) or E (7) groups, respectively (Figs. 4D and 4E). The comparison between C or E fish microbiota suggests that only 4 or 5 genera were abundant in each (Fig. 4F).

Functional metabolic prediction using PICRUSt2

A total of 440 functional pathways were predicted for the 74 different samples (from all components and environments) that were analyzed. Of these, the 60 most abundant and with statistical differences between samples were re-analyzed, with metabolic predictions shown in Fig. S3. It was predicted that certain metabolic pathways would be over-represented in digesta samples relative to others (DLP, D, and DE). The majority of the functional KO pathways predicted to occur in all 9 groups belonged to four main categories: (i) Metabolism, (ii) Genetic Information Processing, (iii) Environmental Information Processing, and (iv) Cellular Processes (Fig. 5).

Within these predicted metabolic pathways, genes associated with the metabolism of energy, cofactors, vitamins, nucleotides, carbohydrates, and amino acids were overrepresented in cultured anterior intestines and digesta samples harvested from all environments (DLP, DC, AC, and DE) (Fig. 5). Predicted differences in metabolic pathways and their regulation are shown in Figs. S3 and S4.