Fluctuations in chromatin state at regulatory loci occur spontaneously under relaxed selection and are associated with epigenetically inherited variation in C . elegans gene expression

Abstract

Some epigenetic information can be transmitted between generations without changes in the underlying DNA sequence. Changes in epigenetic regulators, termed epimutations, can occur spontaneously and be propagated in populations in a manner reminiscent of DNA mutations. Small RNA-based epimutations occur in C. elegans and persist for around 3–5 generations on average. Here, we explored whether chromatin states also undergo spontaneous change and whether this could be a potential alternative mechanism for transgenerational inheritance of gene expression changes. We compared the chromatin and gene expression profiles at matched time points from three independent lineages of C. elegans propagated at minimal population size. Spontaneous changes in chromatin occurred in around 1% of regulatory regions each generation. Some were heritable epimutations and were significantly enriched for heritable changes in expression of nearby protein-coding genes. Most chromatin-based epimutations were short-lived but a subset had longer duration. Genes subject to long-lived epimutations were enriched for multiple components of xenobiotic response pathways. This points to a possible role for epimutations in adaptation to environmental stressors.

Introduction

Epigenetic gene regulation is fundamental to the establishment and maintenance of cell identity during development. Whilst many epigenetic features are erased in germ cell development and the early stages of embryonic development, some transgenerational epigenetic inheritance phenomena have been rigorously established across a wide variety of experimental settings in several different species [1–5]. Several molecular mechanisms capable of transmitting epigenetic information between generations have been demonstrated, including short (18-36nt) non-coding RNA species, histone post-translational modifications, methylation of specific bases within DNA and prions [4,6]. Despite this burgeoning understanding of transgenerational epigenetic inheritance in laboratory settings we still understand little about its significance across evolution [7–9]. Of particular interest is whether transgenerational epigenetic inheritance can contribute to evolutionary processes in populations, including genetic drift and natural selection [5,10].

Heritable changes in epigenetic states may alter gene expression patterns but do not require underlying changes in DNA sequence [11]. When such changes arise in populations they can be referred to as epimutations [12,13]. Epimutations may be induced in response to an environmental stimulus [14–17], but have also been shown to arise spontaneously under controlled conditions [18,19]. Establishing the rate and stability of epimutations requires the use of approaches from evolutionary biology designed to map DNA sequence-based mutations, known as Mutation Accumulation (MA) experiments [20–22]. In the MA approach multiple lines derived from a common ancestor are propagated for many generations at reduced population size. Under these conditions drift dominates over selection so neutral, deleterious and beneficial mutations can be identified [20,23,24]. The MA framework was used to study epimutations in plants, focussing on changes in DNA methylation [25–27]. Spontaneous changes in DNA methylation arise more rapidly than DNA sequence changes, but have limited stability, lasting around 5–10 generations on average [18,28,29].

We recently used the nematode C. elegans as a model to study epimutation rate, spectrum and stability. Transgenerational epigenetic inheritance in C. elegans can be mediated by a class of small RNAs known as 22G-RNAs, which are made by RNA dependent RNA polymerase and map antisense to protein-coding genes [30–34]. We demonstrated that 22G-RNA-based epimutations arise spontaneously in C. elegans MA lines [35]. Epimutations were short-lived, lasting for two to three generations on average, although some epimutations lasted for at least 10 generations. Similar short-term variation was shown to affect metabolite levels, suggesting that short-term epigenetic inheritance could be influential in phenotypic variation [36]. This suggests that epimutations are unlikely to contribute to long-term divergence within populations but could contribute to short-term adaptive processes. However, there was no enrichment of epimutations in 22G-RNAs to target particular categories of genes, thus the functional relevance of this process remains unclear [35].

In addition to 22G-RNAs, other forms of epigenetic regulation could contribute to epimutations in C. elegans. Chromatin is critical to the regulation of gene expression and is structurally dynamic in terms of its susceptibility to undergo remodelling in development and in response to external stimuli [37,38]. The preeminent features of chromatin organisation on C. elegans autosomes are large blocks of constitutively active genes marked by H3K36me3-containing nucleosomes and tissue-specific genes marked by H3K27me3 [37,39]. Additionally, transposable elements are marked by H3K9me2/3 nucleosomes [40]. Chromatin modifying enzymes contribute to transgenerational epigenetic inheritance phenomena either downstream or in parallel to small RNAs [41–46]. Here, we used the MA framework to investigate chromatin-based epimutations using ATAC-seq [47,48] to investigate large-scale changes in chromatin organisation. We demonstrate that heritable epimutations in chromatin structure occur, and are linked to changes in gene expression. Although most epimutations are very short-lived, a subset of epimutations is more durable. Intriguingly, long-lived epimutations are enriched for genes involved in defence response functions, suggesting that epimutations may have consequences for short-term adaptive responses for C. elegans in natural environments.

Results

1. Spontaneous chromatin epimutations occur under relaxed selection

We reasoned that in addition to 22G-RNAs, chromatin-based epimutations could also be implicated in transmission of novel gene expression patterns through the germline (Fig 1A). Chromatin-based epimutations could therefore be a source of potential phenotype variation, which depending on stability and duration could be acted on by selective forces [5,10,49].

[Figure omitted. See PDF.]

Fig 1. Identification of spontaneous chromatin-based epimutations under relaxed selection.

A. Mechanisms of transgenerational epigenetic inheritance. 22G-RNAs have been established in previous work as drivers of epigenetic inheritance through the germ line. Spontaneous changes in chromatin state (epimutations) are an alternate mechanism. B. Experimental design. Three independent lineages of C. elegans nematodes were grown with bottlenecking for a total of 20 generations including a non-restricted pre-mutation-accumulation (PMA) generation (Gen. 0). Paired RNA-seq and ATAC-seq libraries were produced from the PMA generation and at alternate generations from generations 2–20. RNA was further processed to produce small RNA libraries. C. Identification of chromatin-based epimutations. Epimutated loci were defined as regulatory element loci with ATAC-seq counts (rings on plot) that were significantly greater (UP) or lower (DOWN) than that of the same locus in the PMA generation according to a Z-score cut off. UP epimutations correspond to increased opening of chromatin while DOWN epimutations reflect chromatin closing. D. Simulation to identify optimal Z-score cut-off. X-axis shows range of Z-score cut-offs tested. Y-axis shows difference between the percentage of epimutations inherited for two or more generations compared to the predicted percentage if epimutations occurred randomly but were never inherited.

https://doi.org/10.1371/journal.pgen.1010647.g001

We propagated three independent Mutation Accumulation lines of the laboratory strain N2 C. elegans for 20 generations under identical environmental conditions. The population size was limited to two hermaphrodite animals to minimise selective forces [20] (Fig 1B). We collected embryos every two generations and profiled chromatin accessibility using ATAC-seq, small-RNAs using small RNA-seq and protein-coding genes using polyA-RNA-seq.

To identify genes with altered gene expression or antisense 22G-RNA levels we obtained normalized expression and 22G-RNA counts for each protein-coding gene and used a linear model to calculate deviations from the parental line. We then transformed the deviation into a Z-score for each gene for each generation across each lineage (S1–S6 Tables). To identify ATAC-seq changes we obtained accessibility information for regulatory elements on the basis of a previously published set [50]. We used a linear model to compare this to the parental line and transformed this into a Z-score as above (Fig 1C and S7–S9 Tables).

In order to identify putative epimutations from these data, we sought to identify a Z-score that maximised the probability of identifying true epimutations rather than false positives. According to the conventional definition of epigenetic inheritance [1], only inherited changes appearing in multiple successive generations are true epimutations. However, some epigenetic changes that appear in successive generations may have arisen by chance in the same locus over consecutive generations without transmission between generations, which we term false positive epimutations. We cannot distinguish these in any individual case. However, if true epimutations are common then we would expect the proportion of identical epimutations appearing in consecutive generations to be greater than expected given the rate at which epimutations occur overall. At a low threshold for defining an epimutation, there will be more false positives, but increasing the stringency will result in fewer true epimutations being detected. Across all Z-scores tested, the observation of chromatin state changes in two or more consecutive generations was significantly greater in the real data than in the simulated data and was at a maximum Z-score of 2.25 (Fig 1D). We set the threshold for identification of epimutations at 2.25 for subsequent analysis. Therefore, increased ATAC-seq signals (Z > 2.25) were ‘UP’ chromatin-based epimutations, indicating more open chromatin states. Decreased ATAC-seq signals (Z < -2.25) were ‘DOWN’ chromatin-based epimutations, indicating more closed chromatin states. We took the same approach to identify small RNA-based epimutations, for which ‘UP’ indicated increased small RNA levels while ‘DOWN’ indicated decreased small RNA levels.

2. Chromatin accessibility and 22G-RNA diversify at similar rates during evolution of C. elegans in the laboratory

We identified chromatin-based epimutations in three independent lineages propagated under minimal selection and compared their properties to epimutations in small non-coding RNAs and gene expression (S10–S18 Tables). We investigated microRNAs, piwi RNAs and 22G-RNAs separately. Heritable epimutations in miRNA expression levels were very rare, although interestingly slightly above what would be expected by chance. Heritable epimutations in piRNA expression levels were more common but affected fewer loci than changes in 22G-RNAs mapping to protein-coding genes (S1 Fig). On the basis of this and the known role of 22G-RNAs in transgenerational epigenetic inheritance in C. elegans [35], we focussed on 22G-RNA-based epimutations. On average, about 1.6% of protein coding genes exhibited changes in expression level per generation (Fig 2A). The median survival was 5, 5 and 3 generations for gene expression changes, 22G-RNA-based epimutations and chromatin-based epimutations respectively (Fig 2B). All three worm lineages displayed similar rates of epimutations (S2 Fig). This is orders of magnitude higher than the genome-wide rate of spontaneous nucleotide sequence change or insertion-deletion events in MA lines of C. elegans. [51–54].

[Figure omitted. See PDF.]

Fig 2. Rate of onset and survival of gene expression changes, chromatin and 22G-RNA-based epimutations.

A. Boxplots indicate percentage of loci per generation subject to gene expression changes (red), chromatin state changes (blue) or 22G-RNA level changes (gold). Box shows the interquartile range with horizontal line at the median; the whiskers extend to the furthest point no more than 1.5 times the interquartile range B. Survival plot shows step changes in survival of gene expression changes (red), 22G-RNA-based epimutations (gold) and chromatin-based epimutations (blue) persisting over a time frame of 20 generations (generational time points on X-axis). Y-axis showing survival probability. p-value calculated with log rank test.

https://doi.org/10.1371/journal.pgen.1010647.g002

3. Epimutations occur both with and without simultaneous changes in gene expression

We next asked whether changes in gene expression were accompanied by simultaneous changes in chromatin accessibility or 22G-RNA levels. We reasoned that simultaneous epimutations might be mechanistically linked to expression changes as outlined in the diagram (Fig 3A). To test this, we conducted a stepwise analysis of the association between chromatin-based epimutations and expression changes. We first looked at a background of all genes with annotated regulatory elements regardless of whether they had expression changes or not. The proportion of the genes which had expression changes was compared to the proportion of the genes which had chromatin-based epimutations. The overlap, in which genes had both expression changes and epimutations, was calculated (S56 Table). At this level, we found that chromatin-based epimutations rarely accompanied gene expression changes (Tests 1 & 2 in Table 1).

