key: cord-0066264-toz6xd2w
authors: Jia, Xin; Feng, Hui; Bu, Yanhua; Ji, Naizhe; Lyu, Yingmin; Zhao, Shiwei
title: Comparative Transcriptome and Weighted Gene Co-expression Network Analysis Identify Key Transcription Factors of Rosa chinensis ‘Old Blush’ After Exposure to a Gradual Drought Stress Followed by Recovery
date: 2021-07-15
journal: Front Genet
DOI: 10.3389/fgene.2021.690264
sha: ec32afbe5e6ae06c28d78b6366dd40089018773d
doc_id: 66264
cord_uid: toz6xd2w

Rose is one of the most fundamental ornamental crops, but its yield and quality are highly limited by drought. The key transcription factors (TFs) and co-expression networks during rose’s response to drought stress and recovery after drought stress are still limited. In this study, the transcriptomes of leaves of 2-year-old cutting seedlings of Rosa chinensis ‘Old Blush’ from three continuous droughted stages (30, 60, 90 days after full watering) and rewatering were analyzed using RNA sequencing. Weighted gene co-expression network analysis (WGCNA) was used to construct a co-expression network, which was associated with the physiological traits of drought response to discovering the hub TFs involved in drought response. More than 45 million high-quality clean reads were generated from the sample and used for comparison with the rose reference genome. A total of 46433 differentially expressed genes (DEGs) were identified. Gene Ontology (GO) term enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis indicated that drought stress caused significant changes in signal transduction, plant hormones including ABA, auxin, brassinosteroid (BR), cytokinin, ethylene (ET), jasmonic acid (JA) and salicylic acid (SA), primary and secondary metabolism, and a certain degree of recovery after rewatering. Gene co-expression analysis identified 18 modules, in which four modules showed a high degree of correlation with physiological traits. In addition, 42 TFs including members of NACs, WRKYs, MYBs, AP2/ERFs, ARFs, and bHLHs with high connectivity in navajowhite1 and blue modules were screened. This study provides the transcriptome sequencing report of R. chinensis ‘Old Blush’ during drought stress and rewatering process. The study also identifies the response of candidate TFs to drought stress, providing guidelines for improving the drought tolerance of the rose through molecular breeding in the future.

Rose is one of the most fundamental ornamental crops. With an important economic, cultural, and symbolic value, the rose has a significant market share in the floriculture industry. Cut rose flowers account for 32% of the world's total cut flower trade. The rose occupies nearly a quarter of the Netherlands floriculture market (Mordor Intelligence, 2020). However, their growth conditions are closely related to changes in the external environment. Most roses suffer from various degrees of drought, which result in a sizable loss of productivity and quality around the world. As a kind of abiotic stress, drought is becoming an even more serious issue due to global warming. Changes in the earth's climate have an adverse impact on plant growth, crop production, and distribution (Zhu, 2002) . It is therefore pressing to identify the response mechanism in rose under the drought stress.

To resistance drought stress, plants have evolved several adaptive mechanisms depending on the activation of the molecular networks involved in stress perception, signal transcription and transduction, expression of specific stressrelated genes, and physiological reaction of tolerance (Chaves et al., 2003) . Damage to plants under drought stress is manifested in the peroxidation caused by the accumulation of active oxygen in the plant (Niu et al., 2013) . The role of various antioxidant protection enzymes in plants under drought stress has been extensively studied. In addition, malondialdehyde (MDA) is an important indicator of cell stability (Fang et al., 2018) . Studies have shown that drought stress leads to different trends in plant antioxidant systems (Zandalinas et al., 2017; Niu et al., 2018) . The activity of four antioxidant enzymes fluctuate while the MDA content continues to increase in Herbaceous peony with drought stress (Niu et al., 2018) . Drought stress accumulates amounts of osmotic adjustment substances in plants such as proline and soluble sugars. In research on the drought stress of rose, it was found that as the drought stress deepened the content of proline and soluble sugar increased (Gadzinowska et al., 2019; Adamipour et al., 2020; Al-Yasi et al., 2020) . It is crucial to identify the genes and biological pathways involved in drought stress tolerance and improve the understanding of the molecular mechanisms relating to drought resistance in plants (Valliyodan and Nguyen, 2006) . A large number of signaling metabolic pathways and numerous genes have been determined in several model plants exposed to drought environments, such as Arabidopsis (Huang et al., 2018) , rice , and so on. Transcription factors (TFs) control the regulating downstream stress response genes that play essential roles in various abiotic stress response processes (Joshi et al., 2016) . During signal transduction, TFs directly regulate the expression of related genes by acting as a molecular switch. These TFs interact specifically with cis-elements located in the promoter region of the regulator (Franco-Zorrilla et al., 2014) . As the plant drought tolerance network is extremely complex, this study started by exploring the transcriptional regulatory network of TFs for drought resistance. To date, RcLEA , RcMYBPA2 (Li et al., 2019c) , RhMYB96 (Jiang et al., 2018) , RcNAC3 , and RcNAC31 (Ding et al., 2019) genes have been reported to improve the abiotic stress resistance of the transgenic plant. However, information on the mechanism of rose response to drought stress is still limited.

As a deep sequencing technology, transcriptome technology has been successfully used to explore the gene network of nonmodel plant species in response to drought stress. Droughtresponsive transcriptome studies have been conducted in ornamental plants such as chrysanthemum (Xu et al., 2013) , tree peony (Paeonia section Moutan DC.) (Zhao et al., 2019a) , and Dianthus spiculifolius . Previously, only a suppression subtractive hybridization library was constructed for the dehydration stress of cut rose (Rosa hybrida cv. Samantha) petals, which contained 3513 unique expression sequence tags, and its expression profile during the dehydration cycle was analyzed (Dai et al., 2012) . Although the transcriptome of R. chinensis 'Mutabilis' in the three stages of drought stress has been reported (Li et al., 2021a) and provided some basis, drought stress generally happens continuously, and it is more essential to discover how the rose recovers after drought stress. In other words, knowledge about the gene co-expression network involved in rose drought and recovery is very limited. WGCNA has become an important tool to identify gene co-expression related to associating with function, which has been widely used in the research of horticultural plants (Lu et al., 2019; Liu et al., 2021) .

Rosa chinensis 'Old Blush' originated from China and has participated in modern rose hybrid breeding. R. chinensis 'Old Blush' is the diploid species in Rosa, which may represent a convenient genomic model for Rosa research and may be used as a model plant to study rose response to drought stress. Currently, the genome of R. chinensis 'Old Blush' was published, providing a valuable database for further function and genomic research in rose (Raymond et al., 2018) . To obtain unique insights into the molecular mechanisms at the transcriptome level in response to drought and recovery, this study compared transcriptome profiles of R. chinensis 'Old Blush' at continuous drought stages (30, 60, 90 days after full watering) and rewatering, and constructed co-expression modules using WGCNA. The Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis in differentially expressed genes (DEGs) were performed on the modules with significant physiological data and the key TFs in the module were identified. The study aims to obtain an understanding of the drought responsive regulatory network at the molecular level and identify the key TFs. Furthermore, the results of this study can provide the basis for a better understanding of not only other R. hybrida cultivars but also other woody plants responding to drought stress, and odder candidate genes for further rose resistance research.

Two-year-old healthy cutting seedlings of the R. chinensis 'Old Blush' used in the experiments were collected from the Beijing Institute of Landscape Architecture (Beijing, China). The plant material was planted in the pot with an inner size of 26 cm in diameter, 22 cm in height, and 18 cm in diameter at the bottom. It contained sterilized soil with a mixture of nutritive soil, garden soil, perlite, and vermiculite with the ratio of 5:3:1:1 (v/v). The size and height of plant materials used in all experiments were almost the same. Before drought stress, samples were cultured under the following growth conditions: 50% relative humidity, 25 • C/18 • C, day/night temperatures in the same artificial climate chamber in Beijing Forestry University (BJFU) (116.3 • E, 40.0 • N). Then, the samples were subjected to a continuously dry environment (no-watering) for 15, 30, 45, 60, 75, and 90 days. Then, after 7 days with rewatering treatment, the leaves were totally recovered. Fully watered plants were used as a control. Refer to experimental methods of drought treatment from David Gamboa-Tuz et al. (2018) . At each time point, the fourth or fifth fully expanded leaves from the plants were collected, then immediately frozen in liquid nitrogen and stored at −80 • C until the physiological measurement or RNA extraction. Each treatment group contained three replicates, and each sample from the three replicates was used to obtain biological replicates.

