Genome-wide loss of CHH methylation with limited transcriptome changes in Setaria viridis DOMAINS REARRANGED METHYLTRANSFERASE (DRM) mutants.

The DOMAINS REARRANGED METHYLTRANSFERASEs (DRMs) are crucial for RNA-directed DNA methylation (RdDM) in plant species. Setaria viridis is a model monocot species with a relatively compact genome that has limited transposable element (TE) content. CRISPR-based genome editing approaches were used to create loss-of-function alleles for the two putative functional DRM genes in S. viridis to probe the role of RdDM. Double mutant (drm1ab) plants exhibit some morphological abnormalities but are fully viable. Whole-genome methylation profiling provided evidence for the widespread loss of methylation in CHH sequence contexts, particularly in regions with high CHH methylation in wild-type plants. Evidence was also found for the locus-specific loss of CG and CHG methylation, even in some regions that lack CHH methylation. Transcriptome profiling identified genes with altered expression in the drm1ab mutants. However, the majority of genes with high levels of CHH methylation directly surrounding the transcription start site or in nearby promoter regions in wild-type plants do not have altered expression in the drm1ab mutant, even when this methylation is lost, suggesting limited regulation of gene expression by RdDM. Detailed analysis of the expression of TEs identified several transposons that are transcriptionally activated in drm1ab mutants. These transposons are likely to require active RdDM for the maintenance of transcriptional repression.

INTRODUCTION

DNA methylation is a common chromatin modification in many plant genomes. Cytosine methylation is the result of post‐replication modification that adds a methyl group to the 5′ carbon. Although virtually all plants that have been assessed contain DNA methylation, there are differences in the levels and context‐specific patterns of methylation in different species (Niederhuth et al., 2016). The majority of our knowledge about the molecular mechanisms that control DNA methylation and the functions of DNA methylation are based on studies in Arabidopsis thaliana (Arabidopsis), because of the viability of plants with highly reduced DNA methylation (Law & Jacobsen, 2010; Matzke & Mosher, 2014). However, studies in other plants have suggested differences in the patterns and control of DNA methylation (Niederhuth et al., 2016; Springer et al., 2016).

DNA methylation in plant genomes involves several distinct methyltransferases that create or maintain DNA methylation, and these can be distinguished by the local sequence context (Law & Jacobsen, 2010). CG methylation is often present at high levels in plant genomes and is maintained following DNA replication by the preference of MET1 and orthologous genes for hemimethylated sites (Law & Jacobsen, 2010). CHG (where H = A, T or C) methylation is also quite common and is catalyzed by chromomethylase enzymes in a feed‐forward loop with H3K9me2 (Du et al., 2012; Johnson et al., 2007). CHH methylation occurs at non‐symmetrical genomic sites and requires specific targeting mechanisms. Evidence from Arabidopsis suggests that the RNA‐directed DNA methylation (RdDM) pathway is responsible for much of the CHH methylation (Cao & Jacobsen, 2002b; Stroud et al., 2014; Zemach et al., 2013). In this pathway, RNA polymerase IV and V (PolIV and PolV) generate and utilize 24‐nt small RNAs (sRNAs) to initiate or maintain CHH methylation at their target sequences (Matzke & Mosher, 2014). Most of the RdDM activity is focused on either small transposable elements (TEs) or the edges of longer TEs (Zemach et al., 2013). In Zea mays (maize), the RdDM activity seems to be particularly high at the edges of TEs near expressed genes (Gent et al., 2013; Li et al., 2015). There is also evidence that some CHH methylation in Arabidopsis, particularly the CHH methylation found within internal regions of longer TEs, does not depend upon the DRM genes but instead requires the activity of the CMT2 chromomethylase (Stroud et al., 2014; Zemach et al., 2013).

In plants, the DOMAINS REARRANGED METHYLTRANSFERASE (DRM) genes were identified as putative relatives of the mammalian de novo methyltransferase Dnmt3, with a unique rearrangement for the order of the methyltransferase domains (Cao et al., 2000). There are two putative functional DRM genes in Arabidopsis, DRM1 and DRM2, that are present as tandem duplicates. Evidence suggests that DRM2 is responsible for the bulk of RdDM in Arabidopsis, but most studies utilize the drm1 drm2 double mutant to ensure the complete loss of function (Cao & Jacobsen, 2002b; Stroud, Greenberg, et al., 2013). Arabidopsis also encodes a third DRM‐like gene, DRM3, encoding a non‐functional methyltransferase that plays a role in the targeting or regulation of DRM1 and DRM2 (Costa‐Nunes et al., 2014 ; Henderson et al., 2010 ). Most plant species have several putative functional DRM genes as well as orthologs of the non‐catalytic DRM3. Studies in Arabidopsis have indicated that drm1 drm2 mutants had reduced CHH methylation and were compromised for the silencing of some genes and TEs (Cao et al., 2003; Cao & Jacobsen, 2002a; Chan et al., 2004; Stroud, Greenberg, et al., 2013; Tran et al., 2005). However, there are no substantial developmental or morphological abnormalities in Arabidopsis plants that lack DRM1 and DRM2 (Cao & Jacobsen, 2002a; Chan et al., 2006). Combining the drm mutant with a loss of function for CHROMOMETHYLASE 3 (CMT3) results in significant phenotypic impacts, suggesting the partially redundant control of gene silencing and asymmetric methylation by DRM and CMT genes (Chan et al., 2006; Henderson & Jacobsen, 2008; Stroud et al., 2014). In Arabidopsis drm1 drm2 mutants there are substantial reductions of CHH methylation at many loci and there is a partial reduction of CHG in these same regions, which is completely reduced in drm1 drm2 cmt3 mutants, suggesting the combined control of CHH and CHG methylation by DRM and CMT at these sites (Stroud et al., 2014).