[Figure omitted. See PDF.]

Fig 3. Genes with inherited expression changes are enriched for simultaneous chromatin-based epimutations.

A. Chromatin-based epimutations which alter accessibility of regulatory elements to transcription factors (TF) and simultaneous gene expression changes may be inherited in parallel. B. Stacked bar-plot showing proportions of inherited (left bar) and non-inherited (right bar) gene expression changes accompanied by simultaneous chromatin-based epimutations which may themselves be inherited (pink) or non-inherited (blue). C. Stacked bar-plot showing proportions of inherited (left bar) and non-inherited (right bar) gene expression changes accompanied by simultaneous 22G-RNA-based epimutations which may themselves be inherited (blue) or non-inherited (beige). D. Examination of concordance (simultaneous and direction matched) of gene expression changes and chromatin-based epimutations. Y-axis shows associations tested. X-axis shows log2(Odds) of gene expression changes within each association being inherited. Odds ratios and p-values calculated with Fisher’s Exact Test. E. Examination of concordance (simultaneous and direction matched) of gene expression changes and changes in 22G-RNAs. Y-axis shows associations tested. X-axis shows log2(Odds) of gene expression changes within each association being inherited. Odds ratios and p-values calculated with Fisher’s Exact Test. p-value cut off for significance is 0.1.

https://doi.org/10.1371/journal.pgen.1010647.g003

We then looked more closely within the subset of genes which had expression changes. In this background we found a strong association between inherited expression changes and simultaneous chromatin-based epimutations (Fig 3B and Tests 3 & 4 in Table 1).

[Figure omitted. See PDF.]

Table 1. Stepwise analysis of association of gene expression changes and chromatin-based epimutations.

Cells in blue indicate significantly reduced odds of association, cells in grey indicate no significant association, cells in yellow indicate significantly increased odds of association. Odds ratios and p-values calculated with Fisher’s Exact Test.

https://doi.org/10.1371/journal.pgen.1010647.t001

Using the equivalent stepwise approach for 22G-RNAs (S57 Table), we found that changes in 22G-RNAs mapping antisense to protein-coding genes overlapped more strongly with changes in gene expression in a background of all genes (Test 1 in Table 2), and as with chromatin, genes showing heritable changes in gene expression were significantly enriched for simultaneous 22G-RNA-based epimutations (Fig 3C, Tests 3 & 4, Table 2).

[Figure omitted. See PDF.]

Table 2. Stepwise analysis of association of gene expression changes and 22G-RNA-based epimutations.

Cells in grey indicate no significant association, cells in yellow indicate significantly increased odds of association. Odds ratios and p-values calculated with Fisher’s Exact Test.

https://doi.org/10.1371/journal.pgen.1010647.t002

Therefore, although chromatin-based epimutations may occur without concurrent expression changes, on occasions where heritable expression changes are observed, simultaneous chromatin-based epimutations tend to accompany these. Chromatin-based epimutations may produce a favourable environment for changes in gene expression to be inherited but additional factors such as 22G-RNAs may be required to ‘license’ these.

We examined whether chromatin-based epimutations occurring in the same generation (simultaneous) and in the same direction (UP/UP or DOWN/DOWN) as gene expression changes i.e. ‘concordant’ had any effect on the odds that expression changes would be inherited, compared to simultaneous but ‘discordant’ (UP/DOWN, DOWN/UP) changes. Surprisingly, we found that genes with both concordant and discordant chromatin-based epimutations were similarly enriched for inherited expression changes (Fig 3D) which was in keeping with previous findings made from 22G-RNAs described above.

In the case of genes with simultaneous 22G-RNA-based epimutations and gene expression changes, we found that both concordant (UP/UP or DOWN/DOWN) and discordant (UP/DOWN, DOWN/UP) 22G-RNA-based epimutations were similarly enriched for inherited expression changes (Fig 3E). The enrichment for discordant 22G-RNA mediated epimutations was consistent with our previous observations [35] and the function of 22G-RNAs in gene silencing [55,56]

4. Chromatin environment affects the duration of epimutations

Previous work has shown that inherited epigenetic effects are typically short-lived. However, a small proportion of epimutations have been shown to persist over 10 or more generations [35]. In agreement with these results, we found that the duration of the majority of epimutations in chromatin, 22G-RNAs and gene expression was around 3–5 generations (Fig 2B). However, some epimutations lasted considerably longer, in some cases over 10 generations. We used K-means clustering with two clusters to identify a subset of loci that were subject to long-lived epimutations (Methods; S3A–S3C Fig) and showed that these strongly overlapped with loci derived by fitting the distribution to two normal distributions (Methods; S3D–S3F Fig). We compared these to heritable changes that were not in the long-lived set, which comprised loci with “short-lived” changes, and other loci that showed changes which were non-inherited. We examined the chromatin environments of these sets of genes.

Genes with expression changes and 22G-RNA-based epimutations were enriched in Regulated domains (marked by H3K27me3) and depleted from Active domains marked by H3K36me3 [39] (Fig 4A and 4C). The X chromosome, which is subject to H4K20 methylation as a result of dosage compensation and does not show high levels of either H3K27me3 or H3K36me3 [57–59], was also enriched for both non-inherited and heritable short-lived 22G-RNA-based epimutations (Fig 4C).

[Figure omitted. See PDF.]

Fig 4. Enrichment of gene expression changes and epimutations in different constitutive chromatin domains.

A. Bubble plot showing distribution of non-inherited (dark blue), short-lived (light green) and long-lived (light blue) gene expression changes in distinct chromatin domains. B. Bubble plot showing distribution of non-inherited, short-lived and long-lived chromatin-based epimutations in distinct chromatin domains. C. Bubble plot showing distribution of non-inherited, short-lived and long-lived 22G-RNA-based epimutations in distinct chromatin domains. For all plots, Y-axis shows constitutive chromatin domains investigated. X-axis shows log2(Odds) of enrichment. Odds ratios and p-values calculated with Fisher’s Exact Test with Bonferroni correction. p-value cut off for significance is 0.1.

https://doi.org/10.1371/journal.pgen.1010647.g004

Contrastingly, chromatin-based epimutations, regardless of duration, were enriched in Active domains (Fig 4B). Interestingly, we noticed that piRNA clusters, which occupy two broad regions on chromosome IV [60–62] were enriched for non-inherited chromatin-based epimutations (Fig 4B). This was surprising because piRNAs are predominantly produced from H3K27me3-marked regions [63]. Indeed, the proportion of chromatin epimutations in the piRNA clusters compared to the wider genome was greater than that of 22G-RNA-based epimutations and gene expression changes (S4 Fig). Stratification of epimutations by whether they were UP or DOWN did not affect the overall trends, suggesting that both directions contribute to the enrichments observed for all types of epimutation analysed (S5 Fig).

In previous work we investigated stability of 22G-RNA-based epimutations targeting genes associated with different small RNA pathways [35]. We showed that epimutations targeting the CSR-1 pathway, which is associated with gene activation [64], were unstable and short-lived, while those targeting the WAGO-1 pathway, which drives gene silencing [55], were highly stable, at times persisting beyond 10 generations. Having identified 22G-RNA-based epimutations with differing durations in this work, we then checked for differential associations with various small RNA pathways. We found the same enrichment pattern as in our previous work. Targets of non-inherited 22G-RNA-based epimutations strongly overlapped with CSR-1 pathway targets while targets of long-lived 22G-RNA-based epimutations strongly overlapped with WAGO-1 pathway targets (S6 Fig). Epimutations targeting both activating and silencing small RNA pathways is a possible explanation for why gene expression changes were similarly enriched for both concordant and discordant 22G-RNA-based epimutations (Fig 3E).

5. Long-lived epigenetic changes arise in genes functioning in adaptive responses

We next investigated whether epimutated genes with different durations of inheritance were enriched for particular biological functions. Genes with significantly long-lived expression changes were enriched for ontology terms relating to defence responses against bacteria and radiation (Fig 5A). This hints at a role in potentially adaptive processes [65] Compared to non-inherited and short-lived categories, long-lived chromatin-based epimutations at regulatory elements were strongly enriched across an array of pathways (Fig 5B), including regulation of defence response to bacterium. This finding is surprising given the previous observation (Fig 4A and 4B) that gene expression changes and chromatin epimutations are predominantly in distinct chromatin domains. Additionally, we were interested to see enrichment in germ-line stabilising SET domain pathways [66], and other pathways involving translation, protein localization and development. This supports the notion that chromatin-based epimutations establish conducive environments for heritable expression changes to occur which may have wide ranging phenotype effects. Genes with long-lived 22G-RNA mediated epimutations, however, did not show enrichment for terms relating to bacterial defence responses (Fig 5C).

[Figure omitted. See PDF.]

Fig 5. Genes with long-lived epigenetic changes are enriched for nematode defence functions.

A. Bubble plot showing ontology term enrichment of genes with non-inherited (dark blue), short-lived (light green) and long-lived (light blue) expression changes compared to genes without expression changes. Enrichment calculated with χ-squared test, top 10 results shown. X-axis shows log10(χ) for enrichment. Y-axis shows ontology terms. B. Bubble plot showing ontology term enrichment of regulatory element loci with long-lived, short-lived and non-inherited chromatin-based epimutations compared to genes without chromatin-based epimutations. Enrichment calculated with χ-squared test, top 10 results shown. X-axis shows log10(χ) for enrichment. Y-axis shows ontology terms. C. Bubble plot showing ontology term enrichment of genes with long-lived, short-lived and non-inherited 22G-RNA-based epimutations compared to genes without 22G-RNA-based epimutations. Enrichment calculated using Fisher’s Exact Test with Bonferroni correction, top 10 results shown. X-axis shows log10(Odds) of enrichment. Y-axis shows ontology terms. For all plots, p-value cut off for significance is 0.1.

https://doi.org/10.1371/journal.pgen.1010647.g005

We performed an ontology enrichment analysis of the set of genes that had inherited expression changes and simultaneous chromatin epimutations compared to genes which had non-inherited gene expression changes and simultaneous chromatin epimutations. Interestingly, this gene set was enriched for genes implicated in bacterial defence functions (S7 Fig).

6. Four P-glycoprotein genes are in a potential operon and have expression changes which may be explained by simultaneous 22G-RNA-based epimutations

Amongst genes with long-lived changes in expression, we identified 4 P-glycoprotein (pgp) genes. This was particularly interesting as PGP proteins have key functions in defence against xenobiotics, including a key role in drug resistance both in humans and nematodes [67–69]. Pgp genes 5, 6, 7 and 8 had long-lived epigenetic changes in all three lineages. The expression patterns for the 4 pgp genes were strongly similar within each lineage (Fig 6A and 6B). The fact that pgp expression was most often increased relative to the parental lineage could be due to chance, but it could also suggest that this might reflect an unstable epimutation that occurred in the population prior to the onset of the experiment. The organisation of pgp genes 5, 6, 7 & 8 in a 4 gene cluster has been described previously [70] but this cluster has not been annotated as an operon [71].

