Skip to main content
Download PDF

The yellow perch (Perca flavescens) microbiome revealed resistance to colonisation mostly associated with neutralism driven by rare taxa under cadmium disturbance

Animal Microbiome volume 3, Article number: 3 (2021) Cite this article

Abstract

Background

Disentangling the dynamics of microbial interactions within communities improves our comprehension of metacommunity assembly of microbiota during host development and under perturbations. To assess the impact of stochastic variation of neutral processes on microbiota structure and composition under disturbance, two types of microbial habitats, free-living (water), and host-associated (skin and gut) were experimentally exposed to either a constant or gradual selection regime exerted by two sublethal cadmium chloride dosages (CdCl2). Yellow Perch (Perca flavescens) was used as a piscivorous ecotoxicological model. Using 16S rDNA gene based metataxonomics, quantitative diversity metrics of water, skin and gut microbial communities were characterized along with development and across experimental conditions.

Results

After 30 days, constant and gradual selection regimes drove a significant alpha diversity increase for both skin and gut microbiota. In the skin, pervasive negative correlations between taxa in both selection regimes in addition to the taxonomic convergence with the environmental bacterial community, suggest a loss of colonisation resistance resulting in the dysbiosis of yellow perch microbiota. Furthermore, the network connectivity in gut microbiome was exclusively maintained by rare (low abundance) OTUs, while most abundant OTUs were mainly composed of opportunistic invaders such as Mycoplasma and other genera related to fish pathogens such as Flavobacterium. Finally, the mathematical modelling of community assembly using both non-linear least squares models (NLS) based estimates of migration rates and normalized stochasticity ratios (NST) based beta-diversity distances suggested neutral processes drove by taxonomic drift in host and water communities for almost all treatments. The NLS models predicted higher demographic stochasticity in the cadmium-free host and water microbiomes, however, NST models suggested higher ecological stochasticity under perturbations.

Conclusions

Neutral models agree that water and host-microbiota assembly promoted by rare taxa have evolved predominantly under neutral processes with potential involvement of deterministic forces sourced from host filtering and cadmium selection. The early signals of perturbations in the skin microbiome revealed antagonistic interactions by a preponderance of negative correlations in the co-abundance networks. Our findings enhance our understanding of community assembly host-associated and free-living under anthropogenic selective pressure.

Background

Microorganisms drive the biogeochemical cycles of the earth and contribute towards homeostasis, immunity, physiology, behaviour [1, 2] and development [2, 3] across a wide range of metazoan hosts [4]. Host-microbiota symbioses involve complex and dynamic associations between obligate and facultative symbionts [5]. Disentangling the dynamics of microbial interactions within communities improves our comprehension of metacommunity assembly [6]. Ecological processes (i.e., dispersal, selection and ecological drift) shape these interactions and govern the assembly rules of the ecological communities [7, 8]. The impact of ecological processes on community assembly is a long-term debate in macroecology. Stochastic processes related to the neutral theory of biodiversity suggests no impact of ecological interactions on species distribution and abundances. Neutrality reflects that diversity units “Species” are ecologically equivalent, while stochasticity implies random variation in mean demographic rates [9]. In such cases, local communities are randomly connected to a single metacommunity through differing rates of migration, death and birth [10,11,12]. In contrast, deterministic related niche theory considers that environmental conditions and interspecific interactions, including competitive exclusion, determine distribution and abundance of species [13]. The niche processes imply differentiation in mean demographic rates, while ‘determinism’ reflects an absence of random variation in species’ demographic rates. In microbial ecology, the advent of culture-independent approaches such as high-throughput 16S rDNA metabarcoding paved the way for the conceptual framework of the Operational Taxonomic Unit (OTU) or Amplicons Sequences Variants (ASVs) used as units of microbial diversity. Such advancements have opened new mathematical [14,15,16,17,18,19] and network-based [6, 20, 21] models for predicting ecological interactions between microbial communities. These models helped in constructing hypotheses on types of processes driving microbiomes assemblies over evolutionary time.

Models for quantifying the neutral [14, 16, 22, 23] and deterministic [17, 24,25,26] processes in different types of microbial ecosystems continue to provide new comprehensive insights regarding the forces governing microbiome assembly. Both neutral and non-neutral processes have been evidenced as drivers of the microbial metacommunity assembly in many vertebrate microbiomes [27, 28], as well as in the natural environment [29]. For instance, neutral processes were identified as playing a major role during the development of host-associated microbial communities in different domesticated vertebrate and plant models [28, 30, 31]. When focusing on case-control surveys, the influence of selective pressure on microbiota assembly is more salient. In contrasting contexts such as human oral and gut microbiome under antibiotic therapy [32] and euryhaline fish microbiome during salinity acclimation [33], the community assembly was potentially driven by deterministic processes, with little evidence for stochastic colonisation. Nonetheless, there is much work to do to fairly understand the microbial assembly under disruptive conditions across a diverse range of host species. In the present study, we first predict that a selection gradient induced by a concentration gradient of a toxic metal will not only disrupt the host physiology, [34] but will also overwhelm the assembly of symbiont consortia. Second, we predict that the host will lose its filtering capacity to recruit the appropriate symbionts, which in turn, will translate into increased colonisation of opportunistic strains, as predicted by the colonisation resistance model [32]. Third, from the interaction networks of OTUs, we expected that the low abundant taxa might play an important role in the metacommunity assembly [35,36,37,38,39]. To this end, we measured the effect of directional selection along the developmental stages of the host organism -Yellow Perch juveniles (Perca flavescens) were exposed to two selection regimes: a constant (9 ppb) and a gradual (0.8 to 9 ppb) exposure to non-lethal doses of cadmium chloride (CdCl2), over 90 days. The taxonomic compositional dynamic of two types of microbial habitats, free-living (water), and host-associated (skin and gut), were characterised throughout the young developmental stages of the host using a 16S SSU rDNA gene, based metataxonomics approach [40]. Being able to cope with polymetallic gradients generated by acid mine drainages (AMD), the Yellow Perch is a well-established ecotoxicology vertebrate model: many studies have measured the impact of heavy metals on their innate immune system [41], metabolism [42], development [43], parasitism [44,45,46] and gene expression [47]. As these host functions are closely related to gut microbiota composition, the Yellow Perch is a convenient model to test the effects of metal exposure on a vertebrate host-microbiota system. Impact of Cd exposure on microbial communities was documented for AMD bacterioplankton [48], but not within the host-microbiota system. Using non-lethal doses of Cd as a selection pressure, we sought to unravel the neutral and non-neutral processes shaping the microbiota assembly, without triggering significant physiological damage or causing host death.

Materials and methods

Fish rearing

After an acclimation period of 30 days in two containers of 1500 L each, the Yellow Perch (1200 juveniles) were reared in 24 tanks (50 fish per tank) of 36 L and acclimated for the second period of 15 days. Each tank was equipped with an independent filtering system circuit. The fish juveniles were fed daily with the same food from the beginning to the end of the experiment. A second acclimation period of 2 weeks was carried out before the start of Cd exposure (Supplementary file 1).

Cadmium exposure regimes

Control (Control) and Cd (cadmium) treated tanks were randomly distributed in the aquarium facility. The experiment was designed for two Cd exposure regimes (8 tanks per regime), and one negative control regime (8 tanks). In treated tanks, fish were exposed to Cd chloride (CdCl2) provided by Sigma-Aldrich (> 99.9% purity). The Cd was dissolved in ddH2O to target stock concentrations (9 ppb). For the regime of Cd constant concentration (CC), the Cd chloride was initially added at 0.8 ppb (parts per billion), before gradually increasing the concentration every 5 days to reach a maximal concentration at the end of the first month (T1). This maximal concentration was maintained 2 months until the end of treatment (third month, T3). For the regime of Cd variable concentration (CV), the Cd concentration was gradually increased every 5 days to reach the target concentration (9 ppb) at the end of treatment (third month, T3). The maximal CdCl2 concentration was set at 9 μg/L as it is the highest CdCl2 concentration tolerated by Yellow Perch in contaminated Canadian lakes [49].

Host-microbiota and water sampling

A total of 432 mucosa host-microbiota samples were collected for this study, 216 (3 times × 3 regimes × 8 tanks × 3 replicates) skin mucus swabs and 216 (3 times × 3 regimes × 8 tanks × 3 replicates) gut tract samples (Supplementary file 2). Water samples were stored in sterile bottles (Nalgene), 2 l per tank were filtered using a polycarbonate membrane of 0.22 μm. In total, 144 filters (3 times × 3 regimes × 8 tanks × 2 replicates) were conserved in 2 mL sterile microcentrifuge tube and directly stored at − 80 C.

Metal concentration in water and fish liver

Concentrations of metal traces (Cd, Cu and Zn) within the water and liver were determined with the ICPMS (Ionization Coupled Mass spectrometry) technology at the Department of Chemistry, Laval University for T0 and T1, then at INRS (Institut National de la Recherche Scientifique), Quebec, for T1-T3. For water, before ICPMS analysis of Cd ions, the CdCl2 in water samples (10 ML tubes) was fixed by adding 4% of nitric acid. The metals ions (Cd, Cu, Zn) concentrations were measured in water every week until the end of the CdCl2 exposure regimes. For fish, liver samples after lyophilisation were digested with purified nitric acid and kept at room temperature for 5 days. The liver acid digestion protocol was adapted from Pierron et al. [50]. For further details, see the Supplementary file 2. The metal concentrations were analysed using a two-way analysis of variance (ANOVA) of two independent factors: time and treatment. The interactions between time and treatment factors (the interaction means that the effect of treatment depends on time) were analysed using Tukey’s test and Wilcoxon rank test depending on the data (metal concentration) after assessing the normality of data distribution.

DNA extraction to Illumina Miseq sequencing

DNA was extracted from all skin mucus and water samples using the Qiagen DNeasy Blood and tissue kit (Supplementary file 2). For all intestine samples, after an RNA extraction for a transcriptomic project, the DNA was extracted from TRIzol organic phase using BEB (back extraction buffer) and PCI (phenol/chloroform/isoamyl alcohol 25:24:1) solution (Supplementary file 3). The 16S ribosomal DNA was amplified via PCR using universal primers specific to the V3-V4 hypervariable region of the rDNA 16S gene [51]. The purified product of first-round PCR was used as a template for the library preparation by performing second-round PCR. Final amplified DNA was verified by electrophoresis on 2% agarose gel, and finally, DNA concentration of the product was quantified by fluorescence using Quant-iT™ PicoGreen™ dsDNA Assay Kit (Thermo Fischer Scientific).

Analysis of 16S rDNA amplicons

