Skip to main content
Download PDF

Transcriptome and methylome changes in two contrasting mungbean genotypes in response to drought stress

BMC Genomics volume 23, Article number: 80 (2022) Cite this article

Abstract

Background

Due to drought stress, the growth, distribution, and production of mungbean is severely restricted. Previous study combining physiological and transcriptomic data indicated different genotypes of mungbean exhibited variable responses when exposed to drought stress. Aside from the genetic variation, the modifications of environmentally induced epigenetics alterations on mungbean drought-stress responses were still elusive.

Results

In this study, firstly, we compared the drought tolerance capacity at seedling stage by detecting physiological parameters in two contrasting genotypes wild mungbean 61 and cultivar 70 in response to drought stress. We found that wild mungbean 61 showed lower level of MDA and higher levels of POD and CAT, suggesting wild mungbean 61 exhibited stronger drought resistance. Transcriptomic analysis indicated totally 2859 differentially expressed genes (DEGs) were detected when 70 compared with 61 (C70 vs C61), and the number increased to 3121 in the comparison of drought-treated 70 compared with drought-treated 61 (D70 vs D61). In addition, when drought-treated 61 and 70 were compared with their controls, the DEGs were 1117 and 185 respectively, with more down-regulated DEGs than up-regulated in D61 vs C61, which was opposite in D70 vs C70. Interestingly, corresponding to this, after drought stress, more hypermethylated differentially methylated regions (DMRs) in 61 were detected and more hypomethylated DMRs in 70 were detected. Further analysis suggested that the main variations between 61 and 70 existed in CHH methylation in promoter. Moreover, the preference of methylation status alterations in D61 vs C61 and D70 vs C70 also fell in CHH sequence context. Further analysis of the correlation between DMRs and DEGs indicated in both D61 vs C61 and D70 vs C70, the DMRs in gene body was significantly negatively correlated with DEGs.

Conclusions

The physiological parameters in this research suggested that wild mungbean 61 was more resistant to drought stress, with more hypermethylated DMRs and less hypomethylated DMRs after drought stress, corresponding to more down-regulated DEGs than up-regulated DEGs. Among the three DNA methylation contexts CG, CHG, and CHH, asymmetric CHH contexts were more dynamic and prone to be altered by drought stress and genotypic variations.

Peer Review reports

Background

Drought stress is one of the major environmental factors restricting crop growth, production, and distribution with more severe damage than other environmental stresses such as heat, low temperature and salinity stress [1,2,3]. Unlike animals, when subjected to drought stress, plants cannot escape but have to develop complicated defense systems, including a series of cellular, molecular, physiological, biochemical, anatomical and morphological responses [4]. For example, in water-deficient conditions, plants maintain cell turgor through osmotic adjustment to accumulate organic solutes such as glycine betaine, proline, and sugars to adjust water potential [5]. Meanwhile, plant have evolved detoxification systems to scavenge the excessive reactive oxidative species (ROS) caused by drought stress. The antioxidant pathways involve the enzymes superoxide dismutase (SOD), peroxidase (POD), catalase (CAT), and ascorbate peroxidase (APX) as well as the non-enzymatic compounds ascorbate, carotenoids, and glutathione [6].

When plants undergo drought stress, the stimuli induce signals transduction in multiple pathways, resulting in transcriptional changes of drought-responsive genes [7, 8]. With the popularization of next generation sequencing, RNA-seq is widely used to reveal the gene expression changes when exposed to drought stress, including drought-responsive transcription factors, plant hormone related genes, and protein kinases [9,10,11,12]. In the past decades, the genes associated with drought tolerance have been studied in detail, however drought tolerance is a complicated process not only involving transcriptional alterations but also genome-wide DNA methylation changes. DNA methylation, referred to as adding methyl group at the fifth position of cytosine pyrimidine ring, is a well-studied epigenetics mechanism [13, 14]. The methylated cytosines mainly happen in three DNA contexts, CG, CHG and CHH (where H = A, C or T) [15]. The maintenance of symmetrical methylation at CG and CHG occurs by DNA Methyltransferase 1 (MET1) and Chromomethylase 3 (CMT3) during DNA replication respectively [16, 17], while the maintenance of asymmetrical methylation at CHH is not template based but through RNA-directed DNA methylation (RdDM) with Domains Rearranged Methyltransferases 1 and 2 (DRM1, DRM2) [18].

DNA methylation as an important epigenetics marker involves in plant growth and development as well as plant abiotic stress tolerance and adaptations [19,20,21]. DNA methylation patterns in plants undergo dynamic changes depending on the tissues, plant species, and the specific type of stress [22]. Decreased DNA methylation was detected in root tissues of faba bean under drought stress, while increased DNA methylation was observed in root tissues of alfafa under salt stress [23, 24]. It has been clear that epigenetic modifications such as DNA methylation may affect gene expression and further contribute to phenotypic variation in response to environmental stress [25, 26]. Study in popular found genomic alterations of DNA methylation induced by drought stress could influence expression levels of related genes [27]. The obvious correlation between differentially DNA methylated regions and gene expression were detected when under drought stress in apple and mulberry [28, 29].

Mungbean (Vigna radiata L.) is an important fast-growing grain legume crop with rich protein, folate and iron, and is widely distributed in Asian countries especially India, China, Myanmar, and Indonesia [30, 31] . Mungbean can be processed into different food varieties such as sprout, flour, noodles, and porridge, which provide nutrition and tasty flavor for human beings [32]. However, due to abiotic and biotic stress, the yield of mungbean is low, and drought stress restrictions on mungbean production is becoming more severe [33]. Modern mungbean cultivars were derived from domestication and selection of the original mungbean species in India [34]. Normally, wild species contain valuable genes and resources which tend to be disappeared in the process of domestication and selection in breeding [35, 36]. Studies indicated wild soybean possessing multiple valuable candidate genes endowing the plants with stronger tolerance to drought stress [37, 38]. In addition, wild mungbean (TC1996) has shown complete bruchid resistance compared with cultivated mungbean [39].

In this study, we compared the drought tolerance capacity at seedling stage by detecting physiological parameters in two contrasting mungbean genotypes in response to drought stress. The comparative transcriptome analysis integrated with methylome study aimed to reveal the DNA methylation pattern and gene expression variations in control and drought-stressed conditions between the two contrasting genotypes, which might provide clues for.

the potential use of wild germplasm as a drought-tolerant resource in mungbean cultivar breeding.

Results

Comparison of physiological parameters of two mungbean genotypes exposed to drought stress

