Chromatin structure and gene expression changes associated with loss of MOP1 activity in Zea mays.

Though the mechanisms governing nuclear organization are not well understood, it is apparent that epigenetic modifications coordinately modulate chromatin organization as well as transcription. In maize, MEDIATOR OF PARAMUTATION1 (MOP1) is required for 24 nt siRNA-mediated epigenetic regulation and transcriptional gene silencing via a putative Pol IV- RdDM pathway. To elucidate the mechanisms of nuclear chromatin organization, we investigated the relationship between chromatin structure and transcription in response to loss of MOP1 function. We used a microarray based micrococcal nuclease sensitivity assay to identify genome-wide changes in chromatin structure in mop1-1 immature ears and observed an increase in chromatin accessibility at chromosome arms associated with loss of MOP1 function. Within the many genes misregulated in mop1 mutants, we identified one subset likely to be direct targets of epigenetic transcriptional silencing via Pol-IV RdDM. We found that target specificity for MOP1-mediated RdDM activity is governed by multiple signals that include accumulation of 24 nt siRNAs and the presence of specific classes of gene-proximal transposons, but neither of these attributes alone is sufficient to predict transcriptional misregulation in mop1-1 homozygous mutants. Our results suggest a role for MOP1 in regulation of higher-order chromatin organization where loss of MOP1 activity at a subset of loci triggers a broader cascade of transcriptional consequences and genome-wide changes in chromatin structure.

Introduction

Transcriptional regulation of gene expression is important for plant development and genome defense against transposable elements and exogenous DNA, such as viruses and transgenes. Since their discovery in the late 1990s (reviewed in ref. 1), small RNAs (sRNAs) have been demonstrated to play a critical role in transcriptional gene silencing (TGS) pathways including the RNA-directed DNA methylation (RdDM) pathway.

In plants, methylation of cytosines is observed at symmetric (CG and CHG) and asymmetric (CHH) sequence contexts and requires the activity of specific DNA methyltransferases. While symmetric methylation is maintained during DNA replication, CHH methylation is established de novo by CHROMOMETHYLASE3 (CMT3) and DOMAINS REARRANGED METHYLTRANSFERASE1 (DRM1) and DRM2 in association with RdDM, and is often considered an indication of RdDM activity (reviewed in ref. 2).

Relative to total cytosines, low levels of genome-wide asymmetric methyl-cytosines are observed in plants. In Arabidopsis, 1.5% of cytosines are CHH methyl-cytosines (compared with 22.3% CG and 5.9% CHG), whereas 2.2% of cytosines in rice are CHH methyl-cytosines (compared with 59.4% CG and 20.7% CHG methylation). Although higher levels of overall methylation have been reported in maize, the levels of asymmetric methylation are significantly lower than symmetric methylation (approximately 5% CHH compared with 86% CG and 74% CHG methylation). While the genome-wide levels are relatively low, in maize CHH methylation is enriched at intergenic loci of euchromatic chromosomes arms. Certain transposon families preferentially insert within 1 kb of maize genes and are hypothesized to drive enrichment of CHH. Near-gene CHH methylation is guided by 24 nt siRNAs and correlates with TE insertion events.,

Recent evidence in Arabidopsis thaliana suggests the existence of more than one siRNA-dependent DNA methylation TGS pathway in plants that function to establish/maintain silent chromatin states. One such pathway requires a plant-specific DNA-dependent RNA polymerase, RNA polymerase IV (Pol IV), for siRNA biogenesis and is therefore termed Pol IV-RdDM. Pol IV transcripts are used as templates by the Pol IV-associated RNA DEPENDENT RNA POLYMERASE2 (RDR2), and processed into double-stranded RNA (dsRNA). A Dicer-like protein, DICER-LIKE3 (DCL3) cleaves RDR2-generated dsRNA into 24 nt siRNAs that are incorporated into an ARGONAUTE4 (AGO4) or AGO6 complex. This AGO4/6-siRNA complex directs transcripts produced by Pol V, an additional plant-specific DNA-dependent RNA polymerase, to target loci, ultimately acting as a guide for recruitment of chromatin remodelers and DNA methyltransferases (reviewed in ref. 7).

The Pol IV-RdDM pathway was initially described in Arabidopsis thaliana and several maize orthologs have been identified through genetic screens using loss of epigenetic silencing as phenotypes. The b1 and pl1 paramutation systems were used to identify the mediator of paramutation (mop) and required to maintain repression (rmr) genes, respectively, and several MOP and RMR proteins are required for both paramutation and TGS (reviewed in ref. 10). rmr6 appears to encode the largest subunit of Pol IV and rmr1 encodes a Snf2-like ATPase similar to DEFECTIVE IN RNA-DIRECTED DNA METHYLATION1 (DRD1) and CLASSY1 (CSLY1), which are two Arabidopsis proteins implicated in RdDM., mop1 encodes a protein that is orthologous to RDR2, whereas mop2/rmr7, are most similar to Arabidopsis NRPD2/E2, which encodes the second largest subunit of both Pol IV and Pol V.

