Skip to main content

Community recovery dynamics in yellow perch microbiome after gradual and constant metallic perturbations



The eco-evolutionary processes ruling post-disturbance microbial assembly remain poorly studied, particularly in host-microbiome systems. The community recovery depends not only on the type, duration, intensity, and gradient of disturbance, but also on the initial community structure, phylogenetic composition, legacy, and habitat (soil, water, host). In this study, yellow perch (Perca flavescens) juveniles were exposed over 90 days to constant and gradual sublethal doses of cadmium chloride. Afterward, the exposure of aquaria tank system to cadmium was ceased for 60 days. The skin, gut and water tank microbiomes in control and treatment groups, were characterized before, during and after the cadmium exposure using 16s rDNA libraries and high throughput sequencing technology (Illumina, Miseq).


Our data exhibited long-term bioaccumulation of cadmium salts in the liver even after two months since ceasing the exposure. The gradient of cadmium disturbance had differential effects on the perch microbiota recovery, including increases in evenness, taxonomic composition shifts, as well as functional and phylogenetic divergence. The perch microbiome reached an alternative stable state in the skin and nearly complete recovery trajectories in the gut communities. The recovery of skin communities showed a significant proliferation of opportunistic fish pathogens (i.e., Flavobacterium). Our findings provide evidence that neutral processes were a much more significant contributor to microbial community turnover in control treatments than in those treated with cadmium, suggesting the role of selective processes in driving community recovery.


The short-term metallic disturbance of fish development has important long-term implications for host health. The recovery of microbial communities after metallic exposure depends on the magnitude of exposure (constant, gradual), and the nature of the ecological niche (water, skin, and gut). The skin and gut microbiota of fish exposed to constant concentrations of cadmium (CC) were closer to the control negative than those exposed to the gradual concentrations (CV). Overall, our results show that the microbial assembly during the community recovery were both orchestrated by neutral and deterministic processes.

Video Abtract.


Resilience refers to the capacity of a natural ecosystem to maintain a stable state after facing different exogenous disturbances, both in terms of amplitude and frequency [1]. Introduced first by Holling (1973), the concept of resilience was redefined to incorporate the idea of recovery following a temporary disruption [2, 3], not simply the ability to resist this disturbance in the first place [4]. Both ecological concepts, “resistance” and “recovery,” were simultaneously considered as measurable components that together represent resilience [4]. In other microbial studies, the term “resistance” is synonymous with resilience [5] using Holling‘s definition. Notwithstanding, “sensitivity” (inverse of resistance) is also sometimes used to represent the degree to which a community changes in response to disturbance [6]. The recovery rate, time to reach an equilibrium state, and the distance to an alternate stable state are quantitative measures that can be used to compare the resilienc e[4, 7,8,9] and improve our understanding of ecosystem recovery [6, 10]. In this study, we will employ the term “recovery” to describe the pattern of eco-evolutionary change that occurs when a community returns to an alternative stable state.

The recovery of microbial communities depends on the type, duration, intensity, and variability of a disturbance. More importantly, microbial recovery can be impacted by the initial community structure, phylogenetic composition, legacy, and the type of habitat (soil, water, host). After antibiotic treatment, the complete recovery of initial bacterial community composition is rarely achieved, as reported in various host-microbiota systems from honeybees [11] to humans [12]. The incomplete recovery of gut microbiota ecosystems after antibiotic administration results in a shift of the microbial composition to an alternative equilibrium called an “alternate stable state” [6, 13, 14]. This compositional shift occurs when resistance or recovery is weak and/or when the intensity of disturbance is high. Although the understanding of factors that drive such regime shifts to an alternative equilibrium in microbial ecosystems will have tremendous impacts in various fields of application (e.g., personalized medicine, agriculture, bioremediation), this phenomenon is still poorly studied.

The relative roles of ecological and evolutionary processes in the recovery of the structure of microbial communities are still to decipher. Theoretically, the nature of these processes can be neutral (stochastic) [15, 16], or selective (deterministic) [17, 18], the latter being driven either by environmental filtering or competitive exclusion [19, 20], the former by demographic sampling effects alone. In the context of community recovery, a small number of studies revealed that deterministic processes drive bacterial succession dynamics in a soil bacterial community disrupted either by a depletion gradient of nutrients [21], a thermal shock [22], or a rainfall rehydration of dry soil [23].

In the present study, we assessed the relative contribution of neutral and deterministic processes in the recovery of the yellow perch (Perca flavescens) microbiome assembly following an experimental metallic exposure gradient. Polymetallic contamination in aquatic ecosystems mostly results from exposure to acid mine drainages (AMD) occurring around the world [24,25,26,27,28,29,30,31,32,33,34]. For instance, in the natural Canadian lakes, the cadmium (Cd) concentration reaches 9 ppb (parts per billion) in perch liver/water [35, 36], and it has a clear quantitative impact on the perch physiology, gene expression, and genotype diversity [37]. In the same polluted lake system studied by Couture et al. (2008), the microbial communities’ assembly in water has evolved under chronic exposure to a gradient of trace metals due to the AMD expelled in the environment, leaving substantial genotypic signatures of adaption in the taxonomic and functional repertories of AMD communities [34]. Given that yellow perch juveniles can tolerate sublethal doses of cadmium without encountering significant physiological damage or death [38, 39], this host-microbiota model system is well suited to study microbiota recovery following metal exposure stress. In the laboratory, yellow perch juveniles underwent exposure to sublethal doses of cadmium chloride (CdCl2), the accumulation of which was tested in the water and within perch liver. The recovery of community structure and function in water and host microbiome were then studied and compared between constant and variable regimes of metallic stress, which was defined by the levels of Cd detected in liver and water samples. To disentangle the effect of the xenobiotic from host development [40, 41] on bacterial strain recruitment ontogeny, microbiota assembly was also assessed in stable conditions as a control regime. Our expectation was that constant exposure to cadmium chloride, due to its severe implications for host and microbial community physiology, would impede community recovery most severely than in the gradual exposure experimental group.


Fish rearing

The experiment is described in Fig. 1: Schema 1. Briefly, there were two acclimation periods: one in a standard container (1500L) and the second in 24 tanks (36 L) with an independent filtering system circuit for each aquarium. The fish juveniles were reared within the same physicochemical conditions (photoperiod, pH, ammonia, nitrogen dioxide). Throughout the experimental period, to maintain viable conditions for perch in each water tank, fecal and uneaten food particles were removed daily using specific pressing tubes for each set of experimental conditions. A volume of 15 L of water was renewed two times a week for each tank (Fig. 1:Schema 1).

Fig. 1

Schematic illustration of the perch microbiome recovery experiment

Exposure regimes to cadmium

Cd-treated and control (Ctrl) tanks were designed into two cadmium chloride exposure regimes (8 tanks per regime), and one negative control regime (8 tanks) (Fig. 1: Schema 1). The yellow perch in treated tanks were exposed to cadmium chloride (CdCl2) dissolved in water. Under the regime of constant CdCl2 concentration exposure (CC), the cadmium chloride was initially added at 0.8 ppb, then increased to reach a target theoretical concentration of 9 ppb (parts per billion) by the end of the first month (T1). The CdCl2 concentration was adjusted to 9 ppb every 5 days during two additional months until the end of treatment (third month, T3) where the measured concentration reached an average of 5.8 ppb. Under the regime of variable CdCl2 concentration (CV), the CdCl2 was initially added at 0.6 ppb, then the concentration was gradually increased every 5 days to meet the target theoretical concentration of 9 ppb by the end of the third month. The measured concentration reached an average of 6.8 ppb at the end of treatment (third month, T3). The maximal CdCl2 concentration was settled at 9 μg/mL, which was within range of concentrations detected in yellow perch liver in contaminated Canadian lakes [35, 36].

Recovery after the exposure to cadmium

The cadmium administration was stopped after the third month (T3). The experiment was extended 2 months (T5) after T3 to test the recovery of microbiome assembly in water and host.

Host-microbiota and water sampling

