1. Plant Biology
Download icon

Identification and characterisation of hypomethylated DNA loci controlling quantitative resistance in Arabidopsis

  1. Leonardo Furci
  2. Ritushree Jain
  3. Joost Stassen
  4. Oliver Berkowitz
  5. James Whelan
  6. David Roquis
  7. Victoire Baillet
  8. Vincent Colot
  9. Frank Johannes
  10. Jurriaan Ton  Is a corresponding author
  1. University of Sheffield, United Kingdom
  2. La Trobe University, Australia
  3. Technical University of Munich, Germany
  4. Institut de Biologie de l’Ecole Normale Supérieure (IBENS), Ecole Normale Supérieure, Centre National de la Recherche Scientifique (CNRS), Institut National de la Santé et de la Recherche Médicale (INSERM), PSL Université Paris, France
Research Article
  • Cited 0
  • Views 1,536
  • Annotations
Cite this article as: eLife 2019;8:e40655 doi: 10.7554/eLife.40655

Abstract

Variation in DNA methylation enables plants to inherit traits independently of changes to DNA sequence. Here, we have screened an Arabidopsis population of epigenetic recombinant inbred lines (epiRILs) for resistance against Hyaloperonospora arabidopsidis (Hpa). These lines share the same genetic background, but show variation in heritable patterns of DNA methylation. We identified four epigenetic quantitative trait loci (epiQTLs) that provide quantitative resistance without reducing plant growth or resistance to other (a)biotic stresses. Phenotypic characterisation and RNA-sequencing analysis revealed that Hpa-resistant epiRILs are primed to activate defence responses at the relatively early stages of infection. Collectively, our results show that hypomethylation at selected pericentromeric regions is sufficient to provide quantitative disease resistance, which is associated with genome-wide priming of defence-related genes. Based on comparisons of global gene expression and DNA methylation between the wild-type and resistant epiRILs, we discuss mechanisms by which the pericentromeric epiQTLs could regulate the defence-related transcriptome.

https://doi.org/10.7554/eLife.40655.001

eLife digest

In plants, animals and microbes genetic information is encoded by DNA, which are made up of sequences of building blocks, called nucleotide bases. These sequences can be separated into sections known as genes that each encode specific traits. It was previously thought that only changes to the sequence of bases in a DNA molecule could alter the traits passed on to future generations. However, it has recently become clear that some traits can also be inherited through modifications to the DNA that do not alter its sequence.

One such modification is to attach a tag, known as a methyl group, to a nucleotide base known as cytosine. These methyl tags can be added to, or removed from, DNA to create different patterns of methylation. Previous studies have shown that plants whose DNA is less methylated than normal (‘hypo-methylated’) are more resistant to plant diseases. However, the location and identity of the hypo-methylated DNA regions controlling this resistance remained unknown. To address this problem, Furci, Jain et al. studied how DNA methylation in a small weed known as Arabidopsis thaliana affects how well the plants can resist a disease known as downy mildew.

Furci, Jain et al. studied a population of over 100 A. thaliana lines that have the same DNA sequences but different patterns of DNA methylation. The experiments identified four DNA locations that were less methylated in lines with enhanced resistance to downy mildew. Importantly, this form of resistance did not appear to reduce how well the plants grew, or make them less able to resist other diseases or environmental stresses.

The results of further experiments suggested that reduced methylation at the four DNA regions prime the plant’s immune system, enabling a faster and stronger activation of a multitude of defence genes across the genome after attack by downy mildew.

The next steps following on from this work are to investigate exactly how the four DNA regions with reduced methylation can prime so many different defence genes in the plant. Further research is also needed to determine whether it is possible to breed crop plants with lower levels of methylation at specific DNA locations to improve disease resistance, but without decreasing the amount and quality of food produced.

https://doi.org/10.7554/eLife.40655.002

Introduction

Eukaryotic cytosine methylation plays an important role in the regulation of gene expression and genome stability. In plants, this form of DNA methylation occurs at three sequence contexts: CG, CHG and CHH, where H indicates any base except guanine (G) (Vanyushin, 2006; Law and Jacobsen, 2010). Patterns of plant DNA methylation in the plant genome can remain stable over multiple generations and influence heritable phenotypes (Quadrana and Colot, 2016). Recent evidence has suggested that reduced DNA methylation increases the responsiveness of the plant immune system (Espinas et al., 2016) This ‘priming’ of plant defence enables an augmented induction of defence-related genes after pathogen attack, causing increased levels of quantitative resistance (Prime-A-Plant Group et al., 2006; Conrath et al., 2015; Martinez-Medina et al., 2016; Liégard et al., 2018). In some cases, priming of defence-related genes is associated with post-translational histone modifications that mark a more open chromatin structure (Jaskiewicz et al., 2011; Luna et al., 2012). Additional evidence for epigenetic regulation of plant immunity has come from independent studies reporting that disease-exposed Arabidopsis produces progeny that expresses transgenerational acquired resistance (TAR), which is associated with priming of defence-related genes (Luna et al., 2012; Slaughter et al., 2012). Furthermore, Arabidopsis mutants that are impaired in the establishment or maintenance of DNA methylation mimic TAR-related priming without prior priming stimulus (Lopez et al., 2011; Luna and Ton, 2012; López Sánchez et al., 2016). By contrast, the hyper-methylated ros1-4 mutant, which is impaired in active DNA de-methytation, is more susceptible to biotrophic pathogens, affected in defence gene responsiveness, and impaired in TAR (López Sánchez et al., 2016; Yu et al., 2013). Thus, DNA (de)methylation determines quantitative disease resistance by influencing the responsiveness of defence-related genes. However, causal evidence that selected hypomethylated DNA loci are responsible for the meiotic transmission of this form of quantitative disease resistance is lacking.

Epigenetic Recombinant Inbred Lines (epiRILs) have been developed with the aim to study the epigenetic basis of heritable plant traits (Reinders et al., 2009; Johannes et al., 2009). EpiRILs show little differences in DNA sequence, but vary substantially in DNA methylation. A commonly used population of epiRILs is derived from a cross between the Arabidopsis wild-type (Wt) accession Col-0 and the decreased DNA methylation1-2 (ddm1-2) mutant (Johannes et al., 2009). The DDM1 protein is a chromatin re-modelling enzyme that provides DNA methyltransferase enzymes access to heterochromatic transposable elements (TEs) (Jeddeloh et al., 1998; Brzeski and Jerzmanowski, 2003; Zemach et al., 2013). Accordingly, the ddm1-2 mutation causes loss of pericentromeric heterochromatin and reduced DNA methylation in all sequence contexts (Kakutani et al., 1996; Ito et al., 2015). Although the epiRILs from the ddm1-2 x Col-0 cross do not carry the ddm1-2 mutation, they contain stably inherited hypomethylated DNA regions from the ddm1-2 parent, which are maintained up to 16 generations of self-pollination (Johannes et al., 2009; Colomé-Tatché et al., 2012; Latzel et al., 2013). A core set of 123 epiRILs from this population at the eight generation of self-pollination in the wild-type (Wt) background has been characterized for differentially methylated region (DMR) markers, enabling linkage mapping of heritable hypomethylated loci controlling root growth, flowering and abiotic stress tolerance (Liégard et al., 2018; Cortijo et al., 2014; Kooke et al., 2015).

In this study, we have characterised the core set of 123 lines from the ddm1-2 x Col-0 epiRIL population for resistance against the biotrophic downy mildew pathogen Hyaloperonospora arabidopsidis (Hpa) to search for heritable hypomethylated loci controlling disease resistance. We identified four of these epigenetic quantitative trait loci (epiQTLs), accounting for 60% of the variation in disease resistance. None of these epiQTLs were associated with growth impairment, indicating that the resistance does not incur major physiological costs on plant development. Further phenotypic characterisation and transcriptome analysis of selected Hpa-resistant epiRILs revealed that their resistance is associated with genome-wide priming of defence-related genes. Interestingly, bisulfite sequencing did not reveal defence regulatory genes inside the epiQTL regions that were simultaneously primed and hypomethylated, suggesting that DDM1-dependent DNA methylation at the epiQTLs trans-regulates the responsiveness of distant defence genes.

Results

Identification of epiQTLs controlling quantitative resistance against Hpa

To examine the role of DDM1-dependent DNA methylation in heritable disease resistance, 123 epiRILs from the ddm1-2 x Col-0 cross were analysed for Hpa resistance and compared to siblings of the ddm1-2 parent (Figure 1a, red), the Wt parent (Col-0), and five progenies thereof (Figure 1a, green). Leaves of 3-week-old plants were inoculated with Hpa conidiospores and then collected for trypan-blue staining at 6 days post-inoculation (dpi). Microscopic classification of leaves into four classes of Hpa colonisation (Figure 1—figure supplement 1) revealed 51 epiRILs with statistically enhanced levels of resistance compared to each susceptible Wt line (Pearson’s Chi-squared tests, p<0.05). Of these, eight epiRILs showed similar levels of Hpa resistance as the ddm1-2 line (Figure 1a, dark blue triangles; Pearson’s Chi-squared test, p>0.05), whereas 43 epiRILs showed intermediate levels of resistance. To identify the epiQTL(s) responsible for the observed variation in Hpa resistance, the categorical classification of Hpa infection was converted into a single value numerical resistance index (RI; Figure 1a, bottom graph). Using a linkage map of stably inherited DMR markers (Colomé-Tatché et al., 2012) (Supplementary file 1 dataset S1), interval mapping revealed four statistically significant epiQTLs on chromosomes I, II, IV and V (Figure 1b). The epiQTL on chromosome II had the highest logarithm of odds (LOD) value. For all epiQTLs, the DMR markers with the highest LOD scores (‘peak markers’) showed a positive correlation between ddm1-2 haplotype and RI (Figure 1c), indicating that the hypomethylated haplotype from ddm1-2 increases resistance against Hpa. A linear regression model to calculate the percentage of RI variance explained by each peak marker (R (2)(g)) (Cortijo et al., 2014) confirmed that the DMR peak marker of the epiQTL on chromosome II had the strongest contribution to RI variation. Using an additive model, the combined contribution of all epiQTL peak markers to RI variation (R2(G)) (Cortijo et al., 2014) was estimated at 60.0% (Figure 1d).

Figure 1 with 5 supplements see all
Mapping of epigenetic quantitative trait loci (epiQTL) controlling transgenerational resistance against Hyaloperonosopra arabidopsidis (Hpa).

(a) Levels of Hpa resistance in 123 epiRIL lines, the ddm1-2 line (F4; red triangle) and six Wt lines (Col-0; green triangles). Top graph shows distribution of infection classes in each epiRIL; blue triangles pinpoint the eight most resistant epiRILs with statistically similar levels of Hpa colonisation as the ddm1-2 line (Pearson’s Chi-squared test, p>0.05). Bottom graph shows variation in Hpa resistance index (RI). Green bars: Wt lines; red bar: ddm1-2; blue bars eight most resistant epiRILs (n > 100). (b) Linkage analysis of RI (blue line) and green leaf area (GLA) of three-week-old seedlings (green). Green bars at the bottom represent chromosomes. Red line represents the threshold of significance. Peak DMR markers with the highest LOD scores are shown on top. (c) Correlation plots between peak marker haplotype (methylated Wt versus hypomethylated ddm1-2) and RI (blue) or GLA (green). (d) Percentages of resistance variance explained by the peak DMR markers, including covariance between markers (orange).

https://doi.org/10.7554/eLife.40655.003

DNA methylation maintains genome stability by preventing transposition of TEs. In the Col-0 x ddm1-2 epiRIL population, reduced methylation at the ddm1-2 haplotype occurs predominantly at long transposons in heterochromatic pericentromeric regions (Zemach et al., 2013; Colomé-Tatché et al., 2012). Frequent transposition events in the epiRILs are nevertheless rare as most DNA hypomethylation occurs at relic transposons that have lost the ability to transpose, and the occurrence of independent transposition events at similar loci is extremely unlikely (Colomé-Tatché et al., 2012; Matzke and Mosher, 2014). However, it is possible that transposition events originating from the heavily hypomethylated ddm1-2 parent were crossed into the population, resulting into shared transposition events (STEs) between multiple epiRILs, which could have contributed to variation in resistance. To account for this possibility, we compared the genomic DNA sequences of the four epiQTL intervals from 122 epiRILs (LOD drop-off = 2) for the presence of STEs in more than two epiRILs, using TE-tracker software (Gilly et al., 2014). This analysis revealed three STEs in the epiQTL interval on chromosome I (Supplementary file 1 dataset S2), while no STEs could be detected in the other epiQTL intervals. None of the STEs in the epiQTL on chromosome I showed statistically significant linkage with RI (Supplementary file 1 dataset S2). Accordingly, we conclude that the segregating Hpa resistance in the epiRIL population is caused by epigenetic variation in DNA methylation, rather than genetic variation by STEs.

Effects of the resistance epiQTLs on plant growth and resistance against other (a) biotic stresses

