Genome-Wide DNA Methylation Dynamics During Drought Responsiveness in Tibetan Hulless Barley

Differences in drought stress tolerance within diverse grass genotypes have been attributed to epigenetic modifications. DNA methylation is an important epigenetic alteration regulating responses to drought stress. However, its effects on drought tolerance are poorly understood in Tibetan hulless barley. Here, bisulfite sequencing was conducted to profile the DNA methylation patterns of drought-tolerant variety (XL) and drought-sensitive (DQ) under drought and control conditions. A total of 5843 million reads were generated. We found the significant genome-wide changes in CHH methylation rates between XL and DQ, while CG or CHG methylation rates did not. Besides that, the two contrasting varieties do reveal distinct responses to drought stress in antioxidant activities and differentially methylated regions (DMRs). Genes in drought-tolerant varieties XL are rapidly and significantly methylated when exposed to drought stimulus. These DMRs-related genes in XL are significantly enriched in defense response and response to stimuli via gene-ontology enrichment analysis. Then, we focused on 1003 transcription factors and identified 15 specific DMRs-related transcription factors exhibiting specific methylation changes under drought stimuli. Finally, we identified three DMRs-related TFs (HvRR12, HvRR2, and HvCSP41B), where Arabidopsis homologs involve in responses to drought conditions. Altogether, abiotic stresses could be rapidly respond and mediated by methylation of transcription factors in hulless barely.


Introduction
Hulless barley (Hordeum vulgare L. var. nudum Hook. f.), also called naked barley, is an important cereal crop in Tibet Plateau (Xu et al. 2016). It has been served as a healthy food for human consumption and animal feed for over thousands of years. Owing to high altitudes, the naked barley is cultivated in harsh environment such as valleys and higher land on Tibet (Liang et al. 2017). It is also affected by drought and low temperature in March every year when the weather is quite cold and dry (J. B. Du et al. 2011). In order to mitigate adversities such as drought, salinity, and low temperature, various strategies have been evolved in hulless barley (H. Li et al. 2014).
It is well known that drought is the one of the most serious environmental stress that affect crop growth and yield by causing a wide range of physiological and biochemical responses (Flowers 1989;Iqbal et al. 2015). One of inevitable biochemical changes occurring is the accumulation of reactive oxygen species (ROS), which is scavenged by antioxidative enzymes Superoxide dismutase (SOD) Catalase (CAT). Many genes responsible for drought tolerance have been reported for decades (Lenka et al. 2015;Moa et al. 2005). For example, the transcription factor AtHB13 has been reported to act as a positive regulator of drought tolerance in Arabidopsis. Genetically, drought tolerance is extremely complex trait involving in several genetic pathways such as polygenic control and complex morpho-physiological mechanisms (Xuekun Zhang et al. 2014). At the molecular level, drought stress could induce genome-wide changes in gene expression via epigenetic mechanisms like histone modification and DNA methylation in plants Zong et al. 2013). For example, the mutant allele (met1) DNA methyltransferase in one locus could remove methylation at several genomic regions and leading to specific express of 31 stress response-related genes in tobacco (Ku et al. 2006;Zubko et al. 2012).
DNA methylation, an epigenetic modification, plays a crucial role in plant growth and development as well as responses to various abiotic stresses (Arthur et al. 2018;Uthup et al. 2011). Actually, DNA methylation exists in all eukaryotes, and the most common DNA modification methylation of cytosine (C) at the 5-position to produce 5-methylcytosine (5mC) (Kasai & Kawai 2009). Under normal conditions, the proportion of methylcytosine in plants is 20-30% in plants, which usually occurs in three nucleotide sequences: CG, CHG, and CHH (H indicates C, T, or A) (Serre, Lee, & Ting, 2009). CG methylation is produced by the conserved DNA methyltransferase METHYLTRANS-FERASE1 (MET1) (Zubko et al. 2012); CHG methylation is modified by the plant-specific DNA methyltransferase CHROMOMETHYLASE3 (CMT3) (J. Du et al. 2012); and the de novo CHH is methylated by 24-nucleotide small interfering RNA-dependent DNA methylation (RdDM) pathway (Ma et al. 2015). Interestingly, changes in methylation levels contribute greatly to the process of adaptation to stress in plants. The hyper or hypomethylation changes in the hybrids could be an indicator of the expression levels of stress-related gene under drought, when compared to their parents (Boyko et al. 2010;Hai & Zhang 2009). This indicates that methylcytosine and its reversibility may regulate transgenerational response to stresses. An increasing number of studies revealed that the transposon-rich heterochromatic regions show heavy methylation (Melamed-Bessudo & Levy 2012). Over one-third of expressed genes are methylated within transcribed regions, while only 5% genes within their promoter regions in Arabidopsis (Xiaoyu Zhang et al. 2006). Thus, DNA methylation within genes is a common feature of eukaryotic genomes. Intriguingly, methylation of transcribed regions does not usually result in gene silencing (Xiaoyu Zhang et al. 2006). It seems to primarily occur at CG sites and appears to show moderate correlation between the level of gene-body methylation and gene expression (Su et al. 2014). DNA methylation shows the deposition of certain chromatin marks such as differentially modified histones because of the tight link between DNA methylation and histone modifications (Fuks 2005). For instance, there is a relationship between the reduce in DNA methylation levels and the changes in level of H3K4me and H3K9me modifications in DECREASE IN DNA METHYLATION1 (DDM1) mutant, which maintain cytosine DNA methylation within the heterochromatic locus regions (Sasaki et al. 2012;Zhou et al. 2013).
It is the first time that the epigenetic responses to drought between two contrasting varieties were revealed via bisulfite sequencing and RNA sequencing in Tibetan hulless barely. Here, we explored how DNA methylation is involved in drought responsiveness. Changes in the level of DNA methylation affect plant resistance to drought, especially the methylation changes of several transcription factor genes. Therefore, epigenetic changes in genome could be considered as an important regulatory mechanism for plants to adapt to drought and possibly other environmental stresses.

Plant Materials and Treatments
The Tibetan hulless barley (Hordeum vulgare L. var. nudum) varieties XL and DQ were used (Xu et al. 2021). The XL variety is highly resistant to drought stress, while the DQ is sensitive. Healthy seeds were sterilized by soaking in 2% H 2 O 2 for 40 min, and rinsed in sterile water. Finally, seeds were germinated on moistened filter paper at 16-18 °C with 18-h light/10-h dark photoperiod and a relative humidity of 80% in a growth chamber. Artificial water stress was induced with 21% polyethylene glycol (PEG) 6000 solutions to achieve an osmotic potential. The roots and leaves were harvested from 10-day-old seedlings from two genotypes under normal and drought conditions at 0 h, 4 h, and 48 h, respectively; The experimental design consisted of three biological replicates at three time points for each accession. All these tissue samples mentioned above were fast frozen in nitrogen and stored at −80 °C immediately.

Determination of the MDA Content
The level of lipid peroxidation in plant tissues was measured by determination of MDA, which was measured using the thiobarbituric acid (TBA) method (Schmedes & H?lmer, 1989). About 0.3 g leaf sample of hulless barley seedling was homogenized with a mortar and pestle in 5 mL 0.1% Trichloroacetic acid (TCA), then the homogenate was centrifuged at 10,000 g for 15 min. About 5 ml of 20% TCA containing 0.5% TBA was added to 1 ml supernatant aliquot. The mixture was heated at 95 °C for 10 min, cooled immediately, and centrifuged at 10,000 g for 15 min. The reaction products of MDA and TBA show highest absorbance at 532 nm. Finally, MDA concentration was calculated. The measurement was performed with three biological replicate.

Determination of Peroxidase (POD) and Catalase (CAT) Activity
To analyze the activities of plant antioxidant enzyme (POD, CAT), about 0.5 g frozen leaf samples were homogenized 1 3 with a mortar and pestle in ice bath at a 1:10 ratio with 100 mM phosphate buffer saline containing 1 mM EDTA and 1% polyvinylpyrolidone (PVP). The crude homogenates were centrifuged at 10,000 g for 15 min at 4 °C. The supernatant was used to determine the POD and CAT activities.
POD activity was measured by spectrophotometer following change of absorption at 420 nm due to guaiacol oxidation (Tamas et al. 2016). Under these conditions, one unit of POD activity was defined as the amount of POD that was require for 50% inhibition of the enzymatic reaction in 1 ml enzyme extraction of per milligram of protein (Yu et al. 2009). CAT activity was measured, which based on the fact that ammonium molybdate could rapidly terminate the H 2 O 2 degradation reaction catalyzed by CAT and generate a yellow complex. This complex is monitored by absorbance at 405 nm using the spectrophotometer. One unit of CAT activity was defined as 1 mmol H 2 O 2 decomposed in one milligram of protein per second. The measurement was performed with three biological replicate.

Whole Genome Bisulfate Sequencing
Whole Genome Bisulfate Sequencing were performed according to the previously described (Allen et al. 2006) by Wuhan IGENEBOOK Biotechnology Co.,Ltd (http:// www. igene book. com). Briefly, 1ug DNA extracted from Tibetan hulless barley sample were fragmented with an ultrasonic disruptor (Bioruptor), with an average size of about 300-500 bp. In order to evaluate bisulfite conversion rates, 20 ng nonmethylated lambda DNA (251431, Promega) were spiked as controls during library preparation. Fragmented DNA were end-repaired and ligated to a fully methylated adapter by using NEXTflex™ Bisulfite-Seq Barcodes-6 Kit (51191, Bioo Scientific). EZ DNA Methylation-Gold ™ kit (D5005, Zymo Research Corp) was used for bisulfite conversion of adapter-linked DNA.After amplification by PCR, High-quality libraries, having a distribution of DNA fragments centered around 300 bps, were used for sequencing analysis using the Illumina sequencing platform. Three biological replicates were proformed in this study.

BS-Seq Data Analysis
Trimmomatic (version 0.38) was used to filter out low-quality reads. The clean reads were then mapped to reference genome (Xingquan Zeng,2015) by the BSMAP software (version: 2.9). Bisulfite conversion rates was evaluated by lambda DNA. The coverage of Cytosine bases across the genome and Cytosine-base methylation levels in various types were counted by cgmaptools software t (Macisaac et al. 2012). Methylation level was determined by dividing the number of reads covering each mC by the total number of reads covering that cytosine, which was also equal to the mC/C ratio at each reference cytosine. Circlize software was used to do the whole genome methylation distribution analysis. . DMRs were calculated by gmaptools software (Song et al. 2009). DMR-associated genes were defined as genes with the closet DMR located within 2 KB upstream of the transcription start site and 2 KB downstream of the transcription end site. GO and KEGG annotation were done by Eggnog5 and clusterProfiler in R package.

RNA-Seq Data Analysis
The published transcriptome data (Xu et al. 2021) under the same experimental conditions are used for gene expression analysis. The adapter and low-quality reads were filtered out through cutadapt (version 1.11). Clean reads were mapped to the reference genome (Xingquan Zeng, 2015) by Hisat2 (version 2.1.0), allowing up to two mismatches. It comprised RSEM (v1.2.6) for transcript abundance estimation and normalization of expression values as FPKM (Fragments per kilobase of transcript per million fragments mapped). Differentially expressed genes were identified with DESeq2 with a filter threshold of adjusted q-value < 0.05 and |log2FoldChange|> 1.

Physiological and Biochemical Differences Between Two Contrasting Hulless Barley Varieties
Reactive oxygen species (ROS) excessively accumulates in plant tissues exposed to drought stress. To analyze physiological and biochemical changes in drought-tolerant variety (XL) and drought-sensitive one (DQ), leaves were harvested to investigate changes of Malondialdehyde (MDA) concentration, Catalase (CAT), and Peroxidase (POD) activities. During 48-h treatment. Firstly, concentration of MDA was tested in leaves, because MDA could serve as the biomarker for cellular and oxidative damage. MDA concentration significantly decreased in XL under drought stress conditions in leaves, while it enhanced sharply in DQ (Fig. 1A). This suggests XL is more tolerable to drought stress when compared with the DQ variety.
Then, POD and CAT activities are analyzed in Tibetan hulless barley leaves. POD activities obviously reduced in XL under drought conditions, while it rose in DQ (Fig. 1B). On the contrary, activity of CAT showed contrasting trends. Its activities show 1.64-folds rise in XL when compared to normal conditions, whereas it decreased by around 50% in DQ (Fig. 1C). Therefore, CAT might involve in improving drought tolerance in the drought-tolerant variety XL, but POD activities might negatively correlate to the XL in hulless barley. Also, we detected the MDA, POD, and CAT contents in roots to testify the change patterns in different plant tissues. A similar dynamical changes was observed in XL and DQ varieties, suggesting that the plant organs might exhibit identical responses to drought conditions.

DNA Methylomes of Tibetan Hulless Barley Under Drought Conditions
In order to reveal the dynamics of genome-wide DNA methylation in contrasting Tibetan hulless barley varieties, we sequenced DNA methylomes of the hulless barley exposed to drought stresses for 0 h, 4 h, and 48 h via bisulfite sequencing, respectively. Two replicates were conducted to assess biological and experimental reproducibility in DNA methylome analysis. Finally, a total of 5843 million 150 bp paired-end reads were generated from sodium bisulfitetreated DNA with 99% bisulfte conversion rate for each library. On an average, about 583 million 150 bp sequencing reads were aligned to hulless barely genome for each sample with 25-fold coverage of the genome. An average of 52% clean read pairs were aligned to genome for each library (Table S1).
Genome-wide cytosine methylation at CG, CHG, and CHH sequences was determined in all hulless barley samples. An average of 91.71% CG and 67.36% CHG was methylated in these samples, both of which show higher methylation levels. In contrast, only 3.14% CHH was methylated, i.e., the lowest level of methylation (Table S2, Fig. 2A). In addition, we determined DNA methylation patterns in 2 kb upstream of the transcription start site (TSS), gene body, and 2 kb downstream of the transcription termination site (TTS). CG and CHG were methylated at gene bodies and their flanking regions with an obvious drop around TSS and TTS sites, which are consistent with the previous findings (Q. Li et al. 2015). In contrast, CHH methylation levels in both flanking regions were higher than that in gene bodies (Fig. S1).
We further determined and analyzed the CG, CHG, CHH methylation rates in these samples. The results revealed that there were no signifigenome-wide changes in CG or CHG methylation rates between the drought-tolerant variety (XL) and the drought-sensitive one (DQ) at 4 h and 48 h, respectively (Fig. 2B). Interestingly, CHH methylation rates decreased by 6% from 4 to 48 h in drought-tolerant variety (DQ), whereas it significantly increased by 48% in sensitive one (DQ) (Fig. 2C). This suggests that CHH contents methylation positively respond to drought tolerance in hulless barley.

Analysis of Differentially Methylation Regions Within Same Hulless Barley Varieties Between Normal and Drought Conditions
Differentially methylation regions (upstream of TSS, gene body, and downstream of TTS) were selected and identified, when we compared methylation regions within the same variety between drought and normal conditions at different periods. First of all, the CG-type, CHG-type, and CHH-type DMR numbers were compared within droughttolerant genotypes XL or drought-sensitive one DQ, respectively (Table 1). The results exhibited significant differences between XL and DQ. A total of 1840 and 141 DMRs were identified in XL at 4 h and 48 h when compared between the two conditions, respectively. But, there were just about 62 The results indicate that genes in drought-tolerant varieties XL are rapidly and significantly methylated, especially in 4 h after exposed to drought conditions. Then, we analyzed the number of CGtype, CHG-type, and CHH-type DMR in XL and DQ at 4 h. There were 1104 CG-type DMR, 660 CHG-type DMR, and 76 CHH-type DMR in XL, but there were only 15 CG-type DMR, 18 CHG-type DMR, and 29 CHH-type DMR in DQ. Considering the significant difference of DMR numbers between XL and DQ, the two contrasting varieties exhibits distinct responses to drought stress by methylating DNA. Thus, DNA methylation of genes in XL might involve in response to drought stimuli. Then, these genes in XL or DQ were further analyzed by gene ontology (GO). We found that CG-type-related DMR genes in XL were markedly enriched for nucleotide binding, ADP binding, defense response, and response to stimulus (Fig. 3).

Comparative Differentially Methylation Regions Between Two Contrasting Varieties
In order to reveal molecular mechanism of drought tolerance, we compared the levels of methylation and gene expression between two contrasting varieties XL and DQ. Firstly, differentially methylation regions (DMRs) were calculated and compared between XL and DQ under normal and drought conditions at 0 h, 4 h, and 48 h, respectively. As shown in Table 2, the number of CG-type DMR slightly increased by 5.0% at 4 h when compared to 0 h between XL and DQ under normal conditions. However, it sharply increased by 40% at 4 h under drought conditions. The similar trends were observed in the number of CHG-type and CHH-type DMR. The changes of DMR number between XL and DQ under drought conditions are more obvious than that under normal conditions, suggesting drought stress could induce DNA methylation to response to abiotic stimulus quickly. Interestingly, the number of DMR in CG, CHG, and CHH sequences between two varieties reduced from 4 to 48 h under drought stress, but it increased under normal conditions. The contrasting DMR number changes between normal and drought conditions suggested that DNA methylation might regulate responsiveness to environmental stimuli. Furthermore, we conducted gene ontology (GO) analysis by using these DMR-related genes between XL and DQ exposed to drought stress. We found CG-type DMR-related genes were significantly enriched for cellular response to  stimulus, response to stress, nucleotide binding, lyase activity, isomerase activity, and hydrolase activity (Fig. 3). Then, CG-type-, CHG-type-, and CHH-type-related genes were identified and calculated between XL and DQ both in drought and normal conditions. As shown in Table S3, a number of 284 CG-type genes, 205 CHG-type genes, and 51 CHH-type genes between two genotypes were identified at 4 h, while 152 CG-type genes, 97 CHGtype genes, and 35 CHH-type genes were identified at 48 h. The results suggest that Tibetan hulless barley could quickly activate more genes to alleviate drought damage. Besides, we used published transcriptome data (Xu et al. 2021) to explore the relationship between gene expression and methylation.We performed cluster analysis of these DMR-related genes, of which more than 50% were transcribed. The results exhibit markedly opposite expression profiles of CG-type-, CHG-type-, and CHH-type-related genes between two contrasting varieties at 4 and 48 h (Fig. 4).

Fig. 3
Gene ontology analysis of CG-type (A)-and CHH-type (B)-related DMR genes between XL and DQ exposed to drought stress

Identification of Differentially Methylation Regions-Related Transcription Factors in Hulless Barley Exposed to Drought Stress
To identify differentially methylation regions-related transcription factors (TFs) in hulless barley after exposure to drought stress, we compared the bisulfite sequencing profiles of both between two contrasting varieties and within same genotypes under drought conditions and normal conditions at 0 h, 4 h, and 48 h, respectively. We mainly focused on the four groups of TF (AP2/ EREB, bZIP, NAC and MYB), which could involve in response to abiotic stresses in plants such as Arabidopsis thaliana and rice. A total of 659 TFs were identified according to nucleotide sequence analysis of TFs from Tibetan hulless barley, which is homologous to these from rice (Table S4). Of the 659 TFs, 26 unique TFs (11 MYB, 9 NAC, 5 AP2/DREB, and 1 bZIP) appeared to be differentially methylated between different comparative combinations (Table S5). Intriguingly, 46 DMR-related TFs were identified in total, but there were as many as 20 repetitive DMR-related TFs among different combinations (Table S6).
To identify DMR-related TFs that respond to drought conditions, we compared the DMR-related TFs selected from the combinations exposed to drought conditions with these from the combinations exposed to normal conditions. As a result, 15 specific DMR-related TFs that specifically exhibit significant methylation changes under drought stimuli were identified, indicating the 15 TFs are induced by drought stress (Table 3). Among these TFs, three hulless barley TFs whose gene homolog from Arabidopsis thaliana mediate responses to drought stresses (Rohit et al. 2016;Sultana et al. 2016). AtRR12, HVUL6H08680.2 (HvRR12) gene homolog from Arabidopsis thaliana, encoding type B cytokinin response regulators negatively regulate plant responses to drought; AtRR2 that is homologous to HVUL4h39100.2 (HvRR2) exhibits tolerance to drought and salt stresses; AtCSP41B with homology to HVUL2H41931.2 (HvCSP41B) is modulated by miR399f to mediate plant responses to salt, ABA, and drought. Besides that, 4 specific Arabidopsis TFs also could involve in abiotic and biotic stresses except for drought stress. ANAC102 with homolog to both HVUL7H38638.2 and HVUL3H21712.2 mediates response to low oxygen stress; AtTDR1 with homolog to HVUL1H37933.2 enhances tolerance to salt stress via activation of ABA synthesis in Arabidopsis; ATMYB71 with homolog to HVUL4 h24905.2 might regulate plant defense via promoter-based integration. Of those remaining 8 specific DMR TFs, their homologous genes in Arabidopsis encode putative proteins with unknown functions.
To reveal DNA methylation status of the above 15 specific DMR TFs that appear to be differentially methylated by the imposition of drought stress, the sequence corresponding to upstream of TSS, gene body, and downstream of TSS was subjected to analysis by using the modified version of the previous SNV-calling method (Barturen G,2013). We mainly focused on the three TFs (HvRR12, HvRR2, and HvCSP41B) described above. In gene HvRR12, the significantly differential methylation locates in the gene body, which means the cytosines were demethylated after exposure to salinity stress for 48 h in XL (Fig. 5A). In gene HvRR2, the significantly differential methylation exhibits in upstream of TSS, and the methylation levels in XL are notably lower than that in DQ after exposed to drought conditions for 4 h (Fig. 5B). In gene HvCSP41B, the significantly differential methylation locates in downstream of the TSS, which means Fig. 4 The expression profiles of CG-type (A)-, CHG-type (B)-, and CHH-type (C)-related genes between drought-tolerant varieties XL and drought-sensitive one DQ. XL drought-tolerant variety hulless barely, DQ the drought-sensitive hulless barely, CK hulless barely were exposed to normal conditions, S hulless barely were exposed to drought conditions. it is hypermethylated in stressed XL seedlings (Fig. S2). Clearly, DNA methylation in HvRR12, HvRR2, and HvC-SP41B could be affected by drought stress.

Discussion
Reactive oxygen species (ROS) significantly accumulates in plant tissues exposed to drought conditions. We investigated changes of Malondialdehyde (MDA) contents and enzymatic activity (Catalase, Peroxidase) both in leaves and roots after exposure to drought for 48 h in drought-tolerant variety XL and sensitive one DQ, respectively. Interestingly, similar dynamical changes were observed both in leaves and roots in XL and DQ. We only sequenced DNA methylomes of the leaves in hulless barley in this study. It is necessary to analyze and compare the DNA methylomes between leaves and roots, which could be considered as a potential method to uncover the response pathway to drought conditions.
We further determined and analyzed the CG, CHG, CHH methylation rates in XL and DQ. The results revealed that there were no signifcant genome-wide changes in CG or CHG methylation rates between the two varieties, while CHH methylation rates did not. Interestingly, CHH is specifically methylated by 24-nucleotide small interfering RNAdependent DNA methylation (RdDM) pathway. Then, we analyzed the number of CG-type, CHG-type, and CHH-type DMR in XL and DQ at 4 h. There were 1104 CG-type DMR, 660 CHG-type DMR, and 76 CHH-type DMR in XL, but there were only 15 CG-type DMR, 18 CHG-type DMR, and 29 CHH-type DMR in DQ. However, the CG-type DMR and CHG-type DMR sharply decline in drought-sensitive variety DQ when compared to XL, while the CHH type does not. It is possible that the extremely low CHH methylation rates result in the less CHH-type number. CHH methylation is specifically perpetuated by the combined action of siRNAs and DRM2. Thus, it indicates that the siRNAs and DRM2 could play important roles in response to drought stress. the drought-sensitive hulless barely, CK hulless barely were exposed to normal conditions, DR hulless barely were exposed to drought conditions.
In order to identify differentially methylation regionsrelated transcription factors (TF) in hulless barley after exposure to drought stress, 1003 TFs were identified according from Tibetan hulless barley, which is homologous to these from Arabidopsis thaliana and rice (Rohit et al. 2016). Finally, 15 specific DMR-related TFs were identified, which exhibit significant methylation changes under drought stimuli. The Arabidopsis homologs of 3 TFs (HVUL6H08680.2, HVUL4 h39100.2, HVUL2H41931.2) have been proved to be involved in response to drought conditions. Therefore, DNA methylation regulates responsiveness to environmental stimuli by methylating transcription factors in hulless barely.
Altogether, it is the first time that the roles of DNA methylation in drought responsiveness were revealed in Tibetan hulless barley. Significant genome-wide changes were found in CHH methylation rates between two contrasting varieties (XL and DQ). They showed obviously different response to drought stress both in DMR numbers and antioxidant activities. Gene ontology analysis indicates DMR-related genes in XL could involve in defense response and response to stimulus. Of the 1003 TFs, 26 unique TFs appeared to be differentially methylated. Intriguingly, 15 specific DMRrelated TFs exhibiting special methylation changes exposed to drought stress were identified, and 7 TFs of those might involve in response to drought stress or other abiotic stresses. Finally, we found that methylation levels of some transcription factors have changed significantly under drought-stress treatment such as HVUL6H08680.2, HVUL4 h39100.2, and HVUL2H41931.2 whose gene homolog from Arabidopsis thaliana involves in responses to drought stresses. Therefore, our results proved that hulless barely could respond to drought stress via methylation on transcription factors.