Direct and Indirect Transcriptional Effects of Abiotic Stress in Zea mays Plants Defective in RNA-Directed DNA Methylation

Plants respond to abiotic stress stimuli, such as water deprivation, through a hierarchical cascade that includes detection and signaling to mediate transcriptional and physiological changes. The phytohormone abscisic acid (ABA) is well-characterized for its regulatory role in these processes in response to specific environmental cues. ABA-mediated changes in gene expression have been demonstrated to be temporally-dependent, however, the genome-wide timing of these responses are not well-characterized in the agronomically important crop plant Zea mays (maize). ABA-mediated responses are synergistic with other regulatory mechanisms, including the plant-specific RNA-directed DNA methylation (RdDM) epigenetic pathway. Our prior work demonstrated that after relatively long-term ABA induction (8 h), maize plants homozygous for the mop1-1 mutation, defective in a component of the RdDM pathway, exhibit enhanced transcriptional sensitivity to the phytohormone. At this time-point, many hierarchically positioned transcription factors are differentially expressed resulting in primary (direct) and secondary (indirect) transcriptional outcomes. To identify more immediate and direct MOP1-dependent responses to ABA, we conducted a transcriptomic analysis using mop1-1 mutant and wild type plants treated with ABA for 1 h. One h of ABA treatment was sufficient to induce unique categories of differentially expressed genes (DEGs) in mop1-1. A comparative analysis between the two time-points revealed that distinct epigenetically-regulated changes in gene expression occur within the early stages of ABA induction, and that these changes are predicted to influence less immediate, indirect transcriptional responses. Homology with MOP1-dependent siRNAs and a gene regulatory network (GRN) were used to identify putative immediate and indirect targets, respectively. By manipulating two key regulatory networks in a temporal dependent manner, we identified genes and biological processes regulated by RdDM and ABA-mediated stress responses. Consistent with mis-regulation of gene expression, mop1-1 homozygous plants are compromised in their ability to recover from water deprivation. Collectively, these results indicate transcriptionally and physiologically relevant roles for MOP1-mediated regulation of gene expression of plant responses to environmental stress.

Plants respond to abiotic stress stimuli, such as water deprivation, through a hierarchical cascade that includes detection and signaling to mediate transcriptional and physiological changes. The phytohormone abscisic acid (ABA) is well-characterized for its regulatory role in these processes in response to specific environmental cues. ABA-mediated changes in gene expression have been demonstrated to be temporally-dependent, however, the genome-wide timing of these responses are not well-characterized in the agronomically important crop plant Zea mays (maize). ABA-mediated responses are synergistic with other regulatory mechanisms, including the plant-specific RNA-directed DNA methylation (RdDM) epigenetic pathway. Our prior work demonstrated that after relatively long-term ABA induction (8 h), maize plants homozygous for the mop1-1 mutation, defective in a component of the RdDM pathway, exhibit enhanced transcriptional sensitivity to the phytohormone. At this time-point, many hierarchically positioned transcription factors are differentially expressed resulting in primary (direct) and secondary (indirect) transcriptional outcomes. To identify more immediate and direct MOP1-dependent responses to ABA, we conducted a transcriptomic analysis using mop1-1 mutant and wild type plants treated with ABA for 1 h. One h of ABA treatment was sufficient to induce unique categories of differentially expressed genes (DEGs) in mop1-1. A comparative analysis between the two time-points revealed that distinct epigenetically-regulated changes in gene expression occur within the early stages of ABA induction, and that these changes are predicted to influence less immediate, indirect transcriptional responses. Homology with MOP1-dependent siRNAs and a gene regulatory network (GRN) were used to identify putative immediate and indirect targets, respectively. By manipulating two key regulatory networks in a temporal dependent manner, we identified genes and biological processes regulated by RdDM and ABA-mediated stress responses. Consistent with mis-regulation of gene expression, mop1-1 homozygous plants are compromised in their ability to recover from water deprivation. Collectively, these results indicate transcriptionally and physiologically relevant roles for MOP1-mediated regulation of gene expression of plant responses to environmental stress.

