Skip to main content

Transcriptome profiling reveals the effects of drought tolerance in Giant Juncao



Giant Juncao is often used as feed for livestock because of its huge biomass. However, drought stress reduces forage production by affecting the normal growth and development of plants. Therefore, investigating the molecular mechanisms of drought tolerance will provide important information for the improvement of drought tolerance in this grass.


A total of 144.96 Gb of clean data was generated and assembled into 144,806 transcripts and 93,907 unigenes. After 7 and 14 days of drought stress, a total of 16,726 and 46,492 differentially expressed genes (DEGs) were observed, respectively. Compared with normal irrigation, 16,247, 23,503, and 11,598 DEGs were observed in 1, 5, and 9 days following rehydration, respectively. Gene Ontology and Kyoto Encyclopedia of Genes and Genomes pathway analyses revealed abiotic stress-responsive genes and pathways related to catalytic activity, methyltransferase activity, transferase activity, and superoxide metabolic process. We also identified transcription factors belonging to several families, including basic helix-loop-helix (bHLH), WRKY, NAM (no apical meristem), ATAF1/2 and CUC2 (cup-shaped cotyledon) (NAC), fatty acyl-CoA reductase (FAR1), B3, myeloblastosis (MYB)-related, and basic leucine zipper (bZIP) families, which are important drought-rehydration-responsive proteins. Weighted gene co-expression network analysis was also used to analyze the RNA-seq data to predict the interrelationship between genes. Twenty modules were obtained, and four of these modules may be involved in photosynthesis and plant hormone signal transduction that respond to drought and rehydration conditions.


Our research is the first to provide a more comprehensive understanding of DEGs involved in drought stress at the transcriptome level in Giant Juncao with different drought and recovery conditions. These results may reveal insights into the molecular mechanisms of drought tolerance in Giant Juncao and provide diverse genetic resources involved in drought tolerance research.


Drought stress is one of the most threatening environmental constraints that adversely affect plant growth and yield [1]. However, with global climate change, the frequency and intensity of drought have continuously increased [2, 3]. Drought might cause metabolic imbalance in plant cells and influence the optical energy absorption of plant leaves, destroying the photosynthetic organs of plants [4]. Moreover, drought can lead to the accumulation of active oxygen substances in the leaves, which may accelerate the peroxidation of biological membrane lipid to produce toxic products, thereby inhibiting plant growth [5]. For plants exposed to drought stress, the total primary productivity is not only closely related to their resistance and tolerance to drought stress, but also shows an important relationship with the ability of plants to recover from damage after the elimination of the stress [6, 7]. Therefore, the recovery ability after rehydration is important for the successful adaptation of plants to arid environments. Rehydration helps plants recover their physiological functions, and it can offset plant damage from drought stress to a certain extent [7]. However, the compensation of rehydration to plant growth after drought stress is often limited. The recovery degree of plant growth might be related to the degree and duration of drought stress before rehydration and drought resistance of plants [8].

To cope with drought stress, plants have adapted various self-protection and defense mechanisms in the long-term evolution process [9]. Plants decrease the photosynthetic capacity of mesophyll by closing the stomata to adjust the photosynthetic process of leaves [8]. By adjusting the in vivo antioxidation system, plants can eliminate excess active oxygen and maintain in vivo redox equilibrium. The intracellular water potential can be increased by increasing the substances of cell osmotic adjustment to maintain a certain expansion, thereby protecting the continuous growth of plants under drought stress [10]. Plant hormones regulate their own response mechanism through synergistic or antagonistic action in response to arid states. Various expressed genes have been reported in response to drought stress [11]. These genes include stabilizing membrane proteins, heat shock proteins and late embryogenesis abundant proteins, which play an important role in stabilizing protein structure and enhancing cell’s water binding capacity. Early drought-induced proteins protect plants by producing certain metabolic proteins and regulating gene expression through precise signal transduction during drought stress. Dehydrin genes (Dhn), which are among the most frequently observed proteins in plants, protect the cells from water deficit. Additionally, a large number of genes change their expression through the regulation of transcription factors (TFs) [12]. Several TFs also provide response under drought stress, including abscisic acid-responsive element (ABREs) binding factors (ABFs, AREBs, or DPBFs) [13,14,15], dehydration-responsive element binding factors [16], myeloblastosis (MYB), and SNF1-related kinase 2 [17]. Therefore, knowledge about the various genes translated and expressed in response to drought stress conditions will help elucidate the water deficit tolerance mechanisms and facilitate the development of new plant cultivation tools to combat climate change.

Giant Juncao is an ideal Gramineae C4 plant for water-soil conservation, wind prevention, and sand fixation due to its large biomass, disease resistance, and developed root system. This plant has been widely applied for the comprehensive environmental management of regions with vulnerable ecology [18]. However, drought environment influences the yield and restricts its large-scale plantation. To improve its productivity and performance in water deficit conditions, the response mechanisms of Giant Juncao to drought should be elucidated. Recently, the development of molecular biological methods has facilitated the discovery of potential plant response mechanisms to environmental stresses. RNA sequencing (RNA-seq) can quickly and comprehensively obtain the gene expression of a specific cell or tissue in a certain state, so as to determine the molecular mechanism of physiological metabolic responses of plants under abiotic stress conditions [19]. RNA-seq data could provide insights into the discovery of new genes, including annotation genes and differentially expressed genes (DEGs), and molecular markers [20]. Compared with traditional sequencing methods, RNA-seq provides high-throughput sequencing results, is inexpensive, has high sensitivity, and can detect low abundance expressed genes [21]. A large quantity of DEGs associated with drought stress response have been reported in various Gramineae species by RNA-seq, such as wheat [22], maize [23], rice [24], sorghum [25], and foxtail millet [26]. This method can obtain transcripts from different developmental stages, tissues, and organs, so it is a fundamental and efficient method for discovering functional genes.

Herein, we aimed to identify the genes involved in the response of Giant Juncao to drought stress and rehydration treatment using Illumina sequencing. We analyzed DEGs in Giant Juncao seedlings to elucidate its molecular response mechanism to drought stress. Furthermore, we employed the physiological indices to understand the drought response mechanism of Giant Juncao. As only a few studies have identified genes that respond to drought stress in Pennisetum spp., our study provides an important transcriptomic database for further targeted gene modifications in grasses.


Physiological evaluation of Giant Juncao seedling in response to drought and rehydration treatment

Drought stress and rehydration treatments affected the physiological index of Giant Juncao (Fig. 1a). Compared with the control, drought-stressed plants showed a significant decrease in Pn, Tr, and Gs at D1 and D2, whereas the WUE increased at D1 and decreased at D2 (Fig. 1b). After rehydration, all ecological indices except Gs significantly increased compared with those after 14 days of drought stress. With the increase in recovery time, the ecological index of Giant Juncao reached a stable stage during 5 and 9 days of rehydration, and the Gs steadily increased. The Pn was 14.43–15.40 μmol m− 2 s− 1 for nearly 7 days of drought stress. The Tr and Gs were 2.84–3.19 mol m− 2 s− 1 and 55.50–69.50 μmol m− 2 s− 1, respectively, which were higher than those in plants exposed to 7 days of water deficit. The WUE increased to 5.06 mmol mol− 1, which showed no significant difference with the control when the plant was re-watered for 5 days (Fig. 1b).

Fig. 1
figure 1

Effect of drought stress and rehydration condition on the growth of Giant Juncao. a A picture of Giant Juncao under different treatments condition. b physiological and ecological indicators between drought stress and rehydration treatment of Giant Juncao. Different lowercased letters in the same panel indicate statistical significance (P < 0.05). Pn: net photosynthetic rate, Tr: transpiration rate, Gs: stomatal conductance, WUE: water use efficiency, SOD: superoxide dismutase, POD: peroxidase, PRO: proline, MDA: malondialdehyde; P < 0.05

After drought stress, SOD, POD, PRO, and MDA showed an increasing and decreasing trend. Since the re-watering, the SOD and PRO rapidly decreased (Fig. 1b). The POD remained stable, and MDA significantly increased when watering was restored for 1 day. With the increasing rehydration time, POD, PRO, and MDA exhibited a continued declining trend, but SOD showed the opposite tendency. At the end of the experiment, SOD and MDA levels were 4.09 U mg− 1 protein and 17.16 nmol g− 1 FW, respectively, which were significantly higher than those under non drought stress (P < 0.05), whereas the POD and PRO were similar to those in the control.

Transcriptome sequencing and assembly