[Figure omitted. See PDF.]

Fig 6. A four gene cluster of P-glycoprotein genes has long-lived RNA-seq changes and may be co-regulated.

A. Chromatin states of three pgp-6 promoters and pgp-5, pgp-6, pgp-7, pgp-8 gene expression trends over 20 generations in Lineage A. Generational time points on X-axis. Z-score for epimutation status shown on Y-axis with 0 equivalent to PMA state. Horizontal thresholds indicate Z-score cut-offs; > 2.25 = UP epimutation, < - 2.25 = DOWN epimutation. B. Levels of 22G-RNAs targeting pgp-6 and pgp-7 over 20 generations in Lineage A. Generational time points on X-axis. Z-score for epimutation status shown on Y-axis with 0 equivalent to PMA state. Horizontal thresholds indicate Z-score cut-offs; > 2.25 = UP epimutation, < - 2.25 = DOWN epimutation.

https://doi.org/10.1371/journal.pgen.1010647.g006

Although not matching the threshold for epimutations, we observed that chromatin activity at the three regulatory elements for pgp-6 appeared to broadly track with some of the changes in gene expression for pgp genes 5, 6, 7, and 8 (Fig 6A). We wondered whether the genome-wide threshold for qualifying a chromatin state change as an epimutation under-estimated the extent of changes in the setting of the X chromosome due to its unique chromatin environment [59]. To address this, we derived Z-scores for X chromosome gene expression and chromatin state change loci using the same method as before but restricting the distribution of data points to those only derived from the X chromosome in order to produce X chromosome-specific Z-scores. However, the effect of adjusting for the specific context of the X chromosome was minimal (S8 Fig). Interestingly, when we examined the behaviour of 22G-RNA-based epimutations targeting genes within the putative pgp operon, a steep decline in 22G-RNAs occurred at precisely at the same time as the initial steep increase in gene expression (Fig 6B). In subsequent generations, the fluctuations in chromatin state changes appeared to correlate with gene expression changes (Fig 6A). However, it is important to note that neither the 22G-RNAs nor changes in chromatin correlate strongly with gene expression changes, suggesting that additional regulatory factors may be involved in regulating the heritable changes in gene expression that we observed at these loci.

Discussion

Here we have demonstrated that epimutations arise at the level of chromatin during long-term propagation of C. elegans under conditions of relaxed selection. These epimutations are heritable, and, in some cases, correspond to changes in gene expression. This expands our previous work demonstrating that 22G-RNA levels are subject to spontaneous short-lived epimutations in C. elegans [35] It also fits with recent results showing that short-term heritable variation in metabolite levels occurs in C. elegans lines propagated with limited population size [36]. Additionally, we show that epimutated genes display biased localisation across chromatin domains, and are enriched for genes with particular functional characteristics. Below, we discuss the implications of our findings for our understanding of the mechanism of epimutation formation and the potential consequences for evolutionary processes occurring within populations.

Our data shows clearly that spontaneous changes in gene expression arise during long-term propagation of C. elegans lineages. These changes are mostly short-lived, although occasionally long-lived epimutations appearing in more than 10 consecutive generations were observed. Changes in 22G-RNAs may account for around 26% of the heritable changes. Although significantly greater than expected by chance, changes in chromatin accessibility explain only around 12% of the remaining heritable changes. Interestingly, we also observe many changes in chromatin accessibility at regulatory elements that do not correspond to changes in the expression of nearby genes. Thus, changes in chromatin accessibility are to some extent decoupled from gene expression changes [72,73]. We note that in some cases this might reflect the fact that subtle changes, not sufficient to qualify as an epimutation, may be responsible- indeed we observed potential evidence for such changes at clustered pgp gene loci (Fig 6). Nevertheless, it seems likely that considerable changes in chromatin accessibility, which are heritable trans-generationally, can occur without changes in gene expression. Potentially the redundant structure of gene expression regulation in metazoans, whereby several regulatory elements combine to regulate each gene [74,75] could mean that spontaneous changes in one regulatory element are rarely sufficient for gene expression changes to follow. The fact that we find both concordant and discordant changes in chromatin and gene expression further emphasises that the coupling between chromatin accessibility and gene expression can be complex [73].

As we used ATAC-seq to map changes in chromatin accessibility, we are unable to determine the precise molecular events that occur during initiation or propagation of chromatin-based epimutations. However, the ability of the changes we detect to be sustained over multiple generations without obvious changes in transcription are consistent with the proposal that histone modifications can be transmitted by the ability of histone modifying enzymes to recognise the modifications that they themselves introduce [76]. Further investigation using ChIP-seq to identify changes in candidate histone modifications will be required to decipher the histone modifications involved and suggest possible modes of propagation.

An important remaining question is what underpins the remainder of the heritable changes in gene expression that cannot be explained by changes in chromatin or 22G-RNAs. It is possible that chromatin changes affecting specific histone modifications might not manifest as changes in ATAC-seq signal despite being linked to local changes in gene expression. Although C. elegans does not have cytosine DNA methylation [77] and no mechanism to transmit adenine DNA methylation is known [78], other sources of epigenetic information might be important, including prions and modifications to RNA molecules. Another interesting possibility is that RNA itself could act as an epigenetic signal. Random accumulation of additional mRNA molecules due to aberrant transcription during germ cell development might be sufficient to trigger a long-term propagation of a gene expression change, due to the ability of maternally deposited RNAs to “licence” subsequent expression [79].

An important controversy in the field of epigenetic inheritance is whether epigenetic variation produces phenotype variation that could be subject to selection [5,10]. Here we have shown that long-lived epigenetic changes arise in genes with key roles in xenobiotic detoxification. Long-lived chromatin-based epimutations also occurred in regulatory elements acting on xenobiotic response genes (Table 3). These genes showed different changes in different lineages and epimutations in these genes occurred at different time points within individual lineages. This suggests that these epimutations occur spontaneously even under stable environmental conditions, although it remains possible that small, environmental fluctuations outside of our control may be responsible for some of these alterations.

[Figure omitted. See PDF.]

Table 3. Four gene families with critical roles in xenobiotic defence have long-lived epimutations.

Xenobiotic defence genes are represented in the subset of genes with long-lived expression changes and long-lived chromatin-based epimutations but only one such gene had long-lived 22G-RNA-based epimutations.

https://doi.org/10.1371/journal.pgen.1010647.t003

Epimutations arising in these genes in stable environments could generate diversity in xenobiotic response profiles within a population. Xenobiotic defence gene families (Table 3) show a marked propensity for significant up and down fluctuations in expression over 20 generations while no such fluctuations are seen in a set of known housekeeping genes [81] (S9 Fig). Such fluctuations are reminiscent of a ‘bet hedging strategy’ whereby offspring have mixed fitness under different conditions. Stochastic and long-lived epimutations might enable epigenetic ‘anticipation’ of many possible environments [82,83]. The short-lived nature of epimutations may enable populations to partially adapt to environmental stress before adaptive genetic mutations have had time to emerge [5,10]. Further work is required to determine if innate epimutations could be subject to selection, and subsequently if epigenetically driven phenotypes become genomically encoded, through the effect of ‘genetic takeover’ [5,84].

Under conditions of xenobiotic stress, the propagation of pre-existing epimutations which are now advantageous in this setting is enhanced due to selection. On removal of the stressor, epimutations no longer provide a survival advantage, or indeed may be detrimental. Selection then ceases to drive propagation of the epimutation (Fig 7). Alternatively, the mechanism of epimutation inheritance itself may become more robust under stress, for example due to enhanced ‘self-templating’ of the epimutation [76], leading to population expansion of an epimutation which may or may not have adaptive benefit. Enhanced heritability returns to baseline when the stressor is removed. Either or both of these mechanisms might enable epimutations to contribute to adaptation to xenobiotic stress. It will be intriguing to test whether this occurs in nematode populations exposed to stress conditions.

[Figure omitted. See PDF.]

Fig 7. Two potential mechanisms for propagation of spontaneous epimutations under xenobiotic stress.

We propose two scenarios whereby epimutations may be propagated with short term reversibility under conditions of xenobiotic stress and the potential for adaptive consequences.

https://doi.org/10.1371/journal.pgen.1010647.g007

Methods

Construction of mutation accumulation lines

Caenorhabditis elegans N2 worms were grown at 20°C in a temperature-controlled incubator with time out of the incubator kept to a minimum. Worms were grown on Nematode Growth Medium (NGM) plates with an OP50 E. coli lawn. The use of Mutation Accumulation (MA) lines has been described previously [52]. MA lines are produced by restricting the effective population size by selecting the minimal number of individuals required to found each successive generation. Trial studies confirmed that in order to have sufficient embryo material for ATAC-seq, two worms were required to found each new generation. Each of three independent lineages (A, B & C) was therefore founded by picking two N2 L4 hermaphrodite worms from an expanded population. On day 4, two L4 hermaphrodite worms were picked from the founder plate (pre-mutation accumulation, or generation zero) to a new plate to produce the subsequent generation. On day 5, embryos were collected from the original plate. Lineages were propagated in this way for 20 generations. Embryos were collected at generation zero and then at alternate even numbered generations; 0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20. For embryo collection, worms and embryos were washed off plates into 15 ml falcon tubes with 0.1% Triton X. Plates were washed repeatedly to maximise collection of worms and embryos. The falcon tubes were centrifuged at 2000 RPM for 1 minute and the supernatant removed. Remaining embryos on the washed plates were loosened from agar by spraying 0.1% Triton X with the pipette angled against the plate. The volume of liquid containing additional washed off embryos was added to the worm pellet. Three further washes of the pellet were performed until the supernatant was clear to ensure that OP50 bacteria had been removed. Bleaching whole worms destroys the worm, leaving embryos intact. This meant that only embryonic material was retained for epigenetic and RNA-seq analysis. Worms were bleached with hypochlorite treatment in which 5.5 ml of bleach was added to each sample. Samples were vortexed continuously for 5 minutes with intermittent checking to ensure worm carcasses were dissolving and to avoid over exposure of embryos to bleach. Bleach was washed out rapidly once worm carcasses had disappeared with addition of M9 to dilute bleach, centrifugation to pellet undissolved embryo material, complete removal of supernatant and repeating this procedure for a total of three washes. Following the 3 post bleaching washes, the supernatant was removed from the embryo pellet. The pellet was resuspended in 1 ml 0.1% Triton X and then split into 1.5 ml microcentrifuge tubes follows: 25% of each sample was reserved for RNA library preparation and 100 ul TRIzol was added to inhibit RNase and maintain RNA integrity. 75% of each sample was reserved for ATAC-seq and frozen with no additive. Microcentrifuge tube were then submerged in liquid nitrogen before being stored at– 80°C.

Extraction of material for assessment of epimutations

RNA extraction.

