Introduction

Rice is a major staple food crop that contributes to the diets of half of the world’s inhabitants, especially in Asia. However, multiple environmental stresses, including drought, flooding, high temperature, salinity, among others, affect crop productivity and sustainability of agriculture. Drought massively affects the growth and yield of rice, which is a semiaquatic plant and considered as one of the most drought-susceptible plants. Almost all the developmental stages of rice can be affected by drought, especially the reproductive stage, resulting in severe yield penalties. Therefore, it is urgently needed to breed drought-resistant/tolerant rice for improving rice yield under drought to meet with the increasing population and the drought challenge that occurs. Understanding the mechanisms of how rice responds and adapts to drought stress will provide us grounds for optimism.

Rice employs complicated processes and multiple strategies in response to drought, involves various physiological, biochemical, and genetic responses. Many traits under drought have been reported, such as modification of cell wall plasticity, root extension, leaf area, and surface properties (Ganie and Ahammed 2021; Kim et al. 2020a, b), etc. Multiple ways at genetic levels, transcriptomic levels, and proteomic levels were extensively used to reveal the drought stress genes and pathways. Currently, several transcriptional factors (TFs) are known to play crucial roles in regulation of a series of downstream genes (e.g., cell expansion and cell wall-related genes, lignin biosynthetic gens, cellulose synthase) to cope with drought, such as AP2/ERF family TF genes (Jung et al. 2017; Lee et al. 2016), the NAC group of TFs(NAM, ATAF, and CUC) (Hu et al. 2006; Nakashima et al. 2007; Redillas et al. 2012; Zheng et al. 2009), WRKYs (Sahebi et al. 2018; Shen et al. 2012). MicroRNAs (miRNAs), which regulate diverse developmental and stress-related processes, are regarded to modulate these traits as well (Zhang et al. 2018b; Nadarajah and Kumar 2019). Various components of ABA signaling, plant hormones crosstalk (Todaka et al. 2015; Riemann et al. 2015; Santosh Kumar et al. 2020), and metabolic networks (Fàbregas and Fernie 2019) are well described in response to drought stress in different species. Other than these, DNA methylation, which is one of the main types of epigenetic modifications that commonly presented in eukaryotic genomes and significantly regulate gene expression, also play a critical role in response to plant drought stresses (Li et al. 2020).

DNA methylation occurs at sequence contexts of CG, CHG, and CHH (H=A, C, or T). CG sites methylation are usually mediated by methyltranserase 1 (MET1), while CHG sites are mediated by chromomethylase 3 (CMT3), and CHH sites are governed by domains rearranged methyltransferase 1/2 (DRM1/2) and chromomethylase 2 (CMT2) (Kankel et al. 2003; Law and Jacobsen 2010; Lindroth et al. 2001). Using Methylation sensitive amplified polymorphism technique, hypermethylation was found dominant in drought-susceptible rice genotypes IR20 and CO43, while hypomethylation was found to be presented in drought-tolerant rice genotypes PL and PMK3 when under drought condition (Gayacharan 2013). 68 rice accessions of osmotic stress-tolerant and -susceptible groups were further checked for their DNA methylation levels, with similar conclusion that drought-tolerant accessions possessed lower DNA methylation and more de-methylation events under stressed condition (Xia et al. 2017). A latter study using methylated DNA immunoprecipitation sequencing and Affymetrix GeneChip array also reported that differentially methylated regions (DMRs) in drought-tolerant introgression line DK151 was far more less than DMRs in drought-sensitive line IR64, which suggested drought-tolerant rice plants have a more stable methylome (Wang et al. 2016). DNA methylation was positively correlated with spikelet sterility and negatively correlated with yield, and has various effects on gene expression (Gayacharan 2013; Wang et al. 2016). A recent study also reported that methylation in CG and CHH contexts within gene body and distal promoter regions were positively correlated with the expression of abiotic stress response genes (Rajkumar et al. 2020). Hypomethylation of Transposable elements (TEs) in plants were associated with transcriptional/transpositional activities (Hu et al. 2014; Miura et al. 2001). Methylation of TEs nearby genes affect their expression levels. TEs were found methylated mostly (> 80%) in CG sites, and 20–90% in CHG sites, 2–30% in CHH sites (Niederhuth et al. 2016). CG, CHG, and CHH methylation in A. thaliana TEs averages at about 80%, 40%, and 15%, respectively. CG, CHG, and CHH methylation in P. patens TEs averages around 80%, 80%, and 30%, respectively (Domb et al. 2020). DRMs targeted to TEs via RNA-directed DNA methylation pathway, and methylation under all three contexts are important for TE silencing. These evidence supported that DNA methylation is an important epigenetic regulatory mechanism for rice to adapt drought stress.

