Generation of met1 mutants
We used CRISPR-Cas9 mutagenesis [28] to target the MET1 gene (AT5G49160) in 18 early-flowering accessions of A. thaliana, creating frameshift mutations in exon 7 (“Methods,” Additional file 1: Table S1). For each accession, transgene-free lines were genotyped at the MET1 locus and propagated. In the next generation, we obtained homozygous met1 mutants from 17 accessions, with two mutant lines of independent origin for 14 accessions and one mutant line for Cvi-0, Ler-1, and Col-0. For Bl-1, we did not recover a sufficient number of homozygous progeny for in-depth analysis. Because Bl-1 heterozygotes already had morphological defects, we included these in our analyses, along with heterozygous mutants in Col-0, Ler-1, and Bu-0. We also included second-generation met1 homozygotes from Tsu-0 and Tscha-1 (descended from siblings of first-generation homozygotes) to glean first insights into progressive changes at later generations of homozygosity. For one Bs-1 line with bi-allelic mutations in MET1, we only analyzed second-generation progeny that was homozygous for one of the two alleles.
Previous work showed that epigenetic states across the genome can diverge in different lineages of met1 mutants over several generations [5, 11]. To ensure that we could directly link chromatin state and gene expression, we performed paired BS-seq, ATAC-seq, and RNA-seq on leaf tissue of the same plant rosettes, collected as three biological replicates for both wild-type and met1 mutant lines (“Methods”). All together, we obtained 73 BS-seq, 158 ATAC-seq, and 158 RNA-seq libraries that passed quality control.
MET1 can both buffer and increase transcriptomic variation across accessions
Since natural accessions of A. thaliana are known to express diverse transcriptomes in their wild-type state [16], we hypothesized that MET1-induced CG methylation may contribute to generating this diversity and that the transcriptomes would become more similar to each other in met1 mutants. We therefore first grouped all wildtypes (54 samples) and all met1 mutants (104 samples) separately, and analyzed 19,473 genes with sufficient read counts in both groups. Contrary to our hypothesis, UMAP projections of read counts at these genes revealed greater differences between accessions in met1 mutants compared to their respective wild-type samples (Fig. 1a). This suggested an alternate hypothesis where MET1 functions in buffering transcriptomic diversity across accessions, with its absence therefore unmasking larger regulatory differences. Wild-type samples had 1210 differentially expressed genes (DEGs; defined as genes with |log2(fold change)|≥1, FDR adjusted p-value≤0.01) (Additional file 2: Table S2, Additional file 3: Dataset S1) between accessions, fewer than those identified between accessions in met1 mutant samples (1868 DEGs) (Additional file 2: Table S2, Additional file 4: Dataset S2). There were only 406 DEGs that were shared by wildtypes and met1, i.e., which were differentially expressed independently of MET1 activity. These are likely to include genes that are associated with structural variation or major differences in cis-regulatory sequences. There were 804 DEGs that were unique to wild-type samples, i.e., where differences in expression between accessions were greatly reduced upon loss of MET1. Finally, the largest group in this comparison comprised 1462 DEGs that were unique to met1 samples, i.e., which were expressed at similar levels across wildtypes, but became differentially expressed in met1 mutants. These inferences were also apparent from heatmaps (Additional file 5: Fig. S1) and density distributions (Additional file 5: Fig. S2) of Pearson correlation coefficients between accessions. Together, these observations indicate that there is evidence for both of our alternative hypotheses—MET1 reduces transcriptomic diversity across accessions for many genes, but for a smaller group of genes it adds another layer of expression complexity to diversify transcriptomes.
TE- and Non-TE genes are differentially expressed in met1 mutants
A major role of MET1 is the transcriptional silencing of TEs by establishing and maintaining CG methylation [29, 30]. We next asked whether there were any accession-specific differences in the requirement of MET1 in regulating the expression of both TE genes (genes associated with TEs, see “Methods”) and other protein-coding genes. We first analyzed RNA-seq read counts of 21,657 genes in all samples (wildtypes and met1 mutants) and generated a UMAP visualization. Two distinct clusters of samples could be observed, one encompassing wild-type and the other met1 plants (Fig. 1b). In agreement with our initial results, wild-type samples were less spread out than mutant samples, indicating greater gene expression heterogeneity in met1 mutants than in wild-type plants. Comparison of first- and second-generation homozygous mutants in the Tscha-1 and Tsu-0 accessions did not indicate major changes upon propagation of mutant lines.
To obtain first insights into MET1 function across all accessions, we examined protein-coding genes whose expression levels changed significantly in a contrast of all met1 mutants against all wild-type plants, thereby identifying 3479 DEGs. Of these, 1466 genes (42% of all DEGs) were associated with TEs, which we called TE-DEGs (Additional file 2: Table S2, Additional file 6: Dataset S3), and these corresponded to 87% of the 1678 TE-associated genes with sufficient information in our dataset (“Methods”). Almost all of them were upregulated in met1 mutants compared to wild-type plants, and often greatly so (Fig. 1c). This was consistent with previous findings from early- and late-generation homozygous met1 mutants [31,32,33,34,35], including the observation that Class II DNA TEs of the En/Spm superfamily were most often among activated TE-DEGs [36, 37]. Consistent with TE density being highest near the centromeres [38], TE-DEGs were enriched in pericentromeric regions (Fig. 1d). Many of the TE-DEGs were strongly overexpressed in met1 mutants, hundreds of times or more (Fig. 1c).
The picture was different for the remaining 2013 DEGs that were not associated with TEs (Non-TE-DEGs, Additional file 2: Table S2, Additional file 6: Dataset S3), and which corresponded to 10% of the 19,979 Non-TE genes with sufficient information in our dataset. Although the majority was also upregulated in met1 mutants, more than a third, 728, was downregulated (Fig. 1c). Non-TE-DEGs were also more uniformly distributed along the chromosome arms, and overall expression changes for both up- and downregulated genes were much more moderate (Fig. 1d). Gene Ontology (GO) enrichment of Non-TE-DEGs revealed several terms related to abiotic and biotic stresses and stimulus response (Additional file 5: Fig. S3). We conclude that across all accessions, Non-TE-associated genes vary much more in their sensitivity to loss of MET1 than TE-associated genes.
Mis-regulated genes vary among met1 mutants of different accessions
A closer examination of DEGs called by contrasting all met1 mutants against all wildtypes showed that DEGs were not uniformly induced or repressed across accessions (a random subset of such DEGs is shown in Additional file 5: Fig. S4). Hence we individually examined each of the 18 accessions for genes that were differentially expressed between met1 mutants and the corresponding wild-type parents. We found that the number of DEGs in each accession varied substantially (Fig. 1e, Additional file 2: Table S2, Additional file 6: Dataset S3), but that this was much more true for Non-TE-DEGs than TE-DEGs. In addition, while downregulated TE-DEGs were rare in all accessions, the ratio of up- to downregulated Non-TE-DEGs in met1 mutants was much more variable (Fig. 1f, Additional file 5: Fig. S5). The number of Non-TE-DEGs ranged from 278 (26% of all DEGs) in Com-1 to 3409 (77% of all DEGs) in Est. Notably, even though we could only analyze heterozygous Bl-1 met1 mutants, these had more Non-TE-DEGs than homozygous met1 mutants in several other accessions, suggesting that even the removal of only one of the two functional MET1 copies was sufficient to alter expression levels of a large number of genes, at least in this accession.
In total, there were 10,151 Non-TE-DEGs and 1524 TE-DEGs that were differentially expressed in at least one of the accessions or in the all-met1-against-all-wild-type contrast. Expression levels of both DEG sets were significantly more variable across accessions in met1 mutants compared to wildtypes (Fig. 1g, Additional file 5: Fig. S6), once again demonstrating that expression diversity across accessions was overall higher in the absence of MET1 activity.
We used all 10,151 Non-TE-DEGs to build a weighted gene co-expression network (“Methods”), finding nine modules. Genes from module “D” were the most consistent in their expression levels across all accessions, being on average always upregulated in met1 mutants compared to the respective wild-type plants (Additional file 5: Fig. S7, Additional file 7: Dataset S4). GO enrichment analyses revealed that many “D module” genes were associated with nucleic acid metabolic processes, DNA repair and transcription, although only weakly significantly so. This suggests that MET1 likely affects core metabolic machinery genes, which may further impact a different subset of downstream genes in each accession.
We next generated a frequency spectrum to examine the overlap of DEGs between accessions. When we focused on the two extreme accessions, Est and Com-1, a great majority of TE-DEGs were shared across most contrasts, while the picture for Non-TE-DEGs was very different. While almost 30% of the 3409 Non-TE-DEGs in Est were not detected in any other accession, 98% of the 278 Com-1 Non-TE-DEGs were found in at least one other accession (Fig. 2a,b, Additional file 5: Fig. S8). The 983 unique Non-TE-DEGs in Est were enriched for functions with a common theme of RNA and DNA metabolism (Fig. 2c), compatible with a scenario in which mis-regulation of one or a few specific master regulators had led to expression changes in numerous Non-TE-DEGs in Est.
When overlaying DEGs from the 18 accession-specific met1-against-wild-type comparisons and the all-met1-against-all-wild-type comparison, we found 291 universal DEGs (Additional file 8: Dataset S5), albeit with the extent of expression change in met1 mutants differing considerably across accessions (Fig. 1g, Additional file 5: Fig. S6). Only 15 of the universal DEGs were Non-TE-DEGs (Fig. 2d), and eight of these 15 genes are strongly expressed in siliques in Col-0 wild-type plants [39]. We had measured gene expression in rosette leaves, and many of these genes were expressed at low levels in our wild-type samples, but strongly upregulated in met1 mutants. Included among the 15 genes were three maternally imprinted genes, FWA, SDC, and AT1G59930 [40, 41], along with AT5G35120, which we suspect may be maternally imprinted as well, given its sequence similarity with AT1G59930 and selective expression in the endosperm. Ectopic FWA expression is known to delay flowering in both Col-0 and Ler-1 accessions [42, 43], while ectopic activation of SDC has been shown to lead to dwarfism and leaf curling in Col-0 [41, 44]. This is consistent with the whole-plant phenotypes described in detail below. Four consistently observed DEGs encoded an ULP-1 protease family domain. ULP-1 sequences have been found to be associated with Mutator-like TEs in Cucumis melo (“CUMULEs”), although there is no evidence that similar CUMULEs in rice or A. thaliana assist in TE mobilization [45]. Because these genes are not annotated as “TE genes” in TAIR10, they are included in our list of Non-TE-DEGs.
Finally, we assessed to what extent DEGs from one accession changed in the other accessions. Genes identified as up- or downregulated Non-TE-DEGs after loss of MET1 in one accession changed on average almost always in the same direction in each of the other accessions, but did so less strongly (Fig. 2e, f). TE-DEGs on the other hand were similarly upregulated in all accessions, but downregulated TE-DEGs were often variably expressed in other accessions (Additional file 5: Fig. S9). This finding is consistent with the idea that MET1 activity often homogenizes gene expression levels in different genetic backgrounds.
met1 mutants exhibit reduced DNA methylation and increased overall chromatin accessibility
To investigate how much of the variation in DEGs across accessions arose from differences in DNA methylation and chromatin architecture, we characterized the methylomes and chromatin accessibility in our collection of met1 mutants and corresponding wild-type lines. As expected from the knockout of MET1, cytosine methylation in the CG sequence context was drastically and consistently reduced in homozygous met1 mutants, dropping to an average genome-wide level of 0.2%, representing a 98.5% decrease from mean wild-type levels (Additional file 5: Fig. S10). Cytosine methylation in non-CG contexts was moderately affected in first-generation met1 mutants, being on average 6.8% higher in the CHG context and 21% lower in the CHH context. Second-generation met1 mutants of Tsu-0 and Tscha-1 had greatly increased CHG methylation, both relative to first-generation met1 mutants and wild-type parental lines (Additional file 5: Fig. S10). This increase in methylation during successive rounds of inbreeding has previously been described for met1-3 mutants in the Col-0 accession [11].
Because overall methylation patterns were greatly altered in met1 mutants, we wanted to closely examine genomic regions with the most significant changes in methylation. We contrasted 73 methylomes including both met1 mutants and the wild-type parents from all 18 accessions (“Methods”) to identify differentially methylated regions (DMRs) in the CG, CHG, and CHH contexts (Fig. 3a, b, Additional file 5: Fig. S11, Additional file 9: Dataset S6, Additional file 10: Dataset S7, Additional file 11: Dataset S8). The 2388 CG-DMRs overlapped with the vast majority of the 350 CHG-DMRs and 1023 CHH-DMRs. Approximately half of all CG-DMRs overlapped with TE sequences, a quarter overlapped TE genes, while 42% overlapped Non-TE genes (Additional file 5: Fig. S10, Additional file 2: Table S2).
While most CG-DMRs in met1 mutants retained minimal or no CG methylation, the extent of reduction in methylation differed across accessions, consistent with accession-specific methylation patterns in the presence of MET1 [16] (Fig. 3a, b, Additional file 5: Fig. S12). To assess the functional consequences of this variation, if any, we asked how loss of cytosine methylation impacted genome-wide chromatin architecture in each accession. To this end, we identified accessible chromatin regions (ACRs) by ATAC-seq in all met1 mutants and the corresponding wildtypes. With this analysis, which allowed us to define differential ACRs (dACRs) with accessibility changes in at least two accessions (“Methods,” Additional file 12: Dataset S9), we could focus on 9505 highly variable dACRs (HV-dACRs) with particularly stark variation in accessibility across accessions and genotypes (“Methods,” Additional file 13: Dataset S10). We visualized variation in accessibility levels at these HV-dACRs using UMAP (“Methods”), which similarly to the RNA-seq data (Fig. 1b) revealed two distinct clusters of wild-type plants and met1 mutants (Fig. 3c), with wild-type plants from different accessions being more similar to each other than met1 mutants from different accessions. Approximately one third of all HV-dACRs overlapped in position with Non-TE genes, 24% with TE genes, and 62% with TE sequences (Additional file 5: Fig. S13, Additional file 2: Table S2).
Genome wide, the chromatin of met1 mutants was more accessible than chromatin of wild-type plants, and this difference was particularly pronounced for two subgroups, k-groups 1 and 3, of HV-dACRs identified by k-means clustering (Fig. 3d). Across all accessions, approximately one third (31%) of CG-DMRs overlapped with HV-dACRs, but at the same time, most HV-dACRs (88%) did not overlap with DMRs (Additional file 2: Table S2), indicating that the vast majority of HV-dACRs appeared due to trans effects of methylation changes in the genome.
To see whether epigenetic profiles at TEs, which included both TE genes and other TE sequences, differed from Non-TE genes, we averaged methylation levels using genome-wide methylated cytosines in all contexts (“Methods”) and chromatin accessibility levels using 34,993 consensus ACRs (“Methods”) for all met1 mutants and wild-type plants across 31,189 TEs and 29,699 Non-TE genes including 1 kb flanking sequences (Additional file 5: Fig. S14). Genome-wide chromatin accessibility patterns in met1 mutants mirrored cytosine methylation levels, with increases in accessibility accompanied by decreases in methylation for most accessions. This pattern was much more pronounced in magnitude over TEs (Additional file 5: Fig. S14a) than in Non-TE genes (Additional file 5: Fig. S14b). These observations once again demonstrate that TEs are highly sensitive to the absence of MET1, in agreement with TE genes being strongly upregulated in met1 mutants.
Finally, we asked how much of the expression changes at DEGs were explained by these epigenetic alterations. We focused on Est and Com-1, the two accessions with the highest and lowest number of DEGs. While TE-DEGs behaved similarly in the two accessions, Non-TE-DEGs showed contrasting patterns, with accessibility differences between met1 mutants and wildtypes being on average much smaller in Est than in Com-1 (Additional file 5: Fig. S14c), confirming a complex relationship between methylation, chromatin accessibility and gene expression changes at Non-TE-DEGs, which in turn leads to very different numbers of Non-TE-DEGs in different accessions.
Non-TE-DEGs show varying methylation and chromatin accessibility profiles across accessions
Observing that differential gene expression can be accompanied by epigenetic changes, we quantitatively assessed mutual relationships of cytosine methylation and chromatin accessibility with gene expression, having analyzed all of the three factors from the same individuals across our collection of met1 mutants and wildtypes.
To investigate whether the presence or absence of methylation within or in proximity to a gene could influence its expression level, we examined how DEGs were regulated in the presence of a DMR. We first defined a consensus set of 7132 DEGs from all 18 accessions (“Methods,” Additional file 5: Fig. S15, Additional file 14: Dataset S11, Additional file 15: Dataset S12), containing TE-DEGs (1401) and Non-TE-DEGs (5731) and intersected their genomic coordinates with 1569 CG-DMRs with sufficient methylation calls across all samples that were identified from all mutant and wild-type samples (“Methods,” Additional file 16: Dataset S13).
While only a small fraction of consensus DEGs were located close to a CG-DMR (21% of TE-DEGs and 7% of Non-TE-DEGs) (Additional file 2: Table S2), the converse was also true: only a minority of CG-DMRs was found next to DEGs (21% next to TE-DEGs and 27% to Non-TE-DEGs). Local differences in CG methylation are thus neither necessary nor sufficient for differences in gene expression between wild-type and met1 mutant plants. At the loci where both expression and CG methylation were altered, we examined expression counts and methylation levels in homozygous met1 mutants and in corresponding wild-type plants from 17 accessions. As a control, we randomly sampled Non-DEGs (genes that were not differentially expressed; “Methods”) that positionally overlapped with DMRs.
Most TE-DEGs with CG-DMRs in either their extended gene bodies (“Methods,” Additional file 5: Fig. S16, Additional file 5: Fig. S17c) or 1.5 kb up- or downstream sequences (“cis” regions) (“Methods,” Additional file 5: Fig. S16, Additional file 5: Fig. S18c) had lost CG methylation in met1 mutants. These TE-DEGs were upregulated in all met1 mutants, albeit to a different extent in each accession. Non-TE-DEGs, which were already observed to be very variable in their differential expression levels, exhibited a wide gradient of methylation differences between met1 mutants and the respective wild-type parents (ranging from −100 to +10%). This was observed both when DMRs were located in extended gene bodies (Additional file 5: Fig. S17a) and in 1.5 kb up- or downstream sequences (Additional file 5: Fig. S18a). In both cases, one group of Non-TE-DEGs was strongly upregulated and had highly reduced methylation in met1 mutants, suggesting that methylation could have a strong effect on gene regulation at these genes. There was also another group of more moderately up- or downregulated Non-TE-DEGs with negligible methylation changes in met1 mutants. On examining the accession-of-origin of each of these Non-TE-DEGs, we found that the same gene in different accessions could be found in either of the groups identified above, indicating considerable epigenetic plasticity across accessions.
Consequently, we asked whether parental methylation and expression state can predict epigenetic and transcriptional response in met1 mutants. DEGs classified in the top quintiles of wild-type CG methylation levels (“Methods”) were more likely to increase in expression than DEGs from the bottom quintiles, which were similarly likely to be up- or downregulated in met1 mutants (Fig. 4c, g, Additional file 5: Fig. S19c, Fig. S19g). In terms of expression changes, we found that DEGs from the lowest quintile of expression counts in the wild-type parents were the ones that were the most upregulated in met1 mutants (Fig. 4a, e, Additional file 5: Fig. S19a, Fig. S19e). These observations are consistent with high levels of CG methylation in the wild-type state serving to silence genes.
In A. thaliana, many constitutively active genes are marked by gene body CG methylation (gbM) [46] and the methylation levels of these genes are known to vary in tandem with their differential expression across A. thaliana accessions [22, 23]. Among Non-TE-DEGs with gene body CG-DMRs, 91 out of 196 genes (46%) were methylated in the wild-type state (“gbM like” genes; “Methods”) in at least one accession. Most of these “gbM-like” genes were downregulated in met1 mutants (Additional file 5: Fig. S20a). When we examined the larger set of all Non-TE genes (19,979), we found three genes that were consistently gbM-like in wildtypes of 17 accessions, but with highly reduced methylation in homozygous met1 mutants. These were AT5G20490 (XIK, encoding a myosin family protein involved in root hair growth and trichome development), AT5G37130 (encoding a prenylyltransferase superfamily protein), and AT5G44800 (PICKLE, encoding a protein that affects histone methylation levels and impacts floral meristem identity) (Additional file 5: Fig. S20b). Only AT5G37130 was a consensus Non-TE-DEG (being significantly downregulated in met1 mutants in Aa-0, Com-1, Pi-0 and Cvi-0) (Additional file 5: Fig. S20c). Incidentally, genes that were highly methylated and weakly expressed in wildtypes (exhibiting TE methylation characteristics in the CG context or “CG teM-like” genes; “Methods”) were always upregulated in met1 mutants (Additional file 5: Fig. S21), similarly to TE-DEGs. The two groups of genes, with gene body or TE methylation in the CG context, did not overlap.
We next carried out a similar quantitative analysis to examine the relationship between differential expression changes and differential chromatin accessibility at DEGs (“Methods,” Additional file 5: Fig. S22). Approximately 37% of all consensus DEGs (“Methods”) were associated with HV-dACRs within 1.5 kb upstream and downstream of their gene body (Additional file 2: Table S2). While chromatin accessibility increased almost proportionally with expression in met1 mutants for TE-DEGs, increased accessibility could be associated with, but did not necessitate an increase in gene expression for Non-TE-DEGs (Additional file 5: Fig. S23, Additional file 5: Fig. S24). Accessible chromatin is well known to favor gene transcription by facilitating transcription-factor binding [47, 48] although there is also evidence that inaccessible regions can occur in some long and highly transcribed genes [49]. Together, these observations point to complex interactions between gene expression and chromatin accessibility levels in cis, especially at Non-TE-DEGs.
MET1 can have indirect effects on the expression of Non-TE genes
Since variation in the response of gene expression, methylation and chromatin accessibility to loss of MET1 was most apparent for Non-TE-DEGs, we focused on these genes to explore the nature of their epigenetic plasticity. Among 5731 consensus Non-TE-DEGs that varied in their expression response across accessions (Fig. 5a), 21% had differentially accessible chromatin regions (HV-dACRs) in their vicinity, but lacked a cis CG-DMR (“cis” here including the gene body) (Additional file 2: Table S2). Conversely, a minority of Non-TE-DEGs, 255 (5%) had cis CG-DMRs but no nearby HV-dACRs, and a majority of Non-TE-DEGs, 4105 (72%) had neither a nearby HV-dACR nor CG-DMR (Additional file 2: Table S2). Closer inspection of several genes from these categories, AT1G60190, PR1, ROS1 (Fig. 5a), showed that the relationship between altered gene expression, CG methylation, and accessibility in met1 mutants varied substantially, both across genes and accessions.
To focus on close-range interactions between changes in methylation, chromatin accessibility, and gene expression at the same locus, we examined Non-TE genes associated with both cis HV-dACRs and CG-DMRs—comprising 164 DEGs, which we compared with 107 randomly sampled Non-DEGs (“Methods”), by overlaying methylation, accessibility, and gene expression data for each gene across 17 accessions. We observed that 3073 genes (in 17 accessions) with strongly reduced CG methylation (>75% methylation reduction) had a tendency to have increased chromatin accessibility. Of these, 15.6% were strongly upregulated in met1 mutants, while only 1% were strongly downregulated (Fig. 5b). Incidentally, strongly upregulated genes in this group also stood out from the rest because they exhibited moderately increased accessibility in met1 mutants, consistent with methylation changes being directly responsible for induction of expression (Fig. 5b). For both DEGs and Non-DEGs (Fig. 5b, c), there were many other different combinations of methylation and accessibility states, and these did not always cluster by degree of expression change. Together, these observations suggested the presence of multiple epigenetic states, both for different genes in the same accession and for the same gene in different accessions.
We next examined Non-TE genes (271 genes including 164 DEGs and 107 Non-DEGs) as a group for accession-specific epigenetic patterns, first comparing the methylation levels of these genes in wildtypes and mutants. Cvi-0 had the highest fraction of genes with minimal reduction in methylation in met1 mutants (Additional file 5: Fig. S25a). This was explained by Cvi-0 wild-type having already many genes with low methylation levels, limiting the extent of any further reduction in methylation (Additional file 5: Fig. S25a). We asked how genes with low CG methylation in Cvi-0 wild-type, that is, genes with methylation levels that could not be reduced much further by loss of MET1, fared in other accessions. While the average reduction in methylation level was greater in other accessions, as expected, changes in accessibility and expression level after inactivation of MET1 were similar in magnitude when compared to Cvi-0 (Fig. 5d). This observation provides further support for genome-wide hypomethylation indirectly affecting the expression of many genes. Finally, comparing the relationship between MET1-dependent changes in methylation, chromatin accessibility, and gene expression in the reference accession Col-0 with changes in other accessions confirms that Col-0 is not particularly representative of A. thaliana accessions at large (Fig. 5e, f).
met1 mutants express signatures of known epialleles
Many of our met1 mutants had strong methylation and expression changes at well-known loci sensitive to epigenetic regulation, such as FWA [42] (Additional file 5: Fig. S25c,d), SDC [44] (Additional file 5: Fig. S26), the PAI genes [50], IBM1 [51], SNC1 (with alleles similar to the bal variant of SNC1) [52], the ROS1 demethylase, which is known to function as a methylation sensor [11, 53, 54] (Fig. 5a, Additional file 5: Fig. S25b), AG [4] (Additional file 5: Fig. S27), and SUP (with alleles similar to the clark-kent variant of SUP) [55] (Additional file 5: Fig. S28). New epialleles at several of these loci have been reported before in Col-0 and C24 met1 mutants [3, 56]. We observed a variety of methylation patterns at these loci depending on the accession of origin, with considerable differences in chromatin accessibility and gene expression across accessions; several examples at the FWA locus are shown in Additional file 5: Fig. S25c,d. As seen before [11], some epialleles arose only in second-generation met1 mutants (Additional file 5: Fig. S27), suggesting that epialleles continue to accumulate in the absence of MET1 during inbreeding.
met1 mutants vary in phenotypes and segregation distortion
Compared to their wild-type parents, homozygous met1 mutants were dwarfed, flowered late, and had altered rosette leaf architecture—although to different degrees in each accession (Fig. 6a, Additional file 5: Fig. S29-S30). met1 mutants also suffered from silique abnormalities, in some cases affecting fertility (Fig. 6b–d). Additionally, we observed during the preparation of ATAC-seq libraries that met1 mutant nuclei had lower levels of endopolyploidy (Additional file 5: Fig. S31). Homozygous mutants were underrepresented in the progeny of heterozygous parents, consistent with reduced transmission of met1 alleles, as described previously [6, 57]. To estimate the extent of segregation distortion, we grew up to 96 progeny of heterozygous mutants and genotyped all individuals. In all accessions, homozygotes were underrepresented (Fig. 7a, Additional file 17: Table S3), from 2% in Aa-0 to 18% in Uk-1.
Bs-1 Line 2 included two different types of homozygous mutants that were phenotypically distinct, with the frameshift allele, which we used for the genomic analyses described above, causing more severe phenotypic defects and being more strongly underrepresented in a segregating populations (Fig. 7b, Additional file 17: Table S3). Bu-0 Line 2 and Ste-0 Line 2 segregated individuals with different combinations of two mutant alleles along with a wild-type allele, pointing not only a bi-allelic origin of the met1 alleles but also a change in the generative ploidy level. Genotyping by sequencing also revealed heterozygous individuals with a skewed ratio of reads for the wild-type and mutant alleles, especially in Bu-0 Line 1, Pi-0 Line 1, Col-0 Line 2, and Bl-1 Line 1 (Additional file 17: Table S3). Since heterozygotes of these four lines also exhibited phenotypic variation (Fig. 7c, Additional file 5: Fig. S32), we suspected that there was variation in ploidy. Both Bu-0 Line 1 and Line 2 as well as wild-type Bu-0 plants used in this study turned out to be tetraploid, possibly explaining the altered segregation ratios (Fig. 7d, e). However, other lines with skewed heterozygous allele frequencies were derived from diploid parents (Additional file 5: Fig. S32), awaiting a more comprehensive explanation of genotypic and phenotypic variation in these lines.