In this study, 18 cDNA libraries from the non drought stress (0 day), drought stress (7 and 14 days), and rehydration (1, 5, and 9 days) conditions were constructed totally, and then C0, D1, D2, R1, R2, and R3 were referred correspondingly. Each condition had three biological repetitions. An overview of the sequencing is listed in Table S1, and more than 97.18% and approximately 94% of the bases from the over 1 billion raw reads had a q-value ≥20 and ≥ 30 (an error probability of 0.02–0.025%), respectively. The GC-content was between 55.67 and 57.76% (Table S1). All original FASTQ data files were submitted to the NCBI Sequence Read Archive (SRA), under accession numbers is PRJNA632455.

After filtering out low-quality reads, a total of 967 million clean reads were generated from six samples (Table S1). Trinity was used to generate 144,806 transcripts with N50 of 1705 bp and N90 of 645 bp (Table 1). Among them, 93,907 were unigenes, where less than 300 bp has 3630 unigenes, 300–500 bp has 16,003 unigenes, 500–1000 bp has 35,098 unigenes, 1–2 kb has 25,770 unigenes, and > 2 kb has 13,406 unigenes.

Table 1 Summary statistics of Giant Giant Juncao transcriptome assemblies

Functional and pathway annotation of genes

In order to investigate the genes’ function, all the assembled unigenes were annotated into the Nr, Nt, Swiss-Prot, KOG, Pfam, GO, and KEGG pathway databases (see Materials and Methods). A total of 57,941 unigenes were annotated in at least one database, which account for 61.7% of the unigenes (Fig. 2a, Table S2). A total of 6963 (7.4%) unigenes were annotated in the public databases (Fig. 2a). The numbers of unigenes in Nr, Nt, KEGG, Swiss-Prot, Pfam, GO, and KOG databases were 46,952 (50.0%), 47,268 (50.3%), 14,221 (15.1%), 30,901 (32.9%), 40,136 (42.7%), 26,618 (28.4%), and 24,334 (25.9%), respectively (Fig. 2a, Table S2). The unigenes matched sequences from the genomes of Setaria italica, Sorghum bicolor, Zea mays, Oryza sativa Japonica Group, Oryza sativa Indica Group, Brachypodium distachyon, Anthurium amnicola, Aegilops tauschii, Oryza brachyantha, and others (Fig. 2b).

Fig. 2
figure 2

Information of unigenes annotated in different databases (a) and distributed into different species (b)

DEGs under drought and rehydration conditions

The gene expression levels were calculated and normalized using the TMM method; |log2(FoldChange)| > 1 and adjusted p-value < 0.05 were set as the threshold for significant differential expression. We found 62,298 unigenes which showed differential expression in drought-rehydration-treated samples. In total, 16,726 and 46,492 DEGs were detected at days 7 and 14 of drought treatment. A total of 7819 and 22,187 DEGs were upregulated, whereas 8907 and 24,305 DEGs were downregulated at D1 and D2, respectively (Fig. 3a). Compared with control plants, a total of 16,247 (9040 upregulated and 7207 downregulated), 23,503 (11,469 upregulated and 12,034 downregulated), and 11,598 (5048 upregulated and 6550 downregulated) DEGs were detected at R1, R2, and R3, respectively (Fig. 3a).

Fig. 3
figure 3

Up- and downregulated DEGs and venn diagrams showing the numbers of DEGs across five comparisons. a: DEGs were upregulated or downregulated by drought-rehydration treatment in Giant Juncao. D1, D2, R1, R2, and R3 show the comparisons of differential unigenes identified between 7 days of drought vs. control, 14 days of drought vs. control, 1 day of rehydration vs. control, 5 days of rehydration vs. control, and 9 days of rehydration vs. control, respectively; b: upregulated DEGs in venn diagrams; c: downregulated DEGs in venn diagrams

To identify common unigenes in different physiological stages, the overlaps in each comparison were shown in a Venn diagram (Fig. 3b and c). The highest overlap of the up- and downregulated genes was between 14 days of drought and control, which showed > 70% common unigenes, suggesting plant damage after strong drought (Figs. 3b and c). The second highest overlap of the up- and downregulated unigenes was between 5 days of rehydration, suggesting a gradual recovery after re-watering at R2 (Figs. 3b and c). A total of 140 upregulated unigenes (Fig. 3b) and 993 downregulated unigenes (Fig. 3c) overlapped with those of 7 and 14 days of drought and 1, 5, and 9 days of rehydration versus control, respectively, suggesting that a shared group of genes was associated with response to water deficit and rehydration treatment conditions.

GO enrichment in DEGs

To identify the major biological processes that are expressed at the drought-rehydration conditions, we showed each GO enrichment analysis of up- and downregulated DEGs at five physiological stages depending on FDR < 0.05 (Additional file 2). On the basis of the GO enrichment results, we selected the most significant accession to explain the physiological changes. Processes that promote activity were mainly enriched in the downregulated DEGs at D1, whereas only the carbohydrate metabolic process was enriched in the downregulated DEGs at D2 and after re-watering at R1. Drug binding and ATP binding processes were enriched in the downregulated DEGs at D2 and were further enriched at R2, suggesting reduced protein combining ability in response to drought. Ribosome and peptide metabolic processes were enriched in the upregulated DEGs at R1 and R2, whereas the structural constituent of the ribosome was enriched in the upregulated and downregulated DEGs at R3, suggesting the gradual completion of cell repair. DNA-binding TF activity and transcription regulator activity were enriched in the upregulated DEGs, whereas the cellular protein metabolic process was enriched in the downregulated DEGs at R3, suggesting improvement in TF activity and decreased protein degradation after re-watering (Table 2).

Table 2 GO analysis of the most significant up- and downregulated DEGs (FDR < 0.05)

KEGG pathway enrichment analysis

The KEGG pathway enrichment analysis showed that 1972 upregulated and 2337 downregulated DEGs were involved in 117 different pathways in Giant Juncao when comparing 7 days of drought with the control (Tables S3 and Additional file 3). With prolonged drought, 4490 upregulated and 5725 downregulated DEGs were annotated in 119 and 116 different pathways (Tables S3 and Additional file 3). After rehydration for 1, 5, and 9 days, 2166, 3106, and 1126 upregulated DEGs and 2096, 3229, and 2016 downregulated DEGs were annotated in the KEGG pathway database, which involved 111, 122, and 107 different upregulated pathways and 115, 115, and 114 different downregulated pathways (Tables S3 and Additional file 3). According to the KEGG pathway analyses, a heatmap reflecting the different biological processes in the significantly enriched pathways was constructed (Fig. 4). In these pathways, glutathione metabolism was enriched at D1 and R1 in the upregulated DEGs, whereas protein processing in the endoplasmic reticulum was enriched at D1 and in the downregulated DEGs at R2. Pathways related to homologous recombination and mismatch repair were enriched at D1, whereas alpha-linolenic acid metabolism, galactose metabolism, and regulation of autophagy were enriched at D2 in upregulated DEGs, suggesting improved metabolism in response to drought. Ribosome biogenesis in eukaryotes, pantothenate and CoA biosynthesis, and pyrimidine metabolism were enriched in the upregulated DEGs at R1, reflecting recovery after re-watering. The major of pathway related to ribosome was enriched in the upregulated DEGs at R1 and R2 and downregulated DEGs at D2 and R3. Plant pathogen interaction was enriched in the upregulated DEGs at R1 and R3 and downregulated DEGs at D1 and R2. The upregulated DEGs related to photosynthesis were enriched in porphyrin and chlorophyll metabolism, photosynthesis - antenna proteins, and carotenoid biosynthesis pathways at R2 and R3, suggesting that photosynthesis plays an important role after re-watering. Amino sugar and nucleotide sugar metabolism, phenylpropanoid biosynthesis, and starch and sugar metabolism were enriched in the downregulated DEGs at D1 and R1, whereas arginine and PRO metabolism was enriched at D1 and R2. Biosynthesis of secondary metabolism, phenylalanine metabolism, and nucleotide excision repair were enriched in the downregulated DEGs at D1 and D2. The pathways related to linoleic acid metabolism; phagosome; and stilbenoid, diarylheptanoid, and gingerol biosynthesis were enriched in the downregulated DEGs at R1 (Fig. 4).

Fig. 4
figure 4

Heatmap diagram reflecting the dynamics of enriched biological processes in KEGG analysis. Drought-rehydration-related physiological categories over the course of the experiment as obtained from KEGG enrichment (see Additional file 3). Enriched processes with FDR < 0.05 were included in the heatmap. Red and blue colors represent upregulated and downregulated DEGs, respectively. The intensity of the color reflects the number of DEGs as indicated at each physiological stage

Response of key unigenes to photosynthesis, antioxidant, osmoregulation and plant hormone signal transduction processes

