Changes in phenotype and physiology at drought stress in R.chinensis
We applied drought stress by stopping watering on R.chinensis. The plant leaf phenotypes were visually inspected and physiologically measured under control condition (CK), and at 15,30,45,60,75,90 days after stress and rewatering was applied. Under normal circumstances, the leaves were dark green and shiny (Fig. 1a). After 15 days of drought treatment, the leaves had no obvious morphological changes (Fig. 1b). The leaf relative water content is a quantitative indicator of plant water status [17]. Similarly, the relative soil water content and leaf relative water content were not significantly different from the control (Fig. 2a, Fig. 2b). At 30 days the leaves still appeared healthy, but the leaves were slightly wrinkled (Fig. 1c). When the drought stress lasted for 45 days and 60 days, the leaves were all curled (Fig. 1d, Fig. 1e). At this time, the relative water content of the soil decreased significantly, reaching a moderate degree of stress (Fig. 2a). At 60 days, the relative water content of the leaves decreased significantly compared with the control, which was consistent with the phenotype (Fig. 2b). The curl degree of the leaf deepened and the leaf became thinner, and the relative water content of the leaf was significantly reduced during 75 days of drought (Fig. 1f, Fig. 2b). However, after 90 days of drought, both the leaf relative water content and the soil relative water content reached the lowest, reaching the level of severe stress. At the same time, the leaves were severely curled and withered (Fig. 1g, Fig. 2a, Fig. 2b). After rewatering, the leaves were green and thin (Fig. 1h). In order to determine the exact point for transcriptome sequencing, we also measured physiological indicators. Malondialdehyde (MDA) is usually an important physiological indicator to measure oxidative damage under abiotic and biological stress. The MDA content of leaves treated with drought and rewatering was determined. With the deepening of the degree of drought stress, the content of MDA rises in fluctuation and got the highest after 90 days of the drought treatment (Fig. 2d). In addition, the content of superoxide dismutase (SOD) was examined and it was found that SOD increased with the increasing of drought treatment time and reached the maximum value after 90 days (Fig. 2c). We also found that both MDA and SOD content increased significantly during drought treatment for 30 days, 60 days and 90 days. After rewatering, the MDA and SOD contents were all reduced. Based on these results, we selected leaves under normal conditions, drought treatment for 30 days, 60 days, 90 days, and rewatering treatment for transcriptome sequencing.
Transcriptome analysis of R.chinensis
To investigate genes involved in the transcriptome, a cDNA sample was prepared from an equal mixture of total RNA isolated from the leaves of 2-year-old cutting seedlings R.chinensis from the control, drought-treated (30,60,90 d) 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. The clean reads removing adaptor sequences, ambiguous nucleotides, and low-quality reads (Q value < 20) 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–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–96.38%, and the proportion of intron sequences ranged from 1.45–2.66% (Additional file 1: figure S4).
Table 1
RNA sequencing data and corresponding quality control
Sample name
|
Clean reads
|
Raw reads
|
Mapped reads
|
GC Content (%)
|
Q20(%)
|
Q30(%)
|
CK1
|
48126136
|
48181768
|
44510062 (93.27%)
|
48.28
|
97.51
|
92.92
|
CK2
|
49106032
|
49156710
|
45429755 (93.89%)
|
48.66
|
98.03
|
94.17
|
CK3
|
43793126
|
43853234
|
40291282 (93.01%)
|
48.29
|
97.48
|
92.88
|
DT1-1
|
45651270
|
45706946
|
42334367 (93.31%)
|
49.52
|
97.67
|
93.28
|
DT1-2
|
48888834
|
48972378
|
44986749 (92.70%)
|
48.11
|
97.38
|
92.70
|
DT1-3
|
53820314
|
53902378
|
49546849 (92.94%)
|
48.30
|
97.46
|
92.85
|
DT2-1
|
47463436
|
47521292
|
43862640 (93.07%)
|
48.19
|
98.12
|
94.31
|
DT2-2
|
54353944
|
54420772
|
49969726 (93.02%)
|
48.23
|
98.10
|
94.27
|
DT2-3
|
51937470
|
52002482
|
47795161 (92.33%)
|
48.09
|
98.13
|
94.32
|
DT3-1
|
50304612
|
50359618
|
40931525 (93.07%)
|
49.61
|
98.27
|
94.71
|
DT3-2
|
44365148
|
44421624
|
38029228 (93.35%)
|
49.58
|
98.10
|
94.28
|
DT3-3
|
49072308
|
49130228
|
41525793 (93.40%)
|
49.60
|
98.18
|
94.50
|
RW1
|
54756474
|
54838556
|
50746877 (93.11%)
|
48.65
|
97.99
|
94.00
|
RW2
|
41397356
|
41446080
|
38442300 (93.52%)
|
48.94
|
97.91
|
93.88
|
RW3
|
48571772
|
48639152
|
45164692 (93.24%)
|
48.56
|
98.16
|
94.46
|
Comparisons of differentially expressed genes under the different drought stress stages
To analyze the differentially expressed genes (DEGs) between leaves suffered drought stress at four different time duration (30,60,90 d and rewatering) and the control, we conducted comparisons between five groups. The DEGs were filtrated according to an expression level |log2(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 (Fig. 3a).
The number of DEGs is an increase from DT1 to DT3, reached a peak at CK vs. DT3 comparison set, indicated that most of the drought regulated genes are late-response. The results showed DEG enrichment in CK vs. DT1, CK vs. DT2 and CK vs. DT3, illustrating that 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 is the most among the drought stress groups. Meanwhile, only 2503 upregulated and 2367 downregulated between CK and DT1, which suggested that the differentiation of expressed genes between CK and DT3 is larger than that of CK and DT1, this revealed that, transcript abundance changed dramatically at key switched among the drought stages, which drought response genes could be induced and expressed largely. Comparing with CK, slightly more genes upregulated than genes downregulated in RW (Fig. 3b).
GO functional enrichment and KEGG pathway enrichment analysis of DEGs
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’, so we can see, the drought resistance process in R.chinensis is complex. Meanwhile, according to 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) (Fig. 4a, Fig. 4b, Fig. 5b). And CK vs DT3 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.) (Fig. 5a, Additional file 2: Table S1). In addition, the GO analysis of the down-regulated DEGs found that among the top 20 enriched, 13 belong to the category of biological processes, mainly focusing on the response to stimuli in CK vs DT1(Fig. 4a, Additional file 3: Table S2).
The KEGG pathway analysis revealed that it is mainly enriched in metabolism under the different drought degree 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, which suggested that drought stress had been sensed and signaling transduction pathways had been activated (Fig. 6a). 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 (Fig. 6b, Fig. 7a). 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 metabolism processes may act as the predominant process of middle and late stages of drought stress in R.chinensis. It was unsurprising because drought causes changes in the metabolism of plants to improve their drought tolerance. Apart from pathways enriched metabolism, ‘plant-pathogen interaction’(ko04626) was also enriched in CK vs RW, indicating that plants are vulnerable to pathogenic bacteria after recovering from drought stress (Fig. 7b, Additional file 4: Table S3).
Identification of plant hormone and signal transduction-related DEGs
Signal transduction in plants was evidently influenced by drought stress, as verified by both the GO and KEGG enrichment analyses of the DEGs performed in this study. 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 was found to show the highest degree of upregulation in DT3. Mitogen-activated protein kinase (MAPK) occupy an important role in signal transduction induced by abiotic stress [18]. The KEGG pathways analysis of the DEGs, the ‘MAPK signaling pathway’(ko04016) was induced in the drought treatment. For DEGs detected in DT1 and DT2, two MAPKKKs including RchiOBHm_Chr2g0158821 and RchiOBHm_Chr5g0002581 and four MAPKKs including, MPK9(RchiOBHm_Chr7g0178621), MKK10(RchiOBHm_Chr7g0211301), MEKK1(RchiOBHm_Chr7g0240941), MPK4(RchiOBHm_Chr7g0196781) were up-regulated. MKK10(RchiOBHm_Chr7g0211301) was particularly up-regulated 9.18-fold. However, three MAPKKs consist of MPK3(RchiOBHm_Chr5g0061451), MKK4(RchiOBHm_Chr2g0150341) and MEKK1(RchiOBHm_Chr1g0352161) were down-regulated. In addition, MPK4(RchiOBHm_Chr7g0196781) and MPK19(RchiOBHm_Chr5g0028971) were substantially increased in DT3(Additional file 5: Table S4).
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 [19, 20]. Ca2+ signal transmission is achieved through the generation, decryption and transmission of specific calcium signals and the corresponding physiological and biochemical reactions downstream. In this study, we identified DEGs related to calcium ion, encoding calcium-binding proteins, Ca2+-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 (Additional file 6: Table S5).
When plants were subjected to environmental stress, reactive oxygen species generated signals in the plant, triggering a series of changes. Plants also have an antioxidant system protecting against poisonous oxygen [21]. According to GO and KEGG analysis of the DEGs, we found that the GO term ‘hydrogen peroxide metabolic’(GO:0042743) and ‘response to reactive oxygen species’(GO:0000302) were enriched at DT1 and DT2, Peroxisome(ko04146) was induced at DT3. We identified 69 DEGs encoding 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 peroxidase was up-regulated by 4.79-fold at DT1, another gene encoding peroxidase PNC1(RchiOBHm_Chr5g0072281) was no significant change from DT1 to the control, but it increased significantly 9.29-fold at DT2, reaching the peak of expression level at DT3 (Additional file 7: Table S6).
Moreover, the expression of genes involved in hormone biosynthesis or signaling was substantially changed following the drought stress and then rewatering, including ABA, auxin, brassinosteroid (BR), cytokinin, ethylene (ET), jasmonic acid (JA) and salicylic acid (SA). Among these DEGs, most genes were involved in ABA response. Nine-cis-epoxycarotenoid dioxygenase (NCED) is an important enzyme in ABA biosynthesis. For instance, one gene encoding NCED3 (RchiOBHm_Chr5g0014331) was up-regulated by 2.70-fold at DT1. Another gene encoding NCED6(RchiOBHm_Chr4g0397001) was significantly down-regulated 9.40-fold at DT3. Interestingly, the expression of a gene (RchiOBHm_Chr5g0049951) encoding AAO, which is also a vital enzyme in ABA biosynthesis, was up-regulated by 2.78-fold at DT1 and down-regulated by 7.28-fold at DT3. Moreover, some genes involved in 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, was up-regulated (Additional file 8: Table S7).
In total of 24 genes involved in the auxin signaling pathway were also changed. Three genes encoding auxin response factor (ARF3, ARF5, ARF9) were drought induced. Particularly, ARF3 (RchiOBHm_Chr5g0009381) expression continued to decline under drought stress, which was induced for 4.41-fold at DT3, while 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 downed by nearly 10-fold at DT3. In addition, some genes participated in auxin response and signal, such as auxin-binding protein ABP19a-like (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 (Additional file 8: Table S7).
In this study, large numbers of ethylene-responsive transcription factor (ERF) genes displayed a trend of down-regulation under early drought stress. However, one gene encoding ERF023 (RchiOBHm_Chr1g0360021) was up-regulated 3.70-fold at DT1 and then down-regulated 5.02-fold at DT3. Another gene encoding ERF110 (RchiOBHm_Chr6g0274591) was up-regulated 7.34-fold at DT3. Meanwhile, we identified ethylene receptor (ERS1), ethylene-overproduction protein 1 (ETO1) and ethylene insensitive3-like protein (EIL3) (Additional file 8: Table S7).
Most of 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. Furthermore, the majorities of genes related to cytokinin, salicylic acid and brassinosteroid pathways. One gene encoding cytokinin dehydrogenase 3 (RchiOBHm_Chr2g0093491) was down-regulated 7.40-fold at DT3 (Additional file 8: Table S7).
Expression of genes involved in metabolism and biosynthesis
Through the GO and KEGG enrichment analysis, we found 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. Except for a small portion of starch synthase genes that were down-regulated, all were 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, pyrophosphate--fructose 6-phosphate 1-phosphotransferase subunit beta-like and 6-phosphofructo-2-kinase/fructose-2,6-bisphosphatase were down-regulated at DT3 (Additional file 9: Table S8).
In addition, we found that 89 DEGs associated with lipid metabolism were enriched in the ‘lipid metabolic process’(GO:0006629), most of them related to the biosynthesis and metabolic processes of lipoprotein, wax, keratin, lipid, and fatty acid. 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(Additional file 10: Table S9). Moreover, KEGG functional annotation of these differentially expressed genes related to metabolism. Many differentially expressed genes were related to secondary metabolisms. The total of 95 DEGs involved in secondary metabolism were enriched by KEGG function, and the results were shown in Fig.8. Drought stress affected the transcription level of genes related to phenylpropanoid biosynthesis pathways.
Photosynthesis-related genes under drought stress
We observed that many DEGs enriched in the category of ‘photosynthesis’(ko00915) were overexpressed in DEGs down-regulated by DT2 and DT3, and 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. Among them, three (CAB40, LHCB5, psaA) of 16 DEGs encoding the components of photosystem I (PSI) were slightly up-regulated in DT1, and the other genes were significantly down-regulated during the drought treatment period. A total of 33 genes encoding PSⅡ were identified, including 6 up-regulated genes (psbE, 3genes encoding psbC, psbD, PSB27-1) and 27 down-regulated genes. Most of these genes decreased significantly in the middle and late periods of drought, indicating that PSI and PSⅡ 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 6 genes encoding components of chloroplasts were significantly decreased (Additional file 11: Table S10).
Transcription factors responding to drought stress
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 down-regulate [22]. 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 showed large changes under different treatments. The expression of some genes reached the maximum in the middle and late periods of drought. Genes encoding WRKY75 and WRKY71 were all reached the highest expression level at DT3, up-regulated more than 5-fold. WRKY57 (RchiOBHm_Chr5g0083891) was up-regulated regarding to expression responses at different drought stress, especially increased 7.93-fold at DT3. However, over half of the WRKY family has decreased expression levels at different drought levels. Twenty-four of genes encoding WRKYs were down-regulated at DT1 by 1.02- to 3.71-fold (Fig.9c, Additional file 12: Table S11).
More than half of the genes encoding NACs are up-regulated to varying degrees under drought stress. One gene encoding NAC100 (RchiOBHm_Chr2g0167461) continued to increase in the middle and late periods of drought stress, especially by 7.20-fold. A total of 11 genes upregulated during DT1. Some genes increased and then decreased under drought stress, but 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 it was increased by 6.96-fold after rewatering. NAC037 (RchiOBHm_Chr2g0099671) has a similar expression trend, but the up-regulation decrease by a small amount. On the contrary, one gene encoding NAC072 (RchiOBHm_Chr5g0034761) was down-regulated at DT1 and continued to increase afterwards (Fig.9a, Additional file 12: Table S11.).
In addition, under different drought treatment, the expression level of more than three-quarters of the MYBs genes showed different degrees of up-regulation. Among them, the number of up-regulated genes that appeared during the DT1 period was relatively large. For instance, MYB46 (RchiOBHm_Chr1g0315931) and MYB61 (RchiOBHm_Chr1g0315931) showed the highest up-regulation of expression with 7.09 and 7.15-fold. Moreover, some genes were up-regulated significantly at DT3. MYB75 (RchiOBHm_Chr2g0116041) was the most up-regulated among all up-regulated genes at DT3, increased by 7.65-fold. Some genes showed a tendency to increase and then decrease with drought treatment, such as gene encoding MYB61 (RchiOBHm_Chr3g0458721), slightly increased at DT1 and then decreased by 6.10-fold. After rewatering, a small amount of gene expression was up-regulated, and most of the gene expression did not change significantly, comparing with CK (Fig.9d, Additional file 12: Table S11.).
In the three drought treatments of the bHLH families, the number of genes up-regulated in one period was more than the number of genes down-regulated, and the number of down-regulated in two periods was more than that. Most genes of the bHLH family were up-regulated during the early drought period, DT1, while large numbers of genes were significantly down-regulated during the middle and late drought periods (DT2 and DT3). Among these, there was a slight upward adjustment, up-regulated by 1.02- to 4.66-fold at DT1, but there were 37 genes down-regulated at DT2 or DT3. One gene encoding bHLH72 (RchiOBHm_Chr7g0182341) was down-regulated 6.94-fold. A gene bHLH121 (RchiOBHm_Chr6g0245181) belonging to the bHLH family was no significant change in gene expression during drought treatment, but it was up-regulated after rewatering (Fig.9b; Additional file 12: Table S11.).
Among AP2/DREB and bZIP transcription factors, more genes were down-regulated under drought treatment. Nineteen members of bZIPs were found to be responsive to drought stress and 12 were down-regulated, from 1.24- to 4.80-fold (Fig.9e; Additional file 12: Table S11.). Five genes encoding DREBs, all down-regulated exceeded 6.60-fold at DT1. The most significant down-regulation was that the DREB1D (RchiOBHm_Chr7g0199351) was down-regulated 9.64-fold (Fig.9f; Additional file 12: Table S11.).
Co-expression network construction and prediction of regulatory circuitry controlling genes in drought stress response
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 and different colors represented different modules. A total of 18 modules were identified (Fig.10a). The module gene expression pattern in each sample was displayed with the module characteristic value (Fig.10c). Transcription factors are an important class of regulatory proteins in biological processes. This study analyzed the transcription factors in each module. Results showed that the number of transcription factors in the module was 2 (salmon 2)-286 (coral 2) (Fig.10b). There were differences in the distribution of transcription factors in different modules, but mainly concentrated in WRKY, NAC, MYB, ERF, bHLH, bZIP, C2H2 and other transcription factor families. According to previous reports, these transcription factors were involved in the process of plant stress regulation [23]. Furthermore, the module-trait correlation relationships are shown in Fig.11. Four interesting modules, coral 2, lightsteelblue 1, navajowhite 1 and bisque 4, were selected according to the criterion of |r| > 0.5 and P < 0.05. We analyzed the expression pattern of the co-expressed genes in these four selected modules (Fig.12). Meanwhile, further annotated these four modules with GO and KEGG analysis (Additional file 1: Figure S2; Additional file 1: Figure S3; Additional file 1: Figure S4). We found that in the navajowhite 1 module, the significantly rich GO terms indicated the key biological processes of co-expressed gens, and many genes were enriched in various stress responses. 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 enrich (Fig.13a). KEGG enrichment analysis showed that the main enrichment in metabolic pathways (Fig.13b). Therefore, we speculated that this module genes may play an important role in the process of rose drought stress and rewatering.
Weighted gene co-expression network analysis (WGCNA) was used to identify putative transcription factors. The transcription factors of the navajowhite1 module were analyzed, and screened out the highly connected genes, which were used as the core genes. Highly connected nodes in expression networks were defined as hub genes. Many DEGs Identified in the navajowhite1 module were annotated as TFs, including several WRKY, MYB, NAC, ERF, ARF and bHLH TFs (Fig.14a). Based on the hub gene correlation network and their high connectivity, 22 of these TFs were particularly selected: two genes (RchiOBHm_Chr5g0034761 and RchiOBHm_Chr6g0308271) encoding NAC072, two genes (RchiOBHm_Chr2g0167461 and RchiOBHm_Chr7g0181811) encoding NAC100, NAC087 (RchiOBHm_Chr3g0497471), MYB330(RchiOBHm_Chr7g0228621), MYB75 (RchiOBHm_Chr2g0116041), MYB102 (RchiOBHm_Chr5g0006031), MYB78 (RchiOBHm_Chr4g0426181), WRKY3 (RchiOBHm_Chr4g0439041), WRKY4 (RchiOBHm_Chr2g0166991), WRKY6 (RchiOBHm_Chr5g0040801), WRKY27 (RchiOBHm_Chr1g0378621), WRKY28 (RchiOBHm_Chr2g0151681), WRKY71 (RchiOBHm_Chr5g0013131), two genes (RchiOBHm_Chr4g0429851 and RchiOBHm_Chr2g0175911) encoding WRKY75, ERF113 (RchiOBHm_Chr4g0428891), ARF1 (RchiOBHm_Chr7g0219771), ARF5 (RchiOBHm_Chr6g0302551), ARF8 (RchiOBHm_Chr3g0487771) and the bHLH3 (RchiOBHm_Chr6g0288981) (Fig.14b).
Validation of the DEGs by qRT-PCR
To verify the accuracy and reproducibility of the transcriptome analysis, 13 selected DEGs were analyzed the transcript abundance using qRT-PCR, including 7 randomly selected transcription factors from navajowhite1 module, another 4 transcription factors and 2 ABA synthesis-related DEGs. The results showed that the expression profiles detected by qRT-PCR were positively correlated with the RNA-Seq results (Fig.15). Therefore, the reliability of our RNA-seq data was confirmed by the consistency between the qRT-PCR results and RNA-seq analyses.