ARTICLE

https://doi.org/10.1038/s41467-021-27320-6 OPEN

Heterochromatin is a quantitative trait associated

with spontaneous epiallele formation
Yinwen Zhang1, Hosung Jang2, Rui Xiao 1, Ioanna Kakoulidou 3, Robert S. Piecyk3, Frank Johannes 3,4 ✉ &
Robert J. Schmitz 2,4 ✉
1234567890():,;

Epialleles are meiotically heritable variations in expression states that are independent from
changes in DNA sequence. Although they are common in plant genomes, their molecular
origins are unknown. Here we show, using mutant and experimental populations, that
epialleles in Arabidopsis thaliana that result from ectopic hypermethylation are due to feed-
back regulation of pathways that primarily function to maintain DNA methylation at het-
erochromatin. Perturbations to maintenance of heterochromatin methylation leads to
feedback regulation of DNA methylation in genes. Using single base resolution methylomes
from epigenetic recombinant inbred lines (epiRIL), we show that epiallelic variation is
abundant in euchromatin, yet, associates with QTL primarily in heterochromatin regions.
Mapping three-dimensional chromatin contacts shows that genes that are hotspots for
ectopic hypermethylation have increases in contact frequencies with regions possessing
H3K9me2. Altogether, these data show that feedback regulation of pathways that have
evolved to maintain heterochromatin silencing leads to the origins of spontaneous hyper-
methylated epialleles.

1 Institute of Bioinformatics, University of Georgia, Athens, GA, USA. 2 Department of Genetics, University of Georgia, Athens, GA, USA. 3 Department of

Plant Sciences, Technical University of Munich, Freising, Germany. 4 Institute for Advanced Study (IAS), Technical University of Munich, Garching, Germany.
✉email: frank@johanneslab.org; schmitz@uga.edu

NATURE COMMUNICATIONS | (2021)12:6958 | https://doi.org/10.1038/s41467-021-27320-6 | www.nature.com/naturecommunications 1

ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-27320-6

G
enetic variation is the primary driver of phenotypic var- clues from previous published studies. For example, the creation
iation, yet, there are well-documented examples of phe- of hybrids between met1 and wild type reveals widespread
notypes arising independent of changes to the DNA redistribution of DNA methylation29. Wild-type chromosomes
sequence. These variants are referred to as epigenetic alleles experienced hypomethylation, which was complemented by
(epialleles) and they can be transmitted across generations, greater methylation on the met1 chromosomes29. Moreover, a
especially in plants1,2. However, the extent and the contribution comparison of the DNA methylomes between 1st generation
of natural epialleles underlying certain traits is less clear3. ddm1 mutants and a 9th generation ddm1 mutants shows
This lack of knowledge is mostly a consequence of the challenges increasing reductions of DNA methylation over generations30.
to studying epialleles, both in their identification and in Unexpectedly, novel epialleles, in the form of ectopic DNA
demonstrating causality. methylation, are abundant in the 9th generation mutants even
The best characterized epialleles causally linked to phenotypes though these plants are essentially wild type in sequence30. These
were discovered in plants and most were associated with changes data point to a model whereby epigenomes are maintained by
in DNA methylation4–9. For example, colorless non-ripening (cnr) some unknown mechanisms that involves some level of feedback
epiallele results from spontaneous formation of DNA methylation regulation based on the genomes ability to sense its overall
in its upstream regulatory region and a decrease in gene epigenomic state.
expression, resulting in tomato fruits with a mealy discolored skin There are multiple examples of feedback regulation of epige-
compared to wild type5. Once formed, these particular epialleles nomic states in plants. For example, expression of the DNA
are stably inherited in subsequent generations, with rare reversion demethylase REPRESSOR OF SILENCING 1 (ROS1) is sensitive to
events observed along the way. The use of whole genome bisulfite genome-wide levels of DNA methylation31. Similarly, the
sequencing (WGBS)10,11 to identify differentially methylated expression of the full-length transcript encoding an H3K9
regions has rapidly advanced our understanding of epiallele fre- demethylase, named INCREASE IN BONSAI METHYLATION 1
quency, stability and molecular nature. For example, WGBS of a (IBM1)32, is sensitive to H3K9me2/DNA methylation levels
population of mutation accumulation lines of Arabidopsis thali- within an intronic repeat33. Additional clues of epigenomic
ana shows that epialleles are rare, enriched in transcribed regions feedback regulation are provided by increasing phenotypic var-
(genes, repeats and transposons) and are relatively stable once iation with each additional generation of selfing in the Arabi-
formed12–14. Importantly, the single base resolution nature of dopsis thaliana mutants met1, ddm1 and ibm127,30,32,34.
WGBS also reveals which DNA methylation pathways are most Evidence for DNA methylation feedback regulation is even
often associated with epiallele formation10,11. observed in natural populations of Arabidopsis thaliana where a
In plants, cytosine DNA methylation is present at CG, CHG (H negative correlation was found between the abundance of CHG
= A, C or T) and CHH sites15. Multiple independent pathways/ methylation, and the frequency of gene body DNA methylated
enzymes coordinately reinforce DNA methylation, which ensures (gbM) genes35. This leads to a model whereby the abundance of
its stability through mitotic and meiotic cell divisions. For heterochromatin, demarcated by H3K9me2 and CHG methyla-
example, METHYLTRANSFERASE 1 (MET1) is recruited to tion, is linked to the origins of epialleles that transition from
hemimethylated CGs during S-phase of DNA replication to unmethylated to gbM genes35. However, the natural variation in
maintain their methylation16,17. Upon exiting S-phase, CHRO- gbM is dependent on a functioning CMT3 pathway, as without
MOMETHYLASE2 and 3 (CMT2/3) are recruited to DNA CMT3, gbM cannot be established and maintained36–38.
associated with histone 3 lysine 9 dimethylation (H3K9me2) to Although evidence for DNA methylation feedback regulation
methylate CWA (W = A or T) and CHG sites, respectively18–21. between chromosomes and its role in epiallele formation is
After the cell exits the cell cycle, the RNA-directed DNA mounting39, there is still a lack of mechanistic understanding of
methylation (RdDM) pathway, which is guided by small RNAs, this process. How is it that unmethylated sequences become
directs DOMAINS REARRANGED METHYLTRANSFERASE 2 newly methylated spontaneously? What provides the trigger and
(DRM2) to methylate cytosines in any sequence context22. Often, why are some regions more susceptible for epiallele formation
these pathways are found acting upon the same sequences and than others? In this study, we follow the fate of epialleles in genes
their coordinated efforts result in the relatively stable DNA that are either converted from unmethylated to methylated or
methylation patterns that are observed between distinct cell types from gbM to a transposon like methylation (teM) pattern in a
and over generational timescales in plants23. variety of experimental mutant lines and populations. We show
Although epialleles are readily identified and characterized that loss of ibm1 leads to ectopic CMT2 and 3 activity in the form
using WGBS to detect losses and gains of DNA methylation, their of CWA and CHG methylation, respectively, and almost exclu-
spontaneous origins are unknown. This has led to use experi- sively at gbM genes compared to unmethylated genes. This
mental approaches to follow that fate of epialleles and to increase ectopic activity also occur, albeit at a lower frequency, when
the frequency of epiallele formation to study their molecular methylation of these gbM genes is erased using a met1 epiRIL
origins and their association with traits. For example, silenced line. This shows the importance of the CMT3-H3K9me2 feedback
epialleles can be induced at a family of MuDR transposons in loop in the establishment of epialleles. This is further evaluated
maize by exposing them to a silencing trigger Muk24. Although using base resolution methylomes of 169 ddm1-derived epiRILs,
MuDR elements are silenced after exposure to a silencing trigger, where the abundance of heterochromatin is a variable trait in
one specific MuDR element reverts over time in the absence of the these lines. Using these data, we discover that variation in the
trigger. Other examples of experimental induced epialleles genome-wide levels of CHG methylation, which serves as a proxy
includes the creation of epigenetic recombinant inbred lines for heterochromatin, is negatively correlated with the abundance
(epiRILs) in Arabidopsis thaliana25,26, whereby RILs are created of epiallele formation, most of which is targeted to gbM genes.
between wild type and a mutant defective in maintenance of Methylation QTL analysis further confirm the importance of
DNA methylation (met1 or DECREASE IN DNA METHYLA- heterochromatin for epiallelic variation, as pericentromeres are
TION 1—ddm1)27,28. The epialleles in the epiRIL populations are hotspots for QTL. Collectively, this study demonstrates that a
mostly due to losses or gains of parental methylation states, positive feedback loop between H3K9me2 and CMT2/3 is a major
however, epialleles also form at novel regions that were not contributing factor to the origins of spontaneous epialleles
present in either parent. The molecular basis for the origins of and that heterochromatin is a quantitative trait that influences
epialleles is an active area of investigation and there are some epiallele formation.

2 NATURE COMMUNICATIONS | (2021)12:6958 | https://doi.org/10.1038/s41467-021-27320-6 | www.nature.com/naturecommunications

NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-27320-6 ARTICLE

a
1.0 WT
Methylation level 1G ibm1
0.8
3G ibm1
0.6
0.4
0.2
0.0
CG CHG CHH CG CHG CHH CG CHG CHH
Total genes gbM UM

b Chr3
c d
gbM UM gbM UM
21,625,000 21,641,655 21,658,700 (%) WT
AT3G58470.1 AT3G58510.1 AT3G58560.1 WT 1G ddm1
40 40

mCG
AT3G58530.1
AT3G58480.1
1G ibm1 2G ddm1
20 3G ibm1 20 9G ddm1
AT3G58490.1 AT3G58540.1
AT3G58460.2 AT3G58520.1 AT3G58570.1 0 0
AT3G58500.1 AT3G58550.1 5
30

mCHG
WT mCG 3
20
1G ibm1 mCHG 10 1
mCHH 0
3G ibm1
6 3
1G ddm1
mCWA

4 2
2G ddm1 2 1
9G ddm1 0
-1kb TSS TTS +1kb -1kb TSS TTS +1kb -1kb TSS TTS +1kb -1kb TSS TTS +1kb

e 3000 2992 12651 f 5108 12612

2265 10000 4000 10000
2000
130 44
1000 5000 2000 5000
31 75 22
0 57 0 2 0 1 0 6
CWA+ CWA+ CWA+ CWA+
CHG+ CHG+ CHG+ CHG+
4000 0
Total gbM
10000 0
Total UM 4000 0 Total gbM 10000 0 Total UM

Fig. 1 Ectopic CHG methylation of gbM genes increases over generations in ibm1 and ddm1 mutants. a The DNA methylation patterns of WT (Col-0), 1G
ibm1, and 3G ibm1 mutants for all genes, gbM and UM genes (N = 33056, 5314 and 12684, biologically independent samples in each group of total, gbM
and UM genes). Box plots show a median center line, the lower and upper hinges are the first and third quartiles. Whiskers represent 1.5x the interquartile
range. b A genome browser view of genes with ectopic non-CG methylation in ibm1 and ddm1 mutants. c DNA methylation patterns of WT (Col-0), 1G
ibm1 and 3G ibm1 mutants for gbM and UM genes. d mCHG-gain (>0.1) and mCWA-gain (>0.1) in 1G ibm1 and enriched in gbM genes. e DNA methylation
patterns of WT (Col-0), 1G ddm1, 2G ddm1, and 9G ddm1 mutants for gbM and UM genes. f Gain of non-CG methylation in 9G ddm1 mutants is enriched
in gbM genes. Source data underlying Figs. 1a, 1c, and 1d are provided as a Source Data file.

Results methylation in the first-generation mutant, which is highly sta-

Ectopic hypermethylation accumulates in gbM genes of tistically significant (Fisher’s exact test, p-value < 0.00001),
methylation mutants. Phenotypic variation of certain mutants, compared to less than 1% (31/12,684) of the UM genes (Fisher’s
namely met1, ddm1, and ibm1, becomes stronger over repeated exact test, p-value = 1, Fig. 1d; Supplementary Data 4). Fur-
generations of self-crossing27,30,32,34. To understand the possible thermore, the levels of ectopic methylation increased over gen-
reasons for this observation, previous studies examined the DNA erations and was also found at CWA sites, which is catalyzed by
methylomes of ibm1 and ddm1 separated by up to eight gen- CMT2 (Fig. 1b–d, Supplementary Figs. 1a, 2a, b). Together, these
erations of single seed descent, respectively30,40–42. These analyses results are consistent with a proposed model whereby loss of
uncovered non-parental epialleles in the form of ectopic hyper- IBM1 activity leads to ectopic H3K9me2 specifically in gene
methylation of CHG specifically in genes, which was expected for bodies of gbM genes, which recruits CMT2 and CMT3 to
ibm143, but was rather unexpected for ddm1. To explore the fate methylate CWA and CHG sites, respectively43,44. The specificity
of these non-parental epialleles further, we measured CG, CHG to gbM versus UM genes is most likely due to ectopic H3K9me2
and CHH methylation in ibm1 mutants selfed for one and three that results from the activity of the histone lysine nine methyl-
generations using previously published data (Fig. 1a). We cate- transferases, SUV4/5/6, which are targeted to methylated DNA
gorized genes into gene body DNA methylated genes (gbM, N = via their SRA domains45. The activity of this positive feedback
5,314) and unmethylated genes (UM, N = 12,684) (Supplemen- loop specifically at gbM genes is consistent with the increase in
tary Data 3). Although it was known that ectopic CHG methy- DNA methylation and with the increased phenotypic variation of
lation occurs in genes, our analysis reveals that this ectopic ibm1 mutants over time.
methylation is essentially exclusive to gbM genes (Fig. 1a–d). The increase in phenotypic variation over generational time is
Over 40% (2,265/5,314) of the gbM genes acquired ectopic CHG much stronger in ddm1 compared to ibm1, although the ibm1

NATURE COMMUNICATIONS | (2021)12:6958 | https://doi.org/10.1038/s41467-021-27320-6 | www.nature.com/naturecommunications 3

ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-27320-6

mutant used has not been propagated for as many generations as Importantly, this ectopic methylation was not due to disruption
ddm1. Furthermore, much of the phenotypic variation in ddm1 is to maintenance of methylation within the 7th intron of IBM1,
mostly due to the loss of maintenance of DNA methylation over which is known to lead to ectopic mCHG in genes (Supple-
generational time. A similar analysis was performed for data from mentary Fig. 4). In total, ectopic CHG methylation was observed
9th generation selfed ddm1, which revealed a reduction of at 1,384 and 1,059 genes in wild-type and ddm1 haplotypes,
CG methylation in gbM genes (Fig. 1d) and showed ectopic CHG respectively (Fig. 3d; Supplementary Data 7). A deeper inspection
and CWA in genes (Fig. 1b, e, Supplementary Figs. 1b, 2c, d). of these genes revealed a significant enrichment for ectopic
However, although there was a greater probability of gbM genes activity specifically at gbM genes compared to teM genes and UM
acquiring ectopic methylation (2.5%, 130/5,314, Fisher’s exact genes (which actually had a strong depletion compared to back-
test, p-value = 8.85e−10), it was much less than what is observed ground expectations—Fig. 3e). Next, we evaluated genes on a
in ibm1 and the ectopic methylation was also observed in case-by-case basis to determine how many of them were sus-
UM genes albeit at a much lower frequency (0.17%, 22/12,684, ceptible to feedback regulation of DNA methylation. For example,
Fisher’s exact test, p-value = 1) (Fig. 1f; Supplementary Data 4). the gene presented in Fig. 3a shows a positive correlation between
Combined, these results show that loss of IBM1 activity the change of CHG methylation and the hypomethylation index
immediately and directly affects ectopic methylation of gbM regardless of haplotype the gene resided within (Fig. 3f). There-
genes almost exclusively, whereas loss of DDM1 activity leads to fore, we performed a correlation analysis between the difference
preferential accumulation of non-parental epialleles via ectopic in CHG methylation between epiRILs and wild type for each gene
methylation with a preference for gbM genes compared to that has evidence for ectopic CHG methylation in any line. These
UM genes. results revealed, yet again, a strong preference for gbM genes
compared to UM and teM genes (Fig. 3g). Collectively, these
results demonstrate that a greater amount of DNA hypomethy-
Abundance of heterochromatin is associated with the fre-
lation, which is associated with a greater reduction to hetero-
quency of spontaneous epialleles. DNA methylation feedback
chromatin methylation, leads to ectopic activity of CMT3 in
regulation between chromosomes is an emerging theme and has
genes, with a clear preference for gbM loci.
major implications for the formation of spontaneous epialleles39.
One of the best examples of this process is nicely demonstrated
using the ddm1 epiRIL population that had variation in the
Heterochromatin is a hotspot for methylation QTLepi. The fact
number of chromosomal regions inherited from the ddm1 parent
that genes in both wild-type and ddm1 haplotypes are almost
versus the wild-type parent. It was observed that hypomethylated
equally affected indicates that they are regulated in trans (by
chromosomes from ddm1 led to ectopic de novo methylation at
distant loci). To test this hypothesis, we used a methylation
regions inherited from wild type and that the amount of hyper-
QTLepi (meQTLepi) approach to explore the potential causal basis
methylation was negatively correlated with the abundance of
for the observed spontaneous hypermethylation in the ddm1
heterochromatin30.
epiRILs. Traditional genetic markers could not be used given the
The initial observations of the ddm1 epiRIL population were
reduced genetic variations in these lines. Instead, we used dif-
examined using microarray technology30, which does not provide
ferentially methylated regions (DMRs) that segregate in a stable,
the opportunity to evaluate methylation in specific contexts. As a
Mendelian fashion in this population (Supplementary Data 5), as
result, WGBS on a subset of three epiRILs was used to further
markers for QTL analysis9. Using the single base resolution
reveal the type of hypermethylation present30. To further expand
methylomes from 169 ddm1 epiRILs we verified, improved and
on these intriguing initial observations, we performed WGBS on
increased coverage of the existing map compared to previous
169 individuals from the ddm1 epiRIL population (Supplemen-
attempts. We scanned each chromosome for associations with
tary Data 1) that are unique from the previously studied set of 123
global mCHG, average mCHG of genes and the number of
lines using tiling arrays46. Using differentially methylated regions
mCHG-gain genes within each epiRIL line (Fig. 4a; Supplemen-
(Supplementary Data 5) we were able to create a haplotype map
tary Data 8 and 9). There were no obvious associations with
of each line that differentiated whether a chromosomal segment
global levels of mCHG, however, the pericentromeric regions of
was inherited from wild type versus ddm1 (Fig. 2a). Each ddm1
multiple chromosomes were associated with the average mCHG
epiRIL was assigned a hypomethylation index based on the
levels of genes as well as the number of genes with ectopic mCHG
amount of DNA inherited from the ddm1 parent and the
(Fig. 4a). To study this further, we performed a comprehensive
haplotype map (Fig. 2b). The higher the hypomethylation index
QTL analysis for each of the 1,595 genes that had evidence of
indicates increased amounts of DNA from the ddm1 parent is
ectopic mCHG. In total, 701 of these genes were associated with
present in the line. Next, we examined ectopic DNA methylation
meQTLepi and the majority (718/975) of the QTL resided in
in gbM genes (Supplementary Fig. 3) and observed an
heterochromatin regions (Fig. 4b). A genome-wide view of these
enrichment of ectopic CHG and to a certain extent CWA
events shows that distant (trans) meQTLepi are especially enri-
methylation, which was most apparent in lines with a very high
ched in the pericentromeric regions (Fig. 4b, c; Supplementary
hypomethylation index (Fig. 2c). Collectively, these data support
Data 10 and 11). Curiously, chromosome 3 was somewhat devoid
that hypomethylated chromosomes from the ddm1 parent led to
of meQTLepi in contrast to other chromosomes, which was
in imbalance in DNA methylation patterns, resulting in ectopic
unexpected given the presence of well-known repeats such as a 5
activity of CMT2 and CMT3 at gbM genes
kb chloroplast insertion, an rDNA and a telomeric repeat. Of the
associations detected, most were linked with a single meQTLepi,
Heterochromatin impacts spontaneous epiallele formation in although numerous loci that gained mCHG in the gene body were
trans. To further evaluate the impact to the disruption to normal found to be associated with 2–3 meQTLepi (Fig. 4d). Collectively,
DNA methylation states, we analyzed genes based on whether these data show that disruption to DNA methylation at hetero-
they were inherited from wild-type versus ddm1 derived haplo- chromatin that was triggered by the initial loss of ddm1 leads to
types. Ectopic activity of CMT3 and CHG methylation was spontaneous epiallele formation at hundreds of loci across the
observed at gbM genes (Fig. 3a) regardless of which haplotype chromosomes. Importantly, the ectopic CHG methylation is
they were present within and this ectopic activity correlated with conditional on the heterochromatic state and does not segregate
the hypomethylation index (Fig. 3b, c; Supplementary Data 6). independently of it, indicating that there is constant feedback.

4 NATURE COMMUNICATIONS | (2021)12:6958 | https://doi.org/10.1038/s41467-021-27320-6 | www.nature.com/naturecommunications

NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-27320-6 ARTICLE

a WT not covered ddm1 b

1 2 3 4 5 ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
169 ddm1 epiRILs

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

0 10 20 30 5 15 0 10 20 5 15 0 10 20 0.0 0.2 0.4 0.6

Chromosome (M) Hypomethylation index
c d
10 Hypo index

# of mCHG−gain genes
●

0.6
8 0.4 300 ●
other mCHH (%)

0.2 200 ●
●
●

●
●
5 0.0 ● ●
●
mCWA (%)
mCHG (%)

● ●

100 ● ● ●
● ●●
●● ●
●
●
● ● ●● ● ●● ●
● ●
● ● ● ●
●● ● ●●● ● ● ● ● ●●●● ●●●●●● ● ● ●● ●
2 ●
●
● ● ●● ●
● ●● ● ● ● ● ●
●
●●● ● ● ●●●● ●●●●
●●●●●●●●●●●●●●●●●●● ●●●●●●●●●●●●● ●●●●●●● ●● ●● ●●●●●●●●●● ● ● ● ●● ●●●●● ●● ● ● ● ●●● ● ●
● ●
●●

0
0 0.2 0.4 0.6
-1kb TSS TTS +1kb -1kb TSS TTS +1kb -1kb TSS TTS +1kb Hypomethylation index

Fig. 2 Hypomethylated chromosomes from the ddm1 parent lead to ectopic DNA methylation in gbM genes. a A haplotype map for 169 ddm1 epiRIL lines.
Each row represents the genome composition of one epiRIL line which was ordered based on hypomethylation index which indicates amounts of DNA inherited
from the ddm1 parent. Yellow bars represent regions from ddm1 parent, blue bars represent regions from wild-type parent and gray bars show regions not
covered by markers. b Hypomethylation index is calculated as an average of the values on 144 DMR markers (methylated WT marker =0, unmethylated ddm1
marker =1). So higher the value represents a more hypomethylated genome. c Metaplot shows non-CG methylation pattern on gbM genes for 169 ddm1 epiRIL
lines. Methylation patterns for each sample is shown by one line that was colored based on its hypomethylation index. d The scatter plot indicates the number
of gbM genes that gain CHG methylation for each sample. Source data underlying Figs. 2a, 2c, and 2d are provided as a Source Data file.

H3K9me2 is important to the formation of spontaneous CG methylation is lost. Therefore, using these lines we can compare
epialleles at unmethylated loci. The ectopic CHG and CWA genes that maintain gbM versus those that were once gbM, but are
methylation at gbM genes observed in ibm1, ddm1 and ddm1 epiRILs now UM due to inheritance via met1.
demonstrates independent paths for how positive feedback regulation This design enables us to test the hypothesis that H3K9
between H3K9me2 and CMT2/CMT3 is established de novo. methylation can establish itself de novo to recruit CMT2/3
However, understanding the origins of spontaneous ectopic hyper- activity to unmethylated genes. To increase the prevalence of
methylation from these experiments is complicated by pre-existing H3K9me2, we crossed ibm1-6 into the met1 epiRIL-12 and
DNA methylation at gbM genes. Previous research has shown that isolated a homozygous ibm1 line. We produced or used publicly
ectopic CHG accumulates in gene bodies in met1, which is devoid of available WGBS data from Col-0, met1-3, ibm1-6, met1 epiRIL-12
CG methylation11. We re-evaluated met1 WGBS data and found and ibm1;met1 epiRIL-12 (Supplementary Data 236,47). Using CG
1,161 genes that accumulate CHG methylation compared to wild methylation patterns, we were able to identify regions within each
type, with a strong statistical enrichment at gbM genes (Supple- genotype that were derived from the met1 versus the wild-type
mentary Fig. 5). However, loss of MET1 leads to wide range of and/or ibm1-6 parents (Fig. 5a and Supplementary Fig. 6). We
disruption to normal nuclear processes that could lead to indirect identified 9.1 Mbs of sequence in ibm1;met1 epiRIL-12 where
effects that explain these observations. Therefore, to understand how 256 gbM genes had lost all CG methylation on chromosome 2
genes can spontaneously transition from unmethylated to methylated and were homozygous for ibm1-6 (Fig. 5b). This met1-derived
in a wild-type genotype, we took advantage of the met1 epiRIL region also possessed 625 genes that were UM. We compared
population26. In the met1 epiRILs, large segments of chromosomes ibm1-induced ectopic DNA methylation of this collection of 881
have lost gbM because they were derived from the met1 parent where genes, all of which have no CG methylation, yet have different

NATURE COMMUNICATIONS | (2021)12:6958 | https://doi.org/10.1038/s41467-021-27320-6 | www.nature.com/naturecommunications 5

ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-27320-6

Fig. 3 Global hypomethylation leads to ectopic activity of CMT3 at genes present within wild-type and ddm1 haplotypes. a Genome browser view of a
representative gene AT3G59750, in which non-CG methylation was observed in epiRILs with either wild-type or ddm1 derived haplotypes. b Correlation
between hypomethylated genome and the number of gbM genes that gain CHG methylation. Genes inherited from WT (blue dots) and ddm1 (yellow dots)
parents were calculated separately for each epiRIL. Gray area shows 95% confidence level interval for predictions (shown by lines) from a loess regression
model. c Correlation between hypomethylated genomes and the average accumulated CHG methylation level on mCHG-gain genes. Gray area shows 95%
confidence level interval for predictions (showing by lines) from a linear regression model. d A union set of total mCHG-gain genes from all epiRIL lines.
It shows the number of mCHG-gain genes from each haplotype type as well as those shared by both haplotypes. e An enrichment analysis using Fisher’s
Exact test shows that mCHG-gain genes that either from ddm1 (left) or WT haplotypes (right) are enriched in gbM genes. Yellow bars show the ratio of
each type of genes with ectopic mCHG. Gray bars show the ratio of each type of genes in comparison with all coding genes. p- value is based on one-sided
test with alternative hypothesis that odds ratio is greater than 1. f Genome-wide association of methylation level changes (difference of value compared to
WT) at an example gene compared to its hypomethylation index. An example of positive correlation is shown for AT3G59750. For each gene, the epiRILs
were grouped based on haplotypes, either WT (blue) or ddm1 (yellow) derived haplotypes. Gray area shows 95% confidence level interval for predictions
(showing by lines) from a linear regression model. g Correlation test in f was applied for all genes. The bar plot shows the number of genes with CHG
methylation level changes strongly correlated with the hypomethylation index (with Pearson correlation coefficient either >0.35 (above the line) or <−0.35
(below the line)). The number of associated genes is calculated separately based on haplotypes (either WT or ddm1 derived haplotypes) and gene type
(gbM, teM and UM). Source data underlying Figs. 3b, 3e, and 3g are provided as a Source Data file.

6 NATURE COMMUNICATIONS | (2021)12:6958 | https://doi.org/10.1038/s41467-021-27320-6 | www.nature.com/naturecommunications

NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-27320-6 ARTICLE

a global mCHG mCHG number

b
Heterochromatin Euchromatin Heterochromatin
ave. mCHG

4 400

Number of methylQTL
600
LOD value/Threshold

3 300

400
2 200

100 200
1

0 0
0
1 2 3 4 5 E H
Chr1 Chr2 Chr3 Chr4 Chr5 Chromosome Trait position

c
1 2 3 4 5
d
10

500

1
Cis
30

400
5 15 0 10 20

Trans

2
Trait position (Mb)

300

Counts
3
200

100
5

4
15

0
8

1 2 3 4
5

the number of associated QTL

16

10 30 5 15 0 10 20 5 15 8 16
QTL marker position (Mb)

Fig. 4 Methylation QTL for genes with ectopic mCHG are especially enriched in the pericentromeric regions. a Standard interval QTL mapping for genic
mCHG related phenotypes. Methylation status (methylated, unmethylated) on 144 stably inherited DMR makers was used for QTL mapping. The results
for three phenotypes, global mCHG level (global mCHG), the number of mCHG-gain genes (mCHG #), averaged mCHG level on genes (ave. mCHG), were
plotted as likelihood-ratio test statistic (LOD score divided by threshold value of 3) against physical location of markers (Mb). b Barplots shows that
methylQTL are especially enriched in the pericentromeric regions in comparison to euchromatin, whereas affected genes are enriched in euchromatin
(right bar plot). E=Euchromatin and H=Heterochromatin. c The scatter plot shows the location of the 701 genes against the location of their corresponding
methylQTL, respectively. d A comprehensive QTL analysis for each of the1595 genes that had evidence of ectoipc mCHG. Significantly associated QTL loci
were identified for 701 genes. For association detected, half of them were linked with a single methylQTL, and others were found to be associated with 2–4
methylQTL. The majority of methylQTL are distant (trans) to associated genes. Source data underlying Fig. 4b–d are provided as a Source Data file.

histories of CG methylation (either gbM or UM - Fig. 5b). No explanation for the preferential accumulation of mCHG in certain
obvious increase in ectopic CHG or CWA methylation was genes is that they are in nuclear neighborhoods/compartments
observed in ibm1;met1 epiRIL-12 using a meta-analysis of these that have higher concentrations of nucleosomes that possess
genes (Fig. 5b). Therefore, we performed a gene specific analysis, H3K9me2. This increased localized pool of H3K9me2 could lead
which identified a significant enrichment (46/256, Fisher’s exact to spontaneous and rare incorrect incorporation of these
test, p-value = 1.81e−16) of ectopic CHG and/or CWA nucleosomes into gene bodies when chromatin is reassembled
methylation in genes that were at one time gbM compared to after DNA replication. Once present in a gene body, CMT3 would
UM genes (12/625, Fisher’s exact test, p-value = 1) (Fig. 5c–e; bind H3K9me2 to methylate associated DNA18,48. Presumably,
Supplementary Data 12). The genes that acquired CHG H3K9me2 would be removed by IBM1 once the gene is tran-
methylation were on average much longer than genes that did scribed, but mCHG would still be present. To test the hypothesis,
not regardless of their historical gbM versus UM status (Fig. 5f). we used Hi-C data49, which incorporates proximity ligation, to
In summary, these results show that spontaneous epialleles can reveal chromatin interactions. We identified all regions that have
form in unmethylated genes and that there is a preference for at least one edge of the chromatin interaction that overlaps with
genes that were gbM in previous generations. The reason for this H3K9me2 (Fig. 6a). We discovered that gbM and teM genes are
preference is unknown, but it is likely in part due to the length enriched, albeit weakly, for contacts with H3K9me2, whereas UM
and transcriptional activity known to be associated with gbM genes were not enriched (Fig. 6b). The distance between
genes. H3K9me2 regions that contacted gbM and UM genes were
not significantly different from one another, yet the teM
genes were significantly further away (Fig. 6c). What did distin-
Hotspots for epiallele formation have increased chromatin guish gbM genes that were hotspots for spontaneous epiallele
contacts with H3K9me2 regions. Why certain genes are hotspots formation (ectopic mCHG) was that they had a greater contact
for epiallele formation versus others is unknown. One possible frequency with H3K9me2 regions compared to UM genes

NATURE COMMUNICATIONS | (2021)12:6958 | https://doi.org/10.1038/s41467-021-27320-6 | www.nature.com/naturecommunications 7

ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-27320-6

a b gbM genes; UM genes;

Chr2 met1-3 derived region met1-3 derived region
0 100% 60%
mCG mCG ibm1-6
40%
met1-3
Chr (Mb)

20%

Methylation level
mCG
0 Col-0

40% met1
mCHG mCHG
19 0 epiRIL-12
l-0 m1 et1 12 12 20%
C o ib m iRIL- iRIL-
ep 1 ep ibm1;met1
t1 t
me ; me 0 epiRIL-12
m 1
i b -1kb TSS TTS +1kb -1kb TSS TTS +1kb

c gbM gene; AT2G12550 UM gene; AT2G45200

Col-0

ibm1-6
mCG
met1-3 mCHG
mCHH
met1 epiRIL-12

ibm1; met1 epiRIL-12

p=0.0011
216
d e 600
613 f 15
p=0.0081
200

150 10
400

Gene length (kb)

100
200 5
50
32
11 3 12
0 0
mCWA mCHG 0
mCHG
+

;m +

-
UM HG

G
G

CH
UM CH

gb CH
Total UM

C
Total gbM

m
;m

m
600 200
200 0 M
M
gb

met1-3 derived region; gbM genes met1-3 derived region; UM genes

Fig. 5 Mutations in ibm1 lead to ectopic methylation preferentially in genes that previously possessed gbM. a A schematic showing regions that
inherited mCG from the met1 epiRIL-12 parent or the ibm1-6 parent on Chromosome 2 (Col-0, ibm1-6 and met1-3 are shown as controls). Purple to black
colors show the mCG percentage of genes that possess gbM in Col-0 from 100% to 0%, respectively. Red arrows indicate a region in ibm1; met1 epiRIL-12
derived from the original met1-3 parent. b DNA methylation plot of gbM and UM genes including 1kb up and downstream for indicated genotypes and
methylation contexts. Depicted regions shown are derived from met1-3. TSS: Transcription Start Site. TTS: Transcription Termination Site. c Genome
browser examples of gbM and UM genes from met1-3 derived region of the met1 epiRIL-12 genotypes. d and e The number of gbM and UM genes that
gained mCHG, mCWA or non-mCWA for in the met1-3 derived region of ibm1; met1 epiRIL-12. Vertical bar plot shows count of overlapping genes between
different methylation contexts. Horizontal bar plot shows count of specific type of genes. f Gene length distribution of genes that gained mCHG in gbM and
UM genes derived from the met1-3 region of ibm1; met1 epiRIL-12. N = 43, 12, 219, and 613 biologically independent samples for four groups from left to
right. Gain is depicted by ‘+’ and no gain is depicted by ‘-’. p-value is generated by two-sided student t-test. Box plots show a median center line, the lower
and upper hinges are the first and third quartiles. Whiskers represent 1.5x the interquartile range. Source data underlying Figs. 5b and 5f are provided as a
Source Data file.

(Fig. 6d). This result was further confirmed using Hi-C data from possible mechanism by which certain loci are susceptible for
ddm150. Although global hypomethylation in ddm1 leads to a spontaneous epiallele formation could be through association
significant reduction of contact frequency (most obvious for teM with H3K9me2 regions of the genome in three-dimensional
genes) compared to wild-type Col-0, gbM genes still had a greater space.
contact frequency with H3K9me2 regions compared to UM
genes. Even though ectopic non-CG methylation is rare in the
first generation of ddm1 being a mutant, genes with ectopic Discussion
mCHG did show a greater contact frequency with H3K9me2 One possible reason for the prevalence of epialleles in plants is
regions (Supplementary Fig. 8). These results show that one that there is no comprehensive erasure of DNA methylation from

8 NATURE COMMUNICATIONS | (2021)12:6958 | https://doi.org/10.1038/s41467-021-27320-6 | www.nature.com/naturecommunications

NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-27320-6 ARTICLE

H3K9 methyl binding protein SAWADEE HOMEDOMAIN

a b Contact with H3K9me2
Background HOMOLG159 and the DNA methylation readers SUVH2/9
gbM gene UM gene teM gene 60% (Su(var)3-9 homolog) recruit the RdDM pathway to target
sequences for DNA methylation reinforcement45. There are
40% p=4.6e-05
H3K9me2 peak multiple examples in the literature that show the effectiveness of
p=9.36e-54
20% these feedback loops at re-establishing patterns of DNA
c methylation23. For example, multiple proteins have recently been
0%
d used as triggers in epigenome editing whereby the protein is

# of contacts with H3K9me2 peaks

distance to H3K9me2 peaks (kb)

gbM teM UM
targeted to the FWA promoter to establish a DNA
p=5.01e-08 methylation45,60. Once a positive feedback loop is established the
p=0.00191 p=1.549e-05 transgene triggers can be segregated away as they are no longer
10.0 p=1.88e-26
400 p=5.55e -26 needed for maintenance45. Another recent study illustrated how
p=4.88e-53 p=1.16e-53 RNAi-independent pathway(s) effectively re-establishes non-CG
300 7.5
methylation at transposons that had lost all non-CG methylation
200 n=332
5.0 n=332
n=281 due to loss of H3K9me2 and CMT2/3 activity61. Because CG
n=875 n=2424
n=2424 n=875 n=594 methylation was still present at these transposons, it functioned to
100
2.5 recruit SUV4/5/6 to re-establish H3K9me2 and non-CG methy-
0 lation. It is well known how these pathways maintain DNA
gbM teM UM gbM teM UM gbM gbM methylation, but how feedback regulation is established de novo
;m ;m
CH CH at new regions in the genome is challenging to study given their
G+ G-
spontaneous nature.
Fig. 6 Loci that are susceptible to spontaneous epiallele formation have a One possibility that is consistent with the spontaneous nature
greater contact frequency with H3K9me2 region. a Hi-C data was used of epialleles is that the feedback regulation important to main-
for identifying significant three-dimensional contacts between genomic tenance of DNA methylation have a low rate of off-targeting
regions. The schematic plot depicts the interaction from regions possessing activity57. Although their main function is to target repeats and
H3K9me2 and genes. b Enrichment test for genes that have contacts with transposons for silencing, biochemical features of many of these
H3K9me2. Yellow bars show the proportion of each type of genes that components, such as domains that recognize methylated DNA
contacts with H3K9me2. Grey bars show the proportion of each type genes (SRA, MBD) and histones (CHROMO, BAH, SAWADEE) leads
in comparison with all coding genes. Fisher’s exact test was used for to improper establishment of these pathways at unintended
enrichment test. p-value is based on one-sided test with alternative regions in the genome57. Given the function of IBM1, a histone
hypothesis that odds ratio is greater than 1. c Distribution of distance demethylase that removes H3K9me2 from PolII-dependent
between three types of genes (gbM, UM, teM) and H3K9me2 regions. transcribed regions of the genome, it seems likely that plant
Two-sided Wilcoxon tests was used for pairwise comparison of each two epigenomes have evolved mechanisms to reduce off-targeting
groups without adjustment for multiple comparisons. d Distribution of the activity32,43,44,62. In this particular case, pre-existing CG methy-
number of contacts between three types of genes (gbM, UM, teM) and lation at gbM genes serves as a substrate for the SRA domains in
H3K9me2 regions. gbM genes were further classified into a group with SUVH4/5/6 to catalyze H3K9me258, which is counteracted by the
ectopic mCHG (gbM; mCHG+) and a group without ectopic mCHG (gbM; activities of IBM1.
mCHG-) based on ddm1 mutant line data. Two-sided Wilcoxon tests was Our analysis of the ibm1;met1 epiRIL-12, shows that the
used for pairwise comparison of each two groups without adjustment for absence of IBM1 does result in ectopic CMT3 activity at com-
multiple comparisons. Box plots show a median center line, the lower and pletely unmethylated genes, which suggests that H3K9me2 is
upper hinges are the first and third quartiles. Whiskers represent 1.5x the likely a trigger for inducing spontaneous epialleles. Evidence that
interquartile range. Source data underlying Fig. 6b–d are provided as a H3K9me2 can seed de novo methylation independent from
Source Data file. RdDM and pre-existing DNA methylation has also been observed
at the BONSAI locus in A. thaliana42. Therefore, understanding
H3K9me2 dynamics and why certain genes are hotspots for
one generation to the next like there is in mammalian H3K9me2 activity compared to others will be essential to
genomes27,51–53. Instead, extensive reinforcement of DNA understanding the origin of epialleles. The ectopic hypermethy-
methylation occurs during gamete production54, with the lated epialleles identified in the genotypes used in this study show
exception of rare events important for genomic imprinting. This that there is a non-random process that makes certain loci more
feature of flowering plants enables the occurrence of rare spon- susceptible to spontaneous epiallele formation (Supplementary
taneous epialleles in a single cell to propagate in subsequent cell Fig. 7). This has led us to speculate that gbM genes in angios-
divisions and in some cases to become a major allele in the next perms arise from off-targeting activity of positive feedback loops
generation if it is present in a cell that leads to the production of established between CMT3-H3K9me2 and that gbM is evolutio-
gametes. Therefore, once plant epialleles form, they can be stu- narily neutral. These genes are susceptible to this pathway
died to understand their stability over generations, interaction because of their inherent functions. These genes are generally
with alleles with different epigenomic states and association with ‘housekeeping’ genes that are long, moderately expressed in all
phenotypes. cells and typically not regulated by development/environment.
Another likely explanation that is consistent with the epige- Transcription at these loci is likely a prerequisite to becoming a
netic nature of epialleles is the involvement of positive feedback gbM gene, transcription could lead to incorporation of H3K9me2
loops. These positive feedback loops are fundamental to main- nucleosomes at a low rate. Typically, this H3K9me2 is removed
tenance of DNA methylation in plants. Hemimethylated CGs are by IBM1, but over evolutionary time mistakes occur which lead to
recognized by VARIANT IN METHYLATION 155, which ectopic activity of CMT3 leading to CHG methylation and
recruits MET1 to maintain CG methylation. CMT2 and CMT3 eventually CG methylation. How CHG methylation transitions to
are directed to methylate DNA through their binding to CG methylation to provide transgenerationally stability of epial-
regions with H3K9me2 and H3K9 methyltransferases bind leles is currently unknown and was not addressed in this study.
methylated DNA and methylate H3K918–20,48,56–58. Similarly, the There are at least two possible mechanisms by which this

NATURE COMMUNICATIONS | (2021)12:6958 | https://doi.org/10.1038/s41467-021-27320-6 | www.nature.com/naturecommunications 9

ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-27320-6

occurs63. One includes the activity of the RdDM pathway, which 14 °C with reduced light conditions. Following this, they were grown under long
is recruited to regions possessing H3K9me2 via SHH enabling day conditions (16h light, 8h dark) at 20/18 °C, 60–75% humidity and 180-240 μΕ
light intensity. The pots were watered with 55, 30 and 20 mL water the 2, 7, and
DRM2 to methylate CGs. This could be especially prevalent 8 day after sowing (DAS) and then every other day with 55 mL.
during CHH DNA methylation reinforcement that occurs in the
embryo and after fertilization52,54,64–67. The second possibility Library construction. All epiRIL lines were harvested at 27 DAS in a time frame of
includes recruitment of VIM1 to methylated DNA via its SRA 3 h. Flowering stems and roots were removed and all 18 individual plants from
domain (in this case CHG methylation)68, which could lead to each epiRIL line were pooled in 50 mL tubes. They were immediately frozen in
rare de novo methylation of CGs by MET1. Regardless of the liquid nitrogen and stored to −80 °C until processing. Genomic DNA was
exact mechanism of transitioning CHG to CG methylation, the extracted from each pooled sample using the DNAeasy plant mini kit from Qiagen.
169 epiRIL lines with at least 1 μg of DNA were submitted to the Beijing Genome
incorporation of these feedback loops would be prevented at Institute where they were prepared for WGBS libraries. Sequencing was performed
many developmental/environmental regulated genes as they are on an Illumina HiSeq X ten instrument. Clean raw paired-end files were obtained
targeted by the Polycomb Repression Complex and H3K27me3 to from BGI and used for downstream analysis. The ibm1;met1 epiRIL-12 library was
precisely limit their expression69. H3K27me3 is incompatible prepared following the MethylC-seq protocol74. Briefly, genomic DNA was soni-
cated to 200 bp using a Covaris S-series focused ultrasonicator, and end-repaired
with maintenance of DNA methylation in A. thaliana, as it using End-It DNA end-repair kit (Epicentre). End-repaired DNA was subjected to
colocalizes with H2A.Z70. Given H3K9me2 associates with A-tailing using Klenow 3′–5′ exo− (NEB) and ligated to methylated adapters using
H2A.W71, the presence of H2A.Z would prevent the CMT3- T4 DNA ligase (NEB). Ligated DNA was subsequently bisulfite converted using the
H3K9me2 feedback loop from establishing gbM at many genes in EZ DNA methylation-Gold kit as per the manufacturer’s instructions and ampli-
fied using KAPA HiFi uracil + Readymix Polymerase.
the genome.
We hypothesize that spontaneous epialleles form as a bypro-
Methylome mapping. The WGBS data from the ddm1-epiRILs was analyzed using
duct of enzymatic activities that are dedicated to the maintenance the Methylstar v1.4 pipeline75. Region-level methylation calls were obtained with
of heterochromatin and that variation in heterochromatin Methimpute v1.16.0 (200 bp bins, step size = 50 bps, at least 10 cytosines per
abundance and methylation influences the rate of spontaneous bin)76. We used Methimpute’s 2-state Hidden Markov Model to classify a given
epiallele formation35. If genes that maintain heterochromatin region as either homozygous methylated or homozygous unmethylated. These state
calls were used downstream for the construction of an augmented linkage map in
sometimes act on inappropriate targets, causing epialleles, natural the epiRIL panel. Base-resolution methylome analysis of the ddm1-epiRILs, ibm1
selection cannot reduce epialleles by reducing DNA methylation and met1 epiRIL12 data used in this study were all processed by Methylpy v1.3 as
because that would lead to loss of maintenance of hetero- described in77. Quality filtering and adapter trimming were performed using
chromatin and result in genome instability. In this way, the cutadapt v1.9.dev178. Qualified reads were aligned to the A. thaliana TAIR10
evolution of epialleles is similar to the evolution of chromosome reference genome79 (downloaded from https://phytozome.jgi.doe.gov) using bowtie
2.2.480. Only uniquely aligned and nonclonal reads were retained. Chloroplast
rearrangements. In that case, molecular recombination, is favored DNA (which is fully unmethylated) was used as a control to calculate the sodium
to generate gamete diversity (or for some other reason) that has bisulfite reaction non-conversion rate of unmodified cytosines. A binomial test was
an unintended, deleterious consequences, namely the production used to determine the methylation status of cytosines with a minimum coverage of
of genome rearrangements by ectopic recombination. three reads.
Numerous major questions remain, but the most important
one is how is H3K9me2 initially incorporated into an unmethy- Data acquisition. WGBS data of ddm1 and ibm1 mutants and their self-crossing
offspring used in this analysis were obtained from published datasets30. WGBS data
lated region de novo. Transcription coupled incorporation of of met1-3, ibm1-6 and met1 epiRIL-12 were obtained from previously published
H3K9me2 nucleosomes is one possibility given the role of datasets36,47. Hi-C data and H3K9me2 ChIP-seq reads used in this study were
nucleosome eviction and reincorporation during transcription72. obtained from a previously published datasets49,81.
It’s also possible that H3K9me2 is mis-incorporated to certain
regions of the genome during the establishment of chromatin Gene body methylation status classification. To explore the fate of ectopic
upon DNA replication. In this study, we used Hi-C data to show methylation of genes that have different DNA methylation states, we categorized
genes into one of three classes including gbM, teM and UM based on CG, CHG
the spontaneous epialleles we identified interact with H3K9me2 and CHH methylation in wild-type Col-0. The total number of cytosines and the
regions at greater frequencies that unmethylated genes in the methylated cytosines were counted for cytosines in each context (CG, CHG, and
genome. This result supports that the hypothesis that there are CHH) for the coding sequences (CDS) of the primary transcript for each gene. The
sub-nuclear compartments that could have different concentra- percentage of methylated sites for each sequence context in all coding regions were
used as the background probability of having methylation on a single site. Given a
tions of H3K9me2 nucleosome pools and that spontaneous background probability and the total number of cytosines and methylated cyto-
epialleles are more likely to occur in pools with high concentra- sines, a p-value was calculated using a binomial distribution to show the cumulative
tions of H3K9me2 nucleosomes. However, the current evidence probability of having a higher number of methylated cytosines on a given gene82.
to support these conclusions are premature to make stronger Then a q-value was calculated by adjusting p- values by Benjamin–Hochberg FDR
to control the false discovery rate.
conclusions. Regardless, the involvement of CMT3-H3K9me2 Genes were classified as gbM if they had reads mapping to at least 20 CG sites
feedback regulation is a major factor in the origins of spontaneous and had a q-value < 0.05 for mCG and a q-value > 0.05 for mCHG and mCHH.
epialleles and future studies will be required to test these pro- Genes were classified as mCHG if they had reads mapping to at least 20 CHGs, a
posed models. mCHG q-value < 0.05, and a mCHH q-value > 0.05. As mCG is commonly
associated with mCHG, the q-value for mCG could be significant or insignificant in
mCHG genes. Genes were classified as mCHH if they had reads mapping to at least
Methods 20 mCHH sites and a mCHH q-value < 0.05. q-values for mCG and mCHG could
be anything as both types of methylation are associated with mCHH. mCHG and
Plant material. The only new plant material generated for this study was the create
mCHH genes were collectively referred to as teM genes. Genes were classified as
of ibm1;met1 epiRIL-12. These individuals were produced by crossing ibm1-6
unmethylated (UM) if they had reads mapping to at least 20 mCHH sites and had a
(SALK_006042)73 to a met1 epiRIL-1226 and isolating homozygous ibm1-6 lines in
q-value > 0.05 for all sequence contexts. To make sure the selected UM genes are
the F2. The complete set of ddm1-epiRILs from Johannes et al. was obtained from
truly unmethylated, we further remove genes with more than 2 symmetric mCG
the Versailles Arabidopsis Stock center of INRA (http://publiclines.versailles.
from UM genes list. Only genes that fit the definition of gbM, teM and UM above
inra.fr/)25. All epiRIL lines were propagated in a greenhouse at the Leibniz Institute
were used for the downstream analysis in this study.
of Plant Genetics and Crop Plant Research (IPK). The plants were grown in single
seed pots and at a later developmental stage six siliques per plant were left to dry.
Seeds for a selection of 169 epiRILs were sown in the IPK fully automatic phe- Determination of DNA methylation patterns at genes. For ibm1, ddm1 and
notyping facility for small plants (Junker et al., 2015). Three independent experi- ddm1 epiRIL lines, DNA methylation patterns at gbM and UM genes were
ments in three consecutively months were performed. In each of the 3 cultivation explored. Each gene was divided into 20 windows. Additionally, regions 1000 bp
experiments, 6 individual plants from each line were grown in 2 separate trays upstream and downstream were each divided into 20 50-bp windows. Methylation
completely randomized in the chamber. The sown pots were firstly placed for levels (total methylated reads divided by total reads mapped to cytosine sites in a
3 days at 4 °C in darkness and then the plants were acclimated for 2 days under 16/ window) were calculated for each window77, and a mean value for each window

10 NATURE COMMUNICATIONS | (2021)12:6958 | https://doi.org/10.1038/s41467-021-27320-6 | www.nature.com/naturecommunications

NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-27320-6 ARTICLE

was averaged over the methylation level of the same window from genes in the enrichment of observed from the expected contact counts85. Raw reads of H3 and
same group (gbM/UM). The mean methylation levels and their corresponding H3K9me2 ChIP-seq were trimmed with Trim Galore v0.6.5 (https://github.com/
window number was used to generated metaplots using the R package ggplot2. FelixKrueger/TrimGalore) with default parameters. The remaining reads were
aligned to the A. thaliana TAIR10 reference genome using Bowtie2 v2.3.5.1 with
default setting. Aligned reads were sorted using SAMtools v1.1086. Then,
Locating the met1-derived regions in ibm1;met1 epiRIL-12. The met1 epiRIL-12
unmapped reads, duplicated reads and multiple mapped reads were filtered out
line used in this study was carefully selected to ensure the IBM1 locus was derived
using Sambamb v0.7.187. Lastly, the remaining uniquely mapped reads were used
from the Col-0 parent instead of the met1 parent. The met1 epiRIL-12 derived
in MACS2 for peak calling with parameters ‘-g 1.35e+8–broad’88. H3 ChIP-seq
regions of ibm1;met1 epiRIL-12 was identified by comparing CG methylation levels
was used as a control for identification of H3K9me2 enriched regions. Significant
on the exons of gbM genes between the epiRIL line and its parents (ibm1-6 and
contacts that had at least one edge that overlaps with a H3K9me2 enriched region
met1 epiRIL-12). Each chromosome sequence was divided into continuous bins and
and another edge that overlapped with a gbM (n=875), teM (n=332) or UM
each bin included at least 10 exons. A bin with an average mCG level that was
(n=2424) gene was selected for downstream analysis. A Fisher’s Exact test was used
greater than 25% reduced in ibm1;met1 epiRIL-12, as compared to that in ibm1-6
to evaluate the enrichment of genes that have contacts with H3K9me2 peaks in the
was used to differentiated the ibm1-6 and met1 epiRIL-12 derived regions. All
three categories of genes (gbM, UM, teM). The distribution of distance and the
downstream analyses of this line used the same criteria described in ‘Gene body
number of contacts between three categories of genes and H3K9me2 peaks were
methylation status classification’ section of the Methods for defining gbM and UM
compared using a Two-sample Wilcoxon test. The selected gbM genes (n=875)
genes. Regions that possessed met1 epiRIL-12 derived regions in the ibm1-6 mutant
were further classified into an ectopic mCHG-gain group (gbM; mCHG+) and a
background were used to evaluate the consequence of mutant ibm1-6 on genes that
group without ectopic mCHG (gbM; mCHG-) based on the methylome of ddm1
had lost gbM due to previous loss of met1. Ectopic CHG or CWA methylation on
mutant. The distribution of the number of contacts between these two groups and
gbM and UM genes required a minimum of three sites with an average 10% higher
H3K9me2 peaks were also compared using a Two-sample Wilcoxon test.
methylation level when compared to the wild-type Col-0 parent.
We also used previously published Hi-C data to identify chromatin contacts
between two genomic regions in nonadjacent locations along the genome for both
Construction of an augmented epiRIL linkage map. Previously, we used tiling- wild-type Col-0 and ddm1 mutant50. Significant contacts were obtained from Hi-C
array data from 123 epiRILs and their ddm1-2 and Col-wt founder lines and data using HiC-Pro and Fit-Hi-C following the same pipeline and parameters as
identified differentially methylated regions (DMRs) in the founders that were stably the previous analysis (see the first paragraph in ‘Characterizing chromatin
inherited in the epiRILs through at least 10 rounds of meiosis46. Using these ultra- interactions between genes and H3K9me2’ section of the Methods) in Fig. 6.
stable DMRs as physical markers we were able to construct a linkage map involving Significant contacts that had at least one edge that overlaps with a H3K9me2
126 DMRs in this isogenic experimental system46, and later used this map for enriched region and another edge that overlapped with a gbM, teM or UM gene
QTLepi analysis9. The 169 epiRIL methylomes from the present study were were selected. The number of the three types of genes (gbM/teM/UM) used in the
employed to augment the existing linkage map, with the goal to improve marker analysis were shown separately for both wild-type and ddm1 mutant samples in
spacing and mapping resolution. Among the 169 epiRILs measured in the present Supplementary Fig. 8. The gbM genes (n=1775) selected from ddm1 Hi-C data
work, only 37 overlapped with the 123 epiRILs used in the previous our studies9,46. were further classified into an ectopic mCHG-gain group (mCHG+) and a group
We used the tiling-array-derived linkage map of 126 DMRs as a starting point. For without ectopic mCHG (mCHG-) based on methylome of ddm1 1G mutant. The
the 37 overlapping epiRILs, the WGBS derived methylation state calls at these 126 distribution of the number of contacts between these two groups and H3K9me2
DMR position were consistent with the previous tiling-array calls, indicating that peaks were compared by Two-sample Wilcoxon tests.
these DMRs are robust. For a given chromosome, we stepped in new DMRs into
the existing linkage map. DMRs that correlated significantly across chromosomes
Box plots. All box plots presented show a median center line with an upper and
with any of the core DMRs were rejected. For map cleaning we employed Rqtl’s
lower quartile. Whiskers represent 1.5x the interquartile range.
tutorial on map construction83. For DMRs with identical cM positions we selected
those that best optimized overall marker spacing in terms of base pair coverage.
DMRs were kept or rejected by manual inspection in conjunction with the like- Reporting summary. Further information on research design is available in the Nature
lihood score from the drop.one function in Rqtl. The final linkage map contains Research Reporting Summary linked to this article.
144 well-spaced markers, which provide improved coverage in chromosome arms.
Data availability
Methylation analysis for ddm1 epiRILs. For methylation analysis in ddm1 Previously published sequencing data used in this study are available listed in
epiRILs, we needed to determine whether a chromosomal segment was inherited Supplementary Data 2. WGBS sequencing data produced from this study have been
from the wild-type or ddm1 parent for each line. The methylation levels at deposited in the NCBI GEO database under GSE171157 and GSE171414. Source data are
144 selected DMR markers was calculated for each line. DMRs were classified as an provided with this paper.
unmethylated ddm1 marker if their methylation level was less than 0.5. Otherwise,
DMRs were classified as methylated WT makers if their methylation level was
greater or equal to 0.5. Then, the haplotype map in Fig. 2a was generated based on Code availability
the location of 144 DMRs and their methylation status (methylated, unmethylated) Code used in processing and analysis of these data can be found at GitHub [https://
for each line using the R package ggplot2. We also used a measurement of github.com/schmitzlab/Heterochromatin-quantitative-epiallele-formation]89.
hypomethylation index to determine the amount of DNA inherited from the ddm1
parent. The hypomethylation index was calculated as an average of the values at
144 DMR markers (methylated WT marker =0, unmethylated ddm1 marker =1). Received: 30 April 2021; Accepted: 15 November 2021;
A higher value represents a more hypomethylated genome. To further evaluate the
impact to the disruption to normal DNA methylation states, we classified genes
into wild-type versus ddm1 parent-derived regions based on the haplotype map for
each line. Then, we explored how ectopic CHG methylation was distributed at
genes from WT versus ddm1 parents, respectively. A Fisher’s Exact test was used to References
evaluate the enrichment of mCHG-gain genes in three categories of genes (gbM, 1. Heard, E. & Martienssen, R. A. Transgenerational epigenetic inheritance:
UM, teM). A Pearson correlation coefficient was used to evaluate how CHG myths and mechanisms. Cell 157, 95–109 (2014).
methylation changes against hypomethylation index for each gene. 2. Richards, E. J. Inherited epigenetic variation–revisiting soft inheritance. Nat.
Rev. Genet. 7, 395–401 (2006).
QTL mapping analysis. Using the binary methylation status (M/U) at 144 DMR 3. Springer, N. M. & Schmitz, R. J. Exploiting induced and natural epigenetic
markers of each epiRIL line as genotypic data together with the genic mCHG related variation for crop improvement. Nat. Rev. Genet. 18, 563–575 (2017).
phenotypes, including three global genic mCHG related traits and genic mCHG levels 4. Cubas, P., Vincent, C. & Coen, E. An epigenetic mutation responsible for
on each of the 1595 mCHG-gain genes, we performed interval mapping with Rqtl. The natural variation in floral symmetry. Nature 401, 157–161 (1999).
customized R scripts used for the QTL mapping analysis can be obtained from the 5. Manning, K. et al. A naturally occurring epigenetic mutation in a gene
GitHub repository provided in the Code Availability section9. encoding an SBP-box transcription factor inhibits tomato fruit ripening. Nat.
Genet. 38, 948–952 (2006).
Characterizing chromatin interactions between genes and H3K9me2. We used 6. Ong-Abdullah, M. et al. Loss of Karma transposon methylation underlies the
previously published Hi-C data49 to identify chromatin contacts between two mantled somaclonal variant of oil palm. Nature 525, 533–537 (2015).
genomic regions in nonadjacent locations along the genome. The HiC-Pro v2.11.4 7. Soppe, W. J. et al. The late flowering phenotype of fwa mutants is caused by
pipeline was used to process Hi-C data, from raw reads to normalized contact gain-of-function epigenetic alleles of a homeodomain gene. Mol. Cell 6,
matrices for selected windows using the parameter ‘–binsize 2000’, since 2 kb is an 791–802 (2000).
approximate value of the average gene length of A. thaliana84. The contact matrices 8. Henderson, I. & Jacobsen, S. Tandem repeats upstream of the Arabidopsis
were then transformed to Fit-Hi-C v.2.0.7 readable input files with hicpro2fithic.py. endogene SDC recruit non-CG DNA methylation and initiate siRNA
Then, significant contacts were identified using Fit-Hi-C by assessing the spreading. Genes Dev. 22, 1597–1606 (2008).

NATURE COMMUNICATIONS | (2021)12:6958 | https://doi.org/10.1038/s41467-021-27320-6 | www.nature.com/naturecommunications 11

ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-27320-6

9. Cortijo, S. et al. Mapping the epigenetic basis of complex traits. Science 343, 40. Zabet, N. R., Catoni, M., Prischi, F. & Paszkowski, J. Cytosine methylation at
1145–1148 (2014). CpCpG sites triggers accumulation of non-CpG methylation in gene bodies.
10. Cokus, S. J. et al. Shotgun bisulphite sequencing of the Arabidopsis genome Nucleic Acids Res. 45, 3777–3784 (2017).
reveals DNA methylation patterning. Nature 452, 215–219 (2008). 41. Inagaki, S. et al. Autocatalytic differentiation of epigenetic modifications
11. Lister, R. et al. Highly integrated single-base resolution maps of the epigenome within the Arabidopsis genome. EMBO J. 29, 3496–3506 (2010).
in Arabidopsis. Cell 133, 523–536 (2008). 42. Sasaki, T., Kobayashi, A., Saze, H. & Kakutani, T. RNAi-independent de novo
12. Becker, C. et al. Spontaneous epigenetic variation in the Arabidopsis thaliana DNA methylation revealed in Arabidopsis mutants of chromatin remodeling
methylome. Nature 480, 245–249 (2011). gene DDM1. Plant J. 70, 750–758 (2012).
13. Schmitz, R. J. et al. Transgenerational epigenetic instability is a source of novel 43. Miura, A. et al. An Arabidopsis jmjC domain protein protects transcribed
methylation variants. Science 334, 369–373 (2011). genes from DNA methylation at CHG sites. EMBO J. 28, 1078–1086 (2009).
14. Hofmeister, B. T., Lee, K., Rohr, N. A., Hall, D. W. & Schmitz, R. J. Stable 44. Inagaki, S. & Kakutani, T. What triggers differential DNA methylation of
inheritance of DNA methylation allows creation of epigenotype maps and the genes and TEs: contribution of body methylation? Cold Spring Harb. Symp.
study of epiallele inheritance patterns in the absence of genetic variation. Quant. Biol. 77, 155–160 (2012).
Genome Biol. 18, 155 (2017). 45. Johnson, L. M. et al. SRA- and SET-domain-containing proteins link RNA
15. Law, J. A. & Jacobsen, S. E. Establishing, maintaining and modifying DNA polymerase V occupancy to DNA methylation. Nature 507, 124–128 (2014).
methylation patterns in plants and animals. Nat. Rev. Genet 11, 204–220 46. Colome-Tatche, M. et al. Features of the Arabidopsis recombination landscape
(2010). resulting from the combined loss of sequence variation and DNA methylation.
16. Bostick, M. et al. UHRF1 plays a role in maintaining DNA methylation in Proc. Natl Acad. Sci. USA 109, 16240–16245 (2012).
mammalian cells. Science 317, 1760–1764 (2007). 47. Stroud, H., Greenberg, M. V., Feng, S., Bernatavichute, Y. V. & Jacobsen, S. E.
17. Borges, F. et al. Loss of Small-RNA-Directed DNA Methylation in the Plant Comprehensive analysis of silencing mutants reveals complex regulation of
Cell Cycle Promotes Germline Reprogramming and Somaclonal Variation. the Arabidopsis methylome. Cell 152, 352–364 (2013).
Curr. Biol. 31, 591–600 e594 (2021). 48. Du, J. et al. Dual binding of chromomethylase domains to H3K9me2-containing
18. Stroud, H. et al. Non-CG methylation patterns shape the epigenetic landscape nucleosomes directs DNA methylation in plants. Cell 151, 167–180 (2012).
in Arabidopsis. Nat. Struct. Mol. Biol. 21, 64–72 (2014). 49. Liu, C. et al. Genome-wide analysis of chromatin packing in Arabidopsis
19. Lindroth, A. M. et al. Requirement of CHROMOMETHYLASE3 for thaliana at single-gene resolution. Genome Res. 26, 1057–1068 (2016).
maintenance of CpXpG methylation. Science 292, 2077–2080 (2001). 50. Feng, S. et al. Genome-wide Hi-C analyses in wild-type and mutants reveal
20. Jackson, J. P., Lindroth, A. M., Cao, X. & Jacobsen, S. E. Control of CpNpG high-resolution chromatin interactions in Arabidopsis. Mol. Cell 55, 694–707
DNA methylation by the KRYPTONITE histone H3 methyltransferase. (2014).
Nature 416, 556–560 (2002). 51. Ibarra, C. A. et al. Active DNA demethylation in plant companion cells
21. Gouil, Q. & Baulcombe, D. C. DNA methylation signatures of the plant reinforces transposon methylation in gametes. Science 337, 1360–1364 (2012).
chromomethyltransferases. PLoS Genet. 12, e1006526 (2016). 52. Calarco, J. P. et al. Reprogramming of DNA methylation in pollen guides
22. Cao, X. & Jacobsen, S. E. Role of the arabidopsis DRM methyltransferases in de epigenetic inheritance via small RNA. Cell 151, 194–205 (2012).
novo DNA methylation and gene silencing. Curr. Biol. 12, 1138–1144 (2002). 53. Walker, J. et al. Sexual-lineage-specific DNA methylation regulates meiosis in
23. Papareddy, R. K. et al. Chromatin regulates expression of small RNAs to help Arabidopsis. Nat. Genet. 50, 130–137 (2018).
maintain transposon methylome homeostasis in Arabidopsis. Genome Biol. 54. Slotkin, R. K. et al. Epigenetic reprogramming and small RNA silencing of
21, 251 (2020). transposable elements in pollen. Cell 136, 461–472 (2009).
24. Singh, J., Freeling, M. & Lisch, D. A position effect on the heritability of 55. Woo, H. R., Pontes, O., Pikaard, C. S. & Richards, E. J. VIM1, a
epigenetic silencing. PLoS Genet. 4, e1000216 (2008). methylcytosine-binding protein required for centromeric
25. Johannes, F. et al. Assessing the impact of transgenerational epigenetic heterochromatinization. Genes Dev. 21, 267–277 (2007).
variation on complex traits. PLoS Genet. 5, e1000530 (2009). 56. Du, J. et al. Mechanism of DNA methylation-directed histone methylation by
26. Reinders, J. et al. Compromised stability of DNA methylation and transposon KRYPTONITE. Mol. Cell 55, 495–504 (2014).
immobilization in mosaic Arabidopsis epigenomes. Genes Dev. 23, 939–950 57. Wendte, J. M. & Schmitz, R. J. Specifications of targeting heterochromatin
(2009). modifications in plants. Mol Plant, https://doi.org/10.1016/j.molp.2017.10.002
27. Vongs, A., Kakutani, T., Martienssen, R. A. & Richards, E. J. Arabidopsis (2017).
thaliana DNA methylation mutants. Science 260, 1926–1928 (1993). 58. Li, X. et al. Mechanistic insights into plant SUVH family H3K9
28. Finnegan, E. J., Peacock, W. J. & Dennis, E. S. Reduced DNA methylation in methyltransferases and their binding to context-biased non-CG DNA
Arabidopsis thaliana results in abnormal plant development. Proc. Natl Acad. methylation. Proc. Natl Acad. Sci. USA 115, E8793–E8802 (2018).
Sci. USA 93, 8449–8454 (1996). 59. Law, J. A. et al. Polymerase IV occupancy at RNA-directed DNA methylation
29. Rigal, M. et al. Epigenome confrontation triggers immediate reprogramming sites requires SHH1. Nature 498, 385–389 (2013).
of DNA methylation and transposon silencing in Arabidopsis thaliana F1 60. Gallego-Bartolome, J. et al. Co-targeting RNA polymerases IV and V
epihybrids. Proc. Natl Acad. Sci. USA 113, E2083–E2092 (2016). promotes efficient De Novo DNA methylation in Arabidopsis. Cell 176,
30. Ito, T. et al. Genome-wide negative feedback drives transgenerational DNA 1068–1082 e1019 (2019).
methylation dynamics in Arabidopsis. PLoS Genet. 11, e1005154 (2015). 61. To, T. K. et al. RNA interference-independent reprogramming of DNA
31. Williams, B. P., Pignatta, D., Henikoff, S. & Gehring, M. Methylation-sensitive methylation in Arabidopsis. Nat. Plants 6, 1455–1467 (2020).
expression of a DNA demethylase gene serves as an epigenetic rheostat. PLoS 62. Teixeira, F. K. & Colot, V. Gene body DNA methylation in plants: a means to
Genet. 11, e1005142 (2015). an end or an end to a means? EMBO J. 28, 997–998 (2009).
32. Saze, H., Shiraishi, A., Miura, A. & Kakutani, T. Control of genic DNA 63. Wendte, J. M. & Schmitz, R. J. Specifications of targeting heterochromatin
methylation by a jmjC domain-containing protein in Arabidopsis thaliana. modifications in plants. Mol. Plant 11, 381–387 (2018).
Science 319, 462–465 (2008). 64. Hsieh, T. F. et al. Genome-wide demethylation of Arabidopsis endosperm.
33. Rigal, M., Kevei, Z., Pelissier, T. & Mathieu, O. DNA methylation in an intron Science 324, 1451–1454 (2009).
of the IBM1 histone demethylase gene stabilizes chromatin modification 65. Kawakatsu, T., Nery, J. R., Castanon, R. & Ecker, J. R. Dynamic DNA
patterns. EMBO J. 31, 2981–2993 (2012). methylation reconfiguration during seed development and germination.
34. Mathieu, O., Reinders, J., Caikovski, M., Smathajitt, C. & Paszkowski, J. Genome Biol. 18, 171 (2017).
Transgenerational stability of the Arabidopsis epigenome is coordinated by 66. Narsai, R. et al. Extensive transcriptomic and epigenomic remodelling occurs
CG methylation. Cell 130, 851–862 (2007). during Arabidopsis thaliana germination. Genome Biol. 18, 172 (2017).
35. Zhang, Y., Wendte, J. M., Ji, L. & Schmitz, R. J. Natural variation in DNA 67. Bouyer, D. et al. DNA methylation dynamics during early plant life. Genome
methylation homeostasis and the emergence of epialleles. Proc. Natl Acad. Sci. Biol. 18, 179 (2017).
USA 117, 4874–4884 (2020). 68. Woo, H. R., Dittmer, T. A. & Richards, E. J. Three SRA-domain methylcytosine-
36. Bewick, A. J. et al. On the origin and evolutionary consequences of gene body binding proteins cooperate to maintain global CpG methylation and epigenetic
DNA methylation. Proc. Natl Acad. Sci. USA 113, 9111–9116 (2016). silencing in Arabidopsis. PLoS Genet. 4, e1000156 (2008).
37. Kiefer, C. et al. Interspecies association mapping links reduced CG to TG 69. Zhang, X. et al. Whole-genome analysis of histone H3 lysine 27 trimethylation
substitution rates to the loss of gene-body methylation. Nat. Plants 5, 846–855 in Arabidopsis. PLoS Biol. 5, e129 (2007).
(2019). 70. Luo, C. et al. Integrative analysis of chromatin states in Arabidopsis identified
38. Wendte, J. M. et al. Epimutations are associated with potential regulatory mechanisms for natural antisense transcript production.
CHROMOMETHYLASE 3-induced de novo DNA methylation. Elife 8, Plant J, https://doi.org/10.1111/tpj.12017 (2012).
https://doi.org/10.7554/eLife.47891 (2019). 71. Yelagandula, R. et al. The histone variant H2A.W defines heterochromatin
39. Williams, B. P. & Gehring, M. Principles of epigenetic homeostasis shared and promotes chromatin condensation in Arabidopsis. Cell 158, 98–109
between flowering plants and mammals. Trends Genet 36, 751–763 (2020). (2014).

12 NATURE COMMUNICATIONS | (2021)12:6958 | https://doi.org/10.1038/s41467-021-27320-6 | www.nature.com/naturecommunications

NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-27320-6 ARTICLE

72. Kujirai, T. & Kurumizaka, H. Transcription through the nucleosome. Curr. population including Claus Schwechheimer in TUM for letting us perform DNA isola-
Opin. Struct. Biol. 61, 42–49 (2020). tions in his lab. This study was supported by the National Science Foundation (MCB-
73. Alonso, J. M. et al. Genome-wide insertional mutagenesis of Arabidopsis 1856143) and the National Institutes of Health (R01GM134682) to R.J.S. F.J. and R.J.S.
thaliana. Science 301, 653–657 (2003). acknowledge support from the Technical University of Munich-Institute for Advanced
74. Urich, M. A., Nery, J. R., Lister, R., Schmitz, R. J. & Ecker, J. R. MethylC-seq Study funded by the German Excellent Initiative and the European Seventh Framework
library preparation for base-resolution whole-genome bisulfite sequencing. Program under grant agreement no. 291763. F.J., R.S.P., and I.K. were supported by the
Nat. Protoc. 10, 475–483 (2015). SFB Sonderforschungsbereich924 of the Deutsche Forschungsgemeinschaft (DFG).
75. Shahryary, Y. et al. AlphaBeta: computational inference of epimutation rates
and spectra from high-throughput DNA methylation data in plants. Genome
Biol. 21, 260 (2020).
Author contributions
Y.Z., R.J.S. and F.J. conceived the study. I.K. performed experiments. Y.Z., H.J., R.X.,
76. Taudt, A. et al. METHimpute: imputation-guided construction of complete
R.S.P. performed computational analyses. R.J.S. wrote the manuscript with assistance
methylomes from WGBS data. BMC Genomics 19, 444 (2018).
from Y.Z. and F.J.
77. Schultz, M. D. et al. Human body epigenome maps reveal noncanonical DNA
methylation variation. Nature 523, 212–216 (2015).
78. Martin, M. Cutadapt removes adapter sequences from high-throughput Competing interests
sequencing reads. EMBnet. J. 17, 10–12 (2011). The authors declare no competing interests.
79. Berardini, T. Z. et al. The Arabidopsis information resource: Making and
mining the “gold standard” annotated reference plant genome. Genesis 53,
474–485 (2015). Additional information
80. Langmead, B. & Salzberg, S. L. Fast gapped-read alignment with Bowtie 2. Nat. Supplementary information The online version contains supplementary material
Methods 9, 357–359 (2012). available at https://doi.org/10.1038/s41467-021-27320-6.
81. Inagaki, S. et al. Gene-body chromatin modification dynamics mediate
epigenome differentiation in Arabidopsis. EMBO J. 36, 970–980 (2017). Correspondence and requests for materials should be addressed to Frank Johannes or
82. Takuno, S. & Gaut, B. S. Body-methylated genes in Arabidopsis thaliana are Robert J. Schmitz.
functionally important and evolve slowly. Mol. Biol. Evol. 29, 219–227 (2012).
Peer review information Nature Communications thanks Leandro Quadrana and the
83. Broman, K. W., Wu, H., Sen, S. & Churchill, G. A. R/qtl: QTL mapping in
other, anonymous, reviewer(s) for their contribution to the peer review of this work.
experimental crosses. Bioinformatics 19, 889–890 (2003).
84. Servant, N. et al. HiC-Pro: an optimized and flexible pipeline for Hi-C data
Reprints and permission information is available at http://www.nature.com/reprints
processing. Genome Biol. 16, 259 (2015).
85. Ay, F., Bailey, T. L. & Noble, W. S. Statistical confidence estimation for Hi-C
data reveals regulatory chromatin contacts. Genome Res. 24, 999–1011 (2014). Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in
86. Li, H. et al. The sequence alignment/map format and SAMtools. published maps and institutional affiliations.
Bioinformatics 25, 2078–2079 (2009).
87. Tarasov, A., Vilella, A. J., Cuppen, E., Nijman, I. J. & Prins, P. Sambamba: fast
processing of NGS alignment formats. Bioinformatics 31, 2032–2034 (2015). Open Access This article is licensed under a Creative Commons
88. Zhang, Y. et al. Model-based analysis of ChIP-Seq (MACS). Genome Biol. 9, Attribution 4.0 International License, which permits use, sharing,
R137 (2008). adaptation, distribution and reproduction in any medium or format, as long as you give
89. Zhang, Y. et al. Heterochromatin is a quantitative trait associated with appropriate credit to the original author(s) and the source, provide a link to the Creative
spontaneous epiallele formation. Zenodo https://doi.org/10.5281/ Commons license, and indicate if changes were made. The images or other third party
zenodo.5651358 (2021). material in this article are included in the article’s Creative Commons license, unless
indicated otherwise in a credit line to the material. If material is not included in the
article’s Creative Commons license and your intended use is not permitted by statutory
regulation or exceeds the permitted use, you will need to obtain permission directly from
Acknowledgements the copyright holder. To view a copy of this license, visit http://creativecommons.org/
The authors would like to acknowledge David Hall for contributions to the discussion of
licenses/by/4.0/.
evolutionary origins of spontaneous epialleles and William Jordan for help creating
ibm1;met1 epiRIL-12. We also thank Thomas Altmann, Rhonda C. Meyer and the entire
group of Heterosis at the IPK for helping with the cultivation of the ddm1-epiRIL © The Author(s) 2021

NATURE COMMUNICATIONS | (2021)12:6958 | https://doi.org/10.1038/s41467-021-27320-6 | www.nature.com/naturecommunications 13