INTRODUCTION
The sessile nature of plants and their adaptation to terrestrial environments coincided with the evolution of whole plant and molecular responses to fluctuating environmental conditions (reviewed by Gupta et al., 2020). Extreme abiotic environments, including water scarcity, often lead to yield loss in agricultural crop plants across the globe (FAO, 2017). When osmotic stress is first detected, the initial and immediate whole plant response is often the closure of stomata, which allows the plant to conserve water within its tissues, while limiting the energy and resources expended in biological processes such as photosynthesis. More prolonged drought conditions result in responses that often limit plant growth, and are associated with developmental defects in reproductive organs, thus decreasing yield. Indeed, it has long been documented that Zea mays (maize) plants that experience drought stress exhibit reduced yield and the overall effects depend on the specific developmental stage at the time that stress is experienced (Claassen and Shaw, 1970). From a molecular perspective, recent studies demonstrate that whole plant responses are related to the disruption of gene regulatory networks and concomitant changes to stress-response transcriptional programs (Van den Broeck et al., 2017).
Changes in transcription at stress-responsive loci are often associated with genome-wide structural changes to chromatin that affect gene expression and can be detected as alterations in chromatin accessibility (Kim et al., 2015, reviewed by Chang et al., 2020. These changes are strongly influenced by the plant phytohormone abscisic acid (ABA), an important signaling molecule that is responsible for many processes throughout the life cycle of plants such as regulating several important stages of development, including seed germination, ABA synthesis, and signaling, and serves as a critical step in plant response to specific abiotic stress stimuli (reviewed by Ma et al., 2018, andTakahashi et al., 2020). In response to ABA, activation by phosphorylation of trans-acting factors initiates broad scale changes in gene expression, creating a hierarchical response that includes a combination of primary, secondary, and later stage cis and transacting responses at the molecular level (reviewed by Takahashi et al., 2018). It has also been observed that certain transcriptional changes in maize in response to and throughout recovery from drought stress is associated with differential enrichment for specific histone modifications (Forestan et al., 2020), and that differential DNA methylation is associated with water stress response in ABA-deficient maize mutants (Sallam and Moussa, 2021), further suggesting the overlapping regulation of ABA signaling and chromatin-mediated gene expression changes in plant stress responses. The coordinated effect of this multidimensional response can create whole-plant responses that originate from a molecular signal triggered by an environmental or developmental cue.
Activated trans-acting factors differentially regulate target chromosomal sequences, depending in part on the structure of chromatin at cis-regulatory elements (reviewed by Wang and Qiao, 2020). For example, evidence suggests that transcription factor binding is influenced by DNA (cytosine) methylation, although these mechanisms are not completely understood for a broad range of transcription factors (reviewed by Heberle and Bardet, 2019). Our recent investigations in maize seedlings indicates that genotypes defective in RNA-dependent DNA methylation (RdDM), a plant-specific epigenetic regulatory pathway, respond to exogenous ABA at the transcriptional level in a manner distinct from wild type plants (Vendramin et al., 2020). Genotype-specific changes in CHH (H = A, T or C) methylation were also observed at some loci transcriptionally responsive to ABA (Vendramin et al., 2020), which is consistent with prior observations for targets of RdDM (Gent et al., 2014). While this indicates that there is a relationship between transcriptional regulation by RdDM and ABA-mediated responses in maize, this association does not clearly distinguish between causality, dependence or coincidence. Interpretation is confounded by the fact that each regulatory network (ABA and RdDM) has primary and cascading indirect effects influencing to gene expression and methylation.
With regards to hormone signaling in response to environmental stress stimuli, time course experiments are a useful way to elucidate hierarchical relationships in complex regulatory networks, as the primary responses are generally expected to be triggered immediately following the stimulus, and the secondary and other downstream responses may require some time to occur. Time course analysis of ABA-regulatory networks in the model plant Arabidopsis thaliana suggest that ABA responsive changes in gene expression may be spread across an initial response period from 1 to 8 h after exposure to exogenously applied ABA (Song et al., 2016). To better understand the specific regulatory relationships of epigenetic gene regulation and abiotic stress responses in plants, changes in gene expression in maize plants that were either wild type or defective in RdDM were compared after 1-or 8-h exposure to exogenous ABA. Because these early responses can have longterm developmental and physiological effects on stressed plants, we also investigated the whole-plant responses of plants defective in RdDM to a severe drought simulation by withholding water for 14 days.

Plant Materials and Growth Conditions
Maize (Zea mays) plants with the mop1-1 mutation introgressed into the B73 inbred as previously described  were used for this analysis. Homozygous wildtype (Mop1 WT) and homozygous mutant (mop1-1) sibling progeny resulting from the self-pollination of an ear of a heterozygous plant were used. Seedlings were genotyped as previously described .
For abscisic acid treatment of maize seedlings: Seedlings were grown in greenhouse conditions (16 h light period, 25 • C, 50% humidity) in the Department of Biological Science at Florida State University (319 Stadium Drive) until they reached the V3 stage. At the V3 stage, maize seedlings were removed from the soil, roots were rinsed in water, dried, and then submerged in a 1 L beaker with 250 mL of liquid Murashige and Skoog (MS) media (Sigma Aldrich, M6899) with 50 µM ABA [ABA; (Sigma Aldrich, (+/-) Abscisic Acid, A1049)] or without ABA (MS) for 1 h (this study) or 8 h (Vendramin et al., 2020) in greenhouse conditions. After the incubation period, roots where removed and seedlings were immediately flash frozen in liquid nitrogen and stored at −80 • C until use.
For severe drought simulation on maize plants: Plants were grown in greenhouse climate-controlled conditions (25 • C, 50% humidity) at the Florida State University Mission Road Research Facility, Tallahassee, Florida, USA in January of 2017. B73 seeds were sown alongside as a control for drought response. Healthy B73, homozygous Mop1, and homozygous mop1-1 seedlings were then transplanted into 300 size pots and later into 2,000 size pots ∼35 days after sowing (DAS). Plants were randomly assigned to severe drought treatment (water withheld for 14 days) or normally watered groups. B73 (11 plants), Mop1 (11 plants), and mop1-1 (16 plants) individuals were in the normally watered control group and B73 (10 plants), 16 Mop1(16 plants), and mop1-1 (11 plants) individuals were in the drought treated group. The non-uniformity in sample number per category was due to premature death for a few individuals. Drought-treatment began once the individual plant reached the V6 stage to control for variation in development between samples/genotypes. After 14 days, plants entered the recovery phase by application of 7.5 L of water to the soil. After recovery, plants from the drought treated group were normally watered throughout the duration of the experiment. The tip of the V8 leaf (∼4 cm) was dissected from the maize plants, immediately frozen in liquid nitrogen, and stored at −80 • C until use.

Physiological Observations of Drought-Responsive Traits
All observations were made daily between the hours of 10:00 and 14:00. The growth stages of individual plants were determined using the leaf collar method (Nielsen, 2019). Plant growth was monitored beginning after seed germination (VE) and continued through the last collared leaf below the tassel (≤ V18). Daily observations were made to track the emergence of ears, tassels, silks, pollen shed, and the anthesis-silking interval (ASI) which is defined as the number of days between the first pollen shed and silk emergence (Bolanos and Edmeades, 1996). We determined the "effective tassel branch score" by inspecting each tassel and determining the ratio of functional tassel branches (branches with anthers) to total tassel branches (reported as a percentage). We determined the plant height at 90 DAS by measuring the length between the first node above the soil and the tip of the longest tassel branch. We determined the average internodal length by measuring the internodal distance of the three apical internodes above the V4 leaf node.

Total RNA Isolation, RNA Library Preparation and RNA-Sequencing
Total RNA was extracted as previously described (Vendramin et al., 2020). Briefly, frozen tissue was finely ground into powder in liquid nitrogen and homogenized before total RNA extraction was performed using TRI reagent R according to manufacturer's instructions (Molecular Research Center, 18080-051). RNA samples were DNase treated (RQ1 RNase-Free DNase, Promega, M6101) and purified using RNA clean and concentrator TM 25 (Zymo Research, R1018).
Three biological replicates were used for all RNA-seq experiments for each treatment and genotype, for a total of 12 samples per time point: 1 h (this study) and 8 h (Vendramin et al., 2020). The final sample concentration was quantified by Qubit. RNA library preparation (NEBNext R Ultra TM II kit, NEB, E7760) and Illumina paired-end 150 bp (PE 150) sequencing were performed by Novogene Corporation (Sacramento, California). The 1 h samples were sequenced using the Illumina NovaSeq 6000 platform, whereas the previous reported 8 h samples (Vendramin et al., 2020) were sequenced using the Illumina HiSeq 2,500 platform. More than 20 million reads were obtained per library.

Read Alignment, Batch Correction and Differential Gene Expression Calling
Bioinformatics analysis was performed by Linkage Analytics, LLC (Denver, CO). To ensure a consistent and re-producible computation environment, the workflow was containerized using Singularity (3.6.4) (Kurtzer and Sochat, 2017) and the data workflow steps were defined using Snakemake (5.30.1) (Koster and Rahmann, 2018) and read quality control was assessed using fastp (0.21.0) . Reads from the 1 h and 8 h sequencing batches were processed simultaneously. FASTQ adapters were trimmed by Cutadapt 1.8.1 (Martin, 2011) Reads were mapped to the B73 maize genome (AGP B73v4) (Jiao et al., 2017) by HISAT2 v2.2.1 (Pertea et al., 2016). Transcripts were assembled de novo to allow for inclusion of transcripts that are not included in the reference genome annotation and quantified using StringTie v2.1.4 (Pertea et al., 2016). Gene count matrices were generated from this data using the prepDE.py python script available in the StringTie website (http://ccb.jhu. edu/software/stringtie/index.shtml?t=manual). These matrices were used by the Bioconductor package edgeR 3.28.1  in R for differential expression analysis in order to identify upregulated and downregulated genes for the four different genotypes under two treatments. Lowabundance counts of < 0.58 cpm were removed using the DESeq2 filtering method (statquest.org/2017/05/16/statquestfiltering-genes-with-low-read-counts/); (Love and Huber, 2014) incorporated into the edgeR pipeline, and genes with an adjusted p-value of ≤ 0.05 and an absolute log 2 -fold change (FC) value of ≥ 0.95 were considered as differentially expressed for both upregulated and downregulated genes.

Gene Ontology Analysis and Hierarchical Clustering of Significantly Enriched GO Terms and DEGs
Singular Enrichment Analysis (SEA) was performed using the web-based tool agriGO v2.0 (Tian et al., 2017) with the B73 reference version 4 (AGOv4) gene annotations to determine enriched gene ontology terms (GO complete) associated with differentially expressed genes.
Fisher's statistical test, Hochberg (FDR) multi-test adjustment method with a significance level of < 0.05 and minimum number  of mapping entries of 10 genes per GO-term. The GO term enrichment was generated by hierarchically clustering the log10 of the total GO term percentage of a set of genes that were upregulated or downregulated in wildtype or mutant in response to ABA.

RESULTS
Early ABA Treatment Is Sufficient to Induce Unique Categories of Differently Expressed Genes in mop1-1 Mutants To identify genes that are immediately responsive to epigenetic regulation under abiotic stress conditions, RNA from maize seedlings exposed to 1h of abscisic acid (ABA) and nutrient solution without ABA (MS) in mop1 wildtype (WT) and mutant (mop1-1) genotypes was subjected to RNA-sequencing (RNA-seq) and transcriptome analysis as previously described (Vendramin et al., 2020). An average of ∼33 million 150 bp paired end raw reads were obtained per sample ( Table 1) and mapped to the B73 maize genome (AGP B73v4) (Jiao et al., 2017). Significant differentially expressed genes (DEGs) between mop1 genotypes and 1h ABA treatments were categorized into four pairwise comparisons designated "analysis groups" (1h Groups A-D; Table 2) as genes with a two-fold expression change (log 2 FC ≥ 0.95, FDR < 0.05) and identified by making direct comparisons between genotypes and treatments ( Table 2). The DEGs in the four analysis groups were further sub-divided based on gene expression patterns (up-or down-regulated; e.g., 1h A-up and 1 h A-down) ( Table 2; File 1) and subjected to further analysis. The total number of significant DEGs with two-fold change (2FC) in expression identified in the four analysis groups after 1 h of ABA-induction (1 h Groups A-D) included 1,856 genes, where, 737 (40%) of these genes were found to be common to more than one group, resulting in 1,119 (60%) DEGs unique to an individual analysis group (Table 2; Figure 1). After 1 h of ABA induction, only ∼7% of the total 2FC DEGs were differentially expressed in WT and mop1-1 genotypes relative to their own control (1 h Groups A and B). These transcriptional responses are genotype-specific as there was also almost no overlap between the DEGs identified in each of these analysis groups (1 h Groups A and B) ( Figure 1A). The majority of DEGs (63%) were identified in comparisons that included both genotype and treatment (Table 2 Group C; mop1-1 ABA/WT ABA). A comparison between Group C DEGs with control plants of the same genotypes (mop1-1 MS/WT MS from Group D), revealed 417 and 451 up-and down-regulated genes, respectively, that are uniquely responsive to ABA treatment and loss of MOP1 activity ( Figure 1B). These 868 genes were subjected to a more in-depth analysis to identify primary and indirect MOP1 specific targets.
Gene ontology (GO) analysis was used to predict the biological processes of all annotated genes in each of the four analysis groups (1h Groups A-D; FDR < 0.05). As expected, the GO term for response to stimulus (GO:0050896) was highly enriched in all 1 h analysis groups, except for the comparison constituting a genotype control of mutant and wild type plants treated with MS (Group D; Figure 1C). The diversity of enriched DEGs was enhanced in mop1-1 mutants subjected to ABA (1 h Groups A and C) relative to WT plants (Group B) or the genotype control (Group D) (Figure 1C). These mop1-1 ABA unique categories include biological processes associated with cell growth and size ( Figure 1C).
In Mop1-1 Mutants, the Most Distinct Changes in Gene Expression Occur Within the Early Stages of ABA Induction To identify genes that respond to ABA and MOP1 in a temporal manner, the mapped reads from RNA-seq after 1 h (this study) and 8 h (Vendramin et al., 2020) of ABA induction were simultaneously, bioinformatically processed and mapped to the B73 reference genome (AGP B73v4) (Jiao et al., 2017) and used in subsequent analysis (Tables 1-3). Due to the differences in sequencing depth as a result of use of different Illumina sequencing platforms (HiSeq vs. NovaSeq) between the two timepoints, we normalized the read quality score threshold used in HISAT2 ("-score-min") between platforms. Based on consistency between replicates as well as differences in distributions of mapping qualities between sequencing platforms, the HISAT2 "-score-min" parameter was chosen to normalize the number of uniquely mapped reads across datasets (Table 1;  Supplementary Figure 1).
Predictably, the overall number of DEGs increases with increasing time. Eight h of ABA treatment resulted in more genes exhibiting differential expression, but most of the genes were detected in multiple analysis groups, resulting in a lower percentage of DEGs being unique to one analysis group at 8 h (27%) compared to 60% unique DEGs observed after 1 h of ABA-induction (Tables 2, 3). Consistently, while there was almost no overlap between wildtype and mutants DEGs in 1 h Groups A and B (Figure 1), there was more overlap between Groups A and B after 8 h of ABA treatment (Supplementary Figure 2; Supplementary Table 2). For analysis groups C, there were more significant DEGs at 1 h of ABA treatment, compared to the same comparison after 8 h (Tables 2,  3). This suggests that the most distinct changes in gene expression between these genotypes occurs within the early stages of ABA induction.
Genes from the 8 h and 1 h samples were directly compared with each other (8 h/1 h) and categorized into four different pairwise comparisons, designated Groups E-H ( Table 4). The total number of significant 2FC DEGs identified in these four 8 h/1 h analysis groups was 34,147 genes, representative of the magnitude of changes in gene expression that occur over time. However, 29,610 (87%) genes were found to be common to more than one analysis group, where only 4,537 (13%) were found to be unique to an individual group (Table 4; Supplementary Figure 3). This observation is consistent with the similarities in the enriched GO terms per group (Figure 2), with the least diverse biological processes observed in the wildtype control group (Group H). GO terms associated with biological regulation (GO:0065007), regulation of biological process (GO:0050789) and response to stimulus (GO:0050896) were commonly highly represented terms across all groups (Figure 2).

MOP1-Dependent siRNAs and Gene Regulatory Networks (GRNs) Predict Immediate and Indirect Responses to Abiotic Stress
To distinguish between primary and indirect targets of epigenetic regulation under abiotic stress conditions, the 868 genes (451  up-and 417 down-regulated) identified to be uniquely associated with ABA treatment and loss of MOP1 activity (Group C; Figure 1B) were further analyzed for a specific connection with MOP1-mediated RdDM. MOP1 is required for the production of the majority of siRNAs at loci undergoing RdDM (Gent et al., 2014), therefore these genes were compared with a list of genes having promoter homology with MOP1-dependent siRNAs (Vendramin et al., 2020). This comparison identified 97 up-and 76 down-regulated genes from 1 h Group C that are predicted to be direct MOP1-regulatory targets based on homology with siRNAs (Figure 3; Supplementary Table 3), suggesting that MOP1-mediated RdDM is involved in early responses to ABA at these specific genes. It is plausible that these 173 genes are primary targets of MOP1 that in turn influence downstream gene expression in response to ABA. Because these genes are differentially responsive in mop1-1 plants very early after ABA treatment, these genes were designated as MOP1-dependent immediate responsive genes (MIMs). Gene ontology analysis of these genes revealed that there were more significant (FDR < 0.05) enriched GO terms associated with the 97 up-regulated genes compared with the 76 down-regulated genes ( Table 5). This suggests that in   response to ABA, MOP1-dependent activity and siRNAs are directly associated with regulation of specific biological processes, whereas the siRNAs associated with downregulated genes (MOP1-independent) may be indirect, not RdDM targets and/or have less specific biological roles in relation to ABA responses. To understand how MIMs potentially influence downstream transcriptional responses to 1 h ABA treatment in maize, a gene regulatory network (GRN) (Huang et al., 2018) was used to predict targets of these 97 and 76 up-and down-regulated genes, respectively. Twenty one of the 97 (∼22%) upregulated genes and 5 of the 76 (∼7%) downregulated genes had predicted regulatory targets based on the GRN, and the majority of these 26 genes are transcription factors implicated in drought, ABA and stress responses based on homology and phenotypic characterization in other studies ( Table 6). The predicted GRN targets of the 21 Group C 1 h upregulated MIMs included a total of 16,748 genes and the predicted GRN targets of the 5 Group C 1 h downregulated MIMs included a total of 4,221 genes across all tissues types and datasets in the GRN ( Table 6; File 3). Some of these genes (∼14%) were duplicated in the two lists of targets predicted by the GRNs and overall there were 18,014 unique target genes. These targets predicted by the GRN could be considered indirect (secondary or more downstream) targets of MOP1 responsive factors, because they are predicted to be regulated by genes with evidence of direct MOP1-mediated regulation. This group of genes were collectively designated as MOP1-dependent indirectly responsive genes (MINs).
The 18,014 MINs ( Table 6; File 3) were compared with genes in analysis groups E and G (Table 4), representing genes differentially expressed in a temporal manner (8 h/1 h) after treatment with ABA in mutant and wild type, respectively ( Figure 4A). This analysis revealed that in the mutant genotype, 54% of upregulated and 42% of downregulated genes with expression changes from 1 h to 8 h in the presence of ABA were identified as MINs (Figure 4A). A similar comparison in wild type identified that 52% of upregulated and 42% of downregulated were identified as MINs (Figure 4B). This initial observation suggests that over time, there is a similar magnitude of indirect effects in MOP1-regulated targets between mutant and wild type plants subjected to abiotic stress stimuli. However, further analysis revealed qualitative and functional differences in DEGs, as the identity of genes did not completely overlap ( Figure 4C). Specifically, there were 689 (44% from wildtype (G) and 51% from mutant) putative indirect DEGs between 1 h and 8 h of ABA treatment common to both genotypes (Figure 4C;  Supplementary Table 4). A gene ontology analysis was used to identify genes that fall into the three categories ( Figure 4D) and GO terms associated with response to stimulus were conserved in these three categories. It appears that mop1-1 mutants are not expressing some of the developmental genes required for MOP1-mediated responses to ABA. This comparison of DEGs between mutant and wild type may be indicative of the role of MOP1 in maize development in response to some abiotic stress stimuli. This association is also supported by recent work in other labs, indicating a role for RdDM or other chromatin-mediated regulatory events in stress response in maize (Forestan et al., 2016, Forestan et al., 2017, Forestan et al., 2020.

MOP1 Is Required for Recovery From Drought Stress
To determine the role of MOP1 in drought stress response and recovery at the whole-plant level, we characterized the vegetative and reproductive developmental consequences of a severe drought treatment (14-days without watering) on Mop1 WT, mop1-1, and B73 plants. Initially, all plants were watered normally and showed no significant differences in growth rate until reaching the V6 stage (Nielsen, 2019) when drought treatment was applied to randomly selected individuals from each group (Figure 5). We controlled for growth rate differences between individuals by beginning the drought treatment when an individual reached the V6 stage (auricle exposed). The growth rate was significantly delayed among all drought-treated plants compared to normally-watered controls, with mop1-1 droughttreated plants taking the longest time to reach the V7 stage ( Figure 5). After 14 days of drought treatment, water (7.5 L) was given to each plant and plants were normally watered throughout the duration of the experiment. While B73 and Mop1 droughttreated plants recovered rapidly and approached the growth rate curve of normally-watered controls, mop1-1 drought-treated plants significantly lagged behind, with several plants failing to reach reproductive competency ( Figure 5A). Normally-watered mop1-1 plants showed no significant differences in growth rates compared to normally-watered B73 and Mop1 suggesting that MOP1 is required to recover from drought stress.
To determine the effects loss of MOP1 during drought had on reproductive development, we made observations for plant height at maturity, internodal length, ear emergence, number of ears, effective tassel branches, and the anthesissilking interval (ASI) (Figure 5B; Supplementary Figure 4). Because stunted plant height is indicative of severe stress during vegetative development, we measured the average internodal length and the heights of plants at 90 days after sowing (DAS). Drought-treated plants were stunted and had reduced internodal lengths compared to normally-watered  controls across genotypes ( Figure 5B; Supplementary Figure 4). Reproductive development can also be affected by drought-stress and the magnitude of the effect is in some cases dependent on the stage of development in which the plant endures the stress. Because the drought-treatment in our study begins at the V6 stage, which is prior to the transition to reproductive development, we were able to determine the effects of vegetative stress on reproductive traits. To determine how drought affects tassel development, we characterized the effective tassel branches for each individual by measuring the ratio of tassel branches with functional anthers (i.e., anthers shedding pollen) to total tassel branches and found a drought-dependent decrease in effective tassel branches across genotypes, however, these differences were not significant (Supplementary Figure 4). The number of days until ear emergence and the number of ears per plant were also measured. It was found that drought treatment led to a significant delay in ear emergence and that mop1-1 drought treated plants, but not B73 or Mop1, displayed a significant reduction in the number of ears per plant ( Figure 5B;  Supplementary Figure 4). In addition, mop1-1 drought-treated plants displayed a significantly larger anthesis-silking interval compared to B73 and Mop1, suggesting that impaired recovery from drought stress in plants defective in MOP1 function has an effect on reproductive development and competency ( Figure 5B).

DISCUSSION
Understanding the molecular mechanisms that contribute to plant responses to changing environments is essential to ensure that we can develop climate resistant plants that meet the increasing global demands on crop yield. RdDM and ABAsignaling are two critical gene-regulatory pathways that each influence how plants respond to environmental cues at specific developmental stages. The extent of the synergy between these two regulatory systems is largely uncharacterized in agronomically important crop plants, such as maize. To address this gap in knowledge, we recently conducted a transcriptomic analysis which demonstrated that loss of RdDM activity renders maize seedlings more susceptible to transcriptional changes as a result of ABA treatment, and that many genes were responsive to disruption of both regulatory networks after 8 h of phytohormone treatment (Vendramin et al., 2020). The differential response of the RdDM-deficient mutant to treatment with ABA and to water deprivation suggest that stressful growing conditions or exogenous of application of growth hormones like ABA might be sufficient stimuli to alter the epigenome of maize, and could be useful in crop epi-breeding platforms, which may enhance modern breeding efforts (Dalakouras and Vlachostergios, 2021). While this study identified and established synergy between these two networks in maize, interpretation of the results was confounded by the hierarchical nature of cascading transcriptional outcomes for both regulatory pathways, each dependent on varied cis and trans-regulatory elements.
Using an approach based on a temporal response to phytohormone treatment, we have identified immediate and direct MOP1-dependent transcriptional responses to ABA (MIMs) that are predicted to function upstream of genes responsive to longer periods of exposure to abiotic stress stimuli (MINs). These relatively few MIM genes, identified as unique 1 h Group C genes having homology with MOP1-dependent siRNAs, appear to be specific in their biological function. Using a GRN, we were able to establish a hierarchial relationship between predicted MIMs and MINs, where the MIMs identified in this study are predicted to regulate ∼50% of genes differentially expressed after longer exposure to abiotic stress (8 h), suggesting a substantial impact on transcriptional responses by MIMs. The lack of multiple enriched GO terms associated with the 76 down-regulated MIMs is indicative of either a lack of biological specificity of these genes or a reflection of the complexity of regulation of genes in this category that may also be targets of an active demethylaion mechanisms by DNA glycosylases. For this study, only a subsets of possible regulatory features associated with RdDM activity were used to identify MIMs, and yet the predicted regulatory impact of the identified MIMs account for almost half of the MINs, suggesting that this may in fact be an underestimation of the contribution of MOP1 in establishing responsive transcriptional profiles. Additional analysis to include other RdDM regulatory features, such as proximity to specific categories of transposable elements (TEs) Vendramin et al., 2020) and contexts of cytosine methylation that establish boundaries between the TEs and adjacent protein-coding genes in maize (Gent et al., 2014;Li et al., 2015) might identify additional specific ABAinduced MIM genes. It is likely that an extensive genomewide analysis will need to be pursued to elucidate specific examples of direct correlation between DNA methylation, chromatin marks and differential expression in these conditions, because prior work has demonstrated that these coordinated responses are hierarchical and inter-related, and often do not involve simple relationships between differential expression and FIGURE 5 | Effects of drought treatment on growth rate and reproductive development. (A) The gradient on the x-axis represents the approximate vegetative to reproductive developmental transitions, from the growth stage where two vegetative leaves have emerged (V2) through the 8 vegetative leaf stage (V8), and the reproductive stages where tassels and ears are present (VT) through later stages of flowering/seed set (R1 to R2). Individual points indicate mean days after sowing (DAS) when the developmental stage was reached for the indicated genotypes and treatment groups, and error bars show standard deviation within groups. Data is shown for B73 inbred lines under different watering conditions, and wild type (Mop1) and homozygous mop1-1 individuals segregating within a family. Significant differences (p < 0.05) between drought-treated mop1-1 and Mop1 plants are indicated for selected points with an asterisk (*) or NS for no significance. (B) Physiological observations from drought stress experiment showing average internodal length, number of ears per plant, and anthesis-silking interval (ASI). Significant differences (p < 0.05) between groups are indicated by an asterisk (*) or no significance (NS).
hallmarks of RdDM Vendramin et al., 2020). Thus, a locus-specific approach was not attempted in this study.
There is already compelling evidence indicating that RdDM activity in maize has consequential effects on plant growth and development, affecting the male and female inflorescences (Dorweiler et al., 2000;Hultquist and Dorweiler, 2008) and ultimately seed yield (Barber et al., 2012). The study described herein, reveals that, consistent with mop1-1 plants misexpressing genes involved in development (Vendramin et al., 2020), plants defective for RdDM are compromised in their growth rate recovery after water stress. This observation links the differences in transcriptional responses of maize mop1-1 plants to differing abilities to recover from abiotic environmental influences, and highlights the physiological relevance of the gene expression phenotypes of RdDMdeficient plants.
Collectively, this data suggests that MOP1 activity is required for preparedness to respond, early response and later response to ABA signaling at the level of gene expression, and may indicate that MOP1, a component of RdDM in maize, functions in plant response to stressful growth conditions. Future work will include molecular characterization of the MIMs to identify the architecture of upstream cis-regulatory elements of these genes.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. RNA-sequencing data is available at the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (Edgar et al., 2002) through GEO Series accession number GSE179629.

AUTHOR CONTRIBUTIONS
TM, SV, JL, and KM designed the research. SV performed the seedling ABA-induction experiments and RNA-seq library preparation. JL and PL performed the drought stress experiments. TM and KM performed analysis and interpretation of seedling ABA RNA-seq data. KL assisted with bioinformatic analysis and interpretation. TM and KM wrote the manuscript. All authors edited and/or reviewed the original manuscript, and contributed to the article and approved the submitted version.

FUNDING
This work was funded by the National Science Foundation funds to KM (CMB-035919) and start-up funds from the University of Washington Bothell School of STEM to TM.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2021. 694289/full#supplementary-material Supplementary Figure 1 | Threshold normalization for uniquely mapping reads. The mapped reads from RNA-seq after 1 h (this study; Illumina NovaSeq platform) and 8 h (Vendramin et al., 2020; Illumina HiSeq platform) of ABA induction were simultaneously, bioinformatically processed and mapped to the B73 reference genome (AGP B73v4) (Jiao et al., 2017). Read quality score thresholds were normalized using HISAT2 ("-score-min") between sequencing platforms. Based on consistency between replicates as well as differences in distributions of mapping qualities between sequencing platforms, the HISAT2 "-score-min" parameter was chosen to normalize the number of uniquely mapped reads across datasets. HISAT2 filters reads based on a threshold were defined by the slope a linear function between mapping quality score and read length. HISAT2 slope filter threshold for 1 h = −0.2 and for 8 h = −0.6. Genotype-treatment samples include wild type MS (wm), wild type ABA (wa), mutant MS (mm), mutant ABA (ma). Data is shown for B73 inbred lines under different watering conditions, and wild type (Mop1) and homozygous mop1-1 individuals segregating within a family. Significant differences (p < 0.05) between groups are indicated by an asterisk ( * ) or no significance (NS).