C61 and C70, representing control of wild and cultivated mungbean plants, were well-watered, while D61 and D70 plants were drought-stressed. It was obvious after drought treatment, D61 and D70 plants exhibited phenotypic changes, such as small stature and lower height (Fig. 1a). Leaf samples were collected from the vegetative 1 (V1) stage seedlings of the two contrasting mungbean plants for physiological parameters determination. We measured the content of malondialdehyde (MDA), which is the biomarkers of oxidative stress [12], and the antioxidant-related enzymes including SOD, POD, and CAT (Fig. 1b). Compared with control, the MDA content was significantly increased in both drought-stressed genotypes (p < 0.01) (Fig. 1b). However, in C70 and D70, the accumulation of MDA was higher than the levels in C61 and D61. The SOD level was also increased in drought-stressed condition in both genotypes but not statistically significant. Interestingly, both the content of POD and CAT in D61 was significantly higher than C61 (p < 0.01) but there was no significant difference between C70 and D70. These data suggested wild mungbean genotype 61 presented stronger resistance when subjected to drought stress.

Fig. 1
figure 1

Phenotypic and physiological parameters changes in two mungbean genotypes under drought stress. a Phenotypes of two mungbean genotypes plants in control and drought-stressed conditions. C61 and D61 mean well-watered and drought-stressed wild mungbean 61; C70 and D70 mean well-watered and drought-stressed mungbean cultivar 70. b MDA, SOD, POD, and CAT content of two mungbean genotypes plants in control and drought-stressed conditions. ** represents p < 0.01

Full size image

Transcriptomic changes of two mungbean genotypes in response to drought stress