Briefly, we selected 144 mucosa samples of skin (2 times × 3 regimes × 8 tanks × 3 replicates) and 144 gut (2 times × 3 regimes × 8 tanks × 3 replicates) samples corresponding to T0 (no cadmium) and T3 (ultima cadmium treatment). Also, 48 water samples (2 times × 3 regimes × 8 tanks × 1 technical replicate) from T0 and T3 were included. At the end of recovery time (T5), 72 skin mucosa samples (1-time × 3 regimes × 8 tanks × 3 replicates) and 72 gut samples (1-time × 3 regimes × 8 tanks × 3 replicates) were collected from the host. In addition, 240 samples (5 times × 3 regimes × 8 tanks × 2 technical replicates) of water (2 L) microbial filter (0.22 μm) were sampled between T3 and T5 at an interval time of 15 days corresponding to five recovery time points (TR1, TR2, TR3, TR4, and T5).

Metal concentration in water and fish liver

Every week until the end of the CdCl2 exposure regimes, we measured the concentration of trace metals of cadmium (Cd), zinc (Zn), and copper (Cu) in the yellow perch liver and water tanks using the ICPMS (Ionization coupled mass spectrometry) technology available at INRS (Institut National de la Recherche Scientifique). For further details of the measurement of Cd in liver preceded by acid digestion and lyophilization see our under review study Cheaib et al. (2019). Two-way analysis of variance (ANOVA), Tukey’s test, and Wilcoxon rank test were applied to test the significance of cadmium accumulation in liver and water over time, and between treatment groups.

DNA extraction, libraries preparation, and 16S amplicons sequencing

DNA was extracted using the Qiagen DNeasy blood and tissue kit for skin mucosa, and TRIzol organic phase followed by BEB (back extraction buffer) and PCI (phenol/chloroform/isoamyl alcohol 25:24:1) solutions for all gut samples. The V3–V4 hypervariable region of the universal rDNA 16S gene (Werner et al. 2012) was amplified using universal specific primers. The libraries of amplicons were prepared using a set of 384 combinations of adaptors, processed in one sequencing run, on an Illumina Miseq sequencing machine. Reactions of PCR were verified by electrophoresis on 2% agarose gel, purified, and quantified by fluorescence for the double-strand DNA concentration using Quant-iT™ PicoGreen™ dsDNA Assay Kit (Thermo Fischer Scientific).

Bioinformatics and biostatistics analyses

Reads preprocessing and OTUs clustering

Sequence analysis was performed with our bioinformatic pipeline as described previously [42, 43]. In the first instance, we used SICKLE Version 1.2 to trim the reads (>Q30 Phred quality score) followed by utilizing PANDASEQ Version 2.11 [44] assembler for merging paired-end read into a single merged reads (~ 350 bp) corresponding to the amplified 16S rRNA V3–V4 hypervariable region (347 F-805 R). Based on an approach of de novo sequence clustering before the taxonomic assignment, reads were clustered into OTUs at 97% identity with USEARCH Version 9 (Edgar RC. 2010) and filtered out using UNOISE2 algorithm [45] to discard chimeric sequences, putatively produced during PCR amplification cycles using OTUs were annotated using RDP database as previously described in our pipeline [42, 43] to. Community structure and composition of metacommunities were analyzed across time and treatments by richness (OTUs count), evenness (Shannon index), and the Gunifrac phylogenetic distance [46] using the vegan [47] and Rhea [48] packages in R.

We then calculated, the alpha-diversity indices (richness and evenness) and beta-diversity (phylogenetic distance) differences between experimental groups and used rank statistics tests (Kruskal-Wallis/Wilcoxon) to assess their significance. The resulting p values for pairwise comparisons in alpha and beta-diversity were corrected for multiple testing using the Benjamini-Hochberg method (Benjamini and Hochberg, 1995). Note that the Beta-diversity was calculated using the generalized UniFrac metric [49], which considers both the dominant and the rare OTUs. The permutational multivariate analysis of variance (henceforth referred to as PERMANOVA) was applied to the Gunifrac distance matrices to explain the sources of variations including the experimental groups. To test homogeneity of variances, we performed the multivariate homogeneity test which a Multi-Response Permutation Procedure (MRRP) of within versus among group dissimilarities dispersions of Gunifrac distances. The non-metric multi-dimensional scaling (NMDS) was performed to visualize Gunifrac distances in a reduced space with k = 2 dimensions. For statistics comparison of one-dimensional statistics of multiple groups, we used the non-parametric Kruskal-Wallis rank-sum test because of the strong assumption of the normal distribution of OTUs abundance being rarely assumed.

The alpha-diversity variation across time and per treatment was predicted and plotted with linear mixed effect models using the ratio of richness/evenness as a response variable, time, and cadmium concentrations in water and liver as fixed effects, with the categorical variable tank as a random effect.

Using the lmer R package, for water, the model was used as following in R:

Model<- lmer (Richness/Shannon.effective~Time+Cd.Water+(1|Tank), data=mixdata, REML = TRUE)

whereas for each host habitat (Skin, Gut), we employed the following model:

Model_host<- lmer(Richness/Shannon.effective~Time+Cd.Liver+Cd.Water+(1|Tank), data=mixdata, REML = TRUE)

The confidence interval was then predicted using the predict interval() function in R.

Post-OTUs analysis, networks, and function prediction

Structure and diversity measures of groups (control and treatments) were compared with rank statistics tests (Kruskal-Wallis/Wilcoxon) adjusted with BH (Benjamin-Hochberg) test for multiple corrections, and the p value < 0.05 as a threshold of statistical significance. To understand the role of relative abundance of OTUs on the similarities of community structure, correlation networks of communities (samples) were constructed using the Spearman coefficient as a robust approach of correlation detection [50]. Significant positive and negative correlations were filtered and false discovery rate (FDR) was assessed with B-H test for multiple corrections. Next, network visualization and analysis were performed with Cytoscape software [51]. The network centrality was analyzed using the “Network Analyzer” plugin in Cytoscape. The betweenness centrality of a node was calculated as the total number of the shortest paths from all nodes to all other nodes that pass this node [52]. The centrality of the nodes reflected their importance in transmitting information between hubs; it does not depend on the feature of node degree, which describes the total node connectivity. The size of nodes was proportional to the number of OTUs in each sample, and the coefficient of the significant correlation between two nodes was inversely proportional to the size of the edge. Finally, functional profiles of each community type at every time point were predicted using the software TaxforFun [53].

Neutral and deterministic models to asses the recovery of community assembly

In the null hypothesis, the neutral model [16] (Sloan et al. 2006) assumes a beta distribution of OTU abundance. Using the non-linear partial least square method [41], which estimates the migration rate (m) of OTUs from their source to a destination community, the model predicts the frequencies of OTUs. The estimated migration rate (m) is the probability that a random loss (death or emigration) of an OTU in a destination community is replaced by dispersal from the source community. Comparing the predicted versus observed frequencies, we can determine which OTU fits the model in each host and water community, at every time point, across both control and treatment groups. The goodness of fit to the model was measured using the coefficient of determination R2 (R2 > 0.5) within a confidence interval of 95%, where increased strength of goodness of fit to the model suggests an essential role of stochastic processes in the microbiome assembly.


Cadmium concentration bioaccumulation in the fish liver during recovery time

Interestingly, the concentration of cadmium ions measured with ICPMS increased significantly in the fish liver even after 2 months from stopping exposure. The Cd concentration increased from 0.4 ppb to 1 ppb in the variable CdCl2 regime (CV) and from 0.5 ppb to 1.17 ppb in the constant CdCl2 regime (CC). However, in the water, as expected, the Cd concentration significantly decreased from 6.4 ppb to 1.06 ppb in CV, and from 5.8 ppb to 1.34 ppb in CC (Tables 1 and 2). Consequently, the accumulation of Cd in liver and water was always significantly higher in treatments CC and CV compared to the control group (Table 3). Similar Cd concentrations observed among treatment groups CC and CV in the water at times T3 and T5 (expected at maximum Cd concentration added in tanks) were observed in fish liver only at time T5 (Table 3).