To further understand the response of photosynthesis, antioxidant, osmoregulation and plant hormone signal transduction processes to drought and recovery processes, annotated DEGs with fpkm > 2 and verified functions related to six pathways (including carbon fixation in photosynthetic organisms, photosynthesis, photosynthesis - antenna proteins, peroxisome, arginine and PRO metabolism and plant hormone signal transduction) were selected and grouped into clusters (Table 3, 4 and 5). In Giant Juncao under D1, 4, 7, and 2 genes related to the carbon fixation in photosynthetic organisms, photosynthesis, and photosynthesis - antenna proteins were upregulated, whereas 23, 4, and 2 genes were downregulated, respectively (Additional file 4). These genes mainly included encoding alanine aminotransferase 2 (ALT2) (c116433_g3_i1), alanine aminotransferase (ALT) (c111450_g1_i1), phosphoenolpyruvate carboxylase (PPC) (c114318_g2_i1, c119399_g3_i1, c98555_g1_i1), and photosystem II protein M (PsbM) (c47941_g1_i1), which were downregulated, and photosystem I P700 apoprotein A1 (PsaA) (c8283_g1_i1), which was upregulated. With the increasing degree of drought stress, the genes encoding aspartate aminotransferase (AST) (c116017_g2_i1), fructose-bisphosphate aldolase cytoplasmic isozyme (FBA) (c112618_g1_i3), malate dehydrogenase (MDH) (c101919_g2_i1, c112985_g2_i3, c112985_g2_i2), triosephosphate isomerase (TPI1) (c111215_g2_i3), ribose-5-phosphate isomerase (RPI) (c118707_g1_i7), and aldolase C-2 (ALDOC-2) (c75858_g1_i1) were upregulated, and the |log2(FoldChange)| was between 1.5 and 4.1 at D2. However, the |log2(FoldChange)| of the downregulated genes was between 1.1 and 12.5, especially the genes (c111709_g1_i1, c112618_g1_i1, c117517_g1_i1, c116040_g1_i1) whose |log2(FoldChange)| is more than 9; these genes separately encoded chloroplast fructose-1,6-biphosphate aldolase (CpFBA), fructose-bisphosphate aldolase (FDA), P-pyruvate carboxylase (PEPC), and chloroplast ribulose-1,5-bisphosphate carboxylase/oxygenase small subunit (rbcS-Ma5) (Table 3).

Table 3 Selected DEGs involved in the photosynthesis during drought and rehydration of Giant Juncao.
Table 4 Selected DEGs involved in peroxisome and arginine and proline metabolism of Giant Juncao.
Table 5 Selected DEGs involved in hormone signal transduction during drought stress and rehydration.

After rehydration, all genes were downregulated, and the |log2(FoldChange)| was between 1.7 and 2.9 at R1. With the extension of rehydration time, more genes became upregulated at R2, especially those encoding phosphoribulokinase (PRK) (c108665_g1_i1), oxygen-evolving enhancer protein 1 (OEE1) (c94338_g1_i1), oxygen-evolving enhancer protein 2 (OEE2), photosystem I reaction center subunit II (PsaD) (c108324_g1_i2), rbcS-Ma5 (c116040_g1_i1), oxygen-evolving enhancer protein 3–1 (OEE3–1), photosystem II repair protein (PSB27-H1), small subunit of ribulose-1,5-bisphosphate carboxylase (rbcS) (c102727_g1_i1), and TPA:malic enzyme (c112733_g4_i4) (Table 3). At the end of the rehydration, photosystem II CP43 chlorophyll apoprotein (PsbC) and PPC, which were encoded by c117389_g2_i1 and c119399_g3_i1, were downregulated, whereas other genes were upregulated (Table 3).

Antioxidant and osmoregulation processes play an important role in coping with drought stress. With the aggravation of drought stress, peroxisomal (S)-2-hydroxy-acid oxidase GLO1 (GLO1) (c119055_g1_i1) was downregulated at D2 and upregulated at R3. Cu/Zn SOD (c103167_g1_i1) and SOD [Cu-Zn] 4AP (c103687_g1_i3) were upregulated at D1, D2, and R1 but disappeared at R2 and R3 (Table 4). c104659_g2_i1, c112509_g1_i1, and c105587_g1_i1, which encoded protein Mpv17, PXMP2/4 family protein 4, and peroxisomal membrane protein PMP22, respectively, were downregulated at D1 and D2, whereas xanthine dehydrogenase (XDH) (c119859_g1_i1) was only upregulated at D2. After rehydration, long-chain acyl-CoA synthetase 1 (LACS1) (c113157_g1_i1) was only upregulated at rehydration stages, whereas 3-ketoacyl-CoA thiolase (ACAA2) was downregulated at R2 (Table 4). For the arginine and PRO metabolism pathway, AST, cytoplasmic (Got1) (c116017_g2_i1), delta 1-pyrroline-5-carboxylate synthetase 1 (P5CS1) (c118246_g1_i3), and delta-1-pyrroline-5-carboxylate dehydrogenase (P5CDh) (c118038_g1_i2) were upregulated at D2, whereas arginase (c112173_g1_i4) and amine oxidase flavin domain-containing protein (MAO) (c117003_g1_i2, c117003_g1_i1) were downregulated at R1 and R2, respectively (Table 4).

Plant hormones are secondary metabolites in plants. In response to drought stress, a large number of upregulated and downregulated DEGs are expressed in plant hormone signal transduction, especially at 14 days of drought stress (Additional file 5). With the intensification of drought stress, coronatine-insensitive 1-like protein (c117581_g1_i3), serine/threonine-protein kinase SAPK5 (c78738_g1_i2) and type-A response regulator ARR1 (c112083_g3_i1) were downregulated, while the ABF2 (c112944_g1_i3), ATP binding protein (c115628_g1_i3), histidine-containing phosphotransfer protein (c102835_g2_i2), TPA: SAUR56-auxin-responsive SAUR family member (c120116_g2_i4) were upregulated at D2. After rehydration, ethylene-insensitive protein 2 (c120302_g1_i2) and type-A response regulator ARR1 (c112083_g3_i1) were upregulated at R2, while ABF2, ATP binding protein, and TPA: SAUR56-auxin-responsive SAUR family member exhibited the opposite of what is expressed during drought conditions (Table 5).

Classification of TFs in Giant Juncao under drought stress and rehydration treatment

The TFs of 56 families were included in the DEGs during drought and rehydration (Additional file 6). The largest group of TFs was the basic helix-loop-helix (bHLH) family, followed by the WRKY, whereas other TFs belonged to the NAC, FAR1, B3, MYB-related, and basic leucine zipper (bZIP) families (Fig. 5a). During drought and rehydration treatments, 181, 432, 295, 251, and 114 TFs among the 510, 1090, 521, 642, and 301 TFs in the bHLH family were upregulated at D1, D2, R1, R2, and R3, respectively (Fig. 5a). The 91 upregulated and 281 downregulated TFs in the WRKY family increased to 398 upregulated and 497 downregulated TFs during drought stress. After rehydration, the numbers of upregulated TFs in the WRKY family tended to decrease, whereas those of downregulated TFs tended to increase and then decrease (Additional file 6).

Fig. 5
figure 5

TFs that were differentially expressed under drought and rehydration conditions in Giant Juncao. a: major TF families; b: hierarchical clustering heatmap of 181 bHLH which located at high level of expression (FPKM ≥20); c: the most dominant bHLH with upregulated TFs; d: the most dominant bHLH with downregulated TFs; e: qRT-PCR analysis the relative expression level of c104644_g1_i1; f: qRT-PCR analysis the relative expression level of c109045_g3_i1

At a high expression level (FPKM ≥20), 181 bHLH were used for hierarchical clustering analysis in water deficit and rehydration conditions (Fig. 5b). These bHLH were grouped into five cluster from A to E. The 32 bHLH in cluster A were highly expressed at D2 but have low expression in other stages. In cluster B, 6 and 40 bHLH were in B1 and B2, respectively. The gene c1132_g1_i1located in B1 were highly expressed at R1 and R2. Some genes with low expression were found at D2 in cluster B2. Cluster C consists of 13 members that were highly expressed from CK to D2 but had low expression at R1 to R3. In cluster D, 24 genes were expressed from high to low then back to high expression during drought and rehydration conditions. Cluster E was divided into three subclusters, namely, E1, E2, and E3, which had 5, 31, and 30 genes, respectively. In E1, c107976_g3_i2 and c110523_g1_i1 had the highest expression level at CK, R2 and R3 (Fig. 5b). However, the most dominant bHLH, with 68 upregulated and 118 downregulated TFs, was mrna09125.1-v1.0-hybrid (Figs. 5c and d).