Expression of inducible defence mechanisms is often associated with physiological costs, resulting in reduced plant growth (Huot et al., 2014). To determine whether the resistance that is controlled by the four epiQTLs is associated with costs to plant growth, we quantified the green leaf area (GLA) of 12–15 individual plants per line at the stage of Hpa inoculation (Figure 1—figure supplement 2). Subsequent interval mapping revealed one statistically significant epiQTL on chromosome I (Figure 1b). The corresponding peak marker (MM150) showed a negative correlation between GLA and ddm1-2 haplotype (Figure 1c), indicating that the hypomethylated ddm1-2 allele at this locus represses plant growth. The growth epiQTL mapped to a different region than the resistance epiQTL on chromosome I (Figure 1b, inset). Furthermore, none of the eight most resistant epiRILs showed significant growth reduction compared to all Wt lines in the screen (Figure 1—figure supplement 2). Hence, the resistance provided by the four hypomethylated epiQTLs is not associated with major physiological costs to plant growth.

Enhanced defence to one stress can lead to enhanced susceptibility to another stress, which is caused by antagonistic cross-talk between defence signalling pathways (Koornneef and Pieterse, 2008). To examine whether Hpa resistance in the epiRIL population is associated with increased susceptibility to other stresses, we compared the eight most Hpa-resistant epiRILs (Figure 1a; Figure 1—figure supplement 3a) for resistance against the necrotrophic fungus Plectosphaerella cucumerina (Pc) and tolerance to salt (NaCl). At 9 dpi with Pc spores, epiRIL#193 showed a statistically significant reduction in necrotic lesion size compared to the Wt (line #602), indicating enhanced resistance (Figure 1—figure supplement 3b). The seven other epiRILs showed unaffected levels of Pc resistance that were similar to the Wt. Salt tolerance was quantified by the percentage of seedlings with fully developed cotyledons at 6 days after germination on agar medium with increasing NaCl concentrations. Remarkably, all Hpa-resistant epiRILs showed varying degrees of tolerance to the highest NaCl concentration compared to Wt plants (Figure 1—figure supplement 3c). Thus, the quantitative resistance to Hpa in the epiRIL population does not compromise resistance against necrotrophic pathogens or abiotic stress.

Hpa-resistant epiRILs are primed to activate different defence mechanisms

Basal resistance against Hpa involves a combination of salicylic acid (SA)-dependent and SA-independent defence mechanisms (Knoth et al., 2007; Coates and Beynon, 2010). To examine the role of SA-dependent defences, we profiled the expression of the SA-inducible marker gene PR1 at 48 and 72 hr post-inoculation (hpi), which represents a critical time-window for host defence against Hpa (Koch and Slusarenko, 1990; Soylu and Soylu, 2003). None of the epiRILs showed a statistically significant increase in basal PR1 expression after mock inoculation (Figure 2a; Figure 1—figure supplement 4a), indicating that the resistance is not based on constitutive up-regulation of SA-dependent defence signalling. However, in comparison to the Wt line, epiRILs #71, #148, #193, #229 and #508 showed augmented induction of PR1 at 48 and/or 72 hpi with Hpa (Figure 2a; Figure 1—figure supplement 4a), indicating priming of SA-inducible defences (Martinez-Medina et al., 2016). To assess the role of cell wall defence, all lines were analysed for effectiveness of callose deposition, which is a pathogen-inducible defence mechanism that is largely controlled by SA-independent signalling (Luna et al., 2011). Compared to the Wt line, all but one epiRIL (#193) showed a statistically significant increase in the proportion of callose-arrested germ tubes (Figure 2a; Figure 1—figure supplement 4b). Hence, the eight most Hpa-resistant epiRILs are primed to activate differentially regulated defence responses, which explains the lack of major costs on growth and compatibility with other types of (a)biotic stress resistance in the epiRILs (Figures 1b and 2a; Figure 1—figure supplements 24).

Figure 2 with 4 supplements see all
The defence-related transcriptome of Hpa-resistant epiRILs.