Before transcriptome sequencing using leaf tissue, which was collected and stored at −80 • C from plants subjected to 0 to 90 days of drought stress, the relative water content of leaves, superoxide dismutase (SOD) activity, MDA, and relative soil water content (RSWC) was evaluated. RSWC was analyzed using direct drying method. The relative water content of leaves was determined by the initial fresh weight of the leaf samples, referencing from Zhou et al. (2019) with slight modifications, the samples were soaked in distilled water placed on the lab bench with dark cloth for 24 h and finally weighted, The thiobarbituric acid (TBA) reaction was used for evaluation of the MDA content as described by Heath and Packer (1968) with minor adjustments. The SOD activity was measured after the spectrophotometer by determining its ability to inhibit the photochemical reduction of nitro blue tetrazolium (NBT). Each experiment included three biological replicates. The correlation analysis for various physiological parameters used the Pearson correlation coefficient in SPSS software.

According to the manufacturer's instructions, total RNA was isolated with an Easy Spin Plus RNA Extraction Kit (RN53, Aidlab China). RNA integrity was verified using electrophoresis on 1% agarose gel. The 2100 Bioanalyzer (Agilent Technologies, United States) and ND-2000 (NanoDrop Technologies) were employed to test the quality and purity of RNA.

The RNA-Seq library was generated using the NEB#7530 RNA Library Prep Kit (#E7530, New England Biolabs). The total mixed RNA from R. chinensis leaves in all treatments was used for cDNA library construction. The library quality was subsequently detected through the High Sensitivity DNA assay Kit (Agilent Technologies). The three biological replicates from each experiment were sequenced using an Illumina HiSeq TM 4000 platform by Genedenovo Biotechnology, Co., Ltd. (Guangzhou, China).

To investigate genes involved in the transcriptome, a cDNA sample was prepared from an equal mixture of t all RNA. The total RNA was isolated from the leaves of 2-year-old cutting seedlings R. chinensis from the control, drought-treated (30, 60, 90 days) and rehydration groups, designated as CK1, CK2, CK3, DT1-1, DT1-2, DT1-3, DT2-1, DT2-2, DT2-3, DT3-1, DT3-2, DT3-3, RW1, RW2, RW3. Three biological replicates from each group were analyzed on the Illumina HiSeq 4000 platform.

To get high quality clean reads, we were further filtered by fastp (version 0.18.0). The parameters were as follows: first removing reads containing adapters, and then removing reads containing more than 10% of unknown nucleotides (N), finally, removing low quality reads containing more than 50% low quality (Q-value ≤ 20) bases.