A separate forward genetics screen identified the transgene reactivated (tgr) mutants based upon transcriptional reactivation of the silent B1 genomic transgene (BTG-s), which is also reactivated in mop1-1, rmr1-1, rmr2-1, and mop3-1 (McGinnis and Sloan, unpublished data). Reactivation of BTG-s in mop1, rmr1, rmr2, tgr1, and tgr2 mutants is associated with hypomethylation of the 35SCaMV promoter sequence,, and indicates that at specific regulatory loci, TGR1 and TGR2 likely function in the same epigenetic DNA methylation and transcriptional gene silencing pathway as RMR1, RMR2 and MOP1.

Consistent with its putative function in a Pol IV-RdDM pathway, loss of MOP1 results in a dramatic reduction of 24 nt siRNAs., mop1-1 mutants also display pleiotropic developmental defects and phenotypes, which may result from direct and indirect misregulation of genes and siRNAs.,,- These developmental phenotypes are often variably penetrant and expressive, and seem to be influenced by environmental factors,, suggesting an increased potential for variable gene expression in this mutant.

Consistently, many genes and transposable elements were differentially expressed in the shoot apical meristems (SAMs) of mop1-1 homozygous plants. Together, these findings provide evidence that MOP1 not only plays a crucial role in progression of the Pol IV-RdDM pathway but also in maintenance of genome-wide stable epigenetic states that influence gene expression. At least one MOP1-regulated locus also exhibits changes in chromatin structure and nucleosome positioning in association with changes in transcriptional level,, suggesting that RdDM may also influence chromatin structure.

Eukaryotic chromosomes are structurally and functionally organized in the form of chromatin. Nucleosomes, the basic unit of eukaryotic chromatin, are comprised of approximately 147 bp of DNA wrapped 1.65 times around a histone protein octamer (reviewed in ref. 27). In the nucleus, chromatin undergoes several levels of folding, compaction, and organization. The distribution of nucleosomes along DNA, as “beads-on-a-string,” forms the primary structure of chromatin. Nucleosomes interact with each other and give rise to secondary structures, which has been modeled as a 30 nm chromatin fiber. These secondary structures are further condensed to form tertiary (or higher-order) structures of highly condensed chromatin, as seen during the metaphase of mitotic nuclei (reviewed in refs. 27 and 28). The different levels of structural chromatin organization allow for packaging of a significant amount of charged molecules into the confined space of the nucleus while creating a structure that is permissive for essential cellular functions such as replication and transcription.

The mechanisms of higher order chromatin organization and compaction into the nucleus are not well characterized, although several models of this folding have been proposed (reviewed in ref. 27). Epigenetic modifications play a role in chromatin organization. For example, histone acetylation and tri-methylation of H3K4 are often associated with euchromatin, whereas DNA methylation and H3K9 di-methylation are associated with heterochromatin (reviewed in ref. 29). These modifications to DNA and histones are distributed non-randomly along chromosomes to modulate chromatin compaction and regulate transcription by allowing or inhibiting access of DNA to regulatory proteins.

Given the role of epigenetic modifications on higher order chromatin structure and accessibility of DNA to regulatory proteins, it is plausible that genetic mutants of epigenetic factors, such as mop1-1, could also exhibit defects in local and global chromatin organization as an element of a loss-of-function phenotype. Accordingly, mop1-1 homozygous plants were assayed for local nucleosome occupancy at transcription start sites (TSSs) of approximately 400 genes in maize, resulting in the identification of small number of distinct examples where nucleosome occupancy in mutants differed from WT. This suggests that at some loci, MOP1-mediated pathways influence primary chromatin structure.

To date, little is known about how epigenetic mechanisms regulate higher-order chromatin structure. Because chromatin structure is modeled to be hierarchical, disruptions in primary structure may influence changes in higher order chromatin structure. Thus, it is possible that RNA-dependent regulatory mechanisms might also regulate chromatin organization. Herein, we investigated the relationship between higher-order chromatin structure and transcription in response to loss of the epigenetic regulator MOP1. We observed that, on a genome-wide scale, loss of MOP1 resulted in an increase in chromatin accessibility/sensitivity at gene-rich chromosome arms; whereas pericentromeric regions were less accessible/more resistant. We assayed expression, and observed both up and downregulated genes, physically distributed across the genome, indicating more extensive changes in chromatin organization and transcription than would be explained by specific MOP1 epigenetic silencing in localized regions. We identified a subset of differentially expressed genes (DEGs) likely to be direct targets of epigenetic transcriptional silencing. Our results suggest a possible role for MOP1-mediated pathways in regulation of higher-order chromatin organization that results in direct and indirect consequences on gene expression in maize.

Results and Discussion

The mop1-1 mutation results in global changes in nuclease sensitivity

To determine a role for MOP1 on genome-wide chromatin structure, we used a custom maize NimbleGen microarray based platform to assay higher-order chromatin structure similar to that described in humans., Chromatin structure was determined as a measure of sensitivity of formaldehyde-fixed nuclei to enzymatic cleavage by micrococcal nuclease (MNase). MNase cleaves naked/linker DNA (sensitive DNA) more readily than DNA bound in nucleosomes (resistant DNA) (reviewed in ref. 33), and can be used as a probe to assay chromatin structure (reviewed in ref. 34). Nuclei isolated from mop1-1 homozygous mutants and WT immature ears were digested with limiting amounts of MNase and separated by agarose gel-electrophoresis into MNase-resistant and MNase-sensitive. MNase-resistant and MNase-sensitive DNA fractions were differentially labeled and hybridized to a custom NimbleGen microarray containing ~150 000 unique, tiled and isothermal oligonucleotide probes that span the maize genome at ~14 kb spacing. We measured the signal ratio of MNase-resistant DNA/MNase-sensitive DNA for WT and mutant samples at each probe, and then calculated the deviation of the mutant response from WT (Fig. 1).

