Research ArticleDEVELOPMENTAL BIOLOGY

Active DNA demethylation regulates tracheary element differentiation in Arabidopsis

See allHide authors and affiliations

Science Advances  26 Jun 2020:
Vol. 6, no. 26, eaaz2963
DOI: 10.1126/sciadv.aaz2963

Abstract

DNA demethylation is important for the erasure of DNA methylation. The role of DNA demethylation in plant development remains poorly understood. Here, we found extensive DNA demethylation in the CHH context around pericentromeric regions and DNA demethylation in the CG, CHG, and CHH contexts at discrete genomic regions during ectopic xylem tracheary element (TE) differentiation. While loss of pericentromeric methylation occurs passively, DNA demethylation at a subset of regions relies on active DNA demethylation initiated by DNA glycosylases ROS1, DML2, and DML3. The ros1 and rdd mutations impair ectopic TE differentiation and xylem development in the young roots of Arabidopsis seedlings. Active DNA demethylation targets and regulates many genes for TE differentiation. The defect of xylem development in rdd is proposed to be caused by dysregulation of multiple genes. Our study identifies a role of active DNA demethylation in vascular development and reveals an epigenetic mechanism for TE differentiation.

INTRODUCTION

DNA methylation [5-methylcytosine (5mC)] is a conserved epigenetic mechanism used by plants and animals to silence transposons and regulate gene expression. In plants, DNA methylation occurs in CG, CHG, and CHH contexts (H represents A, T, or C) (1). A canonical RNA-directed DNA methylation (RdDM) pathway establishes de novo DNA methylation patterns (1). In RdDM, two plant-specific RNA polymerases, Pol IV and Pol V, generate 24-nucleotide small interfering RNAs (siRNAs) and long scaffold RNAs, respectively. Pol V–dependent long scaffold RNAs recruit the ARGONAUTE 4/siRNA complex and then DOMAIN REARRANGED METHYLTRANSFERASE 2, accomplishing RdDM (1). Once established, CG and CHG methylation are maintained by METHYLTRANSFERASE 1 and CHROMOMETHYLASE 3 (CMT3), respectively, whereas CHH methylation is maintained by RdDM or the chromatin remodeler DECREASE IN DNA METHYLATION 1–dependent CMT2 pathway (1, 2).

DNA methylation patterns can be modified by passive or active DNA demethylation. Passive DNA demethylation refers to loss of DNA methylation during DNA replication owing to inactivation or reduced expression of DNA methyltransferases (3). Active DNA demethylation involves enzyme-catalyzed reactions. In Arabidopsis, active DNA demethylation is initiated by ROS1 family of DNA glycosylases REPRESSOR OF SILENCING 1 (ROS1), DEMETER (DME), DEMETER-LIKE 2 (DML2), and DML3 (46) and completed through a base-excision repair pathway (7, 8).

Over the past decade, evidence has accumulated revealing that DNA demethylation regulates plant development (7, 8). In Arabidopsis, CHH methylation is passively lost during seed germination, suggesting a possible role of passive DNA demethylation in seed germination (9, 10). DME-initiated active DNA demethylation, which occurs in the vegetative cell, the companion cell of sperm (11), promotes pollen germination (12). DME-initiated genome-wide active DNA demethylation, which occurs in the central cell, the companion cell of the egg that develops into the endosperm (11, 13), is required for seed development (14). ROS1-initiated active DNA demethylation, which occurs in somatic cells, regulates the initiation of stomatal lineage cells by preventing the spreading of DNA methylation from a transposon in the EPF2 promoter (15). In rice (Oryza sativa) grains, OsROS1-initiated active DNA demethylation prevents the formation of thickened aleurone (16). In Medicago truncatula, MtDME-initiated active DNA demethylation promotes nodule development (17). In tomato (Solanum lycopersicum), active DNA demethylation, initiated by SlDML2, promotes fruit ripening (18). By contrast, in orange (Citrus sinensis), active DNA demethylation, initiated by CsDME, CsDML1, CsDML4, and CsDML3, prevents fruit ripening (19).

In vascular plants, xylem tissue transports water and provides mechanical support (20). Xylem tissue is constituted of tracheary elements (TEs), parenchyma cells, and fibers. Mature TEs are hollow, elongated, and lignified dead cells connected end-by-end. TEs are differentiated from procambial or cambial cells. The differentiation process is divided into four stages: Cell fate determination, xylem-specific gene activation along with cell expansion/elongation, secondary cell wall (SCW) formation, and programmed cell death (PCD) (21). Hormones, transcription factors, peptides, kinases, and other regulatory molecules play important roles in the orchestration of molecular and cellular events for TE differentiation (table S1). Auxin, cytokinin, and their interplay promote the formation of procambial cells in embryos (22, 23) and transdifferentiation of mesophyll cells into TEs in various xylogenic culture system in vitro (2427). Transcription factors including TARGET OF MONOPTEROS 5, LONESOME HIGHWAY (28, 29), BRI1-EMS SUPPRESSOR 1 (BES1), and VASCULAR-RELATED NAC-DOMAIN (VND) family proteins VND1 to VND7 (30) confer the commitment of procambial cells toward xylem cells. VND6 and VND7 can also activate genes involved in SCW formation (e.g., MYB46 and MYB83) and PCD (e.g., XCP2) (25), and they are the master switches of TE differentiation (24). In kinase cascade, TE DIFFERENTIATION INHIBITORY FACTOR (TDIF), a peptide produced by CLAVATA3/EMBRYO SURROUNDING REGION-RELATED (CLE) genes, negatively regulates the differentiation of procambial cells into TEs in vitro and in vivo (31, 32) by binding to its receptor TDIF RECEPTOR (TDR), liberating GLYCOGEN SYNTHASE KINASE 3 proteins (GSK3s) from TDR and activating the GSK3 protein BRASSINOSTEROID-INSENSITIVE 2 (BIN2), which negatively regulates TE differentiation by inhibiting BES1 (33). Despite these findings, whether epigenetic mechanisms, such as changes in DNA methylation, contribute to TE differentiation remains unknown.

In this study, we investigated genome-wide reprogramming of DNA methylation during TE differentiation using the VISUAL (Vascular cell Induction culture System Using Arabidopsis Leaves) xylogenic culture system (26, 34). We found that the Arabidopsis genome passively loses CHH methylation at pericentromeric regions in the process of TE differentiation. The Arabidopsis genome also undergoes CG and CHG demethylation at discrete loci, which is dependent on ROS1, DML2, and DML3 (RDD)–initiated active DNA demethylation. Our genetic analysis revealed that dysfunction of the RDD trio of DNA demethylases impairs TE differentiation in vitro and causes defective xylem development in vivo. Moreover, integrative analysis of the methylomes and transcriptomes of differentiating TEs of the Col-0 wild type (WT) and rdd mutant backgrounds led us to discover a subset of genes regulated by active DNA demethylation for TE differentiation. Our study highlights a novel role of RDD-initiated active DNA demethylation in plant development and provides insights into the regulation of TE differentiation by an epigenetic mechanism.

RESULTS

Arabidopsis mesophyll cells can differentiate into TEs in vitro

VISUAL, an Arabidopsis tissue culture system for in vitro TE differentiation, has previously been established (26, 34). In this system, mesophyll cells in Arabidopsis leaf discs or cotyledons synchronously differentiate into TEs from procambial cells in the presence of auxin, cytokinin, and the GSK3 inhibitor bikinin. To validate the effectiveness of the strategy used in the system, we detected autofluorescence from lignified SCWs, a hallmark feature of differentiated TEs, every 12 hours after application of auxin, cytokinin, and bikinin to Arabidopsis cotyledons (fig. S1, A and B). We found that SCW lignification appeared 36 hours after induction (hai). The percentage of cells with SCW lignification was greater than 85% at 72 hai in WT cotyledons (fig. S1C). We next detected the expression of genes that mark four stages of TE differentiation [reviewed in (35) and listed in table S1]. Expression of procambial cell marker genes TMO6, ARABIDOPSIS THALIANA HOMEOBOX GENE 8 (AtHB-8), PIN-FORMED 7 (PIN7), and TDR was induced at 12 hai and peaked at 24 hai (fig. S1D). Expression of VND5, VND6, VND7, and LOB DOMAIN-CONTAINING PROTEIN 15 (LBD15), marker genes of the transition from procambial cells to TEs, was induced at 36 hai and continued to increase (fig. S1D). Genes involved in SCW formation and PCD were induced at 36 and 48 hai, respectively (fig. S1D). These results suggest that marker gene expression during the formation of procambial cells and the differentiation of procambial cells into TEs can be recapitulated using the VISUAL system.

Pericentromeric CHH methylation gradually decreases during TE differentiation

To investigate the role of DNA methylation in the regulation of gene expression during TE differentiation, we profiled DNA methylomes at sequential stages of TE differentiation using the VISUAL system. At all stages, the centromeric and pericentromeric regions, which are rich in satellite repeats and transposons and form heterochromatic regions, had high levels of DNA methylation, while gene-rich chromosome arms had low levels of DNA methylation in the CG, CHG, and CHH contexts (Fig. 1A). Genome-wide CG and CHG methylation did not substantially change, but centromeric and pericentromeric CHH methylation decreased prominently during TE differentiation (Fig. 1A). Calculation of the average DNA methylation levels of genes revealed that the CHH methylation level, but not the CG and CHG methylation levels, at gene promoters and 3′ untranslated regions was decreased from 0 to 60 hai (Fig. 1B). The CG methylation level of transposons remained unchanged, but their CHG methylation level was decreased slightly, and their CHH methylation level was decreased markedly, especially from 12 to 24 hai (Fig. 1C).

Fig. 1 Genome-wide DNA methylation dynamics during ectopic TE differentiation.

(A) Circular heatmap depicting the dynamics of DNA methylation patterns during TE differentiation. The distribution of the weighted average of CG, CHG, and CHH methylation levels per 100-kb bins on five chromosomes at six time points (0, 12, 24, 36, 48, and 60 hai) during TE differentiation is displayed. The outer annotation tracks depict the relative density of genes and transposons in the 100-kb equidistant window. (B and C) Metagene plots showing the weighted average of the CG, CHG, and CHH methylation levels of genes (B) and transposons (C) around transcription start sites (TSS) and transcription termination sites (TTS). Weighted average methylation levels for each 100-bp interval are plotted.

We next identified differentially methylated regions (DMRs) in a specific stage relative to the preceding stages. In total, we identified 631 CG DMRs, 1591 CHG DMRs, and 5019 CHH DMRs (table S2). For CG DMRs and CHG DMRs, the numbers of hypermethylated and hypomethylated DMRs were comparable. However, most of the CHH-DMRs were hypomethylated. The average CG methylation level of CG-DMRs was slightly decreased at 24 hai and increased at 60 hai (fig. S2A). The average CHG methylation level of CHG-DMRs was moderately decreased at 24 hai and remained low at 60 hai (fig. S2B). The average CHH methylation level of CHH-DMRs was decreased continuously from 12 to 60 hai (fig. S2C). DMRs were distributed mainly on centromeres and pericentromeres, with the distribution of CHH-DMRs being more widespread (fig. S2D). We further divided DMRs into four categories: genic region, intergenic region, region of transposon overlapping a gene, and region of transposon outside of a gene. About half of the CG-DMRs distributed in genic and intergenic regions, whereas CHH-DMRs tended to be located in regions of transposons outside of genes (fig. S2E). The MuDR and RC/Helitron superfamilies of DNA transposons were overrepresented among the CHH-DMR–associated transposons (fig. S2, F and G).

Passive DNA demethylation mediates pericentromeric CHH demethylation, while RDD-dependent active DNA demethylation mediates CG and CHG demethylation at a subset of loci

To determine whether passive or active DNA demethylation was responsible for the removal of DNA methylation during TE differentiation, we profiled the DNA methylome of rdd-2 (a triple mutant of ros1-4, dml2-2, and dml3-2) at 0 and 48 hai. Forty-eight hours after induction was chosen because important TE differentiation genes (e.g., LBD15, MYB85, and XCP2; fig. S1D) were greatly induced and a great number of hypo-DMRs (38 CG hypo-DMRs, 174 CHG hypo-DMRs, and 1756 CHH hypo-DMRs) were identified (table S2) at this time point. Furthermore, TE differentiation was still ongoing and there was no extensive PCD at this time point. We then compared changes in DNA methylation levels of CG, CHG, and CHH hypo-DMRs in the WT Col-0 and rdd-2. The DNA methylation level of CHH hypo-DMRs in rdd-2, like that in Col-0, was decreased at 48 hai (Fig. 2), suggesting that passive DNA demethylation, instead of active DNA demethylation, was largely responsible for DNA demethylation in the CHH context. However, the decrease in the DNA methylation levels of CG and CHG hypo-DMRs at 48 hai observed in Col-0 was abolished in rdd-2 (Fig. 2), suggesting that active DNA demethylation was required for DNA demethylation at CG and CHG hypo-DMRs.

Fig. 2 Involvement of passive and active DNA demethylation in ectopic TE differentiation.

(A to C) DNA methylation levels of CG (A), CHG (B), and CHH (C) hypo-DMRs (48 hai versus 0 hai) in Col-0 and rdd-2. Upper panel: Heatmaps showing the DNA methylation levels of CG (A), CHG (B), and CHH (C) hypo-DMRs at 0 and 48 hai of ectopic TE differentiation. The y axis of the upper panel is clustered hypo-DMRs. Lower panel: Box plots showing the weighted average DNA methylation levels of different types of hypo-DMRs at 0 and 48 hai of ectopic TE differentiation. The cotyledons of Col-0 and rdd-2 seedlings were used for ectopic TE differentiation. Significant differences between two groups are marked with different letters (P < 0.001, Mann-Whitney U test).

Active DNA demethylation is required for proper TE differentiation

Active DNA demethylation could regulate TE differentiation by mediating DNA demethylation at specific genomic loci. To test the role of active DNA demethylation in TE differentiation, we measured the rates of TE differentiation in rdd, rdd-2, and the single mutants ros1-4, dml2-3, and dml3-2 using the VISUAL system. While the TE differentiation rate of Col-0 was greater than 85% (Fig. 3, A and H), those of rdd and rdd-2 were approximately 30 to 40% (Fig. 3, B, C, and H). Among ros1-4, dml2-2, and dml3-2, the TE differentiation rate of ros1-4 was reduced the most in comparison with that of Col-0 (Fig. 3, D to F and H). The defective TE differentiation of ros1-4 was complemented by the ROS1 transgene under the control of a native promoter (Fig. 3, G and H). Our data suggest that, in the VISUAL system, active DNA demethylation is required for normal TE differentiation and ROS1 is the dominant DNA demethylase that regulates TE differentiation.

Fig. 3 Phenotypic analysis of ectopic TE differentiation in mutants for active DNA demethylation.

(A to G) Images showing ectopic TE differentiation in Col-0 (A), rdd (B), rdd-2 (C), ros1-4 (D), dml2-3 (E), dml3-2 (F), and ros1-4 complemented with the ROS1 transgene. (G) Whole cotyledons from different genotypes at 72 hai were imaged by ultraviolet-excited fluorescence microscopy. Scale bar, 500 μm. (H) Quantification of TE differentiation rates in (A) to (G). Error bars indicate SD (N = 10 to 15). Significant differences between the mutant and Col-0 groups are indicated by asterisks (Tamhane T2 test, a = 0.05).

We next investigated whether DNA demethylase dysfunction causes defective TE differentiation in vivo. The ros1-4, rdd, and rdd-2 mutants did not exhibit obvious abnormal developmental phenotypes. However, while metaxylem and protoxylem vessel formation was normal in the young roots of Col-0 seedlings, protoxylem discontinuities were observed in ros1-4, rdd, and rdd-2 (Fig. 4, A to D). More than 70% of the roots of rdd and rdd-2 seedlings had protoxylem discontinuities. Roots with more than three protoxylem discontinuities accounted for 22 and 17% of the total roots of rdd and rdd-2 plants, respectively (Fig. 4F). Approximately 30% of the roots of ros1-4 seedlings had protoxylem discontinuities, and these roots each had one or two discontinuities (Fig. 4F). The ROS1 transgene complemented the phenotype of protoxylem discontinuities in ros1-4 (Fig. 4, E and F). These results suggest that active DNA demethylation is required for normal TE differentiation in vivo, and dysfunction of ROS1 is sufficient to result in defective protoxylem formation in young roots.

Fig. 4 Phenotypic analysis of xylem development in mutants for active DNA demethylation.

(A to E) Images showing the phenotype of protoxylem discontinuities in 6-day-old roots of Col-0 (A), rdd (B), rdd-2 (C), ros1-4 (D), and ros1-4 plants complemented with the ROS1 transgene. (E) Upper panel: Differential interference contrast (DIC) images of protoxylem. Lower panel: Ultraviolet-excited fluorescence microscopic images of protoxylem. Protoxylem discontinuities are indicated by red arrows. (F) Frequency of protoxylem discontinuities in the roots of different genotypes. Error bars indicate SD (n = 70). Scale bar, 20 μm.

To provide additional genetic evidence that the phenotype of protoxylem discontinuities is caused by dysfunctional DNA demethylation, we generated the quadruple mutants rdd-2 nrpd1 and rdd-2 nrpe1 using CRISPR-Cas9 technology (fig. S3). The nrpd1 and nrpe1 mutations are frameshift mutations that cause dysfunction of the largest subunits of Pol IV and Pol V, respectively. These mutations were introduced into rdd-2 to abolish RdDM. Both the nrpd1 and nrpe1 mutations rescued the phenotype of protoxylem discontinuities in rdd-2 to a great extent (Fig. 4F). These findings suggest that protoxylem discontinuities are a result of defective DNA demethylation, and active DNA demethylation promotes TE differentiation by antagonizing RdDM activities.

Genes with a broad spectrum of functions are differentially expressed during TE differentiation

To obtain a global view of gene expression during TE differentiation, we applied RNA sequencing (RNA-seq) to Col-0 cotyledons at sequential stages of TE differentiation. Principal components analysis (PCA) revealed that dots representing different biological replicates were clustered together (fig. S4A), suggesting that our experiments had good reproducibility. However, the 0 to 60 hai samples were separated in the PCA plot, suggesting that gene expression changed as TE differentiation progressed (fig. S4A). The most marked changes in gene expression occurred from 0 to 12 hai and from 24 to 36 hai (fig. S4B). Identification of differentially expressed genes (DEGs) revealed 3471, 5338, 7156, 7234, and 7784 genes with more than twofold changes in expression levels at 12, 24, 36, 48, and 60 hai relative to 0 hai (fig. S4C and table S3). At each stage, approximately half of the DEGs were up-regulated, while another half were down-regulated (fig. S4C).

To understand how DEGs coordinate TE differentiation, we used Dynamic Regulatory Events Miner (DREM) 2.0 software to construct a gene regulatory network (36). Genes with similar expression patterns during the time course of TE differentiation were grouped together. The resulting network consisted of 12 groups of genes, designated as paths 1 to 12 (Fig. 5A and table S3). The genes in paths 1 to 4 were strongly induced at the early stages of TE differentiation and had higher expression at later stages. The genes in paths 5 to 6 were only slightly induced at the early stages of TE differentiation. The genes in paths 7 to 12 were repressed, and this repression was marked for the genes in paths 10 to 12 (Fig. 5B). Based on previously identified transcription factor–gene interactions (37, 38), DREM 2.0 assigned putative transcription factors to the paths (table S3). The VND transcription factors were identified as master transcription factors that coordinate the induction and repression of genes in different paths (Fig. 5A). Many transcription factors without previously recognized roles in TE differentiation were found to regulate the expression of genes in different paths (fig. S4D). Gene ontology (GO) analysis was conducted to determine the biological processes associated with each path (Fig. 5C and table S3). Path 1 was greatly enriched in genes involved in the synthetic processes of cellulose, xylan, and lignin, which are three major components required for SCW formation. Paths 2 to 4 were enriched in genes involved in microtubule- and cell cycle–dependent processes. Paths 9 to 12 were enriched in genes involved in response to hormones. Path 5 was enriched in genes involved in DNA methylation and demethylation, in addition to microtubule-, cell cycle–, and ubiquitin-dependent processes (Fig. 5C and table S3).

Fig. 5 DREM analysis of DEGs during ectopic TE differentiation.

(A) DREM model showing 12 groups of coexpressed genes designated as paths 1 to 12. The y axis indicates the log2 (fold change) of the expression level at different time points relative to 0 hai during TE differentiation. The x axis indicates the induction time in hours (hai). Major transcription factors assigned to different DREM paths and the number of genes per path are indicated. (B) Heatmap showing the dynamics of the relative expression of genes as log2 (fold change) relative to 0 hai in different paths during TE differentiation. (C) GO analysis of genes in different DREM paths. The color of the circles represents the −log10 (FDR) value. The size of the circles represents the enrichment score of the corresponding GO term. The enrichment score was calculated as follows: (number of input genes in a GO term/number of input genes)/(number of genes in a GO term/number of total genes).

Active DNA demethylation is required for dynamic gene expression during TE differentiation

To determine whether RDD-dependent active DNA demethylation is required for dynamic changes in gene expression during TE differentiation, we performed RNA-seq on rdd-2 cotyledons at 0 and 48 hai of TE differentiation. Fewer than 200 genes were differentially expressed in rdd-2 in comparison with Col-0 at 0 hai (fig. S5A and table S4). However, more than 1000 genes were differentially expressed in rdd-2 in comparison with Col-0 at 48 hai (fig. S5A and table S4), among which 958 (87.9%) were DEGs identified during TE differentiation. The expression levels of most (94.4%) of the rdd-2 versus WT DEGs were lower in rdd-2 at 48 hai (fig. S5, A and B).

To identify genes that are involved in TE differentiation and are targeted and regulated by RDD, we set out to identify TE differentiation–associated DEGs that had significantly different expression levels in rdd-2 in comparison with Col-0 (at 48 hai) and were associated with hyper-DMRs (rdd-2 versus WT at 48 hai; table S4). We found 290 DEGs that met our criteria (Fig. 6A and table S4). Most of the 290 DEGs were induced during TE differentiation and had persistently high expression levels at 60 hai (Fig. 6B). Some DEGs showed decreased expression as TE differentiation progressed (Fig. 6B). Many genes known to be important for TE differentiation were among the 290 DEGs. For example, XYLEM CYSTEINE PROTEASE 2 (XCP2) was greatly induced at 48 hai in Col-0 but induced to a lesser degree in rdd-2 (Fig. 6C). XCP2, which functions redundantly with XCP1, participates in the micro-autolysis of cellular contents and the formation of hollow TEs for water conduction (39). Increased DNA methylation before and during TE differentiation, presumably due to defective DNA demethylation, in the promoter of XCP2 was associated with compromised induction of XCP2 in rdd-2 (Fig. 6C and fig. S5C). The 290 DEGs also include genes with TE differentiation–independent functions and genes with unknown functions (table S4). For instance, AT3G13630 encodes a hypothetical protein. Similar to XCP2, compromised induction of AT3G13630 at 48 hai in rdd-2 was associated with increased DNA methylation in its promoter (Fig. 6D and fig. S5D). YLS9 encodes a late embryogenesis abundant hydroxyproline-rich glycoprotein. YLS9 had a moderate level of expression at 0 hai and got further induced at 48 hai in Col-0. However, YLS9 expression could barely be detected at 0 and 48 hai in rdd-2. This was associated with increased DNA methylation in the YLS9 promoter (fig. S5E).

Fig. 6 Regulation of gene expression by active DNA demethylation for TE differentiation.

(A) Venn diagram showing the overlap of DMR (rdd-2 versus Col-0)–associated genes at 48 hai, DEGs identified during ectopic TE differentiation in Col-0, and DEGs identified in rdd-2 relative to Col-0 at 48 hai. (B) Heatmap showing the dynamics of the relative expression of 290 genes identified in (A). The rows correspond to 290 genes. The columns correspond to the log2 (fold change) of the gene expression level in rdd-2 versus Col-0 at 48 hai and the log2 (fold change) of the gene expression level at each time point relative to 0 hai in Col-0. (C and D) Expression and methylation levels of XCP2 (C) and AT3G13630 (D) in Col-0 and rdd-2 at 0 and 48 hai. Left: Genome browser snapshots of the DNA methylation levels and relative expression levels of the indicated genes. Right: Validation of the expression levels of the indicated genes by reverse transcription quantitative PCR (RT-qPCR). Error bars indicate SD, n = 3 (***P < 0.01, two-tailed Student’s t test). (E) Working model for the role of active DNA demethylation in TE differentiation. ROS1-, DML2-, and DML3-mediated active DNA demethylation antagonizes RdDM to regulate the expression of genes involved in TE differentiation.

Most of the hyper-DMRs (rdd-2 versus WT at 48 hai) were located in the promoter region (fig. S5F), suggesting that changes in the expression of their nearby genes are very likely due to DNA hypermethylation. For hyper-DMRs located in the promoters, we found that 1179 of 2674 (44.1%) DMRs were overlapping with transposons, suggesting that transposons may dictate distribution of DMRs over the promoters.

Notably, 668 DEGs had different expression levels in rdd-2 in comparison with Col-0, but these genes were not associated with DMRs (Fig. 6A). Multiple transcription factors and SCW-related enzymes known to be important for TE differentiation, including LBD15, MYB85, MYB20, TED6 (TE DIFFERENTIATION-RELATED 6), CCoAOMT1 (CAFFEOYL COENZYME A O-METHYLTRANSFERASE 1), and PRX25 (PEROXIDASE 25) (40), were among these 668 DEGs (table S4). Our results suggest that active DNA demethylation is directly and indirectly required for dynamic gene expression during TE differentiation.

DISCUSSION

TE differentiation is precisely regulated to ensure proper development of xylem tissue. In this study, we found that, during TE differentiation, passive DNA demethylation mediates reduction in CHH methylation at pericentromeric regions, while active DNA demethylation initiated by RDD mediates reductions in CG and CHG methylation at a subset of DMRs. Dysfunction of active DNA demethylation led to low TE differentiation rate in vitro. Protoxylem discontinuities were observed in the young roots of rdd mutant plants, suggesting that RDD-dependent active DNA demethylation is required for proper development of xylem tissues. In plants, few discernible phenotypes were previously found to be associated with the absence of ROS1, DML2, DML3, and their homologs (7, 8). The phenotypes linked to ros1 and rdd mutations in Arabidopsis are limited to overproduction of stomatal lineage cells (15), increased susceptibility to the bacteria Pseudomonas syringae (41) and the fungus Fusarium oxysporum (42), and increased seed dormancy and abscisic acid (ABA) sensitivity (43). The phenotype linked to mutation of ROS1 homolog OsROS1 is limited to thickened aleurone (16), and mutation of SlDML2 only was found to impair fruit ripening (18). Here, we identified a role of RDD-initiated active DNA demethylation in vascular development, suggesting that RDD-dependent active DNA demethylation regulates more developmental processes than previously found.

We examined the dynamics of gene expression during TE differentiation to identify genes targeted and regulated by RDD-dependent active DNA demethylation. The DREM model shows that DEGs can be grouped into 12 paths according to their expression patterns (Fig. 5A). NAC (for NAM, ATAF1/2, CUC2) domain transcription factors VND1 to VND7 were identified as master transcription factors for the activation of genes in paths 1 to 6, in agreement with previous studies (30, 44). Unexpectedly, VND1 to VND7 were also identified as transcription factors mediating gene repression in paths 10 to 12 (Fig. 5A). Whether VND1 to VND7 activate or repress genes may be dependent on protein modifications, binding partners, and the local chromatin environment (45). Knowledge of target genes repressed by VND1 to VND7 and the mechanisms dictating gene activation or repression by these transcription factors should provide important insights into the mechanisms underlying the regulation of gene expression during TE differentiation. In addition to transcription factors known to be involved in xylem tissue formation, many other transcription factors appear to regulate TE differentiation. For instance, WRKY transcription factors are involved in the regulation of gene expression in paths 3, 5, and 8 (fig. S4D). WRKY transcription factors have previously been shown to play crucial roles in SCW formation of Arabidopsis stem pith (46). Additional studies of the role of WRKY transcription factors should provide insights into changes in gene expression during TE differentiation.

The ros1 and rdd mutants display defects in TE differentiation. This could be because (i) active DNA demethylation is not operating in ros1 and rdd during TE differentiation. Loci, which undergo active DNA demethylation in Col-0, had no reduction in DNA methylation in ros1 and rdd during TE differentiation (Fig. 2A). (ii) Some loci, which show no change in DNA methylation during TE differentiation in Col-0, had increased DNA methylation in ros1 and rdd before and during TE differentiation (examples in Fig. 6, C and D). To determine whether RDD-mediated active DNA demethylation during TE differentiation contributes to gene regulation, we first identified RDD-mediated DMRs during TE differentiation (identified hypo-DMRs from WT 48 versus 0 hai and chosen those showing no reduction in DNA methylation from 0 to 48 hai in rdd-2). We found 342 unique DMR-associated genes. We then examined changes in their expression levels in rdd-2 versus WT at 48 hai. We found that 17 genes had significantly down-regulated expression in rdd-2 versus WT at 48 hai (fold change > 2; see also table S4). YLS9 (fig. S5F) was among the 17 genes. Thus, active DNA demethylation during TE differentiation regulates gene expression. Among the 17 genes, we did not find any known TE-related genes. Some of these genes may be involved in TE differentiation but need to be tested. This is the first report that genes involved in TE differentiation (including XCP2, LBD15, PRX25, and MYB20) are regulated by RDD-dependent active DNA demethylation. No single mutation of any of these RDD-dependent xylem-related genes was sufficient to result in protoxylem discontinuities. The xyp1 xyp2 double mutation resulted in protoxylem discontinuities (47). However, we did not detect decreased expression of XYP1 in rdd, excluding the possibility that simultaneous down-regulation of XYP1 and XYP2 gave rise to protoxylem discontinuities in the rdd mutant. Because more than 1000 genes were dysregulated in the rdd mutant, we propose that the phenotype of protoxylem discontinuities was a result of multiple dysregulations. In other words, dysregulation of multiple genes (including VNDs, LBDs, MYBs, and XCP2) known to be involved in TE differentiation in rdd caused this phenotype. It is also possible that the phenotype of protoxylem discontinuities was caused by a single mutation of a gene with an unexplored role in xylem development. Because dysfunction of RdDM can rescue the phenotype of protoxylem discontinuities in the rdd mutant plants (Fig. 5F), the genes associated with the phenotype should be RdDM targets. RDD opposes RdDM to ensure the expression of genes involved in TE differentiation (Fig. 6E).

During TE differentiation, the expression levels of ROS1, DML2, and DML3 remain essentially unchanged (fig. S4E), consistent with the observation that RDD-dependent active DNA demethylation occurs at a subset of genomic loci. Given the constant expression of ROS1, DML2, and DML3, they must be specifically activated or specifically recruited to a subset of genomic loci. A protein complex with histone acetyltransferase IDM1 as the core protein regulates ROS1 targeting at a subset of its target loci (8). However, the factors that coordinate ROS1, DML2, and DML3 targeting for xylem development and other biological processes remain unidentified.

Pericentromeric regions gradually lost CHH methylation during TE differentiation (Fig. 1). Passive DNA demethylation could be responsible for the decrease in CHH methylation at pericentromeric regions. We examined the expression levels of key DNA methyltransferases and regulators of DNA methylation to determine the factors that are down-regulated for passive DNA demethylation. However, none of the tested enzymes or regulators had apparent changes in their expression levels (fig. S4E), with the exception of CLASSY1 (CLSY1), which functions upstream of RdDM and showed slightly reduced expression. These findings suggest that passive DNA demethylation is more likely to be achieved through losses of protein activity or targeting to pericentromeric regions. We found that the expression level of DME was higher than that of ROS1, DML2, and DML3 throughout TE differentiation (fig. S4E). It has been reported that DME is specifically expressed in reproductive cells, where it is responsible for active DNA demethylation (14). However, increasing evidence suggests that DME is ubiquitously expressed (48). DME also functions in somatic cells, where it contributes to plant defense against the fungus F. oxysporum (49). Given these lines of evidence, we propose that a loss of pericentromeric CHH methylation could be alternatively mediated by DME-dependent active DNA demethylation. Future investigations are necessary to identify the pathway responsible for CHH demethylation during TE differentiation and explore the role of CHH demethylation in xylem development.

MATERIALS AND METHODS

Plant materials and growth conditions

The T-DNA insertion mutant lines ros1-4 (salk_045303), dml2-2 (salk_113573), and dml3-2 (salk_056440) were obtained from the European Arabidopsis Stock Centre. The rdd (ros1 dml2 dml3 triple) mutant was generated as described previously (6). The rdd-2 mutant is a triple mutant of ros1-4, dml2-2, and dml3-2 in the Col-0 background. To generate the nrpd1 and nrpe1 mutant alleles in the rdd-2 background using CRISPR-Cas9, 20–base pair (bp) single-stranded guide RNAs (sgRNAs) targeting NRPD1 or NRPE1 (table S5) were cloned into a CRISPR-Cas9 system controlled by an egg cell–specific promoter (50). The constructs were transformed into the rdd-2 mutant background using Agrobacterium tumefaciens GV3101 by the standard floral dipping method. Homozygous mutant plants with the Cas9 transgene out-crossed were used for experiments. For complementation of ros1-4, the genomic DNA of ROS1 containing the ~2-kb promoter region was amplified and cloned into the pCAMBIA1305 vector for plant transformation. Primary transformants were selected on 1/2 Murashige-Skoog plates containing hygromycin (25 mg/liter). Homozygous lines were used for phenotypic analysis.

Induction of TE differentiation via VISUAL

VISUAL using cotyledons was performed as previously described (26, 34) with minor modifications. Briefly, preculture was conducted under continuous light (4000 lux) for 5 days at 22°C, and healthy cotyledons with a long axis of approximately 2 mm were used for TE differentiation in induction medium containing 15 mM bikinin (SML0094, Sigma-Aldrich). Approximately 1000 pieces of cotyledons were collected at different time points (0, 12, 24, 36, 48, and 60 hai) for microscopy examination and DNA/RNA extraction.

Real-time reverse transcription polymerase chain reaction

Fresh tissues collected from VISUAL time course experiments were ground into fine powder in liquid nitrogen. Total RNA was extracted using the RNeasy Plant Mini Kit (74904, Qiagen) and quantified using a NanoDrop spectrophotometer (ND2000Claptop, Thermo Fisher Scientific). Complementary DNA (cDNA) was synthesized from 1 μg of RNA template using 5× All-In-One RT MasterMix (with AccuRT Genomic DNA Removal Kit) (G492, ABM) according to the manufacturer’s instructions. Genomic DNA contamination was eliminated before reverse transcription. EvaGreen qPCR Master Mix-S for Roche CFX96 (ABM) was used for real-time polymerase chain reaction (PCR). The expression levels of selected genes were normalized to that of internal control gene TIP41. The primer sets are listed in table S5.

Confocal microscopy observation and image analysis

The cotyledons collected from VISUAL time course experiments or 6-day-old roots were fixed in ethanol/acetic acid (3:1) overnight and mounted between glass slides and coverslips filled with clearing solution (chloral hydrate/water/glycerol; 8:2:1) for at least 4 hours. To observe ectopic xylem cell differentiation, a 405-nm excitation/520- to 550-nm long-pass emission filter set (LSM 710 NLO, Zeiss) was used to visualize lignin autofluorescence in spiral-patterned SCWs. Full-width views of ectopic xylem cell formation at 72 hai were captured by a spinning disc confocal microscope equipped with a 4× objective lens (UltraVIEW VoX, PerkinElmer). The differentiation rate was calculated as the area of ectopic xylem cells in a cotyledon divided by the total area of the cotyledon blade using ImageJ software (National Institutes of Health, Bethesda, MD, USA).

To observe the xylem vessels of 6-day-old roots, every sampled root was carefully scanned from the root tip to the hypocotyl region by bright-field microscopy, and discontinuous sites were counted. Representative discontinuous sites were photographed using a 63× oil immersion objective lens under ultraviolet light (FV3000, Olympus). Vascular tissues were color-coded cyan to improve image contrast.

RNA-seq and data analysis

Total RNA was extracted from the cotyledons in the VISUAL time course experiments using the RNeasy Plant Mini Kit (74904, Qiagen). Poly-A RNA-Seq library preparation and high-throughput sequencing were performed by BGI Genomics (Shenzhen, China). Sequencing was performed on a BGISEQ-500 platform to generate single-end 50-bp reads. After removing the adapters, clean reads were mapped to the reference genome of Arabidopsis (version: TAIR10, www.arabidopsis.org) based on the latest annotation from Araport11 (51) (www.araport.org) with the options --min-intron-length 20, --max-intron-length 5000. The mapped reads were further sorted, indexed, and compressed by samtools (https://github.com/samtools/samtools). Only uniquely mapped reads were retained for further analysis. To visualize and compare the read coverage from different time points, bam files were converted into bigwig files by bam2wig.py from RseQC (http://rseqc.sourceforge.net/) with the option --wigsum = 10,000,000,000. The average values of three biological replicates were calculated by the mean function from wiggletools (https://github.com/Ensembl/WiggleTools) and further transformed into bigwig format by UCSC’s wigToBigWig.

FeatureCounts (52) was used to calculate the read counts of each gene from Araport11. Differential gene expression analysis was performed in a pairwise manner by DESeq2 (53) under strict criteria (fold change ≥ 2 and q ≤ 0.01). To reveal the reproducibility of three independent biological replicates from each time point, PCA was conducted using the expression levels of the top 1000 genes with the greatest variance.

To investigate the gene expression dynamics and key transcription factors responsible for TE differentiation, we applied DREM (36) analysis to our RNA-Seq data. For the DREM input datasets and model, see table S3. To simplify the DREM model, we excluded the data from 60 hai and only considered gene expression dynamics at 0, 12, 24, 36, and 48 hai. The log2 fold change values of DEG expression (fold change ≥ 2 and q ≤ 0.01) were calculated and plotted. Based on DREM analysis, all DEGs were grouped into 12 distinct paths. GO analyses for the genes from each path were conducted by agriGO V2 (http://systemsbiology.cau.edu.cn/agriGOv2/index.php) under default parameters. The heatmap showing the log2 fold change values for the DEG expression data was generated by ComplexHeatmap (https://github.com/jokergoo/ComplexHeatmap) in R. The enrichment score was calculated as follows: (the number of input genes in a GO term/the number of input genes)/(the number of genes in a GO term/the number of total genes). The enrichment score bubble map was created by ggplot2. The screenshots of typical loci were generated by Integrative Genomics Viewer (IGV) (https://software.broadinstitute.org/software/igv/). Transcription factor–gene interactions were obtained from a previous study (38). These interactions were compiled from two sources: the AGRIS AtRegNet database (https://agris-knowledgebase.org/moreNetwork.html) and DAP-Seq peaks from the Plant Cistrome Database (http://neomorph.salk.edu/dap_web/pages/index.php). Key transcriptional regulators were predicted by DREM and marked in the figures.

WGBS and data analysis

Genomic DNA was extracted from the cotyledons in the VISUAL time course experiments using the Hi-DNAsecure Plant Kit (DP350-03, TIANGEN). Bisulfite treatment, library preparation, and sequencing were performed by BGI Genomics (Shenzhen, China). Sequencing was performed on an Illumina HiSeq 4000 instrument to generate 100-bp paired-end reads.

Adapters and low-quality reads were filtered by fastp (https://github.com/OpenGene/fastp) under the option --length_required 70. Read mapping and extraction of methylation information were performed using MethyPy (54) (https://github.com/yupenghe/methylpy). Briefly, reads were mapped to the Arabidopsis reference genome (TAIR10) by bowtie2 (http://bowtie-bio.sourceforge.net/bowtie2/index.shtml) in MethylPy. Only uniquely mapped reads were retained for further analysis. PCR duplicates were marked by Picard (https://broadinstitute.github.io/picard/). Bisulfite conversion rates were calculated using the unmethylated chloroplast genome as a negative control. A binomial test was performed on each cytosine to identify reliable methylated cytosines. Last, the methylation information of all covered cytosines was reported in an allc file. To increase sequencing coverage, we merged the data from two biological replicates.

Wiggle files were converted from the allc file defined in MethylPy by in-house R scripts and visualized to assess local DNA methylation changes in IGV. The weighted average DNA methylation levels of specific regions of interest were calculated by the add-methylation-level function in MethylPy under the restriction (depth ≥ 5) and further visualized by ComplexHeatmap in R. The chromosome-wide DNA methylation level was displayed in 100-kb bins by ShinyCircos (https://github.com/venyao/shinyCircos) in R. The methylation levels among different groups were compared by the pairwise Mann-Whitney U test with false discovery rate (FDR) correction (55).

To identify DMRs from whole-genome bisulfite sequencing (WGBS) datasets, we compared two WGBS samples using the bins (DMRcaller-B) method from DMRcaller. The DMRcaller-B method split the genome into 200-bp bins and identified differentially methylated bins from all reads with methylated cytosines. Score testing was performed to identify significantly differentially methylated bins. To control the FDR, the P values were adjusted using an FDR correction (55). A bin was marked as significantly differentially methylated if it satisfied the following criteria: (i) The difference between two samples was significant based on the score test for each bin; (ii) for CG, CHG, and CHH methylation, the difference between two samples was greater than 0.2, 0.2, and 0.1, respectively; and (iii) the average number of reads per cytosine was more than four. Any significantly differentially methylated bins that were separated by fewer than 100 bp were merged as a DMR. DMR-associated genes include genes overlapping with DMRs and genes close to DMRs (distance <3000 bp).

SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/6/26/eaaz2963/DC1

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.

REFERENCES AND NOTES

Acknowledgments: We thank Y. Li for valuable input on the manuscript. We thank the National Center for Protein Sciences at Peking University in Beijing, China, for assistance with imaging. Most of the WGBS bioinformatic analyses were performed on the High-Performance Computing platform at the Center for Life Sciences (Peking University). Funding: This work was supported by the National Key Research and Development Program (2016YFD0600103 to X.-Q.H. and 2016YFA0500800 to W.Q.) and the National Natural Science Foundation of China (31570581 to X.-Q.H. and 31970614 and 31522005 to W.Q.). Author contributions: X.-Q.H. and W.Q. conceptualized this study. W. Lin, W. Liang, and R.-Z.H. performed experiments. L.S., X.L., H.H., and W.Q. analyzed the data. W. Lin, L.S., H.F., W.Q., and X.-Q.H. wrote the paper. All authors read and approved the final manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials, and with regard to the genomic data generated for this study, they have been deposited in the NCBI under accession number PRJNA547955. Additional data related to this paper may be requested from the authors.
View Abstract

Stay Connected to Science Advances

Navigate This Article