Transcriptomics Analysis Reveals a More Refined Regulation Mechanism of Methylation in a Drought-Tolerant Variety of Potato

Whether DNA methylation modification affects the gene transcription and expression of potatoes under drought stress is still unknown. In this study, we used comparative transcriptomics to explore the expression pattern of related genes of the drought-tolerant variety Qingshu 9 (Q) and the drought-sensitive variety Atlantic (A) under drought stress and DNA methylation inhibitor treatment. The results showed that there was a significant difference in the number of DEGs between the two varieties’ responses to mannitol and 5-azad C, especially when they were co-treated with two reagents, and the gene expression of Q was more sensitive to mannitol after two hours. Furthermore, we found that these differentially expressed genes (DEGs) were significantly enriched in DNA replication, transcription, translation, carbohydrate metabolism, photosynthesis, signal transduction, and glutathione metabolism. These results indicate that the difference in the background of methylation leads to the difference in drought resistance of the two varieties. The complexity of the DNA methylation of variety Q might be higher than that of variety A, and the method of methylation regulation is more refined. This study systematically expands the understanding of the molecular mechanism wherein DNA methylation regulates the response to drought stress.


Introduction
The response of plants to drought stress is a complex biological process that is jointly regulated by a variety of signaling pathways [1]. When suffering from drought stress, plants will activate the downstream functional gene expression through signal transduction, including kinases, hormones, calcium ions, transcription factors, etc. [2,3]. In recent years, more and more studies have shown that epigenetic regulation plays an important role in the response of plants to drought stress [4,5]. DNA methylation is a kind of epigenetic modification, and plant DNA methylation functions in biological processes such as chromatin inhibition, cell differentiation, gene expression regulation, and plant growth and development [6]. In addition, abiotic stress can cause changes in the degree of DNA methylation during plant development, and in turn, DNA methylation is also involved in plants' response to abiotic stress [7,8]. Therefore, it is particularly important to clarify the impact of drought stress on the changes in DNA methylation.
After drought stress, HUW468 showed hypermethylation in the whole genome, while C306 showed hypomethylation in the whole genome [9]. In Populus tomentosa, the proportion of cytosine methylation under drought stress is 10.04% and only 7.75% under watering treatment [10]. The same pattern was noted in drought-sensitive rice varieties, and the DNA methylation level increased under drought [11]. Under drought stress, methylation 100 bp upstream of the transcription start site (TSS) of Populus tomentosa inhibited gene expression [10]. The above results indicate that transcription factors may play an important role in the drought stress response of P. tomentosa through changes in DNA methylation.
The different tolerances to salt and drought among the three rice varieties are related to the existence of differentially methylated regions [12]. Barley has a high and stable level of whole-genome DNA methylation under water-shortage conditions, but there is a demethylation trend in the leaves and a methylation trend in the roots [13]. During the drought-recovery process, DNA methylation and demethylation occurred simultaneously in the leaves of Boea hygrometirica but were dominated by demethylation [14]. In perennial ryegrass, the methylation level was 57.67% under normal conditions and decreased to 47.39% under drought stress [15]. As for Macrotyloma uniflorum, the methylation degree of drought-sensitive HPKC2 (10.1%) was higher than that of drought-tolerant HPK4 (8.6%) [16]. A comparative analysis of the differential methylation sites showed that these genes are mainly involved in plant signal transduction, material and energy metabolism, plant growth and development, stress response, and other physiological reaction processes [17].
Sulfite sequencing analysis showed that compared with the control, the DNA methylation level of the AtGSTF14 promoter region was significantly reduced, and its expression levels increased significantly in the drought-treated Arabidopsis plants [18]. Under drought stress treatment, the level of genomic DNA methylation of pea was significantly increased, thereby inhibiting gene expression from adapting to drought [19]. There is a difference in genome-wide DNA methylation between drought-tolerant (Bachar) and drought-sensitive (F177) Vicia faba genotypes. However, drought stress leads to a decrease in Bachar and F177 genomic DNA methylation levels, while the Bachar demethylation level is higher than that of F177 [20]. Furthermore, six differentially methylated regions related to drought stress were identified, including genes encoding lipoxygenase (LOX), calcium-dependent protein kinase (CDPK), the ABC transporter family (ABC), glycosyl hydrolase (GH), and phosphoenolpyruvate carboxylase (PEPC). Under drought stress, the expression levels of VfLOX, VfCDPK, VfABC, and VfGH genes in Bachar are higher than that of F177. More importantly, the DNA methylation changes caused by drought resistance can be inherited by the offspring, that is, stress memory, which has been confirmed in Arabidopsis [21], rape [22], and rice [23]. These indicate that the DNA methylation of plants is heritable and may contribute to improving their drought resistance, which may have significance in plant breeding.
5-Aza-2 -deoxycytidine (5-azad C) is commonly used in plants to inhibit the degree of DNA methylation [24]. 5-azad C is an analog of cytosine and binds irreversibly to DNA methyltransferase, which makes it difficult for the genome to maintain the methylation status [25,26]. It was found that the DNA demethylation treatment seriously affected the growth and development of the potato variety of Atlantic test-tube seedlings [27]. In wheat, 5-azacytidine can reduce the level of wheat genomic DNA methylation, thereby improving salt tolerance [28]. Treatment of cabbage seedlings and rice seeds showed abnormal phenomena, such as dwarf plants and small leaves after growth and maturity, and the degree of genomic DNA methylation was measured and showed a significant reduction [29]. Therefore, the function of DNA methylation is complex in plants.
The lack of water in potatoes (Solanum tuberosum) seriously affected the normal growth and development of potato plants [30]. At present, research on potato DNA methylation mainly focuses on the changes in DNA methylation under in vitro culture and stress environment [31][32][33]. In our previous work, we found that DNA demethylation under drought stress is involved in regulating the phenotypic traits and physiological and biochemical responses of potatoes [34]. In this study, comparative transcriptomics analysis was used to explore the expression characteristics of related genes in the drought-tolerant variety Qingshu 9 (Q) and drought-sensitive variety Atlantic (A) under drought stress and DNA methylation inhibitors. We hope to clarify the molecular mechanism by which DNA methylation participates in the regulation of potato response to drought stress, and promote potato drought resistance molecular breeding.

