Abstract Background Gymnocarpos przewalskii Bunge ex Maxim. (G. przewalskii) is an endangered xerophytic shrub that plays a crucial role as a source of forage in the Alxa Desert. However, there is limited understanding regarding the forage quality of G. przewalskii and its response to drought. This study aimed to evaluate the forage quality of G. przewalskii and investigate the physiological and transcriptomic changes in G. przewalskii response to drought stress. Results The ash, fat, crude protein, lignin, crude fiber, acid detergent fiber, and neutral detergent fiber contents in G. przewalskii twigs were 10.61%, 1.85%, 5.68%, 7.08%, 21.23%, 42.16%, and 58.42%, respectively. In contrast, these ingredients in its leaves were 20.39%, 0.92%, 11.96%, 2.40%, 17.51%, 14.29% and 20.26%, respectively. Osmotic stress led to a reduction in chlorophyll levels and an increase in malondialdehyde content. Levels of hydrogen peroxide and oxygen free radicals remained relatively stable under osmotic stress. The proline content, SOD and CAT activities, and ·OH scavenging capacity were enhanced in G. przewalskii under osmotic stress. RNA-sequencing of G. przewalskii generated 44.51 Gb clean reads, which were assembled into 102,191 Unigenes and 30,809 Unigenes were successfully annotated. Comparative analysis identified 3,015 differentially expressed genes under osmotic stress. There were 2,134 and 1,739 DEGs enriched in 47 GO secondary categories and 129 KEGG pathways, respectively. 2 up-regulated DEGs were annotated to P5CS, a key enzyme in the biosynthesis of proline. 32 DEGs were annotated to various antioxidases and antioxidants. 81 DEGs were annotated to 8 plant hormone signaling pathways, in which the auxin and ABA signaling pathways exhibited dominant enrichment. 150 DEGs were annotated to 35 transcription factor families with the abundant enrichment of TF families containing WRKY, bZIP, ERF, bHLH, MYB, and NAC. Conclusions High forage quality and drought stress tolerance were observed in G. przewalskii. In response to drought stress, G. przewalskii orchestrates reactive oxygen species scavenging, proline biosynthesis, and other intricate physiological processes, with substantial contributions from plant hormones and transcription factors. This study provides new insights into the forage quality and the mechanisms involved in drought adaptation of G. przewalskii, offering a foundation for its conservation and sustainable utilization. Supplementary Information The online version contains supplementary material available at 10.1186/s12870-025-06185-7. Keywords: Gymnocarpos przewalskii Bunge ex Maxim., Forage quality, Drought stress, RNA-sequencing, Antioxidant system, Plant hormone Introduction Plants encounter multiple abiotic stresses under field conditions due to their immobile nature, and drought is one of the most severe abiotic stresses at a global level [[32]1, [33]2]. According to the Intergovernmental Panel on Climate Change report, temperatures are projected to rise by a minimum of 1.5 degrees by 2040 because of the increasing CO[2] levels in the atmosphere [[34]2]. Due to global temperature escalation, the frequency and intensity of drought have continuously increased, which is threatening environmental biodiversity and crop yields [[35]3]. The report of the United Nations Environment Program showed that the crop yields in semi-arid and arid regions were reduced by approximately 3.6 billion hectares (25% of the world’s highlands) [[36]4]. Therefore, exploring the mechanism of plant adaption to drought stress and improving the utilization of semi-arid and arid regions is closely related to human development. As a common abiotic stress in the natural environment, drought stress causes water deficit in plants [[37]5]. Water deficit will lead to a loss of turgidity and pressure against the cell wall, both essential for maintaining plant structure and growth [[38]6]. Meanwhile, water is the main solvent and material for a lot of physiological and biochemical processes in plant cells, including photosynthesis, protein folding and assembly, signal perception, and material transportation. For plants exposed to drought stress, the water deficit will interfere with these processes in plant cells, resulting in inhibition of enzyme activities, loss of chlorophyll, photosynthesis rates drop, metabolic imbalance, and accumulation of harmful compounds [[39]7–[40]9]. Furthermore, drought stress can cause the excessive accumulation of reactive oxygen species (ROS). Excessive ROS during drought stress can oxidize unsaturated fatty acids in the cell membrane, damage the membrane structure, and cause secondary oxidative stress to plants [[41]6]. To survive in drought stress, plants should maintain enough water and eliminate excessive ROS in plant cells through different physiological pathways, such as reducing stomatal conductance to inhibit transpiration, accumulating more osmolytes to decrease the water potential, and elevating antioxidase activities to maintain ROS balance [[42]6, [43]10]. For example, TaWRKY24 transgenic tobacco promoted the accumulation of proline and soluble sugars to increase the relative water content and elevated the activities of SOD, POD, and CAT to decrease the oxidative membrane damage in cells under drought stress [[44]11]. In rice, ZmGLK transgenic plants had higher stomatal closure than wild-type plants, resulting in ZmGLK transgenic plants having higher drought tolerance [[45]12]. Additionally, when the abiotic stress-related genes ONAC023, PGL3 and OsSF3B1 were knocked out in rice, the drought tolerance of onac023, pgl3 and ossf3b1 mutants were all lower than that of wild-type rice, with the higher H[2]O[2] accumulation were observed in three mutants [[46]10]. Alxa desert (37°25′17.03″~41°50′49.89″ N, 99°21′28.97″~107°4′29.60″ E) is a typical desert, which consists of Badain Jaran, Tengger and Ulanbuh deserts. It locates in the northwest of China and has an average elevation of 1250 m. The total area of Alxa Desert is more than 100,000 km^2, accounting for about one third of the total area of Alxa Plateau [[47]13]. This area has a variety of climate extremes, including hyper-arid conditions (annual average precipitation of 32 ~ 180 mm, annual average evaporation of 2348 ~ 4203 mm), extreme temperature (-34.4℃~44.5℃), intensive UV radiation, high soil salinity, poor soil and sandstorm [[48]14]. However, in this extreme arid desert, there are many rare, relic and endangered plants, such as Zygophyllum xanthoxylon, Nitraria tangutorum, Reaumuria songarica, Ammopiptanthus mongolicus, Haloxylon ammodendron, Amygdalus mongolica, and Gymnocarpos przewalskii Bunge ex Maxim. (G. przewalskii) [[49]14–[50]17]. These plants exhibit high drought resistance and serve as excellent materials for studying the mechanisms of plant response to drought stress. Meanwhile, they provide the necessary herbage for local camels, cattle, sheep and other livestock, which play an important role in improving the economic income of local herdsmen. G. przewalskii (Gymnocarpos genus, Caryophyllaceae family), an endangered xerophytic shrub, is a keystone species in the arid and semi-arid regions of Alxa desert [[51]18, [52]19]. This species has developed distinct morphological characteristics that is characterized by succulent leaves, sunken stomata, thick and waxy cuticle, well-developed palisade tissue, and a high ratio of root to stem [[53]19]. These characteristics make G. przewalskii could tolerate extreme drought and suitable for surviving in the arid and semi-arid deserts. In addition, the twigs and leaves of G. przewalskii contain rich nutrients and have good palatability, making it a highly popular herbage for livestock, such as camels, cattle, and sheep [[54]20]. Because of the rich nutrients and high drought resistance, G. przewalskii has played an important role in controlling desert and promoting economic development in the Alxa desert [[55]19, [56]20]. Previous studies reported that G. przewalskii could adapt to drought stress through various physiological processes. The seeds of G. przewalskii would go into dormancy when the drought stress exceeded its tolerance range, however, about 80% of the dormant seeds began to germinate rapidly after the drought stress was relieved, indicating that G. przewalskii could avoid severe drought stress by seed dormancy [[57]21, [58]22]. In addition, the soluble protein content, soluble sugar content, proline content and antioxidase activity all increased in G. przewalskii under drought stress, indicating that the osmolytes and antioxidases involed in G. przewalskii response to drought stress [[59]23]. However, the physiological processes and the change of genes expressions in G. przewalskii response to drought stress and its forage quality remains unclear. In present study, the forage quality of twigs and leaves, which were collected from G. przewalskii in Alxa desert, were tested. To mimic drought stress, the seedlings of G. przewalskii were treated with 15% polyethylene glycol (PEG) 6000 to induce osmotic stress, and the changes of physiological processes and transcriptome were tested. We had assessed the forage quality of G. przewalskii, and integrated physiological and transcriptome analyses to investigate specialized mechanism of G. przewalskii responses to drought stress. Methods Plant materials and stress treatments The twigs, leaves and seeds were collected from 3 healthy Gymnocarpos przewalskii Bunge ex Maxim. (G. przewalskii) plants with same size (60–70 cm height, 90–100 cm crown diameter) in the natural environment of Yabalai Mountain (39°20′85″N, 101°66′53″E) in Alxa desert, China. The length of twigs, leaves and seeds were about 5 cm, 1 cm and 1 mm, respectively. The ambient natural environmental conditions during collection were 29℃, 85,000 lx, and 17% relative humidity. The intact plump seeds of G. przewalskii with same size were sterilized by immersing in 10% sodium hypochlorite solution for 15 min and then washed 3 times with distilled water to remove the bleach. Subsequently, the seeds were cultivated in a 150 mL glass culture flask containing 40 mL Murashige and Skoog (MS) solid medium at 25 °C under 16-hours (h) light/8-h dark photoperiod for 2 weeks. When seedlings were approximately 5 cm high, they were transferred into a glass tube containing 50 mL MS liquid medium and further cultured for another 4 weeks at same condition. Three healthy seedlings with the same size were selected for further osmotic stress treatments, which were cultivated in MS liquid medium containing 15% PEG 6000 (w/v) for 1 day (d), 3 d, 5 d, and 7 d. The sodium hypochlorite solution and PEG 6000 were bought from Sangon Biotech Co., Ltd (Shanghai, China). The MS medium was bought from PhytoTech Labs, Inc. (Kansas, America). Forage quality measurements The twigs and leaves of G. przewalskii were dried to constant weight at 65 °C. Then its ash content, fat content, crude protein content, lignin content, crude fiber content, acid detergent fiber (ADF) content, and neutral detergent fiber (NDF) were measured based on the dry weight (DW) according to the Chinese standards GB/T 6438 − 2007, GB/T 6433 − 2006, GB/T 24,318 − 2009, GB/T 20,805 − 2006, GB/T 6434 − 2006, NY/T 1459–2007, and GB/T 20,806 − 2006 ([60]https://www.ndls.org.cn/), respectively. Three biological replicates were performed for each sample. All the chemical reagents used for forage quality measurements were bought from Sangon Biotech Co., Ltd (Shanghai, China). Physiological parameter measurements The leaves of G. przewalskii seedlings under normal condition and different osmotic stresses were collected, and its physiological parameters were measured. The assay kits of chlorophyll content, malondialdehyde (MDA) content, hydrogen peroxide (H[2]O[2]) content, oxygen free radical (OFR) content, proline (Pro) content, superoxide dismutase (SOD) activity, catalase (CAT) activity, and hydroxyl radical (·OH) scavenging capacity were bought from Suzhou Comin Biotechnology Co., Ltd (Suzhou, China), and were used to corresponding physiological parameter measurements. Three biological replicates were performed for each sample. RNA‑Sequencing analysis The 6-week-old healthy seedlings of G. przewalskii with the same size were selected for RNA‑Sequencing. The leaves of G. przewalskii seedlings treated by osmotic stress for 3 d were collected and preserved at -80 °C for subsequent RNA‑Sequencing experiments. As a control group, the leaves of G. przewalskii seedlings under normal condition, which were cultivated in MS liquid medium without PEG 6000 for 3 d, were also collected and preserved at -80 °C for subsequent RNA‑Sequencing experiments. Total RNA extraction, RNA-Sequencing library construction, and Illumina sequencing were carried out by the Biomarker Technologies Co., Ltd (Beijing, China). Three biological replicates were performed for each sample. A total of 6 RNA-Sequencing libraries were established. The clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS (Illumia) according to the manufacturer’s instructions. After cluster generation, the library preparations were sequenced on an Illumina Hiseq 2000 platform and paired-end reads were generated. Transcriptome data analysis Raw reads of FASTQ format were firstly processed through in-house Perl scripts. In this step, clean reads were obtained by removing reads containing adapter, reads containing ploy-N and low-quality reads from raw data. At the same time, Q20, Q30, GC-content and sequence duplication level of the clean reads were calculated. All the downstream analyses were based on clean reads with high quality. Transcriptome de novo assembly was accomplished using Trinity [[61]24]. The left files (read1 files) from three biological replicates samples were pooled into one big left.fq file, and right files (read2 files) into one big right.fq file. Transcriptome assembly was accomplished based on the left.fq and right.fq using Trinity with min_kmer_cov set to 2 by default and all other parameters set default. The expression of all genes in each sample was calculated using RSEM [[62]25]. Gene function was annotated based on the following databases: NR (NCBI non-redundant protein sequences); Pfam (Protein family); KOG/COG/eggNOG (Clusters of Orthologous Groups of proteins); Swiss-Prot (A manually annotated and reviewed protein sequence database); KEGG (Kyoto Encyclopedia of Genes and Genomes); GO (Gene Ontology). Differentially expressed gene (DEG) analysis Differential expression analysis of two conditions was performed using the DESeq R package (1.10.1). The gene expression differences with FDR < 0.05,|log[2]FC| ≥ 1 were considered statistically significant. GO enrichment analysis of the DEGs was implemented by the topGO R packages-based Kolmogorov–Smirnov test. KEGG pathways analysis of the DEGs was implemented by the KOBAS [[63]26] software to test the statistical enrichment of differential expression genes. Quantitative real-time PCR (qPCR) analysis To validate the reliability of the transcriptome data, qPCR was random conducted on 9 up-regulated DEGs, 9 down-regulated DEGs, and 9 normal Unigenes. Under normal condition and osmotic stress treatment for 3 d, the total RNA was extracted from the leaves of G. przewalskii seedlings and treated by DNase I to remove gDNA using OminiPlant RNA Kit (DNase I) (CWBio Co., Ltd, Taizhou, China). The cDNA synthesis of the leaves of G. przewalskii seedlings under normal condition and osmotic stress treatment for 3 d were performed using the SweScript All-in-One RT SuperMix for qPCR (One-Step gDNA Remover) (Servicebio Biotechnology Co., Ltd, Wuhan, China). The 2×Universal Blue SYBR Green qPCR Master Mix used for qPCR reaction was bought from Servicebio Biotechnology Co., Ltd (Wuhan, China). The qPCR was conducted on the CFX96 Touch Real-Time PCR Detection System (Bio-Rad Laboratories Co., Ltd, Shanghai, China). The qPCR analysis was performed as previously described [[64]27]. The G. przewalskii β-Actin (GpActin2) gene was used as reference genes in the qPCR reactions. Under normal condition and osmotic stress treatment for 3 d, the mean of GpActin2 Ct from three biological replicates were 16.50 and 16.63, respectively, suggested the expression of GpActin2 was stable under different conditions. Sequences of primers used for qPCR analysis are listed in Supplemental Table [65]S1. Statistical analysis Statistical analysis of data was conducted using GraphPad Prism 10.0 according to one-way ANOVA analysis of variance followed by Duncan’s test. The significance level was shown above the bar graph in the figures * P ≤ 0.05, ** P ≤ 0.01, *** P ≤ 0.001, while ns represents non-significance. The correlation was conducted using GraphPad Prism 10.0 according to Pearson correlation coefficient. Results G. przewalskii is a good forage in the desert region The excellent forages often with the characteristics of high protein, high fat, low lignin and low fiber [[66]28, [67]29]. To evaluate the forage quality of G. przewalskii, its ash content, fat content, crude protein content, lignin content, crude fiber content, ADF content, and NDF content in twigs and leaves were measured. The results showed that the fat contents of twigs and leaves were the lowest compared with other forage quality related indicators, is 1.85% and 0.92% respectively. The ash content, crude protein content, lignin content, crude fiber content, ADF content, and NDF content in twigs were 10.61%, 5.68%, 7.08%, 21.23%, 42.16%, and 58.42%, respectively. The ash content and crude protein content in leaves were higher than that in twigs, were 20.39% and 11.96% respectively. However, the lignin content, crude fiber content, ADF content, and NDF content in leaves were all lower than that in twigs, were 2.40%, 17.51%, 14.29% and 20.26% respectively (Fig. [68]1). These data showed that the twigs and leaves of G. przewalskii have a relatively higher protein content, and lower lignin and fiber content, indicating that G. przewalskii could be utilized as a good forage in the desert region with the leaves exhibiting particularly high forage quality. The leaves of G. przewalskii have higher forage quality even in extremely dry natural environments, inferring it is more important for the forage quality and drought tolerance of G. przewalskii. Therefore, the leaves of G. przewalskii seedlings were selected for the subsequent physiological parameter measurements and RNA‑Sequencing under osmotic stress. Fig. 1. [69]Fig. 1 [70]Open in a new tab The forage quality measurements of twigs and leaves in G. przewalskii G. przewalskii has high tolerance to drought stress To mimic the effect of drought stress on the growth of G. przewalskii, the 6-weeks old seedlings of G. przewalskii were treated by 15% PEG 6000 for 1 d, 3 d, 5 d, and 7 d to induce osmotic stress, and the chlorophyll content, MDA content, H[2]O[2] content, and OFR content of G. przewalskii seedlings under different conditions were subsequently measured. Under 1 d-osmotic stress, the G. przewalskii seedlings had not significantly change in its growth. However, when G. przewalskii seedlings were treated by osmotic stress for 3 d, 5 d, and 7 d, its growth became worse with the prolongation of osmotic stress, which could verify by the yellowing and wilting in the G. przewalskii seedlings (Fig. [71]2). Fig. 2. [72]Fig. 2 [73]Open in a new tab The G. przewalskii seedlings were treated by 15% PEG 6000 for different time. CK denoted the G. przewalskii seedlings under normal condition, which was Murashige and Skoog liquid medium at 25 °C under 60% relative humidity and a 16-hours light/8-hours dark photoperiod. 1d, 3d, 5d, 7d denoted the G. przewalskii seedlings treated by 15% PEG 6000 for 1 day, 3 days, 5 days, and 7 days In addition, the osmotic stress had caused physiological changes in the G. przewalskii seedlings. As shown in Fig. [74]3, the chlorophyll content of the G. przewalskii seedlings under 7 d-osmotic stress was significantly lower than that under normal condition. However, there were no significant differences among the chlorophyll contents of the G. przewalskii seedlings under normal condition, 1 d-osmotic stress, 3 d-osmotic stress, and 5 d-osmotic stress, indicating that G. przewalskii seedling could tolerate longer periods of drought. The MDA contents of G. przewalskii seedlings showed a rising trend under different osmotic stress treatments. The H[2]O[2] contents and OFR contents of G. przewalskii seedlings had not remarkably difference under different osmotic stress treatments, except for a significant reduction in OFR contents under the 7 d-osmotic stress, inferring that G. przewalskii seedlings suffered less secondary oxidative damage under osmotic stress. These data suggested that although drought stress affects the growth of G. przewalskii, it still had a high drought resistance and can tolerate long-term drought stress. Fig. 3. [75]Fig. 3 [76]Open in a new tab Measurement of chlorophyll content, MDA content, H[2]O[2] content, and OFR content of G. przewalskii seedlings under different conditions. CK denoted the G. przewalskii seedlings under normal condition, which was Murashige and Skoog liquid medium at 25 °C under 60% relative humidity and a 16-hours light/8-hours dark photoperiod. 1d, 3d, 5d, 7d denoted the G. przewalskii seedlings treated by 15% PEG 6000 for 1 day, 3 days, 5 days, and 7 days. FW denoted the fresh weight of G. przewalskii seedling The proline and antioxidase involved in G. przewalskii response to osmotic stress To explore the roles of Pro and antioxidase in G. przewalskii response to osmotic stress, the Pro content, SOD activity, CAT activity, and ·OH scavenging capacity of G. przewalskii under different conditions were measured. In contrast with normal condition, the Pro contents of G. przewalskii under osmotic stresses were significantly elevated except under 1 d-osmotic stress. The SOD and CAT activities of G. przewalskii increased under osmotic stresses, reaching peak levels at 7 d of treatment. Meanwhile, the ·OH scavenging capacity of G. przewalskii under osmotic stresses also higher than normal condition (Fig. [77]4). These results demonstrated that Pro, SOD, CAT played an important role in G. przewalskii response to osmotic stress. Fig. 4. [78]Fig. 4 [79]Open in a new tab Measurements of proline content, SOD activity, CAT activity, and ·OH scavenging capacity of G. przewalskii seedlings under different conditions. CK denoted the G. przewalskii seedlings under normal condition, which was Murashige and Skoog liquid medium at 25 °C under 60% relative humidity and a 16-hours light/8-hours dark photoperiod. 1d, 3d, 5d, 7d denotd the G. przewalskii seedlings were treated by 15% PEG 6000 for 1 day, 3 days, 5 days, and 7 days. FW denoted the fresh weight of G. przewalskii seedling RNA-sequencing and assembly Since yellowing and wilting of the G. przewalskii seedlings began after osmotic stress treatment for 3 d, the G. przewalskii seedlings under normal condition and 3 d-osmotic stress were selected to RNA-sequencing. In total, 44.51 clean reads that more than 6.87 Gb per sample were obtained from the RNA-sequencing of G. przewalskii seedlings. The Q20 and Q30 base percentage of each library was > 98% and 94%, and the GC-content of each library was between 44.56% and 45.51% (Table [80]1), indicating that clean reads were obtained and could be used for further analysis. All original FASTQ data files were submitted to the NCBI Sequence Read Archive (SRA) under accession numbers PRJNA1178850. Table 1. Overview of RNA-sequencing Sample Read Number Base Number Q20 (%) Q30 (%) GC content CK1 27,619,668 8,229,615,016 98.33 95.13 44.64% CK2 24,096,410 7,197,556,800 98.33 95.09 45.02% CK3 24,833,679 7,414,634,528 98.30 94.96 45.51% D1 22,954,357 6,870,537,788 98.06 94.03 44.70% D2 24,243,334 7,246,282,412 98.27 94.92 44.94% D3 25,319,502 7,547,457,568 98.30 94.95 44.56% [81]Open in a new tab CK1, CK2 and CK3 respectively denoted 3 biological replicates of the G. przewalskii seedlings under normal condition. D1, D2 and D3 respectively denoted 3 biological replicates of the G. przewalskii seedlings under 3 d-osmotic stress. Read Number denoted the pair-end reads number of clean reads. Base Number denoted the base number of clean reads Using Trinity to de novo assembly, a total of 328,642 transcripts were obtained with an average length of 1488.66 bp and a N50 of 2,310 bp. Taking the longest transcripts as Unigene, there were 102,191 Unigenes were obtained with an average length of 785.15 bp and a N50 of 1,494 bp (Table [82]2). The number of Unigene with length of 200–300 bp, 300–500 bp, 500–1000 bp, 1000–2000 bp and longer than 2000 bp was 35,673, 25,958, 18,910, 11,496, 10,154, respectively (Table [83]3). Table 2. The result of De novo assembly Term Total Number N50 (bp) Total length (bp) Mean Length (bp) Transcript 328,642 2,310 489,235,152 1488.66 Unigene 102,191 1,494 80,235,308 785.15 [84]Open in a new tab Table 3. Transcript and Unigene length distribution Length Range 200–300 bp 300–500 bp 500–1000 bp 1000–2000 bp > 2000 bp Transcript 45,947 43,440 60,053 87,057 92,145 Unigene 35,673 25,958 18,910 11,496 10,154 [85]Open in a new tab Unigene annotation The 102,191 Unigenes showed varying success annotation rates against nine databases: GO (22.67%), KEGG (18.09%), KOG (15.68%), COG (6.28%), Pfam (18.21%), Swiss-Prot (17.12%), TrEMBL (27.84%), eggNOG (23.90%), and NR (28.47%). There were 30,809 Unigenes annotated successfully in at least one of the nine databases, accounting for 30.15% of the total number of Unigenes (Table [86]4). Table 4. Annotation of Unigenes in various databases Database Number of Unigene Ratio (%) 300 ≤ length<1000 length ≥ 1000 Annotated in GO 23,170 22.67 6,599 12,928 Annotated in KEGG 18,484 18.09 4,885 11,026 Annotated in KOG 16,023 15.68 4,071 9,150 Annotated in COG 6,417 6.28 1,033 4,901 Annotated in Pfam 18,611 18.21 4,207 12,388 Annotated in Swiss-Prot 17,500 17.12 4,394 11,050 Annotated in TrEMBL 28,448 27.84 8,420 15,596 Annotated in eggNOG 24,428 23.90 7,117 13,594 Annotated in NR 29,092 28.47 8,729 15,668 Annotated in at least one database 30,809 30.15 9,538 15,844 [87]Open in a new tab Identification and analysis of DEGs in G. przewalskii response to osmotic stress Compared to normal condition, a total of 3,015 DEGs were identified in 3 d-osmotic stress, including 1,637 up-regulated and 1,378 down-regulated DEGs (Fig. [88]5). Among these DEGs, there were 2,607 DEGs were annotated in at least one of the nine databases, including GO, KEGG, KOG, COG, Pfam, Swiss-Prot, TrEMBL, eggNOG, and NR (Table [89]5). However, there were 408 DEGs had no BLAST hits in these nine databases, inferring that these DEGs may be the unique drought-resistant genes of G. przewalskii and they are worth further study in the subsequent experiments. Fig. 5. [90]Fig. 5 [91]Open in a new tab The volcan plot of DEGs in the transcriptome of G. przewalskii. DEGs were determined using a threshold of FDR < 0.05 and|log2FC| ≥ 1. Red spots represented up-regulated DEGs. Green spots indicated down-regulated DEGs. Black spots indicated normal Unigenes Table 5. Annotation of DEGs in various databases Term Annotated GO KEGG KOG COG Pfam Swiss-Prot TrEMBL eggNOG NR DEGs 2,607 2,134 1,739 1,193 811 2,053 1,966 2,571 2,226 2,585 Ratio(%) 86.47 70.78 57.68 39.57 26.90 68.09 65.21 85.27 73.83 85.74 [92]Open in a new tab GO enrichment in DEGs To identify the major biological processes that are expressed under osmotic stress, we performed a GO functional enrichment analysis (P < 0.05) on the all DEGs and other Unigenes. The analysis shown that there were 2,134 DEGs enriched in 3 GO primary categories and 47 GO secondary categories, including 1,130 up-regulated DEGs and 1,004 down-regulated DEGs (Table [93]S2, Fig. [94]6). The 3 GO primary categories were biological process, cellular component, and molecular function. In biological process, “metabolic process” (655 DEGs), “cellular process” (650 DEGs), “single-organism process” (507 DEGs), “response to stimulus” (232 DEGs) and “biological regulation” (231 DEGs) were significantly enriched among DEGs compared with the whole transcriptome background. In cellular component, the most significant secondary categories are related to “membrane” (834 DEGs), “membrane part” (746 DEGs), “cell” (570 DEGs), “cell part” (570 DEGs) and “organelle” (415 DEGs). In molecular function, the prime secondary categories DEGs are associated with the “catalytic activity” (1042 DEGs), “binding” (974 DEGs), “transporter activity” (135 DEGs), “nucleic acid binding transcription factor activity” (93 DEGs) and “structural molecule activity” (18 DEGs) (Table [95]S2, Fig. [96]6). Fig. 6. [97]Fig. 6 [98]Open in a new tab GO enrichment in DEGs. Dark column: the number and percentage of DEGs enriched. Light column: the number and percentage of all Unigenes enriched. X-axis: the names of the primary and secondary categories in GO enrichment. Y-axis (left): the percentage of Unigene numbers. Y-axis (right): above is the number of DEGs, below is the number of all Unigenes To obtain a further comprehension of the roles of these DEGs, we spectively performed a GO functional enrichment analysis (with a significance threshold of P < 0.05) on the 1,130 up-regulated DEGs and 1,004 down-regulated DEGs. Nevertheless, there were less difference in the GO enrichment of up-regulated DEGs and down-regulated DEGs, as evidenced by that the most significant secondary categories of up-regulated DEGs were almost same to down-regulated DEGs and the numbers of up-regulated DEGs and down-regulated DEGs in the different most significant secondary categories were also subequal (Table [99]S2). KEGG pathway enrichment analysis in DEGs The KEGG pathway enrichment analysis showed that a total of 1,739 DEGs were enriched in 129 KEGG pathway (P < 0.05) (Fig. [100]7A). Among these pathways, “Plant hormone signal transduction” (87 DEGs), “Plant-pathogen interaction” (82 DEGs), “MAPK signaling pathway-plant” (64 DEGs), “Phenylpropanoid biosynthesis” (46 DEGs), and “Starch and sucrose metabolism” (39 DEGs) exhibited significant enrichment (Fig. [101]7B). However, there was a little difference between the KEGG pathway enrichments of up-regulated DEGs and down-regulated DEGs. The up-regulated DEGs were primarily enriched in “Plant-pathogen interaction” (50 DEGs), “Plant hormone signal transduction” (36 DEGs), “MAPK signaling pathway-plant” (35 DEGs), “Phenylpropanoid biosynthesis” (23 DEGs), and “Spliceosome” (21 DEGs) (Fig. [102]7C). The down-regulated DEGs were primarily enriched in “Plant hormone signal transduction” (51 DEGs), “Plant-pathogen interaction” (32 DEGs), “MAPK signaling pathway-plant” (29 DEGs), “Starch and sucrose metabolism” (26 DEGs), and “Pentose and glucuronate interconversions” (25 DEGs) (Fig. [103]7D). Fig. 7. [104]Fig. 7 [105]Open in a new tab KEGG pathway enrichment of DEGs. A: Statistical graph of KEGG number. B: The top 20 KEGG-enriched pathways of all DEGs. C: The top 20 KEGG-enriched pathways of up-regulated DEGs. D: The top 20 KEGG-enriched pathways of down-regulated DEGs DEGs related to proline To explore the accumulation mechanism of proline in G. przewalskii response to osmotic stress, we analyzed the Unigenes related to pyrroline-5-carboxylate synthetase (P5CS) protien and proline dehydrogenase (PRODH) protien, which were 2 key enzymes in the metabolism of proline. The results shown that 2 Unigenes were annotated to PRODH and 3 Unigenes were annotated to P5CS. Among these 5 Unigenes, 2 Unigenes related to P5CS were up-regulated DEGs, while the expression of other Unigenes did not change significantly under osmotic stress (Table [106]6). These results demonstrated that G. przewalskii increased the accumulation of proline mainly by promoting the biosynthesis of proline under osmotic stress. Table 6. DEGs related to proline Gene CK-FPKM D-FPKM FDR log2FC Swiss-Prot Annotation c74475.graph_c0 68.07 250.68 5.8 × 10^− 5 1.67 Delta-1-pyrroline-5-carboxylate synthase OS = Mesembryanthemum crystallinum OX = 3544 GN = P5CS PE = 2 SV = 1 c58092.graph_c0 8.27 42.1 6.7 × 10^− 5 2.08 Delta-1-pyrroline-5-carboxylate synthase OS = Mesembryanthemum crystallinum OX = 3544 GN = P5CS PE = 2 SV = 1 c51258.graph_c0 0.79 1.55 Delta-1-pyrroline-5-carboxylate synthase B OS = Arabidopsis thaliana OX = 3702 GN = P5CSB PE = 2 SV = 1 c70093.graph_c0 79.61 170.17 0.1 1.01 Proline dehydrogenase 2, mitochondrial OS = Arabidopsis thaliana OX = 3702 GN = POX2 PE = 2 SV = 1 c93561.graph_c0 0 1.57 Proline dehydrogenase 2, mitochondrial OS = Arabidopsis thaliana OX = 3702 GN = POX2 PE = 2 SV = 1 [107]Open in a new tab CK-FPKM and D-FPKM respectively denoted the average FPKM of 3 biological replicates of the G. przewalskii seedlings under normal condition and 3 d-osmotic stress DEGs related to antioxidant system In G. przewalskii transcriptome, we had identified 32 DEGs related to various antioxidases and antioxidants, including glutaredoxin (GRX) (5 DEGs), ascorbate peroxidase (APX) (1 DEGs), monodehydroascorbate reductase (MDAR) (1 DEGs), glutathione S-transferases (GST) (10 DEGs), peroxidase (POD) (8 DEGs), peroxiredoxin (PRX) (1 DEGs), thioredoxin (TRX) (3 DEGs), SOD (2 DEGs) and CAT (1 DEGs) (Table S3). The expression patterns of DEGs related to antioxidant system were analyzed in G. przewalskii transcriptome. As shown in Fig. [108]8, The number of DEGs in GST was the largest, including 7 up-regulated DEGs and 3 down-regulated DEGs. In POD related DEGs, there were 3 up-regulated DEGs and 5 down-regulated DEGs. In GRX related DEGs, there were 3 up-regulated DEGs and 2 down-regulated DEGs. In TRX related DEGs, there were 2 up-regulated DEGs and 1 down-regulated DEGs. For APX、MDAR、CAT and PRX related DEGs, only the DEG encoding MDAR (Unigene ID: c74365.graph_c0) was up-regulated (Fig. [109]8). These data suggested that antioxidant system played an important role in G. przewalskii response to osmotic stress. Fig. 8. [110]Fig. 8 [111]Open in a new tab Heat map showing the expression patterns of DEGs related to antioxidant system. CK denoted normal condition. D denoted 3 d-drough stress. The colour bar denoted log[2]FPKM values Plant hormone involved in G. przewalskii response to osmotic stress To understand the role of plant hormone playing in G. przewalskii response to osmotic stress, we analyzed the DEGs involved in plant hormone signaling pathways. The results showed that DEGs were annotated to 8 plant hormone signaling pathways: auxin (21 DEGs), cytokinin (6 DEGs), gibberellin (GA) (10 DEGs), abscisic acid (ABA) (16 DEGs), ethylene (5 DEGs), brassinosteroid (15 DEGs), jasmonic acid (JA) (5 DEGs) and salicylic acid (SA) (3 DEGs) (Table S4, Fig. [112]9). Among the 8 plant hormone signaling pathways, auxin and ABA signaling pathways exhibited dominant enrichment. The above data indicated that 8 plant hormones were participated in G. przewalskii response to osmotic stress, with auxin and ABA playing a dominant regulatory role in these plant hormones. Fig. 9. [113]Fig. 9 [114]Open in a new tab The plant hormone signal in G. przewalskii response to drought stress. Green box denoted up-regulated DEGs. Red box denoted down-regulated DEGs. Blue box denoted containing both up-regulated and down-regulated DEGs Classification of transcription factors in G. przewalskii response to osmotic stress Transcription factors (TFs) is a crucial component of stress signalling response in the plant because of regulating the expression of downstream signalling related genes. In G. przewalskii transcriptome, 918 Unigenes were identified to TFs including 150 DEGs and belonging to 35 families (Table S5). Among the 150 DEGs encoding TFs, the abundant enrichment of TF families were WRKY (16 DEGs), bZIP (14 DEGs), ERF (14 DEGs), bHLH (11 DEGs), MYB (10 DEGs), and NAC (9 DEGs) (Fig. [115]10). Fig. 10. [116]Fig. 10 [117]Open in a new tab The number of DEGs encoding TFs. Red columns indicated up-regulated DEGs. Green columns indicated down-regulated DEGs To response drought stress, the expression of plant TF genes must make a timely and correct changes. In the present study, the expression patterns 150 DEGs encoding TFs were analyzed in G. przewalskii transcriptome. As shown in Fig. [118]11, a total of 90 up-regulated DEGs and 60 down-regulated DEGs were identified in the 150 DEGs encoding TFs. In addition, the different DEGs encoding TFs shown various expression patterns, even though the DEGs belong to same family. These results suggested that TFs have a complex regulated network in G. przewalskii response to osmotic stress. Fig. 11. [119]Fig. 11 [120]Open in a new tab Heat map showing the expression patterns of DEGs encoding transcription factors. CK denoted normal condition. D denoted 3 d-drough stress. The colour bar denoted log[2]FPKM values Validation of DEGs by qRT‑PCR To validate the RNA-seq results, we randomly selected 9 up-regulated DEGs, 9 down-regulated DEGs, and 9 normal Unigenes for qRT-PCR to confirm the reliability of the RNA-Seq results. The results showed that the relative expression levels of the 9 up-regulated DEGs were all greater than 1, the relative expression levels of the 9 down-regulated DEGs were all less than − 1, and the relative expression levels of the 9 normal unigenes were all almost equal to 0, indicating that the 9 up-regulated DEGs and 9 down-regulated DEGs reliably changed significantly under osmotic stress, but the 9 normal Unigenes unchanged (Fig. [121]S1). The Pearson correlation coefficient between RNA-seq data and qRT‑PCR data was r = 0.9202 (P < 0.001). The qPCR results for 27 selected Unigenes showed general agreement with their transcript abundance changes determined by RNA-seq, which indicated that the RNA-Seq results used in this study were reliable. Discussion G. przewalskii is an important forage with strong drought resistance in the desert region Desert and semi-desert are important ecosystem of the planet. When plants under deserts, the most threatening abiotic stress they encounter is drought. After a long period of natural selection and biological evolution, some typical desert plants have gradually evolved to be extremely resistant to drought in desert and semi-desert, such as Ammopiptanthus mongolicus, Zygophyllum xanthoxylon, Haloxylon ammodendron and Haloxylon persicum [[122]30–[123]32]. These typical desert plants contain a large number of abiotic stress related genes and are excellent materials for understanding the mechanisms of plant response to drought stress. For example, a non-targeted metabolomics analysis demonstrated that Haloxylon ammodendron responds to drought by increasing the content of organic nitrogen compounds and lignans, neolignans, and related compounds, and reducing the content of alkaloids and derivatives. By contrast, Haloxylon persicum adapts to the dry environment by increasing the content of organic acids and their derivatives and reducing the content of lignans, neolignans, and related compounds [[124]30]. The survival strategies of Ammopiptanthus mongolicus and Zygophyllum xanthoxylon in saline and drought environments showed that seed germination of Ammopiptanthus mongolicus and Zygophyllum xanthoxylon was negatively affected by combined salt and drought stresses. Compared with Zygophyllum xanthoxylon, the higher germination rate and shorter radicle were detected for Ammopiptanthus mongolicus, indicating that Ammopiptanthus mongolicus had less resistance to salt and drought stresses than Zygophyllum xanthoxylon [[125]31]. Furthermore, AmDREB3 [[126]33], ZxZF [[127]34], HaNAC3 [[128]35], ScPIP1 [[129]36], and RtWRKY23 [[130]37], the genes isolated from typical desert plants, were demonstrated to improve the tolerance of plants to abiotic stress. G. przewalskii is a typical desert plant which originated in Tethys Ocean [[131]38]. In the present study, the G. przewalskii seedlings were still alive in 5 d-osmotic stress with a little of yellowing and wilting of the leaves (Fig. [132]2). The chlorophyll, H[2]O[2] and OFR contents of G. przewalskii seedlings had not remarkably difference under 1d, 3d and 5d osmotic stress treatment, inferring that G. przewalskii seedlings suffered lesser oxidative damage and maintained better photosynthesis under 1d, 3d and 5d osmotic stress treatment (Fig. [133]3). Furthermore, there were 3,015 DEGs obtained in the transcriptome of G. przewalskii under osmotic stress (Fig. [134]5). These results demonstrated that G. przewalskii had strong tolerance to drought stress and rich resources of drought stress resistance related genes, similar to other typical desert plants. In addition, one of the important factors affecting the ecological balance and economic development of desert and semi-desert is the availability of forage for herbivores. Previous studies have shown that excellent forages are high in protein and digestibility [[135]29]. However, lignin and fiber contents are inversely related to forage quality due to their higher content reduces the digestibility and energy value of forage for feeding livestock [[136]29, [137]39]. For example, the reductions in NDF, ADF and lignin contents leaded to a significant increase in the forage quality of Medicago sativa, which is consistent with the observed increase in digestibility [[138]40]. Unfortunately, most herbages with good forage quality, such as Fabaceaeand and Poaceae herbages, are difficult to grow and reproduce in extremely dry desert and semi-desert. Therefore, the special forage with strong drought resistance is important for improving economic development in desert and semi-desert regions. On the previous field study, we had found that G. przewalskii grew well in Alxa desert. However, the twigs and leaves of G. przewalskii were more severely chewed by herbivores in Alxa desert than other surrounding plants, such as Zygophyllum xanthoxylon and Reaumuria songarica. These results inferred that G. przewalskii had a high adaptation to desert and high palatability to herbivores (Fig. [139]S2). In the present study, the crude protein content of twigs and leaves in G. przewalskii was determined to be 7.08% and 11.96%, respectively. Furthermore, the lignin and fibre content was found to be low in the twigs and leaves of G. przewalskii. The results suggested that the twigs and leaves of G. przewalskii all possessed high nutrient value, with the leaves exhibiting particularly high forage quality (Fig. [140]1). Furthermore, the seedlings of G. przewalskii shown high tolerance to drought stress because of survival from the osmotic stress for a long time (Fig. [141]2). Therefore, G. przewalskii is a drought-tolerant forage which could be used to explore the mechanism of plants response to drought stress and provide the necessary forage for herbivores in semi-desert and desert. Previous studies suggested that climate change affected the growth and metabolism of forage, which in turn affected their aboveground net primary production and nutritional quality [[142]41, [143]42]. For example, the linear change rates of potential and actual forage crude protein, ether extract, Ash, NDF, ADF and water-soluble carbohydrates contents showed significant nonlinear correlations with annual temperature and annual precipitation in alpine grasslands of the Tibet in 2000–2020 [[144]43]. The comprehensive meta-analysis shown that the co-occurrence of drought stress and elevated CO[2] negatively impacted on the productivity and quality of Fabaceaeand and Poaceae herbages by decreasing crude protein content, nitrogen content and forage dry matter yield [[145]44]. In this study, the leaves of G. przewalski used for forage quality measurements were collected from the natural environment of Yabalai Mountain in Alxa desert, which is an extremely dry natural environment. Nevertheless, the leaves of G. przewalski still had a very high forage quality with 11.96% crude protein content (Fig. [146]1). Therefore, we inferred that the excellent drought tolerance made G. przewalski to survival in desert and maintain a high forage quality. However, how the drought stress effect on the forage quality of G. przewalski remains unclear, which is worth exploring in future studies. Construction of an informative transcriptome dataset for G. Przewalskii The genes expression is the original and crucial process in plants response to drought stress and RNA-sequencing is an effective means of detecting the changes of genes expression [[147]45]. An informative transcriptome dataset is the necessity to explore the mechanism of plants response to drought stress. However, the transcriptome data of G. przewalskii under drought stress are still lacking. In the present study, the leaves of G. przewalskii seedlings under normal condition and 3 d-osmotic stress treatment were selected for RNA-sequencing. A total of 102,191 Unigenes were obtained in G. przewalskii with an average length of 785.15 bp and a N50 of 1,494 bp (Table [148]2). Among the total Unigenes of G. przewalskii, there were 30,809 Unigenes were successfully annotated by BLAST analysis and functional bioinformatics analyses (e.g., GO, KEGG, and Swiss-Prot) (Table [149]4). These results suggested that the sequences of the G. przewalskii Unigenes generated in this study were assembled and annotated correctly and this transcriptome will contribute to uncover the mechanism of G. przewalskii response to drought stress. DEGs are crucial in transcriptome, which represent the significant changes of gene expression in plants [[150]45]. In this study, a total of 3,015 DEGs were identified from the transcriptome of G. przewalskii (Fig. [151]5). In addition, 2,607 DEGs were annotated by BLAST analysis and functional bioinformatics analyses (e.g., GO, KEGG, and Swiss-Prot) (Table [152]5). The DEGs functional bioinformatics analysis results showed that “catalytic activity”, “binding” and “membrane” were the prime enriched GO terms (Fig. [153]6), and “Plant hormone signal transduction”, “Plant-pathogen interaction”, and “MAPK signaling pathway-plant” were the main enriched KEGG pathway (Fig. [154]7A), which provided a new insight into the mechanism of G. przewalskii response to drought stress. However, there were 408 DEGs had no significant BLAST hits in the public databases, inferring that these unannotated DEGs represent a specific drought resistant gene resource of G. przewalskii which warrants further investigation. The proline and antioxidant system maintain the homeostasis of G. przewalskii under drought stress The water deficit and excessive ROS caused by drought stress will destroy the homeostasis and inhibite the growth of plant. Reducing cell water potential to enhance the ability of absorb and maintain water in plant is one of the important ways in which plants resistance to water deficit [[155]46]. To reduce cell water potential, plants have evolved a variety of osmolytes to reduce osmotic potential, and one of the most important osmolytes is proline. Numerous studies have demonstrated that the accumulation of proline significantly promotes plant adaptation to drought stress [[156]47]. In present study, the proline content of G. przewalskii was significantly elevated after osmotic stress treated for 3d, 5d, and 7d (Fig. [157]4), suggesting that proline plays an important role in G. przewalskii response to osmotic stress. Furthermore, 2 P5CS related Unigenes were significantly up-regulated (FDR < 0.05,|log2FC| ≥ 1) (Table [158]6). The P5CS is a key enzyme of the proline biosynthesis in plants [[159]37]. These date suggested that G. przewalskii could elevate the expressions of P5CS related genes and promte the proline biosynthesis to maintain the water homeostasis in G. przewalskii under drought stress, which significantly enhanced the tolerance of G. przewalskii to drought stress. ROS are important signaling molecules in plant, which govern the regular growth, development, and response to abiotic stress [[160]48]. However, the excessive ROS caused by abiotic stress will disrupte the ROS homeostasis and destroy to the cell membrane system in plants. To maintain the ROS homeostasis, plant could elevate the activities of antioxidases to scavenge the excessive ROS caused by external abiotic stress [[161]6, [162]49]. In this study, we had found that G. przewalskii seedlings suffered less oxidative damage under osmotic stress, which inferred from the relatively stable and low MDA, H[2]O[2], and OFR contents in G. przewalskii under osmotic stress (Fig. [163]3). The SOD and CAT activities and ·OH scavenging capacity of G. przewalskii under osmotic stress were significantly increased (Fig. [164]4). Accordingly, 32 DEGs related to various antioxidases and antioxidants were identified in G. przewalskii transcriptome (Fig. [165]8). Based on these data, it was hypothesized that G. przewalskii had an efficient antioxidant system and this system helped G. przewalskii maintain the ROS homeostasis resulting in G. przewalskii suffered less oxidative damage under drought stress. Interestingly, in G. przewalskii transcriptome, nearly half of the antioxidases and antioxidants related DEGs were annotated in GRX (5 DEGs) and GST (10 DEGs). Consequently, glutathione may play an important role in the antioxidant system of G. przewalskii and worth in-depth study. The plant hormone and TFs play a key regulatory function in G. przewalskii response to drought stress When plants are exposed to drought stress, TFs and plant hormone signaling pathways will be integrated a networks to timely regulate the gene expression and physiological processes in response to the abiotic stresses [[166]46]. In G. przewalskii, we discovered that the most enriched KEGG pathway of DEGs was “Plant hormone signal transduction” (87 DEGs) associating with auxin, cytokinin, GA, ABA, ethylene, brassinosteroid, JA and SA mediated signaling pathways (Fig. [167]9), and150 DEGs were identified to 35 TF families including WRKY, bZIP, ERF, bHLH, MYB, and NAC TF families (Fig. [168]10). These results suggested that TFs and plant hormone signaling pathway orchestrated an extremely complex regulating network to regulate G. przewalskii response to drought stress. In G. przewalskii, we have found that the auxin (21 DEGs) and ABA (16 DEGs) signal transduction pathways exhibited dominant enrichment in KEGG pathway enrichment analysis (Table S4), inferring that auxin and ABA signal transduction pathways involved in G. przewalskii response to drought stress. Many studies demonstrated that auxin signal transduction pathways not only promote plant growth and development, but also take part in regulating plant response to drought stress [[169]50]. Auxin mediates plant growth and drought stress responses via a short nuclear pathway that can activate and/or repress various responsive genes by recruiting DNA-binding transcription factors: the AUXIN RESPONSE FACTOR (ARF) proteins. The AUXIN INFLUX CARRIER (AUX1) proteins can transport auxin to different parts of the plant, while the AUXIN/IAA proteins were repressors of auxin signal transduction pathways by modulating the activity of ARFs [[170]51, [171]52]. In our study, the DEGs related to AUX1 and ARF were all down-regulated, but 1 DEG related to AUXIN/IAA (Unigene c70163.graph_c0) were up-regulated (Table S4, Fig. [172]9), inferring that G. przewalskii may inhibit the growth under drought stress by inhibiting the auxin transport and auxin signal pathways transduction. ABA is generally considered to be a stress-hormone due to it involve in plants response to various abiotic stresses. The ABA signal transduction pathways could induce the stomatal closure of plants to decrease the water loss in plants under drought and its related genes, such as PYL, PP2C, SnRK, and ABF, were responsive to drought [[173]6, [174]53]. In current study, 1 DEGs related to ABA RESPONSIVE ELEMENT BINDING FACTOR (ABF) (Unigene c67965.graph_c0) (Table S4, Fig. [175]9), the crucial DNA-binding transcription factor of ABA signal transduction pathways, was up-regulated, inferring that G. przewalskii may regulated the transduction of ABA signal pathways to response to drought. However, there were 7 DEGs related to the negative regulator PP2Cs were up-regulated but all DEGs related to the ABA receptor PYR/PYL were down-regulated (Table S4, Fig. [176]9). The similar results that down-regulation of receptor PYR/PYL, up-regulation of negative regulator PP2Cs, and up-regulation of the downstream transcription factor ABF under drought stress were found in Medicago ruthenica and Casuarina equisetifolia, which demonstrated that the negative feedback of ABA signal transduction pathways could reduce the adverse effects caused by ABA accumulation and increase the tolerance to drought stress in plants [[177]5, [178]7]. Thus, we hypothesized that G. przewalskii regulated the negative feedback of ABA signal transduction pathways to response drought stress. Previous studies have proved that plants can use more material and energy in response to abiotic stresses by inhibiting growth and development [[179]54]. Based on the expression patterns of auxin and ABA signal transduction pathways related DEGs, we inferred that G. przewalskii may reduce the growth by regulating auxin signal pathway and promote the stomatal closure by regulating ABA signal pathway in response to drought stress, which maintained the balance of material distribution between growth and drought stress tolerance. This balance helped G. przewalskii to maximize survival in drought stress. WRKY TF family was one of the most important plant-specific TF families in plants. So far, lots of WRKY TFs have been found to play a key regulating function in plant responses to drought stress by regulating different physiological processes, such as plant hormone signal transduction, photosynthesis, ROS scavenging, root development, and secondary metabolism [[180]54]. For example, Dioscorea composita WRKY5 enhance drought and salt tolerances in by regulating the expressions of AtSOD1 and AtABF2 [[181]49]. In Arabidopsis, WRKY57 and WRKY63 confer drought tolerance to plants by regulating the transduction of ABA signaling pathway [[182]55, [183]56]. In our study, there were 14 up-regulated DEGs and 2 down-regulated DEGs identified to WRKY TFs in G. przewalskii transcriptome, resulting in WRKY TF family was the largest TF family in the DEGs of G. przewalskii transcriptome. These results suggested that WRKY TF family may mainly act as positive regulatory factor to particpate in G. przewalskii response to drought stress. Conclusion In conclusion, we have evaluated the forage quality of G. przewalskii and explored the mechanism of G. przewalskii response to drought stress. The twigs and leaves of G. przewalskii have a good forage quality and a higher palatability to herbivores compared with other plants around it. A total of 102,191 Unigenes were de novo assembled and 30,809 Unigenes were annotated successfully by BLAST analysis and functional bioinformatics analyses in G. przewalskii under osmotic stress. The analyses of DEGs and physiological indexes revealed that G. przewalskii regulated ROS scavenging, proline biosynthesis, and other complex physiological processes to response drought tress, in which auxin, ABA, and WRKY might play a key regulatory role. These results provided new insight into the mechanism of adaption to drought tolerance and laid the foundation to protection and utilization of G. przewalskii. Electronic supplementary material Below is the link to the electronic supplementary material. [184]Supplementary Material 1^ (750.8KB, docx) [185]Supplementary Material 2^ (58.3KB, docx) Acknowledgements