BACKGROUND
an increasing number of documents have evidenced the importance of dna methylation in the process of plant growth and development, plant defense, and the response to adversity stress, including salt, drought and heavy metal stress.
naturally occurring modifications in a single gene locus in plants may yield heritable morphological changes without alteration of the dna sequence  <cit> , and some of these modifications may even be passed down through several generations.
arabidopsis, a model plant, is the first plant for which a whole-genome methylome map was deciphered, for its compact genome and rapid cycle .
cotton has always been known as a primary model plant for studying genome polyploidization  <cit> , cell fate determination, cell elongation and cell wall formation  <cit> .
currently, the genomes of diploid cottons g. arboreum  and g. raimondii  and complex allotetraploid g. hirsutum l.  have all been completed , and high-quality reference genomes and improving sequencing technologies may better ensure the whole-genome methylome analyses of cotton and some other moderate-sized crop genomes, such as the recent discovery of the tomato fruit methylome and of hypomethylation in the rice endosperm  <cit> .
here, we investigate whether whole-genome epigenome reprogramming occurs during the biologically important processes of cotton in response to drought stress.
the response to drought stress in plants is a complicated process, involving several genes and metabolic networks, and dna methylation is one of them.
publications have shown that drought can induce alterations of the dna methylation locus and patterns in plants with variety specificity, tissue specificity and stress specificity .
the establishment and maintenance of dna methylation patterns that result in gene expression modulation is one of the key steps in epigenetic regulation during normal growth and developmental programs.
documents about transcriptome analysis and drought-response mirna identification and functional analysis have been publicated .
recent discoveries have identified long non-coding rnas  as new important players in dna methylation regulation  <cit> .
a recent study published in nature showed a novel lncrna arising from the cebpa gene locus  that is critical for regulating dna methylation at this site through the binding of eccebpa with dna methyltransferase <dig>  dnmt <dig>  <cit> .
interestingly, a robust increase was observed in the levels of dna methylation of the cebpa promoter region following the depletion of this non-coding transcript  <cit> .
however, the alterations of the dna methylome and the possible roles of lncrnas in regulating the occurrence of methylation are still unclear under drought conditions in cotton.
here, we provide a high-resolution and high-coverage map comprising the methylation status of individual cytosines throughout the cotton genome in response to drought stress.
on the basis of the possible regulatory functions of lncrnas, and to investigate the roles of lncrnas in adjusting the epigenome modifications, we integrated the dna methylome and lncrnas.
based on this resource, future studies can use low-coverage methylome sequencing to determine the impact of differentially methylated regions on gene expression, chromosome biology, and transgenerational inheritance.
for example, this will allow breeders to determine the contribution of epigenetic modifications to phenotypic variations  <cit> , and the breeders can utilize a form of “epigenetic selection” to analogously select the individuals with desired epigenomic patterns in breeding programs.
RESULTS
genome-wide patterns of dna methylation
bisulfite treatment of dna, termed the “gold standard” of dna methylation research, can convert unmethylated cytosines into uracils while leaving methylated cytosines unchanged, which allows the generation of “methylomes” at single-base resolution  <cit> .
the bisulfite-converted dna of drought-tolerant cotton variety zhongh <dig> under different treatments, including ck, drought  reached approximately 7% in pots) , and re-watering, was sequenced with illumina hi-seq instruments at a coverage that provided information on the methylation states of individual cytosines with high confidence.
the percentage of methylated cytosines varied depending on the local sequence context  and the external treatments .
the results of wgbs are listed in table  <dig>  we obtained  <dig> , <dig> and  <dig> , <dig> and  <dig> , <dig> reads after ck, drought and re-watering stresses, respectively.
in addition, the mapping rate was  <dig> ,  <dig>  and  <dig> %, respectively, showing that the method and the results have high reliability and accuracy.
the distribution of methylation sites in each chromosome was analyzed .
among the methylated cytosines, we found that more than half of them were located in asymmetric chh contexts .
the total mc in the whole genome under three environments  was  <dig> ,  <dig>  and  <dig> %, respectively, showing an up-regulated methylation level after drought stress and a re-down-regulated trend after recovery from drought stress by re-watering compared with ck.
the duplication rate  increased approximately 2% after drought stress and re-watering , which showed that repetitive sequences may be important in response to drought stress through the alterations of the methylation level and state.
interestingly, we found that methylated cytosines in mcpg, mchg, mchh contexts all showed a hyper-methylation pattern after drought stress compared with ck and re-watering, but the alterations in the chh contexts were more significant than the alterations in other contexts.
this finding suggested that methylation levels in asymmetric chh contexts were dynamically changing with external environments and were mostly correlated with environments.
in symmetric cg and chg contexts, more than half of these cytosines were methylated while only a small proportion were methylated in the chh context .fig.
<dig> epigenome of gossypium hirsutum l. a morphological changes of cotyledons under different treatments.
b density plot of 5-methylcytosine in sequence contexts .
the circles from the outside to inside represent mcg, mchg, mchh, respectively
mc
mcg
mchg
mchh
total reads are the clean sequencing reads after the data filtration.
mapping rate are the reads with unique positions on reference genome after alignments.
mc, mcg, mchg, mchh percent  are the percentages of methylated mc in each context in all methylated cytosines from the whole genome
duplication rate is the percentage of repetitive sequences in all clean sequencing reads.
mc percent  is the percentage of methylated cytosines in all cytosines from the whole genome.
mcg, mchg, mchh percent  are the percentages of methylated mcg, mchg, mchh in all corresponding cytosines context  in genome, respectively
cytosine dna methylation may have sequence preference
studies showed that the occurrence of cytosine methylation was associated with its nearest sequence context  <cit> .
therefore, we used a tool named logo plots, which explores the sequence information of a methylated site or its nearby regions, to help us to understand the preference of methylated sites.
the method is described in the “methods” section.
among all methylated mc sites, the sequence preference around mc sites varied with the cg context and non-cg context .
in the symmetric cg context, mc frequently occurred at tcga sequences, regardless of the presence of low-methylation-level sites or high-methylation-level sites.
in the chg context, we found that the mc of high methylation levels occurred in ctg sequences, whereas at low methylation levels, mc appeared in ctg or cag sequences, which may be caused by drought stress and re-watering treatments.
however, in the chh context, low-methylation-level mc sites often occurred in the caa sequence context, and high-methylation-level mc sites occurred in the cta context.
therefore, we inferred that most high-methylation-level mc sites regularly occurred in the ctn  context, while low-methylation-level mc sites occurred in the can  context.
studies in humans suggest that chg cytosine in the tacag context is more easily methylated, and this type of methylation is always found at splicing sites  <cit> .
along with the alterations of methylation levels, the frequency of cg dinucleotide and its nearest sequences is also changed.
however, this is not always consistent in all creatures.
in humans, the methylation extent of the cpg dinucleotide is not closely related to its nearest sequence context; however, in terms of chh and chg, t and a are always upstream and downstream of cytosines with a higher methylation level  <cit> .fig.
<dig> sequence preference of cytosine dna methylation and differentially methylated cytosines in gene elements.
weblogo analysis of bases around 5m-c in different sequence contexts.
x axis represents bases position and y axis represents the degree of base enrichment.
the position of c was defined as  <dig>  , ,  represent cg, chg and chh context, respectively.
a, b, c represent all mc sites, high methylation sites and low methylation sites, respectively.
in cg context, we defined that methylation level of high methylation level sites was >75% and others were low methylation level sites.
in non-cg context, we defined that methylation level of high methylation level sites was >25% and others were low methylation level sites.
methylation level on different genome elements.
x axis represents different gene elements, and y axis represents methylation levels.
each element of each gene was divided into  <dig> bins and the average of all bins in each gene element was obtained.
different colours represent different sequence contexts
differentially methylated regions  responding to drought stress
based on the annotation of gene elements   <cit> , we compared the methylation levels in each context in each gene element , which can help to identify the functions of methylation alterations in each element in response to drought stress.
we found that promoters displayed significantly higher levels of methylation than did exons and introns, suggesting that the methylation levels of c may be correlated with gene elements.
in each gene element, the frequencies of each context, including cg methylation, chg methylation and chh methylation, were disparate, and cg methylation accounted for approximately half of the total mc, which suggested that cg methylation was the most important methylation pattern in the three gene elements.
previous publications have showed that methylation most frequently occurs in the so-called cpg islands in the 5’ regulatory gene regions   <cit> .
we then compared the methylation levels of different contexts in promoter, exon and intron pre- and post-drought stress, and minuscule changes were found, except for chh methylation .
interestingly, we found that the chh methylation in the promoter region exhibited an up-regulated pattern under drought stress compared with the control and a slightly down-regulated pattern after the re-watering treatment compared with drought stress.
this also suggested that the chh methylation level is dynamic with environments and that chh methylation may be very closely correlated with drought stress.
studies have shown that methylation in the chh context increases compared with the cg and chg context, and this chh methylation can be explained by genes promoting de novo methylation on flanking intergenic chromatin   <cit> .
consequently, we inferred that the mechanism through which methylation is altered in promoter regions can be correlated with genes promoting de novo methylation to regulate gene expression in response to drought stress.
increasing evidence has shown that the methylation and demethylation processes are always correlated with altered gene expression  <cit> .
therefore, we analyzed the influences of methylation changes on gene body and promoter regions, and the results are shown in .
to analyze the relationship between dna methylation and gene expression, the fpkm-log_ratio was used to uniformize the fpkm value of genes and to avoid a large range of gene expression levels.
the results indicated that the range of methylation changes in the gene body region was large, while the promoter was small.
furthermore, the overall methylation level in the promoter was higher than that in the gene body region, which may be closely related to the gene functions.
combined with methylation variations and gene expression changes in the gene body and promoter regions, we found that gene expression could be influenced by dna methylation changes.
in the gene body region, the range of gene expression alterations was narrow compared with those of promoters, and there were some additional genes with a great difference in expression in the promoter region, which indicated that methylation in the gene body may influence gene expression levels to a relatively low and intensive level, while a higher methylation level in a promoter may yield a scattered scope of gene expression.
taken together, our results showed that methylation in the gene body and promoter region can suppress gene expression to some extent.fig.
<dig> differentially methylated regions and go analysis.
a methylation and expression level on different genome elements.
x axis represents different treatments , and y axis represents methylation levels or fpkm-log_ratio of genes.
red was methylation level and green was gene expression level.
the calculation method used in the process was fpkm_log_ratio = log10)/log10).
b number of dmrs under drought stress.
c go enrichment analysis of dmrs.
go analysis contains three branches: cellular component, molecular function and biological process.
red represents biological process and blue represents molecular function.
*represents statistical significance level at p =  <dig>
to assess the role of methylation in response to drought stress, we examined differentially methylated regions  under drought stress compared with controls, which were always considered to participate in the process of gene regulation.
we discovered  <dig>  dmrs under drought stress compared with the control, and the number of dmrs decreased to  <dig>  after the re-watering treatment , which indicated that some dmrs were dynamically changed along with the external environments and that, when the adversity was removed, the methylation sites would return to the normal pattern.
we also performed a dmr length analysis in each chromosome  and found that the length of dmrs varied among  <dig> chromosomes, which may be correlated with the gene numbers, transposon numbers and certain other elements in each chromosome.
to determine the role of dmrs during the response to drought stress, we performed gene ontology term enrichment analysis of differentially methylated genes .
genes involved in protein binding, adp binding, defense response and transporter activity were highly enriched.
moreover, a go analysis of hyper- and hypo- methylated genes showed that these genes with great changes in methylation levels or methylation patterns may be closely correlated with drought , in turn suggesting that methylation alterations may be a mechanism to regulate gene expression in response to drought stress.
the kegg pathway analysis results showed that pathways such as ribosomes, rna degradation, protein processing in the endoplasmic reticulum, and plant hormone signal transduction may be closely correlated with the response to drought .
above all, we inferred that plants can change the methylation levels or patterns of certain genes in some pathways related to stress to regulate gene expression, which is a mechanism to respond to stress.
long non-coding rnas  may mediate the occurrence of dna methylation
dna methylation was the first discovered epigenetic mechanism  <cit>  and was reported to be regulated by lncrnas  <cit> .
to assess the role of lncrnas in regulating dna methylation changes in response to drought stress, we performed whole-genome screening and identification of lncrnas responding to drought, and we obtained  <dig>  lincrnas,  <dig> anti-sense lncrnas and  <dig> intronic lncrnas through five strict filtration steps .
the regulatory potential of lncrnas has never ceased to amaze: from rna catalysis to rna-mediated splicing, rna-based silencing, and rna-directed dna methylation  <cit> .
lncrnas are always the precursors of micrornas and play a role by breaking into one or more micrornas.
therefore, we conducted microrna prediction of these identified lncrnas with a threshold  , and the results indicated that approximately 40% of lncrnas were the precursors of micrornas.fig.
<dig> lncrna results and the prediction of lncrna-dmr related genes.
a the filtration steps and products of lncrnas.
b precursors prediction of lncrnas.
c repetitive sequences and transposons analysis of lncrnas.
d go analysis of lncrna-dmr related genes by prediction
in plants, dna methylation typically occurs by rna-directed dna methylation , which directs transcriptional gene silencing of transposons and endogenous transgenes, and rddm is driven by non-coding rnas   <cit> .
publications have reported that de novo dna methylation can be directed by many long non-coding rnas  and small interfering rnas   <cit> .
we also performed repetitive sequence and transposon analyses of lncrnas , as many 20– <dig> nucleotide sirnas were derived from repetitive sequences and transposons.
consequently, we speculated that it may be a regulatory mechanism through which long non-coding rnas mediate the occurrence by splicing into micrornas during the drought response process.
based on lncrna-seq data, we screened the dmrs with a corrected p-value ≥  <dig> , and  <dig> target genes were identified with high confidence.
the go analysis of these genes is presented in fig.
4d.
dna methylation may participate in the regulation of plant hormones
based on the analysis of go and kegg pathways, we discovered  <dig> differentially methylated genes associated with hormones , of which,  <dig>   <dig>   <dig> and  <dig> were related to ethylene, auxin, gibberellin and cytokinin, respectively.
to further investigate whether dna methylation is involved in the regulation of hormone-related genes, we used methylation-specific pcr  to test the methylation changes pre- and post-drought treatment.
nearly  <dig> %  of genes were induced to change their methylation state by drought stress , of which we found that approximately  <dig> % ,  <dig> % ,  <dig> %  and  <dig> %  of genes related to ethylene, gibberellin, cytokinin and auxin were induced to be de-methylated, which was consistent with the result of the content variations of plant hormones .
to further examine the regulatory roles of dna methylation on the hormone-related genes, we conducted the expression analysis of  <dig> genes with qrt-pcr, and the results showed that  <dig> %  of these genes exhibited a marked change in expression, including  <dig> genes that were up-regulated and  <dig> genes that were down-regulated , which were almost in agreement with the methylation changes in the process.
zhong et al.
showed that differentially methylated cytosines  were mainly concentrated in regions 5′ upstream of genes and were therefore more likely to be associated with promoter regulatory regions  <cit> .
to determine whether the epigenetic variations were associated with promoter regions, we screened the locations of these differentially methylated cytosines  through blast analysis between control and drought-treated samples.
we observed the enticement of dmcs  in promoter regions, which were similar to the methylation variations reported in arabidopsis thaliana leaves  <cit> .
consequently, this finding might indicate that dna methylation in the promoter regions of hormone-related genes is most likely to be involved in the regulation of plant hormones , which is a possible mechanism for controlling plant hormones in the drought response and may contribute to the drought resistance of cotton.fig.
<dig> regulation of plant hormones may be associated with dna methylation.
a venn diagram of differentially expressed genes associated with hormones.
dc, rd and rc represent differentially expressed hormone genes between drought and control, re-watering and drought, and re-watering and control.
b number of genes with methylation-state changed.
c the content variations of plant hormones under drought stress and re-watering treatment.
d expression analysis of plant hormone-related genes with qrt-pcr
methyltransferase inhibition could enhance drought resistance
it has been supposed that dna methylation contributes to the response to drought stress in cotton because the alterations in methylation status would change the transcription and expression of hormone-related genes.
to determine whether the methylation changes in these genes participate in regulating the response to drought , we injected 5-azacytidine, a general inhibitor of dna  methyltransferases, into cotyledons before drought stress treatment and discovered improved resistance in the treated plants.
the cotyledons of the 5-azacytidine-treated cotton plants were all still green, while most of the untreated or dd2o-treated cotyledons turned yellow and even began to fall off under the stress .
gene expression analysis revealed the substantial decrease in the methyltransferase-related lncrnas xloc_ <dig>  xloc_ <dig> and xloc_ <dig> and the sharp increase in their targets, cotad_ <dig>  cotad_ <dig> and cotad_ <dig>  which also confirmed the reliability of the method.
furthermore, the methylation-specific pcr  results confirmed that  <dig> %  of the hormone-related genes identified were induced to methylation changes , of which  <dig> % ,  <dig> % , 50%  and  <dig> %  of the genes related to ethylene, gibberellin, cytokinin and auxin, respectively, were de-methylated.
expression analysis showed that more genes presented an up-regulated trend by a large margin .
interestingly, one gene  associated with auxin was discovered to be regulated through hyper- or hypo-methylation changes, and the position of the gene is shown in fig.
6d.
next, we conducted methylation site prediction for cotad_ <dig>  and two sites were found  .
based on the observations that  more than two-thirds  of hormone-related genes epigenetically changed under drought stress and that  methylation inhibitor treatment could enhance the drought tolerance by delaying the speed of falling and turning yellow of cotyledons, we proposed that the inhibition of dna methyltransferase removes the transcript and expression constraint that prevents massive expression of some hormone-related genes during the drought stress response.fig.
<dig> application of a methylation inhibitor to enhance drought resistance.
a morphological changes of cotyledons between ck and 5-azacytidine-treated under drought stress.
ck represents cotyledons without any treatments; drought represents cotyledons under drought stress; ddh2o and 5-azacytidine represent cotyledons injected with ddh2o and 5-azacytidine, respectively.
two cotyledons from each treatment were selected and all showed a similar symptom.
b number of genes methylation state changed of  <dig> hormone- related genes.
c expression analysis of plant hormone-related genes with qrt-pcr.
d the location of two auxin-related genes.
e prediction of methylation sites of target gene cotad_ <dig>  criteria: size >  <dig>  gc percent >  <dig> , obs/exp >  <dig>
discussion
drought stress, together with other adversities, such as salt, cold and biotic factors, generates serious challenges to plant growth and development.
hence, plants have developed remarkable capabilities to modulate the physiological and molecular machinery through genome-wide gene-expression changes in response to these environmental perturbations  <cit> .
this study provides insights into the potential of a dynamic epigenome in the drought response in gossypium hirsutum l. our data showed that drought stress could induce an up-regulated epigenome, in which three sequence contexts, including mcpg, mchg, mchh contexts, all showed a hyper-methylation pattern after drought stress, which was consistent with previous research  <cit> .
furthermore, the methylation level decreased to some extent but was still slightly higher in re-watered cotton seedlings than in the control, which suggested that some methylation variations could be retained by memory and even inherited by the next few generations and that many methylation variations changed with the treatments; when the stress was removed, these methylation sites could be restored to their original state.
thus, drought-induced epigenetic changes in the cotton genome can be considered a very important regulatory mechanism for cotton plants to adapt to drought and possibly other environmental stresses.
moreover, more significant changes were found in asymmetric chh contexts than in cg and chg contexts, changing dynamically with the external environments, which suggested that chh methylation may be mostly correlated with environments.
publications have shown that chh dna methylation/demethylation may constitute a potential novel epigenetic modification mechanism that regulates growth performance in higher plants over the stress period  <cit> .
thus, we deduced that methylation/demethylation in chh contexts may be correlated with external environments and different growth stages, comprising a complex epigenetic regulatory pathway with other epigenetic modifications, such as histone acetylation and histone methylation.
based on the fact that dna methylation is always associated with nucleotide sequences  <cit> , we found that the preference for dna methylation is possibly associated with the symmetric sequence context of cg and chg and sequence regions with a high or low density of dna methylation.
many lines of evidence have proven that epigenetic modifications, such as dna methylation and histone modification, play a crucial role in regulating gene expression in response to various environmental stresses in plants  <cit> .
dna methylation state alterations in different regions  cause different gene expression changes.
unlike the cg context and the chg context, chh context methylation in promoter regions presents hyper-methylation under drought stress and then recovers to normal levels, which shows that chh context methylation dynamically changes and is most closely related to the environments.
publications have proven that methylation in the chh context promotes de novo methylation on flanking intergenic chromatin   <cit> .
thus, the disproportionally high levels of chh relative to cg and chg that we found in this study suggested that a skewed ratio of de novo methylation near genes would be beneficial for the drought response.
furthermore, methylation alterations in promoter regions would lead to a greater change in expression than those in other gene elements, such as exons and introns, which may be because methylation in promoter regions would influence the transcription of genes with micrornas and rna polymerases.
studies in humans have also shown that dna methylation is an important epigenetic mechanism for gene silencing and cancer progression and that aberrant methylation is mainly found in cpg dinucleotides within promoter regions, which is an important pathway for the repression of gene transcription  <cit> .
on the basis of our observations that  the methylation status of  <dig> selected hormone-related genes changed according to the external environment ,  the hormone content and gene expression level of these  <dig> genes also changed with treatment by drought and re-watering, and  methylation inhibitor injection induced the enhancement of drought tolerance and the alterations of methylation status, we proposed that the methylation status of hormone-related genes is one of the mechanisms responding to drought stress.
studies have shown that hormones  are associated with drought tolerance and that selecting and using cultivars with a higher proline, aba, zr and aux content under drought stress would be a practical approach to improve drought tolerance in plants  <cit> .
expression analysis of these  <dig> hormone-related genes in the investigation also confirmed these results.
combined with the expression data of these hormone-related genes, the long non-coding rnas and micrornas, we inferred that long non-coding rnas would regulate the occurrence of methylation by splicing into micrornas  in hormone content regulation in response to drought stress.
karakulah et al.
have reported that both lncrnas and micrornas could block the binding of specific transcript leading to increase mrna expression by a newly discovered regulatory mechanism called endogenous target mimicry   <cit> .
based on our results and those of previous studies, we propose a four-component model for hormone regulation, in which, through interacting mechanisms that remain unclear, long non-coding rnas can direct dna methylation through splicing into micrornas in response to drought stress .
angiosperms seem to use their epigenomic stability as a supplemental regulatory mechanism controlling a developmental transition that could cause severe energy and substance waste if it were to transit from normal growth to response to various stresses.fig.
<dig> model of regulation of hormone-related genes in drought response
CONCLUSIONS
this study provides insights into the regulatory network of the dynamic epigenome in response to drought.
in addition, previous searches for the “cryptic regulatory role” that promotes responding to stresses and is of significant agricultural value were focused on the expression of stress-tolerance genes, stress-related molecular markers and hormone signaling .
indeed, plant breeding projects depend on dna-base molecular markers  and may overlook much epigenetic variation.
hence, our studies highlight the significance of epigenetic modifications, which are beneficial for future crop-improvement strategies and technologies that consider not only genetic variations but also epigenetic modifications.
during the process of responding to drought stress, the dmrs associated with drought tolerance, as well as the hormone-related genes and transcription factors described here, provide an initial set of targets for analyzing epigenomic variations across the identified drought-tolerant cotton varieties and for assessing variations that may form the basis of future expanded selection strategies.
moreover, other pathways, for example, metabolic pathways, signal transduction and ripening, may involve methylation variations, which will require additional experimental data for confirmation.

 an increasing number of documents have evidenced the importance of dna methylation in the process of plant growth and development, plant defense, and the response to adversity stress, including salt, drought and heavy metal stress. naturally occurring modifications in a single gene locus in plants may yield heritable morphological changes without alteration of the dna sequence, and some of these modifications may even be passed down through several generations. arabidopsis, a model plant, is the first plant for which a whole-genome methylome map was deciphered, for its compact genome and rapid cycle. cotton has always been known as a primary model plant for studying genome polyploidization, cell fate determination, cell elongation and cell wall formation. currently, the genomes of diploid cottons g. arboreum and g. raimondii and complex allotetraploid g. hirsutum l. have all been completed, and high-quality reference genomes and improving sequencing technologies may better ensure the whole-genome methylome analyses of cotton and some other moderate-sized crop genomes, such as the recent discovery of the tomato fruit methylome and of hypomethylation in the rice endosperm. here, we investigate whether whole-genome epigenome reprogramming occurs during the biologically important processes of cotton in response to drought stress. the response to drought stress in plants is a complicated process, involving several genes and metabolic networks, and dna methylation is one of them. publications have shown that drought can induce alterations of the dna methylation locus and patterns in plants with variety specificity, tissue specificity and stress specificity. the establishment and maintenance of dna methylation patterns that result in gene expression modulation is one of the key steps in epigenetic regulation during normal growth and developmental programs. documents about transcriptome analysis and drought-response mirna identification and functional analysis have been publicated. recent discoveries have identified long non-coding rnas as new important players in dna methylation regulation. a recent study published in nature showed a novel lncrna arising from the cebpa gene locus that is critical for regulating dna methylation at this site through the binding of eccebpa with dna methyltransferase dnmt. interestingly, a robust increase was observed in the levels of dna methylation of the cebpa promoter region following the depletion of this non-coding transcript. however, the alterations of the dna methylome and the possible roles of lncrnas in regulating the occurrence of methylation are still unclear under drought conditions in cotton. here, we provide a high-resolution and high-coverage map comprising the methylation status of individual cytosines throughout the cotton genome in response to drought stress. on the basis of the possible regulatory functions of lncrnas, and to investigate the roles of lncrnas in adjusting the epigenome modifications, we integrated the dna methylome and lncrnas. based on this resource, future studies can use low-coverage methylome sequencing to determine the impact of differentially methylated regions on gene expression, chromosome biology, and transgenerational inheritance. for example, this will allow breeders to determine the contribution of epigenetic modifications to phenotypic variations, and the breeders can utilize a form of “epigenetic selection” to analogously select the individuals with desired epigenomic patterns in breeding programs.  genome-wide patterns of dna methylation bisulfite treatment of dna, termed the “gold standard” of dna methylation research, can convert unmethylated cytosines into uracils while leaving methylated cytosines unchanged, which allows the generation of “methylomes” at single-base resolution. the bisulfite-converted dna of drought-tolerant cotton variety zhongh under different treatments, including ck, drought reached approximately 7% in pots), and re-watering, was sequenced with illumina hi-seq instruments at a coverage that provided information on the methylation states of individual cytosines with high confidence. the percentage of methylated cytosines varied depending on the local sequence context and the external treatments. the results of wgbs are listed in table  we obtained, and, and, reads after ck, drought and re-watering stresses, respectively. in addition, the mapping rate was, and %, respectively, showing that the method and the results have high reliability and accuracy. the distribution of methylation sites in each chromosome was analyzed. among the methylated cytosines, we found that more than half of them were located in asymmetric chh contexts. the total mc in the whole genome under three environments was, and %, respectively, showing an up-regulated methylation level after drought stress and a re-down-regulated trend after recovery from drought stress by re-watering compared with ck. the duplication rate increased approximately 2% after drought stress and re-watering, which showed that repetitive sequences may be important in response to drought stress through the alterations of the methylation level and state. interestingly, we found that methylated cytosines in mcpg, mchg, mchh contexts all showed a hyper-methylation pattern after drought stress compared with ck and re-watering, but the alterations in the chh contexts were more significant than the alterations in other contexts. this finding suggested that methylation levels in asymmetric chh contexts were dynamically changing with external environments and were mostly correlated with environments. in symmetric cg and chg contexts, more than half of these cytosines were methylated while only a small proportion were methylated in the chh context.fig. epigenome of gossypium hirsutum l. a morphological changes of cotyledons under different treatments. b density plot of 5-methylcytosine in sequence contexts. the circles from the outside to inside represent mcg, mchg, mchh, respectively mc mcg mchg mchh total reads are the clean sequencing reads after the data filtration. mapping rate are the reads with unique positions on reference genome after alignments. mc, mcg, mchg, mchh percent are the percentages of methylated mc in each context in all methylated cytosines from the whole genome duplication rate is the percentage of repetitive sequences in all clean sequencing reads. mc percent is the percentage of methylated cytosines in all cytosines from the whole genome. mcg, mchg, mchh percent are the percentages of methylated mcg, mchg, mchh in all corresponding cytosines context in genome, respectively cytosine dna methylation may have sequence preference studies showed that the occurrence of cytosine methylation was associated with its nearest sequence context. therefore, we used a tool named logo plots, which explores the sequence information of a methylated site or its nearby regions, to help us to understand the preference of methylated sites. the method is described in the “methods” section. among all methylated mc sites, the sequence preference around mc sites varied with the cg context and non-cg context. in the symmetric cg context, mc frequently occurred at tcga sequences, regardless of the presence of low-methylation-level sites or high-methylation-level sites. in the chg context, we found that the mc of high methylation levels occurred in ctg sequences, whereas at low methylation levels, mc appeared in ctg or cag sequences, which may be caused by drought stress and re-watering treatments. however, in the chh context, low-methylation-level mc sites often occurred in the caa sequence context, and high-methylation-level mc sites occurred in the cta context. therefore, we inferred that most high-methylation-level mc sites regularly occurred in the ctn context, while low-methylation-level mc sites occurred in the can context. studies in humans suggest that chg cytosine in the tacag context is more easily methylated, and this type of methylation is always found at splicing sites. along with the alterations of methylation levels, the frequency of cg dinucleotide and its nearest sequences is also changed. however, this is not always consistent in all creatures. in humans, the methylation extent of the cpg dinucleotide is not closely related to its nearest sequence context; however, in terms of chh and chg, t and a are always upstream and downstream of cytosines with a higher methylation level.fig. sequence preference of cytosine dna methylation and differentially methylated cytosines in gene elements. weblogo analysis of bases around 5m-c in different sequence contexts. x axis represents bases position and y axis represents the degree of base enrichment. the position of c was defined as,, represent cg, chg and chh context, respectively. a, b, c represent all mc sites, high methylation sites and low methylation sites, respectively. in cg context, we defined that methylation level of high methylation level sites was >75% and others were low methylation level sites. in non-cg context, we defined that methylation level of high methylation level sites was >25% and others were low methylation level sites. methylation level on different genome elements. x axis represents different gene elements, and y axis represents methylation levels. each element of each gene was divided into bins and the average of all bins in each gene element was obtained. different colours represent different sequence contexts differentially methylated regions responding to drought stress based on the annotation of gene elements, we compared the methylation levels in each context in each gene element, which can help to identify the functions of methylation alterations in each element in response to drought stress. we found that promoters displayed significantly higher levels of methylation than did exons and introns, suggesting that the methylation levels of c may be correlated with gene elements. in each gene element, the frequencies of each context, including cg methylation, chg methylation and chh methylation, were disparate, and cg methylation accounted for approximately half of the total mc, which suggested that cg methylation was the most important methylation pattern in the three gene elements. previous publications have showed that methylation most frequently occurs in the so-called cpg islands in the 5’ regulatory gene regions. we then compared the methylation levels of different contexts in promoter, exon and intron pre- and post-drought stress, and minuscule changes were found, except for chh methylation. interestingly, we found that the chh methylation in the promoter region exhibited an up-regulated pattern under drought stress compared with the control and a slightly down-regulated pattern after the re-watering treatment compared with drought stress. this also suggested that the chh methylation level is dynamic with environments and that chh methylation may be very closely correlated with drought stress. studies have shown that methylation in the chh context increases compared with the cg and chg context, and this chh methylation can be explained by genes promoting de novo methylation on flanking intergenic chromatin. consequently, we inferred that the mechanism through which methylation is altered in promoter regions can be correlated with genes promoting de novo methylation to regulate gene expression in response to drought stress. increasing evidence has shown that the methylation and demethylation processes are always correlated with altered gene expression. therefore, we analyzed the influences of methylation changes on gene body and promoter regions, and the results are shown in. to analyze the relationship between dna methylation and gene expression, the fpkm-log_ratio was used to uniformize the fpkm value of genes and to avoid a large range of gene expression levels. the results indicated that the range of methylation changes in the gene body region was large, while the promoter was small. furthermore, the overall methylation level in the promoter was higher than that in the gene body region, which may be closely related to the gene functions. combined with methylation variations and gene expression changes in the gene body and promoter regions, we found that gene expression could be influenced by dna methylation changes. in the gene body region, the range of gene expression alterations was narrow compared with those of promoters, and there were some additional genes with a great difference in expression in the promoter region, which indicated that methylation in the gene body may influence gene expression levels to a relatively low and intensive level, while a higher methylation level in a promoter may yield a scattered scope of gene expression. taken together, our results showed that methylation in the gene body and promoter region can suppress gene expression to some extent.fig. differentially methylated regions and go analysis. a methylation and expression level on different genome elements. x axis represents different treatments, and y axis represents methylation levels or fpkm-log_ratio of genes. red was methylation level and green was gene expression level. the calculation method used in the process was fpkm_log_ratio = log10)/log10). b number of dmrs under drought stress. c go enrichment analysis of dmrs. go analysis contains three branches: cellular component, molecular function and biological process. red represents biological process and blue represents molecular function. *represents statistical significance level at p =  to assess the role of methylation in response to drought stress, we examined differentially methylated regions under drought stress compared with controls, which were always considered to participate in the process of gene regulation. we discovered dmrs under drought stress compared with the control, and the number of dmrs decreased to after the re-watering treatment, which indicated that some dmrs were dynamically changed along with the external environments and that, when the adversity was removed, the methylation sites would return to the normal pattern. we also performed a dmr length analysis in each chromosome and found that the length of dmrs varied among chromosomes, which may be correlated with the gene numbers, transposon numbers and certain other elements in each chromosome. to determine the role of dmrs during the response to drought stress, we performed gene ontology term enrichment analysis of differentially methylated genes. genes involved in protein binding, adp binding, defense response and transporter activity were highly enriched. moreover, a go analysis of hyper- and hypo- methylated genes showed that these genes with great changes in methylation levels or methylation patterns may be closely correlated with drought, in turn suggesting that methylation alterations may be a mechanism to regulate gene expression in response to drought stress. the kegg pathway analysis results showed that pathways such as ribosomes, rna degradation, protein processing in the endoplasmic reticulum, and plant hormone signal transduction may be closely correlated with the response to drought. above all, we inferred that plants can change the methylation levels or patterns of certain genes in some pathways related to stress to regulate gene expression, which is a mechanism to respond to stress. long non-coding rnas may mediate the occurrence of dna methylation dna methylation was the first discovered epigenetic mechanism and was reported to be regulated by lncrnas. to assess the role of lncrnas in regulating dna methylation changes in response to drought stress, we performed whole-genome screening and identification of lncrnas responding to drought, and we obtained lincrnas, anti-sense lncrnas and intronic lncrnas through five strict filtration steps. the regulatory potential of lncrnas has never ceased to amaze: from rna catalysis to rna-mediated splicing, rna-based silencing, and rna-directed dna methylation. lncrnas are always the precursors of micrornas and play a role by breaking into one or more micrornas. therefore, we conducted microrna prediction of these identified lncrnas with a threshold, and the results indicated that approximately 40% of lncrnas were the precursors of micrornas.fig. lncrna results and the prediction of lncrna-dmr related genes. a the filtration steps and products of lncrnas. b precursors prediction of lncrnas. c repetitive sequences and transposons analysis of lncrnas. d go analysis of lncrna-dmr related genes by prediction in plants, dna methylation typically occurs by rna-directed dna methylation, which directs transcriptional gene silencing of transposons and endogenous transgenes, and rddm is driven by non-coding rnas. publications have reported that de novo dna methylation can be directed by many long non-coding rnas and small interfering rnas. we also performed repetitive sequence and transposon analyses of lncrnas, as many 20– nucleotide sirnas were derived from repetitive sequences and transposons. consequently, we speculated that it may be a regulatory mechanism through which long non-coding rnas mediate the occurrence by splicing into micrornas during the drought response process. based on lncrna-seq data, we screened the dmrs with a corrected p-value ≥ , and target genes were identified with high confidence. the go analysis of these genes is presented in fig. 4d. dna methylation may participate in the regulation of plant hormones based on the analysis of go and kegg pathways, we discovered differentially methylated genes associated with hormones, of which, and were related to ethylene, auxin, gibberellin and cytokinin, respectively. to further investigate whether dna methylation is involved in the regulation of hormone-related genes, we used methylation-specific pcr to test the methylation changes pre- and post-drought treatment. nearly % of genes were induced to change their methylation state by drought stress, of which we found that approximately %, %, % and % of genes related to ethylene, gibberellin, cytokinin and auxin were induced to be de-methylated, which was consistent with the result of the content variations of plant hormones. to further examine the regulatory roles of dna methylation on the hormone-related genes, we conducted the expression analysis of genes with qrt-pcr, and the results showed that % of these genes exhibited a marked change in expression, including genes that were up-regulated and genes that were down-regulated, which were almost in agreement with the methylation changes in the process. zhong et al. showed that differentially methylated cytosines were mainly concentrated in regions 5′ upstream of genes and were therefore more likely to be associated with promoter regulatory regions. to determine whether the epigenetic variations were associated with promoter regions, we screened the locations of these differentially methylated cytosines through blast analysis between control and drought-treated samples. we observed the enticement of dmcs in promoter regions, which were similar to the methylation variations reported in arabidopsis thaliana leaves. consequently, this finding might indicate that dna methylation in the promoter regions of hormone-related genes is most likely to be involved in the regulation of plant hormones, which is a possible mechanism for controlling plant hormones in the drought response and may contribute to the drought resistance of cotton.fig. regulation of plant hormones may be associated with dna methylation. a venn diagram of differentially expressed genes associated with hormones. dc, rd and rc represent differentially expressed hormone genes between drought and control, re-watering and drought, and re-watering and control. b number of genes with methylation-state changed. c the content variations of plant hormones under drought stress and re-watering treatment. d expression analysis of plant hormone-related genes with qrt-pcr methyltransferase inhibition could enhance drought resistance it has been supposed that dna methylation contributes to the response to drought stress in cotton because the alterations in methylation status would change the transcription and expression of hormone-related genes. to determine whether the methylation changes in these genes participate in regulating the response to drought, we injected 5-azacytidine, a general inhibitor of dna methyltransferases, into cotyledons before drought stress treatment and discovered improved resistance in the treated plants. the cotyledons of the 5-azacytidine-treated cotton plants were all still green, while most of the untreated or dd2o-treated cotyledons turned yellow and even began to fall off under the stress. gene expression analysis revealed the substantial decrease in the methyltransferase-related lncrnas xloc_ xloc_ and xloc_ and the sharp increase in their targets, cotad_ cotad_ and cotad_ which also confirmed the reliability of the method. furthermore, the methylation-specific pcr results confirmed that % of the hormone-related genes identified were induced to methylation changes, of which %, %, 50% and % of the genes related to ethylene, gibberellin, cytokinin and auxin, respectively, were de-methylated. expression analysis showed that more genes presented an up-regulated trend by a large margin. interestingly, one gene associated with auxin was discovered to be regulated through hyper- or hypo-methylation changes, and the position of the gene is shown in fig. 6d. next, we conducted methylation site prediction for cotad_ and two sites were found. based on the observations that more than two-thirds of hormone-related genes epigenetically changed under drought stress and that methylation inhibitor treatment could enhance the drought tolerance by delaying the speed of falling and turning yellow of cotyledons, we proposed that the inhibition of dna methyltransferase removes the transcript and expression constraint that prevents massive expression of some hormone-related genes during the drought stress response.fig. application of a methylation inhibitor to enhance drought resistance. a morphological changes of cotyledons between ck and 5-azacytidine-treated under drought stress. ck represents cotyledons without any treatments; drought represents cotyledons under drought stress; ddh2o and 5-azacytidine represent cotyledons injected with ddh2o and 5-azacytidine, respectively. two cotyledons from each treatment were selected and all showed a similar symptom. b number of genes methylation state changed of hormone- related genes. c expression analysis of plant hormone-related genes with qrt-pcr. d the location of two auxin-related genes. e prediction of methylation sites of target gene cotad_ criteria: size >  gc percent > , obs/exp >  discussion drought stress, together with other adversities, such as salt, cold and biotic factors, generates serious challenges to plant growth and development. hence, plants have developed remarkable capabilities to modulate the physiological and molecular machinery through genome-wide gene-expression changes in response to these environmental perturbations. this study provides insights into the potential of a dynamic epigenome in the drought response in gossypium hirsutum l. our data showed that drought stress could induce an up-regulated epigenome, in which three sequence contexts, including mcpg, mchg, mchh contexts, all showed a hyper-methylation pattern after drought stress, which was consistent with previous research. furthermore, the methylation level decreased to some extent but was still slightly higher in re-watered cotton seedlings than in the control, which suggested that some methylation variations could be retained by memory and even inherited by the next few generations and that many methylation variations changed with the treatments; when the stress was removed, these methylation sites could be restored to their original state. thus, drought-induced epigenetic changes in the cotton genome can be considered a very important regulatory mechanism for cotton plants to adapt to drought and possibly other environmental stresses. moreover, more significant changes were found in asymmetric chh contexts than in cg and chg contexts, changing dynamically with the external environments, which suggested that chh methylation may be mostly correlated with environments. publications have shown that chh dna methylation/demethylation may constitute a potential novel epigenetic modification mechanism that regulates growth performance in higher plants over the stress period. thus, we deduced that methylation/demethylation in chh contexts may be correlated with external environments and different growth stages, comprising a complex epigenetic regulatory pathway with other epigenetic modifications, such as histone acetylation and histone methylation. based on the fact that dna methylation is always associated with nucleotide sequences, we found that the preference for dna methylation is possibly associated with the symmetric sequence context of cg and chg and sequence regions with a high or low density of dna methylation. many lines of evidence have proven that epigenetic modifications, such as dna methylation and histone modification, play a crucial role in regulating gene expression in response to various environmental stresses in plants. dna methylation state alterations in different regions cause different gene expression changes. unlike the cg context and the chg context, chh context methylation in promoter regions presents hyper-methylation under drought stress and then recovers to normal levels, which shows that chh context methylation dynamically changes and is most closely related to the environments. publications have proven that methylation in the chh context promotes de novo methylation on flanking intergenic chromatin. thus, the disproportionally high levels of chh relative to cg and chg that we found in this study suggested that a skewed ratio of de novo methylation near genes would be beneficial for the drought response. furthermore, methylation alterations in promoter regions would lead to a greater change in expression than those in other gene elements, such as exons and introns, which may be because methylation in promoter regions would influence the transcription of genes with micrornas and rna polymerases. studies in humans have also shown that dna methylation is an important epigenetic mechanism for gene silencing and cancer progression and that aberrant methylation is mainly found in cpg dinucleotides within promoter regions, which is an important pathway for the repression of gene transcription. on the basis of our observations that the methylation status of selected hormone-related genes changed according to the external environment, the hormone content and gene expression level of these genes also changed with treatment by drought and re-watering, and methylation inhibitor injection induced the enhancement of drought tolerance and the alterations of methylation status, we proposed that the methylation status of hormone-related genes is one of the mechanisms responding to drought stress. studies have shown that hormones are associated with drought tolerance and that selecting and using cultivars with a higher proline, aba, zr and aux content under drought stress would be a practical approach to improve drought tolerance in plants. expression analysis of these hormone-related genes in the investigation also confirmed these results. combined with the expression data of these hormone-related genes, the long non-coding rnas and micrornas, we inferred that long non-coding rnas would regulate the occurrence of methylation by splicing into micrornas in hormone content regulation in response to drought stress. karakulah et al. have reported that both lncrnas and micrornas could block the binding of specific transcript leading to increase mrna expression by a newly discovered regulatory mechanism called endogenous target mimicry. based on our results and those of previous studies, we propose a four-component model for hormone regulation, in which, through interacting mechanisms that remain unclear, long non-coding rnas can direct dna methylation through splicing into micrornas in response to drought stress. angiosperms seem to use their epigenomic stability as a supplemental regulatory mechanism controlling a developmental transition that could cause severe energy and substance waste if it were to transit from normal growth to response to various stresses.fig. model of regulation of hormone-related genes in drought response  this study provides insights into the regulatory network of the dynamic epigenome in response to drought. in addition, previous searches for the “cryptic regulatory role” that promotes responding to stresses and is of significant agricultural value were focused on the expression of stress-tolerance genes, stress-related molecular markers and hormone signaling. indeed, plant breeding projects depend on dna-base molecular markers and may overlook much epigenetic variation. hence, our studies highlight the significance of epigenetic modifications, which are beneficial for future crop-improvement strategies and technologies that consider not only genetic variations but also epigenetic modifications. during the process of responding to drought stress, the dmrs associated with drought tolerance, as well as the hormone-related genes and transcription factors described here, provide an initial set of targets for analyzing epigenomic variations across the identified drought-tolerant cotton varieties and for assessing variations that may form the basis of future expanded selection strategies. moreover, other pathways, for example, metabolic pathways, signal transduction and ripening, may involve methylation variations, which will require additional experimental data for confirmation. 