Figure 1. Chromatin accessibility (nuclease sensitivity) in mop1-1 and WT immature ears on three representative chromosomes. We measured the signal ratio (log2 [HMW MNase-resistant DNA/MNase-sensitive DNA]) for WT and mutant samples at each probe, then calculated the deviation of the mutant response from WT (Mutant – WT). This deviation was plotted on a genome-browser (genomaize.org) for chromosome 1 (A), chromosome 3 (B) and chromosome 9 (C). WT and mutant are the same at zero. Each plot represents the average of 3 biological replicates with smoothing at 5 pixels. The gene density plot (maizesequence.org) and relative position of the centromere (Jiang/Presting) are also illustrated for each chromosome. All ten chromosomes are included in Figures S1 and S2.

A consistent trend in chromatin response as a result of loss of MOP1 function was observed at all ten chromosomes of maize (Fig. 1; Figs. S1 and S2). The assay provides a low resolution illustration of chromatin structure in which chromatin was more accessible to cleavage by MNase at chromosome arms, whereas pericentromeric regions of chromosomes were more resistant (Fig. 1).

This pattern coincided with gene density distribution (maizesequence.org) across maize chromosomes (Fig. 1). In maize, 24 nt siRNAs are also typically enriched at chromosome arms and correspond with CHH methylation in class I TEs at a high frequency., Relative to other transposon families, most families of DNA elements (except the CATCA superfamily) are found more frequently at chromosome arms within and/or adjacent to genes. It is well documented that the transcriptional silencing and DNA methylation of most Mutator (Mu) elements, and transcriptional silencing of most differentially expressed DNA TEs (with the exception of the CATCA superfamily) in SAMs are reversed in mop1-1 homozygous plants.

Together, this provides evidence that MOP1-mediated activity occurs predominantly at chromosome arms, and less frequently at pericentromeres in a pattern that overlaps with specific types of TEs (e.g., Mu TEs) or their proximity to genes. An increase in chromatin accessibility at chromosome arms as a result of loss of MOP1 is consistent with loss of chromatin repression, conceivably resulting in loss of transcriptional repression at affected loci.

Misregulation of gene expression in mop1–1 immature ears

MOP1 is highly expressed in maize immature ears, and developmental and morphological defects are observed in mop1-1 mutants in this tissue, which include late flowering, and ears that are smaller in size and/or narrower in diameter compared with WT., The developmental phenotypes and changes in chromatin organization we observed in mop1-1 plants could be related through changes in gene expression. While gene expression profiles for mop1-1 SAMs are available, to date, no published genome-wide expression profile exists for mop1-1 immature ears. To identify immature ear-specific DEGs in mop1–1 mutants, we used a NimbleGen expression array platform, as described by Sekhon et al., hybridized with cDNA prepared from the same samples used to assay global chromatin sensitivity to MNase. Because a diagnostic phenotype of mop1-1 homozygous individuals is caused by the release of epigenetic silencing of the b1 gene, we anticipated that one hallmark of MOP1 activity would be upregulation of expression in mop1-1 ears compared with wild type.

We defined differential expression as genes with a 2-fold or greater change in expression in mop1 mutants relative to WT. Using this criterion, we identified ~762 DEGs in mop1-1 mutants in maize immature ears (two biological replicates; 2-fold change cut-off (FC ≥ 2); P < 0.05, one-tail t test). We compared our data with the published data set for expression in mop1-1 SAMs and observed that some genes were differentially expressed in both tissues in mutants compared with wild type (Table S3).

Similar to the variations observed in MNase accessibility, this differential effect in expression in immature ears spanned both arms of all chromosomes. Of these, 349 (~46% of DEGs) genes were upregupregulated in mop1-1 and 413 (~54% of DEGs) genes were downregulated in mop1-1 (Tables S1 and S2, respectively). Upregulation of expression in mop1-1 is consistent with reduced RdDM-mediated transcriptional repression. Downregulation of some loci has been previously reported for mop1-1 homozygous individuals, where approximately two-thirds of DEGs were downregulated in mop1 SAMs, although downregulation is not an anticipated direct consequence of the loss of epigenetic silencing.

Upregulation of transposable elements apparent in mop1-1 ear shoots

Approximately 23% (79/349) of the upregulated loci are highly similar to sequences annotated as transposons in Arabidopsis and/or rice (Table S1). This was much fewer than findings by Jia et al. (2009), where ~175 transposon families were differentially expressed in mop1-1 SAMs. It is likely that our technique might be limited in this regard. Repetitive, conserved sequences, such as DNA and retro-transposons, often do not meet NimbleGen microarray probe design criteria, therefore our experiments had significantly fewer of these types of elements represented. Although the majority of these TEs identified herein were class II elements (DNA transposons), some class I TEs (retrotransposons) were also upregulated in plants deficient in MOP1 activity. MOP1 regulation of Mutator, an abundant group of class II elements, has been reported previously.,

Many genes with apparent regulatory function are upregulated in mop1-1