To validate the accuracy of TF sequencing, the relative expression level of two randomly selected bHLH was examined using quantitative RT-PCR analysis (Fig. 5e and f). After water deficit, the relative expression level of c104644_g1_i1 and c109045_g3_i1 decreased significantly over 14 days of drought and then increased during 1 day of recovery. These results demonstrated similar expression patterns to the results of RNA-Seq.

Weighted gene co-expression network analysis module generation and functional enrichment analysis

The gene co-expression profiles of Giant Juncao in response to water deficit and rehydration conditions were analyzed by using WGCNA to detect the relationship between genes and physiological indexes, as well as inter- or intramodular genes. Twenty co-expression modules and correlation coefficients were identified and obtained (Fig. 6). The physiological indexes were positively correlated with the darkorange2, skyblue, and white modules, and the correlation coefficient was between 0.67–0.87. However, Pn and Tr were negatively correlated with turquoise module, MDA was negatively correlated with darkorange2 module, and the correlation coefficient was 0.71–0.8.

Fig. 6
figure 6

Module-trait relationship with physiological indexes. The number represents the correlation coefficient about modules with physiological indexes. The number in the bracket means p-value

The genes related to the physiological indexes of Giant Juncao were mainly concentrated in the turquoise, darkorange2, skyblue, and white modules. Therefore, the genes with a higher weight in each module were selected for analysis and for drawing the network. The genes, including c118251_g1_i1, c118890_g1_i6, c103721_g1_i1, and c120227_g1_i6, which encoded β-glucosidase, peptide/histidine transporter, LOC103649875, and α/β hydrolase, respectively, had the highest degree in four modules (Fig. 7, Additional file 7).

Fig. 7
figure 7

Net-work analysis the hub genes in four modules which had high correlation coefficient (0.67–0.87) with physiological indexes. a: turquoise module; b: darkorange2 module; c: skyblue module; d: white module. Node size means the degree of genes


Drought tolerance is a complex biological process that is controlled and adjusted by several genes in the plant. In this study, a comprehensive method was used to explore the molecular mechanism of physiological response during drought stress and recovery in Giant Juncao, a species without a reference genome. High-throughput sequencing (RNA-seq) and de novo transcriptome analysis [27] were used to investigate the molecular response of Giant Juncao to drought stress.

Insights into the de novo transcriptome assembly and sequence annotation

Due to the lack of genomic and transcriptome data, the genetic background and functional genes of Giant Juncao are not completely clear. Therefore, a large number of unigenes have not been annotated, and only 61.7% of the total unigenes were annotated in at least one database. This is different to the annotation results of transcriptome sequencing of reported species, such as common vetch, which has approximately 500 million clean reads employed to assemble 174,636 transcripts and 83.48% unigenes annotated in at least one database [28]. In the common buckwheat, 53,404 unigenes were assembled, and 39.33% unigenes were annotated in the GO database [29], whereas 8426 and 12,642 unigenes of Taxus species were annotated in the KEGG and GO database, respectively [30]. This indicates that for those species without a reference genome, a large amount of data in the transcriptome still needs to be further mined.

Response of photosynthesis and related processes to drought and rehydration conditions

Through the analysis of transcriptome data, some genes in the photosynthetic pathway showed difference in abundance, which may cause the photosynthetic indices in Giant Juncao to change after no irrigation and rehydration (Fig. 1b). At the beginning of drought stress (D1), the WUE showed a significant upward trend (Fig. 1b). WUE is a comprehensive reflection of photosynthetic and transpiration characteristics of plants. During water deficit, the decrease rate of Pn was slower than that of Tr, so the WUE exhibited an upward trend. However, according to the analysis of transcriptome data, only PsaA showed an upregulated expression (Table 3). PsaA is a core protein and crucial for the functional assembly of PSI [31]. The upregulation of PsaA will keep PSI stable to delay the decline of the photosynthetic rate. Thus, the gene encoding PsaA could adjust higher WUE to cope with light drought stress. Pn decreased significantly with prolongation of stress, which may be caused by the downregulation of PEPC, CpFBA, FDA, and rbcS-Ma5 proteins to cope with drought stress (Fig. 1b and Table 3). PEPC is an important C4 plant photosynthesis multifunctional enzyme, which can be used to transfer into plants to improve photosynthetic efficiency, and becomes an important way to improve crop yield [32,33,34,35]. Under drought stress, PEPC transgenic rice showed strong drought and light tolerance ability because it can accumulate malic acid or oxaloacetic acid in its guard cells to promote the opening of stomata and enhance photosynthesis [32]. In our study, when the Giant Juncao experienced strong water deficit, the PEPC protein was downregulated, causing Pn, Tr, and Gs to decline. These results showed that downregulated PEPC will cause stomatal closure, which decreases Tr and leads to the decline in Pn with less CO2 into the mesophyll cell. CpFBA is one of the key enzymes controlling the rate of photosynthesis by improving carbon fixation efficiency in the Calvin cycle to enhance the resistance of plants [36, 37]. It plays very important roles in improving the adaptability of Sesuvium portulacastrum under salt and drought stress [38] and maintaining the high photosynthetic rate and biomass of transgenic tobacco [39, 40]. In the present study, drought caused the downregulated expression of several key enzyme regulatory genes in photosynthesis, which was also a major factor in the significant decline in Pn of Giant Juncao.

However, at the beginning of re-watering, Pn, Tr, and WUE significantly increased, whereas Gs remained stable. Meanwhile, the expression of genes was downregulated (Fig. 1b and Table 3). This finding reflects that plants cannot quickly recover from drought stress because the stoma does not open immediately. During rehydration for 5 to 9 days, an increasing number of genes that encoded PRK, OEE1, PSI-LHC4, MDH, CAB, OEE2, psaD, rbcS-Ma5, OEE3–1, FD3, PSB27–1, rbcS, and TPA:malic enzymes were upregulated (Table 3). These protein-encoding genes could be the major regulatory genes related to recovery from drought stress. When Giant Juncao was rehydrated again, WUE could reach the control state (Fig. 1b). According to the transcriptome results, the OEE protein family-encoding genes were upregulated. OEEs are the most important proteins for oxygen evolution and photosystem II stability [41]. In our study, the OEE genes played an important role in the process of drought and rehydration of Giant Juncao. The results were similar to those of cultivar wheat and Norway spruce, wherein OEEs were the key enzymes in the photosynthesis system and enhanced in response to drought stress [42, 43].

Antioxidant enzyme and osmotic adjustment response to drought and re-watering conditions of Giant Juncao

Giant Juncao improved its adaptability to cope with drought environment by regulating the ability of antioxidants and the quantity of osmoregulation substances. In this study, SOD, POD, and PRO initially increased under drought stress and then decreased (Fig. 1b). After rehydration, POD and PRO continued to decrease, whereas SOD showed the opposite trend. SOD is a kind of metal enzyme that widely exists in plants. According to the different metal subgroups, SOD is usually divided into three different types, for example, Cu/Zn-SOD, Mn-SOD, and Fe-SOD [44]. In our study, Cu/Zn SOD was the major regulation protein in coping with drought stress. Cu/Zn-SOD is the first-line antioxidant system to remove reactive oxygen and is closely related to plant stress resistance [45]. Peroxisomes are vital cell organelles and may contain highly variable sets of enzymes that adjust many important cellular processes [46]. MPV17 is one of the import proteins in the peroxisomal membrane in eukaryotes [47]. In our study, MPV17 was downregulated under drought stress and disappeared during recovery (Table 4). This finding is probably because this gene acts as a major regulator in the formation of antioxidant enzymes during water shortage conditions. During drought stress, PRO was probably upregulated through P5CS1, which supports the function of PRO as an osmoregulator [48].

Plant hormone reflect drought and rehydration condition of Giant Juncao

Plant hormones not only widely participate in the various growth and development stages of plants but also play an important role in regulating plant growth to adapt to various kinds of biological or abiotic stress. In our research, ABF2, ATP binding protein, TPA: SAUR 56-auxin-responsive SAUR family member, and type-A response regulator ARR1 l showed opposite expression under drought and recovery processes (Table 5). ABF2, as an ABRE-binding factor, is a principal regulator of the ABA-dependent pathway [49]. When Giant Juncao was subjected to drought stress, ABF2 was upregulated at D1 and D2. However, when the dehydration state disappeared, the ABF2 became downregulated at R2, thereby showing that ABF2, as a hormone regulator, can rapidly regulate plant growth to adapt to environmental stress changes (Table 5). A similar result has been reported with Arabidopsis and rice, in which the overexpression of ABF2 significantly increased the drought tolerance of plants [14, 50]. This finding also confirms that ABF2 plays an important role in plant response to drought stress.