Table 1 Statistics of Cd concentrations in water and fish liver over time and treatments
Table 2 Statistics of cadmium concentration variation over time in water tanks and fish livers
Table 3 Statistics of cadmium concentration variation among treatments in water tanks and fish livers

Genotypic signatures of community recovery

At the alpha-diversity level, to investigate how far diversity metrics could be used as indicators for metacommunity structure recovery, both richness and the evenness were calculated in water and host-microbial communities. In the host microbiome, time had a significant effect on diversity measures within all groups between times T3 and T5. The richness and evenness have significantly increased over time in skin microbiota and significantly decreased in the gut microbiota. Over the five recovery time points (TR1, TR2, TR3, TR4, T5), temporal comparisons in the water microbial communities associated with each experimental group did not show a significant change of evenness for CC and CV, but did during TR2–TR4 within the control group (Ctrl). Within these communities, a significant change in richness during the whole recovery period except TR2–TR3 was also found for CC and CV (Additional file 7: Table S1a-b).

In contrast, both richness and evenness in the control group of skin microbiota significantly fluctuated over the recovery period (T3–T5, after CdCl2 addition had stopped). At time T3, the pairwise comparison of CC and CV against the control group (CC-Ctrl and CV-Ctrl) revealed significant differences in microbial richness in the gut and evenness in the skin. At time T5, statistical tests did not detect any significant change in diversity measures between all groups for gut and skin microbiome (Additional file 7: Table S1-c); however, as at T3, the evenness of skin microbiome at T5 was significantly divergent between cadmium treatments, (p value = 0.0063) (Additional file 7: Table S1-c).

The comparative analysis of richness and evenness among water and host microbiota showed convergent patterns of diversity between the water and the skin communities before the disturbance and after the recovery (Additional file 1: Figure S1).

The predicted alpha-diversity values along with the fitted linear mixed effect model for water communities showed a significant drop in treatment values compared to the control group. On the other hand, for the host communities (skin and gut), they increased under the selection regimes and decreased during the recovery time span (Fig. 2).

Fig. 2

Predicted alpha-diversity plots by linear mixed model. Alpha-diversity in water and host-microbial communities over time and among treatments is predicted using the linear mixed model. The richness/evenness ratio were considered as response variables, the fixed effects were defined by time and cadmium concentration (in water and liver), and tanks were taken as random effects. Over time, the predicted alpha diversity in host microbial communities (skin, gut) highlights stable trends of the Control group compared to the treatments. However, all groups of the water microbial communities decrease overtime. Constant cadmium regime (CC) is in orange, variable cadmium regime (CV) is in yellow, and control (control) is in green

In summary, except for the linear mixed effect model results, the observed patterns of alpha diversity metrics changes across the experiment did not show a clear trend over the course of the experiment. Nonetheless, increasing evenness and richness was a general trend for the skin while decreasing evenness and increasing richness was representative of the gut microbiome community recovery.

Beta diversity (Gunifrac) between samples was compared using a PERMANOVA and a multivariate test for variance homoscedasticity. By T3—at peak cadmium exposure—significant differences (p < 0.05) among treatments were observed in all microbial communities of water and host (Table 4 ; Additional file 2: Figure S2); and by T5, both variable (CV) and constant (CC) cadmium exposure treatments retained differences in skin communities compared to the controls despite the recovery period. Surprisingly, given our expectation that cadmium exposure would have a major impact on community recovery, a high similarity in the community phylogenetic structure between the control and CC groups was detected among gut microbial communities at T5. Beta-diversity between treatments (CC, CV, Ctrl) was always significantly divergent at each time point in water samples except for the observed convergence between CC and CV at recovery time TR2 (Table 4; Additional file 2: Figure S2). The comparison of beta-diversity showed a community structure divergence (p value < 0.001) between water, skin, and the gut microbiota before the disturbance and after the recovery (Additional file 3: Figure S3). The results show that water microbiome at time T3 is not representative enough of fish microbiome (see the blue cluster in the phylogram of CC, page 2 of Additional file 3: Figure S3). However, at the recovery time T5, the water was not representative of the fish microbiome in gradual selection regime (see the blue cluster in the phylogram of CV, page 3 of Additional file 3: Figure S3).

Table 4 Phylogenetic divergence in host and water microbiomes

Microbial taxonomic composition change during recovery

At T5, no significant changes were observed between groups at the phyla level in the water, but Actinobacteria in the gut, and both Euryarchaeota and Tenericutes in the skin, were significantly different between control and treatments (CC and CV). Additional file 8: Table S2 details several taxa that showed significant differential abundance between treatments (Ctrl, CV, CC). Of particular importance, putative pathogenic genus Flavobacterium was significantly enriched in the skin for both groups of CdCl2 exposed fish at T5, despite the recovery period. In the gut microbiome, Syntrophococcus was the only genus to be significantly different between treatments (Fig. 3; Additional file 8: Table S2). In the water, significant differences in taxonomic abundance between CV and Ctrl were restricted to one genus (Kiloniella) at recovery time TR1 and two genera (Marinobacter and Perlucidibaca) at T5. No significant differences in taxonomic composition were detected between CC and Ctrl in the water at T5. Overall, statistical analysis of taxonomic composition dynamics over time within each treatment during the recovery period revealed several minor differences (see Additional file 9: Tables S3 for more details).

Fig. 3

Taxonomic composition dynamics of host communities. Stacked barplots show the most abundant taxa (> 0.5%) overtime in the gut, skin and water microbiomes. The genera that significantly changed among treatments and control at T5 are summarized in Additional file 8: Table S2

On the other hand, the pairwise comparison of taxonomic composition between different type of communities at each time point and for each experimental group showed significant divergence between microbial communities of gut, skin, and water. At the recovery time (T5), Tenericutes, Euryarchaeota, and Firmicutes were inherently associated (significantly upregulated) with the gut microbiome; Actinobacteria and Bacteroidetes, on the other hand, were specific to skin microbiome; with Fibrobacteres and Actinobacteria implicated with the water microbiome. Despite all this, the Proteobacteria were found to be prevalent and common in water and skin microbiome. At the selection time (T3), Fibrobacteria and Actinobacteria were scarcely abundant and were picked up as differentially abundant (Additional file 6: Figure S6). To delineate the most relevant taxa (at the genus level) significantly changing between the communities, we have performed a pairwise test on an overall comparison of skin, gut, and water at each time point for different treatment groups. The results plotted in heatmaps in the Fig. 4 clearly reveal that each community type has an inherent signature and the corresponding proportions of genera differed between control and treatments over time with high similarity between the CC and Ctrl at time T5 (Fig. 4; Additional file 10: Table S4).

Fig. 4

Heatmaps of differential abundance among host and water communities. This figure from left to right includes 9 heatmaps of the significant taxonomic fingerprints at the genus level between gut, the skin and the water at times T0 (first column), T3 (second column), and T5 (third column) in the control (first row), the CV (second row), and the CC (third row) groups. The hierarchical clustering of the relative abundance of phyla which significantly changed over time was performed using Ward’s method and Bray–Curtis dissimilarity distance. Vegan package and pheatmap () function in R were used

Correlational analysis (Additional file 4: Figure S4) revealed a positive relationship between specific genera and concentrations of cadmium in perch liver and water. In aquaria treated with CdCl2 (CV and CC), cadmium concentrations in water and liver showed strong significant positive correlations with seven genera from the gut microbiome, each of which represented a different phylum and had a strong negative correlation with the relative abundance of Mycoplasma. In the skin microbiome of both CC and CV, 15 genera (Sphingomonas; Haloarcula; Legionella; Flavobacterium; Ameyamaea; Dokdonella; Shigella; Massilia; Mycoplasma; Polaromonas; Pseudomonas; Rhodobacter; Rhodococcus; Shewanella; Syntrophococcus) showed significant profiles of positive or negative correlations with the Cd concentrations in the liver (Additional file 4: Figure S4). Divergent profiles between CC and CV were only observed for the correlations of Shewanella and Syntrophococcus with cadmium concentrations. Similar correlation profiles between these groups were observed in the water (Additional file 5: Figure S5).

Correlational networks of host and water microbiome