In addition to TEs, many genes predicted to be protein coding were differentially expressed in mop1-1 homozygous plants. Of the 117 upregulated genes with an identifiable Biological Process Gene Ontology term, 19 are predicted to encode proteins involved in regulation of transcription (Fig. 2). The loss of silencing observed at numerous transcription factors would be expected to have secondary consequences in mop1-1 homozygous individuals. We found that several predicted MADS-box transcription factors and related proteins were misexpressed in mop1-1 (Table 1). MADS-box genes regulate reproductive organ identity during floral development and function through interactions with chromatin-associated proteins and other transcriptional regulators, reviewed by. The upregulated genes included those orthologous to the Arabidopsis SUPERMAN (SUP), SEEDSTICK (STK), AGAMOUS-LIKE 8 (AGL8), and AGL26 proteins. Additionally, expression of a maize gene similar to GIANT KILLER (GIK), a direct target of AGAMOUS (Ng KH, Yu et al. 2009), was also upregulated in mop1–1 compared with wild type. GIK is an AT-hook DNA-binding chromatin modifier and regulates patterning and differentiation of reproductive organs in Arabidopsis. Three genes with orthology to ETTIN (ETT)/AUXIN RESPONSE FACTOR3 (ARF3), a predicted target of GIK, were also misexpressed in maize mop1-1 mutants (Table 1). In Arabidopsis, GIK-mediated repression of ETT is associated with chromatin changes (H3K9me2). Upregulation of maize genes with similarity to both GIK and its target ETT in mop1-1 homozygous plants supports a model where loss of MOP1 triggers primary and secondary changes in gene expression.

Figure 2. Biological process of genes upregulated in mop1-1. Gene ontology was predicted and scored for all mop1-1 differentially expressed genes using BLAST2GO (Conesa, Gotz et al. 2005). The most highly scored categories are listed on the y axis, the bars represent the score for each category.

Table 1. Floral Homeotic and Chromatin-associated Genes Upregulated in mop1–1.

Maize Gene ID2	Avg. FC	P value3	Arabidopsis Homolog	GO Term4
Floral Homeotic-associated Genes
GRMZM2G003927; RAMOSA1 (RA1)	4.162	0.0326	SUPERMAN (SUP)	No Term
GRMZM2G001139	3.074	0.0026	SEEDSTICK (STK)	Regulation of transcription, DNA-dependent
GRMZM2G079727	2.492	0.0025	AGAMOUS-LIKE 8 (AGL8)	Regulation of transcription, DNA-dependent
GRMZM5G854901	7.938	0.0047	AGAMOUS-LIKE 26 (AGL26)	Cellular biosynthetic process; Methylation
GRMZM2G072274; BARREN STALK FASTIGIATED1 (BAF1)	2.391	0.0106	GIANTKILLER (GIK)	No Term
GRMZM2G030710;AUXIN RESPONSE FACTOR24 (ARF24)	2.088	0.0002	ETTIN (ETT)/AUXIN RESPONSE FACTOR3 (ARF3)	Regulation of transcription, DNA-dependent
GRMZM2G441325; AUXIN RESPONSE FACTOR23 (ARF23)	2.366	0.0059	ETTIN (ETT)	Regulation of transcription, DNA-dependent
GRMZM2G056120AUXIN RESPONSE FACTOR11 (ARF11)	2.228	0.0062	ETTIN (ETT)	Regulation of transcription, DNA-dependent
GRMZM2G438438	2.351	0.0063	HOMEOBOX PROTEIN 26 (HB26)	No Term
GRMZM2G317160	2.425	0.0054	AINTEGUMENTA-LIKE 5 (AIL5)	Regulation of transcription, DNA-dependent
GRMZM2G327605	2.031	0.0067	HOMEOBOX PROTEIN 33 (HB33)	No Term
Chromatin-associated Genes
GRMZM2G178435	6.331	0.0050	CHR40	No Term
GRMZM2G126096	4.335	0.0339	REGULATOR OF CHROMOSOME CONDENSATION 1 (RCC1)-related	No Term
GRMZM2G360389	3.162	0.0264	SUVR4	Histone lysine methylation
GRMZM2G149708	3.033	0.0061	HIGH EXPRESSION OF OSMOTICALLY RESPONSIVE GENES 15 (HOS15)	Chromatin silencing
GRMZM2G082538	2.126	0.0200	HISTONE DEACETYLASE 2B (HD2B)	No Term

1 Fold change (FC) ≥ 2; 2Gene name from http://www.gramene.org; 3P value: one-tail t test (4 decimals); 4GO term based on biological process.

As previously reported, several other chromatin associated genes were upregulated in mop1-1 (Table 1). The chromatin remodeling40 (chr40) gene encodes a SNF2 type; DRD1 class protein (www.chromdb.org) similar to REQUIRED TO MAINTAIN REPRESSION1 (RMR1), and was upregulated in mop1-1. Expression of a gene related to chromatin-associated factor REGULATOR OF CHROMOSOME CONDENSATION1 (RCC1) was also upregulated in mop1 mutants. Additionally, we found that expression of a gene highly similar to SUVR4, a SET-domain histone lysine methyltransferase, also increased in mop1-1. Proteins encoded by these genes may, in theory, function to maintain a repressed chromatin structure at specific loci in the absence of MOP1-mediated pathways, inducing downregulation at some loci.