To further reveal the molecular bases accounted for the different performance of two genotypes when exposed to drought stress, the RNA-Seq analysis was conducted in the two genotypes in control and drought stress conditions. Among all eight samples investigated, the raw reads generated were between 40.53 million and 48.09 million, with more than 97% valid bases (Table 1). The Q30 of all eight samples was consistently over 94% (Table 1). These parameters indicated the sequencing data quality was high and could be used for further analysis. From the four pairwise comparisons, it was obvious that the number of differentially expressed genes (DEGs) was the most in D70 vs D61 comparison group, followed by C70 vs C61 (Fig. 2a), among the DEGs in these two groups, there were more upregulated genes than downregulated genes (Fig. 2b, Additional file 1: Fig. S1 ab). Totally 1117 DEGs were identified after drought stress compared with control in wild mungbean 61 (D61 vs C61), with 384 upregulated and 733 downregulated. However, only 185 DEGs were found in D70 vs C70, with 155 upregulated and 30 downregulated (Fig. 2b, Additional file 1: Fig. S1 cd). In order to validate the accuracy of gene expression data generated by RNA-Seq, we selected eight DEGs, most of which were transcription factors related to drought stress, and conducted qRT-PCR assays. The results indicated the trends of relative expression level of upregulated and downregulated genes of qRT-PCR were similar to that calculated by transcriptome sequencing (Additional file 2: Fig. S2), confirming the reliability of transcriptome data. Kyoto Encyclopedia of Genes and Genomes (KEGG: http://www.genome.jp/kegg/) analysis of DEGs in D61 vs C61 and D70 vs C70 was performed [40]. The results indicated that the highest number of transcripts were enriched in carbohydrate metabolism and signal transduction pathways, followed by amino acid metabolism and lipid metabolism (Fig. 2c, d).

Table 1 Statistics of RNA-seq for control and drought-treated samples of wild mungbean 61 and mungbean cultivar 70
Full size table
Fig. 2
figure 2

Differentially expressed genes in different comparisons and the KEGG pathway analysis. a Venn diagrams of DEGs in four pairwise comparisons. b Numbers of upregulated and downregulated DEGs in four pairwise comparisons. KEGG pathway enrichments of DEGs in D61 vs C61 (c) and D70 vs C70 (d). D61 vs C61, drought-stressed wild mungbean 61 versus well-watered; D70 vs C70, drought-stressed cultivar mungbean 70 versus well-watered

Full size image

Methylome profiles in wild and cultivated mungbean

In total, whole-genome bisulfite sequencing (WGBS) results generated 71.9–82.7 million raw reads for each sample (Table 2). After filtering the low-quality data, 70.1–80.6 million clean reads were mapped to the reference genome using Bismark software. The mapping rate ranged from 56.84 to 79.08% (Table 2). For each sodium bisulfite treated library, the unmethylated lambda DNA was used as reference for conversion rate calculation. The results showed that the conversation efficiency was over 99% in all of the samples (Table 2). Genome-wide screening revealed that 18,227,407 methylated cytosines were detected in C61, the proportion of methylated CG, CHG and CHH was 23.8, 28.1, and 48.1% respectively (Fig. 3a). After drought treatment, there was a slight change in the methylated cytosines proportions in three sequence contexts, with CHH methylation increased to 48.9%, CG and CHG methylation decreased to 23.4 and 27.7%. By contrast, the number of methylated cytosines in C70 (30,417,748) and D70 (29,403,037) was more than in wild mungbean, but the proportions of methylated CG, CHG, and CHH was similar to wild mungbean. Interestingly, the percentage of methylated CHH decreased from 47.9 to 46.9% after drought treatment in cultivated mungbean70 (Fig. 3a).

Table 2 Statistics of WGBS for control and drought stress treated samples of wild mungbean 61 and mungbean cultivar 70
Full size table
Fig. 3
figure 3

Methylation profiles in two mungbean genotypes. a The relative proportion of mCGs, mCHGs, and mCHHs in two mungbean genotypes in control and drought-stressed conditions. The level of methylation in different gene features (b) and gene body, upstream 2 K and downstream 2 K regions (c). The comparison of methylation level in different gene features in D61 vs C61 and D70 vs C70 (d). C61 and D61 mean well-watered and drought-stressed wild mungbean 61; C70 and D70 mean well-watered and drought-stressed mungbean cultivar 70

Full size image

DNA methylation patterns in CG, CHG, and CHH were also analyzed in different mungbean genomic regions as well as gene body, promoter, and downstream 2 K region (Fig. 3b,c). It was observed CG methylation was the highest across genomic regions, followed by CHG and CHH. CG and CHG methylation level showed similar trends, for example, in wild mungbean 61 high methylation was observed in promoter and it decreased in 5UTR, increased in intron and decreased again in 3UTR (Fig. 3b). The repeat region showed the highest level of methylation (Fig. 3b). In wild mungbean 61, the highest CG methylation level was observed in upstream 2 K and gene body, followed by downstream 2 K, whereas in mungbean cultivar 70 the highest methylation level was found in upstream 2 K, followed by gene body and downstream 2 K (Fig. 3c). The tends of CHG and CHH methylation changes was similar between mungbean 61 and 70 (Fig. 3c). In addition, it was obvious in D61 vs C61, the increase of CHH was mainly contributed by 5UTR, exon, 3UTR, and repeat regions, while in D70 vs C70, the decrease of CHH was mainly contributed by promoter, intron, and repeat regions (Fig. 3a,d).

Differentially methylated regions in wild and cultivated mungbean

We further compared the differentially methylated regions (DMRs) between wild and cultivated mungbean in control and drought stress conditions. Totally, we identified 12,111 hypermethylated and 6578 hypomethylated DMRs in the wild mungbean D61 vs C61, while in the cultivar mungbean D70 vs C70, the number of hypermethylated DRMs was only 4988 and hypomethylated DMRs was 14,747 (Additional file 3: Fig. S3). After drought stress, increased methylation level of DMRs in wild mungbean 61 were detected in all CG, CHG, and CHH contexts, especially in CHH context (Fig. 4a, Additional file 4: Fig. S4a). On the contrary, the decreased methylation level of DMRs were detected in D70 vs C70 in all CG, CHG, and CHH contexts, especially in CHH context (Fig. 4a, Additional file 4: Fig. S4b). Further detailed comparative analysis related to the genome-wide distribution of DMRs was conducted (Fig. 4b, Additional file 4: Fig. S4c,d). Overall, the hypermethylated DMRs or hypomethylated DMRs in D61 vs C61 and D70 vs C 70 were mainly distributed in promoter, exon, intron and repeat regions (Fig. 4b). Further analysis indicated that the main hypermethylated DMRs after drought stress in wild mungbean 61 were distributed in promoter and intron in CG; promoter, exon and intron in CHG; and promoter, intron and repeat in CHH context (Additional file 4: Fig. S4c). By contrast, the hypomethylated DMRs in mungbean 70 after drought stress mainly distributed in promoter, exon, intron, and repeat regions in all three DNA contexts (Additional file 4: Fig. S4d). The KEGG pathway analysis was conducted in order to investigate the associated biological functions and pathways of the DMRs. The results indicated in D61 vs C61, the DMRs were mainly distributed in pathways such as purine metabolism, RNA transport, pyrimidine metabolism, RNA degradation and carbon metabolism (Fig. 5a). The first two pathways were also observed in the enrichment of D70 vs C70, in addition, protein processing in endoplasmic reticulum and endocytosis were also enriched in D70 vs C70 (Fig. 5b).

Fig. 4
figure 4

Differentially methylated regions distribution in D61 vs C61 and D70 vs C70. a Methylation level distribution of differentially methylated regions (DMRs) by violin boxplots. b Number of DMRs in different regions across genome. D61 vs C61, drought-stressed wild mungbean 61 versus well-watered; D70 vs C70, drought-stressed cultivar mungbean 70 versus well-watered

Full size image
Fig. 5
figure 5

KEGG pathway enrichment of differentially methylated genes. D61 vs C61, drought-stressed wild mungbean 61 versus well-watered; D70 vs C70, drought-stressed cultivar mungbean 70 versus well-watered

Full size image

Relationship between DNA methylation status and gene expression levels

In order to investigate whether DNA methylation is regulated by gene expression, the correlation analysis of gene expression and DNA methylation was conducted. As predicted, the unexpressed genes had the highest methylation level in all gene body, promoter, and downstream 2-kb region in CG and CHG sequence contexts (Fig. 6a,b). While the lowest methylation was detected in the genes showed high expression in all regions of CG methylation and in gene body and downstream 2-kb region of CHG methylation (Fig. 6a,b). In contrast, for CHH methylation, the unexpressed genes showing the highest methylation level was observed in all regions in the wild mungbean 61, but only in gene body and downstream 2-kb region in cultivar mungbean 70 (Fig. 6c). In addition, the low expressed genes had the lowest CHH methylation level in promoter, whereas in gene body and downstream 2-kb region, the high expressed genes had the lowest CHH methylation level (Fig. 6c). We further studied relationship between DNA methylation and gene expression. Based on methylation levels, the methylated genes were divided into five groups with group first the lowest methylation level and group fifth the highest methylation level (Fig. 7, Additional file 5: Fig. S5). We found that in promoter, genes with the highest methylation levels showed the lowest expression levels in all three DNA sequence contexts in the wild mungbean 61, but only in CG and CHG in cultivar mungbean 70 (Fig. 7). In gene body, genes with the highest methylation levels showed the lowest expression levels in all CG, CHG and CHH (Additional file 5: Fig. S5), and moderately CG methylated genes showed the highest level of expression (Additional file 5: Fig. S5a).

Fig. 6
figure 6

Relationship between gene expression and DNA methylation in C61, D61, C70 and D70. DNA methylation levels distributions in upstream 2 K, gene body, and downstream 2 K by different expression levels at CG (a), CHG (b), and CHH (c) DNA contexts. C61 and D61 mean well-watered and drought-stressed wild mungbean 61; C70 and D70 mean well-watered and drought-stressed mungbean cultivar 70

Full size image
Fig. 7
figure 7

Relationship between DNA methylation and gene expression in C61, D61, C70 and D70 in promoter. Expression profiles of different methylated levels at CG (a), CHG (b) and CHH (c) were investigated. The promoter methylation levels were classified into five groups with group.1st the lowest and group.5th the highest. C61 and D61 mean well-watered and drought-stressed wild mungbean 61; C70 and D70 mean well-watered and drought-stressed mungbean cultivar 70

Full size image

Differentially methylated regions and related differentially expressed genes

In order to study the global effect of DRMs on related gene expression, we analyzed the DMRs related genes and promoters. As a result, we found in D61 vs C61 there were 504 and 362 DEGs identified as hyper- and hypomehylated DMR-associated genes, while in D70 vs C70 the corresponding DEGs number were 210 and 606 DEGs respectively (Fig. 8a). Similarly, 482 and 344 DEGs were detected as hyper- and hypomehylated DMR-associated promoters in D61 vs C61. In contrast, in D70 vs C70 there were 210 and 594 hyper- and hypomehylated DMR-associated promoters identified (Fig. 8b). The Spearman rank correlation coefficient was used to test associations between DMRs and DEGs, and Spearman’s rho was used as a measure for correlation. The results indicated in D61 vs C61, the gene body methylation was negatively correlated with gene expression (Spearman rho = − 0.19, p value = 0) (Fig. 8c). Similar result was detected in cultivar mungbean D70 vs C70 (Spearman rho = − 0.18, p value = 0) (Fig. 8d). However, there were no clear correlation between promoter methylation and gene expression in both D61 vs C61 and D70 vs C70 (Fig. 8e,f). Altogether, the results suggested DNA methylation could partially explain the differential transcript abundances of related genes.

Fig. 8
figure 8

Differentially methylated regions and related differentially expressed genes. Differentially expressed genes (DEGs) identified as hyper- and hypomethylated differentially methylated regions (DMRs)-associated genes (a) and promoters (b). Relationship between DMRs in gene body and DEGs in D61 vs C61 (c), and D70 vs C70 (d). Relationship between DMRs in promoter and DEGs in D61 vs C61 (e), and D70 vs C70 (f). D61 vs C61, drought-stressed wild mungbean 61 versus well-watered; D70 vs C70, drought-stressed cultivar mungbean 70 versus well-watered

Full size image

Discussion

In this research, firstly, we compared the responses of two mungbean genotypes after drought stress treatment at seedling stage based on physiological parameters. Recent study indicated variable responses to drought stress among different mungbean varieties according to physiological data and transcriptomic study [33]. Aside from the genetic variation, it was reported environmentally induced epigenetics alterations also could modify stress response and broaden plant phenotypic variation [26]. Therefore, in this research, we integrated physiological parameters with transcriptome and whole genome bisulfite sequencing analysis to reveal the molecular mechanism which might explain the different performance of the two genotypes when exposed to drought stress.

After drought treatment at seedling stage, the two genotypes showed visible differences compared to control. We further found that wild mungbean 61 showed significantly lower level of MDA and higher levels of POD and CAT. As lower level of MDA and increased antioxidant enzymes activity are correlated with cell membrane stability and enhanced antioxidant defense system, which could protect plants from cytotoxic effects [6]. Therefore, our finding indicated wild mungbean 61 exhibited a higher resistance to drought-stress compared with cultivar 70. The distinct phenotypic and physiological responses indicated wild mungbean 61 and cultivar 70 are two contrasting genotypes for drought tolerance. As reported in soybean, the wild germplasm possessed valuable candidate genes which made the plants more drought-tolerant [38]. Thus, we further compared the gene expression changes and DNA methylation patterns alterations from the genome scale.

From the transcriptomic data, when wild mungbean 61 and cultivar 70 were compared, 2859 DEGs were detected. The number increased to 3121 in the comparison of D70 vs D61. However, it was obvious in D61 vs C61 and D70 vs C70, the DEGs were 1117 and 185 respectively. The data indicated the inherent genetic variations existed between these two contrasting genotypes. Interestingly, after drought stress there were more DEGs in D61 vs C61, while less DEGs in D70 vs C70. Study in onion also found that in drought-tolerant genotype more DEGs detected than in drought-sensitive genotype [41]. Among the DEGs, after drought stress, similar level of up-regulated and down-regulated DEGs were found in drought-tolerant onion, and more down-regulated DEGs were found in drought-sensitive genotype [41]. Our study found more down-regulated DEGs than up-regulated in D61 vs C61, whereas opposite response was observed in D70 vs C70. Corresponding to this, after drought stress, more hypermethylated DMRs in wild mungbean 61 were detected and more hypomethylated DMRs in cultivar mungbean 70 were found. Our findings were consistent with the commonly accepted regulative relationship that DNA methylation is negatively associated with gene expression [42, 43]. The pathway enrichment of DEGs suggested that most of them were enriched in carbohydrate metabolism and signal transduction pathways. Carbohydrate metabolism were reported plays important roles in response to drought stress, the starch and sucrose metabolism is correlated with turgor pressure maintenance [44, 45]. In addition, the signal transduction pathway was also found participated in the drought stress response [44, 46].

Our study found that in C61 and C70, the proportion of methylated CG, CHG and CHH was around 24, 28, and 48% respectively, which was similar to the results of previous study using the same mungbean material as cultivar 70 [47]. After drought treatment, genome-wide changes of CHH methylation were relatively bigger than CG and CHG, with increased CHH (from 48.1 to 48.9%) observed in D61 vs C61, but decreased CHH (from 47.9 to 46.9%) was found in D70 vs C70. Similar to our finding in D61 vs C61, in cotton, more significant changes of CHH rather than CG and CHG were found after drought stress, and CHH tended to be hypermethylated [48]. In apple, a slight increased CG and CHG methylation proportions as well as a decreased CHH proportions were revealed after water deficit [29], which was consistent to our report in D70 vs C70. Further investigation confirmed the increase of CHH in D61 vs C61 was mainly contributed by 5UTR, exon, 3UTR, and repeat regions, while the decrease of CHH in D70 vs C70 was mainly contributed by promoter, intron, and repeat regions. The preference of methylation status alterations in CHH suggested asymmetric CHH changes were dynamic and probably associated with external environments [48, 49].

We further compared the epigenetic changes from genome-scale and analyzed the interactions between DNA methylation and gene expression. Our data showed high expression in gene body tended to have lower-level methylation, and non-expressed genes had higher level methylation. Vice versa, the highest methylation levels in gene body showed the lowest expression levels. Earlier report in Arabidopsis thaliana indicated loss of methylation in gene body promoted transcription of genes [50]. However, studies in rice and apple revealed that gene body methylation was commonly positively associated with gene expression [29, 51]. Previous study reported that different DNA sequence context and different genomic regions showing varied effect on gene expression [52]. In our study, in gene body the moderately CG methylated genes showed the highest level of expression, which is consistent with the reports in poplar [53, 54]. As is known, DNA methylation in promoters is likely to impede transcription [50]. In our study, in promoter, the highest methylation levels also showed the lowest expression levels in all three DNA sequence contexts in the wild mungbean 61, however, only in CG and CHG for cultivar mungbean 70. Similarly, CHH methylation levels in apple was found positively assiciated with gene expression, which was different from CG and CHG [29]. In addition, for CHH methylation, the unexpressed genes showing the highest methylation level was observed in promoter in the wild mungbean 61, but not in promoter in cultivar mungbean 70. Altogether, based on the facts that CHH mehylation and gene expression in promoter in mungbean 70 was significantly different from others, it was obvious that main variations between wild mungbean 61 and cultivar mungbean 70 existed in CHH methylation in promoter. Previously, we also found the preference of methylation status alterations in CHH in both D61 vs C61 and D70 vs C70. Taken together, our finding suggests asymmetric CHH contexts were more dynamic and prone to be altered by environmental factor changes and genotypic variations. CHH methylation, which is maintained by CMT2 through RdDM, has been proven to be dynamic and play important roles in regulating gene expression during seed development, germination, and early plant life [55,56,57]. In rice, in response to desiccation and salinity stresses, methylation levels of CHH showed the most variation between different genotypes, suggesting the important role of CHH in abiotic stress response [58]. Further analysis of the correlation between DMRs and DEGs indicated in both D61 vs C61 and D70 vs C70, the DMRs in gene body was significantly negatively correlated with DEGs. However, no significant difference was detected between DMRs and DEGs in promoter. Our results indicated DNA methylation partially contributed to gene expression regulation. In addition, in the poplar salt stress study, few DEGs were identified as different methylation genes, which suggests that the impact of DNA methylation on gene expression is limited [53]..

Conclusions

Compared with cultivar 70, wild mungbean 61 exhibited a higher resistance to drought-stress, reflecting in lower level of MDA and higher levels of SOD, POD, and CAT. Transcriptomic analysis indicated when drought-treated 61 and 70 compared with their controls, more down-regulated DEGs than up-regulated was found, which was opposite in D70 vs C70. Corresponding to this, after drought stress, more hypermethylated DMRs in 61 were detected, with more hypomethylated DMRs in 70. In addition, we found the main variations between the two contrasting genotypes existed in CHH methylation in promoter. Coincidently, the methylation status alterations in D61 vs C61 and D70 vs C70 also fell in CHH sequence context. Further analysis of the correlation between DMRs and DEGs indicated in both D61 vs C61 and D70 vs C70, the DMRs in gene body was significantly negatively correlated with DEGs.

Methods

Plant material and drought treatment

The mungbean cultivar ‘Zhonglu 1’ (germplasm accession no. VC1973A, named 70 in this study) and wild type (germplasm accession no. JP226873, named 61 in this study) were used in this study. The seeds were kindly provided by Dr. Suk-Ha Lee from Department of Plant Science and Research Institute for Agriculture and Life Sciences, Seoul National University. The germinated seeds were grown in pots in growth chamber of Qingdao Agricultural University (Qingdao, Shandong, China) at 24 ± 2 °C day and 17 ± 2 °C night under the photoperiod of 18/6 h day/night. Plants were divided into four groups: a) well-watered 70 (C70); b) drought-stressed 70 (D70); c) well-watered 61 (C61); d) drought-stressed 61 (D61), with three biological replicates in each group for physiological parameters determination, with two biological replicates in each group for transcriptome and methylome study. The well-watered groups were irrigated normally to maintain water capacity, and drought-stressed groups were withheld water since the time planted. Seedlings with the same growth stage were selected for sampling. Leaf samples were collected near to V1 stage (fully developed trifoliate at the second node) when relative water content of soil reached to 39% in drought-stressed groups, which was 69% in well-watered control. The relative water content was calculated by fresh weight subtracting dry weight, and then divided by turgid weight subtracting dry weight according to previous report [59]. The collected leaf samples were stored at − 80 °C until used for RNA and DNA extraction.