In the host and water communities, the network analysis of samples correlations showed a partitioned community distribution between treatment groups, at time T3, and overlapping patterns during the recovery period (Figs. 5 and 6). At time T0, the correlational networks of host microbiome showed unstructured topology with on average a fewer number of edges. An edge can represent significant low correlations (Spearman’s Rho correlation > 0.5) between samples from different groups with the node size proportional to the richness of each sample. The topological distribution of nodes in the network was further analyzed by comparing the betweenness centrality to the eigenvalue centrality (Fig. 7). The results indicate a shift in the mean of the centrality metrics between the control (which is higher) and the cadmium selection regimes. The plots of eigenvalue centrality versus betweenness centrality clearly reveal that these communities shift at times T3 for skin microbiome, as well as at time T3 for gut microbiome. The high betweenness centrality observed in control reflects the efficiency of network centrality measure to predict the effect of perturbation on the community structure during the selection phase, but not during the recovery time as centrality median shifts at T5 was not observed (Fig. 7). The same centrality analysis was obtained for water microbiome networks resulting in similar patterns at time T3 (results not shown).

Fig. 5

Recovery dynamics of the networks of host communities. The networks organization is based on nodes betweenness centrality among treatments and Control. Unstructured patterns in the networks were observed at T0. Node size represents sample richness. The strength of correlation (Spearman correlation from 0.3 to 1) between two nodes is inversely proportional to the size of the edge. This network was built using R and Cytoscape software. Constant Cadmium samples (CC) are in orange, variable cadmium samples (CV) are in yellow, and control (Control) samples are in green

Fig. 6

Recovery dynamics of the network of water communities. The networks organization at every resilience time TR1, TR2, TR3, TR4, and WT5 is based on nodes betweenness centrality among treatments and control. The network modules easliy distinguishable between groups since T0. The size of nodes (sample richness) at the beginning of TR1 and at the end of the time TR4 showed shifts in the community richness. The strength of correlation (Corr. Spearman from 0.5 to 1) between two nodes is inversely proportional to the size of the edge. This network was built using R and Cytoscape software. constant cadmium samples (CC) are in orange, variable cadmium samples (CV) are in yellow, and control (control) samples are in green

Fig. 7

Centrality plots of host microbiome networks. This figure summarizes relationships of betweenness centrality versus eigenvalue centrality of host microbiome networks among treatments and at each time point. The results show evidence of shift in centrality medians between the control regime (which is higher) and the cadmium selection regimes. The plots of eigenvalue centrality versus betweenness centrality clearly reveal that centrality shift at time T3 for skin and gut microbiome

Recovery of microbial functional diversity a time T5

At time T0, an ANOVA of functional richness within the metacommunity showed a significantly higher average of functional diversity in gut and skin microbiomes compared to water microbial communities. Surprisingly stable in water communities, functional diversity did not show any significant divergence among the treatment and control group at T3 regardless of the community type (skin, water, gut). The lack of treatment effect observed may well have been masked by the strong influence of time over microbial diversity (Cheaib et al. 2019 submitted ISMEJ). However, at T5, the functional diversity of skin microbiome was significantly higher in the control group than in treatment groups according to the ANOVA (CC-CV(p value) = 0.04; CV-Ctrl(p value) = 0.0055; CC-Ctrl(p value) = 0.45) (Fig. 8). In the gut microbiota, no significant changes in functional diversity were detected between treatments (CC-CV(p value) = 0.3; CV-Ctrl (p value) = 0.54; CC-Ctrl (p value) = 0.58).

Fig. 8

Function diversity dynamics in host and water microbiome. Boxplots of functions profiles were predicted from the matrices of taxa count using the software Tax4Fun. The statistical significance (p value < 0.05) found using ANOVA followed by FDR (false discovery rate) test are represented with asterisks points (0.001: “***,” 0.01: “**,” 0.05: “*”)

The role of neutral and deterministic processes in the recovery of host microbiota

The goodness of fit of host and water microbial communities to the non-linear partial least square model (NLS) was high (R2 > 0.5), supporting the theory of predominant neutrality (Additional file 11: Table S5). To disentangle gut and skin microbiota ontogeny from the cadmium effect, the NLS model was deployed using the control as a reference. A comparison of observed versus predicted OTU frequencies revealed that the percentage of neutral OTUs in skin and gut microbiota (Fig. 9) at the recovery time T5 is higher in the control group compared to those in treatments at T3 and T5. The same analysis was undertaken in the water and the percentage of neutral versus non-neutral OTUs showed the same trends across Ctrl, CC, and CV at T5. Overall, we noted a preponderance of OTUs that fitted the neutral model in all comparisons. The majority of the OTUs that did not fit the neutral model was assigned to Mycoplasma species (indeed no mycoplasma sp. OTUs fitted the neutral model), which can be seen in Fig. 10 as well as Additional file 12: Table S6 and Additional file 13: Table S7. The neutral process was much more prevalent in the control group at time T3 and T5 as compared to the treatment groups.

Fig. 9

Percentage of neutral OTUs over time and treatment. Using the non-linear least squares model (NLS), the percentage of OTUs that fit the neutral model within a confidence interval of 95% showed variable trends between communities across time and treatments. A goodness of fit R2> 0.5 was considered as the significant threshold of neutrality fit. The cadmium treatment invoked stochasticity in the water communities, while in gut and skin communities, the percentage of neutral OTUs remained higher in the control compared to treatments.

Fig. 10

Demographic variation of metacommunity neutrality across water and host microbiome. This figure summarizes the scatterplots of neutral model fitting the whole metacommunity ( gut skin and water) at times T0 (first column), T3 (second column) and T5 (third column) in the control (first row), the CV (second row), and the CC (third row) groups. Neutral OTUs are shown in black, non-neutral are depicted in grey, while the red is Mycoplasma sp. OTUs. We see no Mycoplasma sp. OTUs that fit the neutral model in the whole metacommunity


Our data clearly show that long-term bioaccumulation of cadmium occurs in the Perca flavescens liver on exposure to aqueous cadmium salts. Our data also showed that the cadmium persists at high concentrations even once the treatment has been stopped for 2 months. We have already shown that cadmium treatment clearly impacts both the skin and gut microbial communities, as compared to controls (Cheaib et al. 2019, under review). Recovery was the focus of the current study, and microbial communities post-exposure showed different routes to (and extents of) recovery in those associated with the skin and gut once cadmium treatment was ceased. In the skin, evenness—the extent to which different microbes in a community share similar abundances—and richness increased during the recovery phase in cadmium-treated fish. Beta-diversity comparisons, meanwhile, revealed significant differences between all experimental cohorts (Ctrl, CC, CV) in water and skin niches. Among gut microbial communities, decreasing richness and increasing evenness were observed over the recovery period. Beta-diversity metrics indicated few significant differences between cadmium and control treatments. Crucially, by the end of the recovery period in the gut, functional richness was comparable between tests and control, a potential signal of full community recovery. We used models to assess the relative roles of microbial assembly in the different groups. We found evidence that neutral processes were a more prevalent contributor to microbial community turnover in control treatments than in those treated with CdCl2—likely indicating the role of selective processes in driving community recovery. Overall our data do not strongly support our prediction that the most extreme cadmium exposure (CC) would lead to the least successful recovery. Instead, CC and CV treatments, especially in the gut, demonstrated a good degree of recovery, both in terms of both alpha and functional diversity.

At the end of the third month of exposure (T3), the cadmium concentration in the liver was significantly higher in CC and CV than in the control group. These concentration differences were still observed two months (T5) after the gradual clearing of Cd started. The liver plays a major role in the accumulation, excretion, and biotransformation of contaminants like metalloids [54, 55], and bioaccumulated metals remain at high concentrations in the liver due to its depuration function of other organs (such as gills and muscles) [56]. Long-term bioaccumulation of cadmium has been documented in perch and other biological systems [39, 57,58,59,60], as has its effect on ecosystem services in soil and water [34, 61] and metazoan gut ecosystems [62,63,64,65]. This study not only confirms the chronic bioaccumulative effect of Cd but also suggests that the sequestered Cd in perch liver presumably cannot predict the regime of exposure (CC, CV), as the concentration did not significantly vary in livers between both regimes, CC and CV, at T5. Over the recovery period, the concentration of Cd in water significantly decreased, but Cd was not completely removed from the tank system since it has a strong affinity to the tank silicone gaskets, and has high competitiveness with Zn for the debris of organic compounds always available in the water aquarium ecosystem [66].

The water microbial communities showed few differential abundances of taxa differences during the recovery period (Fig. 3). Furthermore, the microbial functional diversity in water remained stable throughout the experiment, and no significant differences between treatments were found during the exposure or the recovery periods. However, the community beta-diversity at the phylogenetical level between treatments (CC, CV, Ctrl) showed significant difference at each time point, suggesting a pattern of taxon-function decoupling as an adaptive strategy reported previously in lacustrine water contaminated with cadmium [34].

To assess the yellow perch microbiome recovery, we examined alpha-diversity (richness and evenness), beta-diversity (phylogenetic distance), taxonomic composition, and functional diversity (metabolic functions). Most of these measurements are commonly used as community-wide metrics to assess the recovery of microbial communities, for example, in humans [12, 67], soil [22, 68], and wastewater [69] (Vrieze et al. 2017).

In the skin microbiome, the disturbance intensity (cadmium gradient) had a differential impact on the community recovery trajectories, resulting in a significant difference of evenness (Additional file 7: Table S1c) and functional diversity (Fig. 8) between CC and CV at time T5 (Table 4). During the gradual exposure regime (CV), the cadmium may provoke an endurance effect on the skin microbiota which was progressively adapted to the cadmium accumulating in the tank system, while within the constant exposure regime (CC), abrupt diversity and taxonomic changes might have been triggered. Gradual changes are evident under stress gradients, for example, within bioreactors, the anaerobic microbiome has been shown to gradually adapt following ammonium disturbanc e[70]. Consequently, the significant divergence in the functional diversity between CV-Ctrl and CV-CC, not between CC-Ctrl, perhaps indicates a unique adaptive evolution signature of skin microbiome under CV regime. Therefore, the skin communities from CV and CC may have followed a different recovery trajectory after adaptation. Strikingly, the recovery of skin microbiota of the most extreme exposure (CC) appeared to be the most successful, when considering the convergence of richness, evenness, and functional diversity between CC and Ctrl. However, significant differences among CC, CV, and Ctrl in terms of phylogenetic divergence (Table 4) and taxonomic composition shifts (Fig. 3, Additional file 8: Table S2; Additional file 9: Table S3; Additional file 10: Table S4) suggest this recovery was incomplete. For instance, a significant increase in fish pathogens like Flavobacterium, Legionella and opportunists like Mycoplasma was detected in both cadmium groups (CC and CV) compared to control. The relative abundance of Flavobacterium was significantly lower in the control group with a low percentage (< 0.5 %). Perturbation with cadmium can facilitate the proliferation of opportunistic pathogens, this concern has been found in other studies of fish microbiota recovery after exposure to antibiotic [71] and triclosan biocide [72]. Similar taxonomic changes in both exposure regimes (CC and CV) were expected [73]. Overall, the cadmium disturbance may cause a shift to an alternative stable state, demonstrating differential and incomplete recovery of the skin microbiota in CC and CV.

In the gut microbiome, the recovery routes were different; at time T5, there was only a significant evenness convergence between CC and CV. Overall, the few significant differences in taxonomy, as well as the phylogenetic divergence (Additional file 8: Table S2) between CC-CV and Ctrl-CV, but not between CC-Ctrl, suggests a full recovery of the gut microbiota in CC and gradual recovery in CV. At the level of taxonomic composition, overall the dominance of opportunists Tenericutes was also a feature of farmed Eurasian perch (Perca fluviatilis) gut microbiota studied in a context of stress predation [74], although they were not found in the wild Eurasian perch [75].

In the skin and gut microbiota, the significant increase in diversity (evenness and richness) over the recovery period (T3–T5) was consistent with the diversity increase in other host-associated studies such as the recovery of the fathead minnow gut microbiome from a low-level triclosan exposure [72], the human intestinal microbiota post-infection [67], the murine gut microbiome exposed to antibiotics in early life [76], and the molasses wastewater [69]. Further, the functional redundancy observed in all water and gut microbial communities is a major adaptive strategy behind resistance and recovery [34, 69]. Lastly, the significant divergence of skin and gut microbiota diversity over time within the control group suggests persistent divergence from the initial community structure due to microbiota ontogeny through the developmental stage of fish juveniles [41].

Our findings demonstrate a relative role of neutral processes shaping the bacterial communities’ recovery following exposition to metallic stressors. According to the neutral model fit, the percentage of neutral OTUs in skin and gut microbiota was significantly higher in the control group compared to CdCl2 treated groups, which provides evidence that neutral processes are the major contributor in the microbiota assembly in non-stressed yellow perch, therefore suggesting that selective processes are at play in driving the community recovery in stress-exposed groups. Furthermore, Mycoplasma sp. are a dominant species in perch microbiome, implicated in literature for other fish species [43, 77]. The inability of neutral models to explain the abundance of any OTUs for Mycoplasma sp. in the current study suggests that these bacterial strains can quickly adapt to the host environment. Our study is the first to investigate the relative importance of neutrality and determinism in driving post-disturbance assembly of the host-associated microbiome.


This study not only elucidates the long-term bioaccumulation effect of toxic metals on biological systems but also suggests that the sequestered cadmium in the fish liver will not likely predict the magnitude of exposure regime (constant or variable). The effect of cadmium exposure on microbial communities is also varying and dependent on the nature of the host it is originating from. Surprisingly, after recovery, skin and gut microbiota of fish exposed to constant concentrations of cadmium (CC) were closer to the control group than those exposed to the gradual concentrations (CV). In the skin, the metallic perturbation caused a shift to an alternative stable state, leading to an incomplete recovery and therefore, facilitating the proliferation of opportunistic pathogens (like Flavobacterium). In the gut, the functional and phylogenetic diversity measurements suggest a complete community recovery in the CC group and gradual recovery in the CV group. The selective pressure exerted by cadmium on host and water microbiota may have left adaptive evolution patterns conserving functional diversity at the expense of taxonomic diversity. In both skin and gut microbiota, the recovery was associated with a significant increase of evenness and richness in the skin and vice versa in the gut. In the control group, as expected, the significant divergence from the initial community structure confirms the dynamic of bacterial strains through the developmental stage of fish juveniles. Consequently, community recovery was affected by both cadmium pressure and host development. In addition, our results have shown that the microbial assembly rules during the community recovery were both orchestrated by neutral and deterministic processes. In the water, community recovery was driven by a substantial role of phylogenetic structuring resulting from a combined pattern of stochasticity and cadmium-induced selective pressure, in which the causality remains unknown. Further studies are needed to quantify the interactions of neutrality and determinism in driving post-disturbance assembly of the host-associated microbiome during recovery.

Availability of data and materials

Sequencing data are available in the Sequence Read Archive (SRA) database at NCBI under the BioProject ID PRJNA556617


16 s rDNA:

16S ribosomal DNA


Analysis of variance


Back extraction buffer


Benjamini-Hochberg correction test


Cadmium constant concentration


Cadmium salts


Comités de Protection des Animaux de l’Université Laval


Regime of negative control


Cadmium variable concentration


False discovery rate


Multiple response permutation procedure


Non-linear least squares model


Non-metric multi-dimensional scaling


Operational taxonomic units


Polymerase chain reaction


Permutational analysis of variance


Parts per billion