(a) Defence marker phenotypes and epiQTL haplotypes of 4 Hpa-resistant epiRILs and the Wt (#602), which were analysed by RNA sequencing. Top graph: relative expression of SA-dependent PR1 at 72 hr after inoculation (hpi) with Hpa (red) or water (blue). Middle graph: resistance efficiency of callose deposition in Hpa-inoculated plants. Shown are percentages of arrested (light) and non-arrested (dark) germ tubes at 48 hpi. Bottom panel: epiQTL haplotypes of selected lines. Green: methylated Wt haplotype; yellow: hypomethylated ddm1-2 haplotype. Asterisks indicate statistically significant differences to the Wt. (see Figure 1—figure supplement 4 for statistical information). (b) Principal component analysis of 27,641 genes at 48 (small symbols) and 72 (large symbols) hpi with Hpa (triangles) or water (Mock; circles). Colours indicate different lines. (c) Numbers and expression profiles of Hpa-inducible genes that show constitutively enhanced expression (Group 1) or augmented levels of Hpa-induced expression (Group 2) in the Hpa-resistant epiRILs at 48 or 72 hpi. Heatmaps show normalised standard deviations from the mean (z-scores) for each gene (rows), using rlog-transformed read counts (see Figure 2—figure supplements 3 and 4 for better detail) (d) GO term enrichment of primed and constitutively up-regulated genes. Shown are 469 GO terms (rows), for which one or more epiRIL(s) displayed a statistically significant enrichment in one or more categories (Hypergeometric test, followed by Benjamini-Hochberg FDR correction; q < 0.05). Heatmap-projected values for each GO term (rows) represent percentage of GO-annotated genes in each category relative to all GO-annotated genes in the Arabidopsis genome (TAIR10). Black bar on the top right indicates 111 defence-related GO terms.

https://doi.org/10.7554/eLife.40655.014

Transgenerational stability of the resistance

The 123 epiRILs analysed for Hpa resistance had been self-pollinated for eight generations in a Wt (Col-0) genetic background since the F1 x Col-0 backcross (F9) (Johannes et al., 2009). To examine the transgenerational stability of the resistance phenotype over one more generation, five individuals from the eight most resistant epiRILs and the Wt line (Figure 1a, Figure 1—figure supplement 3a) were selected to generate F10 families, which were then tested for Hpa resistance. Comparing distributions of pooled leaves from all five families per line confirmed that each epiRIL maintained a statistically enhanced level of resistance (Figure 1—figure supplement 5; Pearson’s Chi-squared test, p<0.05; top asterisks). However, when comparing individual F10 families to the Wt, 2 of the 40 F10 families (line #71–2 and line #148–2) exhibited Wt levels of susceptibility, indicating that they had lost Hpa resistance from the F9 to the F10 generation. Furthermore, four of the eight epiRILs tested (#71, #148, #545, and #508) displayed statistically significant variation in Hpa resistance between the 5 F10 families within the epiRIL (Figure 1—figure supplement 5; Pearson’s Chi-squared test, p<0.05; † symbols), suggesting instability of the Hpa resistance.

Hpa-resistant epiRILs show genome-wide priming of defence-related genes

To study the transcriptomic basis of the transgenerational resistance, Wt plants (line #602) and 4 Hpa-resistant epiRILs (#148, #193, #454 and #508), each carrying different combinations of the four epiQTLs, were analysed by RNA sequencing at 48 and 72 hpi (Figure 2a, bottom panel). Principal component analysis (PCA) of biologically replicated samples (n = 3) revealed clear separation between all treatment/time-point/epi-genotype combinations (Figure 2b). The first PCA axis explained 31% of the variation in transcript abundance, separating samples from mock- and Hpa-treated plants, whereas the second PCA axis explained 20% of the variation, mostly separating samples from the different lines (Figure 2b). This PCA pattern indicates that the response to Hpa infection had a bigger effect on global gene expression than epi-genotype. Moreover, samples from Hpa-inoculated epiRILs showed relatively little difference between both time-points (Figure 2b), whereas samples from Hpa-inoculated Wt plants at 48 hpi clustered between samples from mock-inoculated Wt plants and samples from Hpa-inoculated Wt plants at 72 hpi. This pattern suggests a difference in the speed and/or intensity of the transcriptional response to Hpa. To explore this possibility further, we performed three-factorial likelihood ratio tests (q < 0.05) to select differentially expressed genes between all epigenotype/treatment/time-point combinations. This analysis identified 20,569 genes, representing 61% of all annotated RNA-producing genes in the Arabidopsis genome, including transposable elements, non-coding RNA genes and pseudogenes (Supplementary file 1 dataset S3). Of these, 9364 genes were induced by Hpa at 48 and/or 72 hpi in one or more lines (Supplementary file 1 dataset S4). Subsequent hierarchical clustering of this gene selection revealed a large cluster of Hpa-inducible transcripts displaying augmented induction in the epiRILs at the relatively early time-point of 48 hr after Hpa inoculation (Figure 2—figure supplement 1).

To characterize further the pathogen-inducible transcriptome of the resistant epiRILs, we selected Hpa-inducible genes showing elevated levels of expression in the epiRILs during Hpa infection. Within this gene selection, we distinguished two expression profiles. The first group of genes had been selected for constitutively enhanced expression in the resistant epiRILs, using the following criteria (Wald tests, q < 0.05): (i) Hpa-inducible in the Wt, (ii) not inducible by Hpa in the epiRIL and (iii) displaying enhanced accumulation in mock-treated epiRIL that is equal or higher than accumulation in the Hpa-inoculated Wt (‘Group 1’; Figure 2—figure supplement 2a). The second group of genes had been selected for enhanced Hpa-induced expression in the epiRILs, using the following criteria (Wald tests, q < 0.05): (i) Hpa-inducible in the Wt (#602), (ii) Hpa-inducible in the epiRIL(s) and (iii) displaying statistically increased accumulation in Hpa-inoculated epiRILs compared to Hpa-inoculated Wt plants (‘Group 2’; Figure 2—figure supplement 2a). For each epiRIL, we identified more genes in Group 2 than in Group 1 (Figure 2c; Figure 2—figure supplements 2b, 3 and 4; Supplementary file 1 datasets S5 and S6). This difference was most pronounced at 48 hpi, which represents a critical time-point for host defence against Hpa (Koch and Slusarenko, 1990; Soylu and Soylu, 2003). Analysis of a statistical interaction between epi-genotype x Hpa treatment revealed that >92% of all genes in Group 2 are significant for this interaction term (Supplementary file 1 dataset S7), indicating a constitutively primed expression pattern. Visualisation of the expression profiles in heatmaps confirmed this notion, showing that the induction of Group 2 genes by Hpa is strongly augmented in the resistant epiRILs compared to the Wt line (Figure 2c; Figure 2—figure supplement 4), which is consistent with the definition of plant defence priming (Martinez-Medina et al., 2016).

To examine the functional contributions of the Hpa-inducible genes in Groups 1 and 2, we employed gene ontology (GO) term enrichment analysis. After exclusion of redundant GO terms (Jantzen et al., 2011), we identified 469 GO terms, for which one or more of the sets showed statistically significant enrichment. Group 2 genes at 48 hpi displayed dramatically enhanced GO term enrichment compared to all other sets, which was obvious for all epiRILs (Figure 2d). This enrichment was particularly pronounced for 111 GO terms relating to SA-dependent and SA-independent defence mechanisms (Supplementary file 1 dataset S8), which supports our phenotypic characterisation of SA-dependent and SA-independent defence markers (Figure 1—figure supplement 4). Collectively, these results suggest that the quantitative resistance of the epiRILs is based on priming of Hpa-inducible defence genes.

Interestingly, compared to the other gene selections, a relatively large proportion of defence-related genes in Group 2 at 48 hpi was shared between all four epiRILs (Figure 2—figure supplement 2b), pointing to relatively high similarity in the augmented immune response of the epiRILs. Furthermore, only 5% of the genes in the Group 1% and 6.5% of the genes in Group 2 are physically located within the borders of the epiQTL intervals (LOD drop-off = 2). The frequency of Group 1 and 2 genes relative to all other genes was significantly lower for the epiQTL regions compared to the entire Arabidopsis genome (14.6%; Pearson’s Chi-squared test, p<0.05). Thus, the majority of Hpa-inducible Group 1 and 2 genes showing enhanced expression in the more resistant epiRILs are (trans-)regulated by DNA methylation at the four epiQTLs.

The resistance epiQTLs do not contain defence genes that are cis-regulated by DNA methylation, suggesting involvement of trans-regulatory mechanisms

Although 92% of all genes in Group 2 were located outside the physical borders of the four epiQTL intervals (LOD-drop-off=2), we hypothesized that a small set of defence regulatory genes inside the epiQTL regions are directly (cis-)regulated by DNA methylation to mediate augmented levels of defence in response to Hpa infection. Since the Group 2 genes were strongly enriched with defence-related GO terms (Figure 2d), we examined whether their augmented expression during Hpa infection is associated with the hypomethylated ddm1-2 haplotype. To this end, we calculated for each gene in Group 2 the ratio of normalized transcript abundance between Hpa-inoculated epiRIL and the Hpa-inoculated Wt line, which is proportional to their level of augmented expression during Hpa infection. Hierarchical clustering of these ratios enabled us to select for genes that exclusively show augmented expression when associated with the hypomethylated ddm1-2 haplotype of the corresponding epiQTL (Figure 3a; Figure 3—figure supplement 1a). The expression ratios of 279 epiQTL-localised genes did not correlate with the ddm1-2 haplotype (Figure 3a, cluster II; Figure 3—figure supplement 1a; Supplementary file 1 dataset S9), indicating that DNA methylation does not cis-regulate their augmented Hpa-inducible expression. By contrast, 73 epiQTL-localised genes only showed augmented expression when associated with the hypomethylated ddm1-2 haplotype (Figure 3a, cluster I; Figure 3—figure supplement 1a; Supplementary file 1 dataset S10). To confirm the hypomethylated status of these genes, we performed comprehensive bisulfite sequencing analysis of DNA methylation for the four epiRILs and the Wt line. DMR analysis of the gene body (GB), 2 kb promoter region (P) and 1 kb downstream (D) regions confirmed that the levels of augmented gene expression of the 279 genes in cluster II do not correlate positively with the extent of DNA hypomethylation (Figure 3b, Figure 3—figure supplement 1b). This notion was confirmed by linear regression analysis between the augmented expression ratio (48 hpi) and the average level of DNA hypomethylation (Figure 3—figure supplement 2), indicating that the 279 genes in cluster II are regulated indirectly (in trans) by DNA methylation. By contrast, the 73 epiQTL-based genes in cluster I showed a positive correlation between augmented expression ratio (48 hpi) and DNA hypomethylation, which was statistically significant for each epiQTL (p<0.05; Figure 3—figure supplement 2). These results indicate that the 73 genes in cluster I are regulated locally (in cis) by DNA methylation.

Figure 3 with 5 supplements see all
Relationship between augmentation of pathogen-induced expression and DNA methylation for epiQTL-localised genes.

(a) Expression profiles of epiQTL-based genes showing elevated levels of Hpa-induced expression in one or more epiRIL(s) (Group 2). Shown are genes located in the epiQTL interval of chromosome II (epiQTL2; LOD drop-off = 2; see Figure 3—figure supplement 1a for other the epiQTLs). Heatmap shows gene expression ratios between Hpa-inoculated epiRILs and the Wt, representing augmented expression levels during pathogen attack. Hierarchical clustering yielded two distinctly regulated gene clusters (I and II). Coloured bars on the top indicate epiQTL2 haplotypes. Green: methylated Wt haplotype. Yellow: hypomethylated ddm1-2 haplotype. (b) Levels of CG DNA methylation of the same genes in the epiQTL2 interval (see Figure 3—figure supplement 1b for other epiQTLs). Heatmap shows percentages of hypomethylation (blue) or hyper-methylation (brown) relative to the Wt for 2 kb promoter regions (P), gene bodies (GB) and 1 kb downstream regions (D). (c) Distribution of gene annotations of distinctly regulated gene clusters for each epiQTL.

https://doi.org/10.7554/eLife.40655.020

Nearly all cis-regulated genes in cluster I showed a TE-like pattern of DNA methylation in the Wt (teM; methylation at CG, CHG and CHH contexts), whereas most cluster II genes showed either no methylation or a pattern of gene-body methylation in the Wt (gbM; methylation at CG only; Figure 3b and Figure 3—figure supplement 1b). Furthermore, dividing hypomethylation at gene bodies of Group 2 genes by type of DNA methylation (i.e. either teM or gbM) and plotting these values against augmented expression ratio revealed a statistically significant correlation between expression ratio and reduced teM (p=1,06e−8; Figure 3—figure supplement 3), whereas no such correlation was found for reduced gbM (p=0.66; Figure 3—figure supplement 3). These results support the growing notion that reduced teM increases gene expression, whereas changes in gbM have no direct influence on gene expression (Bewick et al., 2017).

The majority of in cis-regulated genes in cluster I genes were annotated as TEs, such as DNA transposons of the CACTA family, retrotransposons of the GYPSY or COPIA families, or TE-related genes, encoding transposases or enzymes necessary for TE function (Supplementary file 1 dataset S10). Only six genes were annotated as protein-coding genes, of which two shared homology to known protein-encoding genes (At2G07240, cysteine-type peptidase; At2G07750, RNA helicase). However, none of these two genes has previously been associated with plant defence. Furthermore, analysis of the genomic context of the six protein-coding genes revealed the presence of overlapping and/or nearby TEs (Figure 3—figure supplement 4), suggesting that their correlation between augmented expression and DNA hypomethylation is determined by association with TEs. Since TE-encoded proteins have no antimicrobial activity or direct defence regulatory function, our results suggest that global defence gene priming by hypomethylated epiQTLs is not based on cis-regulation of defence regulatory genes, but rather on alternative trans-acting mechanisms by DNA methylation of the TE-rich epiQTL.

Discussion

By screening the Col-0 x ddm1-2 epiRIL population for leaf colonisation by the downy mildew pathogen Hpa, we have identified four epiQTLs that provide quantitative disease resistance (Figure 1b). The combined contribution of all 4 DMR peak markers was estimated at 60% of the total variation (Figure 1d), which is higher than previously reported variation in developmental plant traits for this population (Latzel et al., 2013; Cortijo et al., 2014; Kooke et al., 2015). It was previously shown that half of all stably inherited DMRs in the Col-0 x ddm1-2 epiRILs also occur in natural Arabidopsis accessions (Latzel et al., 2013; Schmitz et al., 2013). Considering that the epiRIL population includes heritable variation in a range of ecologically important plant traits, including flowering, root growth, nutrient plasticity and (a)biotic stress resistance (Latzel et al., 2013; Cortijo et al., 2014; Kooke et al., 2015), it is tempting to speculate that variation in DDM1-dependent DNA methylation contributes to natural variation and environmental adaptation of Arabidopsis. Indeed, the phenotypic diversity in the Col-0 x ddm1-2 epiRIL population closely resembles that of natural Arabidopsis accessions (Roux et al., 2011; Mauch-Mani et al., 2017). Furthermore, independent studies have shown that high levels of enduring (a)biotic stress can trigger transgenerational acquired resistance (TAR) in Arabidopsis (Luna et al., 2012; Wibowo et al., 2016; Rasmann et al., 2012). Interestingly, repeated inoculation of 2- to 5 weeks old Arabidopsis seedlings with the hemi-biotrophic leaf pathogen Pseudomonas syringae pv. tomato causes TAR, which is associated with reduced transcription of DDM1 gene in local leaves that is maintained in the apical meristem of paternal plants (Furci and Ton, unpublished results). To what extent this this prolonged repression in DDM1 gene transcription causes heritable reduction in DNA methylation at the epiQTLs requires further study.

Aller et al. (2018) have recently used the same Col-0 x ddm1-2 epiRIL population to map the contribution of heritable variation in DNA methylation to the production of defence-related glucosinolate metabolites (Aller et al., 2018). Interestingly, the resistance epiQTL on chromosome I from our study partially overlaps with an epiQTL that influences basal production of the aliphatic glucosinolate 3-methylthiopropyl (3MTP) (Aller et al., 2018). Glucosinolates contribute to defence against both herbivores and microbes (Ishida et al., 2014). Moreover, myrosinase-dependent breakdown products of indole-derived 4-methoxy-indol-3-ylmethylglucosinolate have been linked to the regulation of callose-mediated cell wall defence in Arabidopsis (Clay et al., 2009; Bednarek et al., 2009). However, the 3MTP-controlling epiQTL identified by Aller et al. (2018) was relatively weak compared to the epiQTL controlling Hpa resistance (Figure 1b), indicating that its contribution to Hpa resistance would at most be marginal. Furthermore, our transcriptome analysis revealed that the largest variation in gene expression between epiRILs and the Wt line comes from the transcriptional response to Hpa, rather than differences in basal gene expression (Figure 2b–c). Moreover, the genes in Group 2, which displayed enhanced Hpa-induced expression in the resistant epiRILs at the critical early time-point of 48 hpi, were strongly enriched with defence-related GO terms (Figure 2d). The majority of these Group 2 genes showed a statistically significant interaction between epi-genotype and Hpa treatment (Supplementary file 1 dataset S7), indicating that these epiRILs were primed to activate defence-related genes. This notion was supported by the actual expression profiles of Group 2 genes (Figure 2c; Figure 2—figure supplement 4), as well as the defence phenotypes of the eight most resistant epiRILs in the population (Figure 2a; Figure 1—figure supplement 4). Furthermore, our epiRIL screen for growth phenotypes demonstrated that the resistance-controlling epiQTLs do not have a major impacts on plant growth (Figure 1b), which is consistent with previous findings that defence priming is a low-cost defence strategy (van Hulten et al., 2006). While we cannot exclude other mechanisms, these independent lines of evidence collectively indicate that genome-wide priming of defence genes is the most plausible mechanism by which the epiQTLs mediate quantitative disease resistance in the population.

Over recent years, various studies have established a link between DNA hypomethylation and plant immune priming (Espinas et al., 2016; Conrath et al., 2015; López Sánchez et al., 2016). However, causal evidence that heritable regions of reduced DNA methylation mediate transgenerational disease resistance is lacking. Our study has shown that heritable regions of hypomethylated DNA are sufficient to mediate resistance in a genetic Wt background. Furthermore, our study is the first to link phenotypic and epigenetic variation of selected epiRILs to profiles of global gene expression, revealing that epigenetically controlled resistance is associated with genome-wide priming of defence-related genes (Figure 2b–d; Figure 2—figure supplement 1; Figure 2—figure supplement 4). The majority of these pathogenesis-related genes showed augmented induction at 48 hpi (Figure 2c), which represents a critical early time-point in the interaction between Arabidopsis and Hpa, during which hyphae from germinating spores start to penetrate the epidermal cell layer and invade the mesophyll (Koch and Slusarenko, 1990; Soylu and Soylu, 2003). Notably, this set of primed genes was substantially more enriched in SA-dependent and SA-independent defence GO terms than the set of Hpa-inducible genes that were constitutively up-regulated in Hpa-resistant epiRILs (Figure 2d), corroborating the analysis of phenotypical defence markers (Figure 2a; Figure 1—figure supplement 4).

DNA methylation of TEs has been reported to cis-regulate expression of nearby genes in Arabidopsis (Soppe et al., 2000; Saze and Kakutani, 2007; Kinoshita et al., 2007; Lei et al., 2015; Williams et al., 2015). By contrast, our study did not find evidence that DNA methylation in the epiQTLs cis-regulates the responsiveness of nearby of defence genes. Firstly, the majority of primed defence genes in the Hpa-resistant epiRILs were located outside the epiQTL intervals (92%). Secondly, of all primed genes within the epiQTLs, only 73 showed augmented induction that coincided with DNA hypomethylation (Figure 3a; Figure 3—figure supplement 1; Figure 3—figure supplement 2; Supplementary file 1 dataset S10). Of these, 67 encoded TEs or TE-related genes, while the six protein-encoding genes were closely associated with one or more TEs and did not have functions related plant defence (Figure 3a; Figure 3—figure supplement 1; Supplementary file 1 dataset S10). Since TEs do not encode defence signalling proteins, we propose that DNA hypomethylation at the TE-rich epiQTLs mediates augmented induction of defence genes across the genome via trans-acting mechanisms. A recent transcriptome study of Hpa-infected Arabidopsis identified 166 defence-related genes that were primed in the hypomethylated nrpe1-11 mutant and/or repressed in hyper-methylated ros1-4 mutant (López Sánchez et al., 2016). The majority of these defence genes were not targeted by NRPE1- and/or ROS1-dependent DNA (de)methylation, indicating that their responsiveness is trans-regulated by DNA methylation. Although NRPE1 and ROS1 target partially different genomic loci than DDM120, this study supports our hypothesis that DNA methylation controls global defence gene responsiveness via trans-acting mechanisms.

There are various mechanisms by which DNA methylation could trans-regulate defence gene expression. It is possible that transcribed TEs in the hypomethylated epiQTLs generate 21-22nt or 24nt small RNAs (sRNAs) that influence distant heterochromatin formation through via RDR6- and DCL3-dependent RdDM pathways (Panda et al., 2016). Support for trans-regulation by sRNAs came from a recent study, which reported that induction and subsequent re-silencing of pericentromeric TEs in Arabidopsis upon Pseudomonas syringae infection is accompanied with accumulation of RdDM-related sRNAs that are complementary to TEs and distal defence genes. Interestingly, while the accumulation of these sRNAs coincided with re-silencing of the complementary TEs, the complementary defence genes remained expressed in the infected tissues (Cambiagno et al., 2018). These findings are supported by another recent study, which demonstrated that AGO1-associated small RNAs can trans-activate distant defence gene expression through interaction with the SWI/SNF chromatin remodelling complex (Liu et al., 2018). Apart from sRNAs, it is also possible that long intergenic noncoding RNAs (lincRNAs) from the hypomethylated epiQTLs regulate pathogen-induced expression of distant defence genes. A recent study revealed that pericentromeric TEs of Arabidopsis can produce DDM1-dependent lincRNAs that are increased by abiotic stress exposure (Wang et al., 2017a). Since lincRNAs can promote euchromatin and heterochromatin formation at distant genomic loci (Heo et al., 2013; Quinodoz and Guttman, 2014), hypomethylated TEs within the epiQTLs could generate priming-inducing lincRNAs. While knowledge about lincRNAs in plants remains limited, like sRNAs, their activity depends on sequence complementary with target loci (Vance and Ponting, 2014). Unlike non-coding RNAs, long-range chromatin interactions can trans-regulate gene expression independently of sequence complementarity (Harmston and Lenhard, 2013; Liu, 2016; Weber et al., 2016; Wang et al., 2017b). Previous high-throughput chromosome conformation capture (Hi-C) analysis revealed that the ddm1-2 mutation has a profound impact on long-range chromatin interactions within and beyond the pericentromeric regions (Feng et al., 2014). Projection of these DDM1-dependent interactions onto the Arabidopsis genome shows extensive coverage of the resistance epiQTLs identified in this study (Figure 3—figure supplement 5). Whether these long-range interactions contribute to trans-regulation of defence gene priming would require further study, including a fully replicated Hi-C analysis of the resistant epiRILs characterised in this study.

In conclusion, our study has shown that heritable DNA hypomethylation at selected pericentromeric regions controls quantitative disease resistance in Arabidopsis, which is associated with genome-wide priming of defence-related genes. This transgenerational resistance is not associated with reductions in plant growth (Figure 1b), nor does it negatively affect resistance to other types of (a)biotic stresses tested in this study (Figure 1—figure supplement 3). However, whether this form of epigenetically controlled resistance can be exploited in crops depends on a variety of factors, including the stability of the disease resistance and potential non-target effects. For instance, our experiments with Arabidopsis revealed that the resistance has limited stability and can erode over one more generation in some epiRILs (Figure 1—figure supplement 5). Furthermore, the genomes of most crop species contain substantially higher numbers of TEs, rendering predictions about the applicability and potentially undesirable side effects on growth and seed production uncertain. Future research will have to point out whether introgression of hypomethylated pericentromeric loci into the background of elite crop varieties allows for selection of meta-stable quantitative disease resistance without side-effects on agronomically important traits.

Materials and methods

Plant material and growth conditions

Epigenetic recombinant inbred lines (epiRILs) seeds of Arabidopsis (Arabidopsis thaliana, accession Col-0) were purchased from Versailles Arabidopsis Stock Centre, INRA, France (http://publiclines.versailles.inra.fr/epirils/index). The epiRIL screen included siblings of the F4 ddm1-2 parental plant of the epiRIL population (IBENS, France). Arabidopsis seeds were stratified in water at 4°C in the dark for 3-5 days. For pathogen bioassays, seeds were sown in a sand:compost mixture (1:3) and grown at short-day conditions for 3 weeks (8.5 hr light/15.5 hr dark, 21°C, 80% relative humidity,~125 µmol s−1 m−1 light intensity). To test transgenerational inheritance and stability of Hpa resistance in the eight most resistant epiRILs (Figure 1—figure supplement 5), five individual F9 plants were cultivated for 4 weeks at short-day conditions and then moved to long-day conditions to initiate flowering (16 hr light/8 hr dark, 21°C, 80% relative humidity,~125 µmol s−1 m−1 light intensity). Seeds of the 40 F10 families were collected for analysis of Hpa resistance (see below).

Screen for variation in disease resistance and seedling growth

Three-week-old seedlings were spray-inoculated with a suspension of asexual conidia from Hyaloperonospora arabidopsidis strain WACO9 (Hpa) at a density of 105 spores/ml. Hpa colonizsation was quantified at 6 days post-inoculation (dpi) by microscopic scoring of leaves, as described previously (López Sánchez et al., 2016). Briefly, trypan blue-stained leaves were analysed with a stereomicroscope (LAB-30, Optika Microscopes) and assigned to 4 Hpa colonisation classes: class I, no hyphal colonisation; class II, ≤50% leaf area colonized by pathogen hyphae without formation of conidiophores; class III,≤75% leaf area colonized by hyphae, presence of conidiophores; class IV, >75% leaf area colonized by the pathogen, abundant conidiophores and sexual oospores (Figure 1—figure supplement 1). At least 100 leaves per (epi)genotype were analysed, not including the cotyledons. Statistically significant differences in frequency distribution of Hpa colonisation classes between lines were determined by Pearson’s Chi-squared tests, using R (v.3.5.1). Growth analysis of the epiRIL population was based on digital photos (Canon 500D, 15MP) of 3-week-old plants, which were taken on the day of Hpa inoculation. Digital image analysis of total green leaf area (GLA) was performed using Adobe Photoshop 6.0. Green pixels corresponding to GLA were selected and converted into mm2 after colour range adjustment, using the magic wand tool.

Mapping of epigenetic quantitative trait loci (epiQTLs)

Mapping of epiQTLs was performed using the ‘scanone’ function of the R/qtl package for R (Broman et al., 2003) (Haley-Knott regression, step size: 2 cM), combining experimental phenotypical data with the recombination map of differentially methylated regions (DMR) generated previously (Colomé-Tatché et al., 2012). For analysis of Hpa resistance, the categorical scoring of Hpa resistance was first converted into a numeric resistance index (RI), using the following formula:

RI=(fclassI4)+(fclassII3)+(fclassIII2)+(fclassIV1)

where f = relative frequency of Hpa colonisation class of each line, multiplied by an arbitrary weight value ranging from four for the most resistant category (class I) to one for most susceptible category (class IV). Mapping of epiQTLs controlling plant growth was based on average GLA values of each line before Hpa infection. A logarithm of odds (LOD) threshold of significance for each trait was determined on the basis of 1000 permutations for each dataset (α = 0.05). The proportion of phenotypic variance R2 (G) explained by the DMR markers with the highest LOD score (peak markers) of all four epiQTLs was calculated with the following formula (Cortijo et al., 2014):

R2G=1-n-1n-k+1in(yi-[β^0+jkβjgij])2inyi-y-2

where n = number of lines analysed, k = number of DMR markers tested; β0 = intercept of the multiple regression model; βj= QTL effect for each QTL j (slopes for each marker in the multiple regression model); gij = (epi) genotype of the jth marker for each individual i (coded as ‘1’ for ddm1-2 epialleles and ‘-1’ for WT epialleles); yi = phenotypic value of individual i; y- = mean of phenotypic values. The contribution of each individual QTL j (R2(g))was calculated, using the following formula:

R2g=1-n-1n-k+1β^j2ingij-g-j2inyi-y-2 ,

as described by (Cortijo et al., 2014), where n= number of lines analysed, = number of markers tested; β= QTL effect for each QTLj (slopes for each peak marker in the multiple regression model); gi= (epi)genotype of the jth marker for each individual i (coded as ‘1’ for ddm1-2 epialleles and ‘-1’ for WT epialleles); g-j = average of the (epi)genotypes values for the jth marker. Covariance was calculated by subtracting the sum of the individual contributions of each QTL j on phenotypical variance (i.e. R2(gQTL1) + R2(gQTL2) + R2(gQTL4) + R2(gQTL5)) from the phenotypical variance explained by the full model (i.e. R2(G)).

Analysis of shared transposition events

TE-tracker software was used to interrogate available Illumina whole-genome sequencing data from 122 epiRILs for the presence of >2 shared transposition evens (STEs) within the epiQTLs intervals (Gilly et al., 2014). STEs were analysed for statistically significant linkage with resistance phenotypes (RIs), using the same linear regression model as described above for DMR linkage analysis.

Plectosphaerella cucumerina pathoassays

Plectosphaerella cucumerina (Pc, strain BMM (Ton and Mauch-Mani, 2004)) was grown from frozen agar plugs (−80° C) on potato dextrose agar (PDA; Difco, UK). Inoculated plates were maintained at room temperature in the dark for at least 2 weeks. Spores were gently scraped from water-inundated plates, after which spore densities were adjusted to 106 spores/ml using a hemocytometer (Improved Neubauer, Hawksley, UK). Four fully expanded leaves of similar age from 5-week-old plants were inoculated by applying 5 µl droplets, minimising variability due to age-related resistance. After inoculation, plants were kept at 100% RH until scoring of lesion diameters. Average lesion diameters at nine dpi were based on four leaves per plant from 12 plants per (epi)genotype (n = 40–48), using a precision caliper (Traceable, Fischer Scientific). Statistically significant differences in necrotic lesions diameter (asterisks) were quantified by two-tailed Student’s t-test (p<0.05) in pairwise comparisons with Wt line (#602), using R (v3.5.1).

Salt stress tolerance assays

Seeds were sterilised by exposure for 4 hr (h) to chlorine vapours from a 200 ml bleach solution containing 10% v/v hydrochloric acid (37% v/v HCl, Fischer Scientific, 7732-18-5). Seeds were air-dried for 1 hr in a sterile laminar flow cabinet and plated on half strength MS plates (Duchefa, M0221;+0.05% w/v MES,+1% w/v sucrose, pH 5.7), containing increasing concentrations of NaCl (0 mM, 50 mM, 75 mM and 100 mM; Fischer Scientific, 7647-14-5). Plates were stratified for 4 days in the dark at 4°C and transferred to short-day growth conditions (8.5 hr light/15.5 hr dark, 21°C, 80% RH, light intensity 100–140 µmol s−1 m−1). Salt tolerance was expressed as percentage of seeds producing fully expanded cotyledons by 6 days after stratification. Germination percentages of epi-genotypes were calculated from >50 seeds per treatment. Statistically significant differences in germination rates (asterisks) were quantified by Fisher’s exact test (p<0.05) in pairwise comparisons with Wt line (#602) at each salt concentration, using R (v3.5.1).

Quantification of callose effectiveness against Hpa infection

Seedlings were collected at three dpi and cleared for >24 hr in 100% ethanol. One day prior to analysis, samples were incubated for 30 min in 0.07 M phosphate buffer (pH 9), followed by 15 min incubation in a 4:1 mixture (v/v) of 0.05% w/v aniline blue (Sigma-Aldrich, 415049) in 0.07M phosphate buffer (pH 9) and 0.025% w/v calcofluor white (Fluorescent brightener 28, Sigma-Aldrich, F3543) in 0.1M Tris-HCL (pH 7.5). After initial staining, samples were incubated overnight in 0.5% w/v aniline blue (Sigma-Aldrich, 415049) in 0.07M phosphate buffer (pH 9) and scored with an epifluorescence microscope (Olympus BX 51) fitted with blue filter (XF02-2; excitation 330 nm, emission 400 nm). Germinated conidia (germ tubes) were divided between in two classes: non-arrested and arrested by callose. In each assay, 10 leaves from different plants for each (epi)genotype were analysed, amounting to >150 conidia-callose interactions. Statistically significant differences in resistance efficiency of callose (asterisks) were analysed using Pearson’s Chi-squared tests (p<0.05) in pairwise comparisons with Wt line (#602), using R (v3.5.1).

Reverse-transcriptase quantitative polymerase chain reactions (RT-qPCR)

Three biologically replicated samples for each genotype/treatment/time-point combination were collected at 48 and 72 hpi, each consisting of six to 12 leaves collected from different plants per pot. Samples were snap-frozen in liquid nitrogen and ground to a fine powder, using a tissue lyser (QIAGEN TissueLyser). Total RNA was extracted using a guanidinium thiocyanate-phenol-chloroform extraction isolation protocol. Frozen powder was vortexed for 30 s in 1 ml Extraction buffer: 1M guanidine thiocyanate (Amresco, 0380), 1M ammonium thiocyanate (Sigma-Aldrich, 1762-95-4), 0.1M sodium acetate (Fisher Scientific, 127-09-3), 38% v/v AquaPhenol (MP Biomedicals, 108-95-2) and 5% v/v glycerol (Fisher Scientific, 56-81-5). Samples were incubated at room temperature (RT) for 1 min and then centrifuged for 5 min at 16,500 g. The supernatant was then transferred to a new tube, mixed with 200 μl chloroform and vortexed for 10–15 s. After centrifuging for 5 min (16,500 g), the aqueous phase was transferred to new tubes, gently mixed by inversion with 350 μl 0.8M sodium citrate (Sigma-Aldrich, 6132-04-3) and 350 μl isopropanol (Fischer Chemicals, 67-63-0) and left at RT for 10 min for RNA precipitation. Samples were centrifuged for 15 min at 16,500 g (4°C), after which pellets were washed twice in 1 ml 70% ethanol, centrifuged at 16,500 g for 1 min, and air-dried before dissolving in 50 μl nuclease-free water. Total RNA was quantified, using a Nanodrop 8000 Spectrophotometer (Thermo Scientific). RNA extracts were treated with DNaseI, using the RQ1 RNase-Free DNase kit (Promega, M6101). First-strand cDNA synthesis was performed from 1 μg RNA, using SuperScript III Reverse Transcriptase (Invitrogen, 18080093) according to the supplier’s recommendations. The qPCR reactions were carried out with a Rotor-Gene Q real-time PCR cycler (Qiagen) and the Rotor-Gene SYBR Green PCR Kit (Qiagen, 204074). Relative PR1 gene expression was calculated, using Livak’s ΔΔCT method (Livak and Schmittgen, 2001) with correction for average PCR efficiencies for each primer pair across experiment samples. Gene expression was normalised against average expression values of At1G13440 (GAPDH), At5G25760 (UBC) and At2G28390 (SAND family protein) (Czechowski et al., 2005). Reactions were performed using previously described primer sequences (López Sánchez et al., 2016). Statistically significant differences in relative expression (asterisks) were quantified by two-tailed Student’s t-test (p<0.05) in pairwise comparisons with Hpa-treated Wt line (#602).

Transcriptome analysis

Samples for RNA sequencing were collected at 48 and 72 hpi of 3-week-old plants. Every epi-genotype/treatment/time-point combination was based on three biologically replicated samples, each consisting of 6–12 shoots from different plants. Initial RNA extraction was performed as described for RT-qPCR reactions. Prior to library preparation, RNA concentration and integrity were measured, using 2100 Bioanalyzer (Agilent) with provided reagents kits and according to manufacturer’s instructions. All RNA samples had RNA integrity numbers (RIN) >7.5. Sequencing libraries were prepared from total RNA, using the TruSeq Stranded Total RNA kit and Ribo-Zero Plant leaf kit (Illumina, RS-122–2401), according to the manufacturer’s instructions. Sequencing runs were performed on a HiSeq1500 platform (Illumina), generating paired-end reads of 125 bp and an average quality score (Q30) >93%. Each sample generated around 35 million paired reads.

Read quality was assessed by FastQC software (Andrews, 2010). Read length and distribution were optimized and adapter sequences were trimmed, using Trimmomatic software (Bolger et al., 2014). Reads were aligned and mapped to the Arabidopsis genome (TAIR10 annotation), using splice site-guided HISAT2 alignment software (John Hopkins University, second iteration of (Kim et al., 2015)). For all samples, more than 95% of reads could successfully be mapped once or more onto the Arabidopsis genome. Number of reads per gene were quantified with the Python package HTseq (Anders et al., 2015). Differential expression analysis was performed using the DESeq2 R package, which applies a negative binomial generalized linear model to estimate mean and dispersion of gene read counts from the average expression strength between samples (Love et al., 2014). Prior to principal component analysis (PCA) by the plotPCA function, gene read counts were subjected to regularized logarithmic transformation, using the rlog function (Love et al., 2014). Likelihood ratio tests of variance within a three-factorial linear model for epigenotype, treatment, time-point and interactions thereof were used to identify genes showing differences in expression across one or more factors (Love et al., 2014). Differentially expressed genes (DEGs) were subjected to hierarchical clustering (Ward method) and presented as a heat map, using the pheatmap R package (Kolde, 2015). For each gene, rlog-normalized read counts of each sample were subtracted from the mean of all samples, and divided by the standard deviation to facilitate heatmap visualisation (z-score). To identify DEGs between two treatment/time-point/epi-genotype combinations, pair-wise comparisons (Wald test; q < 0.05) were performed with the DEGs selection obtained by the lrt test, using the selection criteria illustrated in Figure 2—figure supplement 2a. All Hpa-inducible genes in the Wt and/or epiRILs were selected for elevated expression in the more resistant epiRILs during Hpa infection. Subsequently, these genes were divided between two groups based on their expression profile. Group 1 genes were selected for constitutively enhanced expression in the epiRIL(s) relative to the Wt (Figure 2—figure supplements 2 and 3); Group 2 genes were selected for enhanced levels of Hpa-induced expression in the epiRIL(s) relative to the Wt (Figure 2—figure supplements 2 and 4). To determine the number of Group 2 genes that show a statistically significant interaction between epigenotype x Hpa treatment (all 16,009 genes significant for this interaction were selected from the three-factorial linear model, using the contrast function, and cross-referenced against Group 2 genes.

Gene ontology (GO) term enrichment analysis was performed with the Plant GSEA toolkit (Yi et al., 2013). GO terms were checked for significant enrichment against the whole genome background, using a hypergeometric test and Benjamini-Hochberg false discovery rate correction (q < 0.05). Lists of enriched GO terms in each treatment were analysed by the GO Trimming 2.0 algorithm (Jantzen et al., 2011) to remove redundancy of terms, applying a soft trimming threshold of 0.40. The output list from GO Trimming 2.0 was run through GOSlim Viewer (AgBase) to reduce GO terms according to GO slim ontologies (GO consortium). Enrichment was quantified as the percentage of GO term-annotated genes within a certain selection relative to the total number of Arabidopsis genes in that GO term.

Methylome analysis

For each line, three independent biological replicates were collected, consisting of pooled leaves from six plants of the same developmental stage. High-quality genomic DNA was extracted from leaves of 5-week-old plants, using the GenElute Plant Genomic DNA Miniprep Kit (Sigma-Aldrich). Bisulfite sequencing was performed by GATC Biotech (UK). After quality trimming of read sequences, adapter sequences were removed, and reads were filtered by Cutadapt (version 1.9; Pair end-mode; phred score = 20, min.length = 40). Reads were mapped to an index genome, using of BS-Seeker2 (version 2.0.10, mismatch = 0.05, maximum insert size = 1000 bp). Bowtie2 (version 2.2.2) was used for alignment of reads, as described previously (Lauss et al., 2018). Differential methylation for promoter regions (−2 kb), gene bodies, and downstream regions (+1 kb) relative to the Wt was called using methylkit (version 1.0.0; minimum coverage = 5 x, q = 0.05). Differentially methylated states were visualised as a heat map, using the ‘pheatmap’ R package (version 1.0.8) (Kolde, 2015).

To differentiate Wt methylation states of all epiQTL-based genes in Group 2 (see above), gene bodies of all nuclear genes were categorised between un-methylated, gene body methylated (gbM; CG context only) or TE-like methylated (teM; CHG and/or CHH with or without CG). For each gene containing 20 or more cytosines, methylated and un-methylated cytosine base calls in each context were extracted from the sequence read alignments. Positions with less than 4x coverage were ignored. Methylation patterns were categorised as TE-like if methylated read calls relative to un-methylated read calls in CHG and/or CHH contexts showed a statistically significant increase over average methylation rates of all genes across the genome in the respective context, using the ‘binom.test’ function in R (FDR-adjusted p<0.01). The remaining genes were classified either as gbM if the same test revealed a statistically significant increase in CG context, or as un-methylated if no statistically significant increase in DNA methylation could be detected in any sequence context.

Correlation analysis between gene expression and DNA methylation

Correlations between augmented expression ratio of Group 2 genes (see Transcriptome analysis) and DNA hypomethylation (CG), were determined by plotting augmented gene ratios at 48 hpi against average hypomethylation compared to Wt (%) across promoter region, gene body, and downstream region (see Methylome analysis). To determine which type of DNA hypomethylation correlates with augmented expression in the epiRILs, hypomethylation at gene bodies of Group 2 genes were divided between teM and gbM and plotted against the corresponding expression ratios at 48 hpi. If hypomethylation occurred at CG context only, genes were classified as being reduced in gene body methylation (gbM); if hypomethylation occurred all three sequence contexts (CG, CHG, CHH), genes were classified as being reduced in TE methylation (teM). Values of gbM hypomethylation were expressed as percentage reduction in GC methylation relative to the Wt; values of teM hypomethylation were expressed as percentage reduction in all sequence contexts. Linear regression analyses were performed using R software (v.3.5.1).

Hi-C analysis

HiC sequence libraries SRR1504819 and SRR150482464 were downloaded from NCBI SRA. Sequences were pre-processed and aligned to the TAIR10 Arabidopsis nuclear genome sequence (Berardini et al., 2015), using HiCUP (0.5.9) (Wingett et al., 2015) and Bowtie2 (Langmead and Salzberg, 2012) (2.2.6). Alignments were filtered and de-duplicated as part of the processing by HiCUP, before being further processed in HOMER (Heinz et al., 2010) (4.9.1) at 5 kb resolution. Differential interactions were assessed reciprocally, using each sample as background (analyzeHiC-ped). Interactions were determined to be potentially dependent on genotype if the absolute z-score of the primary versus the secondary experiment was more than 1. Visualisations were generated using Circos (Krzywinski et al., 2009)(0.69–5), based on bundled links (-max_gap 10001).

Data availability

Transcriptome sequencing and bisulfite sequencing reads are available from the European Nucleotide Archive (ENA) under accession code PRJEB26953.

References

  1. 1
  2. 2
  3. 3
    FastQC: A Quality Control Tool for High Throughput Sequence Data
    1. S Andrews
    (2010)
    Babraham Bioinformatics.
  4. 4
  5. 5
  6. 6
  7. 7
  8. 8
  9. 9
  10. 10
  11. 11
  12. 12
  13. 13
  14. 14
  15. 15
  16. 16
  17. 17
  18. 18
  19. 19
  20. 20
  21. 21
  22. 22
  23. 23
  24. 24
  25. 25
  26. 26
  27. 27
  28. 28
  29. 29
  30. 30
  31. 31
  32. 32
  33. 33
  34. 34
  35. 35
    Pheatmap:Pretty Heatmaps
    1. R Kolde
    (2015)
    GitHub.
  36. 36
  37. 37
  38. 38
  39. 39
  40. 40
  41. 41
  42. 42
  43. 43
  44. 44
  45. 45
    Genome-wide analysis of local chromatin packing in arabidopsis thaliana
    1. C Liu
    (2016)
    Genome Research 25:246–256.
  46. 46
  47. 47
  48. 48
    Two distinct Coagulase-dependent barriers protect staphylococcus aureus from neutrophils in a three dimensional in vitro infection model
    1. A Lopez
    2. V Ramirez
    3. J Garcia-Andrade
    4. V Flors
    5. P Vera
    (2011)
    In PLoS Genetics 7:e1002434.
  49. 49
  50. 50
  51. 51
  52. 52
  53. 53
  54. 54
  55. 55
  56. 56
  57. 57
  58. 58
  59. 59
  60. 60
  61. 61
  62. 62
  63. 63
  64. 64
  65. 65
  66. 66
  67. 67
  68. 68
  69. 69
  70. 70
  71. 71
  72. 72
    DNA Methylation: Basic Mechanisms
    1. BF Vanyushin
    (2006)
    W. alter Doerfler, P Böhm, editors. Berlin: Springer.
  73. 73
  74. 74
  75. 75
  76. 76
  77. 77
  78. 78
  79. 79
  80. 80
  81. 81

Decision letter

  1. Daniel J Kliebenstein
    Reviewing Editor; University of California, Davis, United States
  2. Christian S Hardtke
    Senior Editor; University of Lausanne, Switzerland

In the interests of transparency, eLife includes the editorial decision letter and accompanying author responses. A lightly edited version of the letter sent to the authors after peer review is shown, indicating the most substantive concerns; minor comments are not usually included.

Thank you for submitting your article "Identification and characterisation of epigenetic loci controlling transgenerational immune priming in Arabidopsis" for consideration by eLife. Your article has been reviewed by three peer reviewers, one of whom is a member of our Board of Reviewing Editors, and the evaluation has been overseen by Christian Hardtke as the Senior Editor. The reviewers have opted to remain anonymous.

The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.

Summary:

This work utilizes a collection of epiRILs to show that there is potential epigenetic variation in controlling defense responses.

Essential revisions:

After discussion, there were considered two essential components that need addressing to make the manuscript publishable in eLife.

First, it was felt that the work is highly interesting in itself but that at a number of locations, claims were made that were not supported by the data at best and occasionally contradictory information was available in other papers. As such, it was agreed that the authors should use the guidance in the reviews to rewrite the manuscript to more precisely focus on what has and has not been shown in this manuscript and others.

Second, we all agreed that the Hi-C was a potentially interesting idea but that the claims were unsupported. To support the claims would require doing Hi-C in the specific epiRILs or drop it from the manuscript at this point in time.

Reviewer #1:

– If they were just honest about what they found, it would be a really nice paper and I would support it for eLife. But their claims are simply unjustified by their analysis.

The authors map epiQTLs that link to resistance to Hpa, then assess potential epigenetic influences on transcripts and how this may be priming. Finally, the authors look for chromatin interactions in WT/ddm1 and use them to argue that this is the mechanism. The data itself is very nice. but I find that the claims are not justified based on the conducted analysis. The experiments and data themselves will be very useful but the authors would need to focus on what they do show. To illustrate, they nicely map resistance epiQTLs and do some really nice transcriptomics. Just presenting that as plainly as what they found would be interesting.

In contrast, the writing makes extensive claims about priming, transgenerational and chromatin interactions. The transgenerational argument is trying to link in with the work showing that there is one-generation of disease effects on resistance. However, there is no analysis of if their loci really play any role in this phenomena. They just have loci whose ddm1 marks lead to altered resistance. They would have to show that these same loci influence the transgenerational resistance phenomena to make their arguments.

The chromatin interactions are all done in ddm1 v. WT and not in the epiRILs v. WT or v. ddm1 so I'm not sure they can actually make their claims. Wouldn't they need to look at interactions in the epiRILs to make any of these arguments.

Priming concerns –

I should caveat the following in that I have some quant/eco/evo training in my background so I go by a very strict definition of priming; a pre-treatment leads to an elevated response to a secondary treatment that cannot be explained by any effect of the pre-treatment. Similar to the definition I work off of for epistasis. The effect of A x B is not explainable by adding up the effect of A and B separately.

Given this definition, the transcriptomics analysis does not make a direct test of priming. The authors use pairwise analysis which do not allow for a statistical testing of priming directly. To claim priming would require re-analyzing the data with a linear model where Transcript = genotype + treatment + genotype x treatment. And only the transcripts with a significant genotype x treatment term would fit into what might be called priming. Priming is the interaction term rather than the main effect of genotype or treatment.

This concern about what is or is not priming, arises from their own data that argues against priming even when they say it argues for priming. If you notice in Figure 2B, the PCA plot actually shows that the main transcript variance in the epiRILs are not priming but mainly constitutive changes as the change from Mock to Hpa infected is a similar level of variance in each genotype. If it was truly a predominantly priming effect, the uninduced would all be similar to the WT uninduced. Figure 2—figure supplement 2A also shows my worry about priming. If you look at the plots for what a primed transcript looks like, only 2 of those five plots would truly fit what is priming.

The authors argue that priming is the sole thing altered and that other defense traits aren't influenced. There are a couple of potential issues with this exclusionary arguments. The use of the 8 most resistance epiRIL lines to test other interactions is not really evidence of much. There isn't enough statistical power to say anything other than they aren't exactly 100% the same traits. It is possible that there is still a correlation. Additionally, at least one if not two of their loci are linked to variation in defense metabolism in this population based on previous work looking at epigenetic variation in defense metabolites within population. This shows that there are defense traits that are altered in these lines prior to infection. This previous work should be cited and incorporated appropriately when discussing the authors results as it contradicts the central claim about exclusivity of effects.

Reviewer #2:

This study uses the ddm1 epiRIL population to identify QTL controlling resistance to. At least four epiQTL are identified which contribute >50% of the heritable variation. To explore these epiQTL further, transcriptome and methylome analyses were performed. These experiments revealed that four of the epiRIL lines with the highest resistance display a heightened transcriptome response to Hpa indicating these lines are already primed. This enhanced priming explains why the lines are not suffering a growth defect as there is not a constitutive response to Hpa or other a/biotic stresses. To understand how the priming of Hpa induced genes occurs the DNA methylome data were integrated, which generally revealed that the loci within the QTL regions were not directly affected by loss of DNA methylation and that perhaps there is a genome-wide trans affect occurring in these lines due to the global loss of DNA methylation. To provide a mechanistic explanation for how this could occur the authors hypothesize that long range chromatin interactions could be in play. To test this, Hi-C data from WT and ddm1 was analyzed to revealed the existence of long range chromatin interaction between the QTL and the affected defence genes.

Overall, this is an interesting study, especially with the identification of epiRILs that are resistant to Hpa seemingly with no growth defects. However, the causal basis for these phenotypes requires additional evidence. The use of Hi-C is innovative, but it does not prove the strong statement made in the subsection “DDM1-dependent chromatin interactions between the resistance epiQTLs and distant defence genes as a potential trans-priming mechanism”, in the Abstract and in the Discussion.

1) Throughout the manuscript there are references to phrases that do not make sense. What are "epigenetic loci" as stated in the title or "epigenetic resistance". I would refer to these as "loci" and "resistance". Furthermore, why is "transgenerational" invoked throughout the paper? QTL are stably inherited. There are no transgenerational experiments presented in this study. Please remove this phrase from the study unless referring to the literature.

2) There is a potentially interesting observation of a trans effect being causal for the resistance, but this has not been proven as stated in the Abstract and the paper.

3) The epiQTL span regions associated with 4 out of 5 pericentromeric regions indicating that overall depletion of DNA methylation and potentially loss of heterochromatin is associated with Hpa resistance. What is the association between the amount of ddm1 like chromosomes within the epiRILs and the resistance phenotype. Essentially, if you remade Figure 1 to show the amount of Col vs. ddm1 chromosome would it correlate with the observed resistance phenotype?

4) Why was ddm1 not included in the transcriptome study upon Hpa infection? This would be useful to understand if it is stronger than the epiQTL, which seems to be a prediction based on the discussion.

5) The data presented to support a "tight correlation" between hypomethylation and gene expression in Figure 3B needs additional support. While it is expected that loss of methylation and transposon genes could result in their reactivation, this is not as likely to occur at genes given that most genes are unmethylated, a small percent are gene body methylation (which is causally linked to expression) and rare genes show transposon-like methylation profiles. This can be resolved by splitting out the genes into unmethylated, gene body methylation (CG) and TE-like methylation (CG, CHG and CHH). Although some of these genes are cis-regulated by DNA methylation, it does not look like they all are as stated. This could be further tested using a scatter plot between change in DNA methylation versus change in expression.

6) There is no data provided to support the statement that "their transcriptional profile is determined by their associated TEs." Please rephrase this as a hypothesis.

7) The use of Hi-C data is an interesting idea. It is stated that "43 interactions between the epiQTLs…" were reduced or intensified in ddm1 vs. wt. Yet, there is no statistical support to indicate whether this represents an enrichment over background expectations. This is a critical test as this forms the mechanistic basis for how the trans interactions occur.

Reviewer #3:

I think this paper is likely to be of broad interest. The demonstration that demethylation of particular bits of the genome affects a particular phenotype, in this case immunity, is intriguing. The authors have made a serious effort to understand why, and have made substantial progress in that direction. They have demonstrated increased effectiveness of callose, a known factor in resistance to this pathogen, and they have shown increased expression of genes known to be induced in response to infection. They make use of a new technology, Hi-C, and obtain data suggesting that the altered methylation may have long distance effects on expression levels of genes associated with immunity, possibly explaining the resistance phenotype. I found this rather unexpected, as conventional wisdom is that in the small genome of Arabidopsis, control of gene expression is nearly always very local. This work suggests that may not be correct, and a wider view is needed.

The authors were faced with the considerable challenge of working out how the epiQTLs affect expression of immmunity-related genes. In the last experiment, they use Hi-C and find connections between the epiQTL DNA and certain infection-inducible genes outside the QTL regions. Moreover, those genes are up-regulated in the resistant epiRILs carrying the epiQTL alleles. This is temptingly close to showing a mechanism. However, more statistical analysis is needed to make this convincing. A large number of infection-inducible genes are up-regulated in the resistant epi-RILs. Is the fraction of these detected in the Hi-C experiment as being linked to the QTL loci greater than would be expected by chance?

[Editors' note: further revisions were requested prior to acceptance, as described below.]

Thank you for resubmitting your work entitled "Identification and characterisation of hypomethylated DNA loci controlling quantitative resistance in Arabidopsis" for further consideration at eLife. Your revised article has been favorably evaluated by Christian Hardtke as the Senior Editor, and three reviewers.

The manuscript has been improved but there are some remaining issues that need to be addressed before acceptance, as outlined below. There are just a few requests for clarification in terminology and the like that will help to raise the readability and visibility of the manuscript.

Reviewer #1:

The authors have largely responded to my previous concerns. One thing that I think might help this version of the manuscript is to simplify the Discussion and definitions surrounding augmentation and priming. I found myself getting confused about what was being meant. It would help to have discrete explicit definitions being given to allow the reader to fully understand what the authors mean.

Reviewer #2:

The authors have made significant improvements to this study.

Reviewer #3:

This revised version of the manuscript addresses all of my concerns. I agree with the decision to leave out the Hi-C data. The inclusion of additional experiments has strengthened the work. The conclusions are now more conservative.

https://doi.org/10.7554/eLife.40655.034

Author response

Essential revisions:

After discussion, there were considered two essential components that need addressing to make the manuscript publishable in eLife.

First, it was felt that the work is highly interesting in itself but that at a number of locations, claims were made that were not supported by the data at best and occasionally contradictory information was available in other papers. As such, it was agreed that the authors should use the guidance in the reviews to rewrite the manuscript to more precisely focus on what has and has not been shown in this manuscript and others.

Second, we all agreed that the Hi-C was a potentially interesting idea but that the claims were unsupported. To support the claims would require doing Hi-C in the specific epiRILs or drop it from the manuscript at this point in time.

We have edited various sections in our manuscript to make sure that our claims are supported by the evidence provided. To this end, we have toned down some conclusions and/or provided additional experimental and statistical evidence to support our conclusions. In addition, we have included reference to two recent studies about the role of small RNAs in trans-regulation of defence gene expression (Cambiagno et al., 2018 and Liu et al.,). We also make reference to a recent study by Aller et al., 2018, who have used the same epiRIL population to study the contribution of heritable variation in DNA methylation to the production of defence-related metabolites (glucosinolates). Finally, we have removed the preliminary Hi-C analysis from the Results section of our manuscript. Our revised manuscript introduces the possibility of long-range heterochromatic interactions in the Discussion as one of multiple hypotheses by which the epiQTLs could regulate pathogen-induced expression.

Reviewer #1:

If they were just honest about what they found, it would be a really nice paper and I would support it for eLife. But their claims are simply unjustified by their analysis.

The authors map epiQTLs that link to resistance to Hpa, then assess potential epigenetic influences on transcripts and how this may be priming. Finally, the authors look for chromatin interactions in WT/ddm1 and use them to argue that this is the mechanism. The data itself is very nice. but I find that the claims are not justified based on the conducted analysis. The experiments and data themselves will be very useful but the authors would need to focus on what they do show. To illustrate, they nicely map resistance epiQTLs and do some really nice transcriptomics. Just presenting that as plainly as what they found would be interesting.

We thank the reviewer for her/his supportive comments. We agree that some of our assertions were too speculative or not sufficiently supported by experimental evidence. Please, see our below answers for details.

In contrast, the writing makes extensive claims about priming, transgenerational and chromatin interactions. The transgenerational argument is trying to link in with the work showing that there is one-generation of disease effects on resistance. However, there is no analysis of if their loci really play any role in this phenomena. They just have loci whose ddm1 marks lead to altered resistance. They would have to show that these same loci influence the transgenerational resistance phenomena to make their arguments.

To start with this reviewer’s concern about the term ‘transgenerational’, we used this word to refer to the hypo-methylated epiQTLs that are segregating in the epiRIL population. These loci were inherited from the ddm1-2 parent in a wild-type background and can therefore be regarded as ‘transgenerational’. The epiRILs described in our paper were at the 8th generation of self-pollination in a wild-type background (Johannes et al. 2009). In fact, Latzel et al., 2013, have shown that some hypo-methylated DNA regions from the ddm1-2 mutant are inherited for up to 16 generations of self-pollination in a wild-type background. Since we have established a statistically significant link between the segregating resistance in the epiRIL population and four hypo-methylated epiQTLs that are inherited from the ddm1-2 parent, we believe it is justified to refer to the resistance as ‘transgenerational’. To clarify this argument, we have adjusted the following sections in our manuscript:

“A core set of 123 epiRILs from this population at the 8th generation of self-pollination in the wild-type (Wt) background has been characterized for differentially methylated region (DMR) markers, enabling linkage mapping of heritable hypo-methylated loci controlling root growth, flowering and abiotic stress tolerance (Cortijo et al., 2014; Kooke et al., 2015).”

“In this study, we have characterised the core set of 123 lines from the ddm1-2 x Col-0 epiRIL population for resistance against the biotrophic downy mildew pathogen Hyaloperonospora arabidopsidis (Hpa) to search for heritable hypo-methylated loci controlling disease resistance.”

Furthermore, we have carried out an additional experiment to confirm the transgenerational inheritance of the resistance trait and assess the ‘transgenerational’ stability of the trait. To this end. 40 F9 individuals from 8 Hpa-resistant lines were self-pollinated, after which the resulting F10 progenies were tested for Hpa resistance. As is shown in Figure 1—figure supplement 5 of the revised manuscript, the resistance was still evident in the F10 generation for every line tested, confirming that the resistance is indeed ‘transgenerational’ and passed on from the F9 to the F10 generation.

Interestingly, however, 4 epiRILs displayed statistically significant variation between the five F10 families within the line, suggesting that these epiRILs show transgenerational instability of the trait. Moreover, two of the forty F10 families showed similar levels of Hpa resistance as the wild-type, indicating that they had lost the transgenerational resistance completely.

“The 123 epiRILs analysed for Hpa resistance had been self-pollinated for 8 generations in a Wt (Col-0) genetic background since the F1 x Col-0 backcross (F9) (Johannes et al., 2009). […] Furthermore, 4 of the 8 epiRILs tested (#71, #148, #545, and #508) displayed statistically significant variation in Hpa resistance between the 5 F10 families within the epiRIL (Figure 1—figure supplement 5; Pearson’s Chi-squared test, p<0.05; † symbols), suggesting instability of the Hpa resistance.”

The chromatin interactions are all done in ddm1 v. WT and not in the epiRILs v. WT or v. ddm1 so I'm not sure they can actually make their claims. Wouldn't they need to look at interactions in the epiRILs to make any of these arguments.

We agree that the Hi-C analysis of the data from the Feng et al., 2014 study remains too preliminary to conclude that DDM1-dependent chromatin interactions contribute to defence priming in the resistant epiRILs. Based on the collective feedback from all reviewers and the editor, we decided to remove the Hi-C analysis from the Results section of the manuscript. The revised manuscript refers in the Discussion to the results by Feng et al., 2014, who showed that the ddm1-2 mutation has a major effect on heterochromatic interactions in the pericentromeric regions, including the epiQTLs identified in our study. We illustrate this by referring to Figure 3—figure supplement 4, which projects the data from the Feng et al., 2014 paper onto the Arabidopsis genome and shows ample DDM1-dependent chromatin interactions originating from the epiQTLs. Considering that the same ddm1-2 mutant was used to generate the Col-0 x ddm1-2 epiRIL population, this provides support to the hypothesis that DDM1-dependent long-range chromatin interactions with the pericentromeric epiQTLs could play a role in trans-regulation of defence genes. However, we emphasize that this remains a hypothesis until further Hi-C studies with the resistant epiRILs can confirm statistically significant interactions with the constitutively primed defence genes:

“Unlike non-coding RNAs, long-range chromatin interactions can trans-regulate gene expression independently of sequence complementarity (Harmston and Lenhard, 2013; Liu et al., 2016; Weber et al., 2016; Wang et al., 2017). […] Whether these long-range interactions contribute to trans-regulation of defence gene priming would require further study, including a fully replicated Hi-C analysis of the resistant epiRILs characterised in this study.”

Priming concerns –

I should caveat the following in that I have some quant/eco/evo training in my background so I go by a very strict definition of priming; a pre-treatment leads to an elevated response to a secondary treatment that cannot be explained by any effect of the pre-treatment. Similar to the definition I work off of for epistasis. The effect of A x B is not explainable by adding up the effect of A and B separately.

Given this definition, the transcriptomics analysis does not make a direct test of priming. The authors use pairwise analysis which do not allow for a statistical testing of priming directly. To claim priming would require re-analyzing the data with a linear model where Transcript = genotype + treatment + genotype x treatment. And only the transcripts with a significant genotype x treatment term would fit into what might be called priming. Priming is the interaction term rather than the main effect of genotype or treatment.

Firstly, we would like to point out that there is published evidence that basal resistance in plants can be caused by an increased responsiveness of pathogen-inducible defence mechanisms, which is commonly referred to ‘constitutive priming’ (van Hulten, 2006; Ahmad et al., 2010, MPP, 11: 817-827, Ahmad et al., 2011 PCE, 34: 1191-1206, Gamir et al. 2014; Plant J, 78:227-240). This resistance phenotype has also been described for Arabidopsis mutants with reduced DNA methylation (Lopez et al., 2011; Lopez et al., 2016) and progeny from diseased parental plants that express transgenerational acquired resistance (Luna et al., 2012; Rasmann et al., 2012). Thus, the term priming is not strictly limited to describe an induced resistance response to a prior stimulus within the same generation, but can also be used to describe an innate resistance phenotype, as is the case for the epiRILs described in our current study.

Our arguments to support the conclusion that Hpa-resistant epiRILs are constitutively primed for defence-related gene expression is based on three statistical tests, the visualisation of the associated gene expression patterns, and the GO term analysis (which confirms enrichment of defense-related functions). The first statistical selection is based on a comprehensive linear model that includes all factors and interactions (epi-genotype, treatment, time-point and all interactions thereof). Using this model, a likelihood ratio test (LRT) selected for all differentially expressed genes by any factor and/or interaction, yielding a large number of genes. Although hierarchical cluster analysis of this gene set already indicated a large cluster displaying augmented induction in the epiRILs (Figure 2—figure supplement 1), we performed (Wald tests (+ FDR) to select for Hpa-inducible genes that show elevated expression in the more resistant epiRILs during Hpa infection (Figure 2—figure supplement 2). Wald tests are recommended by the DeSEQ2 vignette to extract specific contrasts in gene expression and allowed us e.g. to separate the Hpa-inducible genes from the Hpa-repressed genes. Using a series of Wald tests, we distinguished two Hpa-inducible gene groups that differ in expression profile between the Wt and epiRILs: in addition to Hpa-inducible genes that are constitutively up-regulated in the resistant epiRILs (Group 1), we identified a larger group of genes that show elevated levels of Hpa-induced expression in the resistant epiRILs (Group 2). This latter gene group was strongly enriched with defence-related GO terms at 48 hpi, which represents a critical time-point in the interaction between Arabidospis and Hpa (Figure 2D). Thirdly, on advice of the reviewer and the editor, we have selected all genes from the general linear model that show a statistically significant interaction between epi-genotype x Hpa treatment (16,009 genes), and determined the percentage of genes in Group 2 that are significant for this interaction term. As is presented in Supplementary dataset S7 in Supplementary file 1, the majority of Group 2 genes show a statistically significant interaction term for each epiRIL (> 92%), supporting a primed expression pattern. Finally, we confirm this notion by visualization of the actual expression profiles in heatmaps (Figure 2C and Figure 2—figure supplement 4), which show clearly that the vast majority of Group 2 genes exhibit augmented induction levels by Hpa in the resistant epiRILs.

Thus, to address this reviewer’s concerns about priming, we have applied the following changes to our manuscript:

- We have adjusted our description of the gene expression profiles. The selection for Hpa-inducible genes showing enhanced expression in the more resistant epiRILs now differentiates between genes that show constitutively enhanced expressed in the resistant epiRILs (now referred as “Group 1”) and genes that show elevated levels of Hpa-induced expression in the resistant epiRILs (now referred as “Group 2”). We have removed the notation of “primed” at this stage of the gene selection.

- On advice of the reviewer, we selected from our comprehensive general linear model all genes that show a statistically significant interaction between epi-genotype x Hpa treatment. This criterion alone is not sufficient to select for primed genes, since it also selects for genes that are e.g. differentially repressed by Hpa infection. However, when applied to the genes within Group 2, which had already been selected for enhanced levels of Hpa-induced expression in the epiRILs, we found that the vast majority (> 92%) of genes in this group were statistically significant for this interaction term. The combination of these criteria indicates a constitutively primed expression pattern.

- We have included z-score-normalised heatmaps of the gene expression profiles in Groups 1 and 2. This visualisation of the actual expression patterns supports our conclusion that the majority of Group 2 genes are constitutively primed in the resistant epiRILs: the genes show strongly augmented gene induction by Hpa in the resistant epiRILs, while showing no or only marginally higher levels of basal expression in mock-treated epiRILs (Figure 2C, Figure 2—figure supplements 3 and 4). This expression profile is consistent with the agreed definition of priming, as described by Martinez-Medina et al.,2016).

- We only refer to ‘defence gene priming’ after the analysis of a statistically significant interaction between epi-genotype x Hpa treatment, visualisation of the actual expression patterns, and the GO term analysis. The latter analysis demonstrates that Group 2 genes at the early and critical time-point of 48 hpi are strongly enriched with defence-related GO terms.

This concern about what is or is not priming, arises from their own data that argues against priming even when they say it argues for priming. If you notice in Figure 2B, the PCA plot actually shows that the main transcript variance in the epiRILs are not priming but mainly constitutive changes as the change from Mock to Hpa infected is a similar level of variance in each genotype. If it was truly a predominantly priming effect, the uninduced would all be similar to the WT uninduced. Figure 2—figure supplement 2A also shows my worry about priming. If you look at the plots for what a primed transcript looks like, only 2 of those five plots would truly fit what is priming.

We would like to point out that we did not conclude from the PCA plot that the epiRILs are primed for defence gene expression. What the PCA plot shows is that samples from mock- and Hpa-inoculated plants separate mostly along the first principal component (PC1), which explains 31% of the variation in the experiment, whereas samples from different epi-genotypes separate along the second principal component (PC2), which explains only 20% of the variation. Thus, the main transcript variation in the experiment is caused by the response of the plants to Hpa. Since the samples from Hpa-inoculated Wt plants cluster differently along PC1 than those from the Hpa-inoculated epiRILs, it is justifiable to conclude that the transcriptome of the epiRILs responds differently to Hpa infection than the Wt.

Regarding this reviewer’s concern about the expression profiles depicted in Figure 2—figure supplement 2A, we removed the label ‘Primed’ from the gene selection, which is now referred as ‘Group 2’. As detailed above, we define genes in this Group as displaying enhanced levels of Hpa-induced expression in the resistant epiRILs. Based on further statistical analysis of the interaction between epigenotype x Hpa treatment and visual confirmation of the actual gene expression profiles, we conclude that the vast majority of Group 2 genes are constitutively primed in the resistant epiRILs.

The authors argue that priming is the sole thing altered and that other defense traits aren't influenced. There are a couple of potential issues with this exclusionary arguments. The use of the 8 most resistance epiRIL lines to test other interactions is not really evidence of much. There isn't enough statistical power to say anything other than they aren't exactly 100% the same traits. It is possible that there is still a correlation.

We agree with the reviewer that the callose phenotypes and PR1 gene expression profiles of the 8 epiRILs are insufficient to conclude that the resistance controlled by the epiQTLs is solely based on defence priming. Accordingly, we have removed the exclusionary tone about the causal link between the epiQTLs and defence priming. However, we can conclude that the 8 most resistant epiRILs in the population are constitutively primed for cell wall defence and SA-dependent PR1 induction. This is supported by the genome-wide transcriptomic expression profiles and related GO term analysis, revealing that the majority of Hpa-inducible genes with defence-related functions show a constitutively primed expression pattern (see above). Furthermore, based on our population-wide screen for growth phenotypes, we can conclude that the resistance-controlling epiQTLs have no major impact on plant growth, which is consistent with the fact that defence priming has been shown to provide plants with resistance that is associated with no or minimal effects on plant growth (van Hulten et al., 2006). Thus, while we cannot exclude that other defence mechanisms have a contribution to the segregating resistance in the epiRIL population, independent lines of evidence in our study suggest that defence priming is the most plausible mechanism by which the epiQTLs control Hpa resistance. As is further detailed in our response to comment #7 of this reviewer (see below), we have addressed this point in the second paragraph of the Discussion section.

Additionally, at least one if not two of their loci are linked to variation in defense metabolism in this population based on previous work looking at epigenetic variation in defense metabolites within population. This shows that there are defense traits that are altered in these lines prior to infection. This previous work should be cited and incorporated appropriately when discussing the authors results as it contradicts the central claim about exclusivity of effects.

We thank the reviewer for pointing us to the Aller et al., 2018 paper; “Comparison of the Relative Potential for Epigenetic and Genetic Variation To Contribute to Trait Stability”. The revised version of our manuscript now discusses the potential contribution of pre-existing defence metabolites (glucosinolates) in the Hpa resistance:

“Aller et al., 2018, have recently used the same Col-0 x ddm1-2 epiRIL population to map the contribution of heritable variation in DNA methylation to the production of defence-related glucosinolate metabolites (Wibowo et al., 2016). […] While we cannot exclude other mechanisms, these independent lines of evidence collectively indicate that genome-wide priming of defence genes is the most plausible mechanism by which the epiQTLs mediate quantitative disease resistance in the population.”

Reviewer #2:

[…] Overall, this is an interesting study, especially with the identification of epiRILs that are resistant to Hpa seemingly with no growth defects. However, the causal basis for these phenotypes requires additional evidence. The use of Hi-C is innovative, but it does not prove the strong statement made in the subsection “DDM1-dependent chromatin interactions between the resistance epiQTLs and distant defence genes as a potential trans-priming mechanism”, in the Abstract and in the Discussion.

1) Throughout the manuscript there are references to phrases that do not make sense. What are "epigenetic loci" as stated in the title or "epigenetic resistance". I would refer to these as "loci" and "resistance".

We agree that some statements and phrases in our original manuscript were unclear. While the term ‘epigenetic quantitative trait loci’ (epiQTLs) is a commonly used term to describe heritable regions of DNA hypo-methylated controlling plant complex traits, we have deleted the term ‘epigenetic’ when referring to the resistance phenotype of the epiRILs. To prevent confusion, we have also adjusted the title of manuscript: “Identification and characterisation of hypomethylated DNA loci controlling quantitative resistance in Arabidopsis”.

Furthermore, why is "transgenerational" invoked throughout the paper? QTL are stably inherited. There are no transgenerational experiments presented in this study. Please remove this phrase from the study unless referring to the literature.

We thank the reviewer for this comment. We have now included an extra experiment that addresses the transgenerational nature and stability of the Hpa resistance. For further details, please, see our response to comment #2 from reviewer 1.

2) There is a potentially interesting observation of a trans effect being causal for the resistance, but this has not been proven as stated in the Abstract and the paper.

We agree with the reviewer that there is no hard causal evidence to support this statement. We have therefore removed this statement from the Abstract. In addition, we have toned down the conclusiveness regarding trans-regulation of defence genes in our paper. Instead, we propose different mechanisms of trans-regulation in the Discussion as possible explanations for our finding that we did not find defence-related genes within the epiQTL regions that are cis-regulated by the DNA methylation (i.e. genes that show augmented expression in the epiRILs during Hpa infection and that are simultaneously reduced in DNA methylation).

3) The epiQTL span regions associated with 4 out of 5 pericentromeric regions indicating that overall depletion of DNA methylation and potentially loss of heterochromatin is associated with Hpa resistance. What is the association between the amount of ddm1 like chromosomes within the epiRILs and the resistance phenotype. Essentially, if you remade Figure 1 to show the amount of Col vs. ddm1 chromosome would it correlate with the observed resistance phenotype?

Each of the 123 epiRILs in the population shows different combinations of the 126 DNA methylation markers, which are stably inherited in a Mendelian fashion from the Wt and ddm1-2 parents (Colomé-Tatché et al., 2012) Since the population is derived from a backcross, on average, the epiRIL population contains 25% of the ddm1 haplotype and 75% of the Col-0 haplotype. Interval mapping provides the appropriate method to establish a statistical relationship between the resistance phenotype of each line and the haplotype for each marker (quantified by LOD score).

To address this reviewer’s question, we have included a supplementary dataset (S1 in Supplementary file 1), which details the DNA methylation marker distribution of the Wt and ddm1-2 haplotypes of all epiRILs, ranked by Hpa resistance. The bottom row of this table shows the estimated percentage of ddm1-2 haplotype in each epiRIL, which indeed shows a weak positive correlation with Hpa resistance. We attribute this trend to the size of the epiQTL intervals (particularly on chromosomes I and IV). However, as mentioned above, it is the actual interval mapping that establishes which parts of the ddm1-2 haplotype contribute to the resistance.

4) Why was ddm1 not included in the transcriptome study upon Hpa infection? This would be useful to understand if it is stronger than the epiQTL, which seems to be a prediction based on the discussion.

We agree that including siblings of the F4 ddm1-2 parent of the epiRIL population would have been informative. However, we did not have enough F4 seeds of this line to ensure sufficient replication and plant material for the transcriptome analysis. Because of the nature of the ddm1-2 mutation, propagation of this mutant into a subsequent generation (F5) would increase the level of DNA hypomethylation (Kakutani et al., 1996). As such, using F5 progeny from a F4 sibling of the ddm1-2 parent would not provide an appropriate estimation of the amount of DNA hypomethylation that had been introgressed into the wild-type background for this epiRIL population.

5) The data presented to support a "tight correlation" between hypomethylation and gene expression in Figure 3B needs additional support. While it is expected that loss of methylation and transposon genes could result in their reactivation, this is not as likely to occur at genes given that most genes are unmethylated, a small percent are gene body methylation (which is causally linked to expression) and rare genes show transposon-like methylation profiles. This can be resolved by splitting out the genes into unmethylated, gene body methylation (CG) and TE-like methylation (CG, CHG and CHH). Although some of these genes are cis-regulated by DNA methylation, it does not look like they all are as stated. This could be further tested using a scatter plot between change in DNA methylation versus change in expression.