HISAT2.2.4 software (Kim et al., 2015) was used to carry out comparative analysis based on the rose genome (Raymond et al., 2018) . The comparison results of HISAT2 comparison StringTie (Pertea et al., 2015 (Pertea et al., , 2016 were used to reconstruct the transcript, and the expression of all genes in each sample was calculated.

The Fragments Per Kilobase of Exon Per Million Mapped Fragments (FPKM) method was used to analyze the expression level of each transcript. The reference-based approach was used to map the reads of each sample assembled by StringTie v1.3.1. For each transcription region, an FPKM value was calculated by StringTie software to quantify its expression abundance and variations. Differential expression analysis was performed by DESeq2 software (Love et al., 2014) between two different groups [also by edgeR (Robinson et al., 2010) between two samples]. The genes with the parameter of false discovery rate (FDR) below 0.05 and absolute fold change ≥ 2 were considered DEGs. GO enrichment analyses of DEGs (Ashburner et al., 2000) were able to recognize the main biological functions. Then, KEGG (Ogata et al., 1999) pathway enrichment analysis identified significantly enriched metabolic pathways or signal transduction pathways in DEGs.

Co-expression networks were established using WGCNA (v1.47) package in R (Langfelder and Horvath, 2008) . Gene expression values were imported into WGCNA to construct co-expression modules using automatic network construction function blockwise Modules with default settings, except that the power was 10, TOMType unsigned, and minModuleSize was 50. The correlation coefficient between module eigengenes and physiological data was calculated to find out significant modules, and the hub TFs were selected from those modules . The networks were visualized using Cytoscape_3.3.0 (Shannon et al., 2003) , GO and KEGG enrichment analyses were conducted for genes in each module.

The calculated p-value was subjected to FDR and correction, taking FDR ≤ 0.05 as a threshold.

Total RNAs were isolated from the leaves of the treatment and controlled specimens as described above. According to the manufacturer's instructions, first-strand cDNA synthesis was performed using Prime Script II 1st strand cDNA Synthesis Kit (Takara, Shiga, Japan). Based on the manufacturer's protocol, the qRT-PCR was performed using a Bio-Rad/CFX Connect TM Real-Time PCR Detection System (Bio-Rad, CA, United States) with SYBR R qPCR mix (Takara, Shiga, Japan). Relative mRNA content was calculated using the 2 − Ct method against the internal reference gene RcPP2A (Klie and Debener, 2011) . The primers used in this study were designed with Primer Premier 5 and are listed in Supplementary Table 12 . Three biological replicates were performed for all reactions.

The experiment applied drought stress by stop watering. Plant leaf phenotypes were recorded and physiological indexes were measured under control condition (CK) at 15, 30, 45, 60, 75, 90 days after stress and rewatering. Leaves were dark green and shiny during normal conditions ( Figure 1A) . After 15 days of drought treatment, the leaves has no obvious morphological changes ( Figure 1B) . Similarly, the RSWC and leaf relative water content were not significantly different from the control (Figures 2A,B) . At 30 days the leaves were slightly wrinkled ( Figure 1C ). When drought stress lasted for 45 days and 60 days, the leaves were all curled (Figures 1D,E) . At this time the relative water content of the soil considerably decreased, reaching a moderate degree of stress (Figure 2A) . At 60 days, the relative water content of the leaves decreased significantly compared with the control, which was consistent with the phenotype (Figure 2B) . The curl degree of the leaf deepened and the leaf became thinner, and the relative water content of the leaf has a material decrease during 75 days of drought (Figures 1F, 2B) . Both leaf relative water content and soil relative water content reached the lowest value, and the leaves were severely curled and withered after 90 days of drought (Figures 1G, 2A,B) . After rewatering, the leaves were green and thin ( Figure 1H ). In order to determine the exact point for transcriptome sequencing, physiological indicators were also measured. The MDA content of leaves treated with drought and rewatering was determined. With the deepening degree of drought stress, the content of MDA rises in fluctuation and reached the highest point after 90 days of the drought treatment ( Figure 2D ). In addition, SOD increased with the increasing of drought treatment time and achieved the maximum value after 90 days ( Figure 2C ). Both MDA and SOD content increased significantly during drought treatment for 30 days, 60 days, and 90 days. We conducted correlation analysis on physiological data. MDA content and SOD activity were positively correlated with different treatments and were significantly correlated at the 0.01 level (0.885, 0.979), while the relative leaf water content and soil relative water content were negatively correlated, at 0.01 significantly correlated horizontally (−0.891, −0.987). After rewatering, the MDA and SOD contents were all reduced. Based on these results, leaves were selected under normal conditions, drought treatment for 30 days, 60 days, 90 days, and rewatering treatment for transcriptome sequencing.

The clean reads were obtained by removing adaptor sequences, ambiguous nucleotides, and low-quality reads (Q-value < 20), which ranged from 41397356 to 54756474. The Q20 and Q30 values of all 15 libraries were more than 97% and 92%. The GC content was approximately ranged from 48% to 49%. Between 92.33% and 93.89% of the sequenced reads could be aligned to the rose reference genome ( Table 1) . In this study, the proportion of exon sequences ranged from 95.30% to 96.38%, and the proportion of intron sequences ranged from 1.45% to 2.66% (Supplementary Figure 1) .

To analyze the DEGs between leaves that suffered from drought stress at four different time duration (30, 60, 90 days and rewatering) and the control, comparisons were conducted between five groups. The DEGs were filtrated according to an expression level | log 2 (FC)| > 1 and FDR < 0.05 in each pairwise comparison. 6380 genes were identified exclusively at DT3, while 866 genes were specific in DT1. Especially 6380 (35.14%) genes were identified in DT3, but only 1680 (9.25%) and 1804 (9.94%) genes were identified in DT1 and DT2 ( Figure 3B ).

The number of DEGs is an increase from DT1 to DT3 and reached a peak at the CK vs. DT3 comparison set, indicating most of the drought regulated genes had a late response. The results showed DEG enrichment in CK vs. DT1, CK vs. DT2, and CK vs. DT3, illustrating the number of down-regulated DEGs was significantly higher than the number of upregulated DGEs. Meanwhile, a total of 11545 DEGs including 7886 downregulated and 3659 upregulated between CK and DT3 which contains the largest number of observations among all drought stress groups. Only 2503 upregulated and 2367 downregulated were observed between CK and DT1, suggesting that the differentiation of expressed genes between CK and DT3 is larger than that of CK and DT1. The result indicates that transcript abundance changed dramatically at key switched among the drought stages, which drought response genes could be induced and expressed largely. Compared with CK, there are slightly more genes upregulated than genes downregulated in RW ( Figure 3A ).

The drought resistance process in R. chinensis is complex According to GO functional enrichment analysis, as the drought stress deepened, the most enriched GO category among these DEGs was 'metabolic process, ' followed by 'cellular process, ' 'biological regulation, ' 'regulation of biological process, ' 'response to stimulus, ' 'catalytic activity, ' 'binding, ' 'transport activity, ' 'cell part, ' 'cell, ' 'organelle, ' 'membrane, ' and 'membrane part.' Meanwhile, based on the GO analysis of up-regulated DEGs, the GO terms of up-regulated DEGs in CK vs. DT1, CK vs. DT2, and CK vs. RW were mostly concentrated on membranes, including the cellular components category, such as 'membrane' (GO:0016020), 'membrane part' (GO:0044425), and 'intrinsic component of membrane' (GO:0031224) (Supplementary Table 1 ). CK vs. DT3 also involved more cells and metabolic processes including the biological process category (metabolic process, GO:0008152, cellular, GO:0009987, single-organism process, GO:0044699, cellular metabolic process, GO:0044237) (Supplementary Table 1) . Moreover, the GO analysis of the down-regulated DEGs found that among the top 20 enriched, 13 of them belong to the category of biological processes, mainly focusing on the response to stimuli in CK vs. DT1 (Supplementary Table 2 ). It is worth noting that the GO analysis of DT3 vs. RW mainly focused on cells and cell membranes, indicating that the rose has been restored to a certain extent at the cell level (Supplementary Tables 1, 2) .

The KEGG pathway analysis revealed that it is mainly enriched in the metabolism under the different drought degrees and rewatering treatment. Among the top 10 enriched pathways, 'plant-pathogen interaction' (ko04626), 'biosynthesis of secondary metabolites' (ko01110), 'plant hormone signal transduction' (ko04075), 'MAPK signaling pathway-plant' (ko04016), were regulated in response to DT1, indicating that drought stress was sensed and signaling transduction pathways were activated. However, more pathways related to metabolism were enriched at DT2 and DT3. The 'photosynthesis' (ko00195) and 'photosynthesis-antenna proteins' (ko00196) associated with photosynthesis were all enriched in CK vs. DT2 and CK vs. DT3 (Supplementary Table 3) . Metabolite-related pathways such as 'biosynthesis of secondary metabolites' (ko01110), 'alpha-Linolenic acid metabolism' (ko00592), 'tryptophan metabolism' (ko00380), 'glycosphingolipid biosynthesis-lacto and neolacto series' (ko00601), 'alanine, aspartate and glutamate metabolism' (ko00250), 'galactose metabolism' (ko00052), and 'glyoxylate and dicarboxylate metabolism' (ko00630) were enriched in CK vs. DT2 and CK vs. DT3. These results indicated that the metabolism processes of R. chinensis may act as the predominant process in the middle and later stages of drought stress. This result was unsurprising, as drought causes changes in the metabolism of plants to improve their drought tolerance. Meanwhile, secondary metabolism, primary metabolism, and photosynthesis were focused in DT3 vs. RW, indicating that the metabolism and photosynthesis of the rose after rewatering recovered, and that drought stress did not cause serious damage to the rose (Supplementary Table 3 ). Apart from metabolism enrichment, 'plant-pathogen interaction' (ko04626) was also enriched in CK vs. RW, indicating that plants are vulnerable to pathogenic bacteria after recovering from drought stress (Supplementary Table 3 ).

In this study, signal transduction in plants was influenced by drought stress which was verified by both the GO and KEGG enrichment analysis of the DEGs. More than 240 genes were predicted to encode protein kinases with varying expression levels in R. chinensis. Genes encoding receptor-like protein kinases (RLKs), in which LRR receptor kinases were the majority accounted for the largest proportion of these genes. Twenty-nine of them were up-regulated in DT1 and DT2. Among them, MIK1 (RchiOBHm_Chr6g0286591) was up-regulated 8.02-fold, respectively. More than 66 RLK genes were included in the group of serine/threonine-protein kinases during the drought stress. Furthermore, one gene related to G-type lectin S-receptor-like serine/threonine-protein kinase showed the highest degree of upregulation in DT3. In the KEGG pathway analysis of DEGs, the 'MAPK signaling pathway' (ko04016) was induced during drought treatment. For DEGs detected in DT1 and DT2, two MAPKKKs including RchiOBHm_Chr2g0158821 and RchiOBHm_Chr5g0002581 and four

were up-regulated. MKK10 (RchiOBHm_Chr7g0211301) was particularly up-regulated 9.18-fold. However, three

were down-regulated. Also, MPK4 (RchiOBHm_Chr7g0196781) and MPK19 (RchiOBHm_Chr5g0028971) were substantially increased in DT3 (Supplementary Table 4) .

Ca 2+ signal transmission is achieved through the generation, decryption, and transmission of specific calcium signals and the corresponding physiological and biochemical reactions downstream. The study identified DEGs related to calcium ion, encoding calcium-binding proteins, Ca 2+ -binding protein EF hand, calmodulin-like, calcium-dependent protein kinases, calcium-transporting ATPase, and calcium channel protein. One gene (RchiOBHm_Chr2g0158571) encoding the calcium channel protein sustained high expression under drought stress and was up-regulated 6.44-fold in DT1 vs. CK (Supplementary Table 5 ).

According to the GO and KEGG analyses of the DEGs, the GO term 'hydrogen peroxide metabolic' (GO:0042743) and 'response to reactive oxygen species' (GO:0000302) were enriched at DT1 and DT2, while peroxisome (ko04146) was induced at DT3. 69 DEGs encoded ROS production and scavenging. A total of 55 DEGs encoding enzymes related to ROS scavenging, including peroxidase (POD), ascorbate peroxidase (APX), glutathione S-transferase (GST), glutathione peroxidase (GPX), polyphenol oxidase (PPO), ferritin, glutaredoxin, thioredoxin, and peroxiredoxin. Among these antioxidant enzymes, POD45 (RchiOBHm_Chr4g0400821) encoding POD was up-regulated by 4.79-fold at DT1, another gene encoding POD PNC1 (RchiOBHm_Chr5g0072281) has no significant change from DT1 to the control, but 9.29-fold has a significant at DT2, reaching the peak of expression level at DT3 (Supplementary Table 6 ).

The expression of genes involved in hormone biosynthesis or signaling was substantially changed during the drought stress and rewatering stage. The changed genes included ABA, auxin, BR, cytokinin, ET, JA, and SA. Among these DEGs, most genes were involved in ABA response. Ninecis-epoxycarotenoid dioxygenase (NCED) is a critical enzyme in ABA biosynthesis. For instance, one gene encoding NCED3 (RchiOBHm_Chr5g0014331) was upregulated by 2.70-fold at DT1, while another gene encoding NCED6 (RchiOBHm_Chr4g0397001) was significantly down-regulated 9.40-fold at DT3. The expression of a gene (RchiOBHm_Chr5g0049951) encoding AAO is also a vital enzyme in ABA biosynthesis, which was up-regulated by 2.78-fold at DT1 and down-regulated by 7.28-fold at DT3. Moreover, some genes involved in the ABA mediated signaling pathway were also induced by drought stress. Abscisic acid receptor PYL4 (RchiOBHm_Chr1g0371101) was down-regulated. In contrast, one gene encoding PP2C (RchiOBHm_Chr5g0066401), a negative regulator of ABA, has an up-regulation (Supplementary Table 7) .

A total of 24 genes involved in the auxin signaling pathway changed. Three genes encoding auxin response factors (ARF3, ARF5, ARF9) were drought induced. Particularly, ARF3 (RchiOBHm_Chr5g0009381) expression continued to decline under drought stress, which was induced 4.41-fold at DT3. At the same time, the expression level of ARF5 (RchiOBHm_Chr6g0302551) continued to increase under drought stress. Two genes encoding AUX/IAA proteins were up-regulated at DT1 and down-regulated at DT3. AUX22D (RchiOBHm_Chr4g0389611) was lowered nearly 10-fold at DT3. Furthermore, some genes participated in auxin response and signal, for example of auxin-binding protein ABP19alike (ABP19A), auxin transporter-like protein 2 (LAX2), and auxin efflux carrier (PIN1C). The gene encoding ABP19A (RchiOBHm_Chr2g0094781) was up-regulated 7.41-fold at DT1. Five genes encoding SAUR were identified, among these, SAUR77 (RchiOBHm_Chr5g0026611) was up-regulated 9.78-fold at DT1 (Supplementary Table 7) .

In this study, a large number of ethylene-responsive transcription factor (ERF) genes displayed a down-regulation trend under early drought stress. However, one gene encoding ERF023 (RchiOBHm_Chr1g0360021) had an upregulation 3.70-fold at DT1 and was then down-regulated 5.02-fold at DT3. Another gene encoding example, ERF110 (RchiOBHm_Chr6g0274591) was up-regulated 7.34-fold at DT3. Meanwhile, ethylene receptor (ERS1), ethylene-overproduction protein 1 (ETO1), and ethylene insensitive3-like protein (EIL3) were identified (Supplementary Table 7) .

Most of the genes associated with the biosynthesis or signaling of jasmonic acid (JA) were down-regulated. The gene encoding malate dehydrogenase (MDHG) was down-regulated 4.06-fold at DT3, while 12-oxophytodienoate reductase 3-like (OPR3) was up-regulated 4.13-fold at DT3. Additionally, the majority of genes were related to cytokinin, salicylic acid, and brassinosteroid pathways. One gene encoding cytokinin dehydrogenase 3 (RchiOBHm_Chr2g0093491) was down-regulated 7.40-fold at DT3 (Supplementary Table 7) .

The output of GO and KEGG enrichment analysis indicates that many DEGs were related to metabolism and biosynthesis. The 'polysaccharide catabolic process' (GO:0000272) and the 'Starch and sucrose metabolism' (ko00500) pathway were significantly induced by drought stress. Related to carbohydrate synthase and starch synthase were induced. The expression level of sucrose synthase increased, and one of the genes encoding SUS7 (RchiOBHm_Chr3g0488651) up-regulated 4.03-fold at DT3. Only a small portion of starch synthase genes were downregulated, others appeared to be up-regulated. Additionally, drought treatment significantly induced the expression of enzyme genes involved in starch and sugar metabolism, fructose, and mannose metabolism. Among these, genes encoding fructose-1,6-bisphosphatase, fructose-bisphosphate aldolase, mannose-1-phosphate guanyltransferase alpha, pyrophosphatefructose 6-phosphate 1-phosphotransferase subunit beta-like and 6-phosphofructo-2-kinase/fructose-2,6-bisphosphatase were down-regulated at DT3 (Supplementary Table 8) .

In total, 89 DEGs associated with lipid metabolism were enriched in the 'lipid metabolic process' (GO:0006629). Most of these were related to the biosynthesis and metabolic processes of lipoprotein, wax, keratin, lipid, and fatty acid. Moreover, most genes involved in wax biosynthesis were significantly up-regulated. For example, one gene encoding CYP86B1 (RchiOBHm_Chr4g0432631) up-regulated 9.90-fold at DT3 (Supplementary Table 9 ). The KEGG functional annotation of these differentially expressed genes was related to metabolism and differentially expressed genes were also related to secondary metabolisms. A total of 95 DEGs involved in secondary metabolism were enriched by KEGG function, and the results are shown in Figure 4 . Drought stress affected the transcription level of genes that related to phenylpropanoid biosynthesis pathways.

Many DEGs enriched in the category of 'photosynthesis' (ko00915) were overexpressed in DEGs down-regulated by DT2 and DT3. In the process of DT2 and DT3, many genes involved in photosynthesis were significantly down-regulated. A total of 70 DEGs involved in photosynthesis were identified. In the seam, three (CAB40, LHCB5, psaA) of 16 DEGs encoding the components of photosystem I (PSI) were slightly up-regulated in DT1, and other genes were significantly down-regulated during the drought treatment period. A total of 33 genes encoding PS II were identified, including 6 up-regulated genes (psbE, three genes encoding psbC, psbD, PSB27-1), and 27 down-regulated genes. Most of these genes decreased significantly in the middle and later periods of drought, indicating that PS I and PS II were suppressed under drought stress. Meanwhile, there were 15 genes encoding redox chains, all of which were significantly down-regulated in the middle and late periods of drought stress. Additionally, the expression of six gene encoding components of chloroplasts significantly decreased (Supplementary Table 10 ).

Transcription factors specifically binding to cis-acting regulatory elements in the promoter of target genes are crucial regulatory proteins, which can modulate numbers of genes up or downregulation (Dhriti and Ashverya, 2015) . In this study, WRKY (58), NAC (43), MYB (52), bZIP (19), bHLH (55) families contained large numbers of transcription factors. The expression level of the WRKYs family exhibited large changes under different treatments. The expression of some genes reached a maximum in the middle and later periods of drought. Genes encoding WRKY75 and WRKY71 all reached the highest expression level at DT3 with an up-regulation of more than 5-fold. WRKY57 (RchiOBHm_Chr5g0083891) up-regulated regarding expression responses at different drought stress, especially increased 7.93fold at DT3. However, more than half of the WRKY families have decreased expression levels at different drought levels. Twentyfour genes encoding WRKYs were down-regulated at DT1 by 1.02-to 3.71-fold (Supplementary Table 11) .

More than half of the genes encoding NACs were upregulated to varying degrees under drought stress. One gene encoding NAC100 (RchiOBHm_Chr2g0167461) continued to increase in the middle and later periods of drought stress by 7.20-fold. A total of 11 genes upregulated during DT1. Some genes increased and then decreased under drought stress, but finally increased after rewatering. For example, NAC073 (RchiOBHm_Chr0c19g0500091) was upregulated 5.78-fold at DT1, decreasing slightly at DT2, then down-regulated 2.64-fold at DT3, but had an increase of 6.96-fold after rewatering. NAC037 (RchiOBHm_Chr2g0099671) had a similar expression trend, but the up-regulation decreased a small amount. On the contrary, one gene encoding NAC072 (RchiOBHm_Chr5g0034761) was down-regulated at DT1 and continued to increase afterward (Supplementary Table 11) .

Under different drought treatments, more than threequarters of the MYB genes' expression levels showed different degrees of up-regulation. Among them, a relatively large number of up-regulated genes appeared during the DT1 period. For instance, MYB46 (RchiOBHm_Chr1g0315931) and MYB61 (RchiOBHm_Chr1g0315931) showed the highest upregulation of expression with 7.09 and 7.15-fold. Moreover, some genes up-regulated significantly at DT3. MYB75 (RchiOBHm_Chr2g0116041) displayed the most amount of up-regulation among all up-regulated genes at DT3, with an increase of 7.65-fold. Some genes showed a tendency to increase and then decrease during drought treatment, for example, gene encoding MYB61 (RchiOBHm_Chr3g0458721) slightly increased at DT1 and then decreased by 6.10-fold. After rewatering, most gene expression did not change significantly compared with CK, only a small amount of gene expression up-regulated (Supplementary Table 11) .

Among the three drought treatments of the bHLH families, the numbers of genes up-regulated in one period were more than the numbers of genes down-regulated. The number of downregulated genes in two periods was more than that. Most genes in the bHLH family up-regulated during early drought treatment of DT1. While a large number of genes significantly downregulated during the middle and later drought periods (DT2 and DT3). Among these, there was a slight upward adjustment, up-regulated by 1.02-to 4.66-fold at DT1. Nevertheless, there were 37 genes down-regulated at DT2 or DT3. One gene encoding bHLH72 (RchiOBHm_Chr7g0182341) down-regulated 6.94-fold. A gene bHLH121 (RchiOBHm_Chr6g0245181) that belongs to the bHLH family has no significant change in gene expression during drought treatment but appeared to be upregulated after rewatering (Supplementary Table 11) .

Among AP2/DREB and bZIP transcription factors, more genes exhibited down-regulation under drought treatment. Nineteen members of bZIPs were found to be responsive to drought stress, while 12 down-regulated from 1.24-to 4.80-fold (Supplementary Table 11 ). Five genes encoding DREBs showed down-regulation that exceeded 6.60-fold at DT1. The most significant downregulation occurred by DREB1D (RchiOBHm_Chr7g0199351), which down-regulated 9.64-fold (Supplementary Table 11 ).

Gene co-expression network gene clustering and module cutting combined genes with similar expression patterns on the same branch. Each branch represented a co-expression module with different colors representing different modules. The unexpressed genes in more than half of the samples were filtered, the expression patterns of 30,012 genes obtained from transcriptome sequencing were performed by WGCNA. A total of 18 modules were identified according to the similarity of expression patterns ( Figure 5A ). Cluster analysis was performed to evaluate genes in the modules and were observed through heat maps ( Figure 5B) . Furthermore, to identify modules that were significantly associated with different drought stress, the module-trait correlation relationships were constructed ( Figure 5C ). There were four modules including blue, navajowhite1, salmon4, and coral2 were closely associated, with traits related to four treatments (0.99, 0.99, 0.87, and 0.76, respectively). The expression pattern of the eigengenes represents the gene expression profile of the entire module, and the expression patterns of the four selected modules were analyzed (Figure 6) .

Four modules of blue, navajowhite1, salmon4, and coral2, were annotated through GO and KEGG analysis. The blue module of DT1 was most significant in GO analysis, which was mainly enriched in signal transduction and stress response (Supplementary Table 13 ). The results of KEGG analysis were Table 13 ). The navajowhite1 module showed significant rich GO terms, indicating many genes were enriched in various stress responses. Those that were mainly enriched in 'response to stimulus' (GO:0050896), 'response to chemical' (GO:0042221), 'response to water' (GO:0009415), and 'response to water deprivation' (GO:0009414) were also significantly enriched ( Figure 7A) . KEGG enrichment analysis displayed the main enrichment in metabolic pathways ( Figure 7B ). In the salmon4 module, GO terms and KEGG pathways were enriched in metabolic progress, which indicates that the rose produces metabolites in the later drought stage to resist drought stress (Supplementary Tables 13, 14) . However, the process of compound metabolisms was focused on the coral2 module (Supplementary Tables 13, 14) . The GO terms about stress response are mainly concentrated in navajowhite1 and blue module. As a result, navajowhite1 and blue module genes may play a critical role in the process of rose drought stress and rewatering.

Transcription factors are an essential class of regulatory proteins in biological processes. This study analyzed the TFs in each module. Although there was a difference in the distribution of TFs in different modules, it is mainly concentrated in WRKY, NAC, MYB, ERF, bHLH, bZIP, C2H2, and other transcription factor families. According to previous reports, these TFs were involved in the process of plant stress regulation (Erpen et al., 2018) . TFs of navajowhite1 and the blue module were analyzed, the highly connected examples that were used as core genes were screened out. Many DEGs Identified in the navajowhite1 module were annotated as TFs, including several WRKY, MYB, NAC, ERF, ARF, and bHLH TFs ( Figure 8A) . Based on the hub gene correlation network and their high connectivity, 22 of these TFs were selected: two genes (RchiOBHm_Chr5g0034761 and RchiOBHm_Chr6g0308271) encoding NAC072, two genes (RchiOBHm_Chr2g0167461 and RchiOBHm_Chr7g0181811) encoding NAC100, NAC087 (RchiOBHm_Chr3g0497471), (RchiOBHm_Chr6g0288981) (Figure 8B) . The experiment showed the expression of these genes was significantly upregulated at DT3 except for WRKY6, bHLH3, ARF1, and ARF8 ( Figure 8B ). In the blue module, 20 TFs were identified as hub genes including ERFs, bHLHs, WRKYs, MYBs, and NACs ( Figure 8C ). ERF109 and bHLH162 have higher connectivity among them (Figure 8C) . The expression of all 20 TFs increased in DT2 (Figure 8C) .

To verify the accuracy and reproducibility of the transcriptome analysis, 13 DEGs were selected to analyze the transcript abundance using qRT-PCR, including seven randomly selected transcription factors from the navajowhite1 module, four TFs, and two ABA synthesis-related DEGs. The results showed that the expression profiles detected by qRT-PCR were positively correlated with the RNA-Seq results (Figure 9) . Therefore, the reliability of the RNA-seq data was confirmed by the consistency between the qRT-PCR results and RNA-seq analyses.

Weighted gene co-expression network analysis is an analysis method for analyzing the gene expression patterns of multiple samples. Genes with similar expression patterns can be clustered, moreover, the relationship between modules and specific traits can be analyzed (Zhang and Horvath, 2005) . WGCNA analysis can present the interaction relationship between genes, and identify the hub genes at the center of the regulatory network. It also can use the module eigengenes values to perform correlation analysis with specific traits, more accurately (Silva et al., 2019) . More and more studies use WGCNA to analyze key genes in plants under abiotic stress (Shao et al., 2021) . In this study, through the analysis of WGCNA and transcription regulation modules, the highly connected TFs in modules were selected. Numerous TFs were obtained in the transcriptome sequencing in this study including the NACs, WRKYs, MYBs, AP2/ERFs, bHLHs, and so on. Many studies have shown that NACs play positive roles in the regulation of plant stress resistance. In particular, Arabidopsis (Tran et al., 2004) , wheat (Mao et al., 2014; Zhang et al., 2015) , and rice (Hong et al., 2016; Sung et al., 2018) have received the most research attention on NACs, which found that overexpression of NACs transcription factor family genes can significantly enhance the tolerance of transgenic plants to multiple stresses. In addition, rose (Rosa hybrida) RhNAC3 can improve rose petals dehydration tolerance, and can also improve the drought tolerance of transgenic Arabidopsis . In the study of rose (Rosa chinensis 'Slater's crimson China') heat resistance, it was found that both NAC27 and NAC72 up-regulated . NAC27 and NAC72 were also involved in ABA signal transduction under heat stress of Rhododendron hainanensis . Overexpression of NAC72 in Arabidopsis thaliana can increase the drought resistance of plants, and also proved that it mediates the regulation of ABA-responsive genes . Gene encoding NAC072 with high connectivity in the navajowhite1 module was selected, indicating that it plays an important role in the regulation of rose drought stress. In agreement with previous studies, it can be concluded that NAC072 may be involved in ABA signal transduction under drought stress in R. chinensis. Furthermore, the role of the core gene in the navajowhite1 module NAC087 and NAC100 in drought also requires further consideration.

The WRKY genes play pivotal roles in stress responses. Many researchers have demonstrated that WRKYs have crucial biological functions in the response of plants to different kinds of biotic and abiotic stresses . Studies on Arabidopsis, rice, and soybean have shown the importance of the WRKYs in response to drought stress. For example, AtWRKY57 , AtABO3 (Ren et al., 2010) , OsWRKY47 (Raineri et al., 2015) , and GmWRKY54 (Wei et al., 2019) can improve the drought tolerance of plants. Likewise, VaWRKY14 responds to drought and cold stress, and hence the overexpression of VaWRKY14 can enhance the drought tolerance of transgenic Arabidopsis . However, some WRKYs can act as a negative regulator of abiotic stress in plants. For instance, GhWRKY33 (Wang et al., 2019b) and SlWRKY81 Li et al., 2020) reduce the drought tolerance of transgenic plants. In the analysis, eight WRKY genes including WRKY75, WRKY28, and WRKY27 were found to have high connectivity in the navajowhite1 module, indicating that they may be involved in drought response. Other research has found that PagWRKY75 was down-regulated in the early stages of salt and osmotic stress, transgenic poplar lines overexpressing PagWRKY75 were more sensitive to salt and osmotic (Zhao et al., 2019b) . It can be assumed that WRKY75 plays an essential regulatory role in abiotic stress, but the function of WRKY75 in drought stress screened in this study requires further verification.

APETALA2/ethylene-responsive factor (AP2/ERF) TFs have been found to regulate plants and to enable them to be resistant to abiotic stress (Li et al., 2021b) . In this study, AP2/ERF TFs were observed that may be involved in the regulation of rose drought stress. One gene encoding ERF109 identified in a blue module appeared to play an important role in roses in resisting drought stress, the function of which requires further study. The study showed that PtrERF109 is positively regulated with Poncirus trifoliata to resist the cold (Wang et al., 2019a) . Moreover, DREB2A, DREB3, and DREB1B were found in this study. More and more studies have shown the function of DREBs in resisting abiotic stress in plants. Overexpression of StDREB2 improved the drought resistance of cotton (El-Esawi and Alayafi, 2019). However, RhDREB2B is a negative regulator to resist abiotic stress (Li et al., 2021b) .

The function of bHLHs in resisting abiotic stress in Arabidopsis thaliana has also been studied Qiu et al., 2020) . bHLH162 and bHLH35 have been identified in this study, providing new candidate genes for roses to resist drought stress. Meanwhile, MYBs including MYBS3, MYB75, and MYB78 require further investigation. Particularly for rice, OcMYBS3 is essential for resisting cold stress (Su et al., 2010) .