TFs regulate the drought resistance of Giant Juncao

Transcription regulation is a key step in plant response to stress by temporarily and spatially regulating their target genes of transcription [51, 52]. TFs plays an important role in plant drought resistance via the transcriptional regulation of downstream genes to improve plant stress resistance [53]. Thus, this study provides important information for discovering and separating TFs and further elucidating the molecular mechanisms underlying drought tolerance in Giant Juncao. More than 15,000 DEGs encoding TFs were distributed in the major TF families, including bHLH, WRKY, NAC, MYB-related, FAR1, B3, and bZIP families (Fig. 5a and Additional file 6). These identified TF families are well known in stress tolerance in plants [28].

The bHLH family is one of the largest TFs in plants and play critical roles in light signaling, hormone signaling, wound, and drought stress response [54]. For example, the bHLH family has 183, 231 and 571 members in the rice, maize and wheat genomes respectively [55], 162 members in Arabidopsis [56], 230 members in Chinese cabbage [57], 146 members in Brachypodium distachyon [58], and 175 members in apple [59] that have been identified, and most of them have been characterized with drought stress responses [60]. In addition, a rice domain gene (OsbHLH148) can provide drought tolerance as a jasmonate signaling module component [61]. However, the overexpression of wheat TabHLH39 genes enhances drought resistance, salt tolerance, and frost resistance in transgenic Arabidopsis [62]. In our study, drought stress induced the upregulation of 843 DEGs and downregulation of 923 DEGs belonging to the bHLH TF gene family (Additional file 6). Moreover, 432 and 658 DEGs were upregulated and downregulated, respectively, after 14 days of drought (Additional file 6). After rehydration, the number of bHLH family members declined (Additional file 6). This finding indicates that the genes regulated by bHLH TF will change significantly under aggravated drought conditions, which proves that bHLH TF plays an active role in the regulation of drought stress.

For the cluster analysis, we divided bHLH into 5 categories and 8 subcategories, and the genes in each subcluster showed similar expression trends (Fig. 5b). Thirty-two bHLH located in cluster A were upregulated under severe drought condition (D2) but were downregulated after rehydration (Fig. 5b). The expression of genes in cluster B was reversed in cluster A, thereby indicating that bHLH located in clusters A and B had opposite functions during the corresponding drought stress and rehydration conditions of Giant Juncao. This result might be due to the large number of bHLH genes and their powerful functions, thereby causing them to have different expressions.

Network analysis revealed the key genes of drought resistance of Giant Juncao

WGCNA analysis revealed the hub genes in each module after the four key modules related to drought resistance of Giant Juncao were identified (Figs. 6, 7 and Additional file 7). In the turquoise module, β-glucosidase was the hub gene related to hormone regulation. Moreover, it was upregulated and downregulated when Giant Juncao experienced severe drought stress and was repaired after rewatering (Additional file 7). β-glucosidase recycling, as a key step in determining ABA concentration, affects drought tolerance and photosynthesis. β-glucosidase isoenzyme contributes significantly to cellular ABA pools and plays an important part in stomatal density and aperture size [63]. In Arabidopsis thaliana leaves, inactive ABA was hydrolyzed by β-glucosidase, and then a large amount of ABA accumulated in the leaves to improve its drought resistance [64]. The result was similar to those of this paper because β-glucosidase was highly expressed under drought stress.


In this study, we report the first transcriptome data for the characterization of drought-rehydration-related genes in Giant Juncao. A total of 93,907 unigenes were de novo assembled, and 57,941 genes were annotated by functional databases. The expression of DEGs and the identification of the hub gene encoding β-glucosidase can help elucidate the response mechanism, rehydration recovery, and physiological indices of plants during drought stress. This study provides not only insights into the genomics of adapting drought tolerance in Giant Juncao but also candidate genetic resources for abiotic stress resistance of plant development.


Plant material and drought-rehydration treatment

Giant Juncao seedlings were collected from the plantation base in Fujian Agriculture and Forestry University. Plants with full sprouts and consistent stem thickness were selected and planted in pots filled with 1:1 turfy soil and vermiculite. Each pot contained 2.5 kg of turfy soil and vermiculite with one stem node of Giant Juncao. All materials were planted in a greenhouse at room temperature (27 °C) and well irrigated. One month later, plants with consistent (seven leaves) and healthy growth were selected for drought experiment.

Drought stress was imposed on the plants by stopping irrigation. Seedlings were sampled by collecting their third and fourth leaves at day 0 for control and at days 7 and 14 (D1 and D2) during drought. After 14 days of drought stress, the seedlings were watered to field capacity, and leaf tissue was collected after 1, 5, and 9 days (R1, R2, and R3) following re-watering. To eliminate the influence of the development process on the test results, a unified end time of the experiment was set, and the sampling was unified. Four independent biological replicates were collected for ecological and physiological tests and three replicates were used for RNA-Seq analysis. The tissues were collected immediately into liquid nitrogen and stored at − 80 °C until used.

Measurement of the ecological and physiological indices of Giant Juncao

In this study, a CIRAS-3 portable photosynthetic apparatus (PP-Systems Company, Amesbury, MA01913, USA) was used to observe the photosynthetic parameters of the seedlings. The CO2 concentration in the reference leaf chamber was controlled at a constant value of 390 μmol·mol− 1, and the photosynthetically active radiation and air relative humidity were set to 1200 μmol m− 2 s− 1 and 75%, respectively. The experimental apparatus automatically recorded the net photosynthetic rate (Pn, μmol m− 2 s− 1), transpiration rate (Tr, mol m− 2 s− 1), and stomatal conductance (Gs, μmol m− 2 s− 1). Meanwhile, a formula was used to calculated the water use efficiency (WUE) as follows: WUE = Pn/Tr, mmol mol− 1. The physiological indices, including superoxide dismutase (SOD), peroxidase (POD) activity, and malondialdehyde (MDA) and proline (PRO) contents, were determined in accordance with a previously reported method [65, 66].

RNA extraction

High quality total RNA was isolated and extracted from the Giant Juncao seedling using the TRIzol reagent (TransGen, Beijing, China) according to the product instructions. RNA purity and concentration were tested using a NanoDrop 1000 spectrophotometer (Thermo Fisher Scientific, Wilmington, DE, USA). The integrity of RNA molecules was measured using an Agilent Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA). Total RNA quality was detected using 1% agarose gels.

Sequencing library preparation for transcriptome

The RNA-seq library was constructed with Illumina’s TruSeq RNA Sample Preparation Kit (Illumina Inc., San Diego, CA, USA) according to the product instruction. mRNA connected with poly-T oligo-attached magnetic beads was purified from total RNA, and divalent cations were used to conduct fragmentation under high temperatures in NEBNext First Strand Synthesis Reaction Buffer (5×). Random hexamer primers and M-MuLV Reverse Transcriptase (RNase H-) were used for first-strand cDNAs from the reverse transcribed fragmented mRNA. Then, the DNA Polymerase I and RNase H (Invitrogen, Carlsbad, CA, USA) were used to synthesize the second-strand cDNAs. The remaining moleculars were converted into blunt ends by exonuclease/polymerase activities. After purification with the AMPure XP system (Beckman Coulter, Pasadena, CA, USA), the cDNA fragments were resolved in elution buffer for end reparation and addition of a poly(A) tail and then connected to sequencing adaptors with suitable length fragments. Libraries were sequenced on an Illumina HiSeq 4000 platform to generate paired-end reads of 150 bp.

Sequence read mapping and assembly