Genes downregulated in mop1-1 may not be direct targets of MOP1-mediated regulation

In mop1 immature ears we identified 3 of the 413 (<1%) downregulated loci to encode TEs. As mentioned previously, our platform is limited in its ability to assay repetitive sequences, and so we were unable to conclude if immature ears exhibit a similar pattern of downregulation of retrotransposons as observed in shoot apical meristems. Most downregulated loci observed in our analysis were predicted to be protein coding.

One of the genes downregulated in mop1-1 immature ears is homologous to REPRESSOR OF SILENCING 1 (ROS1)/DEMETER-LIKE1 (DML1) (FC = 2.0485; P = 0.0001; Table S2), an Arabidopsis DNA glycosylase protein that functions in active demethylation.- ROS1 is downregulated in several Arabidopsis RdDM mutants which include pol IVa, IVb, rdr2, and drd1., Hypermethylation and reduced expression can occur when ROS1 function is lost in ros1 mutants. It is possible that reduced ROS1 function leads to the downregulation of a subset of genes in mop1-1 homozygous individuals. Consistent with this, a recent report suggests that RdDM based silencing of some loci might be enhanced in plants deficient for ROS1 activity.

Another possible explanation for downregulated genes in mop1-1 homozygous individuals would be the expanded activity of an alternative silencing pathway in plants. A small RNA-dependent silencing pathway that operates independently of RDR2 activity was recently described in Arabidopsis, and is mediated by 21–22 nt siRNAs.,- In mop1-1 homozygous individuals, while 24 nt siRNAs are substantially reduced, 22 nt siRNAs are abundant. In the absence of the 24 nt siRNA regulatory pathway, perhaps the persisting 22 nt siRNAs downregulate a subset of genes. However, we were unable to detect an immediate correlation between an increase in homologous 22 nt siRNAs and decreased transcript abundance in mop1-1 individuals (data not shown).

Identification of putative targets of epigenetic transcriptional gene silencing

Based upon our observations, we reasoned that loss of MOP1 activity has both primary and secondary consequences, meaning that mop1-1 DEGs would include direct targets of an epigenetic transcriptional gene silencing (TGS) pathway and secondary consequences resulting from genome-wide changes triggered indirectly by the loss of MOP1. To identify genes that are likely direct targets of epigenetic TGS, we used the transgene reactivated1 (tgr1) mutant which releases silencing of some of the same loci as mop1-1., We identified a subset of the 349 genes upregulated in mop1-1 that were also upregulated through loss of TGR1 activity (Table S4) (log2FC ≥ 0.8) using next-generation sequencing of transcripts (RNA-seq) from immature ears of tgr1 plants. Similarly, to identify genes repressed due to loss of epigenetic TGS activity, of the 413 genes downregulated in mop1-1, we identified genes also repressed in tgr1-1 homozygous mutants (Table S5). We found 81 genes that were upregulated (Fig. 3A; Table S4) and 44 genes that were downregulated (Fig. 3B; Table S5) in both genotypes, and considered these to be genes most likely directly regulated by a MOP1-TGR1 epigenetic TGS pathway.

Figure 3. Identification of unifying features in mop1-1 DEG subsets. A subset of genes that were upregulated in mop1-1 were also misregulated in tgr1-1 (A). A different subset of genes were downregulated in both mop1-1 and tgr1 (B). We further categorized each subset to identify those with 24 nt siRNAs within 2 kb upstream of the transcription start site (TSS) (C), a class II transposable element within 1 kb up or downstream of the TSS (D), or both 24 nt siRNAs within 2 kb upstream of the TSS and a class II transposable element within 1 kb up or downstream of the TSS (E).

Upregulation of TGS targets associated with near-gene transposons and 24 nt siRNAs

Analysis of promoter-proximal 24 nt siRNAs at DEGs

In maize, the majority of 24 nt siRNAs are depleted in mop1-1 homozygous mutants, indicating that they are produced in a MOP1-dependent manner. Based on the current model for RdDM, 24 nt siRNAs guide de novo methylation to target loci (reviewed in ref. 7). At maize genes, 24 nt siRNAs peak within 1 kb of transcription start sites and are depleted in mop1-1 mutants (Fig. S3).We predicted that the regulatory regions, including the promoter-proximal regions, of genes directly regulated by MOP1 activity might be expected to share homology with 24 nt siRNAs produced in a MOP1-dependent manner.

Using publically available mop1 siRNA data sets generated in the same tissue as our data, we identified mop1 DEGs with a minimum of 5 mapped 23–24 nt siRNAs detected within 2kb upstream of transcription start sites (TSSs) in either wild type or mutant genotypes; then measured the siRNA fold change between WT and mutant. 142/349 mop1-1 upregulated genes and 195/413 downregulated genes had a ≥2-fold change in 23–24 nt siRNAs within –2kb of the TSS (Figs. 3C and 4A; Tables S6 and S7). 30/81 upregulated and 22/44 genes repressed in mop1-1 and tgr1 had a > 2-fold change in 24 nt siRNAs within –2kb of the TSS (Fig. 3C and 4B; Tables S8 and S9). At genes differentially expressed in mop1-1 immature ears, our analysis indicated that MOP1-dependent siRNAs were consistently more abundant in WT at almost all mop1-1 up- and down-regulated genes (Fig. 3A). The average expression fold change of the subset of MOP1 upregulated genes with 24 nt siRNAs was 4.07 and the P value (0.012) was lower than the P value (0.015) for the average expression fold change of 2.99 observed for MOP1 downregulated genes with 24 nt siRNAs (Fig. 5). This suggests that on average, this subset of genes was more strongly affected by loss of MOP1 activity.