Plant Cultivation and Treatment
The drought-tolerant potato variety used was Q, and the drought-sensitive variety was A. We can use the paper bridge method to inoculate the stems of two varieties of potato tissue cultured seedlings in MS liquid medium (without agar); each stem has an axillary bud, and the test-tube seedling culture environment condition is 22 (±2) • C. The light intensity is 25-37.5 µmol m −2 s −1 and the light duration is 16 h every day. After 24 days of cultivation, the potato test-tube plantlets were treated with three treatment methods; namely, 200 mmolL-1 mannitol treatment (simulated drought treatment and referred to as M in the following), 60 µmol L-1 5-azadC treatment (demethylation and referred to as A in the following), and the two reagents were treated simultaneously (referred to as D in the following). After the treatment, the stems and leaves of the test tube seedlings were treated for 2 h, 6 h, 12 h, and untreated (control, 0 h), quickly frozen with liquid nitrogen, and stored in an ultra-low temperature refrigerator at −80 • C for later use. Three biological replicates were used for each treatment. 5-azadC (5-aza-2 -deoxycytidine) was purchased from Sigma-Aldrich.

RNA Extraction, Library Construction, and Transcriptome Sequencing
First, the total RNA of 60 samples (2 varieties × (control + (3 treatments × 3 time points)) × 3 biological replicates) were extracted using Trizol reagent (Invitgen, Waltham, MA, USA). Then, total RNA was treated using Turbo DNase I (Ambion, Austin, TX, USA) for 30 min and purified using RNeasy ® Plant Mini Kit (QIAGEN, Hilden, Germany). Then, the TruSeq RNA sample Prep V2 kit (Illumina, San Diego, CA, USA) was used to construct a transcriptome-sequencing cDNA library and the Agilent 2100 TapeStation (Agilent, Santa Clara, CA, USA) system was used to detect the quality of the cDNA library. Finally, pairedend sequencing was performed on the Illumina NovaSeq 6000 sequencer by Beijing Biomics Biotech Co., Ltd. (Beijing, China). After the sequencing was completed, we uploaded the data to the NCBI database under the BioProject number PRJNA783478.