We thank the reviewer for this helpful suggestion. The revised manuscript now includes additional scatterplots to visualise and quantify the correlation between augmented expression ratio and DNA hypomethylation of the corresponding Group 2 genes (Figure 3—figure supplement 2). For each epiQTL, linear regression analysis confirmed the positive correlations for genes in expression cluster I, which were statistically significant (Figure 3A and Figure 3—figure supplement 1). By contrast, no positive correlations could be detected for genes in expression cluster II (Figure 3A and Figure 3—figure supplement 1). This extra analysis provides statistical support to our conclusion that cluster I genes are cis-regulated by DNA methylation, whereas the cluster II genes are trans-regulated by DNA methylation.

6) There is no data provided to support the statement that "their transcriptional profile is determined by their associated TEs." Please rephrase this as a hypothesis.

We have toned done the conclusiveness of this statement and modified the sentence as follows: Furthermore, analysis of the genomic context of the sixprotein-coding genes revealed the presence of overlapping and/or nearby TEs (Figure 3—figure supplement 3), suggesting that the correlation between augmented expression and DNA hypomethylation for these genes is determined by association with TEs.”

7) The use of Hi-C data is an interesting idea. It is stated that "43 interactions between the epiQTLs…" were reduced or intensified in ddm1 vs. wt. Yet, there is no statistical support to indicate whether this represents an enrichment over background expectations. This is a critical test as this forms the mechanistic basis for how the trans interactions occur.