Physiological parameters determination

The oxidative stress biomarker and antioxidant-related indicators MDA (Solarbio, BC0025), SOD (Solarbio, BC0175), POD (BC0095), and CAT (Solarbio, BC0205) was detected by using assay kits and with a BioTek Cytation 1 cell imaging multimode reader (BioTek, Winooski, VT, USA). Fresh mungbean leaf tissue was collected and the measurement was performed following the manufacturer’s instructions of Beijing Solarbio Science & Technology Co., Ltd. (Beijing, China).

RNA isolation, RNA-sequencing and data analysis

Total RNA was extracted using RNAprep Pure Plant Kit (DP441, TIANGEN Biotech). The integrity of isolated RNA was assayed through the RNA 6000 Nano labchip on 2100 Agilent Bioanalyzer (Agilent Technologies, Santa Clara, CA) before RNA-Seq libraries preparation. RNA-Seq libraries were constructed using TruSeq Stranded mRNA LTSample Prep Kit (Illumina, San Diego, CA, USA) according to the manufacturer’s instructions. After quality inspection through 2100 Agilent Bioanalyzer, the prepared RNA-Seq libraries were sequenced on Illumina HiSeq X Ten platform by OE biotech Co., Ltd. (Shanghai, China). The raw reads were filtered and trimmed using Trimmomatic v0.32 [60]. The obtained clean reads were aligned to reference genome using HISAT2 [61]. Fragments Per Kilobase of transcript per Million mapped fragments (FPKM) with Cufflinks were used to calculate the expression levels of each gene, and the read count of each gene were generated by HTSeq [62]. DESeq was used to determine DEGs [63], with p-value < 0.05 and |log2 Fold change (logFC)| > 1 setting as the cutoff for significantly DEGs. KEGG pathway enrichment analysis was performed to investigate the biological functions of DEGs using R based on the hypergeometric distribution [40].

