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

doi:10.21203/rs.3.rs-52321/v1

Download PDF

Research Article

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

https://doi.org/10.21203/rs.3.rs-52321/v1

This work is licensed under a CC BY 4.0 License

Journal Publication

published 05 Jan, 2021

Read the published version in Animal Microbiome →

You are reading this older preprint version

Read the latest preprint version →

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 (CdCl₂). Yellow Perch (Perca flavescens) was used as a piscivorous ecotoxicological model. Using 16S rDNA gene based metataxonomics, quantitative 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) 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 stochasticity in the cadmium-free host and water communities, however, NST models suggested higher stochasticity ratios under disturbance.

Conclusions: Neutral models agree that water and host-microbiota assembly promoted by the occurrence frequence of 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 preponderance of negative correlations in the co-abundace networks. Our findings enhance our understanding of community assembly host-associated and free-living under anthropogenic selective pressure.

General Microbiology

CdCl2

skin and gut

potential involvement

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 (macroscopic organisms). 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 emigration, death and birth [10–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–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, 22, 23, 16] and deterministic [17, 24–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 rare taxa might play an important role in the metacommunity assembly [35–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–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.

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 two 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 (Schema 1). 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 (CdCl₂) 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 five days to reach a maximal concentration at the end of the first month (T1). This maximal concentration was maintained two months until the end of treatment (third month, T3). For the regime of Cd variable concentration (CV), the Cd concentration was gradually increased every five days to reach the target concentration (9 ppb) at the end of treatment (third month, T3). The maximal CdCl₂ concentration was set at 9 mg/L as it is the highest CdCl₂ 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 x 3 regimes x 8 tanks x 3 replicates) skin mucus swabs and 216 (3 times x 3 regimes x 8 tanks x 3 replicates) gut tract samples (Supplementary file 2). Water samples were stored in sterile bottles (Nalgene), 2 litres per tank were filtered using a polycarbonate membrane of 0.22 μm. In total, 144 filters (3 times x 3 regimes x 8 tanks x 2 replicates) were conserved in 2mL 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 CdCl₂ 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 CdCl₂ exposure regimes. For fish, liver samples after lyophilisation was digested with purified nitric acid and kept at room temperature for five days. The liver acid digestion protocol was adapted from Pierron et al. (2009). 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 were analysed using Tukey’s test and Wilcoxon rank test depending on the data (metal concentration), after assessing data distribution normality with the Shapiro test.

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 (Supp. 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 (Werner et al. 2012). 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) (Supplementary file 2).

Analysis of 16S rDNA amplicons

Sequence analysis was performed with our bioinformatic pipeline as described previously [27, 50]. The analysis of the alpha and beta diversity of metacommunities was proceeded using the Rhea package [51]. 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 [52], which considers both dominant and rare OTUs. P-values of pairwise comparisons in alpha and beta-diversity were validated with multiple correction tests [53] using the B-H (Benjamini-Hochberg) for avoiding the Type I errors (false positives). Additional details on amplicons analysis are documented in Supplementary file 2.

Correlational networks analyses

The Spearman coefficients supported with multiple corrections tests 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 (Weiss et al. 2016). 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 (Shannon et al. 2003) 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 (Weiss et al. 2016) are comparable, possible bias of indirect associations of nodes can be generated by spearman coefficient [54]. To deal with this bias, the SPIEC-EASI (SParse Inverse Covariance Estimation for Ecological Association Inference) method [54] 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 [55]. 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 [56]

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. 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 relative contributions of both deterministic and stochastic processes in community assembly, the new proposed index, normalized stochasticity ratio (NST) was computed to estimate the ecological stochasticity based on beta diversity [57]. The NST is an index to estimate average stochasticity within a group of samples. 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.

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 (Supp.Table 1). For alpha-diversity, statistical comparisons using both diversity indices for evenness (Shannon effective) and richness (Chao1) demonstrate the importance of time as a driver in microbial community richness, evenness (Figure 1-b, Supp. figure 1, Supp. Table 2), and taxonomic composition (Supp. figure 2&3) rather than treatment. The treatment reveals a significant effect only on the evenness in the skin mucosal communities at T1, and on the richness in the gut communities at T3 (Figure 1b-c & Figure 2b-c). 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 (Supp.Table 3, Figure 1c-d); however, by T3, both cadmium exposure regimes indicate significant changes in both skin and gut microbial communities compared to the control (Supp.Table 3, Figure 2c-d). 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 (Supp Table 3; Supp. figure 4).

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 [58], suggesting low overall connectivity (Figure 2a). 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 (Figure 1a).

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 (Figure 2a). 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 1). In term of nodes connectivity, the average degree was higher in the control group compared to Cd treatment groups at T1 and T3 (Figure 3).

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) and a significantly high average degree in the Cd treatment groups (Figure 3). However, at time T3, like in the gut microbiome, the skin microbiome networks reveal 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 Figure 1a) 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) (Supp. 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 (Supp. Table 5), not for control, showing a significantly higher abundance in skin compared to water and gut microbial communities for Cd groups (Supp. Figure 5). Similar to 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 (Figure 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 (Figure 5). Overall, the structure of the water networks encompassed small disconnected worlds of independent hubs.

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 module 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 (Figure 6).

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 6 ).

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 (CC_T1:T3= 93:85% ; CV_T1:T3= 94:88% ; Ctrl_T1:T3= 99.6:99.4% ), the skin (CC_T1:T3=83:91% ; CV T1:T3=86:88% ; Ctrl_T1:T3=99.4:99.3%) and the water (CC_T1:T3=81:99% ; CV_T1:T3=75:88% ; Ctrl _T1:T3=97:93%) was confidently predicted by the NLS model (Supp Table 6, Figure 1e &&2e). According to this model, the percentage of neutral OTUs and the goodness of fit (R²) are always higher in the control group compared to both Cd-treated regimes in water and host communities (Supp 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 (Supp 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 indicate that stochasticity was higher in the water microbiome compared to host microbial communities in the Control (Gut_T1=0.172 ; Skin_T1=0.200 ; Water_T1=0.326), CV (Gut_T1=0.118 ; Skin_T1=0.204 ; Water_T1=0.423) and CC (Gut_T1=0.129 ; Skin_T1=0.242 ; Water_T1=0.335). Similarly, at T3, the emigration rate indicates higher stochasticity in water compared to the host-microbial communities for the Control (Gut_T3=0.269; Skin_T3=0.177 ; Water_T3=0.488), and CV (Gut_T3=0.351 ; Skin_T3=0.138 ; Water_T3=0.424) groups, but not for the CC (Gut_T3=0.564 ; Skin_T3=0.169 ; Water_T3=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, Figure 1e && 2e, Supp. Figure 7 , Supp Table 6).

To investigate the magnitude of stochastic processes on 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 assembly 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 (NST_avg.) for almost all communities is higher than 0.5 (Figure 7). Comparing stochasticity ratios between water and host microbiota, the NST modelling lead to same results of NLS models. The NST_avg. 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, the NST and NLS models agree in terms of the dominance of neutral processes, according to the stochasticity ratio (NST_avg>0.5), and goodness of fit (R²>0.5). However, in NLS models we always observe that neutral OTUs are higher in the control group, whereas NST models indicate significantly higher stochasticity in Cd treatment groups early at T1 in the skin microbiome and lately at T3 in the gut microbiome.

Our study evidenced salient differential changes in community assembly across three community types (environmental or host-associated), time and xenobiotic exposure regimes. This highlights for the first time 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 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. 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. Lastly, to predict the nature of evolutionary processes driving the metacommunity changes, application of both non-linear least squares model (NLS), and normalized stochastic ratio model (NST) to the data (R²> 0.5; NST_avg> 0.5) agree that neutral processes drove taxonomic drift in host and water microbial communities in the control and Cd-treated groups, suggesting high stochasticity in the community assembly under disturbance [57]. However, NLS and NST disagree on the extent of stochastic processes on the microbiota assembly in the skin at T1 and the gut at T3. This conflicting results would result from the differences in terms of null hypotheses underlying NLS and NST models. The NST model was developed and applied to estimate stochasticity in the succession of a groundwater microbial community in response to organic carbon (vegetable oil) injection [57]. The authors of NST have also compared their model, with not normalized stochasticity strength measure (an early version of NST) [59] and the neutral OTUS predicted by the NLS model of Sloan et al. [16]. They hastily concluded that the accuracy and precision of NST are better than NLS for the majority of similarity metrics examined. They often demonstrated higher stochastic values with NST, especially under disturbance conditions. So expectedly, our results (in the skin at T1 and the gut at T3 figure 7) with NST also demonstrate higher stochasticity under disturbance using the same null-model (proportional fixed) and abundance-based metrics (Figure 7). The NLS model has different assumptions based on the estimation of the emigration rates from a community source but is not based on similarity distance calculation. With the NST model, we could not distinguish between neutral and non-neutral OTUs, although the authors missed considering the emigration rate estimated in NLS in their comparison, which is more accurate in assessing the demographic stochasticity than the percentage of neutral OTUS fitting the Sloan’s model within the 95% confidence interval. Furthermore, in contrast to NLS model, the NST null-modelling is not based on an empirical community source, but on a family of null-models assuming that occurrence frequency is i “Equiprobable”, which means that all taxa have equal probability to occur; ii “Proportional”, which means that the occurrence probability of a taxon is proportional to its observed occurrence frequency; iii “Fixed”, which means that the occurrence frequency of a taxon is fixed as observed. From our point of view the comparison of NST with NLS is not informative on 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, and it is difficult to distinguish ecological stochasticity from the noise caused by deterministic environmental factors, as shown in their simulation [57]. Despite these limitations, NST agree with NLS model on the dominance of demographic stochasticity, a conclusion based on arbitrary threshold NST_avg> 0.5 for the former and NLS goodness of fit > 0.5 for the latter.

Furthermore, we would stress that NLS model, with our null hypothesis based on a real bacterial community (i.e. T0, 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: 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 in the control group would approach specialized, partitioned microbiome networks in the gut and skin, reflecting niche differentiation and microbiota maturation along with perch juvenile’s development. 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 [60] and lesser interactions and more stochastic changes under disturbance.

Concerning 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” [61]; 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 [61], the rise of negative correlations suggests a dysbiosis state of skin-associated microbiota [62]. 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 [63, 64]. 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 shift 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, 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) [65] and the long jaw mudsucker (Gillichythys mirabilis) [66, 67], 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 [62]. 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).

Low abundant or rare OTUs have been demonstrated to play a pivotal role in community assembly [35–39] either in promoting homeostasis [37] or dysbiosis [68]. Therefore, we applied non-linear least squares model (NLS) 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 (Supp. Table 6) suggested a neutral role from 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].

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. In view of the modelling frameworks, neutral processes were thus the major forces of the community assembly in the environmental microbiome, 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 colonisation 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 colonisation 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 [69–71].

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.

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;

PCR: Polymerase chain reaction;

PERMANOVA: Permutational analysis of variance;

RDP: Ribosomal Database Project;

R²: goodness of fit

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

Availability of data and materials

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

Competing interests

The authors declare that they have no competing interests.

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 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 wet-lab steps including libraries preparations. P-L.M. quantified DNA of libraries, B.C. produced and analysed the results. B.C. wrote the manuscript. M.L and N. D. corrected the and revised all the version of this manuscript. All authors reviewed the manuscript.

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.

Johnson KV-A, Foster KR. Why does the microbiome affect behaviour? Nat Rev Microbiol 2018; 16: 647–655.
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–654.
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.
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-18.
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.
Faust K, Raes J. Microbial interactions: from networks to models. Nat Rev Microbiol 2012; 10: 538–550.
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–356.
Vellend M. Conceptual synthesis in community ecology. Q Rev Biol 2010; 85: 183–206.
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–789.
Bell G. Neutral Macroecology. Science 2001; 293: 2413–2418.
Hubbell Stephen P. Neutral theory in community ecology and the hypothesis of functional equivalence. Funct Ecol 2005; 19: 166–172.
Hubbell SP. Neutral theory and the evolution of ecological equivalence. Ecology 2006; 87: 1387–1398.
Chase JM, Leibold MA. Ecological niches: linking classical and contemporary approaches. 2003. University of Chicago Press, Chicago.
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–529.
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.
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–740.
Stegen JC, Lin X, Konopka AE, Fredrickson JK. Stochastic and deterministic assembly processes in subsurface microbial communities. ISME J 2012; 6: 1653–1664.
Zeng Q, Sukumaran J, Wu S, Rodrigo A. Neutral Models of Microbiome Evolution. PLOS Comput Biol 2015; 11: e1004365.
Zeng Q, Wu S, Sukumaran J, Rodrigo A. Models of microbiome evolution incorporating host and microbial selection. Microbiome 2017; 5: 127.
Foster JA, Krone SM, Forney LJ. Application of Ecological Network Theory to the Human Microbiome. Interdisciplinary Perspectives on Infectious Diseases. https://www.hindawi.com/journals/ipid/2008/839501/. Accessed 16 Sep 2018.
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–9612.
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.
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–15350.
Morrison-Whittle P, Goddard MR. Quantifying the relative roles of selective and neutral processes in defining eukaryotic microbial communities. ISME J 2015; 9: 2003–2011.
O’Dwyer JP, Kembel SW, Sharpton TJ. Backbones of evolutionary history test biodiversity theory for microbes. Proc Natl Acad Sci 2015; 112: 8356–8361.
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.
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.
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–9698.
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–1094.
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.
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–2125.
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–1262.
Schmidt VT, Smith KF, Melvin DW, Amaral‐Zettler LA. Community assembly of a euryhaline fish microbiome during salinity acclimation. Mol Ecol 2015; 24: 2537–2550.
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.
Banerjee S, Schlaeppi K, Heijden MGA van der. Keystone taxa as drivers of microbiome structure and functioning. Nat Rev Microbiol 2018; 16: 567–576.
Fierer N. Embracing the unknown: disentangling the complexities of the soil microbiome. Nat Rev Microbiol 2017; 15: 579–590.
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–862.
Lynch MDJ, Neufeld JD. Ecology and exploration of the rare biosphere. Nat Rev Microbiol 2015; 13: 217–229.
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–1602.
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.
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–1064.
Couture P, Rajender Kumar P. Impairment of metabolic capacities in copper and cadmium contaminated wild yellow perch (Perca flavescens). Aquat Toxicol 2003; 64: 107–120.
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–245.
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 Organ 2005; 63: 77–84.
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–258.
Ryman JE, Walleghem JLAV, Blanchfield PJ. Methylmercury levels in a parasite (<Emphasis Type="Italic">Apophallus brevis</Emphasis> metacercariae) and its host, yellow perch (<Emphasis Type="Italic">Perca flavescens</Emphasis>). Aquat Ecol 2008; 42: 495–501.
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–364.
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.
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.
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.
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.
Chen J, Bittinger K, Charlson ES, Hoffmann C, Lewis J, Wu GD, et al. Associating microbiome composition with environmental covariates using generalized UniFrac distances. Bioinforma Oxf Engl 2012; 28: 2106–2113.
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.
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.
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.
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–1936.
Ning D, Deng Y, Tiedje JM, Zhou J. A general framework for quantitatively assessing ecological stochasticity. Proc Natl Acad Sci 2019; 116: 16892–16898.
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–391.
Zhou J, Deng Y, Zhang P, Xue K, Liang Y, Van Nostrand JD, et al. Stochasticity, succession, and environmental perturbations in a fluidic ecosystem. Proc Natl Acad Sci 2014; 111: E836–E845.
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.
Buffie CG, Pamer EG. Microbiota-mediated colonization resistance against intestinal pathogens. Nat Rev Immunol 2013; 13: 790–801.
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.
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.
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–185.
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–1284.
Bano N, DeRae Smith A, 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–2641.
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–223.
Hajishengallis G, Darveau RP, Curtis MA. The keystone-pathogen hypothesis. Nat Rev Microbiol 2012; 10: 717–725.
Snedeker SM, Hay AG. Do interactions between gut ecology and environmental chemicals contribute to obesity and diabetes? Environ Health Perspect 2012; 120: 332–339.
Tilg H, Kaser A. Gut microbiome, obesity, and metabolic dysfunction. J Clin Invest 2011; 121: 2126–2132.
Tilg H, Moschen AR. Microbiota and diabetes: an evolving relationship. Gut 2014; 63: 1513–1521.

Table 1. Host microbiome networks modularity.

Skin microbiome network	T0	T1	T3
Ctrl	78	47	42
CV	51	104	28
CC	56	149	30

Gut microbiome network	T0	T1	T3
Ctrl	375	12	94
CV	315	41	50
CC	445	44	65

SuppFigure1.pdf
Box plots representation of alpha-diversity variations over time and between treatments in 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 time3. 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 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 (Control) is in green.
SuppFigure2.pdf
representing significant overtime changes of taxonomic composition at the phylum level in the gut (2a.), in the 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.
SuppFigure3.pdf
Heatmaps of significant taxonomic variation at the genus level. This figure represented with three heatmaps the significant overtime changes of taxonomic composition at the genus Level in GMC (Gut Microbial Community) (2a.), in 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.
Suppfigure4.pdf
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 host-microbial samples based on the phylogenetic content of their microbiota. The plot shows a significant separation of sample groups according to treatment (see in Supp. 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.
SuppFigure5.pdf
Boxplots of Bacteroidetes variation over time and between treatments. At T0, the relative abundance of Bacteroidetes was significantly lower in Skin compared to the water and the gut microbial communities. However, at T1, Bacteroidetes abundance was significantly higher (Wilcox ‘s P-value <0.05) in Skin communities only in treatment groups (CC and CV), while in the Ctrl group they showed any significant variation.
SuppFigure6.pdf
Skin microbiome networks built with SPIEC method.
The SPIEC-EASI (SParse Inverse Covariance Estimation for Ecological Association Inference) method [54] implemented in R was applied using meinshausen-buhlmann's neighborhood 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.
SuppFigure7.pdf
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.
SupplementaryTable1.xlsx
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 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 between Zinc and Cooper between treatments.
SupplementaryTable2.xlsx
Statistical summary of alpha-diversity changes over time and treatments. Over time, the richness and evenness have significantly changed in treatments and Control, in all types of communities. The significant changes of alpha-diversity indexes difference between treatments and Control were statistically compared (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.
SupplementaryTable3.xlsx
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).
SupplementaryTable4.xlsx
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, 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 connectivity higher in CV and CC at T1 was proportional to the higher 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 displayed variables features over time and between regimes (see sheet 2).
SupplementaryTable5.xlsx
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, 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.
SupplementaryTable6.xlsx
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.

Download PDF

Journal Publication

published 05 Jan, 2021

Read the published version in Animal Microbiome →

Editorial decision: Minor revision
14 Sep, 2020
Review #2 received at journal
13 Sep, 2020
Review #1 received at journal
27 Aug, 2020
Reviewer #2 agreed at journal
17 Aug, 2020
Reviewer #1 agreed at journal
07 Aug, 2020
Reviewers invited by journal
06 Aug, 2020
Submission checks completed at journal
03 Aug, 2020
Editor assigned by journal
03 Aug, 2020
Editor invited by journal
02 Aug, 2020
First submitted to journal
31 Jul, 2020

You are reading this older preprint version

Read the latest preprint version →

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

Status:

Journal Publication

Version 1

Abstract

Figures

Background

Materials And Methods

Results

Discussion

Conclusions

Abbreviations

Declarations

References

Tables

Supplementary Files

Status:

Journal Publication

Version 1