Figure 4. Relative expression and siRNA abundance in -2 kb region of mop1-1 DEGs. Relative expression is expressed as log2 fold change in mop1-1 mutant/wild type (black bars). Log2 fold change of 23–24nt siRNAs (gray bars) -2kb upstream of predicted transcript start site in mutant/wild type for genes differentially expressed in mop1-1 mutants (A) and both mop1-1 and tgr1-1 mutants (B).

Figure 5. Average fold change of expression in mop1–1 DEG subsets. We measured the magnitude of average expression fold change for mop1-1 DEGs within different subsets. One subset included gene differentially expressed in mop1-1 alone and the other included genes differentially expressed in mop1-1 and tgr1. Within each subset, we measured the average expression of upregulated (Up) and downregulated (Down) genes; those Up and Down genes with 24 nt siRNAs within -2 kb of TSSs (+siRNA) and those that have a class II element within + or -1 kb of the TSS (+TIR) (A). The distribution of p-values for each of the genes falling within a particular subset (B), and the averages of these p-values for upregulated (black line) and downregulated (gray line) genes were calculated (C).

We predicted that direct targets of a MOP1-TGR1 epigenetic silencing pathway might be more likely than indirect targets to exhibit a correlation between accumulation of homologous promoter-proximal 24 nt siRNAs and gene expression. However, as observed at genes differentially expressed in mop1-1 alone, the presence or absence of siRNAs at epigenetic (MOP1 and TGR1) targets showed no direct correlation with the direction of transcriptional activity (Fig. 4B). These results are also consistent with observations made at genes upregulated in Arabidopsis RdDM mutants (rdr2-1 and ddc); 24 nt siRNAs were detected and absent upstream of RdDM targets. If we confined our analysis to specific subset of genes using multiple defining features, we did detect that the average expression fold change of genes upregulated in both mop1-1 and tgr1-1 and with promoter-proximal 24 nt siRNAs was 6.5 at a P value of 0.008. This P value was much lower than a P value of 0.013 for the average expression of 3.01 which was observed for the 22 TGS downregulated genes with 24 nt siRNAs (Fig. 5). Genes downregulated in both mutants, with or without 24 nt siRNAs, exhibit a level of expression fold change similar to genes repressed in mop1-1 mutants alone (Fig. 5A).

While the expression of the strongly upregulated genes that appear to be direct targets of a MOP1-mediated pathway varied significantly, the distribution of confidence of expression was higher and had a narrower range than the entire set of mop1-1 DEGs (Fig. 5B) and the average P value decreased when genes were categorized using multiple features (Fig. 5C). This suggests that in mop1-1, expression of the subset of genes considered to be direct targets of TGS was the most strongly and consistently affected. Based on findings reported herein, and those reviewed in ref. 56, 24 nt siRNAs alone are not sufficient to predict gene silencing by RdDM and other potential signals and features must be considered.

Higher changes in expression fold-change observed in TGS targets with near-gene transposons and 24 nt siRNAs

Pol IV-RdDM has been hypothesized to function at gene-proximal TEs to regulate expression of the gene. Based on our observation of enhanced expression changes at TGS targets with 24 nt siRNAs, we hypothesized that specific classes of TEs might be enriched in the regulatory regions of MOP1-regulated genes. We investigated the relationship between the presence of class II elements (Maize Transposable Element Consortium) within 2 kb flanking TSSs (+ or – 1 kb) and expression at targets of MOP1-mediated TGS. Of the genes differentially expressed in mop1-1 alone, we identified 153 upregulated and 177 downregulated genes with proximity to the conserved terminal inverted repeat sequence (TIRs) of class II transposons (Fig. 3D; Tables S6 and S7). Of these, 76 upregulated and 98 downregulated genes also had proximal 24 nt siRNAs (Fig. 3E). Downregulated genes displayed an expression level similar to that observed in the 413 mop1–1 downregulated genes. In comparison to downregulated genes, the average fold change expression of upregulated genes in these categories was higher, yet similar to upregulated genes without proximal RdDM features (Fig. 5A).

At putative TGS targets, we identified 29 upregulated and 18 repressed genes with proximal TIRs (Fig. 3D; Tables S8 and S9) and 16 upregulated and 11 downregulated with proximal TIRs and 24 nt siRNAs (Fig. 3E). We found that putative epigenetic targets with 24 nt siRNAs and proximal DNA transposons (TIRs) include genes with the highest average expression fold change in mop1 mutants (Fig. 5A). The average P value also decreased for genes in this subset (Fig. 5B and C), suggesting that these genes respond strongly and consistently to MOP1 activity. Notably, this group of 16 genes that we expect to be directly targeted by MOP1 was not represented in our prior analyses of nucleosome position in mop1-1 mutants, and may explain the small number of genes that we observed to be effected in that study.