However, throughout its whole life, plants can constantly suffer from environmental fluctuations such as drought, which might occur in a short time (transient) or over long periods. Plants submitted to repeated stress treatments responded differently than to a single stress event, which called “stress memory”. Stress memory can trigger physiological and molecular processes and help plants to enhance their resistance to stress. Many plants were found to exhibit memory behavior, such as switchgrass (Zhang et al. 2018a, b), Zea mays (Virlouvet et al. 2018), and rice (Li et al. 2019), Brazilian savanna (Alves et al. 2020), Glycine max L. (Kim et al. 2020a, b), Pinus pinaster (Fernández de Simón et al. 2020). Exposure to drought resulted in several changes in plants, and these changes in plants submitted to recurrent drought stress differed when compared in a single drought event. This helps in acclimation to drought stress for plant survival, involving several metabolic pathways, gene expression changes, and physiological and phytohormones changes. Proline has been considered to be a critical component of stress memory, including drought (Li et al. 2019). The free proline levels increased by the first drought stress and remained stable throughout the subsequent drought treatment, in corresponding to the gene expression of proline biosynthetic enzyme ∆1-pyrroline-5-carboxylate synthetase 1 (P5CS1) in rice (Li et al. 2019). Besides, DNA methylation was also believed to participate in rice short-term drought memory, although the mechanism was not clear. Multiple ways were employed to study the plant epigenetic stress memory (Godwin Farrona 2020). By methylation-sensitive amplification polymorphism (MSAP) method, both drought-sensitive and -resistant varieties showed a cumulative effect on DNA methylation pattern from the original generation to the sixth generation, which suggested that some DNA methylation status induced by drought stress could inherit and transmit to the next generation (Zheng et al. 2013). Wheat seeds from terminal drought stressed plants showed better growth than the progenies of well-watered crop, and seed priming improved drought tolerance (Tabassum et al. 2018). Meanwhile, rice also was found exhibited short-term drought stress memory in response to cycles of mild drought and re-watering, and several genes involved in rice drought memory response were identified through a whole-transcriptome strand-specific RNA sequencing (Li et al. 2019). A linkage was found between drought memory transcripts and DNA methylation (Li et al. 2019), which suggested that DNA methylation participated in plant drought memory. However, how DNA methylation changes respond to drought and the subsequently drought stress, and how it participated in drought memory is still not clearly identified.

Previously, we established a system that could induce rice drought memory via cycles of mild drought and re-watering treatment and reported 6885 transcripts that involved in the drought memory response. Moreover, we also found some DNA methylation changed located in regions of drought memory genes we found. In this study, we further studied the DNA methylation patterns after cycles of drought stress and re-watering process, and uncovered the contribution of DNA methylation to plant short-term drought memory formation.

Methods

Plant Materials and Growth Conditions

Seeds of rice (Oryza sativa L. ssp. Japonica cv. Zhonghua 11) were germinated under sterilized condition and placed on four layers of moistened filter paper as described previously (Li et al. 2019). Seedlings were grown in darkness for four days at 30 °C and then transferred to a hydroponic growth system using 1/4 modified Hoagland solution at 28 °C for 4 weeks. The light/dark cycle was 12 h/12 h and light intensity was 180 μmol m−2 s−1.

Drought and Re-watering Treatment

4 weeks old seedlings (non-treatment control, R0) were air-dried for 80 min (first drought stress treatment, S1) and fully re-watered for 22 h (first re-watering samples, R1). To monitor the memory of drought stress, seedlings were sequentially dried and re-watered twice. Leaf samples from R0, S1, R3, and S4 treated seedlings were harvested for analysis (Supplemental Fig. 1). Samples were immediately frozen in liquid nitrogen and stored at − 80 °C before processing. All the treatments were performed at 28 °C. Relative water content (RWC) was used for rapid estimation of drought stress: RWC = (FW − DW)/(RW − DW) × 100% (Ding et al. 2014), where FW corresponds to the fresh weight, DW and RW to drought-treated weight and re-watered weight, respectively.