We agree that additional analysis to validate a statistically significant enrichment of heterochromatic interactions between the epiQTLs and distant defence genes identified in our RNA-seq analysis would be necessary to conclude that long-range chromatin interactions contribute to priming of distant defence genes. Sadly, the publicly available dataset from Feng et al., 2014, was based on only one replicate per genotype, thereby limiting the statistical stringency of such analysis. In addition, as pointed out by the reviewer #1 and the editor, even if we would be able to provide a statistically significant justification for these specific interactions, the analysis would still only be based on the ddm1-2 mutant, and not on the epiRILs that were analysed for genome-wide gene expression. For these reasons, we have come to the conclusion that the Hi-C analysis remains too preliminary at the current state to be presented as a main result in the manuscript.

In its current form, the Discussion of our revised manuscript raises the possibility that long-range heterochromatic interactions with peri-centromeric epiQTLs could have a contribution to trans-regulation of distant defence genes. We support this hypothesis by referring to the study by Feng et al., 2014, who demonstrated that the ddm1-2 mutation has a profound effect on long-range heterochromatic interactions with pericentromeric regions, including those that are covered by our epiQTLs. To illustrate this point, we refer in the Discussion to Figure 3—figure supplement 4, which presents the interaction maps of the Wt (Col-0) and the ddm1-2 mutant from the Feng et al., 2014 paper, and highlights the inferred DDM1-dependent interactions with our epiQTL intervals. For the reasons outlined above, we do not take this analysis further.