In all categories examined, downregulated genes exhibited similar average expression changes irrespective of the presence of siRNAs and/or proximal TIRs (Fig. 5A). The observed increment in fold change expression with additive RdDM-associated features at putative epigenetic TGS targets demonstrates that multiple signals are utilized to confer MOP1 target specificity. To date, no consensus sequence exists for Pol IV and Pol V transcriptional targets (reviewed in reference 57) so it is difficult to identify direct targets of Pol IV-RdDM. Herein, we demonstrate that the effect of MOP1 on gene expression is enhanced at genes with proximal TIRs and homologous siRNAs. This trend may apply to only a subset of MOP1-mediated targets, but provides evidence that direct targets of MOP1 are best predicted by overlaying multiple indicating factors.

In summary, the chromatin accessibility profile established in mop1-1 immature ears is consistent with direct and secondary consequences of misregulation of gene expression, related to alterations in global chromatin structure in mop1-1. Consistently, we identified genes where expression increased and other genes where transcription was repressed in mop1-1 mutants. Although secondary consequences in expression are expected, a portion of upregulated genes are likely to be direct targets a MOP1-mediated transcriptional gene silencing pathway. We found that, on average, genes likely to be epigenetic TGS targets with homologous 24 nt siRNAs and proximal class II transposons included genes whose expression was more dramatically and consistently affected in mop1-1. We were unable to identify a similar subset of genes by limiting features in downregulated genes, suggesting that transcriptional repression in mop1-1 individuals may result from reduced ROS1 expression, or silencing by an alternative RdDM-TGS pathway.

Materials and Methods

Plant materials/genetic stocks

Genetic stocks were obtained from J. Dorweiler and are the result of the mop1-1 mutation (Dorweiler et al. 2000) introgressed in the B73 reference genome, and backcrossed for seven generations, at which stocks are considered near isogenic to the recurrent parent. The genotypes of the parental lines were as mop1-1 parent (K55: B’, Pl, r-g, mop1-1, and B73 recurrent parent: b, pl-sr, r-r, Mop1. 3–5 cm ear shoots were harvested and immediately frozen in liquid nitrogen and stored at –80 °C until use. Homozygous wild type (WT) and homozygous mop1 mutant (mop1-1) individuals were used. Plants were genotyped under the following PCR conditions: 94 °C for 5 min (1×); 95 °C for 30 s, 56 °C for 45 s, 72 °C for 45 s (30×); 72 °C for 10 min (1×); hold for 4 °C (as needed); using the following primers: KM384: (5′-TCTCCACCGC CCACTTGAT-3′); KM385: (5′-CCCAAGAGCT GTCTCGTATC CGT-3′); KM386: (5′-CTTCATCTCG AAGTAGCGCT TGTTGTCC-3′).

Nuclei isolation and cleavage with micrococcal nuclease (MNase)

Maize 3–5cm ear shoots were harvested and immediately flash frozen in liquid nitrogen. Flash frozen tissue was ground into a fine powder in liquid nitrogen. Nuclei were isolated from mop1–1 mutants and WT immature ears using a nuclei isolation buffer modified from Steinmüller and Apel, 1986. Frozen, ground tissue was added to Nuclei Isolation Fix Buffer (NIFB: 20 mM TRIS-HCl pH 7.8; 250 mM Sucrose; 5 mM MgCl2; 5 mM KCl; 40% Glycerol, 0.1% β-mercaptoethanol [βME]; 0.1 mM PMSF and 1% Formaldehyde). Nuclei were fixed for 10 min at room temperature (RT) after which formaldehyde was quenched by adding Glycine to 125 mM for 5 min at RT. Tissue was collected by centrifugation (1000 × g, 10 min, 4 °C) and resuspended in Nuclei Isolation Buffer with Triton (NITB: 20 mM TRIS-HCl pH 7.8; 250 mM Sucrose; 5 mM MgCl2; 5 mM KCl; 40% Glycerol, 0.1% BME; 0.1 mM PMSF; 1% Triton X-100). The suspension was filtered through 2 layers of nylon mesh and nuclei were collected by centrifugation (1000 × g, 10 min, 4 °C). Nuclei were gently resuspended in ice-cold Nuclei Resuspension Buffer/MNase Digestion Buffer (NRBT/MDB: 50 mM TRIS-HCl pH 7.5, 320 mM sucrose, 4 mM MgCl2, 1 mMCaCl2, 1% Triton, 0.1 mM PMSF), flash frozen, and stored until use. Nuclei (equivalent to 25 µg DNA) were subjected to digestion with limiting amounts of micrococcal nuclease (MNase: 0.1, 0.75, and 0.25 units) (Worthington Biochemical Corp) for 5mins at 37 °C and the reaction stopped with EGTA. Digested nuclei were decrosslinked at 65 °C overnight in the presence of 1% SDS and 0.2 mg/mL Proteinase K. DNA was extracted via phase extraction using phenol: chloroform and precipitated with NaOAc and ethanol and treated with RNaseA. DNA fragments were separated on a 1% agarose TAE gel. DNA was isolated from the respective gel fragments: MNase-resistant (high-molecular weight [HMW] DNA, undigested by MNase) and MNase-sensitive (low-molecular weight MNase digested DNA ladder, 150 bp to ~4000 bp representing DNA fragments readily accessible to digestion by limiting amounts of MNase). Each biological replicate consists of two ear shoots pooled together. The data represents the average of three biological replicates.

Microarray design, labeling and hybridization/MNase-chip

NimbleGen HD2 12-plex custom designed maize genome microarrays were used in these experiments. Each subarray contains ~134 300 unique, tiled and isothermal oligonucleotide probes that span the maize genome at a ~14 kb spacing designed to span all 10 maize chromosomes. The NimbleGen 12-plex slide ID is: 100428_ZMB73_JD_CGH_HX12. 500–1000 ng of isolated DNA was labeled using the NimbleGen Dual-Color Labeling Kit (NimbleGen) according to manufacturer’s protocol. Resistant/Inaccessible DNA (test) was labeled with Cy3 and Sensitive/Accessible DNA (reference) was labeled with Cy5. Twenty micrograms each (Cy3 and Cy5) labeled DNA were combined and used for hybridization. Arrays were hybridized for ~72 h at 42 °C after which slides were washed according to manufacturer’s protocol. Slides were scanned using a NimbleGen MS 200 Microarray Scanner and images were aligned using NimbleGen software.

Expression microarrays

Total RNA was isolated from maize ear shoots harvested and immediately flash frozen in liquid nitrogen and stored at –80 °C until use. Frozen tissue was ground into a fine powder in liquid nitrogen. Total RNA isolation was performed using TRI reagent according to manufacturer’s instructions (Molecular Research Center, Cat No. TR118). RNA was purified using Zymo RNA Clean and ConcentratorTM Kit (Zymo; Cat. No. R1017), according to manufacturer’s instructions. cDNA was synthesized using SuperScript TM III Reverse Transcriptase (Invitrogen Cat. No.18080-051). Briefly, 10 µg of RNA was reverse transcribed according to the manufacturer’s instructions using Oligo(dT)20 primers. cDNA synthesis was performed at 50 °C for 50 min. The reaction was terminated at 85 °C for 5 min. RNA was removed using 10 units RNaseH (Invitrogen) and 4 µg RNaseA (Qiagen) at 37 °C for 30 min. cDNA was purified using the QIAquick PCR Purification Kit (Qiagen). cDNA was labeled using the NimbleGen Dual-Color Labeling Kit (NimbleGen) according to manufacturer’s protocol. WT cDNA was labeled with Cy3 and mutant cDNA was labeled with Cy5. Six micrograms each labeled cDNA were combined and hybridized to a 385 K array. Arrays we hybridized overnight at 42 °C after which slides were washed according to manufacturer’s protocol. Slides were scanned using a NimbleGen MS 200 Microarray Scanner and images were aligned using NimbleGen software. NimbleGen gene expression microarrays designed by S. Kaeppler and R. Buell as described by Sekhon RS et al. The NimbleGen Array ID is: 090319_Zea_KR_ExpTil. Normalized using RMA (Robust Multi-chip Average). RMA normalized expression values were log2 transformed and fold change was analyzed using DNASTAR Arraystar® 4 software (http://www.dnastar.com/t-sub-products-genomics-arraystar.aspx). P values were calculated using the Student t test, one-tail.

Maize microarray expression annotation and gene ontology analysis

Maize gene annotations were obtained using NCBI standalone BLAST (http://blast.ncbi.nlm.nih.gov/Blast.cgi) by searching Arabidopsis TAIR10 (http://www.arabidopsis.org) and Rice Genome Annotation Project version 7.0 (http://rice.plantbiology.msu.edu/) databases. Top results were used as annotations for maize genes and e-value was set to 0.001. For upregulated and downregulated genes (2-fold cutoff), gene ontology (GO) analysis performed using BLAST2GO. Briefly, sequences were identified through a local Blast search of the NCBI nr database and imported into BLAST2GO. This program was used to extract GO terms associated to each of the obtained hits and returned an evaluated GO annotation for the query sequences. If there were no meaningful BLAST searches, the program returned N/A.

mop1 siRNA alignment

mop1 siRNA data from Nobuta et al. (2009) was downloaded from NCBI using the following GEO accession numbers: GSM306487 (WT) and GSM306488 (mop1-1). siRNA sequences were aligned to the B73 reference genome (v2) with Bowtie2 using default parameters. Small RNAs were grouped into four categories: ≤20 nt, 21–22 nt, 23–24 nt, and ≥25 nt using R. Herein, we used the 23–24 nt class of siRNAs for further analysis. BEDTools was used to calculate siRNA counts for the maize filtered gene set (FGS) and at transcription start sites (TSSs) of mop1-1 differentially expressed genes (DEGs). In-house Python and R scripts were applied to organize results and make a pipeline for analysis.

TGR1 gene expression analysis

Total RNA was isolated and pooled from immature ears of 20 Tgr1 non-mutant and 20 tgr1 mutant plants. RNA isolation was performed using TRI reagent according to manufacturer’s instructions (Molecular Research Center, Cat No. TR118). The libraries were prepared using NEBNext® Ultra™ Directional RNA Library Prep Kit for Illumina® (New England Biolabs). The quality of total RNA were measured by a Bioanalyzer and then subjected to RNA-Sequencing using an Illumina HiSeq2000 single-end sequencing instrument. Raw reads were trimmed by FASTX-Toolkit (http://hannonlab.cshl.edu/fastx_toolkit/) and aligned to B73 reference genome (AGPv2) using Tophat/Bowtie. Differential gene expression analysis was performed using the Cufflink v2.1.1 software package adopting default settings.