DNA Library Construction and Whole-Genome Bisulfite Sequencing

Total genomic DNA of rice seedlings were extracted using the CTAB method. The amount and quality of DNA was examined with Nanodrop-2000 (Thermo Fisher Scientific) and agarose gel electrophoresis. DNA library construction and sequencing was performed by Biomarker technologies. In brief, DNA was fragmented via sonication, and then the fragments were ligated with adapters and purified, then bisulfite converted using EpiTech Bisulfite kit (Qiagen, United states). A pair-end sequencing was performed on an illumina Hiseq ™ 2500 Platform.

Bioinformatics Analysis

Raw reads were quality controlled and trimmed and reads were then aligned to the rice reference genome on phytozome (https://phytozome.jgi.doe.gov/pz/portal.html#!info?alias=Org_Osativa) using BSMAP (version 2.74). The individual methylation ratio was extracted using python from the output of BSMAP mapping results. The methylation rate of mCG, mCHG, and mCHH in the whole genome was determined as Methylated cytosine/total cytosine. DMRs were calculated by sliding windows with 200 bps slide window width and 50 bps step size. Cytosine sites covered by at least six reads in each bin were considered. Fisher test and False discovery rate (FDR) were used for P value (P value < 0.05) calculation and correction. A minimum difference of 40% methylation level, and also the methylation rate was greater than 10% in at least one time point was identified as DMRs. For DMRs that only found in the subsequent time points but not the initial one, the threshold of methylation difference for DMRs was set at 30% for mCG and mCHG, and 20% for mCHH. Circos software was used to construct Circos Plots (Krzywinski et al. 2009). Short Time-series Expression Miner program (STEM, version 1.3.11) (Ernst and Bar-Joseph 2006) was used to visualize and compare the methylation rate among each sample to obtain drought memory profiles.

Gene ontology (GO) enrichment analysis was performed using WEGO website (http://wego.genomics.org.cn). Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis was performed by KEGG orthology based annotation system (KOBAS, version 3.0).

Results

Differential DNA Methylation in Rice Genome Under Drought Stresses

According to our previous research, rice seedlings display drought stress memory after cycles of appropriate water deficient treatment. To analyze and visualize the global DNA cytosine methylation pattern after drought stress and the re-watering process, BS-seq was used in this study. Samples were taken from four-week-old rice leaves before (R0) and after the first 2-h drought stress (S1), and also at the time points after the third 22-h re-watering process (R3) and after the fourth drought stress (S4) (Supplemental Fig. 1).

After performing quality control and trimming process, 63–78 millions of high-quality reads were obtained and 88–94% of them were matched to the rice reference genome with genome coverage of 81.3–91.0% (Supplemental Table 1). Moreover, all the bisulfite conversion rates were higher than 99.5%, which indicating the accuracy and integrity of whole-genome bisulfite sequencing of this study.

The genome-wide cytosine sites (Supplemental Table 2) and their methylation levels at different time points (R0, S1, R3, S4) under all three sequence contexts (CG, CHG, and CHH; H represent A, T, or C) were calculated. The results showed that mCG was the most abundant in three contexts in all conditions we tested, while cytosine methylation level under CHG and CHH showed comparable levels, which were 68.56% and 63.12% of the numbers under CG at R0, respectively (Fig. 1A). The relative DNA methylation rate of mCG is calculated as mCG/(mCG + mCHG + mCHH) × 100%, the methylation rate of mCG was 43.18% at R0 without any treatment (Fig. 1B), while the genomic C sites rate under CG sequence context was about 18.96% in the whole genome (Supplemental Table 2). On the other hand, mCHG and mCHH methylation rates were 29.60% and 27.21% (Fig. 1B), while C site rates under CHG and CHH in the entire genome were 16.84% and 61.2% (Supplemental Table 2), respectively. These results suggested that large numbers of CG sites were cytosine methylated, while only a small part of CHH sites was methylated even the context of CHH is most abundant in rice genome.

Fig. 1
figure 1

DNA methylation levels during drought training treatment. A Total number of methylation events under contexts of CG, CHG, and CGG. B Proportion of methylation under each context. C Circos plots of DNA methylation under contexts of CG, CHG, and CHH in chromosomes. Rings of J01, J02, J03, and J04 show genomic positions of DNA methylation at R0, S1, R3, and S4 time points. Rings of Gene and TE represent the density of genes and TEs on chromosomes

Full size image

The genome-wide cytosine methylation levels in each condition (R0, S1, R3, S4) under CG, CHG, and CHH were further compared. DNA methylation level increased after S1 (first drought stress) compare to R0 (without any treatment), among which cytosine methylation levels under CG showed 10.86% of increasing, while under CHG and CHH only 7.52% and 2.01% increasing was found compared to their R0, respectively. After subsequently drought stresses at R3, the methylation was reduced to an even lower level compared to R0 (R3 vs R0), which decreased 7.34%, 5.45%, and 0.60% under CG, CHG, and CHH, respectively. A second induction after subsequent drought stress (S4) was observed under all three sequences, the increase of methylation level (S4 vs R3) was 22.80%, 16.87%, and 14.34% under CG, CHG, and CHH, respectively (Fig. 1A). The corresponding methylation rates were also analyzed, and the results showed that mCG and mCHH methylation rate are sensitive to drought stress and re-watering process, while mCHG rate did not change much (range from 29.01 to 29.62%) in all four conditions we tested (Fig. 1B).

To visualize the global pattern in each sample, a heat map on all 12 chromosomes of rice was performed including the information about the DNA methylation levels under different stress stages, the distribution of genes and TEs in all sequence contexts (Fig. 1C). For all the three contexts, the methylated cytosines enriched in the pericentromeric regions where there are low levels of genes in intensity. CHH methylation level was around 10–30%, which was strikingly lower than CG and CHG contexts (Fig. 1C, Supplemental Figure S2), which were about 70–100% and 10–100%, respectively.

DNA Methylation Involved in Rice Drought Responses and Drought Memory

DMRs were identified to compare methylation levels after drought stress using sliding-windows approach. Under all contexts of DMRs we detected, CHH DMRs were most abundant (Table 1). We found 44,415 and 45,205 DMRs under CHH sequence in the first (S1) and subsequent (S4) drought stress, while less than 10,000 DMRs under CG (9067 for the first drought stress and 7831 for the second drought stress) or CHG (9249 for the first drought stress and 8138 for the second drought stress) were identified in these two cycles of treatment. Moreover, the number of DMRs under CHH sequence contexts increased 1.78% during the subsequence drought treatment (S4) comparing to first drought stress (S1), while the numbers of DMRs under CG and CHG decreased 13.63% and 12.01%, respectively (Table 1). The number of both hypermethylation and hypomethylation in CG and CHG DMRs were decreased in subsequent drought stress (S4) than first stress treatment (S1), although we found that most of the regions in genome of DNA methylation changes were different (Table 1). Interestingly, the number of hypermethylation in CHH DMRs was increased from 19,669 to 39,388 (54.50%) dramatically. This result suggests that cytosine methylation participated in rice short-term drought memory and CHH methylation is very sensitive to drought stress.

Table 1 Information of DMRs under drought stress
Full size table

Next, we examined the distribution of DMRs in genomic features such as TEs, 2-kb upstream of genes (upstream), gene bodies (body), 2-kb downstream of genes (downstream), and intergenic regions (IG) (Fig. 2). Most CG-related DMRs were located at IG and 2-kb upstream (promoter region) of genes, while most CHG- and CHH-related DMRs were located in TE region and IG. Overall, more hypermethylation were found than hypomethylation were found in both R0–S1 and R3–S4, regardless CG-, CHG-, or CHH-related DMRs. Comparing R0–S1 and R3–S4, CHH-related hypermethylation increased dramatically (54.50%), which TE located DMRs contributed the most (83.11% of the total). On the contrary, in CG- and CHG-related hypermethylation, TE located DMRs decreased slightly (4.28% and 6.45%) (Fig. 2). Less hypomethylated DMRs were found in subsequent drought stress than the first drought treatment regardless sequence contexts, with 21.39%, 29.58%, and 40.12% decrease under CG, CHG, and CHH, respectively (Fig. 2). Meanwhile, hypermethylated DMRs showed complicated patterns. These results showed that CHG-related methylation did not respond much comparing to CG- and CHH-related hypermethylation. Most interestingly, TE located methylation was found to be most sensitive in responding to the sequential drought stress (Fig. 2).

Fig. 2
figure 2

Distribution of drought DMRs under initial and subsequently drought stresses. A Numbers of hypermethylation DMRs under the contexts of CG, CHG, and CHH among TE, 2-kbs upstream of a gene, gene body, and 2-kbs downstream of a gene, and IGs under two drought stresses. B Numbers of hypomethylation DMRs under the contexts of CG, CHG, and CHH among TE, 2-kbs upstream of a gene, gene body, and 2-kbs downstream of a gene, and IGs under two drought stresses

Full size image

To better compare the DNA methylation under the initial and subsequently drought responses, we performed GO analyses of DMR-related genes obtained from both treatments. In total, 22,605 and 22,506 DMR-related genes were identified in the first and the subsequent drought stress, respectively. DMR-related genes were categorized into three groups: biological process, cellular component, and molecular function. Our results showed that in both stress cycles, DMR-related genes were associated with catalytic activity, component binding, metabolic process, cellular process, and other processes (Supplemental Fig. S3). Notably, although we did not observe a significant difference between two drought stress treatments in GO analysis, a lot of DMR-related genes enrich in membrane-bounded organelle and in biological process of response to stimulus including abiotic stimulus.

Rice Drought Stress Memory-Related DMRs Showed Dynamic and Distinct Patterns

Previously, we defined memory genes as those transcript levels in subsequent stress were significantly different from their levels during the initial drought stress period. Among 10,124 memory genes we identified, 5373 memory transcripts were identified to be possibly regulated by DNA methylation (Li et al. 2019). To further explore DMRs patterns in rice related to drought stress memory, we defined memory DMRs based on the following criteria: among regions that were responsive to the first drought stress in genome, those regions for which DNA methylation levels in subsequent stress periods were significantly different from their levels during the first stress period (R0–S1). STEM program was utilized to cluster the memory DMRs according their methylation status. In all, 10 out of 26 clusters of DMRs represented memory DMRs (Profile 0, 1, 3, 4, 7, 18, 21, 22, 24, 25), which memory patterns were classified as lineage, stable, accumulated, dosage dependent, and Initial memory (Fig. 3A). The rest DMRs do not follow the memory patterns which identified in R0–S1, R0–R3, R0–R4 were classified as non-memory DMRs.

Fig. 3
figure 3

Drought-Memory DMR profiles in rice. A CG methylation drought-memory DMR profiles under different drought and re-watering treatment stages. B Distribution of memory DMRs in each profile

Full size image

Our results showed that 27.73% and 26% DMRs under sequence contexts CG and CHH were clustered as memory DMRs, respectively. While only about 22.82% DMRs under sequence contexts CHG were memory DMRs. Most of memory DMRs gathered in profile4, profile7, profile18, and profile21 (Fig. 3B), which indicated that memory DMRs initially changed (both hypomethylation and hypermethylation) in R0–S1 stage, maintained at a stable level despite subsequent drought treatment (R3–S4). We further examined the distribution of memory DMRs, comparing with non-memory DMRs identified in genomic features. In memory MDRs under CHG context, TE-related distribution is lower than in non-memory DMRs (Fig. 4A), with 21.63% in memory DMRs and 34.29% in non-memory DMRs. Meanwhile, the distribution of DMRs under CHG contexts showed similar pattern related to gene body, upstream, downstream, and IG region between memory (16.45%) and non-memory DMRs (15.22%) (Fig. 4B). Interestingly, we also noticed that under contexts of CG and CHH located at IG, the proportion of memory DMRs were 12.45% and 36.84% higher than non-memory DMRs, while the proportion of upstream and downstream of genes were less (8.70% and 14.02% less under CG and CHH upstream of genes, 10.11% and 7.20% less under CG and CHH downstream of genes) (Fig. 4B).

Fig. 4
figure 4

Drought-memory DMRs distinct from drought non-memory DMRs. A Distribution of drought memory and drought non-memory DMRs under the contexts of CG, CHG, and CHH among TE, 2-kbs upstream of a gene, gene body, and 2-kbs downstream of a gene, and IGs. B Proportion of DMRs identified in (A) under each context. M memory DMRs; N non-memory DMRs. Colors represent DMRs located in 2-kb upstream of a gene, gene body, 2-kb downstream of a gene, or IGs. C Density map of the distances of memory DMRs and non-memory DMRs to TEs (Color figure online)

Full size image

Distribution of TEs over DMRs was examined and the distances of DMRs between their nearest TEs were analyzed too. We found the overlap rates (DMRs overlap with TEs) were comparable between CG-related memory DMRs and non-memory DMRs, while the rates in CHG- and CHH-related non-memory DMRs showed higher ratios compare to their memory DMRs, respectively (38.1% in memory vs 39.3% in non-memory under CG contexts, 43.9% in memory vs 49.2% in non-memory under CHG contexts, and 67.2% vs 74.6% under CHH contexts). It is also noteworthy that there are significant differences between memory DMRs and non-memory DMRs regarding distances to TE under all sequence contexts, as peaks shifted noticeably in density plots comparing the distance between memory/non-memory DMRs and TEs (Fig. 4C). Distances between DMRs and TEs were found globally smaller in memory DMRs which suggested that memory DMR might affect TEs globally in drought memory.

GO enrichment and KEGG were performed comparing the memory and non-memory DMR-related genes to obtain a better idea of drought memory mechanism. Similar patterns of GO enrichment were found in both memory DMRs and non-memory DMRs (Supplemental Figure S4). Results of KEGG enrichment showed that both memory DMRs and non-memory DMRs regulated stilbenoid, diarylheptanoid, and gingerol biosynthesis, biosynthesis of secondary metabolites, metabolic pathways, limonene, and pinene degradation pathways (Fig. 5). Interestingly, we found alpha-Linolenic acid metabolism, linoleic acid metabolism, biosynthesis of amino acids, Glycerophospholipid metabolism, Cysteine and methionine metabolism, and Lysine biosynthesis pathways were regulated significantly by memory DMRs specifically, while non-memory DMRs specifically regulated phenylpropanoid biosynthesis, starch and sucrose metabolism, zeatin biosynthesis, phosphatidylinostiol signaling system, diterpendoid biosynthesis, and plant-pathogen interaction pathways. Overall, our results suggested that memory DMRs may play important roles in drought memory and regulate pathways for plant coping with repeated drought stress.

Fig. 5
figure 5

KEGG enrichment analysis of DMR-related genes under recurrent drought stresses. A KEGG enrichment analysis of drought-memory DMRs-related genes. B KEGG enrichment of non-memory DMRs-related genes

Full size image

Memory DMRs Could Directly Regulate Rice Drought Memory Genes

Our previous results showed that 5373 drought memory transcripts might be regulated by DNA methylation by linkage analysis on differentially DNA methylated regions and expression of memory genes (Li et al. 2019). The gene encoding the proline biosynthetic enzyme Δ1-pyrroline-5-carboxylate synthetase 1 (P5CS1) was reported to be a critical gene in rice drought stress memory (Li et al. 2019). The expression of two P5CS1 homologous (LOC_Os01g62900 and LOC_Os05g38150) were found rapidly induced after the first drought stress (S1) and reach a peak at R2, then remained stable throughout the subsequent drought stress treatment, in correspond with the level of free proline content(Li et al. 2019). In the current study, the methylation status of these two drought memory P5CS1 genes in rice were explored to test if the drought memory genes were directly regulated by memory DMRs. Our results showed that memory DMRs were located in promoter region of LOC_Os05g38150 and in gene body of LOC_Os01g62900, and both the memory DMR pattern belonged to profile 21 (Figs. 3A and 6A).

Fig. 6
figure 6

Drought-memory DMRs directly regulate drought memory genes. A Genome browser views of DMRs in two homologous genes of P5CS1 (LOC_Os05g38150 and LOC_Os01g62900). B KEGG enrichment analysis of memory genes regulated by memory DMRs. The size of each point indicates numbers of genes. Color indicates the range of P values. Rich Factor indicates the ratio of the number of memory genes regulated by memory DMRs in the pathway relative to all genes in these pathways (Color figure online)

Full size image

Although DNA methylation was found involved in drought memory genes in our previous work (Li et al. 2019), the relationship between memory genes and memory DMRs was not known. Therefore, we analyzed 663 drought memory genes which identified to have memory DMRs distribution in gene body region or 2 kb flanking region of the gene. The numbers of memory DMRs under the sequence contexts CG, CHG, and CHH were found to be 151, 94, and 484 among these drought memory genes, respectively. Next, we performed GO enrichment and KEGG pathway analysis of these 663 genes to have a better idea of how these drought-memory DMRs regulate memory genes (Supplemental Fig. S5 and Fig. 6B). Most of the memory genes were found involving in metabolic, cellular process, binding, and catalytic activity. Results of KEGG enrichment analysis showed that memory DMRs related memory genes mainly through sesquiterpenoid and triterpenoid biosynthesis, phenylpropanoid biosynthesis, arginine and proline metabolism pathways. These evidences supported memory DMRs involved in rice drought memory by directly regulate rice drought memory genes.

Discussion

Growing evidence of DNA methylation play important roles in rice response to drought has been reported in recent years. However, how DNA methylation involved in drought memory is still unclear. In this study, we assessed the dynamic patterns of DNA methylation patterns in rice cultivar under recurrent drought stresses and recovery treatments using genome-wide bisulphite sequencing at a single base resolution methylome profiling level. This method enables us to observe and explain how DNA methylation change dynamically to regulate rice response to the drought stresses that happens repeatedly in a memory way.

Different plant species may have diverse methylation patterns that occur under the contexts of CG, CHG, and CHH. In our results, most of CG sequences were methylated (80–100%) while only a small portion of CHH sequences were methylated (10–30%) of the genome. The methylation rate of mCG was 43.18%, and mCHG and mCHH methylation rates were 29.60% and 27.21% at R0 without any treatment (Fig. 1B). This result is similar with previous study showed that the percentage of methylcytosines was 46.70–49.52%, 28.36–30.72%, and 19.76–24.30% under CG, CHG, and CHH context (Rajkumar et al. 2020). The slight differences might because the differences in rice cultivar, development stages, growth, and stress treatment conditions. Similar levels of methylation were observed in forward and reverse strands under all contexts in all the samples, which are consistent with a previous report (Garg et al. 2015).

Under drought stress, many previous studies reported DNA methylation level decreased or not change (Rajkumar et al. 2020; Wang et al. 2011; Zheng et al. 2013), while other studies reported strongly increased when rice seedlings encountered drought stress (Gayacharan 2013; Xia et al. 2017). In our study, we observed the general methylation level increased in response to drought stress (R0–S1), which in coincidence with the latter reports. Furthermore, we also observed that the methylation level increased again in the subsequently drought stress (R3–S4) when the methylation level recovered after the initial drought stress. This indicates that the methylation status of rice in responses to drought stress is high dynamic throughout the duration of stress and the recovery time, and also may vary among different rice genotypes growing/stressing under different conditions. This is also supported by another study showed that the methylome divergence among the cultivars is higher than it changes under drought stress within one cultivar (Rajkumar et al. 2020). Our results also suggested that methylation under CG and CHH contexts are more sensitive and responsive to drought stress than that under CHG context.

DMRs in response to drought stress were used to study the methylation dynamics. In rice, more DMRs in sensitive genotype than in tolerant line, and more hypermethylated DMRs were reported compare hypomethylated ones in both drought-sensitive and drought-tolerant lines under drought stress (Wang et al. 2016). This was also observed in our line (Zhonghua 11), which is a drought-sensitive genotype. In our results, much higher numbers of DMRs were identified under CHH context, which confirmed the important role of CHH methylation under drought stress as reported previously (Rajkumar et al. 2020). DMRs identified in our report was significantly higher than the previous reports, maybe due to the cultivar and the way of detecting its methylation status and data analysis (whole-genome bisulphite sequencing vs. MSAPs, MeDIP-seq), also might because that the time points before and after both stress treatments were included in our study. Our results shown that subsequently drought stress (S4) resulted less hypomethylated DMRs compare to the initial drought treatment under contexts of all sequences (CG, CHG, and CHH), regardless of their locations (TEs, gene body, or Igs). These results suggested that hypomethylation is more stable in constantly drought stressed plant. More interestingly, the hypermethylated DMRs in our study showed different patterns: with increasing numbers after subsequently drought stress under CHH context while the numbers of hypermethylated DMRs under CG and CHG were only slightly decreased or stable. Together with previous report showed that more DMRs under CHH context were found in response to drought (Rajkumar et al. 2020), our results confirmed that methylation occur under different contexts were differentially influenced by constantly drought stress, and CHH-related DNA methylation might affect TEs globally, and take important roles in drought memory.

Although drought memory genes in rice were well described, the mechanism on how they were regulated still needs to be explored and drought-memory DMRs were rarely studied (Auler et al. 2017; Li et al. 2019; Liu et al. 2016). In our study, drought stress DMRs showed dynamic and distinctive patterns when treat with recurrent drought stresses. We also defined patterns of memory DMRs to better reflect the ability and effects of DNA methylation on rice drought memory, 10 out of 26 clusters of DMRs were identified as memory DMRs (Profile 0, 1, 3, 4, 7, 18, 21, 22, 24, 25), according to their methylated status which classified as lineage, stable, accumulated, dosage dependent, and initial memory patterns. It is known that methylation occurs at diverse sites could affect different changes in gene regulation. Methylations in promoter regions are generally considered negatively associated with gene expression, while DNA methylation on genes body can cause both silencing or activation of genes. Meanwhile, methylation located on Igs may regulate alternative or antisense transcription. In our results, most DMRs were found targeted to TEs and least were targeted to gene bodies, which showed the important roles of DMRs on TEs to regulate drought memory. Comparing memory DMRs with non-memory DMRs, similar distribution patterns were found under the CG, CHG, and CHG contexts. However, distances of drought-memory DMRs to TEs were found significantly different when comparing their target sites. Under all context, the distances between memory DMRs to TEs were found significantly smaller than non-memory DMRs, with CHH-related ones showed the nearest comparing with CG and CHG. Similar distribution patterns were found between CG and CHG related Non-memory DMRs, while CHH-related DMRs distributed much nearer. Our results revealed the close relationship of TEs and drought-memory DMRs, and support the idea that epigenetic factors play a crucial role in drought memory formation. It would be interesting to reveal how they specifically regulate the methylation in drought memory in the future as CHH sites are governed by DRM1/2 and CMT2.

Previously, we detected 6885 transcripts involved in drought memory response during recurring drought stresses. Other than transcriptional memory, short-term epigenetic memory was discovered in the current report. GO terms and KEGG pathway analysis of the memory DMRs- and non-memory DMRs-related genes showed that they enriched in the same pathway biosynthesis of secondary metabolites, metabolic pathways, limonene and pinene degradation pathways, while memory DMRs specifically enriched in alpha-Linolenic acid metabolism, linoleic acid metabolism, biosynthesis of amino acids, Glycerophospholipid metabolism, Cysteine and methionine metabolism, and Lysine biosynthesis pathways. Previously, we found that 5373 memory transcripts might be regulated by DNA methylation. In this report, we cross checked the memory DMRs-related genes and find 663 drought memory genes have memory DMRs distribution, which enriched in sesquiterpeoid and triterpenoid biosynthesis, phenylpropanoid biosynthesis, arginine and proline metabolism pathways. The two proline biosynthetic enzymes P5CS1 in rice also gave solid evidence that memory DMRs could directly regulate rice drought memory genes. These results further confirmed that short-term epigenetic memory DNA methylation participated in rice seedlings drought stress memory by regulating genes, in addition to affecting TEs. Epigenetic memory is known could inherit and transmit to subsequent generations after exposed to drought stress conditions and increase their tolerance to drought and survival rate. How is the relationship between short-term DNA methylation memory and epigenetic memory which inherited to the next generations? How DNA methylation and histone acetylation work in concert to regulate memory formation? Those would be interesting questions for further studies.

Conclusion

Sessile plants experience constant drought challenges throughout their lives. Acute responses to drought stress have been studied worldwide, however, the drought memory mechanisms which are critical for plant coping with recurring and repetitive drought stress are less well understood. Through genome-wide bisulphite sequencing, our study provided new insight on how DNA methylation dynamic changed and showed that short-term epigenetic memory probably regulated the global expression of genes and TEs in rice drought memory.