Phytohormone plays a key role in abiotic stress responses and coordinates various signal transduction pathways. In this study, the gene expression involved in several plant hormone-related signal transduction pathways changed significantly after drought treatment, indicating that plant hormones may play a crucial role in the response of the rose to drought stress. ABA acts as an endogenous messenger in a plant's abiotic stress responses (Teng et al., 2014) . Drought results in a substantial increase in plant ABA levels, accompanied by a major changes in gene expression and adaptive physiological responses. NCED is a key rate-limiting enzyme in the ABA synthesis pathway (Schwartz et al., 1997; Iuchi, 2002) . Currently, some NCED3 genes have been isolated and identified, corresponding gene function studies have also shown that NCED3 genes play an important role in enhancing plant drought resistance (Li et al., 2019a) . The study found that an NCED3 gene was up-regulated by 2.70fold at DT1. Interestingly, NCED6 down-regulated by 9.40fold in DT3. Lefebvre et al. (2006) state that AtNCED6 is involved in ABA biosynthesis during seed development. The role of NCED6 in drought stress still requires further study. Another ABA synthesis related to gene aldehyde oxidase (AAO3) and zeaxanthin cyclooxygenase (ZEP) were identified, both of which were significantly down-regulated in DT3. These results suggested that the synthesis of ABA may increase during the early drought stage and decrease during the later drought stage. Therefore, ABA-mediated signal transduction may be involved in the drought response of R. chinensis.