Ribosomal Database Project


  1. 1.

    Holling CS. Resilience and stability of ecological systems. Annu Rev Ecol Syst. 1973;4:1–23.

    Article  Google Scholar 

  2. 2.

    Pimm SL. The complexity and stability of ecosystems. Nature. 1984;307:321.

    Article  Google Scholar 

  3. 3.

    Grimm V, Wissel C. Babel, or the ecological stability discussions: an inventory and analysis of terminology and a guide for avoiding confusion. Oecologia. 1997;109:323–34.

    CAS  Article  PubMed  Google Scholar 

  4. 4.

    Hodgson D, McDonald JL, Hosken DJ. What do you mean, ‘resilient’? Trends Ecol Evol. 2015;30:503–6.

    Article  PubMed  Google Scholar 

  5. 5.

    Ziegler M, Seneca FO, Yum LK, et al. Bacterial community dynamics are linked to patterns of coral heat tolerance. Nat Commun. 2017;8:14213.

  6. 6.

    Shade A, Peter H, Allison SD, et al. Fundamentals of microbial community resistance and resilience. Front Microbiol. 2012;3.

  7. 7.

    Ingrisch J, Bahn M. Towards a comparable quantification of resilience. Trends Ecol Evol. 2018;33:251–9.

    Article  PubMed  Google Scholar 

  8. 8.

    Cabrol L, Poly F, Malhautier L, et al. Management of microbial communities through transient disturbances enhances the functional resilience of nitrifying gas-biofilters to future disturbances. Environ Sci Technol. 2016;50:338–48.

    CAS  Article  PubMed  Google Scholar 

  9. 9.

    Scheffer M, Carpenter SR, Dakos V, van Nes EH. Generic indicators of ecological resilience: inferring the chance of a critical transition. Annu Rev Ecol Evol Syst. 2015;46:145–67.

    Article  Google Scholar 

  10. 10.

    Allison SD, Martiny JBH. Resistance, resilience, and redundancy in microbial communities. Proc Natl Acad Sci. 2008;105:11512–9.

    Article  PubMed  Google Scholar 

  11. 11.

    Raymann K, Shaffer Z, Moran NA. Antibiotic exposure perturbs the gut microbiota and elevates mortality in honeybees. PLOS Biol. 2017;15:e2001861.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  12. 12.

    Dethlefsen L, Relman DA. Incomplete recovery and individualized responses of the human distal gut microbiota to repeated antibiotic perturbation. Proc Natl Acad Sci. 2011;108:4554–61.

    Article  PubMed  Google Scholar 

  13. 13.

    Lozupone CA, Stombaugh JI, Gordon JI, et al. Diversity, stability and resilience of the human gut microbiota. Nature. 2012;489:220–30.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  14. 14.

    Sommer F, Anderson JM, Bharti R, et al. The resilience of the intestinal microbiota influences health and disease. Nat Rev Microbiol. 2017;15:630–8.

    CAS  Article  PubMed  Google Scholar 

  15. 15.

    Hubbell SP. Neutral theory in community ecology and the hypothesis of functional equivalence. Funct Ecol. 2005;19:166–72.

    Article  Google Scholar 

  16. 16.

    Sloan WT, Lunn M, Woodcock S, et al. Quantifying the roles of immigration and chance in shaping prokaryote community structure. Environ Microbiol. 2006;8:732–40.

    Article  PubMed  PubMed Central  Google Scholar 

  17. 17.

    Chase JM. Community assembly: when should history matter? Oecologia. 2003;136:489–98.

    Article  PubMed  Google Scholar 

  18. 18.

    Webb CO, Ackerly DD, McPeek MA, Donoghue MJ. Phylogenies and community ecology. Annu Rev Ecol Syst. 2002;33:475–505.

    Article  Google Scholar 

  19. 19.

    Cadotte MW, Davies TJ, Regetz J, et al. Phylogenetic diversity metrics for ecological communities: integrating species richness, abundance and evolutionary history. Ecol Lett. 2010;13:96–105.

    Article  PubMed  Google Scholar 

  20. 20.

    Stegen JC, Lin X, Fredrickson JK, et al. Quantifying community assembly processes and identifying features that impose them. ISME J. 2013;7:2069–79.

    Article  PubMed  PubMed Central  Google Scholar 

  21. 21.

    Song H-S, Renslow RS, Fredrickson JK, Lindemann SR. Integrating ecological and engineering concepts of resilience in microbial communities. Front Microbiol. 2015;6.

  22. 22.

    Jurburg SD, Nunes I, Stegen JC, et al. Autogenic succession and deterministic recovery following disturbance in soil bacterial communities. Sci Rep. 2017;7.

  23. 23.

    Placella SA, Brodie EL, Firestone MK. Rainfall-induced carbon dioxide pulses result from sequential resuscitation of phylogenetically clustered microbial groups. Proc Natl Acad Sci. 2012;109:10931–6.

    Article  PubMed  Google Scholar 

  24. 24.

    Wendt-Potthoff K, Koschorreck M. Functional groups and activities of bacteria in a highly acidic volcanic mountain stream and lake in Patagonia, Argentina. Microb Ecol. 2002;43:92–106.

    CAS  Article  PubMed  Google Scholar 

  25. 25.

    Wang F, Liu C, Liang X, Wei Z. Remobilization of trace metals induced by microbiological activities near sediment-water interface, Aha Lake, Guiyang. Chin Sci Bull. 2003;48:2352–6.

    CAS  Article  Google Scholar 

  26. 26.

    Gough HL, Dahl AL, Nolan MA, et al. Metal impacts on microbial biomass in the anoxic sediments of a contaminated lake. J Geophys Res Biogeosciences. 2008;113:G02017.

    CAS  Article  Google Scholar 

  27. 27.

    Hudson-Edwards KA, Jamieson HE, Lottermoser BG. Mine wastes: past, present, future. Elements. 2011;7:375–80.

    Article  Google Scholar 

  28. 28.

    Moser M, Weisse T. The most acidified Austrian lake in comparison to a neutralized mining lake. Limnologica. 2011;41:303–15.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  29. 29.

    Urbieta MS, Toril EG, Aguilera A, et al. First prokaryotic biodiversity assessment using molecular techniques of an acidic river in Neuquén, Argentina. Microb Ecol. 2012;64:91–104.

    Article  PubMed  Google Scholar 

  30. 30.

    Douglas TA, Loseto LL, Macdonald RW, et al. The fate of mercury in Arctic terrestrial and aquatic ecosystems, a review. Environ Chem. 2012;9:321.

    CAS  Article  Google Scholar 

  31. 31.

    Masmoudi S, Nguyen-Deroche N, Caruso A, et al. Cadmium, Copper, Sodium and Zinc Effects on Diatoms: from Heaven to Hell — a Review. Cryptogam Algol. 2013;34:185–225.

    Article  Google Scholar 

  32. 32.

    Stankovic S, Moric I, Pavic A, et al. Investigation of the microbial diversity of an extremely acidic, metal-rich water body (Lake Robule, Bor, Serbia). J Serbian Chem Soc. 2014;79:729–41.

    CAS  Article  Google Scholar 

  33. 33.

    Valente T, Rivera MJ, Almeida SFP, et al. Characterization of water reservoirs affected by acid mine drainage: geochemical, mineralogical, and biological (diatoms) properties of the water. Environ Sci Pollut Res. 2015;23:6002–11.

    CAS  Article  Google Scholar 

  34. 34.

    Cheaib B, Le Boulch M, Mercier P-L, Derome N. Taxon-function decoupling as an adaptive signature of lake microbial metacommunities under a chronic polymetallic pollution gradient. Front Microbiol. 2018;9.

  35. 35.

    Pyle GG, Rajotte JW, Couture P. Effects of industrial metals on wild fish populations along a metal contamination gradient. Ecotoxicol Environ Saf. 2005;61:287–312.

    CAS  Article  PubMed  Google Scholar 

  36. 36.

    Couture P, Rajotte JW, Pyle G. Seasonal and regional variations in metal contamination and condition indicators in yellow perch (Perca flavescens) along two polymetallic gradients. III. Energetic and physiological indicators; 2008.

    Google Scholar 

  37. 37.

    Bougas B, Normandeau E, Pierron F, et al. How does exposure to nickel and cadmium affect the transcriptome of yellow perch (Perca flavescens) – Results from a 1000 candidate-gene microarray. Aquat Toxicol. 2013;142–143:355–64.

    CAS  Article  PubMed  Google Scholar 

  38. 38.

    Giguère A, Campbell PG, Hare L, et al. Influence of lake chemistry and fish age on cadmium, copper, and zinc concentrations in various organs of indigenous yellow perch (Perca flavescens). Can J Fish Aquat Sci. 2004;61:1702–16.

    Article  Google Scholar 

  39. 39.

    Campbell PGC, Giguère A, Bonneris E, Hare L. Cadmium-handling strategies in two chronically exposed indigenous freshwater organisms—the yellow perch (Perca flavescens) and the floater mollusc (Pyganodon grandis). Aquat Toxicol. 2005;72:83–97.

    CAS  Article  PubMed  Google Scholar 

  40. 40.

    Sylvain F-É, Derome N. Vertically and horizontally transmitted microbial symbionts shape the gut microbiota ontogenesis of a skin-mucus feeding discus fish progeny. Sci Rep. 2017;7:5263.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  41. 41.

    Burns AR, Stephens WZ, Stagaman K, et al. Contribution of neutral processes to the assembly of gut microbial communities in the zebrafish over host development. ISME J. 2016;10:655–64.

    CAS  Article  PubMed  Google Scholar 

  42. 42.

    Sylvain F-É, Cheaib B, Llewellyn M, et al. pH drop impacts differentially skin and gut microbiota of the Amazonian fish tambaqui (Colossoma macropomum). Sci Rep. 2016;6:32032.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  43. 43.

    Llewellyn MS, McGinnity P, Dionne M, et al. The biogeography of the atlantic salmon (Salmo salar) gut microbiome. ISME J. 2016;10:1280–4.

    Article  PubMed  Google Scholar 

  44. 44.

    Masella AP, Bartram AK, Truszkowski JM, et al. PANDAseq: paired-end assembler for illumina sequences. BMC Bioinformatics. 2012;13:31.

    CAS  Article  Google Scholar 

  45. 45.

    Edgar RC. UNOISE2: improved error-correction for Illumina 16S and ITS amplicon sequencing. bioRxiv. 2016:081257.

  46. 46.

    Chen J, Bittinger K, Charlson ES, et al. Associating microbiome composition with environmental covariates using generalized UniFrac distances. Bioinforma Oxf Engl. 2012;28:2106–13.

    CAS  Article  Google Scholar 

  47. 47.

    Oksanen J, Blanchet FG, Friendly M, et al. Community ecology package. R Packag. 2016:2.3–3.

  48. 48.

    Lagkouvardos I, Fischer S, Kumar N, Clavel T. Rhea: a transparent and modular R pipeline for microbial profiling based on 16S rRNA gene amplicons. PeerJ. 2017;5:e2836.

    Article  PubMed  PubMed Central  Google Scholar 

  49. 49.

    Chen W, Zhang CK, Cheng Y, et al. A comparison of methods for clustering 16S rRNA sequences into OTUs. PLoS One. 2013;8:e70837

    CAS  Article  Google Scholar 

  50. 50.

    Weiss S, Van Treuren W, Lozupone C, et al. Correlation detection strategies in microbial data sets vary widely in sensitivity and precision. ISME J. 2016;10:1669–81.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  51. 51.

    Shannon P, Markiel A, Ozier O, et al. Cytoscape: A Software Environment for Integrated Models of Biomolecular Interaction Networks. Genome Res. 2003;13:2498–504.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  52. 52.

    Röttjers L, Faust K. From hairballs to hypotheses-biological insights from microbial networks. FEMS Microbiol Rev. 2018;42:761–80.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  53. 53.

    Aßhauer KP, Wemheuer B, Daniel R, Meinicke P. Tax4Fun: predicting functional profiles from metagenomic 16S rRNA data. Bioinformatics. 2015;31:2882–4.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  54. 54.

    Silver S, Phung LT. BACTERIAL HEAVY METAL RESISTANCE: New Surprises. Annu Rev Microbiol. 1996;50:753–89.

    CAS  Article  PubMed  Google Scholar 

  55. 55.

    Defo MA, Spear PA, Couture P. Consequences of metal exposure on retinoid metabolism in vertebrates: A review. Toxicol Lett. 2014;225:1–11.

    CAS  Article  PubMed  Google Scholar 

  56. 56.

    Jezierska B, Witeska M. The metal uptake and accumulation in fish living in polluted waters. In: Twardowska I, Allen HE, Häggblom MM, Stefaniak S, editors. Soil and Water Pollution Monitoring, Protection and Remediation. Springer Netherlands; 2006. p. 107–14.

  57. 57.

    Giguère A, Campbell PGC, Hare L, Couture P. Sub-cellular partitioning of cadmium, copper, nickel and zinc in indigenous yellow perch (Perca flavescens) sampled along a polymetallic gradient. Aquat Toxicol. 2006;77:178–89.

    CAS  Article  PubMed  Google Scholar 

  58. 58.

    Klinck JS, Green WW, Mirza RS, et al. Branchial cadmium and copper binding and intestinal cadmium uptake in wild yellow perch (Perca flavescens) from clean and metal-contaminated lakes. Aquat Toxicol. 2007;84:198–207.

    CAS  Article  PubMed  Google Scholar 

  59. 59.

    Xie L, Lambert D, Martin C, et al. Cadmium biodynamics in the oligochaete Lumbriculus variegatus and its implications for trophic transfer. Aquat Toxicol. 2008;86:265–71.

    CAS  Article  PubMed  Google Scholar 

  60. 60.

    Nirola R, Megharaj M, Saint C, et al. Metal bioavailability to Eisenia fetida through copper mine dwelling animal and plant litter, a new challenge on contaminated environment remediation. Int Biodeterior Biodegrad. 2016;113:208–16.

    CAS  Article  Google Scholar 

  61. 61.

    Chen YP, Liu Q, Liu YJ, et al. Responses of soil microbial activity to cadmium pollution and elevated CO2. Sci Rep. 2014;4:4287.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  62. 62.

    Liu Y, Li Y, Liu K, Shen J. Exposing to cadmium stress cause profound toxic effect on microbiota of the mice intestinal tract. PLOS ONE. 2014;9:e85323.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  63. 63.

    Zhang S, Jin Y, Zeng Z, et al. Subchronic exposure of mice to cadmium perturbs their hepatic energy metabolism and gut microbiome. Chem Res Toxicol. 2015;28:2000–9.

    CAS  Article  PubMed  Google Scholar 

  64. 64.

    Qian B, Mian L, Peizhan C, et al. Sex-dependent effects of cadmium exposure in early life on gut microbiota and fat accumulation in mice. Environ Health Perspect. 2017;125:437–46.

    Article  Google Scholar 

  65. 65.

    Šrut M, Menke S, Höckner M, Sommer S. Earthworms and cadmium – heavy metal resistant gut bacteria as indicators for heavy metal pollution in soils? Ecotoxicol Environ Saf. 2019;171:843–53.

    CAS  Article  PubMed  Google Scholar 

  66. 66.

    Pinter TBJ, Stillman MJ. Kinetics of zinc and cadmium exchanges between metallothionein and carbonic anhydrase. Biochemistry. 2015;54:6284–93.

    CAS  Article  PubMed  Google Scholar 

  67. 67.

    Singh P, Teal TK, Marsh TL, et al. Intestinal microbial communities associated with acute enteric infections and disease recovery. Microbiome. 2015:3.

  68. 68.

    Griffiths BS, Philippot L. Insights into the resistance and resilience of the soil microbial community. FEMS Microbiol Rev. 2013;37:112–29.

    CAS  Article  PubMed  Google Scholar 

  69. 69.

    De Vrieze J, Christiaens MER, Walraedt D, et al. Microbial community redundancy in anaerobic digestion drives process recovery after salinity exposure. Water Res. 2017;111:109–17.

    CAS  Article  PubMed  Google Scholar 

  70. 70.

    Regueiro L, Carballa M, Lema JM. Microbiome response to controlled shifts in ammonium and LCFA levels in co-digestion systems. J Biotechnol. 2016;220:35–44.

    CAS  Article  PubMed  Google Scholar 

  71. 71.

    Navarrete P, Mardones P, Opazo R, et al. Oxytetracycline treatment reduces bacterial diversity of intestinal microbiota of atlantic salmon. J Aquat Anim Health. 2008;20:177–83.

    Article  PubMed  Google Scholar 

  72. 72.

    Narrowe AB, Albuthi-Lantz M, Smith EP, et al. Perturbation and restoration of the fathead minnow gut microbiome after low-level triclosan exposure. Microbiome. 2015;3:6.

    Article  PubMed  PubMed Central  Google Scholar 

  73. 73.

    Regueiro L, Carballa M, Lema JM. Outlining microbial community dynamics during temperature drop and subsequent recovery period in anaerobic co-digestion systems. J Biotechnol. 2014;192:179–86.

    CAS  Article  PubMed  Google Scholar 

  74. 74.

    Zha Y, Eiler A, Johansson F, Svanbäck R. Effects of predation stress and food ration on perch gut microbiota. Microbiome. 2018;6:28.

    Article  PubMed  PubMed Central  Google Scholar 

  75. 75.

    Bolnick DI, Snowberg LK, Hirsch PE, et al. Individual diet has sex-dependent effects on vertebrate gut microbiota. Nat Commun. 2014;5:4500.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  76. 76.

    Cho I, Yamanishi S, Cox L, et al. Antibiotics in early life alter the murine colonic microbiome and adiposity. Nature. 2012;488:621–6.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  77. 77.

    Holben WE, Williams P, Saarinen M, et al. Phylogenetic analysis of intestinal microflora indicates a novel mycoplasma phylotype in farmed and wild salmon. Microb Ecol. 2002;44:175–85.

    CAS  Article  PubMed  Google Scholar 