DNA extraction, WGBS, and data analysis

Genomic DNA was extracted using modified CTAB method [64]. DNA concentration was quantified using Qubit DNA BR Assay Kits (Invitrogen, Eugene, OR, USA) according to the manufacturer’s instructions. Totally, 100 ng genomic DNA spiked with 9 ng lambda DNA were sonicated into 200–300 bp fragments with Covaris S220, and then treated with sodium bisulfite using the EZ DNA Methylation-Gold Kit (Zymo Research). The spiked lambda DNA was used as an unmethylated reference for conversion efficiency calculation. The prepared libraries were sequenced on Illumina Novaseq platform by OE biotech Co., Ltd. (Shanghai, China) after quality assessment on the 2100 Agilent Bioanalyzer. Bismark software (version 0.16.3) was used for alignments of reads to a reference genome [65]. Bioconductor package DSS software was used for identification of DMRs [66]. The genes related to DMRs was defined as the genes with gene body region or promoter region have an overlap with the DMRs. GOseq R package was used for Gene Ontology (GO) enrichment analysis of genes related to DMRs [67], and KOBAS software was used to determine the statistical enrichment of DMR-related genes in KEGG pathways [68].

Quantitative real-time PCR analysis

A total of 1.5 μg RNA was used for reverse transcription to obtain complementary DNA (cDNA) using 5X All-In-One RT MasterMix (abm, China). The primers used for qRT-PCR were list in Additional file 6: Table S1. The reactions for qRT-PCR were performed based on the protocol of ChamQ SYBR Color qPCR Master Mix (Vazyme, Shanghai, China), with two biological and three technical replicates using a CFX96 instrument (Bio-Rad, Hercules, CA, USA). The total amplification volume was 10 μL per reaction, and the conditions for PCR reaction were as follows: 95 °C for 30 s, 35 cycles of 95 °C for 10 s, 53 °C for 30 s, and 72 °C for 30 s, then followed by 65 °C for 5 s and 95 °C for 5 min. The relative gene expression data was calculated using 2-ΔΔCt method.

Availability of data and materials

All related sequencing data is deposited in NCBI Sequence Read Archive (SRA) database with the link of https://www.ncbi.nlm.nih.gov/Traces/study/?acc=PRJNA771920. The bioProject accession is PRJNA771920 and BioSample accessions are from SAMN22349573 to SAMN22349580.