Chloroform Isopropanol extraction of RNA was performed following standard protocols which had been optimised for our lab [35,85] as follows; embryos were collected in TRIzol as described above. Embryos were lysed through 5 freeze-cracking cycles in which embryos were frozen in liquid nitrogen and then thawed in a 37°C water bath. Tubes were vortexed for 5 minutes with 30 second pauses at 30 second intervals, and then incubated at room temperature for 5 minutes. This causes disruption of RNA-protein complexes. 200 μl chloroform per ml of TRIzol were added to samples followed by vigorous shaking. Samples were then incubated for two to three minutes at room temperature followed by centrifugation at a maximum speed at 4°C for 10 minutes. The top aqueous layer was aspirated and transferred to a new tube into which 1 μl glycogen and an equal volume of isopropanol was added. RNA was precipitated overnight at -20°C. After overnight precipitation, the samples were centrifuged for 1 hour at 4°C. The supernatant was removed and 500 μl of 75% ETOH added to the pellet. The samples were centrifuged at max speed for 10 minutes at 4°C. All ETOH was removed and the pellet allowed to air dry. The RNA pellet was resuspended in 10 μl H20. RNA sample concentration and quality was quantified on a 2200 TapeStation instrument using Agilent RNA screen tapes. An RNA Intengrity (RIN) score was derived. Samples were additionally quantified using Nanodrop.

Assay for Tranposase Accessible Chromatin (ATAC).

The protocol was adapted from Daugherty [48] after the original method for ATAC-seq [47]. Nuclei were extracted from embryo samples through cycles of freeze-cracking as follows. 200 μl of Nuclear Preparation Buffer was added to each microcentrifuge tube containing frozen embryos. The sample were submerged in liquid nitrogen for 90 seconds and then transferred to a 37°C water bath for 90 seconds. This was repeated three times. Embryo material was then transferred to a Wheaton glass tissue homogenizer in ice and was ground (3.5 grinds). The vessels containing the embryo material were covered and centrifuged for two minutes at 200 RCF at 4°C. The supernatant containing nuclei was transferred to a microcentrifuge tube on ice. The remaining embryo material was ground again, spun down and the supernatant again transferred to the collection tube. This cycle was repeated 4 times. The collection tubes were spun at 1000 RCF at 4°C for 10 minutes. The supernatant was discarded, leaving the nuclei pellet behind. To each sample nuclei pellet, Tagmentation enzyme and Tagmentation buffer were added (Illumina Tagment DNA enzyme and buffer). The samples were incubated for 30 minutes at 37°C on a Thermoshaker set to 580 RPM. Samples then underwent DNA clean up with a Qiagen MiniElute Reaction Cleanup Kit. The resulting DNA was eluted in 10 μl of EB buffer. PCR adapters from the Illumina Nextera DNA prep kit (Illumina) were selected and added to each sample along with PCR master mix. PCR was run with the following cycle parameters: 72°C for 5 minutes, 98°C for 30 seconds, 14 cycles (98°C for 10 seconds, 63°C for 30 seconds, 72°C for 1 minute). 4°C hold. Size selection isolates the desired DNA fragments and was achieved with magnetic AMP X beads (Beckman Coulter). First, excessively large DNA fragments were removed as follows; 25 μl AMP X bead solution was added to each sample. Samples were incubated at room temperature for 10 minutes. The PCR tubes containing the bead and sample mix were then put onto a magnetic rack for 5 minutes until the sample appeared clear. The samples were transferred without disturbing the magnetic beads to a new set of PCR tubes. To remove excessively small DNA fragments, the same procedure was followed, but the volume of AMP X bead solution added was adjusted to 65 μl. After incubating on the magnetic rack, the clear liquid was carefully removed and discarded. At this point, the beads had the desired library bound. The bead pellet was washed twice with 80% EtOH and then allowed to air dry at room temperature for two minutes. Excess EtOH droplets were removed with a pipette tip from inside each PCR tube. 22 μl of nuclease free water was added to each sample and the beads were washed down and dispersed into the liquid. The PCR tubes were returned to the magnetic rack, the DNA having eluted into the nuclease free water. The liquid was removed without disturbing the magnetic beads to a final collection tube which could then be frozen at -80°C. Validation of ATAC libraries was performed by quantification on the 2200 TapeStation using Agilent D1000 Screen tapes. ATAC libraries were also quantified on Qubit using the Qubit dsDNA HS Assay kit from Invitrogen.

Small RNA extraction.

The TruSeq Small RNA Prep Kit (Illumina) was used to prepare small RNA libraries from RNA pyrophosphohydrolase (Rpph) treated RNA. Rpph increases efficiency for 22G-RNA sequencing by converting 22G 5-prime triphosphate to monophosphate [86]. RNA was mixed with 1.5 μl RPPH enzyme and 2 μl 10X NEB Buffer 2 in a final volume of 20 μl. This was incubated at 37°C for 1 hour. The number of PCR cycles was increased from 11 to 15 as per protocol optimised previously [35] otherwise all steps were done according to manufacturer’s instructions. Size selection was performed using gel extraction with a 6% TBE gel (Invitrogen). Validation of 22G-RNA libraries was done on a 2200 TapeStation (Agilent).

Library construction

RNA library construction.

cDNA library preparation and PolyA-selected RNA sequencing was carried out in the MRC LMS Genomics Facility to generate PE100bp reads. A cDNA library from 31 samples was prepared using the TruSeq RNA Sample Prep kit (Illumina) following the manufacturer’s instructions. The cDNA library was then hybridised to the flowcell of the Illumina HiSeq 2500 and sequenced. The library was run on 3 separate lanes in order to achieve good coverage depth. The sequencing data were processed by the instrument’s Real Time Analysis (RTA) software application, version 1.18.64. De-multiplexing was done within the sequencing facility with CASAVA.

ATAC library construction.

For submission to LMS Genomics Facility, multiplexing of samples was performed, with 10–12 samples grouped into a pool. Samples were combined as 10 nM dilutions. Paired end sequencing of ATAC samples was done in the LMS Genomics Facility on a HiSeq instrument. The sequencing data were processed by the instrument’s Real Time Analysis (RTA) software application, version 1.18.66.3 using default filter and quality settings. Demultiplexing was done within the sequencing facility using CASAVA- 2.17, allowing zero mismatches.

Small RNA library construction.

For submission to LMS Genomics Facility, 3 nM dilutions were prepared from each small RNA sample. 3 nM sample dilutions were pooled at ~ 30 samples per pool with unique indexes. Sample pools were validated on TapeStation and quantified on Qubit using dsDNA reagents. MiSEQ Single Read sequencing (50 bp read length) was used to validate sample pools. Balancing or preparation of new pools was undertaken to improve the balance of samples in the pool as needed. High Output single read sequencing was performed on NExtSeq 2000 instrument in the LMS Genomics Facility. The max read length was 60 bpm. The sequencing data were processed by the instrument’s Real Time Analysis (RTA) software application, version 2.11.3.

Pre-processing and alignments

RNA-seq libraries.

RNA fastq files were aligned to the C. elegans genome using Bowtie2 [87] to produce a file containing the raw counts data with the following command;

bowtie2 \

-1 ${file1} \

-2 ${file2} \

-x (path to genome)/cel_genomeWS252 \

-S ${Prefix1}.sam \

-p 16

Sam files were sorted and converted to bam and then bed files. Sorted counts were intersected with coordinates of coding genes derived from the UCSC genome browser [88] using the following BEDTools commands [89].

intersectBed -c -a C.elegans_gene_names2.bed -b ${Prefix1}_2.bed \

-wa > ${Prefix1}_intersected_with_genes.bed

The files were simplified so that only the longest transcript was represented and 21U-RNAs were removed. In addition, all genes for which no expression was seen in any of the samples were removed. The read values were normalised using DESeq2 [90].

ATAC-seq libraries.

Trimming of reads to remove adapters was performed using FASTX Toolkit from the Hannon Lab (RRID:SCR_005534). Sequence alignment was carried out using Bowtie2. [87] to produce a file containing the raw counts data with the following command;

bowtie2 \

-1 ${file1}.cut.paired.fq

-2 ${file2}.cut.paired.fq \

-x (path to genome)/cel_genomeWS252 \

-S ${file1%.fastq}.sam \

-p 8

Sam files were sorted and converted to bam and then bed files. Sorted counts were intersected with coordinates for regulatory elements obtained from data produced by the Ahringer Group [50]. Using the following BEDTools commands [89]:

bedtools intersect -c -a Promoter_enhancer.gff \

-b ${Prefix1}_2.sorted.bed > ${Prefix1}_promoter_enhancer.bed

ATAC-seq across regulatory element loci was then annotated with C. elegans gene names obtained from the UCSC genome browser [88]. This resulted in quantification of the number of ATAC-seq reads mapping to each regulatory element. These counts were normalized using the DESeq2 pipeline [90].

Small RNA.

Secondary processing and demultiplexing of data was done within the sequencing facility with Bcl2fastq 2_2.20.0 (Illumina). Adapters were removed using FASTX Toolkit from the Hannon Lab (RRID:SCR_005534). For 22G-RNA mapping the un-collapsed reads were aligned to the Ce11 genome using Bowtie [91] allowing 0 mismatches. 22G-RNAs were selected using a perl script and antisense reads mapping to transcribed regions corresponding to protein-coding genes were selected using Bedtools [89]. Reads mapping to each transcribed gene were summed. Data was normalized using DESeq2. For piRNAs and miRNAs, reads were collapsed using FASTX Toolkit from the Hannon Lab (RRID:SCR_005534) and exact matching to annotated piRNAs [63] and miRNAs [92] was performed using the seqinr package. Simulations to find the expected number of epimutations with a cutoff of a Z-score of 2.25 were performed by randomly populating a table with epimutations at the observed average rate and counting the number that were inherited for 1 generation.

Identification of epimutations

Following normalisation of count values, a linear model was created in which log2(fold changes) were plotted against the log2(mean counts) for each gene in line with the approach taken in previous work to identify epimutations in small RNAs [35]. For each data point, a Z-score was derived by subtracting the mean of the residuals from the linear model from each individual residual value and dividing the output by the standard deviation of the residuals. Regions with a Z-score greater or less than 2.25 were defined as showing differential accessibility. Regions with a Z-score greater than 2.25 were annotated as ‘Up’ epimutations and regions with a Z-score less than 2.25 were annotated as ‘Down’ epimutations.

Test for epimutation inheritance

Simulation of data with random distribution of epimutations across a range of Z-scores was performed as described above. Briefly, we applied each value from a range of Z-score cut offs (1–3) to the ATAC-seq Z-score data set to identify epimutations and produce a corresponding set of binarized data tables in which either 1 (Up epimutation) or -1 (Down epimutation) indicated an epimutated locus at a specific generational time point. As would be expected, sequentially raising the threshold for data points to qualify as epimutations imposed increased stringency thereby reducing the number of points which qualified as epimutations, while the lowest and therefore most permissive threshold qualified a larger number of points as epimutations. For each binarized data table, we then produced a large number of simulated data tables for which the rate of epimutations arising per generation in the real data set was preserved, but the genomic sites undergoing epimutation was determined through random re-assortment of epimutations across loci. Random re-assortment was restricted to loci deemed to be ‘epimutable’, i.e. those with epimutations in the real data set. This was to prevent the confounding effect of incorporating loci which are never epimutated. Therefore, we examined whether epimutations arose in consecutive generational time points in the same locus, as would be expected if a mechanism for inheritance is active. We were able to reject the null hypothesis that epimutations were no more likely to arise in consecutive generations in real data compared to simulated data.

Removal of any epimutations which potentially could have been the result of environmental factors

We recognised that despite tightly controlled conditions, minor changes in the experiment environment could occur as the mutation accumulation experiments ran over an approximate twelve-week period. As the experiment was conducted in a standard laboratory setting, slight changes to UV light levels and ambient temperature over the course of the experiment could potentially affect worms when outside of the incubator. These could constitute external stimuli which could potentially have induced epimutations. Therefore, these would not represent spontaneous changes. We reasoned that such effects would feature in all generations similarly as the mutation accumulation lines were run in parallel and handled outside of the incubator at the same time points. In order to remove epimutations which may have been induced, we excluded any epimutations which had the same onset time and duration in all three lineages.

Quantification of epimutation duration

Our data derived from a pre-mutation founder generation and alternate generation time points (even numbers; generations 2, 4, 6, 8, 10, 12, 14, 16, 18 and 20) from a 20 generation lineage. This required us to form certain assumptions regarding the epimutation status of the intervening generations for which no data was collected. In addition, certain generations were missing from several of the data sets. RNA-seq and small RNA-seq data sets for MA line B lacked generation 4, and for MA line C lacked generation 14. ATAC-seq data set for MA line B lacked generation 4, and for MA line C lacked generation 8. Both ATAC-seq and RNA-seq data sets were complete for MA line A. The reason for missing generations was insufficient yield of RNA or cDNA from these generations, therefore these samples were omitted from sequencing.

When examining the length of epimutations, we assumed that intervening odd numbered generations (1, 3, 5, 7, 9, 11, 13, 15, 17, 19) were not epimutated unless the even numbered generations either side were epimutated. In order to be considered to be inherited, we required that two consecutive even numbered generations displayed an epimutation. This is therefore a conservative approach to detecting inherited epimutations. In measuring the duration of epimutations we accounted for the intervening odd number generations which were assumed to be epimutated. For example, for the germline coding promoter which maps to gene Y65B4BL.7 (coordinates: chr1:510778:510928), the epimutation status over the course of Lineage A was as follows (Table 4):

[Figure omitted. See PDF.]

Table 4. Y65B4BL.7 Chromatin-based epimutation status over 20 generations in Lineage A.

https://doi.org/10.1371/journal.pgen.1010647.t004

This locus is recorded as having 3 separate epimutations. 1 putative epimutation (down) which arose in Generation 2 and had a length of 1, i.e. non-inherited. 1 epimutation (down) which arose in Generation 10 and had a length of 5 (1 + 4 additional generations until disappeared in Generation 14), and finally a non-inherited putative epimutation (down) which arose in Generation 18 and had a length of 1. Therefore, for the epimutation running across Generations 10, 12 and 14, we assumed that epimutations are present in the intervening odd numbered generations 11 and 13.

To be considered as within a single epimutation run, we required that epimutation data points have the same directionality, i.e. all -1 or all 1. A transition from -1 to 1 or vice versa was considered to be a termination of the preceding epimutation and commencement of a new epimutation run. We reasoned that epimutation inheritance would also include inheritance of direction in which chromatin state changed as molecular processes establishing heterochromatin are antagonistic to those dismantling heterochromatin.

We applied a linear model using the lme4 package in RStudio [93] to obtain estimates for the effect of each gene on epimutation duration. We used K-means clustering of the estimates to identify the subset of genes predicted to cause ‘long-lived’ changes in RNA-seq, ATAC-seq and small RNA data sets. The remaining loci in each data set were separated into those with short-lived inherited changes and those with non-inherited changes. The three ‘test’ lists of genes (long-lived, short-lived and non-inherited) were obtained for RNA-seq, ATAC-seq and small RNA data sets and compared to genes without epimutations in subsequent gene set enrichment analyses (below).

Integration of gene expression change and chromatin epimutation data

In order to match up gene expression changes with chromatin states at associated regulatory elements we were required to deal with the unequal ratio of genes to regulatory elements. As would be expected, our data comprised a larger number of regulatory elements than genes. Certain genes have multiple associated coding promoters and putative enhancers as defined previously [50,94–96]. The C. elegans lifespan regulating gene sesn-1 [97], for example, has 3 annotated coding promoters and 12 annotated putative enhancers.

We integrated ATAC-seq data with RNA-seq data in a gene-centric manner. For each gene we recorded its expression state and if expression state change was inherited and the presence of any epimutated regulatory elements. Inheritance of epimutation was recorded if an epimutation was present in at least two consecutive generations for any single regulatory element. We also looked for evidence of simultaneous epimutations in associated regulatory elements. For example, for C. elegans detoxification gene ugt-26 [69] the expression status (Table 5) and associated regulatory element epimutation status (Table 6) over the course of Lineage A were as follows:

[Figure omitted. See PDF.]

Table 5. Gene expression status for ugt-26.

https://doi.org/10.1371/journal.pgen.1010647.t005

[Figure omitted. See PDF.]

Table 6. Chromatin-based epimutation status for ugt-26.

https://doi.org/10.1371/journal.pgen.1010647.t006

Therefore, ugt-26 was recorded as having both inherited expression changes (onset Generation 8 running to at least Generation 20) which were synchronised with chromatin state changes in an associated regulatory element, which are themselves inherited (onset Generation 8, offset Generation 10). Importantly, we did not require that simultaneous epimutations persist for the same duration as inherited expression changes. This is so as to take into account the potential that once established, epigenetic changes may persist in the absence of the initiating mechanism [41].

Gene set enrichment analyses

Lists of epimutated or expression changed genes (test lists) were uploaded to WormEnrichr [98,99] and the number of genes annotated with each identified ontology term was compared between the test lists and background sets of genes which had no epimutations or expression changes. The background gene set was too large to be uploaded in entirety to WormEnrichr, therefore samples were taken from the list (50 samples each of 6000 genes), analysed for ontology enrichments and compared to the test list. We required a minimum of 5 genes to be present in either the test or background sample gene list for the ontology term to be included in the analysis. Enrichment of test list genes for ontology terms compared to genes sampled from the background list was calculated using a Chi Squared Test. When background lists were small enough to be uploaded in entirety to WormEnrichr, enrichment was calculated with Fisher’s Exact Test.

Production of plots

Plots were made in R Studio (Rstudio Team (2020). Rstudio: Integrated Development for R. Rstudio, PBC, Boston, MA.) with aesthetic editing in Adobe Illustrator (Adobe Illustrator software). Schematic diagrams were made in BioRender (BioRender.com) and Adobe Illustrator.

Supporting information

S1 Fig. 22G-RNA class of small RNAs have greatest proportion of heritable epimutations compared to miRNAs and piRNAs, and this is greater than what would be expected by chance.

For each small RNA class (22G-RNAs, mi-RNAs and piRNAS), comparison of proportion of changes exceeding Z-score threshold (2.25, - 2.25) that are inherited to at least two subsequent generations between observed data (dot) and observations derived from 10000 random simulations (Methods).

https://doi.org/10.1371/journal.pgen.1010647.s001

(TIF)

S2 Fig. Epimutation rates are similar across three independent worm lineages.

A & B. Gene expression changes. Percentage of genes showing significantly increased (A) and decreased (B) expression per generation. C & D. Chromatin state changes. Percentage of regulatory elements showing significantly increased (C) and decreased (D) chromatin accessibility per generation. E & F. 22G-RNA level changes. Percentage of genes in each lineage showing significantly increased (E) and decreased (F) antisense 22G-RNA levels per generation. For all plots, box shows the interquartile range with horizontal line at the median; the whiskers extend to the furthest point no more than 1.5 times the interquartile range. Kruskal-Wallis rank sum test gives non-significant p-values > 0.1 for comparison of worm lineages A, B and C within the categories of data type (Gene expression, Chromatin state, 22G-RNA level) and epimutation direction (Up transitions or Down transitions).

https://doi.org/10.1371/journal.pgen.1010647.s002

(TIF)

S3 Fig. Deriving subsets of genes with long-lived expression changes and epimutations.

K-means clustering was done to identify the subset of genes predicted from a linear model to have significantly long-lived expression changes (A), chromatin-based epimutations (B), and 22G-RNA-based epimutations (C). The long-lived gene sets were derived independently using an alternative method based on expectation maximisation (EM) clustering using MixTools package in R (‘MixEM method’). The overlap between long-lived gene sets derived through both methods is shown for genes with expression changes (D), loci with chromatin-based epimutations (E), and loci with 22G-RNA-based epimutations (F).

https://doi.org/10.1371/journal.pgen.1010647.s003

(TIF)

S4 Fig. Small RNA pathways in chromatin and 22G-RNA-based epimutations.

From left to right. Percentage of genes with expression changes (red bar), percentage of regulatory loci with chromatin epimutations (green bar) and percentage of genes with antisense 22G-RNA epimutations (gold bar) which overlap with piRNA cluster genes.

https://doi.org/10.1371/journal.pgen.1010647.s004

(TIF)

S5 Fig. Enrichment of direction stratified gene expression changes and epimutations in different constitutive chromatin domains.

A. Bubble plot showing distribution of UP and DOWN non-inherited short-lived and long-lived gene expression changes in distinct chromatin domains. B. Bubble plot showing distribution of UP and DOWN non-inherited short-lived and long-lived chromatin-based epimutations in distinct chromatin domains. C. Bubble plot showing distribution of UP and DOWN non-inherited short-lived and long-lived 22G-RNA-based epimutations in distinct chromatin domains. For all plots, Y-axis shows constitutive chromatin domains investigated. X-axis shows log2(Odds) of enrichment. Odds ratios and p-values calculated with Fisher’s Exact Test with Bonferroni correction. p-value cut off for significance is 0.1.

https://doi.org/10.1371/journal.pgen.1010647.s005

(TIF)

S6 Fig. Differential enrichment of 22G-RNA-based epimutations in small-RNA pathways.

Bubble plot showing enrichment of different small RNA pathways for 22G-RNA-based epimutations of different durations. Y-axis shows specific small-RNA pathway associated proteins. X-axis shows log2(Odds) of enrichment. Odds ratios and p-values are calculated using Fisher’s Exact Test with Bonferroni Correction. p-value cut off for significance is 0.1.

https://doi.org/10.1371/journal.pgen.1010647.s006

(PDF)

S7 Fig. Genes with inherited expression changes and simultaneous chromatin-based epimutations are enriched for functions involved in defence response to bacteria.

Y-axis shows ontology terms. X-axis shows log10(Odds) of enrichment. Top 10 results shown. Odds ratios and p-values are calculated using Fisher’s Exact Test with Bonferroni Correction. p-value cut off for significance is 0.1.

https://doi.org/10.1371/journal.pgen.1010647.s007

(TIF)

S8 Fig. Chromatin fluctuations at X chromosome are not altered by adjusting normalisation thresholds.

pgp-5, pgp-6, pgp-7, pgp-8 expression and regulatory element chromatin state over 20 generations. Normalisation of ATAC-seq counts is restricted to X chromosome. Data are for worm lineage A. Generational time points on X-axis. Z-score for epimutation status shown on Y-axis with 0 equivalent to PMA state. Horizontal thresholds indicate Z-score cut-offs; > 2.25 = UP epimutation and < - 2.25 = DOWN epimutation.