The CHH methylation triggered by RdDM appears to play a limited role in regulating expression (Cao & Jacobsen, 2002b; Stroud et al., 2014). There are certainly some endogenous loci and transgenes that are silenced by DRM or other components of the RdDM machinery. However, the number of genes or transposons that are activated in a drm1 drm2 mutant is limited, and it seems that there is substantial redundancy between DRM and CMT pathways to maintain silencing in Arabidopsis (Stroud et al., 2014). In contrast, Oryza sativa (rice) plants with loss of function in DRM genes exhibit pleiotropic phenotypes. Loss of function for the rice DRM ortholog OsDRM2 results in pleiotropic phenotypes as well as the aberrant expression of some transposons and genes (Hu et al., 2021; Moritoh et al., 2012; Tan et al., 2016). Combining loss of function for OsDRM2 with CMT genes results in drastic changes in methylation and phenotype in rice (Hu et al., 2021).

In this study we sought to understand the roles of the DRM genes in Setaria viridis. Setaria viridis is an emerging model C4 grass that has a relatively small genome and a short generation time (Bennetzen et al., 2012; Brutnell et al., 2010;Mamidi et al., 2020; Thielen et al., 2020). Significant progress has been made to develop genetic and genomic resources for S. viridis, including two high‐quality reference genomes (A10 and ME034V accessions) and pangenome sequence resources for 598 diverse genotypes (Mamidi et al., 2020; Thielen et al., 2020). TEs account for 46% of the genome and could lead to structural variation that impacts traits (Mamidi et al., 2020; Thielen et al., 2020). The proportion of the genome and distribution of the TEs in S. viridis is similar to that observed for rice (40%) but is much less than that observed in the maize genome (85%) (Akakpo et al., 2020; Schnable et al., 2009). To date, epigenomic resources are largely lacking for S. viridis. In this work, we first develop a whole‐genome DNA methylation map for S. viridis accession ME034V and use CRISPR‐Cas9 to create loss‐of‐function alleles for the two putative functional DRM orthologs, DRM1a and DRM1b. Our results provide insight into how DNA methylation is controlled in monocot species, and present evidence of subtle differences in DNA methylation mechanisms in dicot and monocot species.

RESULTS

Isolation of loss‐of‐function alleles

The S. viridis genome has two genes that encode putatively functional orthologs of the Arabidopsis DRM1 and DRM2 genes: Drm1a, Sevir.9G574800 (A10)/Svm9G0069770 (ME034V); and Drm1b, Sevir.9G496200 (A10)/Svm9G0060600 (ME034V). These two genes are located approximately 5 Mb apart on chromosome 9 and represent a duplication event in Setaria that is not observed in maize or rice (Figure S1). The two protein sequences have 67.1% identity and 76.5% similarity. Expression atlas data for both accession A10 and ME034v suggest that Drm1a is much more highly expressed than Drm1b in leaf tissue (Figure S1). This suggests that Drm1a may have more functional relevance but that there is potential redundancy for these two genes. There is also a DRM3‐like gene, Sevir.3G052500 (A10)/Svm3G0006370 (ME034v), that lacks critical residues in the catalytic domain and is unlikely to provide functional methyltransferase activity. This is likely an ortholog of the Arabidopsis DRM3, which encodes a catalytically inert protein that appears to be required as a cofactor for proper CHH methylation at some loci (Costa‐Nunes et al., 2014; Henderson et al., 2010).

A total of three guide RNAs (gRNAs) targeting Drm1a and Drm1b were designed as described previously (Weiss et al., 2020). To generate S. viridis mutant plants with double gene knockouts, the T‐DNA construct (pTW045) expressing Cas9_Trex2 with all three gRNA sequences was transformed through Agrobacterium‐mediated transformation into the transformable S. viridis genotype ME034V (Weiss et al., 2020) (Figure S2; for details, see Experimental procedures). T0 plants with edits at both targeted genes were identified and selected for self‐pollination. Progeny were screened as described by Weiss et al. (2020), and two transgene‐free T1 plants were selected for further propagation: one containing edits at drm1a and drm1b (T1_12‐9) and one plant with wild‐type (WT) alleles (T1_84‐27). To identify homozygous progeny with frame‐shift mutations at both genes, two generations of self‐pollination were performed. T3 plants with edits at both genes were identified (hereafter referred to as drm1ab), which included a 3‐bp deletion in Drm1a that introduces an early stop codon, as well as a 2‐bp and a 6‐bp deletion in Drm1b that results in a frameshift mutation (Figures 1a and S2). The predicted proteins produced by these mutant alleles both lack the critical domains necessary for methyltransferase activity. The drm1ab plants are fully viable (Figure 1b). The drm1ab plants are reduced in stature and have reduced leaf length (Figure 1b,c). In addition, the drm1ab plants exhibit delayed flowering relative to the WT. The severity of the change in stature and flowering time were variable under different growth conditions.

Figure 1

Isolation of loss‐of‐function alleles for DRM genes in Setaria viridis.

(a) Sequencing of transgene‐free plants derived from transgenic parents expressing gRNAs targeting the Drm1a and Drm1b genes identified individuals that are homozygous for mutations at both target genes. The schematic indicates the position and sequence change at each locus.

(b) Images showing wild‐type (WT) ME034V and drm1ab double mutant plants grown in growth‐chamber conditions. Scale bar: 9 cm.

(c) Leaf length is reduced in glasshouse‐grown drm1ab mutant plants relative to WT controls or the progeny of plants derived from tissue culture. Letters indicate significant differences (P < 0.05). [Colour figure can be viewed at wileyonlinelibrary.com]

Characterization of methylation domains within the S. viridis genome

Whole‐genome DNA methylation profiles were generated for a single replicate sample of WT S. viridis ME034V as well as three biological replicates from plants whose parent (T1_84‐27) was regenerated from tissue culture and three biological replicates of transgene‐free drm1ab plants. All samples were collected from seedling leaf tissue at a developmental stage in which there are no phenotypic differences between the mutant and WT plants. Enzymatic conversion rates (based on alignments with the chloroplast genome) for all samples ranged from 99.43 to 99.81%. The genome‐wide DNA methylation levels for WT ME034V plants (Figures 2a and S3) are quite similar to the reported levels for S. viridis accession A10 (Niederhuth et al., 2016). Prior studies have suggested some changes in DNA methylation induced through tissue culture in rice and maize (Han et al., 2018; Stroud, Ding, et al., 2013). We did not observe significant differences in the overall DNA methylation levels between WT and tissue culture‐derived samples (Figure S3). The genome was divided into 100‐bp tiles, each of which was classified based on the levels of CG, CHG and CHH methylation (Figure S3; for details, see Experimental procedures). The WT and tissue culture‐derived plants had very similar proportions of the genome classified as high CG and CHG (approx. 28% of genome), CG only (approx. 15% of genome) or high CHH (approx. 1.4% of genome) (Figure S3). For the analysis of drm1ab we focused on contrasts between the three biological replicates of tissue culture‐derived plants and the loss‐of‐function lines to ensure that the differences detected would not be solely linked to tissue culture‐induced changes in methylation.

Figure 2

DNA methylation changes in drm1ab plants.

(a) Genome‐wide mean mCG, mCHG and mCHH levels were assessed in three biological replicates of non‐edited tissue‐culture control plants and drm1ab plants. Asterisk indicates significantly lower levels of mCHH methylation in drm1ab.

(b) Metaplots of mCG, mCHG or mCHH levels in genic regions. Solid lines show the profile in tissue‐culture control plants, dashed lines show the levels in drm1ab. A different scale is used for mCHH because of the overall lower methylation levels in this context. Genes are all oriented 5′ → 3′ and the dashed lines indicate the genic region normalized to the same length. The region to the left or right of the dashed vertical lines include 1 kb of upstream or downstream sequences.

(c) Similar metaprofiles of mC levels within and surrounding structurally annotated transposable elements. Genome viewer snapshots showing an example of (d) DRM‐dependent mC loss and (e) DRM‐independent maintenance of mC. [Colour figure can be viewed at wileyonlinelibrary.com]

Widespread loss of CHH methylation in drm1ab mutant plants

Context‐specific DNA methylation levels were evaluated genome‐wide (Figure 2a) and using metaprofiles over genes or TEs (Figure 2b,c). This revealed a substantial loss of CHH methylation in the drm1ab mutant relative to the control. CHH methylation is lost in regions that exhibit elevated levels of methylation, including the regions surrounding genes as well as within TEs (Figure 2b,c). However, CHH methylation still remains in the drm1ab plants, especially within TEs. A visualization of several genomic regions revealed that regions of high CHH methylation in the tissue‐culture control plants can be divided into regions that require DRM (DRM dependent) and regions that have CHH methylation that is not dependent upon DRM (DRM independent) (Figure 2d). The proportion of CHH methylation that is lost in drm1ab plants was assessed for genomic regions with varying levels of CHH methylation in the control. The majority (79.7%) of genomic loci with high (>20%) CHH methylation in the control exhibit at least 80% loss of this methylation in drm1ab plants. In contrast, regions with intermediate (5–20%) or low (2–5%) levels of CHH methylation only lose methylation at 43.3 or 11.8% of loci, respectively. Thus, the loss of DRM activity results in a loss of CHH methylation at regions with high CHH methylation but rarely affects the large number of regions with low levels of CHH methylation.

In order to gain a better understanding of CHH methylation changes, and how these are related to CG and CHG methylation, we identified all 100‐bp tiles with ≥20% CHH methylation in the control sample. The methylation levels of these tiles were assessed in drm1ab to identify DRM‐dependent tiles (≥80% methylation loss in drm1ab), DRM‐intermediate tiles (20–80% methylation loss in drm1ab) and DRM‐independent CHH tiles (<20% methylation loss in drm1ab) (Table 1). The CG, CHG and CHH methylation levels in both control and drm1ab plants were evaluated at these regions (Figure 3). The majority (79%) of regions with CHH levels ≥20% in the WT are DRM dependent, with only 4% that are DRM independent (Table 1). We observed that DRM‐dependent CHH methylation is often accompanied by high levels of CG and CHG methylation (Figure 3). In drm1ab, the CHG methylation is lost in the vast majority of these regions and the CG methylation is reduced at some loci but not at others (Figure 3). The regions of the genome where we observed DRM‐independent CHH methylation generally have high levels of CG and CHG methylation that are not dependent upon DRM (Figure 3).

Table 1

Count of DRM‐dependent, DRM‐intermediate and DRM‐independent methylated tiles in all contexts

	DRM dependent(≥80% loss in drm1ab)	DRM intermediate(20–80% loss in drm1ab)	DRM independent(<20% loss in drm1ab)	Total
CG ≥40%	33 146 (2%)	93 268 (6%)	1 520 749 (92%)	1 647 163
CHG ≥40%	43 603 (4%)	148 452 (13%)	986 556 (92%)	1 178 611
CHH ≥20%	22 883 (79%)	5016 (17.25%)	1174 (4%)	29 073

Figure 3

Comparisons of context‐specific mC levels at DRM‐dependent and DRM‐independent loci. The 100‐bp tiles classified as methylated (≥40% for CG or CHG, ≥20% for CHH) were examined to see whether methylation was lost in drm1ab plants. DRM‐dependent tiles lost >80% mC, whereas DRM‐independent tiles lost <20% mC in drm1ab. Each violin plot includes context‐specific DNA methylation levels in both mutant and control plants for the set of DRM‐dependent or ‐independent tiles. Shaded boxes highlight the mC context used to select the subset of tiles included in each plot. Details of tile classification can be found in Table 1. [Colour figure can be viewed at wileyonlinelibrary.com]

DRM is expected to function in the RdDM pathway and to primarily contribute to the maintenance of CHH methylation. However, genome‐wide levels of CG and CHG methylation show 2–6% reductions (Figure 2a,b). Similar findings were reported for rice DRM mutants (Hu et al., 2021). Methylation levels in all contexts were determined at each 100‐bp genomic tile in the control and drm1ab samples. Differentially methylated tiles were classified as either hypermethylated (higher methylation in drm1ab) or hypomethylated (lower methylation in drm1ab) (Figure S4). In all three contexts there are many more examples of hypomethylated tiles in drm1ab relative to the control (Figure S4). All genomic regions with high (≥40%) CG or CHG methylation levels in the control sample were evaluated to determine what proportion are DRM dependent (Table 1). In contrast to regions with high CHH methylation, the majority of which are dependent on DRM, only a small proportion of these high CG or CHG methylated regions are affected in drm1ab. However, as the number of genomic regions with high CG or CHG methylation vastly outnumbers the regions with elevated CHH, there are more total CG or CHG hypomethylated tiles genome wide (Table 1).

The drm1ab mutant might be expected to affect CG and CHG methylation in regions with high CHH methylation, but we did not expect substantial changes in CG and CHG methylation at regions without CHH methylation. We sought to determine whether the changes in CG and CHG methylation in drm1ab co‐occurred with changes in CHH methylation, or whether some CG and CHG methylation losses occurred in regions without CHH methylation. We found many examples of CG or CHG DRM‐dependent hypomethylation at tiles with low or no CHH methylation (Figure 3). A comparison of the regions with CG and CHG methylation loss found some examples of dual loss in both contexts as well as many with specific loss in CG or CHG (Figure 4a). The CG, CHG or CG/CHG methylation losses were then evaluated to assess what proportion overlap, are within 300 bp of a tile or are farther than 300 bp from a tile with reduced CHH methylation (including both DRM‐dependent and CHH DRM‐intermediate tiles). A portion of the DRM‐dependent CG (9.7%) or CHG (24%) methylation is found in regions that have moderate or high levels of CHH that require DRM (Figure 4a). In addition, another 4–10% of the CG‐ or CHG‐dependent methylation occurs within 300 bp of a region with CHH methylation loss in drm1ab (Figure 4b). This suggests that the loss of a small region of CHH methylation can be associated with a broader loss of CG and/or CHG methylation at some loci (example in Figure S5). The mechanisms leading to losses of CG or CHG methylation in drm1ab in regions distal to CHH methylated regions are unclear.

Figure 4

Relationship of DRM‐dependent mCG and mCHG losses to mCHH hypomethylated tiles.

(a) Venn diagram of the overlap of tiles with mCG and/or mCHG DRM‐dependent methylation with tiles with mCHH DRM‐dependent or ‐intermediate loss.

(b) The proportion of mCG, mCHG or mCG/CHG DRM‐dependent tiles that overlap, are proximal to (within 300 bp) or distal to (greater than 300 bp) mCHH hypomethylated tiles. The number of hypomethylated tiles in each mC context can be found in Figure S4. [Colour figure can be viewed at wileyonlinelibrary.com]

Limited changes in gene expression in mutants

Transcriptome profiling was performed using RNA‐seq on the same seedling tissue samples used for whole‐genome bisulfite sequencing (WGBS), with the addition of multiple replicates for WT ME034V. Principal component analysis (PCA) showed limited variation between WT ME034V and tissue culture‐derived ME034V plants, and, as expected based on the PCA result, there are very few differentially expressed genes between these samples (Figures 5a and S6). However, PCA clustering and differential gene expression analysis finds evidence for hundreds of gene expression changes in drm1ab plants (Figure 5a,b). The drm1ab plants have more genes that are upregulated compared with the control samples, including 136 genes with a greater than 10‐fold upregulation (Figure 5a,b). The observed changes in expression in drm1ab plants may represent direct effects of changes in DNA methylation on gene expression or could represent secondary effects arising from a small number of direct targets that influence the expression of other genes. In order to identify potential direct effects of the loss of RdDM, we initially focused on the subset of genes with CHH methylation immediately over the transcription start site (TSS) region. We identified 1043 genes with CHH methylation (>20%) in the region immediately surrounding the annotated TSS in tissue‐culture control plants. Many (529) of these genes that contain high CHH methylation immediately at or surrounding the TSS are expressed in the control, suggesting that the presence of CHH at or near the TSS is not necessarily silencing gene expression. Although the vast majority of these genes (98%) are hypomethylated in drm1ab, only 3.5% are differentially expressed (24 upregulated in drm1ab and 12 downregulated in drm1ab) (Figure 5c). Over 95% of the genes with elevated CHH methylation surrounding the TSS do not exhibit changes in expression in drm1ab. These observations suggest that there are relatively few genes that are direct targets for silencing by RdDM in seedling leaf tissue of S. viridis.

Figure 5

Transcriptome changes in drm1ab plants.

(a) The number of genes with significant differences in expression (P adj < 0.05; >twofold‐change) was determined for all contrasts.

(b) A volcano plot showing magnitude and P adj for differentially expressed genes in the drm1ab:tissue‐culture comparison. Significant differences are indicated using blue and orange data points.

(c) The proportion of genes that are up‐ (orange) or downregulated (blue) is shown for all genes, the subset of 5424 genes that contain a tile of >20% CHH within 1 kb of the TSS (CHH‐Island genes) or the 1043 genes with >20% CHH methylation in the 100‐bp tile directly over the TSS or the two adjacent tiles (CHH‐TSS genes). [Colour figure can be viewed at wileyonlinelibrary.com]

Previous analysis of DNA methylation in monocots has shown that methylated CHH regions (mCHH islands) are often found upstream of highly expressed genes (Gent et al., 2013; Li et al., 2015; Niederhuth et al., 2016). It is unclear whether the mCHH islands influence nearby gene expression or whether open chromatin, associated with gene expression, enables the DRM‐dependent methylation of these regions. We sought to determine whether losses of CHH methylation at mCHH islands in drm1ab resulted in changes in expression of these genes. We classified 5424 genes as having an mCHH island (requires at least one 100‐bp tile with >20% CHH methylation in the 1‐kb promoter region). Many (3171) of these genes are expressed in control conditions and these genes tend to be more highly expressed than genes without mCHH islands, as previously observed in maize (Gent et al., 2013; Li et al., 2015). The genes containing mCHH islands exhibit only slightly higher proportions of DEGs than in all genes and there are similar proportions of up‐ and downregulated genes in drm1ab relative to the control (Figure 5c). Over 95% of the genes with an mCHH island do not show altered expression in the drm1ab plants, suggesting that the presence of mCHH islands in gene promoters has limited functional significance for gene expression levels in seedling leaf tissue.

Identification of TEs that are upregulated in drm1ab mutants

RNA‐directed DNA methylation (RdDM) has been shown to play important roles in maintaining the silencing of TEs in Arabidopsis (Stroud, Greenberg, et al., 2013; Cao et al., 2003; Cao & Jacobsen, 2002a; Chan et al., 2004; Tran et al., 2005). However, in Arabidopsis and other plants it has been shown that there is often redundant control of TE silencing through multiple DNA methylation pathways and the sole loss of CHH methylation only results in limited TE activation (Chan et al., 2006; Henderson & Jacobsen, 2008; Stroud et al., 2014). We sought to investigate whether there are Setaria TEs that are transcriptionally activated in the drm1ab plants. The Extensive de novo TE Annotator (EDTA) (Ou et al., 2019) pipeline was used to perform both a structural and a homology‐based annotation of TEs in the S. viridis ME034V genome. Two distinct approaches were used to monitor the expression of TEs. The first approach assessed uniquely aligned RNA‐seq reads that align with regions annotated as TEs using the structural annotation of intact ME034V TEs. There were 454 TEs with detectable expression (more than five uniquely mapping reads) in either tissue‐culture or drm1ab samples and 33 of these were differentially expressed (P adj < 0.05 and greater than twofold change) between tissue‐culture and drm1ab samples (Figure 6a). The majority (25/33) of the differentially expressed TEs were upregulated in drm1ab, and 14 of these had little or no expression in the control plants, indicating the requirement for RdDM to maintain the effective silencing of these TEs. The upregulated TEs include 11 class‐I retrotransposons as well as 14 class‐II terminal inverted repeat and helitron DNA transposons (Table S1). It is worth noting that this approach was useful to identify TEs that are upregulated but has two significant limitations. First, the reliance on uniquely mapping reads limits the detection of TE families with multiple highly similar elements, which might be common in TE families regulated by RdDM. Second, in many cases the expression that was detected for TEs is likely to reflect partial transcripts rather than expression of the full‐length TE.

Figure 6

Expression of structurally intact transposable elements (TEs) and analysis of DST‐1‐like TEs.

(a) A volcano plot showing magnitude and P adj for differentially expressed structurally annotated TEs in the tissue‐culture control:drm1ab comparison. Significant differences are indicated using blue and orange data points.

(b) A maximum‐likelihood relatedness tree generated from an alignment of the putative transcripts of the 14 DST1‐like TEs. DSTs that overlap with a structurally annotated TE, have expression in control and/or drm1ab plants are indicated.

(c) Genome browser snapshot of the DST1‐1 region showing mC levels in tissue‐culture and drm1ab plants, (d) RNA‐seq coverage from tissue‐culture control and drm1ab plants, and (e) the location of the predicted DST1‐1 transcript and TE‐associated domains identified by CD‐BLAST. [Colour figure can be viewed at wileyonlinelibrary.com]

The second approach to monitor TE expression was implemented using a de novo transcriptome assembly of the RNA‐seq reads from the WT and drm1ab plants. This enables the identification of TEs that generate potential full‐length transcripts. De novo transcriptome assembly does not rely upon alignment with the reference genome and therefore can potentially identify transcripts that arise from repetitive sequences. The RNA‐seq reads from the WT, tissue‐culture WT and drm1ab samples were aligned with de novo assembled transcripts that were greater than 1 kb in length to identify 103 upregulated transcripts in drm1ab (minimum twofold change and P adj < 0.05). These transcripts include both gene and TE sequences. To focus on putative TEs, we removed any of the 103 drm1ab upregulated transcripts with >50% overlap of an annotated gene based on the alignment of the transcripts with the genome. There were 29 upregulated transcripts that do not align with annotated genes. The analysis of conserved domains within the putative open reading frames (ORFs) of these transcripts identified five of these transcripts that contain TE‐associated domains and three additional transcripts that overlap with structurally annotated TEs. We refer to these eight transcripts as DRM‐silenced TEs (DSTs) (Table S2). The eight DSTs were aligned with the genome to identify the best matching genomic sequence and we assessed the presence of long terminal repeats (LTRs)/terminal inverted repeats (TIRs) and target site duplications. Only DST1 contained the intact structural features that would be necessary for active transposition, and we focused on further characterization of this element and related family members.

The DST1 de novo assembled transcript is 8.7 kb in length and is upregulated by more than eightfold in drm1ab (Table S2). Alignment of the de novo assembled DST1 transcript with the genome assembly of ME034V revealed 15 highly similar sequences (Table S3). Four of the DST1 elements are identified as TEs in the structural annotation and another 10 are at least partially annotated as TEs in the homology‐based annotation. These annotated TEs were not identified as differentially expressed, based on alignments of RNA‐seq reads with the genome, probably because of the repetitive nature of the family and the limited number of unique mapping reads. The internal (non‐LTR) sequences typically have 96–98% identity between members of this family with one element (DST1‐12) having lower similarity, indicating that this may be the oldest member of the family. A phylogeny of the DST1 family based on the internal alignable sequences does not reveal any subgroups with very close relationships that would indicate the continuing movement of this set of elements (Figure 6b). The majority of these (14/15) have intact LTR sequences on the 5′ and 3′ end with 90–96% identity of the two LTRs, suggesting that these elements do not represent particularly recent transposition events. However, there is substantial variability in the LTR sequences between different elements of the family. Only DST1‐2/DST1‐11 and DST1‐4/DST1‐6 share similar LTRs. A comparison among all other elements revealed that there is conservation (90–95%) of identity in the first 250–400 bp and in the last 1.2 kb of the approximately 2‐kb LTRs. The middle region of 300–800 bp is highly variable and cannot be aligned among family members. The DST1‐1 element only contains homology for the first 250‐bp region. DST1 is an intriguing LTR family with highly conserved internal sequences but with a variable region within the LTR, reminiscent of observations for the Tnt1 TE family in Nicotiana tabacum (tobacco) (Casacuberta et al., 1997).

Alignment of RNA‐seq reads with the DST1 transcript reveals an approximately eightfold increase in expression in the drm1ab plants. However, many of the aligned reads contain SNPs relative to the DST1 assembled transcript, potentially reflecting expression of multiple family members with slight sequence variation. The transcripts that arise from DST1 elements were assessed to determine whether the activation in drm1ab occurs at a single element or reflects coordinate activation of several members of the family. Based on SNPs within the DST1 sequence there is evidence for the expression of 11 members of the DST1 family and five of the DST1 elements are only detected in the drm1ab mutant (Figure 6b). The DST1 family members are often highly methylated with high CHH levels (>75%) at, or near, the LTRs. Six of the DST1 elements have a loss of CHH methylation at or within 1 kb of the annotated element in drm1ab (Figure 6c), and several of these are the elements that exhibit strong transcriptional activation in the mutant.

DISCUSSION

The DRM methyltransferases are critical for RdDM in plants. However, there are variable consequences for the loss of functional RdDM among different plant species. Although there are some gene and TE expression differences in A. thaliana, there are limited phenotypic differences (Cao & Jacobsen, 2002a; Chan et al., 2006). In contrast, rice and maize mutants lacking functional RdDM exhibit developmental abnormalities (Alleman et al., 2006; Moritoh et al., 2012; Sidorenko et al., 2009). In this study we isolated S. viridis plants with loss‐of‐function mutations in both catalytically active DRM genes. Although there are some phenotypic differences, such as reduced stature, the plants are fully viable with normal inflorescences. This suggests that a functional RdDM pathway is dispensable under standard growth conditions. However, it is worth noting that we have only maintained the drm1ab mutants in the homozygous mutant condition for several generations. One important function of RdDM might be to ensure the faithful maintenance and inheritance of DNA methylation patterns. Loss of fidelity in maintaining heterochromatin may have growing phenotypic consequences after many generations in the absence of RdDM activity.

A detailed analysis of the seedling leaf methylome in drm1ab plants revealed significant changes in CHH methylation, as expected. In particular, we find a loss of CHH methylation at the vast majority of genomic loci with high (>20%) CHH methylation in WT plants. In contrast, many genomic regions that had lower, but detectable, levels of CHH methylation (5–10%) are not changed in the drm1ab mutant. It seems that the RdDM pathway is responsible for the high levels of CHH methylation that are often found at the edges of TEs, especially near genes. In contrast, the lower levels of CHH methylation are often found within larger TEs and this methylation is likely to be the result of CMT2 or other chromomethylases, similar to observations in other plant species (Zemach et al., 2013). We noted that the distribution of CHH methylation in the drm1ab mutant is quite similar to the distribution of CHG methylation, potentially supporting a role of CMT genes in contributing to the remaining CHH methylation.

The analysis of CG and CHG methylation patterns in drm1ab plants revealed some unexpected findings. First, we found that CHG methylation was rarely maintained at loci that had lost CHH methylation. The regions with high CHH methylation typically had CHG methylation in WT plants and both CHH and CHG methylation were lost in the drm1ab mutant. This suggests a widespread failure to maintain CHG methylation at RdDM targets in the absence of functional DRM. Second, we found numerous examples of CG and/or CHG methylation loss in regions that were not located at or near loci with high CHH methylation. It is not clear why functional DRM is necessary for the maintenance of CG/CHG methylation at these loci, but similar results have been reported for the rice OsDRM2 mutant (Hu et al., 2021). One possible explanation is that these regions have elevated CHH at other developmental stages and the loss of active RdDM at that stage results in the loss of CG/CHG methylation, which is observed in leaf tissue. It is also possible that these sites are targets of active DNA demethylation and require RdDM for maintenance.

There are relatively few changes in gene expression in drm1ab plants. We hypothesized that the subset of genes with high CHH methylation over the TSS would be upregulated in drm1ab. However, we found that many of these genes are already expressed in WT plants that contain methylated TSS regions and very few of these genes have altered expression when the methylation is lost. This observation suggests a complex relationship between RdDM activity and gene expression. It is possible that many of these genes are redundantly silenced by RdDM activity as well as CMT2/3 and MET1 maintenance activities. It is possible that the loss of CHH methylation may destabilize silencing but does not lead to activation.

A small set of TEs that are silenced in WT plants exhibit transcriptional activity in drm1ab plants. These TEs appear to require RdDM for full silencing. In particular, we identified a novel TE family that exhibits the coordinated activation of multiple loci in drm1ab mutants. This family will be of particular interest in future studies of potentially active TEs in the Setaria genome.

EXPERIMENTAL PROCEDURES

Plant material and growth conditions

Setaria viridis variety ME034V was used in this study. The tissue‐culture WT control plants and the drm1ab plants were derived from the same T0 transgenic plant as described in the previous study (Weiss et al., 2020). In the process of mutant screening, seed dormancy was broken by incubating freshly harvested seeds at 29°C for 24 h in a solution of 1.4 mm gibberellic acid and 30 mm potassium nitrate (Sebastian et al., 2014). Seeds were then sterilized with 50% bleach for 10 min, rinsed five times with water and then planted on germination media (0.5× MS, 0.5% sucrose, 0.4% Phytagel, pH 5.7; Sigma‐Aldrich, https://www.sigmaaldrich.com). At 6 days after germination, seedlings were transplanted to soil and grown under a 16‐h light/8‐h dark photocycle at 26°C/22°C (day/night) and 30% relative humidity, according to a protocol modified from Mamidi et al. (2020).

Guide RNA design and vector construction

The genomic sequences of Drm1a and Drm1b were obtained prior to the publication of the ME034V genome and were identified by searching the S. viridis A10.1 reference genome with BLAST (Mamidi et al., 2020) using the Phytozome database (https://phytozome.jgi.doe.gov). CRISPR gRNAs were designed to target the conserved domains in each gene using CRISPR (Haeussler et al., 2016). Conserved domains were identified by aligning the coding sequences from S. viridis with orthologs from brachypodium, maize and Arabidopsis. Construction of the T‐DNA construct, pTW045, was described previously using the Golden Gate assembly method (Čermák et al., 2017; Weiss et al., 2020).

T‐DNA transformation and tissue culture

Agrobacterium tumefaciens‐mediated transformation of S. viridis ME034V was performed as described previously with a few modifications (Van Eck et al., 2017; Weiss et al., 2020). Callus initiation was first performed by removing the seed coats and sterilizing seeds with a solution of 10% bleach plus 0.1% Tween for 5–10 min under gentle agitation. Seeds were then placed on callus induction media with the embryos facing upwards at 24°C in light for 1 week and then moved to dark for callus initiation. Embryogenic calli were collected after 4–7 weeks and inoculated with the AGL1 strain harboring the T‐DNA construct pTW045 (Weiss et al., 2020). Inoculated calli were placed on a co‐culture medium and incubated in the dark at 20°C for 5–7 days. Transformed calli were transferred to selection medium with 50 mg L−1 hygromycin for 4 weeks at 24°C. Selected calli were subcultured on plant regeneration media with 20 mg L–1 hygromycin with 16 h of light to allow the growth of transformed shoots. Elongated shoots were transferred to rooting medium with 20 mg L−1 hygromycin. Shoots were transplanted to soil and grown to maturity.

Genotyping and drm1ab identification

The drm1ab and WT tissue‐culture control plants were identified using genomic PCR with restriction enzyme digestion (cleaved amplified polymorphic sequence (CAPS) assay), followed by Sanger sequencing. PCR was performed with GoTaq Green Master Mix (Promega, https://www.promega.com), in accordance with the manufacturer's instructions, with an annealing temperature of 58°C and an extension time of 1 min. Amplicons were then subjected to restriction enzyme digestion using an enzyme that overlaps with the CRISPR‐Cas9 cleavage site. PCR amplicons made with the corresponding primers were subjected to Sanger sequencing. T‐DNA transgene detection was conducted using two methods: genomic PCR amplification of the hygromycin gene that is close to the T‐DNA left border and a luciferase assay to detect the expression of the luciferase reporter gene that is next to the T‐DNA right border. The luciferase assay procedure was conducted using the Bio‐GloTM Luciferase Assay System (Promega), in accordance with the manufacturer’s instructions. All the primer sequences used in the present study can be found in Table S4.

Methylome profiling

DNA was extracted from tissue collected from the third and fourth leaf of 2.5‐week‐old S. viridis plants grown in a growth chamber with 31°C/21°C, 12‐h/12‐h, day/night conditions. For each sample, tissue from three or four plants were pooled prior to cetyl trimethylammonium bromide (CTAB) DNA extraction. In total, seven samples from pooled tissue were converted and analyzed: one from WT plants, three biological replicates of unedited plants regenerated from tissue‐culture and three biological replicates of drm1ab edited plants. The samples were converted for sequencing with the NEBNext® enzymatic methyl‐seq kit (NEB, https://international.neb.com) and sequenced at the University of Minnesota Genomics Center. All samples were multiplexed in a full Novaseq S1 lane with 150‐bp paired‐end sequencing. Sequencing reads were trimmed with trim galore! 0.4.3, powered by cutadapt 1.8.1 (Martin, 2011) and fastqc 0.11.5, and aligned with the ME034V reference genome (Thielen et al., 2020) with bsmap 2.74 using the following parameters: ‐v 5 ‐r 0 ‐p 8 ‐q 20 (Xi & Li, 2009). Sample conversion rates were calculated for each sample based on the ratio of predicted unconverted to converted cytosines in the unmethylated chloroplast genome. Methyl‐seq alignment metrics for each sample are summarized in Table S5.

The bsmap alignments for the three replicates from the tissue culture and drm1ab were each merged using samtools merge (Danecek et al., 2021). The genome was divided into adjoining 100‐bp tiles and, using the merged data, each tile was classified as: ‘missing data’ (including ‘no data’ and ‘no sites’), ‘CHH > 15%’, ‘CG/CHG’, ‘CG only’, ‘unmethylated’ or ‘intermediate’ following the classifiers and hierarchy outlined by Crisp et al. (2020). The ratio of tiles in each category is displayed in Figure S2. Tissue culture and drm1ab replicates showed little variance, with R 2 of CG methylation values greater than 0.975 within each treatment. Significantly hypomethylated and hypermethylated 100‐bp tiles were identified from the merged replicates using the DSS library (Park & Hao, 2016) in r using a false‐positive adjusted P‐value of 0.05 and a twofold change as cut‐offs (Figure S4).

We further classified tiles as DRM dependent, DRM intermediate or DRM independent. Tiles classified as ‘missing data’ in either the tissue‐culture or the drm1ab mutant samples were omitted from this analysis. For all remaining tiles, we first asked whether the tile was methylated in one or more contexts (≥40% for mCG or mCHG, ≥20% for mCHH) in the tissue‐culture control plant. If yes, the percentage of methylation loss in drm1ab was determined ((mC tissue control − mC drm1ab)/mC tissue control) × 100. A loss of >80% methylation was classified as DRM dependent, a loss of 20–80% was classified as DRM intermediate and a loss of <20% was classified as DRM independent (Table 1).

In order to determine the relationship of CHH methylation with gene expression, all genes were classified as either CHH‐TSS genes and/or CHH‐island genes. CHH‐TSS genes have a tile with >20% mCHH in the tissue‐culture samples that is overlapping, or within one tile, of the annotated TSS. CHH‐island genes have at least one tile with >20% mCHH in the tissue‐culture sample that is between 100‐ and 1000‐bp upstream of the TSS. It is possible for a gene to be both a CHH‐TSS gene and a CHH‐island gene. In addition to association with gene expression, we were interested in whether DRM‐dependent CG and/or CHG methylated tiles were often found at or near CHH DRM‐dependent and CHH DRM‐intermediate tiles. For this, we classified each CG, CHG and CG/CHG DRM‐dependent hypomethylated tile as overlapping, within 300 bp (proximal) or more than 300 bp (distal) from a CHH DRM‐dependent/intermediate tile (Figure 4).

Transcriptome profiling

Prior to DNA extraction, a portion of the ground tissue used for methylome profiling was saved for RNA extraction with the QIAGEN RNeasy mini kit (QIAGEN, https://www.qiagen.com). For RNA‐seq, two additional WT biological replicates were included. RNA was submitted to the UMGC facility for 150‐bp cDNA paired‐end library preparation and run on the Illumina NovaSeq 6000 (Illumina, https://www.illumina.com). Sequencing reads were trimmed as described in the above methylation profiling section. The ME034V genome was indexed using the –runMode geneomeGenerate command of star 2.7.1 (Dobin et al., 2013) using an annotation file that included the primary transcript for each gene as well as all structural TEs. The trimmed reads were aligned with the indexed genome using the –quantMode GeneCounts feature of star 2.7.1 (Dobin et al., 2013). RNA‐seq metrics for each sample are summarized in Table S6.

Read counts were imported into r 4.1.1 (R Core Team, 2021). Normalization (median of ratios) and differential expression was determined using deseq2 (Love et al., 2014). A gene or structural TE was determined to be differentially expressed if the absolute value of the log2 fold change was greater than 1 and the adjusted P‐value was less than 0.05. enhancedvolcano (https://github.com/kevinblighe/EnhancedVolcano) was used to visualize differentially expressed genes and TEs.

Additional RNA‐seq data sets were used to compare our observed expression ratios of Drm1a to Drm1b and Drm3 with previously published data (Figure S1). The data for S. viridis cultivar A10 were downloaded from the Phytozome 13 portal (Goodstein et al., 2012). The data for ME034v are from Thielen et al. (2020).

De novo transcript assembly and analysis

A de novo transcriptome assembly was generated using pooled trimmed RNA‐seq reads from all nine samples with trinity 2.10.0 (Haas et al., 2013). The minimum contig length was set to 200 bp. The de novo transcripts were indexed with the ‘gmap_build’ command of gmap 2015‐09‐26 (Wu & Watanabe, 2005). Default ‘gmap’ parameters were used to map the de novo transcripts back to the ME034V reference genome.

Next, the RNA‐seq reads from each sample were mapped to the transcriptome assembly with salmon 1.2.1 (Patro et al., 2017) using default parameters in order to determine transcripts per million (TPM) in individual biological samples. Differential transcript expression was determined as previously described using the salmon TPM data as the input.

To find transcripts from the transcriptome assembly that may represent TEs that are upregulated in the drm1ab mutant, we first filtered our DE transcript list to only include transcripts greater than 1 kb, reasoning that any shorter transcripts would not encode functional TEs. Next, we set a DE threshold of twofold upregulation in the mutant with an adjusted P‐value of <0.05. Transcripts with at least 50% of their length annotated as genes based on overlap of coordinates from the gmap alignment and the current ME034V gene annotation were removed. The remaining transcripts were subjected to conserved domain BLAST (Lu et al., 2020). Transcripts without TE‐associated domains (e‐value < 0.01) were omitted from further analyses. Finally, the transcript plus the 3‐kb flanking sequence on either side was submitted for a self versus self BLAST to determine whether putative LTR or TIR sequences could be observed.

Of the five transcripts that had evidence for upregulation, one had detectable LTR sequences (now referred to as DST1). DST1 family members were identified in the ME034V genome using BLAST (Altschul et al., 1990) with parameters ‘blastn ‐perc_identity 75 ‐qcov_hsp_perc 75’. A total of 15 similar sequences were identified and classified as a TE family using the 80–80–80 rule (Seberg & Petersen, 2009). To determine the phylogenetic relationship between the DST1 copies, the internal element sequence of the DST1 elements was first trimmed with trimal (Capella‐Gutiérrez et al., 2009) with parameter ‘‐automated1’. Trimmed internal sequences were aligned with muscle (Edgar, 2004), using default settings with manual inspection, and a tree was generated with raxml (Stamatakis, 2014) with settings ‘‐m GTRGAMMA ‐p 12345 ‐x 12345 ‐# autoMRE’. To examine the expression of individual DST1 elements, diagnostic SNPs were identified for each member of the DST1 family from the multiple sequence alignment and expression per individual element was quantified by requiring a minimum of four RNA‐seq reads supported by a diagnostic SNP for each sample.

Annotation of TEs

The ME034V repetitive elements were previously identified using a homology‐based repeat masking approach (Thielen et al., 2020). We were interested in monitoring potentially active TEs that have intact structural elements, and therefore performed a TE annotation using edta (Ou et al., 2019). This approach was implemented using edta 1.9.6 with ‘‐‐species others ‐‐sensitive 1 ‐‐anno 1’ and all remaining parameters set as the default. This produced an initial structural annotation of 9459 intact elements. Simple repeats (‘target_site_duplication, ‘repeat_region’, ‘long_terminal_repeat’) were removed, resulting in a filtered structural annotation of 6369 intact elements accounting for 16.9 Mb (referred to as ‘structural annotation’). These structural elements were then used for a homology search, which identified an additional 188 707 elements (115 Mb) that have similarity to structural TEs but lack intact structural features.

CONFLICT OF INTEREST

The authors declare that they have no conflicts of interest associated with this work.

Table S1. Significantly differentially expressed structurally annotated TEs.

Table S2. TRINITY de novo transcripts that passed filters for: length, expression and overlap % of annotated genes.

Table S3. DST1‐like elements identified in the ME034V genome assembly.

Table S4. Oligos used for genotyping drm1ab plants.

Table S5. Enzymatic methyl‐sequencing metrics.

Table S6. RNA sequencing metrics.

Click here for additional data file.

Figure S1. DRM genes in Setaria viridis.

Figure S2. Genome editing reagents and genotyping data.

Figure S3. Comparisons of DNA methylation levels in drm1ab and unedited plants.

Figure S4. Methylation changes in the drm1ab edited line.

Figure S5. DRM‐dependent loss of mCG and mCHG at regions demarcating edges between high and low mC levels.

Figure S6. A principal component analysis was used to compare expression profiles of the drm1ab mutant with the wild type and tissue‐culture controls.

Click here for additional data file.