Abbreviations

DEGs:

Differentially expressed genes

DMRs:

Differentially methylated regions

KEGG:

Kyoto encyclopedia of genes and genomes

SOD:

Superoxide dismutase

POD:

Peroxidase

CAT:

Catalase

MDA:

Malondialdehyde

RdDM:

RNA-directed DNA methylation

MET1:

DNA Methyltransferase 1

CMT3:

Chromomethylase 3

DRM1:

Domains Rearranged Methyltransferases 1

DRM2:

Domains Rearranged Methyltransferases 2

References

  1. Beck EH, Fettig S, Knake C, Hartig K, Bhattarai T. Specific and unspecific responses of plants to cold and drought stress. J Biosci. 2007;32(3):501–10.

    Article  CAS  PubMed  Google Scholar 

  2. Boyer JS. Plant productivity and environment. Science. 1982;218(4571):443–8.

    Article  CAS  PubMed  Google Scholar 

  3. Lambers H, Chapin FS, Pons TL. Plant physiological ecology. New York: Springer New York; 2008. p.183.

  4. Joshi R, Wani SH, Singh B, Bohra A, Dar ZA, Lone AA, et al. Transcription factors and plants response to drought stress: current understanding and future directions. Front Plant Sci. 2016;7:1029.

    Article  PubMed  PubMed Central  Google Scholar 

  5. Turner NC, Abbo S, Berger JD, Chaturvedi S, French RJ, Ludwig C, et al. Osmotic adjustment in chickpea (Cicer arietinum L.) results in no yield benefit under terminal drought. J Exp Bot. 2007;58(2):187–94.

    Article  CAS  PubMed  Google Scholar 

  6. Gill SS, Khan NA, Tuteja N. Differential cadmium stress tolerance in five Indian mustard (Brassica juncea L.) cultivars: an evaluation of the role of antioxidant machinery. Plant Signal Behav. 2011;6(2):293–300.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  7. Zhu J-K. Salt and drought stress signal transduction in plants. Annu Rev Plant Biol. 2002;53(1):247–73.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  8. Hirayama T, Shinozaki K. Research on plant abiotic stress responses in the post-genome era: past, present and future. Plant J. 2010;61(6):1041–52.

    Article  CAS  PubMed  Google Scholar 

  9. Muthusamy M, Uma S, Backiyarani S, Saraswathi MS, Chandrasekar A. Transcriptomic changes of drought-tolerant and sensitive Banana cultivars exposed to drought stress. Front Plant Sci. 2016;7:1609.

    Article  PubMed  PubMed Central  Google Scholar 

  10. Singh D, Singh CK, Taunk J, Tomar RS, Chaturvedi AK, Gaikwad K, et al. Transcriptome analysis of lentil (Lens culinaris Medikus) in response to seedling drought stress. BMC Genomics. 2017;18(1):206.

    Article  PubMed  PubMed Central  Google Scholar 

  11. Shi D, Wang J, Bai Y, Liu Y. Transcriptome sequencing of okra (Abelmoschus esculentus L. Moench) uncovers differently expressed genes responding to drought stress. J Plant Biochem Biotechnol. 2019;29(2):155–70.

    Article  Google Scholar 

  12. You J, Zhang Y, Liu A, Li D, Wang X, Dossa K, et al. Transcriptomic and metabolomic profiling of drought-tolerant and susceptible sesame genotypes in response to drought stress. BMC Plant Biol. 2019;19(1):267.

    Article  PubMed  PubMed Central  Google Scholar 

  13. Bird A. DNA methylation patterns and epigenetic memory. Genes Dev. 2002;16(1):6–21.

    Article  CAS  PubMed  Google Scholar 

  14. Reinders J, Wulff BB, Mirouze M, Marí-Ordóñez A, Dapp M, Rozhon W, et al. Compromised stability of DNA methylation and transposon immobilization in mosaic Arabidopsis epigenomes. Genes Dev. 2009;23(8):939–50.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  15. Oakeley EJ, Jost J-P. Non-symmetrical cytosine methylation in tobacco pollen DNA. Plant Mol Biol. 1996;31(4):927–30.

    Article  CAS  PubMed  Google Scholar 

  16. Lindroth AM, Cao X, Jackson JP, Zilberman D, McCallum CM, Henikoff S, et al. Requirement of CHROMOMETHYLASE3 for maintenance of CpXpG methylation. Science. 2001;292(5524):2077–80.

    Article  CAS  PubMed  Google Scholar 

  17. Kankel MW, Ramsey DE, Stokes TL, Flowers SK, Haag JR, Jeddeloh JA, et al. Arabidopsis MET1 cytosine methyltransferase mutants. Genetics. 2003;163(3):1109–22.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Wassenegger M, Heimes S, Riedel L, Sänger HL. RNA-directed de novo methylation of genomic sequences in plants. Cell. 1994;76(3):567–76.

    Article  CAS  PubMed  Google Scholar 

  19. Law JA, Jacobsen SE. Establishing, maintaining and modifying DNA methylation patterns in plants and animals. Nat Rev Genet. 2010;11(3):204–20.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. Huang W, Xian Z, Hu G, Li Z. SlAGO4A, a core factor of RNA-directed DNA methylation (RdDM) pathway, plays an important role under salt and drought stress in tomato. Mol Breed. 2016;36(3):28.

    Article  Google Scholar 

  21. Sun L, Miao X, Cui J, Deng J, Wang X, Wang Y, et al. Genome-wide high-resolution mapping of DNA methylation identifies epigenetic variation across different salt stress in maize (Zea mays L.). Euphytica. 2018;214(2):1–15.

    Article  Google Scholar 

  22. Komivi D, Marie AM, Rong Z, Qi Z, Mei Y, Ndiaga C, et al. The contrasting response to drought and waterlogging is underpinned by divergent DNA methylation programs associated with transcript accumulation in sesame. Plant Sci. 2018;277:207–17.

    Article  CAS  PubMed  Google Scholar 

  23. Abid G, Mingeot D, Muhovski Y, Mergeai G, Aouida M, Abdelkarim S, et al. Analysis of DNA methylation patterns associated with drought stress response in faba bean ( Vicia faba L.) using methylation-sensitive amplification polymorphism (MSAP). Environ Exp Bot. 2017;142:34–44.

    Article  CAS  Google Scholar 

  24. Al-Lawati A, Al-Bahry S, Victor R, Al-Lawati AH, Yaish MW. Salt stress alters DNA methylation levels in alfalfa (Medicago spp). Genet Mol Res. 2016;15(1):15018299.

    Article  CAS  PubMed  Google Scholar 

  25. Angers B, Castonguay E, Massicotte R. Environmentally induced phenotypes and DNA methylation: how to deal with unpredictable conditions until the next generation and after. Mol Ecol. 2010;19(7):1283–95.

    Article  CAS  PubMed  Google Scholar 

  26. Mirouze M, Paszkowski J. Epigenetic contribution to stress adaptation in plants. Curr Opin Plant Biol. 2011;14(3):267–74.

    Article  CAS  PubMed  Google Scholar 

  27. Bräutigam K, Vining KJ, Lafon-Placette C, Fossdal CG, Mirouze M, Marcos JG, et al. Epigenetic regulation of adaptive responses of forest tree species to the environment. Ecol Evol. 2013;3(2):399–415.

    Article  PubMed  PubMed Central  Google Scholar 

  28. Li R, Hu F, Li B, Zhang Y, Chen M, Fan T, et al. Whole genome bisulfite sequencing methylome analysis of mulberry (Morus alba) reveals epigenome modifications in response to drought stress. Sci Rep. 2020;10(1):8013.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  29. Xu J, Zhou S, Gong X, Song Y, van Nocker S, Ma F, et al. Single-base methylome analysis reveals dynamic epigenomic differences associated with water deficit in apple. Plant Biotechnol J. 2018;16(2):672–87.

    Article  CAS  PubMed  Google Scholar 

  30. Day L. Proteins from land plants–potential resources for human nutrition and food security. Trends Food Sci Technol. 2013;32(1):25–42.

    Article  CAS  Google Scholar 

  31. Kang YJ, Kim SK, Kim MY, Lestari P, Kim KH, Ha B-K, et al. Genome sequence of mungbean and insights into evolution within Vigna species. Nat Commun. 2014;5(1):1–9.

    Article  Google Scholar 

  32. Kim SK, Nair RM, Lee J, Lee SH. Genomic resources in mungbean for future breeding programs. Front Plant Sci. 2015;6:626.

    Article  PubMed  PubMed Central  Google Scholar 

  33. Kumar S, Ayachit G, Sahoo L. Screening of mungbean for drought tolerance and transcriptome profiling between drought-tolerant and susceptible genotype in response to drought stress. Plant Physiol Biochem. 2020;157:229–38.

    Article  CAS  PubMed  Google Scholar 

  34. Fuller DQ, Harvey EL. The archaeobotany of Indian pulses: identification, processing and evidence for cultivation. Environ Archaeol. 2006;11(2):219–46.

    Article  Google Scholar 

  35. Hawkes J. The importance of wild germplasm in plant breeding. Euphytica. 1977;26(3):615–21.

    Article  Google Scholar 

  36. Tanksley SD, McCouch SR. Seed banks and molecular maps: unlocking genetic potential from the wild. Science. 1997;277(5329):1063–6.

    Article  CAS  PubMed  Google Scholar 

  37. Ning W, Zhai H, Yu J, Liang S, Yang X, Xing X, et al. Overexpression of Glycine soja WRKY20 enhances drought tolerance and improves plant yields under drought stress in transgenic soybean. Mol Breed. 2017;37(2):19.

    Article  Google Scholar 

  38. James A, Lawn R, Cooper M. Genotypic variation for drought stress response traits in soybean. II. Inter-relations between epidermal conductance, osmotic potential, relative water content, and plant survival. Aust J Agric Res. 2008;59(7):670–8.

    Article  Google Scholar 

  39. Chen H-M, Ku H-M, Schafleitner R, Bains TS, Kuo CG, Liu C-A, et al. The major quantitative trait locus for mungbean yellow mosaic Indian virus resistance is tightly linked in repulsion phase to the major bruchid resistance locus in a cross between mungbean [Vigna radiata (L.) Wilczek] and its wild relative Vigna radiata ssp. sublobata. Euphytica. 2013;192(2):205–16.

    Article  Google Scholar 

  40. Kanehisa M, Araki M, Goto S, Hattori M, Hirakawa M, Itoh M, et al. KEGG for linking genomes to life and the environment. Nucleic Acids Res. 2007;36:D480–4.

    Article  PubMed  PubMed Central  Google Scholar 

  41. Ghodke P, Khandagale K, Thangasamy A, Kulkarni A, Narwade N, Shirsat D, et al. Comparative transcriptome analyses in contrasting onion (Allium cepa L.) genotypes for drought stress. PLoS One. 2020;15(8):e0237457.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  42. Rauluseviciute I, Drablos F, Rye MB. DNA hypermethylation associated with upregulated gene expression in prostate cancer demonstrates the diversity of epigenetic regulation. BMC Med Genet. 2020;13(1):6.

    Google Scholar 

  43. Wang X, Duan CG, Tang K, Wang B, Zhang H, Lei M, et al. RNA-binding protein regulates plant DNA methylation by controlling mRNA processing at the intronic heterochromatin-containing gene IBM1. Proc Natl Acad Sci U S A. 2013;110(38):15467–72.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  44. Min X, Lin X, Ndayambaza B, Wang Y, Liu W. Coordinated mechanisms of leaves and roots in response to drought stress underlying full-length transcriptome profiling in Vicia sativa L. BMC Plant Biol. 2020;20(1):165.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  45. Basu P, Sharma A, Garg I, Sukumaran N. Tuber sink modifies photosynthetic response in potato under water stress. Environ Exp Bot. 1999;42(1):25–39.

    Article  Google Scholar 

  46. Zhu Y, Liu Q, Xu W, Zhang J, Wang X, Nie G, et al. De novo assembly and discovery of genes that involved in drought tolerance in the common vetch. Int J Mol Sci. 2019;20(2):328.

    Article  PubMed Central  Google Scholar 

  47. Kang YJ, Bae A, Shim S, Lee T, Lee J, Satyawan D, et al. Genome-wide DNA methylation profile in mungbean. Sci Rep. 2017;7:40503.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  48. Lu X, Wang X, Chen X, Shu N, Wang J, Wang D, et al. Single-base resolution methylomes of upland cotton (Gossypium hirsutum L.) reveal epigenome modifications in response to drought stress. BMC Genomics. 2017;18(1):297.

    Article  PubMed  PubMed Central  Google Scholar 

  49. Ng DW, Miller M, Yu HH, Huang TY, Kim ED, Lu J, et al. A role for CHH methylation in the parent-of-origin effect on altered circadian rhythms and biomass Heterosis in Arabidopsis intraspecific hybrids. Plant Cell. 2014;26(6):2430–40.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  50. Zilberman D, Gehring M, Tran RK, Ballinger T, Henikoff S. Genome-wide analysis of Arabidopsis thaliana DNA methylation uncovers an interdependence between methylation and transcription. Nat Genet. 2007;39(1):61–9.

    Article  CAS  PubMed  Google Scholar 

  51. Li X, Zhu J, Hu F, Ge S, Ye M, Xiang H, et al. Single-base resolution maps of cultivated and wild rice methylomes and regulatory roles of DNA methylation in plant gene expression. BMC Genomics. 2012;13(1):1–15.

    Article  Google Scholar 

  52. Zhang X, Yazaki J, Sundaresan A, Cokus S, Chan SW, Chen H, et al. Genome-wide high-resolution mapping and functional analysis of DNA methylation in arabidopsis. Cell. 2006;126(6):1189–201.

    Article  CAS  PubMed  Google Scholar 

  53. Su Y, Bai X, Yang W, Wang W, Chen Z, Ma J, et al. Single-base-resolution methylomes of Populus euphratica reveal the association between DNA methylation and salt stress. Tree Genet Genomes. 2018;14(6):1-11.

  54. Song Y, Tian M, Ci D, Zhang D. Methylation of microRNA genes regulates gene expression in bisexual flower development in andromonoecious poplar. J Exp Bot. 2015;66(7):1891–905.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  55. Kawakatsu T, Nery JR, Castanon R, Ecker JR. Dynamic DNA methylation reconfiguration during seed development and germination. Genome Biol. 2017;18(1):171.

    Article  PubMed  PubMed Central  Google Scholar 

  56. An YC, Goettel W, Han Q, Bartels A, Liu Z, Xiao W. Dynamic changes of genome-wide DNA methylation during soybean seed development. Sci Rep. 2017;7(1):12263.

    Article  PubMed  PubMed Central  Google Scholar 

  57. Bouyer D, Kramdi A, Kassam M, Heese M, Schnittger A, Roudier F, et al. DNA methylation dynamics during early plant life. Genome Biol. 2017;18(1):179.

    Article  PubMed  PubMed Central  Google Scholar 

  58. Rajkumar MS, Shankar R, Garg R, Jain M. Bisulphite sequencing reveals dynamic DNA methylation under desiccation and salinity stresses in rice cultivars. Genomics. 2020;112(5):3537–48.

    Article  CAS  PubMed  Google Scholar 

  59. Zhao Q, Du Y, Wang H, Rogers HJ, Yu C, Liu W, et al. 5-Azacytidine promotes shoot regeneration during Agrobacterium-mediated soybean transformation. Plant Physiol Biochem. 2019;141:40–50.

    Article  CAS  PubMed  Google Scholar 

  60. Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30(15):2114–20.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  61. Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods. 2015;12(4):357–60.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  62. Anders S, Pyl PT, Huber W. HTSeq--a Python framework to work with high-throughput sequencing data. Bioinformatics. 2015;31(2):166–9.

    Article  CAS  PubMed  Google Scholar 

  63. Anders S, Huber W. Differential expression of RNA-Seq data at the gene level–the DESeq package. Heidelberg: European Molecular Biology Laboratory (EMBL); 2012. p. 10. f1000research

    Google Scholar 

  64. Murray M, Thompson WF. Rapid isolation of high molecular weight plant DNA. Nucleic Acids Res. 1980;8(19):4321–6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  65. Krueger F, Andrews SR. Bismark: a flexible aligner and methylation caller for bisulfite-Seq applications. Bioinformatics. 2011;27(11):1571–2.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  66. Feng H, Conneely KN, Wu H. A Bayesian hierarchical model to detect differentially methylated loci from single nucleotide resolution sequencing data. Nucleic Acids Res. 2014;42(8):e69.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  67. Young MD, Wakefield MJ, Smyth GK, Oshlack A. Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biol. 2010;11(2):1–12.

    Article  Google Scholar 

  68. Mao X, Cai T, Olyarchuk JG, Wei L. Automated genome annotation and pathway identification using the KEGG Orthology (KO) as a controlled vocabulary. Bioinformatics. 2005;21(19):3787–93.

    Article  CAS  PubMed  Google Scholar 