https://doi.org/10.1371/journal.pgen.1010647.s008

(TIF)

S9 Fig. Xenobiotic defence genes fluctuate under relaxed conditions while housekeeping genes do not.

A. Expression over 20 generations across 3 independent worm lineages in housekeeping genes was compared to that of families of xenobiotic defence genes identified to contain long-lived epimutations. Housekeeping genes n = 13, pgp genes n = 14, nhr genes n = 268, ugt genes n = 66, cyp genes n = 74. B. As in (A) but the 4 gene pgp cluster (pgp-5, pgp-6, pgp-7, pgp-8) is removed showing that removing this gene cluster removes bias of pgp gene family towards Up epimutations as seen in (A).

https://doi.org/10.1371/journal.pgen.1010647.s009

(TIF)

S1 Table. Worm Lineage A gene expression, table of Z scores.

File name: S1_Table.csv.

https://doi.org/10.1371/journal.pgen.1010647.s010

(CSV)

S2 Table. Worm Lineage B gene expression, table of Z scores.

File name: S2_Table.csv.

https://doi.org/10.1371/journal.pgen.1010647.s011

(CSV)

S3 Table. Worm Lineage C gene expression, table of Z scores.

File name: S3_Table.csv.

https://doi.org/10.1371/journal.pgen.1010647.s012

(CSV)

S4 Table. Worm Lineage A 22G-RNA levels, table of Z scores.

File name: S4_Table.csv.

https://doi.org/10.1371/journal.pgen.1010647.s013

(CSV)

S5 Table. Worm Lineage B 22G-RNA levels, table of Z scores.

File name: S5_Table.csv.

https://doi.org/10.1371/journal.pgen.1010647.s014

(CSV)

S6 Table. Worm Lineage C 22G-RNA levels, table of Z scores.

File name: S6 Table.csv.

https://doi.org/10.1371/journal.pgen.1010647.s015

(CSV)

S7 Table. Worm Lineage A ATAC counts, table of Z scores.

File name: S7 Table.csv.

https://doi.org/10.1371/journal.pgen.1010647.s016

(CSV)

S8 Table. Worm Lineage B ATAC counts, table of Z scores.

File name: S8 Table.csv.

https://doi.org/10.1371/journal.pgen.1010647.s017

(CSV)

S9 Table. Worm Lineage C ATAC counts, table of Z scores.

File name: S9_Table.csv.

https://doi.org/10.1371/journal.pgen.1010647.s018

(CSV)

S10 Table. Worm Lineage A gene expression Table of Up (1) and Down (-1) epimutations.

File name: S10_Table.csv.

https://doi.org/10.1371/journal.pgen.1010647.s019

(CSV)

S11 Table. Worm Lineage B gene expression Table of Up (1) and Down (-1) epimutations.

File name: S11_Table.csv.

https://doi.org/10.1371/journal.pgen.1010647.s020

(CSV)

S12 Table. Worm Lineage C gene expression Table of Up (1) and Down (-1) epimutations.

File name: S12_Table.csv.

https://doi.org/10.1371/journal.pgen.1010647.s021

(CSV)

S13 Table. Worm Lineage A ATAC Table of Up (1) and Down (-1) epimutations.

File name: S13_Table.csv.

https://doi.org/10.1371/journal.pgen.1010647.s022

(CSV)

S14 Table. Worm Lineage B ATAC Table of Up (1) and Down (-1) epimutations.

File name: S14_Table.csv.

https://doi.org/10.1371/journal.pgen.1010647.s023

(CSV)

S15 Table. Worm Lineage C ATAC Table of Up (1) and Down (-1) epimutations.

File name: S15_Table.csv.

https://doi.org/10.1371/journal.pgen.1010647.s024

(CSV)

S16 Table. Worm Lineage A 22G-RNA Table of Up (1) and Down (-1) epimutations.

File name: S16_Table.csv.

https://doi.org/10.1371/journal.pgen.1010647.s025

(CSV)

S17 Table. Worm Lineage B 22G-RNA Table of Up (1) and Down (-1) epimutations.

File name: S17_Table.csv.

https://doi.org/10.1371/journal.pgen.1010647.s026

(CSV)

S18 Table. Worm Lineage C 22G-RNA Table of Up (1) and Down (-1) epimutations.

File name: S18_Table.csv.

https://doi.org/10.1371/journal.pgen.1010647.s027

(CSV)

S19 Table. Table of values used in Fig 1C.

File name: S19_Table.csv.

https://doi.org/10.1371/journal.pgen.1010647.s028

(CSV)

S20 Table. Table of values used in Fig 2A.

File name: S20_Table.csv.

https://doi.org/10.1371/journal.pgen.1010647.s029

(CSV)

S21 Table. Table of values used in Fig 2B.

File name: S21_Table.csv.

https://doi.org/10.1371/journal.pgen.1010647.s030

(CSV)

S22 Table. Table of values used in Fig 3B.

File name: S22_Table.csv.

https://doi.org/10.1371/journal.pgen.1010647.s031

(CSV)

S23 Table. Table of values used in Fig 3C.

File name: S23_Table.csv.

https://doi.org/10.1371/journal.pgen.1010647.s032

(CSV)

S24 Table. Table of values used in Fig 3D.

File name: S24_Table.csv.

https://doi.org/10.1371/journal.pgen.1010647.s033

(CSV)

S25 Table. Table of values used in Fig 3E.

File name: S25_Table.csv.

https://doi.org/10.1371/journal.pgen.1010647.s034

(CSV)

S26 Table. Table of values used in Fig 4A.

File name: S26_Table.csv.

https://doi.org/10.1371/journal.pgen.1010647.s035

(CSV)

S27 Table. Table of values used in Fig 4B.

File name: S27_Table.csv.

https://doi.org/10.1371/journal.pgen.1010647.s036

(CSV)

S28 Table. Table of values used in Fig 4C.

File name: S28_Table.csv.

https://doi.org/10.1371/journal.pgen.1010647.s037

(CSV)

S29 Table. Table of values used in Fig 5A.

File name: S29_Table.csv.

https://doi.org/10.1371/journal.pgen.1010647.s038

(CSV)

S30 Table. Table of values used in Fig 5B.

File name: S30_Table.csv.

https://doi.org/10.1371/journal.pgen.1010647.s039

(CSV)

S31 Table. Table of values used in Fig 5C.

File name: S31_Table.csv.

https://doi.org/10.1371/journal.pgen.1010647.s040

(CSV)

S32 Table. Table of values for gene expression used in Fig 6A and 6B.

File name: S32_Table.csv.

https://doi.org/10.1371/journal.pgen.1010647.s041

(CSV)

S33 Table. Table of values for chromatin state used in Fig 6A.

File name: S33_Table.csv.

https://doi.org/10.1371/journal.pgen.1010647.s042

(CSV)

S34 Table. Table of values for 22G-RNA level used in Fig 6B.

File name: S34_Table.csv.

https://doi.org/10.1371/journal.pgen.1010647.s043

(CSV)

S35 Table. Table of values used in S2A Fig File name: S35_Table.csv.

https://doi.org/10.1371/journal.pgen.1010647.s044

(CSV)

S36 Table. Table of values used in S2B Fig File name: S36_Table.csv.

https://doi.org/10.1371/journal.pgen.1010647.s045

(CSV)

S37 Table. Table of values used in S2C Fig File name: S37_Table.csv.

https://doi.org/10.1371/journal.pgen.1010647.s046

(CSV)

S38 Table. Table of values used in S2D Fig File name: S38_Table.csv.

https://doi.org/10.1371/journal.pgen.1010647.s047

(CSV)

S39 Table. Table of values used in S2E Fig File name: S39_Table.csv.

https://doi.org/10.1371/journal.pgen.1010647.s048

(CSV)

S40 Table. Table of values used in S2F Fig File name: S40_Table.csv.

https://doi.org/10.1371/journal.pgen.1010647.s049

(CSV)

S41 Table. Table of values used in S3A Fig File name: S41_Table.csv.

https://doi.org/10.1371/journal.pgen.1010647.s050

(CSV)

S42 Table. Table of values used in S3B Fig File name: S42_Table.csv.

https://doi.org/10.1371/journal.pgen.1010647.s051

(CSV)

S43 Table. Table of values used in S3C Fig File name: S43_Table.csv.

https://doi.org/10.1371/journal.pgen.1010647.s052

(CSV)

S44 Table. Table of values used in S3D Fig File name: S44_Table.csv.

https://doi.org/10.1371/journal.pgen.1010647.s053

(CSV)

S45 Table. Table of values used in S3E Fig File name: S45_Table.csv.

https://doi.org/10.1371/journal.pgen.1010647.s054

(CSV)

S46 Table. Table of values used in S3F Fig File name: S46_Table.csv.

https://doi.org/10.1371/journal.pgen.1010647.s055

(CSV)

S47 Table. Table of values used in S4 Fig File name: S47_Table.csv.

https://doi.org/10.1371/journal.pgen.1010647.s056

(CSV)

S48 Table. Table of values used in S5A Fig File name: S48_Table.csv.

https://doi.org/10.1371/journal.pgen.1010647.s057

(CSV)

S49 Table. Table of values used in S5B Fig File name: S49_Table.csv.

https://doi.org/10.1371/journal.pgen.1010647.s058

(CSV)

S50 Table. Table of values used in S5C Fig File name: S50_Table.csv.

https://doi.org/10.1371/journal.pgen.1010647.s059

(CSV)

S51 Table. Table of values used in S6 Fig File name: S51_Table.csv.

https://doi.org/10.1371/journal.pgen.1010647.s060

(CSV)

S52 Table. Table of values used in S7 Fig File name: S52_Table.csv.

https://doi.org/10.1371/journal.pgen.1010647.s061

(CSV)

S53 Table. Table of values for gene expression used in S8 Fig File name: S53_Table.csv.

https://doi.org/10.1371/journal.pgen.1010647.s062

(CSV)

S54 Table. Table of values for chromatin state used in S8 Fig File name: S54_Table.csv.

https://doi.org/10.1371/journal.pgen.1010647.s063

(CSV)

S55 Table. Table of values used in S9 Fig File name: S55_Table.csv.

https://doi.org/10.1371/journal.pgen.1010647.s064

(CSV)

S56 Table. Table of values used in Table 1.

File name: S56_Table.csv.

https://doi.org/10.1371/journal.pgen.1010647.s065

(CSV)

S57 Table. Table of values used in Table 2.

File name: S57_Table.csv.

https://doi.org/10.1371/journal.pgen.1010647.s066

(CSV)

Acknowledgments

We thank Prof. Vahid Shahrezaei (Department of Biomathematics, Imperial College London) for guidance on the statistical approach used to identify epimutations, Dr. Subhanita Ghosh for assistance in developing the ATAC-seq protocol and members of the Epigenetics and Evolution group for helpful comments on the manuscript.