Transcriptome Data Analysis
First, we used FASTQC (https://www.bioinformatics.babraham.ac.uk/projects/fas tqc/, accessed on 6 September 2021) to perform quality control on the off-machine data and calculate the Q20, Q30, and GC content. The potato reference genome sequence was downloaded from the ENSEMBL plant database (http://solanaceae.plantbiology.msu.edu/, accessed on 6 September 2021). We used TopHat2 [35] to compare Clean Reads with the reference genome to obtain the position information on the reference genome or gene, as well as the unique sequence feature information of the sequenced sample. We used the Cuffquant and Cuffnorm components of the Cufflinks (https://github.com/cole-trapne ll-lab/cufflinks, accessed on 6 September 2021) software to quantify the transcript and gene-expression level through the location information of Mapped Reads on the gene.
Then, we used HISAT2 (http://daehwankimlab.github.io/hisat2/, accessed on 6 September 2021) to compare sequencing reads to the potato reference genome to obtain sam files. We used Samtools (https://github.com/samtools/samtools, accessed on 6 September 2021) to convert sam files into bam files and sort them, used Cufflinks to assemble and quantify the sorted bam files, used HTseq-count to calculate the counts of each sample transcript, used Python script to calculate the length of each gene in the reference genome, and used the R package Dseq2 to calculate the FPKM (fragments per kilobase of transcript per million fragments mapped) of each gene and analyze the differential expression between samples. Corrected p < 0.01 and a multiple of change greater than 2 were defined as indicating differentially expressed genes (DEGs). Pearson's correlation coefficient was used as the evaluation index of biological repetition correlation, and the correlation between samples was statistically mapped [36]. The closer r was to 1, the stronger the correlation between the two repetitive samples.
Next, we used the online tool:Profiler (https://biit.cs.ut.ee/gprofiler_beta/gost, accessed on 12 October 2021) to perform COG, GO, and KEGG functional annotation and enrichment analysis of differentially expressed genes, and p < 0.05 as the significant enrichment GO term. Additionally, we used KEGG pathway screening conditions. We used the interaction relationship in the STRING protein interaction database (http://string-db.org/, accessed on 12 October 2021) to analyze the differential gene protein interaction network.
Based on the TopHat2 (http://ccb.jhu.edu/software/tophat/index.shtml, accessed on 12 October 2021) alignment results of the reads of each sample and the reference genome sequence, GATK [37] software was used to identify single-base mismatches between the sequenced samples and the reference genome and potential SNP sites and InDel (insertiondeletion).

Promoter Methylation CpG Island Prediction and Cis-Elements Analysis of Important DEGs
We first use the perl script to extract the 1500 bp promoter sequence upstream of the analyzed gene in batches and then the online software MethPrimer (http://www.urog ene.org/cgi-bin/methprimer/MethPrimer.cgi, accessed on 12 October 2021) to regulate the differentially enriched gene. 5 CpG island prediction was performed in the area, and the prediction parameters were set as: CpG island length > 100 bp, CPG detection content/CpG expected content (Obs/Exp) > 0.60, and GC content > 50%. We used Plant Care (http://bioinformatics.psb.ug ent.be/webtools/plantcare/html/, accessed on 12 October 2021) to predict the cis-acting elements of gene promoters.

Real-Time qPCR Verification
Primer Express 3.0.1 was used to design primers for the 20 selected DEGs (Supplementary Table S1), and then we performed real-time fluorescent quantitative PCR analysis. The RNA of three biological replicates of the selected sample (after 6 h of three kinds of treatment) was reverse-transcribed into cDNA using a cDNA synthesis kit (TaKaRa, Kyoto, Japan), actin was used as an internal reference, and qPCR was performed using the QuantStudio5 real-time fluorescent quantitative PCR system (ABI, California, USA)-PCR detection. Three biological replicates of each gene were used in each sample, and the relative expression level was calculated according to the 2 −∆∆CT method.

Transcriptome Sequencing and Functional Annotation and Identification of the DEGs
As shown in the Supplementary Materials, Table S2, after data filtering and quality evaluation, 60 samples produced approximately 410 G of clean data, and each sample produced approximately 6.8 G of effective data on average. The GC content of each sample was in the range of 42.15~47.82%, the Q20 value was above 97%, and the Q30 value was above 91.98%. At the same time, the unique comparison rate between the reads of each sample and the reference genome sequence was over 98%.  A Venn diagram was used to count the common DEGs and specific DEGs under different treatment conditions. After 2 h of treatment (Figure 2a), there were 583 DEGs shared by the three kinds of treatment, 1585 unique DEGs for methylation inhibitors, 109 unique DEGs for simulated drought, and 3740 unique DEGs for dual treatments. After 6 h of treatment, there were 641 DEGs shared by three kinds of treatment, 486 unique DEGs for methylation inhibitors, 1339 unique DEGs for drought treatment, and 2588 DEGs for double dry treatment. After 12 h of treatment, there were 1184 DEGs shared by three kinds of treatment, 402 unique DEGs for methylation inhibitors, 858 unique DEGs for simulated drought treatment, and 3648 DEGs for dual treatments. Compared with treatment for 6 h and 2 h, there were 35 DEGs shared by three kinds of treatment, 192 specific DEGs for methylation inhibitors, 215 specific DEGs for simulated drought treatment, and 3620 DEGs for dual treatments. Compared with treatments for 12 h and 6 h, there were only eight DEGs shared by three kinds of treatment; there were no specific DEGs for the methylation inhibitor and dual treatment, and the number of specific DEGs of the simulated drought treatment was 12. Compared with the simulated drought treatment (Figure 1d), under the dual treatment condition, after 2 h, there were 3983 DEGs, of which 1637 were upregulated and 2346 were downregulated. After 6 h, there were 6951 DEGs, of which 3344 were upregulated and 3607 were downregulated. After 12 h, there were 5445 DEGs, of which 2316 were upregulated and 3129 were downregulated.
A Venn diagram counted the common and specific DEGs under different treatment conditions (Figure 2b). After 2 h, there were 659 DEGs shared by three kinds of treatment, 1759 DEGs were specific for treatment with methylation inhibitors, 3218 DEGs were specific for dual treatments, and 58 DEGs were specific for simulated drought treatment. After 6 h, there were 233 DEGs shared by three kinds of treatment, 1017 unique DEGs for treatment with methylation inhibitors, and 848 unique DEGs for treatment with dual treatments. There were 2284 DEGs unique to simulated drought treatment. After 12 h, there were 858 DEGs shared by three kinds of treatment, 943 DEGs specific for treatment with methylation inhibitors, and 3731 DEGs specific for dual treatments.
There were 405 DEGs unique to simulated drought treatment. Compared with treatment for 6 h and 2 h, there were 12 DEGs shared by three kinds of treatment, 2142 specific DEGs treated with methylation inhibitors, and 989 DEGs treated with dual treatments. There were 145 DEGs unique to simulated drought treatment. Compared with treatment for 12 h and 6 h, there were 5 DEGs shared by three kinds of treatment, 522 specific DEGs treated with methylation inhibitors, and 1944 DEGs treated with dual treatments. There were 66 DEGs unique to simulated drought treatment. Under the double-treatment condition (Figure 1e), after 2 h, there were 4668 DEGs in common in the two varieties, of which 2457 were upregulated and 2211 were downregulated. After 6 h of treatment, there were 4363 DEGs in common in the two varieties, of which 2525 were upregulated and 1838 were downregulated. After 12 h of treatment, there were 2950 DEGs in common in the two varieties, of which 1500 were upregulated and 1450 were downregulated.
A Venn diagram was used to count the differentially expressed genes and specific differentially expressed genes under different treatment conditions. After 2 h of treatment (Figure 2c), there were 432 DEGs shared by three kinds of treatment, 994 specific DEGs treated with methylation inhibitors, and 3626 DEGs treated with dual treatments. There were 172 DEGs unique to simulated drought treatment. After 6 h of treatment, there were 442 DEGs shared by three kinds of treatment, 595 unique DEGs treated with methylation inhibitors, and 3288 DEGs treated with dual treatments. There were 546 DEGs unique to simulated drought treatment. After 12 h of treatment, there were 329 DEGs shared by three kinds of treatment, 756 unique DEGs treated with methylation inhibitors, and 2325 DEGs treated with dual treatments. There were 189 DEGs unique to simulated drought treatment.

Venn Diagram Statistics of Each Treatment's DEGs in the Two Varieties
Under the methylation inhibitor treatment condition ( Figure 3a) for 2 h, there were 1660 DEGs in common between the two varieties, 1836 DEGs specific to species A, and 2491 DEGs specific to species Q. After 6 h, there were 872 DEGs in common between the two varieties, 1545 DEGs specific to variety A, and 1515 DEGs specific to variety Q. After 12 h of treatment, there were 1555 DEGs shared between the two varieties, 1081 unique DEGs of variety A, and 1641 unique DEGs of variety Q. Compared with the treatment for 6 h and 2 h, there were 184 DEGs in common between the two varieties, 239 DEGs specific to variety A, and 2434 DEGs specific to variety Q. Compared with the treatment for 12 h and 6 h, there were 16 DEGs in common between the two varieties, 50 DEGs unique to variety A, and 765 unique DEGs to variety Q.

Protein-Protein Interaction Network Analysis
We constructed an interaction network based on the DEGs retrieved from known PPI databases (Figures 7-9, Supplementary Material, Figures S12-S17). In terms of comparisons between variety A and Q under control conditions, 689 DEGs can be linked via known relationships, including reaction, catalysis, binding, and inhibition (Supplementary Materials, Figure S12). Significantly, the Soltu.DM.08G017320 (17 links and was mainly upregulated), Soltu.DM.07G003860 (14 links and was mainly downregulated), and Soltu.DM.10G024790 (13 links and was mainly upregulated) encoding aldehyde dehydrogenase, ribosomal protein L30/L7 family protein, and ribosome production factor 2 were proposed as potential hubs.

Analysis of Methylation Sites in the Promoter of Key DEGs
CpG island prediction and cis-element analysis were performed on the promoter of DEGs enriched by KEGG pathways related to plant drought resistance. The results ( Figure 10) showed that Soltu.DM.09G006730 related to the plant glutathione metabolism pathway was predicted, and its CpG island was located between −736and −866 in the promoter sequence of this gene, with a size of 131 bp ( Figure 10). Next, we used the online database Plant Care to analyze the cis action element of the gene promoter sequence. As shown in Figure 10, there were ABRE, CAAT-box, and A-box elements in the CpG island region of the Soltu.DM.09G006730 gene, of which ABRE and CAAT -box are related to the response of plants to drought stress. The above results indicate that when potatoes respond to drought stress, the ABRE and CAAT-box elements in the promoter region of the glutathione transferase gene activate the expression of the gene through DNA demethylation to cope with stress.

Verification of qPCR
Twenty genes were randomly selected from DEGs, and qPCR was used to detect the relative expression of these genes in Atlantic and Qingshu 9 under different treatments (Supplementary Materials, Figure S18). The results show that the qPCR results and the transcriptome sequencing results have basically the same change trends. This demonstrates the reliability of the transcriptome sequencing data and further validates the response of these genes to drought and DNA methylation inhibitors.

Discussion
Various epigenetic modifications before transcription, such as DNA methylation, histone, and chromatin remodeling modifications, also play an important role in regulating plants' response to drought stress [4,5]. Under drought-stress conditions, plants start or stop a series of cascade reactions (response pathways) caused by genes related to stress response through methylation and demethylation of genomic DNA and increase the activity of enzymes in response to stress, thereby resisting drought, as stress is known to impact plant growth and development [35]. Therefore, understanding the dynamic changes in plant epigenetic DNA methylation modification under drought stress will help to study functional gene expression and further clarify the molecular mechanisms of plants' responses to drought. This result also fully implies epigenetic modifications such as DNA methylation. The potential importance and application value of our research are in the fields of plant drought resistance and breeding. Epigenetic modification with DNA methylation as the core is involved in the balance of gene expression, genome stability, and adaptability, which is very important for improving crop traits [38]. At present, there is little research on the relationship between potato drought resistance and DNA methylation. In a previous study, the responses of different potato cultivars to mannitol and 5-azadC were similar. With treatments of mannitol or 5-azadC, the dry/fresh weight, shoot height, leaf number, and chlorophyll content decreased significantly, while the activities of SOD, POD, CAT, and contents of proline and MDA increased significantly. However, the number of branches, root length, and average root thickness did not change significantly, indicating that the regulatory pathways may be different [39].
Analysis of DNA methylation polymorphism in response to drought stress in rice has shown that drought-sensitive genotypes are more likely to have hypermethylation than drought-tolerant genotypes [11]. In this study, a large batch of RNA-seq was used to further study the effect of methylation between different varieties on the drought stress response. This study found that under untreated conditions, there were 1295 DEGs between the two varieties ( Figure 1e). Through functional annotation analysis, we found that these differential genes mainly belong to plant hormone signal transduction, defense response, plant-pathogen interaction, protein processing in the endoplasmic reticulum, and starch and sucrose metabolism (Supplementary Materials, Figure S5).
In both drought-tolerant variety Q and drought-sensitive variety A, many genes were upregulated or downregulated after 6 h of drought stress ( Figure 1B), which lasted for up to 12 h. Compared with drought-sensitive variety A, the number of DEGs in drought-tolerant variety Q decreased sharply after 12 h of drought treatment, which indicated that variety Q could quickly adapt to drought stress by regulating gene expression. At each treatment time point under the dual-treatment condition (Figure 1c), the number of DEGs in variety Q was lower than that in variety A, and the number of highly expressed genes in variety Q was lower than that in variety A regardless of drought or dual-treatment conditions. A (Figure 1d). These results also indicated that variety Q only needs to change relatively few gene expressions to cope with drought stress.
After 2 h of methylation inhibitor treatment, many genes began to be upregulated or downregulated (Figure 1a). In drought-sensitive variety A, the difference between treatment for 6 h and treatment for 2 h was not significant, and the difference between treatment for 6 h and treatment for 12 h was not obvious, indicating that treatment for 2 h can inhibit or induce gene expression. In drought-tolerant variety Q, the difference between treatment for 6 h and treatment for 2 h was larger, and there were 783 DEGs between treatment for 6 h and treatment for 12 h, indicating that the methylation inhibitor could continue to act on variety Q. These results indicate that, to some extent, the background methylation complexity of the genomic DNA of variety Q may be higher than that of variety A.
We can also infer the same result from the Venn diagram (Supplementary Materials, Figures S3 and S4). In cultivar A (Figure 2a), the number of specific DEGs under dual treatment conditions is always much greater than that of methylation inhibitors and drought treatments. For variety Q, the number of DEGs unique to the dual treatment dropped sharply when the drought treatment was 6 h. Regardless of variety A or variety Q, the number of DEGs shared by drought treatment and methylation inhibitors was far more than that of unique DEGs after 2 h of treatment, which indicated that potato genome DNA methylation was deeply involved in the drought stress response. In addition, under the conditions of methylation inhibition and drought treatment for 2 h (Figure 3a,b), the number of DEGs unique to drought-tolerant variety Q is more than that of drought-sensitive variety A, which also indicates that the level of DNA methylation in variety Q is higher, and the reaction is more sensitive and faster. However, under the dual-treatment conditions, the number of DEGs specific to drought-tolerant variety Q was slightly less than that of drought-sensitive variety A, which indicated that there were certain differences in the regulation of methylation levels between the two varieties.
According to the functional annotation enrichment analysis of DEGs, in droughtsensitive variety A (Supplementary Materials, Figure S6), under the conditions of methylation inhibitor treatment, light stress response, photosynthesis, and other pathways were first affected (2 h), followed by transcription regulation and calcium signal transduction, etc. (6 h). Under drought conditions, DNA transcriptional regulation and other pathways first reacted in large numbers (2 h), and after 6 h of treatment, genes in the photosynthesis pathway began to be differentially expressed in large numbers. The situation is different in drought-tolerant variety Q (Supplementary Materials, Figure S7). Under methylation inhibitors and drought treatment, they react quickly (2 h) in the pathways of transcription regulation and calcium signal transduction, but at the same time, under treatment with methylation inhibitors, membrane system components are also regulated. As the treatment time increases, the metabolic pathways related to photosynthesis are gradually affected.
In our previous work, 5-azadC was used to treat different drought-sensitive potato varieties [34] and we found that DNA demethylation under drought stress is involved in regulating the phenotypic traits and physiological and biochemical responses of potatoes. Similar results were also observed by application of DNA methylation inhibitor 5azacytidine on watermelon (Citrullus lanatus), especially in those involved in the regulation of inositol-related genes [40], providing a preliminary insight into the altered methylation levels of watermelon cells in response to osmotic stress. Recent 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 downregulated DEGs than upregulated DEGs, which proved that the differentially methylated regions (DMRs) in the gene body were significantly negatively correlated with DEGs [41]. In this study, the functional annotation of DEGs between the two varieties ( Figure 4) also verified this phenomenon. Under treatment with methylation inhibitors, DEGs were first significantly enriched in photosynthesis, and under drought conditions, DEGs were first significantly enriched in defense reactions and signal transduction. Therefore, these results indicate that the response process of different varieties to drought stress is different, and this process is affected by genome methylation.
As we learned from the results of PPI analysis (Figures 7-9 and S12-S17), under methylation treatment, the key regulatory pathways for the differences between the two varieties focus on protein translation and ATP synthesis. During drought treatment, in addition to protein translation, the key regulatory pathways for differences between the two varieties focused on DNA replication and transcription and lignin synthesis. Under dualtreatment conditions, differences in DNA replication were more significant between the two species, and the MCM complex became a key regulatory node, followed by photosynthesis, protein translation, defense response, and secondary metabolism. These results indicate that the differences in the regulation of genomic DNA methylation levels between the two varieties are mainly based on the replication and transcription of genomic DNA, protein translation, photosynthesis, and secondary metabolism.

Conclusions
Our research shows that compared with drought-sensitive variety A, drought-tolerant variety Q is more sensitive and milder in response to drought stress. Potato genome DNA methylation is deeply involved in the drought stress response. The DNA methylation levels between the two varieties are significantly different in DNA replication and transcription, protein translation, photosynthesis, defense response, and secondary metabolism regula-tion. Variety Q's complexity of background methylation of genomic DNA may be higher than that of variety A, and methylation regulation is more refined. More importantly, we also identify some hub proteins that could be used as key candidate genes in molecular breeding.

Conflicts of Interest:
The authors declare no conflict of interest.