Download references

Acknowledgements

Not applicable.

Funding

This study was financially supported by National Natural Science Foundation of China (No. 31801431), the Natural Science Foundation of Shandong (ZR2018LC023), and the Qingdao Agricultural University Scientific Research Foundation (663/1120075).

Author information

Author notes

  1. Peilei Zhao and Bao Ma contributed equally to this work.

Authors and Affiliations

  1. College of Life Sciences, Qingdao Agricultural University, Qingdao, China

    Peilei Zhao, Bao Ma, Chunmei Cai & Jihua Xu

  2. College of Horticulture, Qingdao Agricultural University, Qingdao, China

    Peilei Zhao

Authors
  1. Peilei Zhao
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Bao Ma
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Chunmei Cai
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Jihua Xu
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

JX wrote the manuscript. JX designed all the experiments. PZ, BM, CC, and JX performed all the experiments and analyzed the data. All the authors have approved the final manuscript.

Corresponding author

Correspondence to Jihua Xu.

Ethics declarations

Ethics approval and consent to participate

The study complied with relevant institutional, national, and international guidelines and legislation.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Figure S1.

Heatmap of DEGs in four pairwise comparisons. a C70 vs C61. b D70 vs D61. c D70 vs C70. d D61 vs C61.

Additional file 2: Figure S2.