Citation: Wilson R, Le Bourgeois M, Perez M, Sarkies P (2023) Fluctuations in chromatin state at regulatory loci occur spontaneously under relaxed selection and are associated with epigenetically inherited variation in C. elegans gene expression. PLoS Genet 19(3): e1010647. https://doi.org/10.1371/journal.pgen.1010647

References

1. Heard E, Martienssen RA. Transgenerational Epigenetic Inheritance: Myths and Mechanisms. Cell. 2014;157: 95–109. pmid:24679529

2. Jablonka E, Raz G. Transgenerational Epigenetic Inheritance: Prevalence, Mechanisms, and Implications for the Study of Heredity and Evolution. Q Rev Biol. 2009;84: 131–176. pmid:19606595

3. Lacal I, Ventura R. Epigenetic Inheritance: Concepts, Mechanisms and Perspectives. Front Mol Neurosci. 2018;11: 292. pmid:30323739

4. Miska EA, Ferguson-Smith AC. Transgenerational inheritance: Models and mechanisms of non–DNA sequence–based inheritance. Science. 2016;354: 59–63. pmid:27846492

5. Sarkies P. Molecular mechanisms of epigenetic inheritance: Possible evolutionary implications. Semin Cell Dev Biol. 2020;97: 106–115. pmid:31228598

6. Harvey ZH, Chakravarty AK, Futia RA, Jarosz DF. A Prion Epigenetic Switch Establishes an Active Chromatin State. Cell. 2020;180: 928–940.e14. pmid:32109413

7. Klironomos FD, Berg J, Collins S. How epigenetic mutations can affect genetic evolution: Model and mechanism: Problems & Paradigms. BioEssays. 2013;35: 571–578. pmid:23580343

8. Nilsson EE, Maamar MB, Skinner MK. Environmentally Induced Epigenetic Transgenerational Inheritance and the Weismann Barrier: The Dawn of Neo-Lamarckian Theory. J Dev Biol. 2020;8: 28. pmid:33291540

9. Richards EJ. Inherited epigenetic variation—revisiting soft inheritance. Nat Rev Genet. 2006;7: 395–401. pmid:16534512

10. Burggren W. Epigenetic Inheritance and Its Role in Evolutionary Biology: Re-Evaluation and New Perspectives. Biology. 2016;5: 24. pmid:27231949

11. Holliday R. Epigenetics: A Historical Overview. Epigenetics. 2006;1: 76–80. pmid:17998809

12. Chong S, Vickaryous N, Ashe A, Zamudio N, Youngson N, Hemley S, et al. Modifiers of epigenetic reprogramming show paternal effects in the mouse. Nat Genet. 2007;39: 614–622. pmid:17450140

13. Horsthemke B. Heritable germline epimutations in humans. Nat Genet. 2007;39: 573–574. pmid:17460680

14. Calo S, Shertz-Wall C, Lee SC, Bastidas RJ, Nicolás FE, Granek JA, et al. Antifungal drug resistance evoked via RNAi-dependent epimutations. Nature. 2014;513: 555–558. pmid:25079329

15. Gutierrez H, Taghizada B, Meneghini MD. Nutritional and Meiotic Induction of Heritable Stress Resistant States in Budding Yeast. Microb Cell. 2018;5: 511–521. pmid:30483522

16. Jirtle RL, Skinner MK. Environmental epigenomics and disease susceptibility. Nat Rev Genet. 2007;8: 253–262. pmid:17363974

17. Torres-Garcia S, Yaseen I, Shukla M, Audergon PNCB, White SA, Pidoux AL, et al. Epigenetic gene silencing by heterochromatin primes fungal resistance. Nature. 2020;585: 453–458. pmid:32908306

18. van der Graaf A, Wardenaar R, Neumann DA, Taudt A, Shaw RG, Jansen RC, et al. Rate, spectrum, and evolutionary dynamics of spontaneous epimutations. Proc Natl Acad Sci. 2015;112: 6676–6681. pmid:25964364

19. Yao N, Schmitz RJ, Johannes F. Epimutations Define a Fast-Ticking Molecular Clock in Plants. Trends Genet. 2021;37: 699–710. pmid:34016450

20. Katju V, Bergthorsson U. Old Trade, New Tricks: Insights into the Spontaneous Mutation Process from the Partnering of Classical Mutation Accumulation Experiments with High-Throughput Genomic Approaches. Makova K, editor. Genome Biol Evol. 2019;11: 136–165. pmid:30476040

21. Ohnishi O. 529.pdf. Genetics. 1977.

22. Teotónio H, Estes S, Phillips PC, Baer CF. Experimental Evolution with Caenorhabditis Nematodes. Genetics. 2017;206: 691–716. pmid:28592504

23. Lynch M, Ackerman MS, Gout J-F, Long H, Sung W, Thomas WK, et al. Genetic drift, selection and the evolution of the mutation rate. Nat Rev Genet. 2016;17: 704–714. pmid:27739533

24. Saxena AS, Salomon MP, Matsuba C, Yeh S-D, Baer CF. Evolution of the Mutational Process under Relaxed Selection in Caenorhabditis elegans. Ruvinsky I, editor. Mol Biol Evol. 2019;36: 239–251. pmid:30445510

25. Becker C, Hagmann J, Müller J, Koenig D, Stegle O, Borgwardt K, et al. Spontaneous epigenetic variation in the Arabidopsis thaliana methylome. Nature. 2011;480: 245–249. pmid:22057020

26. Hazarika RR, Serra M, Zhang Z, Zhang Y, Schmitz RJ, Johannes F. Molecular properties of epimutation hotspots. Nat Plants. 2022;8: 146–156. pmid:35087209

27. Schmitz RJ, Schultz MD, Lewsey MG, O’Malley RC, Urich MA, Libiger O, et al. Transgenerational Epigenetic Instability Is a Source of Novel Methylation Variants. Science. 2011;334: 369–373. pmid:21921155

28. Drinnenberg IA, Berger F, Elsässer SJ, Andersen PR, Ausió J, Bickmore WA, et al. EvoChromo: towards a synthesis of chromatin biology and evolution. Development. 2019;146: dev178962. pmid:31558570

29. Johannes F, Schmitz RJ. Spontaneous epimutations in plants. New Phytol. 2019;221: 1253–1259. pmid:30216456

30. Cecere G. Small RNAs in epigenetic inheritance: from mechanisms to trait transmission. FEBS Lett. 2021;595: 2953–2977. pmid:34671979

31. Frolows N, Ashe A. Small RNAs and chromatin in the multigenerational epigenetic landscape of Caenorhabditis elegans. Philos Trans R Soc B Biol Sci. 2021;376: rstb.2020.0112, 20200112. pmid:33866817

32. Houri-Zeevi L, Korem Kohanim Y, Antonova O, Rechavi O. Three Rules Explain Transgenerational Small RNA Inheritance in C. elegans. Cell. 2020;182: 1186–1197.e12. pmid:32841602

33. Posner R, Toker IA, Antonova O, Star E, Anava S, Azmon E, et al. Neuronal Small RNAs Control Behavior Transgenerationally. Cell. 2019;177: 1814–1826.e15. pmid:31178120

34. Rechavi O, Lev I. Principles of Transgenerational Small RNA Inheritance in Caenorhabditis elegans. Curr Biol. 2017;27: R720–R730. pmid:28743023

35. Beltran T, Shahrezaei V, Katju V, Sarkies P. Epimutations driven by small RNAs arise frequently but most have limited duration in Caenorhabditis elegans. Nat Ecol Evol. 2020;4: 1539–1548. pmid:32868918

36. Johnson LM, Smith OJ, Hahn D, Baer CF. Short-term heritable variation overwhelms two hundred generations of mutational variance for metabolic traits in Caenorhabditis elegans. Evolution. 2020;74: 2451–2464.

37. Carelli FN, Sharma G, Ahringer J. Broad Chromatin Domains: An Important Facet of Genome Regulation. BioEssays. 2017;39: 1700124. pmid:29058338

38. Klemm SL, Shipony Z, Greenleaf WJ. Chromatin accessibility and the regulatory epigenome. Nat Rev Genet. 2019;20: 207–220. pmid:30675018

39. Evans KJ, Huang N, Stempor P, Chesney MA, Down TA, Ahringer J. Stable Caenorhabditis elegans chromatin domains separate broadly expressed and developmentally regulated genes. Proc Natl Acad Sci. 2016;113: E7020–E7029. pmid:27791097

40. McMurchy AN, Stempor P, Gaarenstroom T, Wysolmerski B, Dong Y, Aussianikava D, et al. A team of heterochromatin factors collaborates with small RNA pathways to combat repetitive elements and germline stress. eLife. 2017;6: e21666. pmid:28294943

41. Ashe A, Sapetschnig A, Weick E-M, Mitchell J, Bagijn MP, Cording AC, et al. piRNAs Can Trigger a Multigenerational Epigenetic Memory in the Germline of C. elegans. Cell. 2012;150: 88–99. pmid:22738725

42. Greer EL, Maures TJ, Ucar D, Hauswirth AG, Mancini E, Lim JP, et al. Transgenerational epigenetic inheritance of longevity in Caenorhabditis elegans. Nature. 2011;479: 365–371. pmid:22012258

43. Klosin A, Casas E, Hidalgo-Carcedo C, Vavouri T, Lehner B. Transgenerational transmission of environmental information in C. elegans. Science. 2017;356: 320–323. pmid:28428426

44. Lev I, Seroussi U, Gingold H, Bril R, Anava S, Rechavi O. MET-2-Dependent H3K9 Methylation Suppresses Transgenerational Small RNA Inheritance. Curr Biol. 2017;27: 1138–1147. pmid:28343968

45. Shirayama M, Seth M, Lee H-C, Gu W, Ishidate T, Conte D, et al. piRNAs Initiate an Epigenetic Memory of Nonself RNA in the C. elegans Germline. Cell. 2012;150: 65–77. pmid:22738726

46. Woodhouse RM, Buchmann G, Hoe M, Harney DJ, Low JKK, Larance M, et al. Chromatin Modifiers SET-25 and SET-32 Are Required for Establishment but Not Long-Term Maintenance of Transgenerational Epigenetic Inheritance. Cell Rep. 2018;25: 2259–2272.e5. pmid:30463020

47. Buenrostro JD, Wu B, Chang HY, Greenleaf WJ. ATAC-seq: A Method for Assaying Chromatin Accessibility Genome-Wide. Curr Protoc Mol Biol. 2015;109. pmid:25559105

48. Daugherty AC, Yeo RW, Buenrostro JD, Greenleaf WJ, Kundaje A, Brunet A. Chromatin accessibility dynamics reveal novel functional enhancers in C. elegans. Genome Res. 2017;27: 2096–2107. pmid:29141961

49. Zhang Y-Y, Latzel V, Fischer M, Bossdorf O. Understanding the evolutionary potential of epigenetic variation: a comparison of heritable phenotypic variation in epiRILs, RILs, and natural ecotypes of Arabidopsis thaliana. Heredity. 2018;121: 257–265. pmid:29875373

50. Jänes J, Dong Y, Schoof M, Serizay J, Appert A, Cerrato C, et al. Chromatin accessibility dynamics across C. elegans development and ageing. eLife. 2018;7: e37344. pmid:30362940

51. Denver DR, Dolan PC, Wilhelm LJ, Sung W, Lucas-Lledo JI, Howe DK, et al. A genome-wide view of Caenorhabditis elegans base-substitution mutation processes. Proc Natl Acad Sci. 2009;106: 16310–16314. pmid:19805298

52. Konrad A, Brady MJ, Bergthorsson U, Katju V. Mutational Landscape of Spontaneous Base Substitutions and Small Indels in Experimental Caenorhabditis elegans Populations of Differing Size. Genetics. 2019;212: 837–854. pmid:31110155

53. Meier B, Cooke SL, Weiss J, Bailly AP, Alexandrov LB, Marshall J, et al. C. elegans whole-genome sequencing reveals mutational signatures related to carcinogens and DNA repair deficiency. Genome Res. 2014;24: 1624–1636. pmid:25030888

54. Volkova NV, Meier B, González-Huici V, Bertolini S, Gonzalez S, Vöhringer H, et al. Mutational signatures are jointly shaped by DNA damage and repair. Nat Commun. 2020;11: 2169. pmid:32358516

55. Gu W, Shirayama M, Conte D, Vasale J, Batista PJ, Claycomb JM, et al. Distinct Argonaute-Mediated 22G-RNA Pathways Direct Genome Surveillance in the C. elegans Germline. Mol Cell. 2009;36: 231–244. pmid:19800275

56. Hoogstrate SW, Volkers RJ, Sterken MG, Kammenga JE, Snoek LB. Nematode endogenous small RNA pathways. Worm. 2014;3: e28234. pmid:25340013

57. Lau AC, Nabeshima K, Csankovszki G. The C. elegans dosage compensation complex mediates interphase X chromosome compaction. Epigenetics Chromatin. 2014;7: 31. pmid:25400696

58. Lau AC, Csankovszki G. Balancing up and downregulation of the C. elegans X chromosomes. Curr Opin Genet Dev. 2015;31: 50–56. pmid:25966908

59. Meyer BJ. The X chromosome in C. elegans sex determination and dosage compensation. Curr Opin Genet Dev. 2022;74: 101912. pmid:35490475

60. Batista PJ, Ruby JG, Claycomb JM, Chiang R, Fahlgren N, Kasschau KD, et al. PRG-1 and 21U-RNAs Interact to Form the piRNA Complex Required for Fertility in C. elegans. Mol Cell. 2008;31: 67–78. pmid:18571452

61. Das PP, Bagijn MP, Goldstein LD, Woolford JR, Lehrbach NJ, Sapetschnig A, et al. Piwi and piRNAs Act Upstream of an Endogenous siRNA Pathway to Suppress Tc3 Transposon Mobility in the Caenorhabditis elegans Germline. Mol Cell. 2008;31: 79–90. pmid:18571451

62. Ruby JG, Jan C, Player C, Axtell MJ, Lee W, Nusbaum C, et al. Large-Scale Sequencing Reveals 21U-RNAs and Additional MicroRNAs and Endogenous siRNAs in C. elegans. Cell. 2006;127: 1193–1207. pmid:17174894

63. Beltran T, Barroso C, Birkle TY, Stevens L, Schwartz HT, Sternberg PW, et al. Comparative Epigenomics Reveals that RNA Polymerase II Pausing and Chromatin Domain Organization Control Nematode piRNA Biogenesis. Dev Cell. 2019;48: 793–810.e6. pmid:30713076

64. Seth M, Shirayama M, Gu W, Ishidate T, Conte D, Mello CC. The C. elegans CSR-1 Argonaute Pathway Counteracts Epigenetic Silencing to Promote Germline Gene Expression. Dev Cell. 2013;27: 656–663. pmid:24360782

65. Jones LM, Rayson SJ, Flemming AJ, Urwin PE. Adaptive and Specialised Transcriptional Responses to Xenobiotic Stress in Caenorhabditis elegans Are Regulated by Nuclear Hormone Receptors. Chatterjee B, editor. PLoS ONE. 2013;8: e69956. pmid:23922869

66. Herbette M, Mercier MG, Michal F, Cluet D, Burny C, Yvert G, et al. The C. elegans SET-2/SET1 histone H3 Lys4 (H3K4) methyltransferase preserves genome stability in the germline. DNA Repair. 2017;57: 139–150. pmid:28779964

67. Harlow PH, Perry SJ, Stevens AJ, Flemming AJ. Comparative metabolism of xenobiotic chemicals by cytochrome P450s in the nematode Caenorhabditis elegans. Sci Rep. 2018;8: 13333. pmid:30190484

68. Hartman JH, Widmayer SJ, Bergemann CM, King DE, Morton KS, Romersi RF, et al. Xenobiotic metabolism and transport in Caenorhabditis elegans. J Toxicol Environ Health Part B. 2021;24: 51–94. pmid:33616007

69. Lindblom TH, Dodd AK. Xenobiotic Detoxification in the Nematode Caenorhabditis elegans. J Exp Zool. 2006;1: 720–730. pmid:16902959

70. Kurz CL, Shapira M, Chen K, Baillie DL, Tan M-W. C. elegans pgp-5 IS INVOLVED IN RESISTANCE TO BACTERIAL. Biochem Biophys Res Commun. 2007;16: 438–443.

71. Allen MA, Hillier LW, Waterston RH, Blumenthal T. A global analysis of C. elegans trans -splicing. Genome Res. 2011;21: 255–264. pmid:21177958

72. Ing-Simmons E, Vaid R, Bing XY, Levine M, Mannervik M, Vaquerizas JM. Independence of chromatin conformation and gene regulation during Drosophila dorsoventral patterning. Nat Genet. 2021;53: 487–499. pmid:33795866

73. Kiani K, Sanford EM, Goyal Y, Raj A. Changes in chromatin accessibility are not concordant with transcriptional changes for single-factor perturbations. Bioinformatics; 2022 Feb.

74. Cannavò E, Khoueiry P, Garfield DA, Geeleher P, Zichner T, Gustafson EH, et al. Shadow Enhancers Are Pervasive Features of Developmental Regulatory Networks. Curr Biol. 2016;26: 38–51. pmid:26687625

75. Kvon EZ, Waymack R, Gad M, Wunderlich Z. Enhancer redundancy in development and disease. Nat Rev Genet. 2021;22: 324–336. pmid:33442000

76. Zhang T, Cooper S, Brockdorff N. The interplay of histone modifications–writers that read. EMBO Rep. 2015;16: 1467–1481. pmid:26474904

77. Simpson VJ, Johnson TE, Hammen RF. Caenorhabds elegans DNA does not contain 5-methylcytosine at any time during development or aging. Nucleic Acids Res. 1986;14: 6711–6719.

78. Greer EL, Blanco MA, Gu L, Sendinc E, Liu J, Aristizábal-Corrales D, et al. DNA Methylation on N6-Adenine in C. elegans. Cell. 2015;161: 868–878. pmid:25936839

79. Johnson CL, Spence AM. Epigenetic Licensing of Germline Gene Expression by Maternal RNA in C. elegans. Science. 2011;333: 1311–1314. pmid:21885785

80. Zhao Z, Sheps JA, Ling V, Fang LL, Baillie DL. Expression Analysis of ABC Transporters Reveals Differential Functions of Tandemly Duplicated Genes in Caenorhabditis elegans. J Mol Biol. 2004;344: 409–417. pmid:15522294

81. Tao J, Hao Y, Li X, Yin H, Nie X, Zhang J, et al. Systematic Identification of Housekeeping Genes Possibly Used as References in Caenorhabditis elegans by Large-Scale Data Integration. Cells. 2020;9: 786. pmid:32213971

82. Angers B, Castonguay E, Massicotte R. Environmentally induced phenotypes and DNA methylation: how to deal with unpredictable conditions until the next generation and after. Mol Ecol. 2010;19: 1283–1295. pmid:20298470

83. Leung C, Angers B, Bergeron P. Epigenetic anticipation for food and reproduction. Faulk C, editor. Environ Epigenetics. 2020;6: dvz026. pmid:32015901

84. Danchin E, Pocheville A, Rey O, Pujol B, Blanchet S. Epigenetically facilitated mutational assimilation: epigenetics as a hub within the inclusive evolutionary synthesis: Epigenetics as a hub for genetic assimilation. Biol Rev. 2019;94: 259–282.

85. Green MR, Sambrook J. Precipitation of RNA with Ethanol. Cold Spring Harb Protoc. 2020;2020: pdb.prot101717. pmid:32123016

86. Almeida MV, de Jesus Domingues AM, Lukas H, Mendez-Lago M, Ketting RF. RppH can faithfully replace TAP to allow cloning of 5′-triphosphate carrying small RNAs. MethodsX. 2019;6: 265–272. pmid:30788220

87. Langmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nat Methods. 2012;9: 357–359. pmid:22388286

88. Kent WJ, Sugnet CW, Furey TS, Roskin KM, Pringle TH, Zahler AM, et al. The Human Genome Browser at UCSC. Genome Res. 2002;12: 996–1006. pmid:12045153

89. Quinlan AR, Hall IM. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics. 2010;26: 841–842. pmid:20110278

90. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15: 550. pmid:25516281

91. Langmead B, Trapnell C, Pop M, Salzberg SL. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 2009;10: R25. pmid:19261174

92. Kozomara A, Birgaoanu M, Griffiths-Jones S. miRBase: from microRNA sequences to function. Nucleic Acids Res. 2019;47: D155–D162. pmid:30423142

93. Bates D, Mächler M, Bolker B, Walker S. Fitting Linear Mixed-Effects Models Using lme4. J Stat Softw. 2015;67.

94. Hwang SB, Lee J. Neuron Cell Type-specific SNAP-25 Expression Driven by Multiple Regulatory Elements in the Nematode Caenorhabditis elegans. J Mol Biol. 2003;333: 237–247. pmid:14529613

95. Landmann F, Quintin S, Labouesse M. Multiple regulatory elements with spatially and temporally distinct activities control the expression of the epithelial differentiation gene lin-26 in C. elegans. Dev Biol. 2004;265: 478–490. pmid:14732406

96. Zhao Z, Fang L, Chen N, Johnsen RC, Stein L, Baillie DL. Distinct Regulatory Elements Mediate Similar Expression Patterns in the Excretory Cell of Caenorhabditis elegans. J Biol Chem. 2005;280: 38787–38794. pmid:16159881

97. Yang Y-L, Loh K-S, Liou B-Y, Chu I-H, Kuo C-J, Chen H-D, et al. SESN-1 is a positive regulator of lifespan in Caenorhabditis elegans. Exp Gerontol. 2013;48: 371–379. pmid:23318476

98. Chen EY, Tan CM, Kou Y, Duan Q, Wang Z, Meirelles GV, et al. Enrichr: interactive and collaborative HTML5 gene list enrichment analysis tool. BMC Bioinformatics. 2013;14: 128. pmid:23586463

99. Kuleshov MV, Jones MR, Rouillard AD, Fernandez NF, Duan Q, Wang Z, et al. Enrichr: a comprehensive gene set enrichment analysis web server 2016 update. Nucleic Acids Res. 2016;44: W90–W97. pmid:27141961