Exposing oyster larvae to Microorganism-Enriched seawater shifts their bacterial microbiota throughout their lifespan and in the next generation
To investigate whether an environmental microbial exposure early in life could influence the trajectory of the oyster microbiome, we developed an experimental setup to compare the effects of control or Microbially-Enriched seawater environments during early larval development. Pathogen-free larvae (F1 generation) were produced in a bio-secured (filtered and UV-treated seawater) hatchery. A subset of F1 larvae were exposed from two hours to ten days post-fertilization to a microorganism-enriched environment by cohabitation with oysters transferred from a natural environment during a POMS-free period (Microbial Enriched condition, ME seawater, Fig. 1). As a control, a subset of F1 larvae were raised in bio-secured conditions with no cohabitation or no microbial exposure (Control seawater Fig. 1). From 10 days onward, both ME-exposed and control oysters were raised in the same bio-secured conditions. A part of these two oyster subsets were maintained in bio-secured conditions and reproduced, one year later, to generate the F2 generation (Fig. 1). Between the ME and control oysters, we observed equivalent developmental success and survival rate in the F1 generation (Table 1 of the additional file 3). Moreover, the absence of OsHV-1 was confirmed for ME and control oysters during exposure time. The nature of the seawater treatments was evaluated by analyzing the bacterial load and composition of the ME and control seawater by qPCR targeting 16S rRNA genes and 16S barcoding at day 2 post-fertilization. As expected, the ME contained 5-fold more total bacteria than the control seawater (unpaired t-test with Welch’s correction, p = 0.001) (Fig. 2a) and carried a more diverse microbiota as evidenced by the Chao1 index (ANOVA, p < 0.05) (Fig. 2b). This trend was confirmed by plating the seawater sampled in each tank containing the recipient larvae on marine agar, revealing that ME tanks contained 3.4 times more cultivable bacteria (16544 CFU/ml) than the control seawater tanks (4916 CFU/ml) (Wilcoxon test: p < 0.05) (Fig. 2c). Altogether, the ME condition was considered as an exposure to a safe, Microorganism-Enriched environment.
To test the immediate and long-term impact of early ME exposure on the oyster microbiota, we analyzed the bacterial community composition by 16S amplicon sequencing in both F1 and F2 whole body oysters (Additional file 4). Differences in composition and diversity were evidenced between ME-exposed and control oyster larvae during the ME seawater exposure (Fig. 3a, Fig. 1 of the additional file 3). Among the 41 genera that had higher relative abundance in the ME-exposed larvae compared to control larvae, 29 (70%) were also more highly represented in the ME seawater (Fig. 4, Additional file 5). Conversely, among the 33 genera that make up a higher proportion in the control compared to ME-exposed oysters, 18 (54.5%) were also more represented in the control seawater, strongly suggesting that the microorganisms from the ME seawater colonized the oyster larvae during the exposure.
Dissimilarity analysis, based on the Bray Curtis index, showed that the oyster microbiota profiles clustered first by developmental stage (Fig. 2 of the additional file 3) and then by treatment (Fig. 3a to d). This analysis indicated that the microbiota composition of ME-exposed vs control oysters was significantly different, not only during the exposure (Fig. 3a, permutation test p-val = 0.002), but also several months later during the F1 juvenile stage (Fig. 3b permutation test p-val = 0.005) as well as during larval and juvenile stages of the F2 generation (Fig. 3c and d, permutation test p-val = 0.001 and 0.005 respectively).
Taken together, this barcoding analysis clearly indicated that the oyster microbiota significantly shifts across developmental stages, but despite this strong developmental effect, the ME seawater exposure during larval stages induced a persistent modification of the oysters’ bacterial microbiota composition that even persisted in the subsequent generation.
Early life microbial exposure primes intergenerational immunity against Pacific Oyster Mortality Syndrome
To test whether ME exposure of oyster larvae can produce a long-term impact on their resistance to disease, we conducted an ecologically realistic experimental infection mimicking the Pacific Oyster Mortality Syndrome (POMS) disease on juvenile oysters from F1 and F2 generations (Fig. 1 of the additionnal file 1)[28, 33]. The number of surviving oysters was monitored for 300 hours while oyster OsHV-1 load was measured before the onset of the mortalities (Fig. 5). The increase in virus load during the first 48 hours confirmed successful infection. The viral load was significantly lower in the ME-exposed oysters or their offspring compared to the control lineage (p-val of two-way ANOVA with Bonferroni’s correction for multiple comparisons test: p-val < 0.01) (Fig. 5a and b). Consistent with these results, we observed that ME-exposed oysters had a better survival rate compared to controls in both F1 (66.3% vs. 57.4%, Log rank test, p-val < 0.05) and F2 generations (31,4% vs. 18.4%, Log rank test, p-val < 0.0001) (Fig. 5c and d, respectively). These results were confirmed in a parallel field infection test conducted with oysters from both F1 and F2 generations (F1: 16.4% vs. 14%, log-rank test, p-val < 0.001, F2: 8.5% vs. 1%, log-rank test, p-val < 0.0001) (Fig. 3 of the additional file 3).
No evidence for genetic selection as the mechanism of increased immune capacity
We investigated if a genetic selection could have occurred through ME exposure and would have selected more resistant oysters based on specific allele associations. To this end, we evaluated genome-wide SNP allele frequencies in juvenile oyster samples using whole genome sequencing (WGS). Principal components analysis (PCA) of these data showed little genetic divergence between the ME-exposed vs. control oysters for the F1 and F2 generations (Fig. 6a). Next, we conducted a genome scan comparing allele frequencies in ME-exposed vs. control oysters using the FLK test to interrogate any signals of positive selection. The FLK statistic considers genome wide allele frequency data in a set of populations and aims at detecting positions where genetic differentiation between these populations is higher than expected under neutral evolution. It returns for each SNP a p-value allowing to reject or accept neutrality. In the case of genetic selection on some SNPs, an excess of low p-values is expected. No such excess was detected here, revealing an absence of genetic selection between exposed and control lines for the F1 and F2 generations (Fig. 6b). Furthermore, no significant SNPs could be detected based on a FDR value below 0.05 (even below 0.15 for F2, Table 2 of the additional file 3). This absence of genetic selection is consistent with the fact that the survival rate of ME-exposed larvae was not significantly lower than of control larvae (Table 1 of the additional file 3). Altogether, these findings indicate that genetic alterations are not responsible for the increased resistance among the ME-exposed oyster lineage.
Upregulation of immune-related and other transcripts in microbially-exposed oyster lineages
Next, we asked whether ME exposure impacted oyster gene expression by performing transcriptomic analyses on larvae and on juveniles just before and during the POMS disease breakouts for both F1 and F2 oysters. During ME exposure in the F1 larval stages, we observed a large shift in gene expression (3410 and 1100 DEGs at day 2 and 10) (Table 3 of the additional file 3 and Additional file 6). However, the difference in gene expression between ME-exposed and control oysters is much more nuanced at juvenile stages (35 DEGs at day 120). These observations were similar to what we observed in the F2 generation (6 029 DEGs at day 10 and 120 DEGs at day 120).
To investigate which biological processes are modulated by ME exposure, we performed a rank-based gene ontology analysis (RBGOA; false-discovery rate [FDR] < 0.01) (Additional file 7). The broad gene expression shifts in larval oysters encompassed many functional annotations, including general cellular process, metabolism, response to environmental stimulus, infection and immune response, transcription and gene expression, development, cell fate, RNA process, translation and protein processing, signal transduction and transport. In juvenile oysters of both generations, upregulation of genes involved in responses to external stimuli and immunity persisted from the larval stage, suggesting a potential role for these genes in mediating resistance to POMS at the time of infection. During POMS disease onset at the juvenile life stage, we observed a strong over-representation of immune functions, especially in the F2 generation (Fig. 7b). Analysis of the individual genes driving this enrichment revealed gene families typically involved in microbial-associated molecular pattern (MAMP), recognition (PGRP, lectins, scavenger receptors, TLR, RLR, Macrophage receptor), innate immune pathways (components of IFN- TLR- JAK/STAT pathways as MyD88, IRF2, STING), interaction with bacteria (dual oxidase), and antimicrobial effectors (TNF, proteinases, SOD, Interferon- stimulated genes) (Additional file 6). These immunity-linked families were found differentially expressed in both generations, especially at larval stages, although the individual genes encoding for these immune functions were generally different in F1 compared to F2 generation (different CGI numbers). In addition, a closer look at antimicrobial peptides or proteins (AMP) expression revealed a significant over-expression in ME-exposed compared to control oysters, either during the exposure period at larval stages in F1 (Big-Def1 and BPI at day 2 and day 10) (Fig. 7c), or in F2 (Big-Def2 and DefH at day 10) (Fig. 7d).
Apart from immune functions, our transcriptomic analysis highlighted that ME exposure during larval stages also affected key metabolic pathways. The expression of genes encoding for enzymes involved in glycolysis and the TCA cycle was lower in both generations during the larval stages, whereas the oxidative phosphorylation pathway and folate metabolism enzymes were downregulated at day 10 of the F2 generation only (Additional file 8). We also observed that functions linked to chromatin structure (RBGOA analysis, Additional file 7) and genes encoding for DNA methylation machinery enzymes were repressed at day 10 of the F2 generation (Additional file 8).
Taken together, these transcriptomic analyses showed that the ME seawater larval exposure of C. gigas resulted in modification of the immune response of the oysters. This immunomodulatory effect was maintained up to the juvenile stages and in the subsequent generation. These results support the idea that transcriptional changes may be responsible for the increased immune capacity that we observed in the survival assay.
Differences in DNA methylation may explain the transcriptional changes observed in microbially-exposed oyster lineages
The observed multi-generational impact of the ME exposure on oyster survival capacities and transcriptomic response, as well as the absence of genetic selection between conditions, led us to investigate the impact on epigenetic information through analysis of differentially methylated regions between ME and control oyster lineages. Whole Genome Bisulfite sequencing (WGBS) analysis was performed on oysters sampled at day 10 and 120 of the F1 and F2 generations (Table 2 of the additionnal file 1). We found that C. gigas DNA was mainly methylated in a CpG context with a mosaic-type cytosine methylation pattern as previously described [57]. The PCA results of the global pattern of cytosine methylation data showed a clustering according to developmental stages but not according to the treatment (ME-exposed or control oysters) (Fig. 4 of the additional file 3). This result suggested that the cytosine methylation pattern changed during development as previously observed by others [58] and this observation was confirmed by a global decrease in the cytosine methylation level observed from larval to juvenile stages (1.77–1.58% for F1 and 1.82–1.56% for F2, respectively, p-val from Wilcox test < 0.01 for both generations) (Fig. 5 of the additional file 3). Although ME exposure did not appear to strongly affect the level of cytosine methylation at the genome wide scale (Fig. 5 of the additional file 3), a trend toward a hyper-methylation was observed in ME-exposed compared to control larvae in the F1 generation. In contrast, an opposite trend toward a hypo-methylation was observed in F1 juveniles and in F2 larvae (Fig. 5 of the additional file 3).
To gain deeper insights into the impact of the ME exposure on methylation patterns of oysters, we used DMR-Seq software [55] to identify Differentially Methylated Regions (DMRs) between ME-exposed and control oysters for each generation. The differential methylation analysis led to the detection of 4 325 and 5 531 DMRs for larvae (day 10) of the F1 and F2 generation, respectively, and 4 985 and 5 207 for juveniles (day 120) of the F1 and F2 generation, respectively (Table 4 of the additional file 3, Additional file 9). Hyper-methylated DMRs in ME-exposed compared to control oysters were more frequent than hypomethylated DMRs at day 10 of the F1 generation (57.4% hyper-methylated vs. 42.6% hypo-methylated DMRs). However, hypo-methylated DMRs in ME-exposed oysters were more frequent at day 120 of F1 generation (40.9% hyper-methylated vs. 59.1% hypo-methylated DMRs) and at day 10 of the F2 generation (22.2% hyper-methylated vs. 77.8% hypo-methylated DMRs) (Table 4 of the additional file 3). These observations were in agreement with the previous trend observed at the genome wide level (Fig. 5 of the additional file 3) and corroborated the transcriptional analysis which depicted a strong repression of genes involved in the DNA methylation pathway in the larvae of the F2 generation (Additional file 8).
Next, we analyzed DMRs that intersected with gene positions, defining these regions as Differentially Methylated Genes (DMGs), and asked which functional annotations are overrepresented among DMGs. According to this analysis, the functions mostly impacted by DNA methylation changes were related to general cellular process, metabolism, response to environmental stimulus, signal transduction, translation and protein processing and development (Additional file 10). Similar functions were found to be modified in the oyster transcriptome in response to the ME exposure (Additional file 7 and 10). Although immune functions were not statistically highlighted by the RBGOA analysis, 128 DMGs were found in genes encoding for immune functions. Genes coding for the Interferon pathway, immune signaling pathway, viral production and ubiquitin modification displayed changes in their cytosine methylation profiles (Fig. 6 of the additional file 3). However, we did not observe a canonical association between expression levels and methylation changes when analyzing correlations between methylation and transcriptome profiles (Fig. 7 of the additional file 3).
Some DMRs were meiotically inherited from the F1 to the F2 generations (Fig. 8). 48 hyper-methylated DMRs and 120 hypo-methylated DMRs were conserved from F1 to F2 at the larval stage, and in the juvenile stage we detected 147 hyper-methylated DMRs and 252 hypo-methylated DMRs conserved from F1 to F2 (Table 1). To test whether this number of meiotically heritable DMRs was higher than would be expected by chance, we randomly generated 5 000 files containing artificial DMRs of identical size and number as the DMRs detected in the real dataset of the F1 generation and intersected these with the real dataset of the F2 DMRs. Mean values of the number of intersections between F2 DMRs and these randomized regions were always significantly lower than the number of intersections between F1 and F2 DMRs from the real dataset (Table 1 and Additional file 11), indicating that the similarity of DMRs in F1 and F2 did not occur by chance.
Taken together, our epigenetic analysis shows that microbial exposure during larval stages impacts the DNA methylation pattern in both the directly exposed oysters as well as their offspring. The DNA methylation profile of genes involved in immune functions were clearly impacted in both generations although a functional consequence on their expression was not evidenced. We showed that inheritance in the DMRs between F1 and F2 generation was not obtained by chance which suggests that DNA methylation changes can be inherited via epigenetic memory from the F1 to F2 generation.
Table 1: Shared DMRs between F1 and F2 generations are likely the result of intergenerational epigenetic inheritance
The number of genomic coordinate intersections between either (a) F1 and F2 DMRs induced after the microbial exposure (ME induced DMRs), or (b) mock randomly-generated F1 DMRs (randomly generated DMRs) and F2 DMRs (ME induced DMRs) was compared by a T-test. Hyper- and hypo-methylated DMRs were tested separately within each developmental stage. See also Additional file 11.
Day 10 hyper-methylation |
| F1 ME induced DMRs (n = 2482) | F1 randomly generated DMRs (n = 2482) |
F2 ME induced DMRs (n = 1230) | 48.0 ± 17.9 | 7.20 ± 2.7 |
Day 10 hypo-methylation |
| F1 ME induced DMRs (n = 1843) | F1 randomly generated DMRs (n = 1843) |
F2 ME induced hypo DMRs (n = 4301) | 120.0 ± 27.8 | 18.2 ± 4.2 |
Day 120 hyper-methylation |
| F1 ME induced DMRs (n = 2040) | F1 randomly generated DMRs (n = 2040) |
F2 ME induced DMRs (n = 2550) | 147.0 ± 38.5 | 14.8 ± 3.9 |
Day 120 hypo-methylation |
| F1 ME induced DMRs (n = 2945) | F1 randomly generated DMRs (n = 2945) |
F2 ME induced DMRs (n = 2657) | 252.0 ± 48.1 | 27.6 ± 5.3 |