Reviewer #3:

I think this paper is likely to be of broad interest. The demonstration that demethylation of particular bits of the genome affects a particular phenotype, in this case immunity, is intriguing. The authors have made a serious effort to understand why, and have made substantial progress in that direction. They have demonstrated increased effectiveness of callose, a known factor in resistance to this pathogen, and they have shown increased expression of genes known to be induced in response to infection. They make use of a new technology, Hi-C, and obtain data suggesting that the altered methylation may have long distance effects on expression levels of genes associated with immunity, possibly explaining the resistance phenotype. I found this rather unexpected, as conventional wisdom is that in the small genome of Arabidopsis, control of gene expression is nearly always very local. This work suggests that may not be correct, and a wider view is needed.

The authors were faced with the considerable challenge of working out how the epiQTLs affect expression of immmunity-related genes. In the last experiment, they use Hi-C and find connections between the epiQTL DNA and certain infection-inducible genes outside the QTL regions. Moreover, those genes are up-regulated in the resistant epiRILs carrying the epiQTL alleles. This is temptingly close to showing a mechanism. However, more statistical analysis is needed to make this convincing. A large number of infection-inducible genes are up-regulated in the resistant epi-RILs. Is the fraction of these detected in the Hi-C experiment as being linked to the QTL loci greater than would be expected by chance?