There is more and more evidence suggesting that the phytohormone ethylene is essential for regulating various developmental processes and stress responses of plants (Pei et al., 2017) . Mcmichael et al. (1972) reported that the release of ethylene may increase under drought conditions. Studies on pineapples have shown that drought-stressed plants produce significantly less ethylene in leaf and stem tissues compared to control plants (Min and Bartholomew, 2005) . These conflicting results may be due to the differences in distinct species. In general, the production of ethylene can be affected by changing ACC biosynthesis under drought stress (Larrainzar et al., 2014) . In plants, ACC synthase (ACS) and ACC oxidase (ACO) are the two most essential enzymes in the ethylene biosynthetic pathway (Kende, 1989) . In the study, the expression of three ACS genes was suppressed under drought stress, the expression of one ACS gene up-regulated, while the expression of two ACO genes was also up-regulated. It is not yet possible to speculate whether the ethylene content has increased or decreased under drought conditions. Ethylene response factors (ERFs) are downstream components of the ethylene signaling pathway (Guo and Ecker, 2004) . Experiments showed that ERFs play important roles in plant abiotic stress response (Pan et al., 2012; Yu et al., 2017) . A large number of ERFs exhibit different expressions under different drought levels. Their role under drought stress needs further study.

Jasmonates (JAs), a class of oxygenated lipid derivatives, are phytohormones necessary for plant growth and environmental adaptation (Wasternack and Hause, 2013) . MYCs proteins are one of the JASMONATE-ZIM (JAZ) target proteins and a key transcription factor regulating the corresponding genes downstream of JA. According to the reports, MYC2, MYC3, and MYC4 are the central nodes of the JA signal regulation (Chico et al., 2020) . Lipoxygenase (LOX), allene oxide synthase (AOS), allene oxide cyclase (AOC), and oxophytodienoate reductase (OPR) are important enzyme genes in the JA synthesis pathway (Devoto and Turner, 2005) . The study has observed that RhHB1 binds to the RhLOX4 promoter to suppress its expression in the cut rose (Rosa hybrida). As a result, dehydration tolerance decreases (Fan et al., 2020) . These genes were identified, while most genes were up-regulated under drought stress. Other studies have shown that JA was observed to increase in Pinus pinaster (Pedranzani et al., 2007) plants and Oryza sativa (rice) (Du et al., 2013) leaves exposed to drought conditions. Based on these results, it can be conjectured that the JA pathway plays a role in the response to drought stress in R. chinensis.