Initial processing of original data in FASTQ file was achieved. Bowtie 2 was used to filter out rRNA [67]. Then, the raw short reads were processed through in-house Perl scripts. In this step, clean reads and more reliable results were obtained by removing reads that contain adapter, low-quality reads with less than 50% bases with quality scores lower than 5 and unknown bases that were larger than 10% N bases. At the same time, GC-content, Q20, Q30 and sequence repetition level of clean data were calculated. Transcriptome de novo assembly was achieved using Trinity software with parameters of Kmer = 25 [68]. The short reads were assembled into contigs on the basis of their overlap regions. Then, CD-HIT ( with 95% global sequence identity was used to obtain unigenes by clustering the sequences of the de novo assembled transcriptome to remove any redundancy.

Functional annotation of gene transcripts

For functional annotation, the assembled unigenes were analyzed by BLAST to by using the Nucleotide database (Nt) with 10–5 E-value threshold [69]. The unigenes were annotated using the BLASTx tool with E-value < 10− 5 against the NCBI non-redundant protein sequences database (Nr), Swiss-Prot protein database, and eukaryotic Cluster of Orthologous Groups of proteins (KOG) database [69]. The search results were imported into Blast2GO version b2g4pipe_v2.5 for gene ontology (GO) assignments [70]. The KOBAS tool (Kegg Automatic Annotation Server) was used for KEGG Orthology and KEGG pathway assignments. In addition, the unigenes were used as query sequences for searching the Pfam (Protein family).

Gene expression quantification and differential expression analysis

Gene expression quantification was performed using the Bowtie aligner and expectation–maximization method (RSEM) to obtain the number of read counts via the Perl script with –est_method RSEM from the Trinity protocol [71, 72]. The expression levels of unigenes were normalized and calculated as the values of fragments per kilobase of transcripts per million mapped fragments (FPKM) during the assembly and clustering process [73]. After standardizing the read count data with TMM, the p-values of edgeR package analyses were adjusted via Benjamini-Hochberg’s method to determine the false discovery rate (FDR) and identify DEGs [74]. Then, the p-value obtained from the test was corrected to get the q-value. The standard of differential gene expression screening is |log2(FoldChange)| > 1 and q-value < 0.05.

GO function and KEGG pathway enrichment analyses

To study the biological significance of DEGs, the GO database was employed to exhibit GO enrichment analysis of DEGs during the different treatments of Giant Juncao using the GOseq (v1.22) software [75]. Parameter settings was that the corrected p- value was less than 0.05. The various metabolic pathways of DEGs were analyzed by using the KEGG database. The statistical enrichment of DEGs was tested using the KOBAS 2.0 web server and the corrected p-value < 0.05 was considered to be significantly enriched in KEGG [76].

TF analysis of DEGs

All identified DEGs were blasted with PlantTFDB (Plant Transcription Factor Database) 4.0 (, and the threshold was set as 1 × e− 5.

Weighted gene co-expression network analysis

WGCNA was used to explore the relationship between genes and physiological indexes, as well as that between genes and genes. The genes with FPKM less than 0.5 were filtered out, and the remaining genes were input into WGCNA network construction (WGCNA v1.69 package in R) [77]. Pearson correlation matrix and network topology analysis were used to calculate the gene correlation and soft thresholding power, respectively. Then, the adjacency was converted to a topological overlap matrix. In standard WGCNA networks, power, minModuleSize, and mergeCutHeight value were set to 7, 30, and 0.25, respectively. The networks were visualized using Cytoscape v3.7.2.

Quantitative RT-PCR validation

To validate the accuracy of the RNA-seq results, quantitative RT-PCR analysis was conducted on a CFX Connect qPCR detection system. Two bHLH transcription factors were randomly selected, and PgACT gene was used as a reference gene (the primers are shown in Table S4 as supplementary data). The 2-∆∆CT method was used to calculate the relative expression levels of genes.

Statistical analysis of physiological data

Statistical analysis of physiological data was conducted using SPSS20.0 (SPSS Inc., Chicago, IL, USA). The significance of differences among every treatment was tested by one-way ANOVA and Duncan’s multiple comparative analysis (P < 0.05).

Availability of data and materials

The transcriptome sequences have been deposited in NCBI under BioProject ID: PRJNA632455 and the URL is . All data generated or analyzed during this study are included in this published article and its supplementary information files (Additional files 1, 2, 3, 4, 5).



Differentially expressed genes


Transcription factors


Basic helix-loop-helix


NAM (no apical meristem), ATAF1/2 and CUC2 (cup-shaped cotyledon);


Fatty acyl-CoA reductase




Basic leucine zipper


Superoxide dismutase








Net photosynthetic rate


Transpiration rate


Stomatal conductance


Water use efficiency


Weighted gene co-expression network analysis


  1. Luo LJ, Xia H, Lu BR. Editorial: crop breeding for drought resistance. Front Plant Sci. 2019;10:314.

    Article  PubMed  PubMed Central  Google Scholar 

  2. Banks JM, Percival GC, Rose G. Variations in seasonal drought tolerance rankings. Trees. 2019;33(4):1063–72.

    Article  Google Scholar 

  3. Pflug EE, Buchmann N, Siegwolf RTW, Schaub M, Rigling A, Arend M. Resilient leaf physiological response of european beech (Fagus sylvatica L.) to summer drought and drought release. Front. Plant Sci. 2018;9:187.

    Google Scholar 

  4. Hein JA, Sherrard ME, Manfredi KP, Abebe T. The fifth leaf and spike organs of barley (Hordeum vulgare L.) display different physiological and metabolic responses to drought stress. BMC Plant Biol. 2016;16(1):248.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  5. Torres-Franklin ML, Gigon A, Melo DFD, Zuily-Fodil Y, Pham-Thi AT. Drought stress and rehydration affect the balance between MGDG and DGDG synthesis in cowpea leaves. Physiol Plant. 2007;131(2):201–10.

    CAS  PubMed  Google Scholar 

  6. Furlan A, Bianucci E, Tordable MDC, Kleinert A, Valentine AJ, Castro S. Dynamic responses of photosynthesis and the antioxidant system during a drought and rehydration cycle in peanut plants. Funct Plant Biol. 2016;43(4):337–45.

    Article  CAS  PubMed  Google Scholar 

  7. Tina RR, Shan XR, Wang Y, Guo SY, Mao B, Wang W, et al. Response of antioxidant system to drought stress and rewatering in alfalfa during branching. IOP Conf Ser Earth Environ Sci. 2017;94:14.

    Article  Google Scholar 

  8. Xu ZZ, Zhou GS, Shimizu H. Are plant growth and photosynthesis limited by pre-drought following rewatering in grass? J Exp Bot. 2009;60(13):3737–49.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  9. Farooq M, Wahid A, Kobayashi N, Fujita D, Basra SMA. Plant drought stress: effects, mechanisms and management. Agron Sustain Dev. 2009;29(1):185–212.

    Article  Google Scholar 

  10. Guo YY, Yu HY, Yang MM, Kong DS, Zhang YJ. Effect of drought stress on lipid peroxidation, osmotic adjustment and antioxidant enzyme activity of leaves and roots of Lycium ruthenicum Murr. Seedling. Russ J Plant Physiol. 2018;65(2):244–50.

    Article  CAS  Google Scholar 

  11. Gong ZZ, Xiong LM, Shi HZ, Yang SH, Herrera-Estrella LR, Xu GH, et al. Plant abiotic stress response and nutrient use efficiency. Sci China Life Sci. 2020;63(5):635–74.

    Article  PubMed  Google Scholar 

  12. Maruyama K, Urano K, Yoshiwara K, Morishita Y, Sakurai N, Suzuki H, et al. Integrated analysis of the effects of cold and dehydration on rice metabolites, phytohormones, and gene transcripts. Plant Physiol. 2014;164(4):1759–71.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  13. Fujita Y, Fujita M, Satoh R, Maruyama K, Parvez MM, Seki M, et al. AREB1 is a transcription activator of novel ABRE-dependent ABA signaling that enhances drought stress tolerance in Arabidopsis. Plant Cell. 2005;17(12):3470–88.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  14. Hossain MA, Cho JI, Han M, Ahn CH, Jeon JS, An G, et al. The ABRE-binding bZIP transcription factor OsABF2 is a positive regulator of abiotic stress and ABA signaling in rice. J Plant Physiol. 2010;167(17):1512–20.

    Article  CAS  PubMed  Google Scholar 

  15. Kim S, Kang JY, Cho DI, Park JH, Kim SY. ABF2, an ABRE-binding bZIP factor, is an essential component of glucose signaling and its overexpression affects multiple stress tolerance. Plant J. 2004;40(1):75–87.

    Article  CAS  PubMed  Google Scholar 

  16. Agarwal PK, Gupta K, Lopato S, Agarwal P. Dehydration responsive element binding transcription factors and their applications for the engineering of stress tolerance. J Exp Bot. 2017;68(9):2135–48.

    Article  CAS  PubMed  Google Scholar 

  17. Fox H, Doron-Faigenboim A, Kelly G, Zhou J, Moshe Y, David-Schwartz R, et al. Transcriptome analysis of Pinus halepensis under drought stress and during recovery. Tree Physiol. 2018;38(3):423–41.

    Article  CAS  PubMed  Google Scholar 

  18. Liu F, Lin D, Lin H, Liu H, Lin Z. Physiological and photosynthetic responses of Giant JunCao (Pennisetum Giganteum) to drought stress. Fresenius Environ Bull. 2017;26(6):3868–71.

    CAS  Google Scholar 

  19. Martin LBB, Fei ZJ, Giovannoni JJ, Rose JKC. Catalyzing plant science research with RNA-seq. Front Plant Sci. 2013;4:1.

    Article  Google Scholar 

  20. Zhang H, He L, Cai L. Transcriptome sequencing: RNA-Seq. Methods Mol Biol. 1754;2018:15–27.

    Google Scholar 

  21. Li W, Jiang T. Transcriptome assembly and isoform expression level estimation from biased RNA-Seq reads. Bioinform. 2012;28(22):2914–21.

    Article  CAS  Google Scholar 

  22. Camilios-Neto D, Bonato P, Wassem R, Tadra-Sfeir MZ, Brusamarello-Santos LCC, Valdameri G, et al. Dual RNA-seq transcriptional analysis of wheat roots colonized by Azospirillum brasilense reveals up-regulation of nutrient acquisition and cell cycle genes. BMC Genet. 2014;15:378.

    Article  CAS  Google Scholar 

  23. Zhang LM, Liu XG, Qu XN, Yu Y, Han SP, Dou Y, et al. Early transcriptomic adaptation to Na2CO3 stress altered the expression of a quarter of the total genes in the maize genome and exhibited shared and distinctive profiles with NaCl and high pH stresses. J Integr Plant Biol. 2013;55(11):1147–65.

    Article  CAS  PubMed  Google Scholar 

  24. Zhou Y, Yang P, Cui F, Zhang FT, Luo XD, Xie JK. Transcriptome analysis of salt stress responsiveness in the seedlings of Dongxiang wild rice (Oryza rufipogon Griff.). Plos One. 2016;11(1):e0146242.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  25. Fracasso A, Trindade LM, Amaducci S. Drought stress tolerance strategies revealed by RNA-Seq in two sorghum genotypes with contrasting WUE. BMC Plant Biol. 2016;16(1):115.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  26. Wang YQ, Li L, Tang S, Liu JG, Zhang HS, Zhi H, et al. Combined small RNA and degradome sequencing to identify miRNAs and their targets in response to drought in foxtail millet. BMC Genet. 2016;17:57.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  27. Chen YJ, Chen Y, Shi ZJ, Jin YF, Sun HS, Xie FC, et al. Biosynthesis and signal transduction of ABA, JA, and BRs in response to drought stress of Kentucky bluegrass. Int J Mol Sci. 2019;20(6):1289.

    Article  CAS  PubMed Central  Google Scholar 

  28. Zhu YQ, Liu QX, Xu WZ, Zhang JH, Wang X, Nie G, et al. De Novo assembly and discovery of genes that involved in drought tolerance in the common vetch. Int J Mol Sci. 2019;20(2):328.

    Article  PubMed Central  CAS  Google Scholar 

  29. Wu Q, Zhao G, Bai X, Zhao W, Xiang DB, Wan Y, et al. Characterization of the transcriptional profiles in common buckwheat (Fagopyrum esculentum) under PEG-mediated drought stress. Electron J Biotechnol. 2019;39:42–51.

    Article  CAS  Google Scholar 

  30. Zhou T, Luo XJ, Yu CN, Zhang CC, Zhang L, Song YB, et al. Transcriptome analyses provide insights into the expression pattern and sequence similarity of several taxol biosynthesis-related genes in three Taxus species. BMC Plant Biol. 2019;19:33.

    Article  PubMed  PubMed Central  Google Scholar 

  31. Leelavathi S, Bhardwaj A, Kumar S, Dass A, Pathak R, Pandey SS, et al. Genome-wide transcriptome and proteome analyses of tobacco psaA and psbA deletion mutants. Plant Mol Biol. 2011;76:407–23.

    Article  CAS  PubMed  Google Scholar 

  32. Ku MSB, Agarie S, Nomura M, Fukayama H, Matsuoka M, Ono K, et al. High-level expression of maize phosphoenolpyruvate carboxylase in transgenic rice plants. Nat Biotechnol. 1999;17(1):76–80.

    Article  CAS  PubMed  Google Scholar 

  33. Lebouteiller B, Rolin D. Pierre J-Nl, Bleton J, Tchapla a, Maucourt M, et al. physiological impacts of modulating phosphoenolpyruvate carboxylase levels in leaves and seeds of Arabidopsis thaliana. Plant Sci. 2007;172(2):265–72.

    Article  CAS  Google Scholar 

  34. Ort DR, Merchant SS, Alric J, Barkan A, Blankenship RE, Bock R, et al. Redesigning photosynthesis to sustainably meet global food and bioenergy demand. Proc Natl Acad Sci Unit States Am. 2015;112(28):8529–36.

    Article  CAS  Google Scholar 

  35. Sánchez R, Flores A, Cejudo FJ. Arabidopsis phosphoenolpyruvate carboxylase genes encode immunologically unrelated polypeptides and are differentially expressed in response to drought and salt stress. Planta. 2006;223(5):901–9.

    Article  PubMed  CAS  Google Scholar 

  36. Bowler C, Allen AE, Moustafa A, Montsant A, Eckert A, Kroth PG. Evolution and functional diversification of fructose bisphosphate aldolase genes in photosynthetic marine diatoms. Mol Biol Evol. 2011;29(1):367–79.

    PubMed  PubMed Central  Google Scholar 

  37. Ma WM, Wei LZ, Wang QX, Shi DJ, Chen HB. Increased activity of the tandem fructose-1,6-bisphosphate aldolase, triosephosphate isomerase and fructose-1,6-bisphosphatase enzymes in Anabaena sp. strain PCC 7120 stimulates photosynthetic yield. J Appl Phycol. 2008;20(4):437–43.

    Article  CAS  Google Scholar 

  38. Fan W, Zhang ZL, Zhang YL. Cloning and molecular characterization of fructose-1,6-bisphosphate aldolase gene regulated by high-salinity and drought in Sesuvium portulacastrum. Plant Cell Rep. 2009;28(6):975–84.

    Article  CAS  PubMed  Google Scholar 

  39. Simkin AJ, McAusland L, Headland LR, Lawson T, Raines CA. Multigene manipulation of photosynthetic carbon assimilation increases CO2 fixation and biomass yield in tobacco. J Exp Bot. 2015;66:13.

    Article  CAS  Google Scholar 

  40. Uematsu K, Suzuki N, Iwamae T, Inui M, Yukawa H. Increased fructose 1,6-bisphosphate aldolase in plastids enhances growth and photosynthesis of tobacco plants. J Exp Bot. 2012;63(8):3001–9.

    Article  CAS  PubMed  Google Scholar 

  41. Li T, Liu LN, Jiang CD, Liu YJ, Shi L. Effects of mutual shading on the regulation of photosynthesis in field-grown sorghum. J Photochem Photobiol B. 2014;137:31–8.

    Article  PubMed  CAS  Google Scholar 

  42. Blodner C, Majcherczyk A, Kues U, Polle A. Early drought-induced changes to the needle proteome of Norway spruce. Tree Physiol. 2007;27(10):1423–31.

    Article  PubMed  Google Scholar 

  43. Cheng ZW, Dong K, Gei P, Bian YW, Dong LW, Deng X, et al. Identification of leaf proteins differentially accumulated between wheat cultivars distinct in their levels of drought tolerance. PLoS One. 2015;10(5):e0125302.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  44. Ogawa K, Kanematsu S, Asada K. Intra- and extra-cellular localization of “cytosolic” CuZn-superoxide dismutase in spinach leaf and hypocotyl. Plant Cell Physiol. 1996;37:6.

    Article  Google Scholar 

  45. Marikovsky M, Ziv V, Nevo N, Harris-Cerruti C, Mahler O. Cu/Zn superoxide dismutase plays important role in immune response. J Immunol. 2003;170(6):2993–3001.

    Article  CAS  PubMed  Google Scholar 

  46. Kiel JAKW, Emmrich K, Meyer HE, Kunau WH. Ubiquitination of the peroxisomal targeting signal type 1 receptor, pex5p, suggests the presence of a quality control mechanism during peroxisomal matrix protein import. J Biol Chem. 2005;280(3):1921–30.

    Article  CAS  PubMed  Google Scholar 

  47. Charton L, Plett A, Linka N. Plant peroxisomal solute transporter proteins. J Integr Plant Biol. 2019;61(7):817–35.

    CAS  PubMed  PubMed Central  Google Scholar 

  48. Han YS, Yang H, Wu MY, Yi HL. Enhanced drought tolerance of foxtail millet seedlings by sulfur dioxide fumigation. Ecotox Environ Safe. 2019;178:9–16.

    Article  CAS  Google Scholar 

  49. Wu JJ, Yan GB, Duan ZQ, Wang ZJ, Kang CY, Guo L, et al. Roles of the Brassica napus DELLA protein BnaA6.RGA, in modulating drought tolerance by interacting with the ABA signaling component BnaA10.ABF2. Front. Plant Sci. 2020;11:577.

    Google Scholar 

  50. Li XY, Liu X, Yao Y, Li YH, Liu S, He CY, et al. Overexpression of Arachis hypogaea AREB1 gene enhances drought tolerance by modulating ROS scavenging and maintaining endogenous ABA content. Int J Mol Sci. 2013;14:12827–42.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  51. Wei C, Yao QM, Li S, Agarwal G, Wang B, Li L, et al. Identification and comparative analysis of differential gene expression in soybean leaf tissue under drought and flooding stress revealed by RNA-Seq. Front Plant Sci. 2016;7:1044.

    Google Scholar 

  52. Jin JP, Zhang H, Kong L, Gao G, Luo JC. PlantTFDB 3.0: a portal for the functional and evolutionary study of plant transcription factors. Nucleic Acids Res. 2014;42(1):1182–7.

    Article  CAS  Google Scholar 

  53. Gonçalves LP, Boscariol Camargo RL, Takita MA, Machado MA, dos Soares Filho WS, Costa MGC. Rootstock-induced molecular responses associated with drought tolerance in sweet orange as revealed by RNA-Seq. BMC Genomics. 2019;20:110.

    Article  PubMed  PubMed Central  Google Scholar 

  54. Lorenzo CP, Anahit G, Irma RV, Martinez-Garcia JF, Bilbao-Castro JR, Robertson DL. Genome-wide classification and evolutionary analysis of the bHLH family of transcription factors in Arabidopsis, poplar, rice, moss, and algae. Plant Physiol. 2010;153:1398–412.

    Article  CAS  Google Scholar 

  55. Wei KF, Chen HQ. Comparative functional genomics analysis of bHLH gene family in rice, maize and wheat. BMC Plant Biol. 2018;18:309.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  56. Bailey PC, Martin C, Toledo-Ortiz G, Quail P, Huq E, Heim M, et al. Update on the basic helix-loop-helix transcription factor gene family in Arabidopsis thaliana. Plant Cell. 2003;15(11):2497–501.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  57. Song XM, Huang ZN, Duan WK, Ren J, Liu TK, Li Y, et al. Genome-wide analysis of the bHLH transcription factor family in Chinese cabbage (Brassica rapa ssp.pekinensis) mol. Genet Genomics. 2014;289(1):77–91.

    Article  CAS  Google Scholar 

  58. Niu X, Guan YX, Chen SK, Li HF. Genome-wide analysis of basic helix-loop-helix (bHLH) transcription factors in Brachypodium distachyon. BMC Genomics. 2017;18(1):619.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  59. Mao K, Dong QL, Li C, Liu CH, Ma FW. Genome wide identification and characterization of apple bHLH transcription factors and expression analysis in response to drought and salt stress. Front Plant Sci. 2017;8:1.

    Article  Google Scholar 

  60. Cui X, Wang YX, Liu ZW, Wang WL, Li H, Zhuang J. Transcriptome-wide identification and expression profile analysis of the bHLH family genes in Camellia sinensis. Funct Integr Genomics. 2018;18:10.

    Article  CAS  Google Scholar 

  61. Seo JS, Joo J, Kim MJ, Kim YK, Nahm BH, Song SI, et al. OsbHLH148, a basic helix-loop-helix protein, interacts with OsJAZ proteins in a jasmonate signaling pathway leading to drought tolerance in rice. Plant J. 2011;65(6):907–21.

    Article  CAS  PubMed  Google Scholar 

  62. Zhai YQ, Zhang LC, Xia C, Fu SL, Zhao GY, Jia JZ, et al. The wheat transcription factor, TabHLH39, improves tolerance to multiple abiotic stressors in transgenic plants. Biochem Biophys Res Commun. 2016;473:4.

    Article  Google Scholar 

  63. Wang CL, Chen S, Dong YP, Ren RJ, Chen DF, Chen XW. Chloroplastic Os3BGlu6 contributes significantly to cellular ABA pools and impacts drought tolerance and photosynthesis in rice. New Phytol. 2020;226:1042–54.

    Article  CAS  PubMed  Google Scholar 

  64. Wang Z, Wang FX, Hong YC, Yao JJ, Ren ZZ, Shi HZ, et al. The flowering repressor SVP confers drought resistance in Arabidopsis by regulating abscisic acid catabolism. Mol Plant. 2018;11:1184–97.

    Article  CAS  Google Scholar 

  65. Jiang AL, Tian SP, Xu YZ. Effects of controlled atmospheres with high_O_2 or high_CO_2 concentrations on postharvest physiology and storability of "Napoleon" sweet cherry. J Integr Plant Biol. 2002;44(8):925–30.

    CAS  Google Scholar 

  66. Yingsanga P, Srilaong V, Kanlayanarat S, Noichinda S, McGlasson WB. Relationship between browning and related enzymes (PAL, PPO and POD) in rambutan fruit (Nephelium lappaceum Linn.) cvs. Rongrien and see-Chompoo. Postharvest Biol Technol. 2008;50(2):164–8.

    Article  CAS  Google Scholar 

  67. Langmead B, Salzberg SL. Fast gapped-read alignment with bowtie 2. Nat Methods. 2012;9(4):357–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  68. Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, et al. Full-length transcriptome assembly from RNA-seq data without a reference genome. Nat Biotech. 2011;29(7):644–52.

    Article  CAS  Google Scholar 

  69. Altschul SF, Madden TL, Schäffer AA, Zhang J, Zhang Z, Miller W, et al. Gapped BLAST and PSI-BLAST: a new generation of protein database search programs. Nucleic Acids Res. 1997;25(17):3389–402.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  70. Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ. Basic local alignment search tool. J Mol Biol. 1990;215:403–10.

    Article  CAS  PubMed  Google Scholar 

  71. Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinform. 2011;12(1):323.

    Article  CAS  Google Scholar 

  72. Haas BJ, Papanicolaou A, Yassour M, Grabherr M, Blood PD, Bowden J, et al. De novo transcript sequence reconstruction from RNA-Seq: reference generation and analysis with trinity. Nat Protoc. 2013;8(8):1–43.

    Article  CAS  Google Scholar 

  73. Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, Marijke J. Van Baren, et al. transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat Biotech. 2010;28(5):511–5.

    Article  CAS  Google Scholar 

  74. Smyth GK, Robinson MD, McCarthy DJ. edgeR: a bioconductor package for differential expression analysis of digital gene expression data. Bioinform. 2010;26(1):139–40.

    Article  CAS  Google Scholar 

  75. Young MD, Wakefield MJ, Smyth GK, Oshlack A. Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biol. 2010;11:R14.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

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

    Article  CAS  PubMed  Google Scholar 

  77. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008;9:559.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

Download references


Not applicable.


This work was supported by the Natural Science Foundation of Fujian Province (2019J05053), Education department of Fujian Province (JT180120) and National Natural Science Foundation of China (41775105). The funders had no role in study design, collection, analysis, and interpretation of data and in writing the manuscript, but just provide the financial support.

Author information

Authors and Affiliations



JZ conceived, designed and performed the experiments and wrote the manuscript. SQC and WJS performed parts of the experiments. RDS revised and edited the manuscript. STL and FLY made substantial contributions to the data analysis. ZXL conceived the experiments and provided the Giant Juncao seedling. All authors have read and approved the final manuscript.

Corresponding author

Correspondence to Jing Zhou.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

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

Supplementary Information

Additional file 1: Table S1-S4

. Table S1: Overview of the sequencing. Table S2: Unigene information annotated in different databases. Table S3: Number of differential genes annotated by KEGG pathway. Table S4: The primers of two bHLH TF and PgACT.

Additional file 2.

GO enrichment in DEGs when FDR<0.01.

Additional file 3.

KEGG enrichment in DEGs at different drought stress and rehydration conditions.

Additional file 4.

DEGs associated with photosynthesis at different drought stress and rehydration conditions.

Additional file 5.

DEGs associated with plant hormone signal transduction pathway at different drought stress and rehydration condition.

Additional file 6.

TF analysis results.

Additional file 7.

Degree and log2(FoldChange) value of hub genes under different treatments in the four modules.

Rights and permissions

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

Reprints and Permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Zhou, J., Chen, S., Shi, W. et al. Transcriptome profiling reveals the effects of drought tolerance in Giant Juncao. BMC Plant Biol 21, 2 (2021).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: