Integrated physiological and transcriptomic analysis reveals key drought tolerance mechanisms in Festuca rubra

In recent years, droughts caused by global climate change have become more frequent and intense, posing severe challenges to agricultural production and ecosystems^[1]. Arid regions cover around 41% of the Earth's land surface and support over 38% of the global population, thereby playing a vital role in maintaining the stability of terrestrial ecosystems^[2,3]. In grassland ecosystems and urban landscapes, grassland degradation caused by drought stress has become a more prominent problem, resulting in annual economic losses amounting to billions of dollars^[4]. Given these challenges, the cultivation and screening of drought-resistant grass seeds have become some of the most important areas of current agricultural research.

Festuca rubra L., a perennial herbaceous plant belonging to the Poaceae family, is widely used as a high-quality lawn grass and forage in temperate regions because of its notable grazing tolerance, prolonged green period, exceptional stress resistance, and high degree of adaptability^[5−7]. Compared with commonly cultivated species such as Festuca arundinacea and Lolium perenne, F. rubra exhibits superior drought tolerance and water use efficiency^[8]. As reported by Taleb et al.^[8], under identical drought conditions, F. rubra maintains a significantly lower decline in biomass than other common turfgrass species, highlighting its potential as an ideal choice for turf establishing in arid regions.

Under drought stress, plants activate a series of morphological, physiological, and biochemical response mechanisms^[9]. The plants' response to this stress involves the suppression of photosynthesis, restriction of growth, and reduction in dry matter accumulation^[10]. When antioxidant defenses are overwhelmed, the ensuing accumulation of reactive oxygen species (ROS) triggers morphological damage and reductions in photosynthetic electron transport efficiency^[11]. Plants' antioxidant systems are broadly categorized into enzymatic and nonenzymatic types, which function collectively to eliminate the accumulatinon of peroxide induced by drought^[12]. Simultaneously, plants accumulate a range of compatible solutes, such as proline, trehalose, fructose, and mannitol, to maintain cellular homeostasis under prolonged drought^[13]. The plant hormone abscisic acid (ABA) plays a pivotal role in mediating drought and other abiotic stress responses. It triggers transcriptional reprogramming, leading to adaptive changessuch as osmoprotectant enrichment and closure of the stomata^[14].

The advent of next-generation sequencing technology (NGS) has established RNA sequencing (RNA-Seq) as a potent instrument for the analysis of plants' gene expression profiles, particularly in drought tolerance research^[15]. In addition, various transcription factors, such as NAC, WRKY, and C2H2, are known to play an important role in regulating the expression of multiple genes in plants under drought stress^[16,17].

Although significant progress has been made in understanding the mechanisms of plants' responses to drought stress, including recent studies on related forage grasses like tall fescue (Festuca arundinacea) that have elucidated key physiological and genetic adaptations^[18,19], the specific mechanisms underlying the response of F. rubra to drought stress remain unclear, particularly its comprehensive regulatory network and key factors^[20]. Previous studies have mainly focused on crops or model plants, with insufficient attention paid to turfgrass, especially F. rubra^[21]. Moreover, most existing studies have been limited to single-level (physiological or molecular) analyses^[22], thereby lacking an integrated perspective. Therefore, this study aimed to systematically characterize the physiological and molecular changes in F. rubra under drought stress through an integrated analysis of physiological phenotypes and transcriptome dynamics. Our objectives are to elucidate its drought resistance mechanisms and to identify key differentially expressed genes (DEGs) and signaling pathways. This study will reveal the regulatory mechanisms of F. rubra under drought stress, providing an important theoretical basis and genetic resources for a deeper understanding of drought resistance in Festuca plants and for breeding drought-resistant turfgrass varieties.

Seeds of Festuca rubra L. 'Mengshen' were obtained from Beijing Zhengdao Seed Co., Ltd. The experiment was conducted in the greenhouse of the Plant Science Center at Beijing Forestry University. Uniform-sized seeds with an intact morphology were selected and sown in plastic pots filled with a 3:1:1 (v/v/v) mixture of nutrient soil, vermiculite, and perlite. Approximately 20–30 seeds were sown per pot at a depth of 1–2 cm. Plants were grown in a controlled environment chamber set to a 16-h photoperiod with a light intensity of 500 μmol·m⁻²·s⁻¹, and day/night temperatures of 25/20 °C. The substrate's moisture was maintained at approximately 70% of field capacity by watering every 2 d. After germination, seedlings were irrigated every 3 d with half-strength Hoagland nutrient solution. Approximately 2 months later, vigorous and uniform seedlings were thinned and retained for the stress experiment. The experimental design comprised a control group (CK) and a drought-stressed group (D), with each group consisting of 15 pots. The control group was regularly irrigated throughout the experiment to maintain normal growth conditions. Guided by both preliminary trials and existing studies^[23,24], the drought-stressed group was subjected to 21 consecutive days when water was withheld to simulate drought stress. After that, leaf samples from both the control and drought-stressed groups were collected. For each treatment, three biological replicates (0.3 g per sample) were collected. Each biological replicate comprised pooled material from five individual plants per treatment. All samples were immediately flash-frozen in liquid nitrogen and stored at −80 °C for subsequent physiological measurements and transcriptome sequencing.

Physiological and biochemical analysis methods

Relative electrical conductivity (REC) was determined by the conductivity method^[25]. Chlorophyll was extracted by the acetone method and quantified spectrophotometrically following Lichtenthaler^[26]. Chlorophyll fluorescence parameters, including initial fluorescence (F₀), maximum fluorescence (F_m) and nonphotochemical quenching (NPQ), were measured using a compact multifunctional plant phenotyping system (Plant Explorer) (Pheno Vation B.V.). On the basis of F₀ and F_m, we calculated the variable fluorescence (F_v = F_m − F₀) and the maximum quantum yield of Photosystem II (PSII) (F_v/F_m). Leaf photosynthetic characteristics including photosynthetic rate, stomatal conductance, intercellular CO₂ concentration (Ci), and transpiration rate were measured between 09:00 and 11:00 h under saturating light and stable temperature conditions using a LI-6400 photosynthesis system (Li-Cor, Inc., Lincoln, NE, USA). For all measurements, healthy, fully expanded leaves from plants with uniform growth were selected, with three to five biological replicates per treatment.

Osmotic and oxidative balance

Crude enzymes extracts were prepared from ground plant samples by homogenization in a 50 mM sodium phosphate buffer (pH 7.0) supplemented with 1% (w/v) polyvinylpyrrolidone (PVP) and 0.2 mM ethylenediaminetetraacetic acid disodium salt (EDTA-Na₂)^[27]. Following centrifugation, the resulting supernatant was utilized for subsequent enzymatic assays of catalase (CAT), superoxide dismutase (SOD), peroxidase (POD), and ascorbate peroxidase (APX).

CAT activity was assayed by monitoring the decomposition of H₂O₂ at 240 nm (A₂₄₀) for 60 s^[28]. The reaction mixture consisted of a 50 mM phosphate buffer (pH 7.0) and 150 mM H₂O₂. One unit of CAT activity was defined as a decrease in absorbance of 0.1 per minute. SOD activity was measured using the nitroblue tetrazolium (NBT) method^[29]. The reaction mixture, containing a 50 mM phosphate buffer (pH 7.8), 39 mM methionine, 1 mM EDTA-Na₂, 7.5 mM NBT, 0.1 mM riboflavin, and the enzyme extract, was illuminated at 4,000 lx for 20 min. The absorbance at 560 nm (A₅₆₀) was then measured. One unit of SOD activity was defined as the amount of enzyme causing 50% inhibition of NBT photoreduction. POD activity was determined via the guaiacol colorimetric method^[30]. The reaction mixture contained a 50 mM phosphate buffer (pH 7.0), 0.3% (v/v) H₂O₂, 0.2% (v/v) guaiacol, and the enzyme extract. The increase in absorbance at 470 nm (A₄₇₀) was for 1 min. One unit of POD activity was defined as an increase of 0.01 in absorbance per minute. APX activity was evaluated by tracking the oxidation of ascorbate at 290 nm (A₂₉₀)^[31]. The reaction system contained a 50 mM phosphate buffer (pH 7.0) with 50 mM ascorbate, 10 mM EDTA-Na₂, 10 mM H₂O₂, and the enzyme extract. One unit of APX activity was defined as a decrease in absorbance of 0.01 per minute.

Proline (Pro) content was determined using the sulfonylsalicylic acid method of Shen^[32]. Soluble sugar content was measured via the anthrone colorimetric method^[33]. Soluble protein content was quantified using the Bradford method (Coomassie Brilliant Blue)^[34]. Malondialdehyde (MDA) content was analyzed with the thiobarbituric acid (TBA) method^[35].

Transcriptome sequencing analysis

Total RNA was extracted from the leaves of F. rubra. The RNA's purity and concentration were assessed using a NanoDrop 2000 spectrophotometer, and the RNA's integrity was evaluated with an Agilent 2100 Bioanalyzer. After passing quality control, transcriptome libraries were constructed. Following library quality inspection, paired-end sequencing (150 bp) was performed on the Illumina NovaSeq 6000 platform. Transcriptome sequencing and data analysis was carried out by Shanghai OE Biotech Co., Ltd. To obtain high-quality reads, raw sequencing data were processed for quality filtering. Trimmomatic was used to remove adapters, trim low-quality bases, and eliminate reads containing ambiguous nucleotides (N bases), resulting in clean reads^[36]. The clean reads were then assembled into transcript sequences using Trinity (version 2.4) with a paired-end approach^[37]. On the basis of sequence similarity and length, the longest transcript was selected as a representative unigene for subsequent analysis. For differential expression analysis, unigenes with mean counts of ≤ 2 were filtered out. DESeq2 was used to normalize the raw count data (using the BaseMean value for quantifying expression), calculate the fold changes, and perform statistical significance tests via the negative binomial distribution (NB) model^[38]. DEGs were identified using a threshold of q < 0.05 (adjusted p-value) and fold change (FC) > 2.

Functional annotation of unigenes

Functional annotation of unigenes was performed using DIAMOND^[39] to align sequences against the Nonredundant Protein Sequence Database (NR), Clusters of Orthologous Groups for Eukaryotic Complete Genomes (KOG), Gene Ontology (GO), manually annotated protein database (Swiss-Prot), Evolutionary Genealogy of Genes: Nonsupervised Orthologous Groups (eggNOG), and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases. Additionally, HMMER^[40] was used to align sequences against the Pfam (Protein Family) database for domain-based functional characterization.

GO enrichment analysis of DEGs

After obtaining DEGs, GO enrichment analysis was performed using the GO database (http://geneontology.org) to functionally characterize the DEGs according to their annotations. The GO enrichment method involved counting the number of DEGs associated with each GO term and calculating the statistical significance of enrichment using the hypergeometric distribution algorithm. Fisher's exact test was applied to compute the p-value for each GO term within the three categories: biological process (BP), cellular component (CC), and molecular function (MF). Terms with lower p-values were considered to be more significant. DEGs were prioritized for further study by integrating the GO enrichment results with biological relevance. For visualization, the top 30 enriched GO terms (with PopHits ≥ 5 in BP, CC, or MF) were selected and ranked by –log₁₀ (p-value) (top 10 terms per category). Enrichment chord diagrams were generated to display the 10 most significant terms (lowest q-value or p-value).

KEGG enrichment analysis of DEGs

Pathway analysis of DEGs was conducted using the KEGG database (www.genome.jp/kegg), based on KEGG annotations. The significance of enrichment for each pathway was calculated via the hypergeometric distribution test. The top 20 enriched pathways (with PopHits ≥ 5) were filtered and ranked by –log₁₀ (p-value) in descending order.

Transcription factor database annotation

The distribution of transcription factors (TFs) across all genes and DEGs (upregulated and downregulated) was compared to identify TFs with significant proportional differences. TFs showing marked distribution shifts were further analyzed to determine their association with differential expression.

Weighted gene co-expression network analysis

To systematically identify key gene modules co-varying with physiological traits under drought stress, a weighted gene co-expression network analysis (WGCNA) was performed using the WGCNA package (v1.72) in R. The gene expression matrix from six samples [three control (CK) and three drought-treated (D)] was preprocessed by filtering genes with low expression variance (standard deviation ≤ 0.8), resulting in high-variance genes retained for constructing the network. To meet the assumption of a scale-free topology, the soft-thresholding power was determined by analyzing the scale-free topology fit index (R²) and mean connectivity for powers ranging from 1 to 30. A power of 30 was selected (R² > 0.9) to construct the weighted adjacency matrix, which was subsequently transformed into a topological overlap matrix (TOM). A hierarchical clustering algorithm using dynamic tree cutting was applied to the TOM-based dissimilarity matrix to identify co-expression modules, ultimately classifying all genes into 40 distinct modules.

To identify modules significantly associated with specific physiological traits, Pearson's correlation coefficients and their corresponding p-values were calculated between modules' eigengenes and the measured traits. Modules with an absolute correlation coefficient ≥ 0.3 and p < 0.05 were considered to be significantly trait-associated. For each significant module, gene significance (GS), defined as the correlation between the expression of each gene and the trait, and module membership (MM), defined as the correlation between the expression of each gene and the module's eigengenes, were calculated. Scatter plots of GS versus MM were used to validate the concordance between trait association and intramodular connectivity for genes within the module.

Leveraging the scale-free property of the WGCNA network, where most nodes (genes) have low connectivity but a few hub genes are highly connected, the top 50 genes with the highest intramodular connectivity in each module were defined as hub genes, which are presumed to play pivotal roles in maintaining the module's structure and function.

Data processing and analysis

Statistical data organization was performed using Microsoft WPS. Data analysis was conducted in SPSS 23.0, with independent sample t-tests applied to compare significant differences between the control and drought-stress treatment groups. The results were visualized as bar charts.

As shown in Fig. 1a, drought-stressed plants exhibited severe morphological alterations compared with the controls, including growth inhibition and leaf yellowing, curling, and wilting. Consistent with these phenotypic changes, the REC was significantly increased under drought stress (Fig. 1b). These results demonstrate that drought stress profoundly affects the growth and physiology of F. rubra.

Figure 1. 

(a) Morphological changes in F. rubra under control (CK) and drought stress (D). (b) Relative electrical conductivity (REC). Data are presented as mean ± standard deviation (SD). Statistical significance between groups was determined by analysis of variance (ANOVA). *** p < 0.001.

Determination of photosynthetic and fluorescence parameters

As shown in Fig. 2, under drought stress, the fluorescence parameters F_v/F_m, F_q'/F_m', and ETR (Fig. 2a, c, d) of F. rubra were significantly reduced (*** p < 0.001) compared with the control. F_q'/F_m' decreased by a highly significant 58.7% (from 0.46 to 0.19). However, F₀ and NPQ (Fig. 2b, e) showed no significant changes (p > 0.05). In contrast, the fluorescence parameter 1–qP (Fig. 2f) rose markedly under drought stress, increasing from 0.2 to 0.42, which is a significant increase of 110% (p < 0.01). Among the photosynthetic parameters, the photosynthetic rate, stomatal conductance, intercellular CO₂ concentration (Ci), and transpiration rate (Fig. 2g−j) of F. rubra under drought stress were also severely inhibited by drought (*** p < 0.001). Specifically, the photosynthetic rate dropped sharply from 18.44 to 1.59 μmol·m⁻²·s⁻¹, a highly significant decrease of 91.37% (Fig. 2g). The transpiration rate also decreased from 4.34 to 1.37 mmol H₂O m⁻²·s⁻¹, a highly significant reduction of 68.43% (Fig. 2j). These results indicate that drought stress severely suppresses the physiological processes related to light energy utilization and photosynthetic capacity in F. rubra.

Figure 2. 

Photosynthetic changes in F. rubra under drought stress (D) compared with the control (CK): (a) F_v/F_m, (b) F₀, (c) F_q'/F_m', (d) ETR, (e) NPQ, (f) 1–qP, (g) photosynthetic rate, (h) stomatal conductance, (i) intercellular CO₂ concentration (Ci), and (j) transpiration rate. All data are presented as the mean ± SD. Two groups were analyzed by ANOVA. ** 0.001 < p < 0.01; *** p < 0.001.

Osmotic and oxidative balance

As shown in Fig. 3, drought stress significantly induced the accumulation of osmoregulatory substances and oxidative markers in F. rubra. The contents of MDA, soluble protein, soluble sugar, and proline (Fig. 3a–d) increased substantially (*** p < 0.001) in comparison with the control. MDA content exhibited an increase from 9.29 µmol·g⁻¹ in the control to 24.27 µmol·g⁻¹ under drought stress (Fig. 3a), representing a 161.25% rise. As shown in Fig. 3d, there was a substantial increase in proline content from 116.12 to 768.30 µmol·g⁻¹, representing a 561.64% rise.

Figure 3. 

Osmoregulatory and antioxidant changes in F. rubra under drought stress (D) compared with the control (CK): (a) Malondialdehyde (MDA); (b) soluble protein; (c) peroxidase (POD); (d) catalase (CAT); (e) soluble sugar content; (f) proline; (g) ascorbate peroxidase (APX); and (h) superoxide dismutase (SOD). All data are presented as the mean ± SD. Differences between the two groups were analyzed by ANOVA. ** 0.001 < p < 0.01; *** p < 0.001.

The antioxidant enzyme activities of F. rubra were significantly enhanced under drought stress. POD, CAT, APX, and SOD activities exhibited an increase from 4.50 u·g⁻¹·min⁻¹, 22.30 u·g⁻¹·min⁻¹, 334.27 u·g⁻¹·min⁻¹, and 109.13 u·g⁻¹ in the control to 9.45 u·g⁻¹·min⁻¹, 189.09 u·g⁻¹·min⁻¹, 1,389.63 u·g⁻¹·min⁻¹, and 296.85 u·g⁻¹ under drought stress, respectively. POD activity exhibited a marked increase (Fig. 3e) (** p < 0.01), whereas CAT, APX, and SOD activities demonstrated a significant increase (Fig. 3f–h) (*** p < 0.001). These findings indicate that F. rubra possesses the capacity to enhance its antioxidant capacity by increasing the activity of antioxidant enzymes in response to drought stress.

Evaluation of transcriptome sequencing data

Six cDNA libraries were constructed via RNA-Seq analysis, with the raw sequencing data for each library detailed in Supplementary Table S1. All six samples exhibited Q30 percentages exceeding 95.81%, confirming high-precision measurements and superior data quality suitable for subsequent bioinformatic analysis.

DEGs screening

Differential expression analysis under drought stress identified 25,063 DEGs in F. rubra, comprising 14,773 upregulated genes and 10,290 downregulated genes (Fig. 4).

Figure 4. 

Volcano plot of DEGs in F. rubra under drought stress vs. the control group. Gray indicates non-DEGs, red indicates significantly upregulated DEGs, and blue indicates significantly downregulated DEGs.

GO functional enrichment analysis of DEGs

In order to explore the biological functions of DEGs in F. rubra under drought stress (Supplementary Table S2), GO enrichment analysis was performed (p < 0.05). In the BP category, the most significantly enriched terms were primarily associated with 'response to water', 'photosynthesis', and 'reductive pentose-phosphate cycle' (Fig. 5a). The expression of genes associated with the response to water were predominantly upregulated (Fig. 5b), whereas those involved in photosynthesis and the reductive pentose phosphate cycle were markedly downregulated (Fig. 5c). Within the CC category, the top four subcategories were 'chloroplast thylakoid membrane', 'Photosystem II', 'Photosystem I', and 'thylakoid' (Fig. 5a). It is evident from Fig. 5c that all of these were downregulated. In the MF category, the primary subcategories encompassed 'chlorophyll II binding', 'ribulose-bisphosphate carboxylase activity', 'peroxidase activity', and 'electron transfer activity' (Fig. 5a). Among these, the first three activities were downregulated (Fig. 5c), whereas 'electron transfer activity' was upregulated (Fig. 5b).

Figure 5. 

Top 30 bar charts of GO functional enrichment analysis of DEGs in F. rubra under drought stress compared with the control: (a) Total, (b) upregulated, and (c) downregulated. Screened GO terms with PopHits ≥ 5 in the biological process (BP), cellular component (CC), and molecular function (MF) categories, sorted by –log₁₀ p-value.

Transcription factor analysis

Transcriptome analysis identified 630 differentially expressed transcription factors (TFs), distributed across 64 families, in F. rubra under drought treatment (Fig. 6). The top five TFs families are AP2/ERF-ERF (58), NAC (51), bHLH (44), bZIP (43), and C2H2 (41). Notably, the AP2/ERF-ERF family has been found to contain the highest number of both up- and downregulated genes, suggesting that it plays a key role in the response to drought stress. Genes belonging to the MADS-M-type, Tify, DBP, and E2F-DP families have been observed to be upregulated, implying their potential function in promoting drought tolerance. Conversely, genes belonging to the AP2/ERF-AP2, HB-BELL, C2C2-YABBY, and BES1 families were found to be downregulated, suggesting that they may be subject to suppression under drought stress.

Figure 6. 

Differentially expressed transcription factors of F. rubra in response to drought stress. Dark blue indicates all unigenes, yellow indicates differentially expressed unigenes, light blue indicates downregulated unigenes, and pink indicates upregulated unigenes.

KEGG enrichment analysis of DEGs

To explore the functions and metabolic pathways of DEGs in F. rubra under drought stress, KEGG enrichment analysis was conducted on these DEGs. The top 20 most significantly enriched pathways are presented in Fig. 7. Among these, 'plant hormone signal transduction' and 'mitogen-activated protein kinase (MAPK) signaling pathway–plant' were the most active physiological processes in leaves. The 'plant hormone signal transduction' pathway revealed 21 upregulated and 28 downregulated DEGs, whereas the 'MAPK signaling pathway–plant' pathway revealed 25 upregulated and 8 downregulated DEGs. 'Arginine and proline metabolism' was significantly enriched in upregulated pathways. 'Photosynthesis–antenna proteins' and 'Photosynthesis' were significantly enriched in downregulated pathways (Fig. 7). These findings indicate that drought has a significant inhibitory effect on photosynthesis and, consequently, on plants' energy conversion and growth.

Figure 7. 

Top 20 KEGG enrichment bubble charts of DEGs in F. rubra under drought stress: (a) Total, (b) upregulated, and (c) downregulated DEGs. Pathways with PopHits ≥ 5 are screened and sorted by –log₁₀p-value. The x-axis shows the enrichment score. Larger bubbles indicate more DEGs, and bubble colors range from blue through white and yellow to red, with redder colors indicating smaller p-values and greater significance.

DEGs analysis related to key response pathways to drought stress

Analysis of the key metabolic pathways revealed distinct expression patterns under drought stress (Fig. 8). In the photosynthesis pathway, the expression of most genes related to Photosystem I (PSI), PSII, Cytochrome b₆f complex (Cytb₆f) and Photosynthetic Electron Transport (PET) were downregulated (blue) under drought stress (Fig. 8a). Conversely, in the peroxisome pathway, the expression of genes such as HACL1, ACAA1, ECH1, IDH1, and CAT were found to be significantly upregulated (red) under drought stress, whereas FAR-related genes were predominantly downregulated (blue) (Fig. 8b). In the proline biosynthesis pathway, the expression of key genes, such as P5CS, rocF, rocD and pip, were significantly upregulated (red) under drought stress (Fig. 8c), indicating accelerated conversion of glutamate to L-glutamate-5-semialdehyde and enhanced metabolic pathways from ornithine and glutamate to proline.

Figure 8. 

Heatmap of DEGs in F. rubra under control (CK) and drought stress (D) conditions: (a) Photosynthesis; (b) peroxisome; and (c) proline biosynthesis. Rows represent genes, and columns represent treatments (CK-1, CK-2, CK-3, D-1, D-2, D-3), with colors indicating the relative gene expression levels (red for upregulation, blue for downregulation). The heatmap was created using gene_id and fragments per kilobase per million mapped fragments (FPKM) values (Supplementary Table S2).

The plant hormone signal transduction pathway plays a crucial role in mediating stress responses. A substantial number of DEGs related to plant hormone signal transduction were found in F. rubra under drought stress, including genes related to ABA, auxin, ethylene, jasmonic acid (JA), and salicylic acid, with the most pronounced changes in ABA and auxin. Within the ABA signaling pathway, 12 PP2C genes were found to be upregulated, and 6 of the 7 SnRK2 genes were upregulated in F. rubra under drought stress, suggesting the accumulation of ABA (Fig. 9a). In the auxin pathway, key genes associated with auxin biosynthesis, including TIR1, IAA, ARF, GH3, and SAUR, were differentially expressed. Among these, two TIR1 genes were found to be upregulated, and there were both upregulated and downregulated genes in IAA, GH3, and SAUR (Fig. 9b), indicating complex auxin-mediated regulation of growth and stress adaptation. The plant MAPK signaling pathway has been shown to be associated with stress tolerance and wounding responses. As illustrated in Fig. 9c, the analysis revealed 29 DEGs associated with ABA biosynthesis, along with PYL, PP2C, and SnRK2 genes that were also present in the plant hormone signal transduction pathway. This pathway encompasses MAPKKK17_18, MKK3, MPK6, and CAT genes. It has been highlighted that MAPKKK17_18, MKK3, and CAT were upregulated in ABA biosynthesis. Furthermore, the MAPK signalling pathway showed significant crosstalk with ethylene and JA biosynthesis, with 11 and 5 DEGs identified, respectively (Fig. 9d). Critically, MPK6 was enriched in ABA, ethylene, and JA biosynthesis, whereas MKK3 was enriched in ABA and JA biosynthesis. This finding indicates that MPK6 and MKK3 play a regulatory role in the biosynthesis of multiple hormones, thereby enhancing plants' drought resistance.

Figure 9. 

Heatmap of DEGs related to (a, b) plant hormone signal transduction and (c, d) the MAPK signaling pathway in F. rubra under control (CK) and drought stress (D) conditions. Red boxes indicate upregulated genes, blue boxes indicate downregulated genes, and green boxes indicate genes with both upregulation and downregulation. Rows represent different genes, and columns represent different treatments (CK-1, CK-2, CK-3, D-1, D-2, D-3). Colors show relative gene expression levels (red for upregulation, blue for downregulation). The heatmap was created using gene_id and FPKM values (Supplementary Table S2).

Weighted gene co-expression network analysis

A weighted gene co-expression network was successfully constructed using 25,305 genes that passed the variance filter. By setting a soft-thresholding power of 30 to achieve a scale-free topology, all genes were grouped into 40 co-expression modules (Supplementary Fig. S1). Subsequent module–trait association analysis identified two modules, ivory and blue2 (Supplementary Fig. S1), which showed the most significant correlations with drought-responsive traits (e.g., proline content and photosynthetic rate), with correlation coefficients exceeding |0.8| (p < 0.01). Hub genes, defined as the most highly connected nodes within each module, were then identified. The top five hub genes in the ivory module were TRINITY_DN1227_c0_g1_i6_2, TRINITY_DN36541_c0_g1_i6_2, TRINITY_DN46909_c0_g1_i2_1, TRINITY_DN73423_c0_g1_i1_2, and TRINITY_DN4884_c0_g1_i22_2, and the top five in the blue2 module were TRINITY_DN1282_c0_g1_i14_3, TRINITY_DN13667_c0_g1_i1_3, TRINITY_DN1193_c0_g1_i14_1, TRINITY_DN7229_c0_g1_i7_1, and TRINITY_DN16159_c0_g1_i6_1, suggesting their potential pivotal roles in the drought response network.

Drought stress, as a critical abiotic factor, severely restricts plant growth and the ecosystem's stability, particularly in arid and semi-arid regions where water scarcity directly threatens agricultural sustainability and turfgrass quality. To elucidate the physiological and molecular mechanisms underlying drought tolerance in F. rubra, we subjected plants to a sustained 21 d drought treatment and systematically evaluated their responses, including measurements of REC, osmotic adjustment substances (proline and soluble sugars), photosynthetic efficiency (F_v/F_m, ETR), antioxidant enzyme activities, and transcriptome profiles.

Physiological mechanism of F. rubra under drought stress

Drought stress markedly inhibited the growth of F. rubra, as evidenced by leaf yellowing and curling (Fig. 1a), with the latter being a known drought avoidance strategy through the regulation of leaf water potential^[41].The significant increase in REC (Fig. 1b) suggests cell membrane damage under stress^[42]. At the photosynthetic level, drought stress reduced the fluorescence parameters (F_v/F_m, F_q'/F_m', and ETR) (Fig. 2a, c, e) but increased 1–qP (Fig. 2f), reflecting enhanced nonphotochemical quenching as a protective response to excess light^[43], though at the cost of reduced light-use efficiency. Additionally, the photosynthetic rate, stomatal conductance, and transpiration rate were all substantially reduced under drought (Fig. 2g–j), likely as a result of CO₂ limitation through efficient stomatal closure^[44]. These findings suggest that F. rubra adopts a typical water-conserving strategy through the suppression of photosynthesis and gas exchange, a pattern consistent with reports in species such as Festuca arundinacea and Lolium perenne^[45]. Notably, despite the observed decline in photosynthetic rate and the morphological appearance of leaf yellowing, the chlorophyll content in plants subjected to drought stress actually exhibited an increasing trend (Supplementary Fig. S2). This finding appears to contrast with the common conclusion that drought leads to chlorophyll's degradation^[46], but it may point to a potential unique stress strategy of F. rubra.

Under drought stress, plants partly maintain membrane stability through osmotic regulation. ROS disrupt redox homeostasis, inducing senescence and apoptosis and inhibiting growth^[47]. Excess ROS leads to lipid peroxidation and tissue damage, with MDA as a key marker^[48]. In this study, the MDA content in F. rubra increased by 161.25% under drought (Fig. 3a), indicating significant oxidative damage to membrane lipids. Conversely, osmolytes such as proline, soluble sugars, and soluble proteins contribute to cellular osmotic adjustment and stress protection. Their accumulation under drought conditions (Fig. 3b–d) aligns with previous findings linking their levels to stress resistance^[49]. Notably, the marked rise in proline underscores its central role in the drought response of F. rubra. Concurrently, the activities of antioxidant enzymes, including SOD, CAT, POD, and APX, which mitigate ROS-induced damage by catalyzing ROS conversion and preventing ·OH formation^[50], were significantly increased under drought (Fig. 3e–h). This coordinated upregulation demonstrates that F. rubra activates its enzymatic antioxidant defense system as a key adaptation to oxidative stress^[50]. Interestingly, differing results in other species, such as decreased CAT activity in Sorghum bicolor under drought^[51], may reflect the influence of species-specific responses, stress intensity, or experimental settings. Further investigation into phase-specific stress-regulating mechanisms under field conditions is therefore warranted.

Transcriptomic analysis of F. rubra under drought stress

Transcriptome analysis identified 25,063 DEGs in F. rubra under drought stress, including 14,773 upregulated and 10,290 downregulated genes, thus demonstrating extensive transcriptional reprogramming. Among these, TFs are key regulators of stress responses^[52]. In this study, the most enriched TF families were AP2/ERF-ERF (58), NAC (51), and bHLH (44) in response to drought stress in F. rubra (Fig. 6). The AP2/ERF family, which is unique to plants, plays a central role in ethylene-responsive transcription and adaptation to abiotic stress^[53]. Our results revealed the extensive activation of AP2/ERF-ERF family members, consistent with the enrichment of hormone signal transduction pathways. In rice (Oryza sativa), the NAC family, known to regulate development and stress responses, has been shown to enhance drought tolerance through root system remodeling and activation of downstream stress-related genes^[54]. Similarly, bHLH TFs are involved in ABA signaling and stomatal regulation. For instance, overexpression of AhbHLH112 in peanut (Arachis hypogaea) increased ABA biosynthesis and enhanced drought tolerance^[55]. The prominent enrichment of these three TF families implies a synergistic role in drought adaptation in F. rubra, possibly through stomatal regulation, osmolyte accumulation, and hormone signaling. These TFs thus represent potential targets for further functional validation in studies of drought resistance in turfgrass.

GO and KEGG enrichment analyses revealed that drought stress significantly affected photosynthesis, proline metabolism, and peroxisome-related pathways in F. rubra (Figs 5 and 7). Consistent with the decline in physiological photosynthetic parameters (Fig. 2), most genes related to PSI, PSII, Cytb₆f, and PET complexes were markedly downregulated (Fig. 8a), indicating transcriptional suppression of photosynthetic capacity. Conversely, the upregulation of peroxisome-related genes (Fig. 8b) suggests enhanced peroxisome activity, which may contribute to ROS scavenging and cellular protection. Additionally, proline metabolism-related genes such as P5CS, rocF, rocD, and pip were significantly upregulated (Fig. 8c). In particular, upregulation of P5CS promoted the conversion of glutamate and proline synthesis, corroborating the physiological accumulation of proline (Fig. 3f). This suggests that F. rubra enhances proline biosynthesis at the transcriptional level to improve osmotic adjustment under drought, confirming proline metabolism as a central drought-response pathway^[56]. Overall, F. rubra copes with drought by downregulating photosynthetic genes while upregulating genes related to proline metabolism and peroxisomal metabolism, forming a coordinated molecular response to alleviate drought-induced damage.

KEGG enrichment revealed that 'plant hormone signal transduction' and 'MAPK signaling pathway – plant' were among the most significantly activated pathways in F. rubra leaves under drought stress. ABA signaling is a key regulator of stomatal closure and water retention^[57]. In this study, the expression levels of 12 PP2C genes and 6 SnRK2 genes were increased (Fig. 9a). These findings suggest activation of the core ABA signaling module (PYL–PP2C–SnRK2), consistent with known ABA-mediated stomatal regulation^[58]. This transcriptional upregulation of PP2C and SnRK2 genes is consistent with the established Arabidopsis thaliana model, where the ABA binding to PYL receptors relieves the repression of SnRK2 kinases by PP2C phosphatases, leading to the activation of stress-responsive gene expression^[59]. This supports the conclusion that F. rubra enhances drought tolerance by modulating stomatal aperture and water conservation via ABA signaling^[60]. Auxin signaling, another key drought-responsive pathway, showed altered expression of TIR1, IAA, ARF, GH3, and SAUR family genes (Fig. 9b). Although TIR1 was upregulated, downstream IAA showed mixed expression patterns, possibly caused by interference with ubiquitin-mediated regulation^[61]. The downregulation of ARF may affect auxin–DNA signal transduction, whereas the differential expression of SAUR and GH3 points to their roles in calcium/calmodulin signaling and auxin homeostasis^[62]. These findings indicate that hormone signal pathways, particularly ABA and auxin signaling, play crucial roles in mediating the transcriptional adaptive response of F. rubra.

ABA regulates stomatal conductance and induces drought-responsive genes, serving as a key signal in plants' drought response^[63]. In Arabidopsis, a complete ABA-activated MAPK cascade has been identified, comprising MAP3K17/18, MKK3, and C-group MAPKs (MPK1/2/7/14), which is transcriptionally upregulated by ABA and is essential for regulating abiotic stress responses^[64]. One of its downstream targets is the MAPK cascade (MAPKKK–MAPKK–MAPK), which regulates the expression of antioxidant enzymes and alleviates ROS damage^[65]. Although our data suggest a conserved role for the MAPKKK17_18-MKK3 module in ABA signaling between Arabidopsis^[64] and F. rubra, the specific downstream MAPK diverges (MPK1/2/7/14 in Arabidopsis^[64] versus MPK6 in F. rubra), indicating potential species-specific rewiring of the pathway's output. In F. rubra, transcriptome analysis identified the differential expression of MAPKKK17_18, MKK3, and MPK6 under drought stress (Fig. 9c). Upregulation of MAPKKK17 may trigger downstream MKK3, which, in turn, activates MPK6, a kinase implicated in abiotic stress responses and programmed cell death^[66]. Interestingly, MKK3 and MPK6 were also upregulated in the ABA and JA signaling pathways (Fig. 9c,d), and MPK6 was notably expressed in the ABA, ethylene, and JA pathways. This suggests that MPK6 functions as a signal integration node, coordinating multiple hormone signals to regulate drought resistance^[67]. These hormones are involved in stomatal closure (ABA), antioxidant activation (JA), and defence responses (ethylene). The coordinated expression of MKK3 and MPK6 highlights the crosstalk between MAPK and hormone pathways in F. rubra's stress response. Their roles in fine-tuning drought adaptation underscore the complexity of MAPK signaling. Future studies exploring their tissue-specific and stage-specific expression could further clarify their functional significance.

This study comprehensively elucidated the drought resistance mechanisms of F. rubra by integrating physiological and transcriptomic analyses. Under drought stress, plants exhibited photosynthetic inhibition, enhanced antioxidant capacity, and the accumulation of osmoregulatory substances. Transcriptome profiling revealed that these responses are orchestrated by the central roles of ABA signaling, MAPK cascades, and transcription factors such as AP2/ERF, NAC, and bHLH. The coordinated regulation of proline metabolism and peroxisome function also contributed to stress mitigation. These findings deepen the theoretical understanding of stress biology in turfgrass and provide a valuable foundation for breeding drought-tolerant F. rubra varieties. The identified key TFs and signaling pathways will provide candidate resources for future functional validation and genetic engineering of drought-tolerant varieties.