We agree with the reviewer that our initial Hi-C analysis remains too preliminary to conclude a contribution of DDM1-dependent heterochromatic interactions to defence gene regulation in the resistant epiRILs. As outlined in our response to comment #3 from the first reviewer and comment #7 from the second reviewer, additional statistical analysis would be needed to support this claim. Sadly, the publicly available dataset from Col-0 and ddm1-2 plants is based on one replicate per genotype, so lacks statistical power to make this analysis more conclusive. Moreover, the analysis would still only be based on the ddm1-2 mutant, and not on the epiRILs that were analysed for genome-wide gene expression. Ultimately, we believe that a fully replicated Hi-C analysis on the Hpa-resistant epiRILs compared to Wt would be needed to provide more conclusive statistical evidence for long-range chromatin interactions with the genes identified in from RNA-seq analysis.

For these reasons, we have omitted the Hi-C analysis from the main Results section of the manuscript. However, the Discussion of our revised manuscript addresses the possibility that long-range chromatin interactions could play a role in the trans-regulation of global defence gene expression by the epiQTLs. We support this hypothesis by referring to the results in the Feng et al., 2014 study, which demonstrate that the ddm1-2 mutation has a profound impact on the long-range heterochromatic interactions with pericentromeric regions, including those that are covered by our epiRIL regions. For details, see our response to comment #3 of reviewer 1.