our data also showed that three sequence contexts, including mcpg, mchg, mchh, all presented a hyper-methylation pattern under drought stress and were nearly restored to normal levels after the re-watering treatment. furthermore, the many hormone-related genes with methylation variations suggested that plant hormones might be a potential mechanism in the drought response. dna methylation, with a cryptic role in genome stability, gene transcription and expression, is involved in the drought response process in plants, but the complex regulatory mechanism is still largely unknown. combined with the targets of long non-coding rnas, we found that long non-coding rnas may mediate variations in methylation patterns by splicing into micrornas. future crop-improvement strategies may benefit by taking into account not only the dna genetic variations in cotton varieties but also the epigenetic modifications of the genome. among all the methylation variations, asymmetric chh methylation was the most consistent with external environments, suggesting that methylation/demethylation in a chh context may constitute a novel epigenetic modification in response to drought stress. here, we performed whole-genome bisulfite sequencing and identified long non-coding rnas on cotton leaves under drought stress and re-watering treatments. we obtained and differentially methylated regions after drought stress and re-watering treatments, respectively. methylomesupland cottonepigenome modificationslong non-coding rnasdrought stressthe national 13th five-year major science and technology projects 2016yfd0101006ye wuwei issue-copyright-statement© the author 2017
 dna methylation, with a cryptic role in genome stability, gene transcription and expression, is involved in the drought response process in plants, but the complex regulatory mechanism is still largely unknown.  here, we performed whole-genome bisulfite sequencing and identified long non-coding rnas on cotton leaves under drought stress and re-watering treatments. we obtained and differentially methylated regions after drought stress and re-watering treatments, respectively. our data also showed that three sequence contexts, including mcpg, mchg, mchh, all presented a hyper-methylation pattern under drought stress and were nearly restored to normal levels after the re-watering treatment. among all the methylation variations, asymmetric chh methylation was the most consistent with external environments, suggesting that methylation/demethylation in a chh context may constitute a novel epigenetic modification in response to drought stress. combined with the targets of long non-coding rnas, we found that long non-coding rnas may mediate variations in methylation patterns by splicing into micrornas. furthermore, the many hormone-related genes with methylation variations suggested that plant hormones might be a potential mechanism in the drought response.  future crop-improvement strategies may benefit by taking into account not only the dna genetic variations in cotton varieties but also the epigenetic modifications of the genome. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords methylomesupland cottonepigenome modificationslong non-coding rnasdrought stressthe national 13th five-year major science and technology projects 2016yfd0101006ye wuwei issue-copyright-statement© the author 2017 