Download references


We thank Martin Llewellyn and Eleanor Lindsay for comments on the earlier version of this manuscript. Thanks to the Derome laboratory team for their help in sampling. Thanks to Brian Boyle for his precious advice regarding library preparation, and to all his team members running the sequencing platform at IBIS. We also thank Michel Lavoie for his advice on Cadmium manipulation. Umer Zeeshan Ijaz is supported by NERC Independent Research Fellowship NERC NE/L011956/1 as well as Lord Kelvin Adam Smith Leadership Fellowship (Glasgow).


Funding for this project was provided by the National Sciences and Engineering Council of Canada (NSERC), a Canadian Research Chair in genomics and conservation of aquatic resources to ND.

Author information




ND and BC conceived the experiment. BC and HS reared fish tanks and performed sampling. BC performed the DNA extraction. BC performed molecular biology experiments including libraries preparation and quantification. BC produced and analyzed the results. BC and UZI performed functional diversity analysis. BC and ND wrote the manuscript. All authors reviewed the manuscript. All authors edited the manuscript and approved the final draft.

Corresponding author

Correspondence to Bachar Cheaib.

Ethics declarations

Ethics approval

All animals used in this study were treated in accordance with the approved University of Laval’s CPAUL (Comités de Protection des Animaux de l’Université Laval).