[Editors' note: further revisions were requested prior to acceptance, as described below.]

Reviewer #1:

The authors have largely responded to my previous concerns. One thing that I think might help this version of the manuscript is to simplify the Discussion and definitions surrounding augmentation and priming. I found myself getting confused about what was being meant. It would help to have discrete explicit definitions being given to allow the reader to fully understand what the authors mean.

The Introduction of our revised manuscript provides a better definition of priming, and how priming relates to augmentation of defence-related gene induction and quantitative resistance:

“Recent evidence has suggested that reduced DNA methylation increases the responsiveness of the plant immune system (Espinas, Saze and Saijo, 2016. […] In some cases, priming of defense-related genes is associated with post-translational histone modifications that mark a more open chromatin structure (Jaskiewicz, Conrath and Peterhänsel, 2011; Luna et al., 2012).”

https://doi.org/10.7554/eLife.40655.035

Article and author information

Author details

  1. Leonardo Furci

    P3 Centre for Plant and Soil Biology, Department of Animal and Plant Sciences, University of Sheffield, Sheffield, United Kingdom
    Contribution
    Data curation, Formal analysis, Investigation, Visualization, Methodology, Writing—original draft, Writing—review and editing, Performed the experiments, Wrote the manuscript
    Contributed equally with
    Ritushree Jain
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-2155-6289
  2. Ritushree Jain

    P3 Centre for Plant and Soil Biology, Department of Animal and Plant Sciences, University of Sheffield, Sheffield, United Kingdom
    Contribution
    Data curation, Software, Formal analysis, Visualization, Methodology, Writing—original draft, Writing—review and editing, Coordinated and performed data analysis of RNA-seq and bisulfite-seq experiments, Wrote the manuscript
    Contributed equally with
    Leonardo Furci
    Competing interests
    No competing interests declared
  3. Joost Stassen

    P3 Centre for Plant and Soil Biology, Department of Animal and Plant Sciences, University of Sheffield, Sheffield, United Kingdom
    Contribution
    Data curation, Software, Formal analysis, Methodology, Analysed Hi-C data and helped in analysing RNA-seq and bisulfite-seq experiments
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-5483-325X
  4. Oliver Berkowitz

    Department of Animal, Plant and Soil Science, ARC Centre of Excellence in Plant Energy Biology, La Trobe University, Melbourne, Australia
    Contribution
    Investigation, Prepared RNA-seq libraries and conducted the sequencing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-7671-6983
  5. James Whelan

    Department of Animal, Plant and Soil Science, ARC Centre of Excellence in Plant Energy Biology, La Trobe University, Melbourne, Australia
    Contribution
    Supervision, Investigation, Prepared RNA-seq libraries and conducted the sequencing
    Competing interests
    No competing interests declared
  6. David Roquis

    1. Department of Plant Sciences, Technical University of Munich, Freising, Germany
    2. Institute for Advanced Study, Technical University of Munich, Garching, Germany
    Contribution
    Software, Formal analysis, Contributed to data analysis of bisulfite-seq experiment
    Competing interests
    No competing interests declared
  7. Victoire Baillet

    Institut de Biologie de l’Ecole Normale Supérieure (IBENS), Ecole Normale Supérieure, Centre National de la Recherche Scientifique (CNRS), Institut National de la Santé et de la Recherche Médicale (INSERM), PSL Université Paris, Paris, France
    Contribution
    Software, Formal analysis, Performed genomic analysis of shared TE insertion events in the epiQTL intervals
    Competing interests
    No competing interests declared
  8. Vincent Colot

    Institut de Biologie de l’Ecole Normale Supérieure (IBENS), Ecole Normale Supérieure, Centre National de la Recherche Scientifique (CNRS), Institut National de la Santé et de la Recherche Médicale (INSERM), PSL Université Paris, Paris, France
    Contribution
    Resources, Software, Formal analysis, Supervision, Performed genomic analysis of shared TE insertion events in the epiQTL intervals
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-6382-1610
  9. Frank Johannes

    1. Department of Plant Sciences, Technical University of Munich, Freising, Germany
    2. Institute for Advanced Study, Technical University of Munich, Garching, Germany
    Contribution
    Software, Formal analysis, Supervision, Contributed to data analysis of bisulfite-seq experiment
    Competing interests
    No competing interests declared
  10. Jurriaan Ton

    P3 Centre for Plant and Soil Biology, Department of Animal and Plant Sciences, University of Sheffield, Sheffield, United Kingdom
    Contribution
    Conceptualization, Resources, Formal analysis, Supervision, Funding acquisition, Validation, Investigation, Visualization, Methodology, Writing—original draft, Project administration, Writing—review and editing, Conceived and supervised the project, Wrote the draft and revised manuscript
    For correspondence
    j.ton@sheffield.ac.uk
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-8512-2802

Funding

Agence Nationale de la Recherche (ANR-09-BLAN-0237EPIMOBILE)

  • Vincent Colot

European Commission Seventh Framework Programme (291763)

  • Frank Johannes

Deutsche Forschungsgemeinschaft (SFB924)

  • Frank Johannes

Leverhulme Trust (RL-2012-042)

  • Jurriaan Ton

H2020 European Research Council (309944)

  • Jurriaan Ton

Biotechnology and Biological Sciences Research Council (BB/L008939/1)

  • Jurriaan Ton

Biotechnology and Biological Sciences Research Council (BB/P006698/1)

  • Jurriaan Ton

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We thank David Pardo and Ana Lopez for technical assistance in the lab. We thank the La Trobe University’s Genomics Platform for technical support. The research was supported by a consolidator grant from the European Research Council (ERC; no. 309944 ‘Prime-A-Plant’) to JT, a Research Leadership Award from the Leverhulme Trust (no. RL-2012–042) to JT and two BBSRC grants to JT (BB/L008939/1 and BB/P006698/1). Work in the group of VC was supported by the Agence Nationale de la Recherche (ANR-09-BLAN-0237 EPIMOBILE). FJ acknowledges support from the Technical University of Munich - Institute for Advanced Study funded by the German Excellent Initiative and the European Seventh Framework Programme under grant agreement no. 291763. FJ is also supported by the SFB/Sonderforschungsbereich924 of the Deutsche Forschungsgemeinschaft (DFG).

Senior Editor

  1. Christian S Hardtke, University of Lausanne, Switzerland

Reviewing Editor

  1. Daniel J Kliebenstein, University of California, Davis, United States

Publication history

  1. Received: July 31, 2018
  2. Accepted: January 3, 2019
  3. Accepted Manuscript published: January 4, 2019 (version 1)
  4. Version of Record published: January 22, 2019 (version 2)

Copyright

© 2019, Furci et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 1,536
    Page views
  • 344
    Downloads
  • 0
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Download citations (links to download the citations from this article in formats compatible with various reference manager tools)

Open citations (links to open the citations from this article in various online reference manager services)