Moreover, recent studies have indicated that the role of brassinosteroid (BR) in plant response to drought stress (Nolan et al., 2017a (Nolan et al., ,b, 2020 . In Fabregas et al.'s (2018) study, the overexpression of BRL3, a vascular-rich member of the BR family, can confer Arabidopsis drought resistance without compromising normal plant growth. In this study, the expression of two BRL genes increased under drought stress, and hence BRLs may play positive parts in the response to drought. In addition, BR and ABA regulated the drought resistance of plants. Indeed, the molecular basis of the antagonistic effect of the BR-ABA pathway has been defined Nolan et al., 2020) . The role of ABA in promoting the response to drought stress suggests that BRs will inhibit drought stress responses. Yet, only the role of hormone series resistance has been explored in model plants. This study also laid the foundation for the interaction of hormones in the drought resistance of rose.

When plants are subjected to adversity stress due to changes in the external environment, they sense adversity stimuli from outside, recognize and transduce between cells and cells, and then transmit them to downstream genes, thereby causing downstream affect gene expression to resist adversity. Calcium ions are the second messengers of the cell, which mediates the calcium signal plays an essential role in the plant response to external stimuli (Kudla et al., 2010; Reddy et al., 2011) . As a second messenger, calcium ion can not only maintain the stability of cell structure, but also play a vital role in signal transmission and participate in plant stress response (Xiong et al., 2002; Shinozaki et al., 2003) . The three main families of calcium sensors in Ca 2+ signaling in plants are calmodulin (CaM), calcineurin B-like (CBL), and calcium-dependent protein kinases (CDPKs) (Ranty et al., 2006; Boudsocq and Sheen, 2013) . Analysis of the CaM/CMLs identified in this study revealed that five genes up-regulated, two up-regulated and then downregulated, and fifteen down-regulated under drought stress. AtCML8, AtCML13, AtCML18 and AtCML25 down-regulated by salt damage and drought (McCormack et al., 2005) . ShCML44 up-regulated under cold, drought, and salt damage (Munir et al., 2016) . In contrast, the gene encoding CML44 down-regulated and the gene encoding CML8 up-regulated under drought induction in this study. CDPKs/CPKs are involved in multiple stress signaling pathways. In Arabidopsis, Overexpressed AtCPK6 confers drought tolerance (Xu et al., 2010) , while AtCPK21 is a negative regulator of osmotic response (Franz et al., 2011) . In this study, six genes encoding putative Ca 2+ -dependent protein kinases up-regulated during drought stress, and four of them showed high expression levels during the DT2 period. These results indicate that the CDPK pathway was significantly activated in DT2. CBLs communicate with downstream CIPKs for signal transmission. The CBL-CIPK pathways also participate in plant abiotic stress response (Xiang et al., 2007) . In the current research, three genes encoding CBLs-CIPKs significantly down-regulated in DT3, which indicates that the CBL-CIPK pathway may contribute to regulating a later drought response in R. chinensis. In summary, the above results suggest that Ca 2+ mediated signaling pathway has key functions in the response of R. chinensis to drought stress.

Plant's sense and mediate various intracellular and extracellular signals through cell surface receptor kinases after stimulating by the external environment. Studies have shown that receptor-like protein kinases (RLKs) play a vital role in plant growth and stress adaptation (Wang et al., 2007) . Leucine-rich repeat receptor protein kinases (LRR-RLKs), a class of single transmembrane proteins, are the largest family of plant receptor protein kinases (Gou et al., 2010) . The LRR-RLKs genes identified by the study were 1 up-regulated and 10 down-regulated under drought stress, which was similar to the expression of LRR-RLKs genes (3 up-regulated, 38 downregulated) in chrysanthemum under dehydration stress (Xu et al., 2013) . Additionally, multiple MAPKs were observed during the plant response to abiotic stresses such as salt, drought, cold, and heat, mainly rapid activation of MPK3, MPK4, and MPK6 (Zhu, 2016) . The results in this study are similar, as gene encoding MPK3 significantly down-regulated on DT1, indicating that MPK3 responded quickly when drought stress occurred.

MAPK cascade signals can also transmit reactive oxygen species (ROS) signals to downstream targets. When plants were subjected to environmental stress, ROS generated signals in the plant, triggering a series of changes. Plants also have an antioxidant system protecting against poisonous oxygen (Jaspers and KangasjcTrvi, 2010) . Abiotic stress can cause plants to produce excessive ROS which can cause damage to plants. Plants have efficient enzymes including SOD, POD, APX, GPX, and so on. The antioxidant defense mechanisms can protect plants from oxidative stress (Gill and Tuteja, 2010) . In this study, two genes encoding APX were up-regulated and four were down-regulated. Similarly, in the study of Citrullus lanatus, the changes of APX genes are different under drought stress (Goitseone et al., 2018) . The post-harvest petals of cut rose (Rosa hybrida cv. Samantha) under water deficit stress reported that the Transcription level of RhAPX may be involved in the response of water deficit stress (Jin et al., 2006) . Interestingly, three genes encoding SOD were identified, all of which were reduced under drought stress. When SOD activity was measured, the level of SOD activity increased as the degree of drought stress increased. The expression pattern of SOD genes in this study is similar to the expression pattern of the SmFSD2 in drought. The SmFSD2 did not change significantly over 3 days, but decreased to the lowest value at 9 days (Shafi et al., 2019) . These results indicate that the ROS signaling pathway may be involved in the regulation of R. chinensis response to drought stress.

Primary and secondary metabolites act as signaling molecules or protective agents when plants respond to adversity stress. The primary metabolic processes were involved in the response to the drought stress of R. chinensis, producing a series of osmotic protective agents. Such as sugars, starches, amino acids, and lipids. Accumulation of sugars in various plants is related to high tolerance to drought stress (Ali et al., 2019) , which is consistent with the observation of the experiment. Most genes related to carbohydrate synthesis including fructose, glucose, mannose, sucrose, and trehalose were induced and significantly up-regulated under drought stress. The gene encoding trehalose phosphate synthase (TPP1) down-regulated under drought stress. Among the 11 OsTPS genes in the rice (Oryza sativa) genome, only OsTPS1 has TPS activity . Overexpression of OsTPS1 can improve the tolerance of rice seedlings to drought stress. At the same time, it also causes the expression of some stress-related genes to be up-regulated and the phenotype has no obvious change .

Trehalose may involve in the process of rose's resistance to drought. The expression levels of the four genes encoding starch synthase significantly reduced at DT3, with no significant change in expression levels in the remaining periods compared with the control. It is speculated that the starch content may decrease in a later drought.

Lipids are essential components of cells and organelles. A series of genes that synthesize lipid biosynthesis have been identified in drought treated R. chinensis plants. Long-chain acyl-CoA synthetases (LACS) are indispensable in the pathway of lipids synthesis and degradation of higher plants. There are nine genes in the LACSs family in Arabidopsis, among which LACS1 and LACS2 are involved in the synthesis of wax and cutin together, with an overlap in function (Lü et al., 2009; Weng et al., 2010) . Drought stress can cause rapid accumulation of wax on plant surface . In this study, six genes encoding LACS2 were induced under drought, while five of them significantly down-regulated at the DT3. Studies on bananas, Agave sisalana (Sarwar et al., 2019) , and Poa pratensis (Ni et al., 2016) have reported that waxy synthetic transcripts significantly down-regulated under drought stress. These results confirm that wax may play key roles in regulating the drought resistance of R. chinensis. In brief, it can be speculated that the wax content may be reduced in the later period of drought stress. In addition, the Frontiers in Genetics | www.frontiersin.org gene encoding GDSL esterase/lipase was significantly induced at the beginning of drought stress. Studies on the GDSL esterase/lipase gene of Arabidopsis thaliana revealed that this gene helps plants resist abiotic stress (Lai et al., 2017) , and is induced in M. wufengensis by cold stress (Deng et al., 2019) . In summary, lipid metabolism is crucial for R. chinensis adaptation to drought stress.

Furthermore, as part of the adaptation mechanism to the environment of plants, secondary metabolism is sensitive to both biological and abiotic stresses (Hartmann, 2007) . It is known that the phenylpropane pathway has high responsiveness to different abiotic stresses such as injury, drought, and low or high temperature. White grapes (Vitis vinifera L.) cope with drought by stimulating the phenylpropane pathway, which reflects the result of this study. The experiment showed that most of the secondary metabolites are derived from phenylpropane biosynthesis and flavonoid biosynthesis pathways. The genes that encode phenylalanine ammonia lyase (PAL) and coumaric acid-CoA ligase (4CL) were identified. They are all related to enzyme genes in the phenylpropane metabolism pathway. They significantly up-regulated at DT1, indicating that the phenylpropane metabolism pathway was induced under drought stress. The transcription of PAL of Caragana korshinskii increased under field and laboratory drought conditions . These results also confirm the results of the present report. Another major derivative of the phenylpropane pathway is flavonoids, which play a key role in adversity stress. Flavonoids with free radical scavenging activity alleviate oxidative and drought stress in Arabidopsis thaliana (Nakabayashi et al., 2014) . This study reported that genes encoding flavonoid synthase (FLS), chalcone isomerase (CHS), and anthocyanin synthase (ANS) were induced under drought. Flavonoid products may affect the regulation of R. chinensis by drought stress.

This study aimed to enhance understanding of the mechanisms through which R. chinensis responds to drought and identify the TFs related to drought tolerance. Figure 10 summarizes the model of R. chinensis drought response. According to the transcriptome data, signal receptors, such as ion channel proteins sense external drought stimuli and transmit through Ca 2+ , ROS, and phytohormone signaling transduction. The signal transmits protein kinases and protein phosphorylases, which can activate the transcription factors. The activation of TFs triggers downstream drought-responsive gene transcription such as lipid metabolism, carbohydrate metabolism, and secondary metabolism to regulate cell homeostasis. The hub TFs involved in roses' drought response and rewatering were identified through WGCNA, providing resources for further research on the drought control network of roses. The new findings of this study, on the relationship between unreported TFs and drought resistance regulation, still require further experimental explore.

The datasets presented in this study can be found in the NCBI repositories (BioProject ID: PRJNA722055).

YL, SZ, XJ, and HF conceived and designed the research. XJ performed the experiments, analyzed the data, and wrote the manuscript. YB and NJ cultivated and provided the rose plants. All authors have read and agreed to the published version of the manuscript. 

Metabolic and genes expression analyses involved in proline metabolism of two rose species under drought stress

Tomato WRKY81 acts as a negative regulator for drought tolerance by modulating guard cell H2O2-mediated stomatal closure

Trehalose Metabolism in Plants under Abiotic Stresses: Profiling and CounterAction

Impact of drought on growth, photosynthesis, osmotic adjustment, and cell wall elasticity in Damask rose

Gene ontology: tool for the unification of biology

CDPKs in immune and stress signaling

Understanding plant responses to drought -From genes to the whole plant

CUL3(BPM) E3 ubiquitin ligases regulate MYC2, MYC3, and MYC4 stability and JA responses

RhNAC2 and RhEXPA4 are involved in the regulation of dehydration tolerance during the expansion of rose petals

Transcriptomics and co-expression networks reveal tissue-specific responses and regulatory hubs under mild and severe drought in papaya (Carica papaya L.)

De novo transcriptome sequencing and gene expression profiling of Magnolia wufengensis in response to cold stress

Jasmonate-regulated Arabidopsis stress signalling network

Transcriptional regulation of drought response: a tortuous network of transcriptional factors

RhNAC31, a novel rose NAC transcription factor, enhances tolerance to multiple abiotic stresses in Arabidopsis

Endogenous auxin and jasmonic acid levels are differentially modulated by abiotic stresses in rice

Overexpression of StDREB2 transcription factor enhances drought stress tolerance in cotton (Gossypium barbadense L.)

Potential use of the DREB/ERF, MYB, NAC and WRKY transcription factors to improve abiotic and biotic stress in transgenic plants

Overexpression of the vascular brassinosteroid receptor BRL3 confers drought resistance without penalizing plant growth

The RhHB1/RhLOX4 module affects the dehydration tolerance of rose flowers (Rosa hybrida) by fine-tuning jasmonic acid levels

PeCHYR1, a ubiquitin E3 ligase from Populus euphratica, enhances drought tolerance via ABA-induced stomatal closure by ROS production in Populus

DNA-binding specificities of plant transcription factors and their potential to define target genes

Calcium-dependent protein kinase CPK21 functions in abiotic stress response in Arabidopsis thaliana

Physiological traits determining high adaptation potential of sweet briar (Rosa rubiginosa L.) at early stage of growth to dry lands

Reactive oxygen species and antioxidant machinery in abiotic stress tolerance in crop plants

The cDNA structures and expression profile of the ascorbate peroxidase gene family during drought stress in wild watermelon

Genomewide cloning and sequence analysis of leucine-rich repeat receptor-like protein kinase genes in Arabidopsis thaliana

The ethylene signaling pathway: new insights

From waste products to ecochemicals: fifty years research of plant secondary metabolism

Leaf senescense: correlated with increased levels of membrane permeability and lipid peroxidation, and decreased levels of superoxide dismutase and catalase

Overexpression of a stressresponsive NAC transcription factor gene ONAC022 improves drought and salt tolerance in rice

The antagonistic action of abscisic acid and cytokinin signaling mediates drought stress response in Arabidopsis

Transcriptome analysis and weighted gene co-expression network reveals potential genes responses to heat stress in turbot Scophthalmus maximus

Regulation of drought tolerance by gene manipulation of 9-cisepoxycarotenoid dioxygenase

Reactive oxygen species in abiotic stress signaling

WRKY transcription factors in plant responses to stresses

The novel rose MYB transcription factor RhMYB96 enhances salt tolerance in transgenic Arabidopsis

RhNAC3, a stress-associated NAC transcription factor, has a role in dehydration tolerance through regulating osmotic stress-related genes in rose petals

Activated expression of WRKY57 confers drought tolerance in Arabidopsis

Regulation of ascorbate peroxidase at the transcript level is involved in tolerance to postharvest water deficit stress in the cut rose

Transcription factors and plants response to drought stress: current understanding and future directions

Enzymes of ethylene biosynthesis

HISAT: a fast spliced aligner with low memory requirements

Identification of superior reference genes for data normalisation of expression studies via quantitative PCR in hybrid roses (Rosa hybrida)

The impact of water deficiency on leaf cuticle lipids of Arabidopsis

Calcium signals: the lead currency of plant information processing

Genome-wide analysis of GDSL-type esterases/lipases in Arabidopsis

WGCNA: an R package for weighted correlation network analysis

Drought stress provokes the down-regulation of methionine and ethylene biosynthesis pathways in Medicago truncatula roots and nodules

Functional analysis of Arabidopsis NCED6 and NCED9 genes indicates that ABA synthesized in the endosperm is involved in the induction of seed dormancy

Overexpression of the trehalose-6-phosphate synthase gene OsTPS1 enhances abiotic stress tolerance in rice

Cloning and expression analysis of the FvNCED3 gene and its promoter from ash (Fraxinus velutina)

Comparative transcriptome analysis of Rosa chinensis 'Slaters' crimson China' provides insights into the crucial factors and signaling pathways in heat stress response

RcMYBPA2 of Rosa chinensis functions in proanthocyanidin biosynthesis and enhances abiotic stress tolerance in transgenic tobacco

Physiological characteristic changes and full-length transcriptome of rose (Rosa chinensis) roots and leaves in response to drought stress

Whole-genome characterization of Rosa chinensis AP2/ERF transcription factors and analysis of negative regulator RcDREB2B in Arabidopsis

Dual function of NAC072 in ABF3-mediated ABA-responsive gene regulation in Arabidopsis

SlWRKY81 reduces drought tolerance by attenuating proline biosynthesis in tomato

Correction to: OsbZIP71, a bZIP transcription factor, confers salinity and drought tolerance in rice

Caragana korshinskii phenylalanine ammonialyase is up-regulated in the phenylpropanoid biosynthesis pathway in response to drought stress

Global transcriptome and weighted gene co-expression network analyses of growthstage-specific drought stress responses in maize

bHLH122 is important for drought and osmotic stress resistance in Arabidopsis and in the repression of ABA catabolism

Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2

Comparative transcriptomics and weighted gene co-expression correlation network analysis (WGCNA) reveal potential regulation mechanism of carotenoid accumulation in Chrysanthemum x morifolium

Arabidopsis CER8 encodes LONG-CHAIN ACYL-COA SYNTHETASE 1 (LACS1) that has overlapping functions with LACS2 in plant wax and cutin synthesis

Novel NAC transcription factor TaNAC67 confers enhanced multi-abiotic stress tolerances in Arabidopsis

Handling calcium signaling: Arabidopsis CaMs and CMLs

An effect of water stress on ethylene production by intact cotton petioles

Effects of flooding and drought on ethylene metabolism, titratable acidity and fruiting of pineapple

Mordor Intelligence (2020). Netherlands Floriculture Market-Growth, Trends, COVID-19 Impact, and Forecasts

Overexpression of calmodulin-like (ShCML44) stress-responsive gene from Solanum habrochaites enhances tolerance to multiple abiotic stresses

Enhancement of oxidative and drought tolerance in Arabidopsis by overaccumulation of antioxidant flavonoids

Identification of candidate genes involved in wax deposition in Poa pratensis by RNA-seq

The compensation effects of physiology and yield in cotton after drought stress

Drought stress induces oxidative stress and the antioxidant defense system in ascorbatedeficient vtc1 mutants of Arabidopsis thaliana

Cross-talk of brassinosteroid signaling in controlling growth and stress responses

Selective autophagy of BES1 Mediated by DSK2 balances plant growth and survival

Brassinosteroids: multidimensional regulators of plant growth, development, and stress responses

KEGG: kyoto encyclopedia of genes and genomes

An ethylene response factor (ERF5) promoting adaptation to drought and salt tolerance in tomato

Cold and water stresses produce changes in endogenous jasmonates in two populations of Pinus pinaster Ait

Regulatory Function of Ethylene in Plant Responses to Drought, Cold, and Salt Stresses

Transcriptlevel expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown

StringTie enables improved reconstruction of a transcriptome from RNA-seq reads

MfbHLH38, a Myrothamnus flabellifolia bHLH transcription factor, confers tolerance to drought and salinity stresses in Arabidopsis

The rice transcription factor OsWRKY47 is a positive regulator of the response to water deficit stress

Plant calmodulins and calmodulinrelated proteins: multifaceted relays to decode calcium signals

The Rosa genome provides new insights into the domestication of modern roses

Coping with stresses: roles of calcium-and calcium/calmodulin-regulated gene expression

ABO3, a WRKY transcription factor, mediates plant responses to abscisic acid and drought tolerance in Arabidopsis

edgeR: a bioconductor package for differential expression analysis of digital gene expression data

De novo assembly of Agave sisalana transcriptome in response to drought stress provides insight into the tolerance mechanisms

Specific oxidative cleavage of carotenoids by VP14 of maize

Ectopic expression of SOD and APX genes in Arabidopsis alters metabolic pools and genes related to secondary cell wall cellulose biosynthesis and improve salt tolerance

Cytoscape: a software environment for integrated models of biomolecular interaction networks

Comprehensive transcriptional analysis reveals salt stress-regulated key pathways, hub genes and time-specific responsive gene categories in common bermudagrass (Cynodon dactylon (L.) Pers.) roots

Regulatory network of gene expression in the drought and cold stress responses

Differential gene regulatory pathways and co-expression networks associated with fire blight infection in apple (Malus × domestica)

A Novel MYBS3-Dependent Pathway Confers Cold Tolerance in Rice

Overexpression of OsNAC14 improves drought tolerance in rice

Exogenous ABA induces drought tolerance in upland rice: the role of chloroplast and ABA biosynthesis-related gene expression on photosystem II during PEG stress

Isolation and functional analysis of Arabidopsis stress-inducible NAC transcription factors that bind to a drought-responsive cis-element in the early responsive to dehydration stress 1 promoter

Understanding regulatory networks and engineering for enhanced drought tolerance in plants

The protein phosphatases and protein kinases of Arabidopsis thaliana

Abscisic acid signaling inhibits brassinosteroid signaling through dampening the dephosphorylation of BIN2 by ABI1 and ABI2

ERF109 of trifoliate orange (Poncirus trifoliata (L.) Raf.) contributes to cold tolerance by directly regulating expression of Prx1 involved in antioxidative process

The cotton WRKY transcription factor (GhWRKY33) reduces transgenic Arabidopsis resistance to drought stress

Jasmonates: biosynthesis, perception, signal transduction and action in plant stress response, growth and development. an update to the 2007 review in annals of botany

Gm WRKY 54 improves drought tolerance through activating genes in ABA and Ca 2+ signaling pathways in transgenic soybean

Organ fusion and defective cuticle function in a lacs1lacs2 double mutant of Arabidopsis

Characterization of stress-responsive CIPK genes in rice for stress tolerance improvement

Cell signaling during cold. Drought, and Salt Stress

AtCPK6, a functionally redundant and positive regulator involved in salt/drought stress tolerance in Arabidopsis

Transcriptome sequencing and whole genome expression profiling of chrysanthemum under dehydration stress

The ethylene response factor OsERF109 negatively affects ethylene biosynthesis and drought tolerance in rice

Modulation of antioxidant defense system is associated with combined drought and heat stress tolerance in citrus

Analysis of trehalose-6-phosphate synthase (TPS) gene family suggests the formation of TPS complexes in rice

A general framework for weighted gene co-expression network analysis

Overexpression of VaWRKY14 increases drought tolerance in Arabidopsis by modulating the expression of stress-related genes

The novel wheat transcription factor TaNAC47 enhances multiple abiotic stress tolerances in transgenic plants

RcLEA, a late embryogenesis abundant protein gene isolated from Rosa chinensis, confers tolerance to Escherichia coli and Arabidopsis thaliana and stabilizes enzyme activity under diverse stresses

Physiological and transcriptomic analysis of tree peony (Paeonia section Moutan DC.) in response to drought stress

Functional characterization of poplar WRKY75 in salt and osmotic tolerance

Physiological and transcriptomic analysis revealed the involvement of crucial factors in heat stress response of Rhododendron hainanense

Transcriptome sequencing of dianthus spiculifolius and analysis of the genes involved in responses to combined cold and drought stress

Physiological analysis and transcriptome sequencing reveal the effects of combined cold and drought on tomato leaf

Salt and drought stress signal transduction in plants

Abiotic stress signaling and responses in plants