Consent for publication

Not applicable

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information

Additional file 1: Figure S1. Dynamic of alpha-diversity divergence between host and water communities. The significant ANOVA results of alpha diversity between water (W), Skin(S) and Gut (Gut) communities in Control, CV and CC groups before and during disturbance, and after recovery period are represented with asterisks on the boxplots (0.001 : “***”, 0.01 : “**”, 0.05 : “*”).

Additional file 2: Figure S2. Beta-diversity divergence at the treatment level. This file combines all the NMDS (non-metric Multi-Dimensional Scaling) plots showing first two dimensions in the ordination of when using generalized Unifrac distance measure of water and host-microbial communities. The NMDS plots and PERMANOVA revealed a significant separation between different treatments and control (for the pairwise, see Table 4 for adjusted p-values after Benjamini-Hochberg correction in PERMANOVA and MRPP tests) at T0, T3, and T5 for skin and gut microbiota, and at T0, T3, TR1-TR4, and T5 for water microbial communities.

Additional file 3: Figure S3. Beta-diversity divergence at the community level. This file combines all the NMDS (non-metric Multi-Dimensional Scaling) plots and phylograms based on generalized Unifrac distances between water and host-microbial communities. The NMDS plots and PERMANOVA revealed a significant separation of among all type of communities per time (T0, T3, and T5) and treatment (Control, CC, CV).

Additional file 4: Figure S4. Heatmaps of cadmium with taxa diversity and composition in host and water communities. The correlations indicate a gradient from positive (blue) to negative (red) along a colour gradient, with rows representing diversity measures (richness, evenness) as well as cadmium concentrations, and columns indicating taxonomic levels. The gut microbiome in the constant CdCl2 (CC) and variable CdCl2 (CV) regimes showed a negative correlation between Mycoplasma and diversity indices. A strong positive correlation between Actinomycetales is noticeable in the CV. For the skin microbiome, not only Actinomycetales, but also Burkholderiales, and Chromatiales showed strong positive correlations with CV. This figure was produced using the Rhea package.

Additional file 5: Figure S5. Heatmaps of cadmium with taxa diversity and composition in water during the recovery period. The correlations indicate a gradient from positive (blue) to negative (red) along a colour gradient, with rows representing diversity measures (richness, evenness) as well as cadmium concentrations, and columns indicating taxonomic levels. In the constant CdCl2 (CC) and variable CdCl2 (CV) regimes, correlations of cadmium with taxa abundance showed variable profiles over time. This figure was produced using the Rhea package.

Additional file 6: Figure S6. Heat trees and stacked bar plots of water and host microbiome structure. This figure summarizes pairwise comparison of the community composition of water and each of the host communities for different treatments (Ctrl, CC and CV). Additionally, stacked bar plots of relative abundance at phylum level are provided for each community (water, skin, gut). The non-grey coloring (which category the branches are upregulated in) indicates significant differences in terms of log median ratios for samples from different habitats (Gut, Skin and Water) as determined by a Wilcox rank-sum test followed by a Benjamini-Hochberg (FDR) correction for multiple testing. The heat trees were built using metacoder and stacked barplots were produced using the Rhea package.

Additional file 7: Table S1. Alpha-diversity dynamics over time and treatments. This statistical summary reveals richness or evenness changes over time (Tables 1 and 2) and between control and treatments (Table 3) in water and host communities. The significant changes of alpha-diversity indices between treatments and control were statistically tested using Kruskal-Wallis and Wilcoxon tests by applying Benjamini-Hochberg correction. The same statistics were used to compare alpha-diversity over time.

Additional file 8: Table S2. Statistical summary of taxa divergence between treatments and control after recovery. This table summarises significant taxonomic changes in water and host-microbial communities between control and treatments using the Fisher test, and by applying Benjamini-Hochberg correction.

Additional file 9: Tables S3. Statistical summary of taxa divergence over time in host and water microbial communities after recovery. These tables summarise significant taxonomic changes over time in water and host-microbial communities using Kruskal-Wallis and Wilcoxon tests, and by applying Benjamini-Hochberg correction.

Additional file 10: Tables S4. Statistical summary of differential abundance between water and host microbial communities at all taxonomic levels. These tables summarise significant taxonomic changes between water and each of the host communities (skin and gut) at times T0 (sheet1), T3 (sheet2), and T5 (Sheet3) using Kruskal-Wallis and Wilcoxon tests, and by applying Benjamini-Hochberg correction.

Additional file 11: Table S5. Statistics of the neutral model in host and water microbiomes. This table summarises the neutral model fit based on the following parameters; the migration rate ( within 95% of the confidence interval, the goodness of fit(R2), number of samples, richness, abundance cutoff, percentage of % neutral OTUs and non-neutral OTUs.

Additional file 12: Table S6. List of OTUs that accounted for those that did not fit the neutral model.

Additional file 13: Table S7. List of OTUs that accounted for those that did not fit the neutral model and were assigned to Mycoplasma species.

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Cheaib, B., Seghouani, H., Ijaz, U.Z. et al. Community recovery dynamics in yellow perch microbiome after gradual and constant metallic perturbations. Microbiome 8, 14 (2020).

Download citation


  • Fish microbiome
  • Disturbance
  • Recovery
  • Stress gradient
  • Neutrality
  • Evolutionary forces
  • Community assembly
  • Pathogens
  • Metagenomics