Validation of the reliability of RNA-seq data by qRT-PCR. The vertical axis indicates the fold change when drought stressed D61 compared with control C61 (a), and D70 compared with C70 (b); the horizontal axis shows the eight DEGs selected.

Additional file 3: Figure S3.

Number of differentially methylated regions in D61 vs C61 and D70 vs C70.

Additional file 4: Figure S4.

DNA methylation levels of DMRs in all CG, CHG, and CHH contexts displayed by violin boxplots in D61 vs C61 (a) and D70 vs C70 (b). Number of DMRs in different regions of the genome in D61 vs C61 (c) and D70 vs C70 (d).

Additional file 5: Figure S5.

Relationship between DNA methylation and gene expression in C61, D61, C70 and D70 in gene body. Expression profiles of different methylated levels at CG (a), CHG (b) and CHH (c) were investigated. The gene body methylation levels were classified into five groups with group.1st the lowest and group.5th the highest.

Additional file 6: Table S1.

Primers used for qRT-PCR analysis.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Zhao, P., Ma, B., Cai, C. et al. Transcriptome and methylome changes in two contrasting mungbean genotypes in response to drought stress. BMC Genomics 23, 80 (2022). https://doi.org/10.1186/s12864-022-08315-z

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-022-08315-z

Keywords

BMC Genomics

ISSN: 1471-2164

Contact us