Sequence analysis was performed with our bioinformatic pipeline as described previously [27, 40]. The code of our pipeline is available on Github (https://github.com/BachBioinformatics/MicrobiomePipelines). After the construction of OTUs table, a decontamination step was performed using BWA mapper [52] implemented in DeConseq tool [53], which consisted of mapping the OTUs sequences against the draft genome of Perca flavescens available on NCBI (98% identity threshold). The analysis of the alpha and beta diversity of metacommunities was proceeded using the Rhea package [54]. Afterwards, the significance of variations in alpha-diversity indexes (richness and evenness) and beta-diversity (phylogenetic distance) divergence between experimental groups was assessed using pairwise and multiple rank statistics tests (Wilcoxon /Kruskal-Wallis). Beta-diversity was measured using generalised UniFrac distance [55], which considers both dominant and rare OTUs. P-values of pairwise comparisons in alpha and beta-diversity were validated with multiple correction tests [56] using the B-H (Benjamini-Hochberg) for avoiding the Type I errors (false positives).

Correlational networks analyses

The Spearman coefficients supported with multiple corrections test Benjamini-Hochberg (BH) were computed to measure the correlations between OTUs. This coefficient was recently demonstrated as a robust approach in terms of sensitivity and precision of correlation detection [57]. Only strong correlations (p-value BH < 0.05) positive (Corr > 0.6), and negative (Corr < − 0.6) were visualized in the OTU networks. Cytoscape software [58] was used to perform network visualization and analysis. The number of components indicated the fragmentation (number of subnetworks) in each community per condition. Each node size in the network was proportional to the average OTU’ relative abundance in all the samples per group. Although Spearman and new methods of networks inference [57] are comparable, possible bias of indirect associations of nodes can be generated by spearman coefficient [59]. To deal with this bias, the SPIEC-EASI (SParse Inverse Covariance Estimation for Ecological Association Inference) method [59] was used to supplement the OTUs networks built with spearman coefficient. In SPIEC-EASI networks, to avoid false-negative correlations, the OTUs which occur in less than three samples were combined in one synthetic OTUs. To compare the connectivity of networks between different groups, the average of nodes degree was compared using the Kruskal-Wallis test. Also, the modularity of different networks was assessed with the Markov Clustering Algorithm (MCL) using a value of 2.5 as a reasonable inflation parameter (granularity) of clustering [60]. Modularity is a measure of network structure that was designed to measure the strength of division of a network into modules. Networks with high modularity have dense connections between the nodes within modules but sparse connections between nodes in different modules [61].

Metacommunity assembly modelling

To investigate the role of neutral processes in community assembly, we fit the distribution of OTUs to a neutral model of microbial assembly [16] using a non-linear least-squares approach and beta distributions, which has recently been implemented by Burns and colleagues [22]. The neutral model compares the frequency of OTU occurrence to their abundance in the metacommunity by estimating a parameter (m), which represents the migration rate and which can be interpreted as a measure of dispersal limitation (low migration rate means high dispersal limitation) [22]. The estimated migration rate (m) is the probability that a random loss (death or emigration) of an OTU in a local community is replaced by dispersal from the metacommunity source [22]. The temporal comparisons (T0-T1; T1-T3) of predicted versus observed OTU frequencies from the neutral model were used to highlight the percentage of OTUs fitting the model with a confidence interval of 95%. The goodness of fit to the neutral model was assessed using R-square as the coefficient of determination. R-squared equal or higher to a value = 0.5 is an acceptable threshold of the goodness of fit to Sloan’s neutral model [16]. In addition to Sloan’s model, to quantify the ecological stochasticity in the community structure, the new proposed index, normalized stochasticity ratio (NST) was computed to estimate the ecological stochasticity based on beta diversity [62]. The NST index was used here to estimate the average stochasticity between local communities (treatments and control) using Jaccard and Ruzicka metrics (recommended by the authors of NST model) and “fixed proportion” as null-model which means that the occurrence probability of a taxon is proportional to its observed occurrence frequency. The NST values were computed using the function “tNST”, then because pairwise comparisons of observed/null dissimilarity values are not independent, bootstrapping analysis based random comparisons (N = 1000) for stochasticity ratio was also performed.

Results

Metal concentrations, diversity measures and taxonomic composition

In water, the difference in Cd concentration is significant between all treatments (CC, CV or Control) at T1 and T3. In the fish liver, the difference is only significant at T3 (Supplementary Table 1). For alpha-diversity, the treatment reveals a significant effect only on the evenness (Shannon effective) in the skin mucosal communities at T1, and on the richness in the gut communities at T3 (Figs. 1a-b and 2a-b). The Cd treatment reveals also a significant imbalance in the taxonomic composition in the skin microbiome at T1 (Table 1). The statistical comparisons using both diversity indices for evenness and richness demonstrate the importance of time as a driver in microbial community richness, evenness (Fig. 1a, Supplemetal Figure 1, Supplemetal Table 2), and taxonomic composition (Supplemetal Figures 2 and 3) rather than treatment. For beta-diversity (GUniFrac distance), by T1, significant differences (BH p-value < 0.05 of significant PERMANOVA tests) among treatments are only observed in the mucosal skin communities (Supplemetal Table 3, Fig. 1b-c); however, by T3, both cadmium exposure regimes indicate significant changes in both skin and gut microbial communities compared to the control (Supplementary Table 3, Fig. 2b-c). Among water microbial communities, the phylogenetic distance between treatments and control is significantly different at every time points (BH’s p-value of PERMANOVA tests for CC-Ctrl T1: 0.0135; CC-Ctrl T3: 0.002; CC-CV T1: 0.136; CC-CV T3: 0.0015; CV-Ctrl T1: 0.006; CV-Ctrl T3:0.0015). Interestingly, the phylogenetic distance between CV and CC treatments becomes significant at T3 (Supplementary Table 3; Supplementary Figure 4).

Fig. 1
figure1

Diversity, structure, and assembly of the skin microbiome. a highlights the negative correlations in the co-abundance networks of skin microbial community. Each node size in the network is proportional to the average of the OTUs relative abundance in all samples. These networks are based on significant Spearman coefficients and were constructed using R scripts and Cytoscape software. b shows the boxplots of the significant difference in alpha-diversity (Shannon effective). c summarises the statistical tests of alpha-diversity and beta-diversity (Gunifrac distance). d represents the PCOA 3D plot of the microbial skin communities of all treatment groups based on the significant difference of generalized Unifrac distances tested with PERMANOVA, MRPP and multiple correction test (See Supplementary Table 3) e reports the distribution of neutrality versus abundance cut-off and goodness of fit. The plot shows the variation of neutral OTUs percentage (Y-axis) along with the goodness of fit predicted by NLS models using 12 cut-offs of relative abundance averages (facet panels) in the skin metacommunity at T1 and T3

Full size image
Fig. 2
figure2

Diversity, structure, and assembly of the gut microbiome. a highlights the peripheral location of the most abundant OTUs in the correlation networks of the gut microbial community. Each node size in the network is proportional to the average of the OTUs relative abundance in all samples. These networks are based on significant Spearman coefficients and were constructed using R scripts and Cytoscape software. b shows the boxplots of the significant difference in alpha-diversity (Shannon effective). c summarizes the statistical tests of alpha-diversity and beta-diversity (Gunifrac distance). d represents the PCOA 3D plot of the microbial skin communities of all treatment groups based on the significant difference of generalized Unifrac distances tested with PERMANOVA, MRPP and multiple correction test (See Materials and Methods and Supplementary Table 3) (e) reports the distribution of neutrality versus abundance cut-off and goodness of fit. The plot shows the variation of neutral OTUs percentage (Y-axis) along with the goodness of fit predicted by NLS models using 12 cut-offs of relative abundance percentages (facet panels) in the gut metacommunity at T1 and T3

Full size image
Table 1 Statistical summary of taxonomic significant changes between treatments
Full size table

Overall, the effect of treatment on the host-microbiota for beta diversity is more significant than time (i.e. microbiota ontogeny) for skin (T1, T3) and gut (T3) microbiota. Contrastingly, for alpha diversity, the effect of time is more significant than treatment. To explore the effect of Cd treatment on the interactions within and between host and water microbial communities we analyzed the networks parameters mainly degree and modularity.

Substantial role of rare taxa in microbiome network connectivity

Within gut microbial communities, the high abundant OTUs are peripheral (with minimal interactions) or disassortative (not connected) to the giant network component, see [63], suggesting low overall connectivity (Fig. 2d). Within skin and gut microbiome networks, no correlation was found between the degree (number of connections per node) and the average relative abundance of OTUs (nodes size in the network). However, given that the high abundance is a feature of few OTUs, most of the connections within host-microbiome networks occur among rare or low abundant OTUs (relative abundance < 0.1%) – especially in skin communities early in both Cd exposure regimes a T1 in the networks of Cd treatment groups (Fig. 1d).

Poor connectivity in gut microbiome network reflects cd selection regimes

The exposure to Cd has also a significant impact on network connectivity and integrity in the gut microbial communities. In the control group, most abundant OTUs were connected to a central hub (Fig. 2d). However, in the Cd-treated groups, the most abundant OTUs (> 1%) were gradually disconnected from the main network in small independent hubs or sub-networks. Considering the MCL clustering of gut microbiome networks into modules, the modularity indicates higher values at time T1 in Cd treatment groups (Ctrl modules = 12; CV modules = 41; CC modules = 44). However, at time T3 the network’s modularity is higher in the control group (Ctrl modules = 94; CV modules = 50; CC modules = 65) (Table 2; Supplementary Table 4). In term of nodes connectivity, the average degree was higher in the control group compared to Cd treatment groups at T1 and T3 (Fig. 3; Supplementary Table 4).

Table 2 Host microbiome networks modularity
Full size table
Fig. 3
figure3

The average degree of host-microbiome networks over time and between treatments. The connectivity represented with violin plots are significantly higher in average for the control groups in the gut and the skin at time T3, whilst at T1, they were significantly higher for cadmium-treated groups only in the skin. The average of nodes’ degree computed with Network Analyzer was compared using the Kruskal-Wallis test followed by Benjamini-Hochberg test. The value of 0.05 is the threshold of B-H p-value significance. Only the significant Dunn test p-values for pairwise comparisons are displayed on this figure, however in order to improve visibility, the significant p-values of Kruskal-wallis test for multiple groups comparisons were not plotted

Full size image

Negative correlations in the skin microbiome network suggest dysbiosis under disturbance

In skin microbial community, the correlational networks at time T1 are characterised by high modularity (Ctrl modules = 12; CV modules = 41; CC modules = 42) (Table 1; Supplementary Table 4) and a significantly high average degree in the Cd treatment groups (Fig. 3). However, at time T3, like in the gut microbiome, the skin microbiome networks indicate higher modularity in the control group (Ctrl modules = 42; CV modules = 28; CC modules = 30), and significantly high average degree. The higher average degree connectivity observed in CV and CC relative to the Control treatment at T1 is manifest in the significantly high percentage of significant negative correlations (red edges in networks of Fig. 1d) observed in those groups (CV neg. corr.: 6.9%; CC neg. corr.: 6.3%) compared to the control group (Control neg. corr.: 2.2%). A significant increase over time (T0-T1) in the abundance of Tenericutes (T), and Bacteroidetes (B), and decrease of Firmicutes (F),and Proteobacteria (F) (Supplementary Table 5) is associated with nodes implicated with negative correlations in the control group Ctrl (T:6% B: 3% F: 6% P: 63%), CV (T:6% B:7% F: 19% P: 57%) and CC (T:8% B: % 7 F: 13% P: 61%). In fact, at the time T1, the Bacteroidetes have significantly increased over time in the skin for CC & CV groups (Supplementary Table 5), not for control, showing a significantly higher abundance in skin compared to water and gut microbial communities for Cd groups (Supplementary Figure 5). Like Bacteroidetes, few nodes are exclusively implicated in negative correlations in CC and CV for Synergistetes (~ 1%), Acidobacteria (~ 1%), and Deinococcus-Thermus (~ 1%). Interestingly, however, Euryarchaeota (E) is more implicated in negative correlations within the control group Ctrl (10%), than CV (2%) and CC (1%).

Fragmentation of water microbiome network explains community dynamics

The water microbiome networks are fragmented, both over time and among treatments (Fig. 4; Supplementary Table 4). There was a notable lack of Tenericutes (Mycoplasma) in comparison to the skin and intestinal microbial communities’ networks. The degree average is significantly high in the control group of water microbiome networks at T1 and T3 (Fig. 5). Overall, the structure of the water networks encompassed small, disconnected worlds of independent hubs.

Fig. 4
figure4

Correlational co-abundance networks of water microbial community. Water microbial community networks displayed fragmented interactions, over time and between treatments — the topology of water microbial encompassed sub-networks disconnected in small independent hubs. The number of independent hubs was the highest in Control network at T0 and T3, at T3 in CC, and it was always intermediate in CV network (see the Supplementary Table 3). These networks are based on significant Spearman coefficients and were constructed using R scripts and Cytoscape software

Full size image
Fig. 5
figure5

The average degree of water microbial networks over time and between treatments. The connectivity represented with violin plots is significantly higher on average for the control groups in the gut and the skin at times T3 and T1. The average of nodes’ degree computed with Network Analyzer was compared using the Kruskal-Wallis test followed by Benjamini-Hochberg test. The value of 0.05 is the threshold of B-H p-value significance. Only the significant Dunn test p-values for pairwise comparisons are displayed on this figure, however in order to improve visibility, the significant p-values of Kruskal-wallis test for multiple groups comparisons were not plotted

Full size image

Metacommunity dynamic indicates niche specialisation and time effect

The dynamics of interactions of water microbial communities with the different host (skin and gut) microbiomes reveal the similar structure of networks between treatments, but not overtime. At T1, the subnetwork assembling gut nodes is weakly connected to the modules of the water and skin microbial subnetworks. A T3, the three microbial subnetworks of water, gut and skin are almost disconnected within the Cd treatment groups and weakly connected in the control group (Fig. 6). The comparison of connectivity assessed with nodes degrees and represented with violin plots (Supplementary Figure 6) indicate that the average of connections is significantly higher in the skin compared to the water and gut microbiomes in all treatments and at all time points. However, at T1 the connectivity converged (which means not significantly different) between the water and skin microbiome for cadmium-treated groups. The comparison of centrality assessed with nodes closeness centrality and represented with violin plots (Supplementary Figure 7) indicates that the average of connections is significantly higher in the skin compared to the water and gut microbiomes in all treatments at T0 before disturbance. However, at T1 the average node centrality converged (which means not significantly different) between the water and skin microbiome for cadmium-treated groups. At T3, the centrality of node converged between the gut and skin microbiome in the control group. Finally, from a methodology point of view, the networks built with SPIEC-EASI method and Spearman coefficient converge to similar topologies (see the gut and skin microbiome networks built with SPEIC method in Supplementary Figure 8).

Fig. 6
figure6

Statistical analysis of the node degree in the water and host-microbiome networks. The comparison of connectivity assessed with nodes degrees and represented with violin plots indicate that the average of connections is significantly higher in the skin compared to the water and gut microbiomes in all treatments and at all time points. However, at T1 the connectivity converged (which means not significantly different) between the water and skin microbiome for cadmium-treated groups. The average of nodes’ degree computed with Network Analyzer was compared using the Kruskal-Wallis test followed by Benjamini-Hochberg test. The value of 0.05 is the threshold of B-H p-value significance. Only the significant Dunn test p-values for pairwise comparisons are displayed on this figure, however in order to improve visibility, the significant p-values of Kruskal-wallis test for multiple groups comparisons were not plotted

Full size image

Stochasticity and determinism involved in water and host-associated bacterial community assembly

In host and water bacterial communities, considering all the OTUs without any filtration of abundance, the occurrence frequency of the majority of OTUs in the Gut (CCT1:T3 = 93:85%; CVT1:T3 = 94:88%; CtrlT1:T3 = 99.6:99.4%), the skin (CCT1:T3 = 83:91%; CVT1:T3 = 86:88%; CtrlT1:T3 = 99.4:99.3%) and the water (CCT1:T3 = 81:99%; CVT1:T3 = 75:88%; Ctrl T1:T3 = 97:93%) was confidently predicted by the NLS model (Supplementary Table 6, Figs. 1e and 2e). According to this model, the percentage of neutral OTUs and the goodness of fit (R2) are always higher in the control group compared to both Cd-treated regimes in water and host communities (Supplementary Table 6). The percentage comparison of neutral OTUs between different models could not quantify the relative impact of the non-neutral OTUs (overfit or underfit the model) because the goodness of fit and the estimated emigration rate with maximum likelihood (m.mle) varied between water and host models (Supplementary Table 6). The emigration rate in the NLS model predicts the stochastic demography in community assembly. So, at T1, the comparison of the estimated emigration rates indicates that stochasticity was higher in the water microbiome compared to host microbial communities in the Control (GutT1 = 0.172; SkinT1 = 0.200; WaterT1 = 0.326), CV (GutT1 = 0.118; SkinT1 = 0.204; WaterT1 = 0.423) and CC (GutT1 = 0.129; SkinT1 = 0.242; WaterT1 = 0.335). Similarly, at T3, the migration rate indicates higher stochasticity in water compared to the host-microbial communities for the Control (GutT3 = 0.269; SkinT3 = 0.177; WaterT3 = 0.488), and CV (GutT3 = 0.351; SkinT3 = 0.138; WaterT3 = 0.424) groups, but not for the CC (GutT3 = 0.564; SkinT3 = 0.169; WaterT3 = 0.362). Investigating the relationship of OTU abundance and neutrality, our analysis indicates that the percentage of neutral OTUs decreased in the Ctrl groups compared to Cd-treated groups, in host and water microbial communities at T1 and T3, especially when the low abundant OTUs were discarded from the NLS models (see the neutral OTUs percentage and goodness of fit across OTUs abundance cut-offs, Figs. 1e and 2e, Supplementary Figure 9, Supplementary Table 6).

To quantify the ecological stochasticity in the community structure, we applied a null modelling approach using normalized stochastic ratio (NST) index (see Materials and methods). When NST is higher than 0.5, the community structure is more likely driven by the neutral processes, whilst when the NST is lower, the main driving processes are deterministic. According to the results, at each time point, and in all treatments, the average of NST (NSTavg.) for almost all communities is higher than 0.5 (Fig. 7). Comparing stochasticity ratios between water and host microbiota, the NST modelling lead to same results of NLS models. The NSTavg. values indicate stochasticity equal (at time T0) or higher (at T1 and T3) in water compared to the gut and skin microbiota. When it comes to comparing the magnitude of stochasticity between treatments and control in each community type structure, both NST and NLS models agree in terms of the dominance of neutral processes in the community assembly (goodness of fit R2 > 0.5), and structure (NSTavg > 0.5), although giving further insights in both demographical and ecological aspects, respectively. With NLS model, the percentage of neutral OTUs and the goodness of fit are always higher in the control group, suggesting that a higher proportion of non-neutral OTUs influenced demographical dynamics in Cd-treated groups relatively to control group. With NST models, the ecological stochasticity in the community structure was significantly higher (NSTavg > 0.5) that the null hypothesis in all groups. Furthermore, higher averages of ecological stochasticity were observed under perturbations in Cd-treated groups early at T1 in the skin microbiome and lately at T3 in the gut microbiome (Fig. 7).

Fig. 7
figure7

Comparative analysis of stochasticity ratios averages between treatments. The NST values were computed using the NST model available in the R package. With this model the” Proportional fixed” null-model, Jaccard similarity coefficient (also known then as Ruzicka similarity) and n = 1000 for random permutations were used to compute the normalized stochasticity ratios. The NST in each group represented in violin plots were compared for their average using Kruskal-Wallis test followed by Benjamini-Hochberg as a multiple correction test. Only the significant Dunn test p-values for pairwise comparisons are displayed on this figure, however in order to improve visibility, the significant p-values of Kruskal-wallis test for multiple groups comparisons were not plotted

Full size image

Discussion

Our study evidenced salient differential changes in community assembly across three community types (environmental or host-associated), time and xenobiotic exposure regimes. This study for the first time highlights the relative contribution of neutral and non-neutral factors in shaping the microbiota during the early life-stages of an ecotoxicology vertebrate model. First, alpha diversity increased significantly in the skin and gut microbiota in both constant and gradual selection regimes. Then, at the community-level, significant phylogenetic divergence was observed between the control and treatment groups in the three community types at two specific time points, T1 and T3. These two key time points were investigated further with co-abundance network analysis and community assembly modelling. Frequent significant negative correlations between taxa in both selection regimes in the skin and the increasing richness of environmental bacterial strains suggest a dysbiosis in the mucosal host-associated microbiota [64, 65]. Negative correlations between OTUs occur is all groups at time T1, but they are remarkably more frequent in both selection regimes within the skin microbiome networks. In many studies, the increase of negative correlations in networks was associated to microbiota dysbiosis in several contexts of human lung infectious diseases [66,67,68], in dog’s gut (Vázquez-Baeza et al. [69]) and macaque model of tuberculosis for instance. Furthermore, network connectivity under stress was maintained by rare OTUs, while abundant OTUs were mainly composed of opportunistic invaders such as Mycoplasma and other genera related to fish pathogens like Aeromonas, Pseudomonas and Flavobacterium [70, 71].

Lastly, to predict the nature of evolutionary processes driving the metacommunity assembly under Cd perturbations, the fit of the neutral model, based non-linear least squares (NLS), indicates that neutral processes predominantly drove taxonomic drift during the assembly of host and water microbial communities in all experimental groups (R2 > 0.5). However, regarding the percentage of neutral OTUs between groups suggests significantly lower stochasticity of community assembly in Cd-treated groups. At this point, from the higher percentage of non-neutral OTUs, we propose that Cadmium perturbations have an effect not only on the community assembly but also on the community structure. Previous studies quantified high ecological stochasticity in the community structure under perturbations by suggesting specific [72] and general indexes [62] of stochastic strength. The application of NST model to our data (NSTavg > 0.5) confirms high ecological stochasticity in the community structure changes under disturbance, as suggested by Ning and co-authors [62], and this result agrees with the dysbiosis and the model of resistance to colonisation. However, we disagree with the Ning et al. (2019) [62] on the comparison of the NLS model with NST indexes, because this comparison is meaningless. The general NST null model was developed for quantitatively assessing ecological stochasticity, while the NLS model Sloan’s neutral model was developed for quantitatively assessing the demographic stochasticity. Initially, the NST model was developed and applied to estimate the stochasticity strength during the succession of a groundwater microbial community in response to perturbations caused by organic carbon (vegetable oil) injection [62]. The authors of NST have also compared their model index with not normalized stochasticity strength measure (an early version of NST) [72] and the percentage of neutral OTUs predicted by the NLS of Sloan’ model [22]. NST demonstrated high stochasticity under perturbation (oil injection), and in our host-microbiome system, we found the same trends. On the other side, the NLS model has different assumptions based on the estimation of the emigration rates (can be interpreted as dispersal limitation index) from a community source but not based on an ecological similarity/dissimilarity distance. To our knowledge, with the NST model, we could not distinguish between neutral and non-neutral OTUs or quantify a migration rate to assess the dispersal ability from the community structure. From our point of view, the comparison of NST with NLS could not infer the relative contribution of deterministic processes. The NST authors ‘model also acknowledged in their discussion that the operational distinction of stochasticity and determinism can appear somewhat arbitrary with their NST index, and it is difficult to distinguish ecological stochasticity from the noise caused by deterministic environmental factors, as shown in their simulation [62].

Last not least, we would stress that NLS model, with our null hypothesis based on a real bacterial community (i.e. the control community), detected non-negligible proportion of non-neutral OTUs, potentially suggesting that host-microbiota assembly partly resulted from deterministic processes. Although moderate, the percentage of OTUs not fitting the NLS model in both Cd regimes may suggest patterns of deterministic effects on microbiota recruitment [27] host filtering and Cd treatment. The host filtering can be interpreted lately at T3 throughout the network’s modularity and average connectivity in the control group. The high modularity and degree observed at T3 would approach specialized, partitioned microbiome networks in the gut and skin [73, 74], reflecting niche differentiation [75] and microbiota maturation along with perch juvenile’s development [3]. However, in an opposite direction, the gut and skin networks modularity and degree decreased over time to reach a lower average in CC & CV groups, suggesting broader niche partitioning [74] and lesser interactions and more stochastic changes under disturbance [62].

The dominant ecological stochasticity quantified under disturbance left interesting patterns in the community structure of the skin microbiota. At T1, significant phylogenetic divergence occurred between treatments (CC & CV) and control (Control). Most importantly, a taxonomic convergence between treatments (CC & CV) not only in the skin but also within the water community occurred. This convergence mainly resulted from an invasion of environmental bacterial strains in the skin. The significant increase of Bacteroidetes only in treatments groups (CC &CV) in skin communities compared to water communities could strongly support this hypothesis. Such gain or loss of tissue-specific community type suggests a disruption of the host’s ability to control the assembly of skin microbiota, which is correlated with Cd exposure. This phenomenon is termed “direct colonisation resistance” [76]; however, we do not exclude that this colonisation failure also resulted, at least partly, from a host immune system failure (termed as “immune colonisation resistance”). This compositional disruption translated into many negative correlations between taxa in both selection regimes in the skin-associated microbiota at T1. Furthermore, the impact of Cd exposure on skin community structure was also observed at T3, where the phylogenetic distance became significantly divergent, even between both selection regimes, where many negative correlations were detected between taxa. In addition to the increasing invasion of environmental bacterial strains in the skin (i.e. failure of colonisation resistance, see [76]), and the rise of negative correlations suggest a dysbiosis state of skin-associated microbiota [66,67,68,69, 77]. This dysbiosis might be associated not only with an increase in evenness and phylogenetic convergence with the water bacterial community but also with the rise of antagonism among OTU co-abundance networks in both selection regimes (CC & CV). Most of antagonism was mediated through rare and abundant Tenericutes (Mycoplasma) and Proteobacteria. Depending upon the strain, Mycoplasmas are thought to be either fish opportunists, or innocuous commensals in fish [78, 79]. In comparison to the skin, the significant divergence between control (Control) and treatments (CC, CV), and the rise of negative correlations, appeared later in the gut community: at T3. This delayed pattern of dysbiosis strongly suggests that the physiological impact of cadmium exposure was mitigated more effectively within the gut. In fish (and other vertebrates), the liver is the main organ to accumulate xenobiotics including cadmium. Therefore, the late compositional change in the gut microbiota potentially occurred when bioaccumulation of cadmium within the liver reached its maximum carrying capacity. Another noticeable compositional change was the disconnection of abundant taxa from the main gut interacting network, which was proportional to the stress intensity. By T3, the overall taxonomic network connectivity was formed exclusively from rare OTUs. Contrastingly, abundant OTUs were peripherals and disconnected from central hubs [63], mainly composed of putative opportunistic invaders such as Mycoplasma and other genera encompassing strains associated to fish pathogens, like Bacillus (> 6% in CC and CV). As observed in other fish species such as Atlantic salmon (Salmo salar) [80] and the long jaw mudsucker (Gillichythys mirabilis) [81, 82], Yellow Perch have intestinal microflora dominated by Tenericutes (Mycoplasma sp.). It is therefore difficult to conclude whether the increase of several Mycoplasma strains is beneficial or not to the host. Concerning Bacillus, a similar increase was associated with irritable bowel disease (IBD) in dogs and negatively correlated with bacterial strains associated with healthy individuals [69]. Interestingly, rare OTUs and negative correlations did not play an essential role in water co-abundance networks, which were highly fragmented in Cd treatment regimes, and the average connectivity was significantly higher in the control regime. It is worthy to notice that the taxonomic composition of water bacterial communities was characterized by a very low occurrence of Tenericutes (Mycoplasma) compared to host communities. The metacommunity networks combining host and water bacterial communities highlights the niche divergence and dynamics between the water bacterial communities and each of the host microbiota. A T0, all these community types showed interactions within a distinctive structure of subnetworks, which differentiated over time in independent hubs weakly connected (better connected in the control group).

Finally, low abundant or rare OTUs have been demonstrated to play a pivotal role in community assembly [35,36,37,38,39] either in promoting homeostasis [37] or dysbiosis [83]. Therefore, we applied the NLS model to disentangle neutral and non-neutral evolutionary processes that were at play for rare OTUs. The variation of the presence of neutral OTUs across different abundance cut-offs (Supplementary Table 6) suggested a neutral role around 35 to 50% of rare OTUs in host water communities. Overall, our data demonstrated that stochastic processes mainly drove taxonomic drift even under Cd disturbance. However, host-microbiota assembly evolved by involving non-neutral processes in cadmium treatment groups, although this trend was less salient in both experimental groups at T3. The majority of the OTUs that did not fit the neutral model was assigned to Mycoplasma genus under and after the exposure during the recovery period [40]. The assembly of gut microbial communities may have evolved under non-neutral processes due not only to the cadmium as a disrupting factor but also due to the selection imposed by the host development [22].

Conclusions

In this experimental evolution study, our findings demonstrate the extensive involvement of low abundant (rare) taxa throughout community assembly and interacting network connectivity under perturbations. We have unearthed a niche-specific response to cadmium disturbance. In the skin microbiome network, we detected early signals of antagonistic interactions by a preponderance of negative correlations, whilst at the same time, the overall connectivity in the gut microbiome was degraded over time. The network topology in the control group suggests specialized, partitioned microbiome interactions in the gut and skin, reflecting niche differentiation and microbiota maturation along with perch juveniles’ development. On the other hand, the network topology for Cadmium-treated groups suggests a broader niche, fewer interactions. Considering the modelling frameworks, neutral processes were thus the major forces of the community assembly in the environmental microbiome, although involving a low percentage of deterministic changes in the host-microbiota under disturbance. The taxonomic convergence between water- and skin-associated bacterial communities across both cadmium exposure groups highlights the loss of the colonization resistance capacity of the host. This was due to physiological stress experienced by the host: cadmium bioaccumulation in Perch’s liver has already been documented to disrupt host physiology. By highlighting the link between a loss of colonization resistance and dysbiosis within the host (which in turn is known to induce an inflammatory response), our results will be useful not only for the field of microbial ecology but also for biomedical research, as dysbiosis of gut microbiome composition has been shown to result in the onset of various inflammatory diseases such as diabetes, IBD, Crohn disease, cancer, and obesity [84,85,86].

Perspective

The patterns of neutral and non-neutral assembly in contrasting types of bacterial communities (i.e. one environmental and two host-associated) described here provide novel key insights regarding our understanding of evolutionary forces that are at play in shaping the host-microbiota when facing sublethal environmental stress. Living organisms are currently facing unprecedented levels of environmental stressors that impact their capacity to cope with natural pathogens, essentially by altering their overall immune defence. Therefore, there is an urgent need to accurately decipher the early warning signals occurring at the first stages of xenobiotic exposure.

Availability of data and materials

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

Abbreviations

ANOVA:

Analysis of variance

16 s rDNA:

16S Ribosomal DNA

BEB:

Back extraction buffer

BH:

Benjamini-Hochberg correction test

CdCl2:

Cadmium Salts

CPAUL:

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

Ctrl:

The regime of negative Control

CV:

The regime of Cadmium Variable Concentration

CC:

The regime of Cadmium Constant Concentration

IBD:

Irritable bowel disease

m:

Emigration rate

m.mle:

Maximum likelihood estimation of emigration rate

NLS:

Non-linear least-squares model

NST:

Normalized stochasticity ratios (NST)

MRPP:

Multiple Response Permutation Procedure

NMDS:

Non-metric Multi-Dimensional Scaling

OTUs:

Operational Taxonomic Units

ppb:

Parts per billion

PCoA:

Principal Coordinates Analysis

PCR:

Polymerase chain reaction

PERMANOVA:

Permutational analysis of variance

RDP:

Ribosomal Database Project

R2 :

Goodness of fit

References

  1. 1.

    Johnson KV-A, Foster KR. Why does the microbiome affect behaviour? Nat Rev Microbiol. 2018;16:647–55.

    CAS  PubMed  Article  Google Scholar 

  2. 2.

    Stephens WZ, Burns AR, Stagaman K, Wong S, Rawls JF, Guillemin K, et al. The composition of the zebrafish intestinal microbial community varies across development. ISME J. 2016;10:644–54.

    PubMed  Article  CAS  Google Scholar 

  3. 3.

    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.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  4. 4.

    Sherrill-Mix S, McCormick K, Lauder A, Bailey A, Zimmerman L, Li Y, et al. Allometry and ecology of the bilaterian gut microbiome. mBio. 2018;9:e00319–8.

    PubMed  PubMed Central  Article  Google Scholar 

  5. 5.

    Theis KR, Dheilly NM, Klassen JL, Brucker RM, Baines JF, Bosch TCG, et al. Getting the hologenome concept right: an eco-evolutionary framework for hosts and their microbiomes. mSystems. 2016;1-2:e00028-16.

  6. 6.

    Faust K, Raes J. Microbial interactions: from networks to models. Nat Rev Microbiol. 2012;10:538–50.

    CAS  PubMed  Article  Google Scholar 

  7. 7.

    Nemergut DR, Schmidt SK, Fukami T, O’Neill SP, Bilinski TM, Stanish LF, et al. Patterns and processes of microbial community assembly. Microbiol Mol Biol Rev. 2013;77:342–56.

    PubMed  PubMed Central  Article  Google Scholar 

  8. 8.

    Vellend M. Conceptual synthesis in community ecology. Q Rev Biol. 2010;85:183–206.

    PubMed  Article  Google Scholar 

  9. 9.

    Tucker CM, Shoemaker LG, Davies KF, Nemergut DR, Melbourne BA. Differentiating between niche and neutral assembly in metacommunities using null models of β-diversity. Oikos. 2016;125:778–89.

    Article  Google Scholar 

  10. 10.

    Bell G. Neutral macroecology. Science. 2001;293:2413–8.

    CAS  PubMed  Article  Google Scholar 

  11. 11.

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

    Article  Google Scholar 

  12. 12.

    Hubbell SP. Neutral theory and the evolution of ecological equivalence. Ecology. 2006;87:1387–98.

    PubMed  Article  Google Scholar 

  13. 13.

    Chase JM, Leibold MA. Ecological niches: linking classical and contemporary approaches. Chicago: University of Chicago Press; 2003.

    Book  Google Scholar 

  14. 14.

    Harris K, Parsons TL, Ijaz UZ, Lahti L, Holmes I, Quince C. Linking statistical and ecological theory: Hubbell #x0027;s unified neutral theory of biodiversity as a hierarchical Dirichlet process. Proc IEEE. 2017;105:516–29.

    Article  Google Scholar 

  15. 15.

    Jayathilake PG, Gupta P, Li B, Madsen C, Oyebamiji O, González-Cabaleiro R, et al. A mechanistic individual-based model of microbial communities. PLoS One. 2017;12:e0181965.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  16. 16.

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

    PubMed  Article  Google Scholar 

  17. 17.

    Stegen JC, Lin X, Konopka AE, Fredrickson JK. Stochastic and deterministic assembly processes in subsurface microbial communities. ISME J. 2012;6:1653–64.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  18. 18.

    Zeng Q, Sukumaran J, Wu S, Rodrigo A. Neutral models of microbiome evolution. PLoS Comput Biol. 2015;11:e1004365.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  19. 19.

    Zeng Q, Wu S, Sukumaran J, Rodrigo A. Models of microbiome evolution incorporating host and microbial selection. Microbiome. 2017;5:127.

    PubMed  PubMed Central  Article  Google Scholar 

  20. 20.

    Foster JA, Krone SM, Forney LJ. Application of ecological network theory to the human microbiome. Interdiscip Perspect Infect Dis. 2008;2008:839501 https://www.hindawi.com/journals/ipid/2008/839501/. Accessed 17 Sep 2018.

    PubMed  PubMed Central  Article  Google Scholar 

  21. 21.

    Taxis TM, Wolff S, Gregg SJ, Minton NO, Zhang C, Dai J, et al. The players may change but the game remains: network analyses of ruminal microbiomes suggest taxonomic differences mask functional similarity. Nucleic Acids Res. 2015;43:9600–12.

    CAS  PubMed  PubMed Central  Google Scholar 

  22. 22.

    Burns AR, Stephens WZ, Stagaman K, Wong S, Rawls JF, Guillemin 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.

    CAS  PubMed  Article  Google Scholar 

  23. 23.

    Ofiţeru ID, Lunn M, Curtis TP, Wells GF, Criddle CS, Francis CA, et al. Combined niche and neutral effects in a microbial wastewater treatment community. Proc Natl Acad Sci. 2010;107:15345–50.

    PubMed  Article  Google Scholar 

  24. 24.

    Morrison-Whittle P, Goddard MR. Quantifying the relative roles of selective and neutral processes in defining eukaryotic microbial communities. ISME J. 2015;9:2003–11.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  25. 25.

    O’Dwyer JP, Kembel SW, Sharpton TJ. Backbones of evolutionary history test biodiversity theory for microbes. Proc Natl Acad Sci. 2015;112:8356–61.

    PubMed  Article  CAS  Google Scholar 

  26. 26.

    Yeh Y-C, Peres-Neto PR, Huang S-W, Lai Y-C, Tu C-Y, Shiah F-K, et al. Determinism of bacterial metacommunity dynamics in the southern East China Sea varies depending on hydrography. Ecography. 2015;38:198–212.

    Article  Google Scholar 

  27. 27.

    Heys C, Cheaib B, Busetti A, Kazlauskaite R, Maier L, Sloan WT, et al. Neutral Processes Dominate Microbial Community Assembly in Atlantic Salmon, Salmo salar. Appl Environ Microbiol 2020;86:e02283-19. .

  28. 28.

    Jeraldo P, Sipos M, Chia N, Brulc JM, Dhillon AS, Konkel ME, et al. Quantification of the relative roles of niche and neutral processes in structuring gastrointestinal microbiomes. Proc Natl Acad Sci. 2012;109:9692–8.

    CAS  PubMed  Article  Google Scholar 

  29. 29.

    Langenheder S, Székely AJ. Species sorting and neutral processes are both important during the initial assembly of bacterial communities. ISME J. 2011;5:1086–94.

    PubMed  PubMed Central  Article  Google Scholar 

  30. 30.

    Maignien L, DeForce EA, Chafee ME, Eren AM, Simmons SL. Ecological succession and stochastic variation in the assembly of Arabidopsis thaliana phyllosphere communities. mBio. 2014;5-1:e00682-13.

  31. 31.

    McCafferty J, Mühlbauer M, Gharaibeh RZ, Arthur JC, Perez-Chanona E, Sha W, et al. Stochastic changes over time and not founder effects drive cage effects in microbial community assembly in a mouse model. ISME J. 2013;7:2116–25.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  32. 32.

    Costello EK, Stagaman K, Dethlefsen L, Bohannan BJM, Relman DA. The application of ecological theory toward an understanding of the human microbiome. Science. 2012;336:1255–62.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  33. 33.

    Schmidt VT, Smith KF, Melvin DW, Amaral-Zettler LA. Community assembly of a euryhaline fish microbiome during salinity acclimation. Mol Ecol. 2015;24:2537–50.

    PubMed  Article  Google Scholar 

  34. 34.

    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  PubMed  Article  Google Scholar 

  35. 35.

    Banerjee S, Schlaeppi K, van der Heijden MGA. Keystone taxa as drivers of microbiome structure and functioning. Nat Rev Microbiol. 2018;16:567–76.

    CAS  PubMed  Article  Google Scholar 

  36. 36.

    Fierer N. Embracing the unknown: disentangling the complexities of the soil microbiome. Nat Rev Microbiol. 2017;15:579–90.

    CAS  PubMed  Article  Google Scholar 

  37. 37.

    Jousset A, Bienhold C, Chatzinotas A, Gallien L, Gobet A, Kurm V, et al. Where less may be more: how the rare biosphere pulls ecosystems strings. ISME J. 2017;11:853–62.

    PubMed  PubMed Central  Article  Google Scholar 

  38. 38.

    Lynch MDJ, Neufeld JD. Ecology and exploration of the rare biosphere. Nat Rev Microbiol. 2015;13:217–29.

    CAS  PubMed  Article  Google Scholar 

  39. 39.

    Pester M, Bittner N, Deevong P, Wagner M, Loy A. A ‘rare biosphere’ microorganism contributes to sulfate reduction in a peatland. ISME J. 2010;4:1591–602.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  40. 40.

    Cheaib B, Seghouani H, Ijaz UZ, Derome N. Community recovery dynamics in yellow perch microbiome after gradual and constant metallic perturbations. Microbiome. 2020;8:14.

    PubMed  PubMed Central  Article  Google Scholar 

  41. 41.

    Dautremepuits C, Marcogliese DJ, Gendron AD, Fournier M. Gill and head kidney antioxidant processes and innate immune system responses of yellow perch (Perca flavescens) exposed to different contaminants in the St. Lawrence River, Canada. Sci Total Environ. 2009;407:1055–64.

    CAS  PubMed  Article  Google Scholar 

  42. 42.

    Couture P, Rajender Kumar P. Impairment of metabolic capacities in copper and cadmium contaminated wild yellow perch (Perca flavescens). Aquat Toxicol. 2003;64:107–20.

    CAS  PubMed  Article  Google Scholar 

  43. 43.

    Azizishirazi A, Dew WA, Bougas B, Dashtban M, Bernatchez L, Pyle GG. Chemosensory mediated behaviors and gene transcription profiles in wild yellow perch (Perca flavescens) from metal contaminated lakes. Ecotoxicol Environ Saf. 2014;106:239–45.

    CAS  PubMed  Article  Google Scholar 

  44. 44.

    Marcogliese DJ, Brambilla LG, Gagné F, Gendron AD. Joint effects of parasitism and pollution on oxidative stress biomarkers in yellow perch Perca flavescens. Dis Aquat Org. 2005;63:77–84.

    CAS  Article  Google Scholar 

  45. 45.

    Marcogliese DJ, Dautremepuits C, Gendron AD, Fournier M. Interactions between parasites and pollutants in yellow perch (Perca flavescens) in the St. Lawrence River, Canada: implications for resistance and tolerance to parasites. Can J Zool. 2010;88:247–58.

    CAS  Article  Google Scholar 

  46. 46.

    Ryman JE, Walleghem JLAV, Blanchfield PJ. Methylmercury levels in a parasite (Apophallus brevis metacercariae) and its host, yellow perch (Perca flavescens). Aquat Ecol. 2008;42:495–501.

    CAS  Article  Google Scholar 

  47. 47.

    Bougas B, Normandeau E, Pierron F, Campbell PGC, Bernatchez L, Couture P. 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.

    PubMed  Article  CAS  Google Scholar 

  48. 48.

    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:869.

  49. 49.

    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. Factors Influencing Tissue Metal Concentrations, Human and Ecological Risk Assessment: An International Journal. 2008;14(1):97-125.

  50. 50.

    Pierron F, Bourret V, St-Cyr J, Campbell PGC, Bernatchez L, Couture P. Transcriptional responses to environmental metal exposure in wild yellow perch (Perca flavescens) collected in lakes with differing environmental metal concentrations (Cd, Cu, Ni). Ecotoxicology. 2009;18:620–31.

  51. 51.

    Werner JJ, Zhou D, Caporaso JG, Knight R, Angenent LT. Comparison of Illumina paired-end and single-direction sequencing for microbial 16S rRNA gene amplicon surveys. ISME J. 2012;6:1273–76.

  52. 52.

    Li H. Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. ArXiv13033997 Q-Bio 2013.

  53. 53.

    Schmieder R, Edwards R. Fast identification and removal of sequence contamination from genomic and metagenomic datasets. PLoS One. 2011;6:e17288.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  54. 54.

    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.

    PubMed  PubMed Central  Article  Google Scholar 

  55. 55.

    Chen J, Bittinger K, Charlson ES, Hoffmann C, Lewis J, Wu GD, et al. Associating microbiome composition with environmental covariates using generalized UniFrac distances. Bioinformatics. 2012;28:2106–13.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  56. 56.

    Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc Ser B Methodol. 1995;57:289–300.

    Google Scholar 

  57. 57.

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

  58. 58.

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

  59. 59.

    Kurtz ZD, Müller CL, Miraldi ER, Littman DR, Blaser MJ, Bonneau RA. Sparse and compositionally robust inference of microbial ecological networks. PLoS Comput Biol. 2015;11:e1004226.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  60. 60.

    Morris JH, Apeltsin L, Newman AM, Baumbach J, Wittkop T, Su G, et al. clusterMaker: a multi-algorithm clustering plugin for cytoscape. BMC Bioinformatics. 2011;12:436.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  61. 61.

    Ji X, Machiraju R, Ritter A, Yen P-Y. Examining the distribution, modularity, and community structure in article networks for systematic reviews. AMIA Annu Symp Proc. 2015;2015:1927–36.

    PubMed  PubMed Central  Google Scholar 

  62. 62.

    Ning D, Deng Y, Tiedje JM, Zhou J. A general framework for quantitatively assessing ecological stochasticity. Proc Natl Acad Sci. 2019;116:16892–8.

    CAS  PubMed  Article  Google Scholar 

  63. 63.

    Hall CV, Lord A, Betzel R, Zakrzewski M, Simms LA, Zalesky A, et al. Co-existence of network architectures supporting the human gut microbiome. iScience. 2019;22:380–91.

    PubMed  PubMed Central  Article  Google Scholar 

  64. 64.

    Johnson KV-A, Burnet PWJ. Microbiome: should we diversify from diversity? Gut Microbes. 2016;7:455–8.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  65. 65.

    Hooks KB, O’Malley MA. Dysbiosis and its discontents. mBio. 2017;8:e01492–17.

    PubMed  PubMed Central  Article  Google Scholar 

  66. 66.

    Bao Y, Li Y, Qiu C, Wang W, Yang Z, Huang L, et al. Bronchoalveolar lavage fluid microbiota dysbiosis in infants with protracted bacterial bronchitis. J Thorac Dis. 2018;10:168–74.

    PubMed  PubMed Central  Article  Google Scholar 

  67. 67.

    Sai Prashanthi G, Jayasudha R, Chakravarthy SK, Padakandla SR, SaiAbhilash CR, Sharma S, et al. Alterations in the ocular surface fungal microbiome in fungal keratitis patients. Microorganisms. 2019;7:309.

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  68. 68.

    Wang Z, Bafadhel M, Haldar K, Spivak A, Mayhew D, Miller BE, et al. Lung microbiome dynamics in COPD exacerbations. Eur Respir J. 2016;47:1082–92.

    PubMed  Article  Google Scholar 

  69. 69.

    Vázquez-Baeza Y, Hyde ER, Suchodolski JS, Knight R. Dog and human inflammatory bowel disease rely on overlapping yet distinct dysbiosis networks. Nat Microbiol. 2016;1:16177.

    PubMed  Article  CAS  Google Scholar 

  70. 70.

    Llewellyn MS, Boutin S, Hoseinifar SH, Derome N. Teleost microbiomes: the state of the art in their characterization, manipulation and importance in aquaculture and fisheries. Front Microbiol. 2014;5-207.

  71. 71.

    de Bruijn I, Liu Y, Wiegertjes GF, Raaijmakers JM. Exploring fish microbial communities to mitigate emerging diseases in aquaculture. FEMS Microbiol Ecol. 2018;94:fix161.

  72. 72.

    Zhou J, Deng Y, Zhang P, Xue K, Liang Y, Nostrand JDV, et al. Stochasticity, succession, and environmental perturbations in a fluidic ecosystem. Proc Natl Acad Sci. 2014;111:E836–45.

    CAS  PubMed  Article  Google Scholar 

  73. 73.

    Quinn RA, Comstock W, Zhang T, Morton JT, da Silva R, Tran A, et al. Niche partitioning of a pathogenic microbiome driven by chemical gradients. Sci Adv. 2018;4:eaau1908.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  74. 74.

    Ponisio LC, Valdovinos FS, Allhoff KT, Gaiarsa MP, Barner A, Guimarães PRJ, et al. A network perspective for community assembly. Front Ecol Evol. 2019;7:103.

  75. 75.

    Liu Z, Cichocki N, Hübschmann T, Süring C, Ofiţeru ID, Sloan WT, et al. Neutral mechanisms and niche differentiation in steady-state insular microbial communities revealed by single cell analysis. Environ Microbiol. 2019;21:164–81.

    CAS  PubMed  Article  Google Scholar 

  76. 76.

    Buffie CG, Pamer EG. Microbiota-mediated colonization resistance against intestinal pathogens. Nat Rev Immunol. 2013;13:790–801.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  77. 77.

    Cadena AM, Ma Y, Ding T, Bryant M, Maiello P, Geber A, et al. Profiling the airway in the macaque model of tuberculosis reveals variable microbial dysbiosis and alteration of community structure. Microbiome. 2018;6:180.

    PubMed  PubMed Central  Article  Google Scholar 

  78. 78.

    Brown RM, Wiens GD, Salinas I. Analysis of the gut and gill microbiome of resistant and susceptible lines of rainbow trout (Oncorhynchus mykiss). Fish Shellfish Immunol. 2018;86:497–506.

    PubMed  Article  CAS  Google Scholar 

  79. 79.

    Holben WE, Williams P, Saarinen M, Särkilahti LK, Apajalahti JHA. Phylogenetic analysis of intestinal microflora indicates a novel mycoplasma Phylotype in farmed and wild Salmon. Microb Ecol. 2002;44:175–85.

    CAS  PubMed  Article  Google Scholar 

  80. 80.

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

    PubMed  Article  Google Scholar 

  81. 81.

    Bano N, DeRae SA, Bennett W, Vasquez L, Hollibaugh JT. Dominance of mycoplasma in the guts of the long-jawed Mudsucker, Gillichthys mirabilis, from five California salt marshes. Environ Microbiol. 2007;9:2636–41.

    CAS  PubMed  Article  Google Scholar 

  82. 82.

    Givens C, Ransom B, Bano N, Hollibaugh J. Comparison of the gut microbiomes of 12 bony fish and 3 shark species. Mar Ecol Prog Ser. 2015;518:209–23.

    Article  Google Scholar 

  83. 83.

    Hajishengallis G, Darveau RP, Curtis MA. The keystone-pathogen hypothesis. Nat Rev Microbiol. 2012;10:717–25.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  84. 84.

    Snedeker SM, Hay AG. Do interactions between gut ecology and environmental chemicals contribute to obesity and diabetes? Environ Health Perspect. 2012;120:332–9.

    CAS  PubMed  Article  Google Scholar 

  85. 85.

    Tilg H, Kaser A. Gut microbiome, obesity, and metabolic dysfunction. J Clin Invest. 2011;121:2126–32.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  86. 86.

    Tilg H, Moschen AR. Microbiota and diabetes: an evolving relationship. Gut. 2014;63:1513–21.

    CAS  PubMed  Article  Google Scholar 

Download references

Acknowledgements

We thank all of the Derome laboratory team for their help during long days of sampling. Thank you to M. Pierre-Luc Mercier for Qubit PCR’s product quantification, and Dr. Brian Boyle for his precious advice in libraries preparations, thank you for all his team members at sequencing platform of IBIS. We thank Louis Bernatchez and Connie Lovejoy laboratories for their help in materials in DNA extraction and water filtering. We also thank Dr. Michel Lavoie for his advice in Cd manipulation. Thanks for my BSc student Sara-Jane Bedard for her help in adjusting PCR protocols. Thanks for Dr. Ciara Keating for her precious comments and suggestions.

Funding

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 Pr. Nicolas Derome.

Author information

Affiliations

  1. Institut de Biologie Intégrative et des Systèmes (IBIS), Pavillon Charles-Eugène Marchand, Université Laval, 1030, avenue de la Médecine, Québec, QC, G1V 0A6, Canada

    Bachar Cheaib, Hamza Seghouani, Katherine Vandal-Lenghan, Pierre-Luc Mercier & Nicolas Derome

  2. Institute of Biodiversity, Animal Health and Comparative Medicine (BACHM), Glasgow, University of Glasgow, Glasgow, UK

    Bachar Cheaib & Martin Llewellyn

  3. School of Engineering, University of Glasgow, Glasgow, G12 8QQ, UK

    Bachar Cheaib

Authors
  1. Bachar Cheaib
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Hamza Seghouani
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Martin Llewellyn
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Katherine Vandal-Lenghan
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Pierre-Luc Mercier
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Nicolas Derome
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

N.D. and B.C. conceived the experiment, B.C and H.S. reared fish tanks and sampling, B.C. and K.V-L. Performed the DNA extraction, B.C performed molecular biology experiments including PCR and libraries preparations. P-L.M. quantified the DNA in libraries, B.C. produced and analysed the results. B.C. wrote the manuscript. M.L and N.D. corrected and revised all the version of this manuscript. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Bachar Cheaib.

Ethics declarations

Ethics approval and consent to participate

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: Table S1.

ANOVA summary of metals concentrations in water tanks and Perca flavescens fish livers. Tukey and Wilcoxon’s tests showed that the Cadmium concentration in fish liver and water has significantly changed between treatments and Control; Constant Cadmium regime (CC), Variable Cadmium regime (CV), and Control (Control) and overtime. No significant changes were observed in Zinc and Cooper between treatments.

Additional file 2: Table S2.

Statistical summary of alpha-diversity changes over time and treatments. Over time, the richness and evenness have significantly changed in treatments and Control for all types of communities. The significant changes of alpha-diversity between treatments and Control were statistically demonstrated (Table 2b) using rank statistics tests (Kruskal-Wallis/Wilcoxon). The same statistics were used to compare alpha-diversity overtime. Overall, the evenness in microbial skin communities has significantly diverged between treatments and Control at T1, and richness in Gut microbial communities at T3.

Additional file 3: Table S3.

Phylogenetic divergence in host and water microbiomes. The phylogenetic distances between OTUs were computed using Gunifrac distance (see Materials and Methods). The divergence between treatments and control was assessed using PERMANOVA, and the homogeneity for group dispersions (distance from centroid) was evaluated using two multivariate tests, BETADISPER and Multi-Response Permutation Procedure (MRRP) of within versus among group dissimilarities. The significance of divergence between groups was measured by applying multiple correction tests with Benjamini-Hochberg BH (p-value< 0.05).

Additional file 4: Table S4.

Summary of correlation network features in host and water microbial communities. The networks of gut microbial communities are composed of 30, 46, 65 of CC (connected components or hubs), respectively in Control, CV, and CC groups (see sheet 1). Therefore, the nodes number (NN) was lower in the Control group (ControlNN: 410) compared to Cd regimes networks (CVNN:488; CCNN:526). On the contrary, an average of neighbours (AN) was higher in the Control network (Control AN: 4) than in Cd groups networks (CVAN:2.98; CCAN:2.82). Skin microbial networks were composed of 47, 10, 3 of CC (connected components), respectively in Control, CV, and CC groups (see sheet 1). The number of CC was inversely proportional to nodes number, and low in Control group network (ControlNN: 651) in comparison with Cd networks (CV NN: 661; CC NN:759). Similarly, the average of neighbours (AN) was lower in the Control network (Control AN: 4) compared to Cd group networks (CVAN: 10.7; CCAN: 24.6). This high connectivity in CV and CC at T1 was proportional to the high percentage of strong negative correlations (R neg. Corr. < − 0.6; B-H p-value < 0.05), observed in those groups (CV neg. Corr.: 6.3%; CC neg. Corr.: 6.9%) compared to Control group (Control neg. Corr.: 2.88%). Water microbial networks indicate variables features over time and between regimes (see sheet 2).

Additional file 5: Table S5.

Significant taxonomic changes in water and host-microbial communities over time and between treatments. In the gut, the significant changes in the relative abundance between T0 and T1 were detected for Synergistetes in CC and Tenericutes in CC and CV. Later, between T1 and T3, the relative abundance of four phyla, Euryarchaeota, Firmicutes, Proteobacteria, and Tenericutes has significantly changed in all groups. In the skin, between T0 and T1, the significant changes were observed for Firmicutes and Fibrobacteres in the CV, and Bacteroidetes and Fusobacteria in CC and CV. However, between T1 and T3, the significant changes of relative abundance concerned Actinobacteria in CC and CV; Proteobacteria, Tenericutes, Bacteroidetes, and Firmicutes in all groups. In water, between T0 and T1, the significant changes were observed for Bacteroidetes in Control; Proteobacteria in the CV; Firmicutes and Tenericutes in all groups. However, between T1 and T3, the significant changes concerned Proteobacteria and Tenericutes in CV, and Firmicutes in all groups. This table only summarized significant taxa variation with multiple test corrections (Benjamini-Hochberg p-value BH < 0.05) performed on significant Wilcoxon rank-sum test.

Additional file 6: Table S6.

This table summarises the variation of the goodness of fit (R2) to the neutral model, the percentage of neutral and non-neutral OTUs through different abundance cut-offs in the different treatment groups of host and water microbial communities at times T1 and T3.

Additional file 7: Figure S1.

Box plots of alpha-diversity variations over time and between treatments in the host and water microbial communities. The boxplots of richness and evenness variations showed different trends between treatments and Control. In the gut, the alpha-diversity showed the same tendency in all groups, except at time T3. In the skin, the evenness at T1 was higher in Cadmium treatments compared to Control while the opposite produced for richness at T3. In water, the evenness and richness were intermediate in the Control group compared to variable and constant Cadmium selection treatments, except for the evenness which was the highest in the Control group at T1. Constant Cadmium regime (CC) is in orange, variable Cadmium regime (CV) is in Yellow, and Control (Ctrl) is in green.

Additional file 8: Figure S2.

Heatmaps of significant taxonomic variation at the phylum level. This figure includes three heatmaps representing significant overtime changes of taxonomic composition at the phylum level in the gut (2a.), skin (2b.) and water (2c.). 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 heatmap () function in R were used to produce these heatmaps.

Additional file 9: Figure S3.

Heatmaps of significant taxonomic variation at the genus level. This figure indicates with three heatmaps the significant overtime changes of taxonomic composition at the genus Level in GMC (Gut Microbial Community) (2a.), SMC (Skin Microbial Community) (2b.) and WMC (Water Microbial Community) (2c.). 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 heatmap () function in R were used to produce these heatmaps.

Additional file 10: Figure S4.

Phylogenetic divergence at the community-level in the water microbiome. NMDS (non-metric Multi-Dimensional Scaling) plot of generalized Unifrac distances showing the distribution of the water samples based on the phylogenetic content of their microbiota. The plot shows a significant separation of sample groups according to treatment (see in Supplementary Table 5-b, p-values of the PERMANOVA test indicating the significance of group separations) at each time point, but treatment samples are closer to each other than the Control group at T1 and T3.

Additional file 11: Figure S5.

Boxplots of Bacteroidetes variation over time and between treatments. At T0, the relative abundance of Bacteroidetes was significantly lower in the skin compared to the water and the gut microbial communities. However, at T1, Bacteroidetes abundance was significantly higher (Wilcoxon ‘s P-value < 0.05) in the skin microbial communities only in treatment groups (CC and CV), while in the Ctrl group they showed any significant variation.

Additional file 12: Figure S6

. Water and host-microbial interactions network over time and between communities. This figure summarises the dynamic of interactions of water with the host microbiome networks. Each node size in the network is proportional to the average of the OTUs relative abundance in all samples. These networks are based on significant Spearman coefficients and were constructed using R scripts and Cytoscape software.

Additional file 13: Figure S7

. Statistical analysis of the closeness centrality in the water and host-microbiome networks. The comparison of centrality assessed with nodes closeness centrality and represented with violin plots indicate that the average of connections is significantly higher in the skin compared to the water and gut microbiomes in all treatments at T0 before disturbance. However, at T1 the centrality of node converged (which means not significantly different) between the water and skin microbiome for cadmium-treated groups. At T3, the centrality of node converged between the gut and skin microbiome in the control group. The average of nodes’ degree computed with Network Analyzer was compared using the Kruskal-Wallis test followed by Benjamini-Hochberg test. The value of 0.05 is the threshold of B-H p-value significance.

Additional file 14: Figure S8

. Skin microbiome networks built with SPIEC method. The SPIEC-EASI (SParse Inverse Covariance Estimation for Ecological Association Inference) method [55] implemented in R was applied using Meinshausen-Buhlmann’s neighbourhood selection (MB) method to estimate the inverse covariance matrix. The OTUs having low frequency occurrence (occurrence <=3) were dummied in one synthetic OTUs (black node). Red and grey edges represent negative and positive regression coefficients of the inverse covariance matrix. The node size in the network is proportional to the average of an OTU relative abundance in all samples. These networks were visualized using Cytoscape software.

Additional file 15: Figure S9.

Distribution of neutrality versus abundance cut-off and goodness of fit. This figure shows the variation of neutral OTUs percentage (Y-axis) and goodness of fit predicted by NLS models using 12 cut-offs thresholds of relative abundance percentages (gg facet panels) in the entire metacommunity at T1 and T3.

Additional file 16: Supplementary file 1

.

Additional file 17: Supplementary file 2

.

Additional file 18: Supplementary file 3

.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Cheaib, B., Seghouani, H., Llewellyn, M. et al. The yellow perch (Perca flavescens) microbiome revealed resistance to colonisation mostly associated with neutralism driven by rare taxa under cadmium disturbance. anim microbiome 3, 3 (2021). https://doi.org/10.1186/s42523-020-00063-3

Download citation

Animal Microbiome

ISSN: 2524-4671