Abstract Background Sustainable production of high-quality feedstock has been of great interest in bioenergy research. Despite the economic importance, high temperatures and water deficit are limiting factors for the successful cultivation of switchgrass in semi-arid areas. There are limited reports on the molecular basis of combined abiotic stress tolerance in switchgrass, particularly the combination of drought and heat stress. We used transcriptomic approaches to elucidate the changes in the response of switchgrass to drought and high temperature simultaneously. Results We conducted solely drought treatment in switchgrass plant Alamo AP13 by withholding water after 45 days of growing. For the combination of drought and heat effect, heat treatment (35 °C/25 °C day/night) was imposed after 72 h of the initiation of drought. Samples were collected at 0 h, 72 h, 96 h, 120 h, 144 h, and 168 h after treatment imposition, total RNA was extracted, and RNA-Seq conducted. Out of a total of 32,190 genes, we identified 3912, as drought (DT) responsive genes, 2339 and 4635 as, heat (HT) and drought and heat (DTHT) responsive genes, respectively. There were 209, 106, and 220 transcription factors (TFs) differentially expressed under DT, HT and DTHT respectively. Gene ontology annotation identified the metabolic process as the significant term enriched in DTHT genes. Other biological processes identified in DTHT responsive genes included: response to water, photosynthesis, oxidation-reduction processes, and response to stress. KEGG pathway enrichment analysis on DT and DTHT responsive genes revealed that TFs and genes controlling phenylpropanoid pathways were important for individual as well as combined stress response. For example, hydroxycinnamoyl-CoA shikimate/quinate hydroxycinnamoyl transferase (HCT) from the phenylpropanoid pathway was induced by single DT and combinations of DTHT stress. Conclusion Through RNA-Seq analysis, we have identified unique and overlapping genes in response to DT and combined DTHT stress in switchgrass. The combination of DT and HT stress may affect the photosynthetic machinery and phenylpropanoid pathway of switchgrass which negatively impacts lignin synthesis and biomass production of switchgrass. The biological function of genes identified particularly in response to DTHT stress could further be confirmed by techniques such as single point mutation or RNAi. Supplementary Information The online version contains supplementary material available at 10.1186/s12870-022-03477-0. Keywords: Panicum virgatum, Transcriptome, Drought stress, Heat stress, Transcription factors, Gene ontology, Differential gene expression Background Plants in the field are exposed to various environmental stresses which affect production and yield. These environmental stresses include abiotic factors such as DT, HT, and salinity and biotic stresses like pathogens, and insect pests, [[39]1]. Abiotic stresses are reported to reduce about 50% of crop production [[40]2]. Stress tolerance research has primarily focused on the response of plants to individual stress with limited information on plants’ adaptability to combined stresses such as HT and DT and salinity and DT [[41]3–[42]5]. Moreover, plants exhibit a unique expression pattern when exposed to multiple stresses [[43]6]. Hence to bridge the knowledge gap, we have compared the transcriptional response of switchgrass when exposed to individual DT stress or a combination of DT and HT stresses. The combined effect of DT and HT stresses has been shown to cause more damage to plants than when these stresses occur at separate times [[44]7, [45]8]. The mechanisms used by plants to adapt to multiple stresses can be complex. It has been shown that the effect of one stress could have a synergistic or antagonizing effect on other stress. DT, salinity, high and low temperature have been shown to promote the occurrence of pathogens and pests [[46]5]. In addition, the antagonizing effect of cold stress on osmotic stress during the induction of dehydration-responsive gene RD29A has been reported [[47]9]. Abscisic acid (ABA) was found to antagonize jasmonate-ethylene signaling pathways and mediates defense gene expression and disease resistance in Arabidopsis [[48]10]. Multiple stress in plants led to the expression of common overlapping genes due to a cross-talk of a signaling pathway. A previous study identified 22 genes that were induced commonly during DT, cold, and NaCl treatment [[49]11]. Some of the molecular mechanisms adopted by plants to combat stress include the release of Heat shock proteins or chaperons that are expressed during exposure to environmental cues [[50]12]. Transcriptome analysis of Arabidopsis showed that HT resistance is conferred by HT stress-responsive genes, plant hormones, and antioxidant enzymes [[51]13]. The importance of transcriptional gene regulation in plants under DT and HT stresses has been previously reported [[52]13]. RNA sequencing (RNA-Seq) has been commonly used to identify genome-wide transcript profiles in plants. Stress-responsive genes have been identified in tobacco and Arabidopsis when exposed to combined DT and HT stress by RNA-Seq technology [[53]14, [54]15]. Plant responses to single stress treatment of cold, high light, salt, HT, and flagellin have been compared to various combinations of these six pair of stresses (cold and high light, salt and HT, salt and high light, HT and high light, HT and flagellin respectively). The outcome of this study revealed how plants have evolved to withstand combination of these stresses [[55]4]. The combined effect of DT and HT stress has been studied in wheat [[56]16]. The effect of combined abiotic stress signaling such as DT, salinity, and metal in rice was found to be complex with the involvement of multiple genes, differential expression patterns in different developmental tissues, and protein-protein interaction [[57]17]. Furthermore, the separate impact of DT and HT and their combined effect on grain filling, physiological, vegetative, and yield traits were investigated in wheat [[58]8]. Switchgrass (Panicum virgatum L.) is a C4 warm-season perennial grass identified as a potential bioenergy crop [[59]18, [60]19]. It has been investigated for lignocellulosic ethanol production in the US, Canada, and Europe [[61]20] due to its high biomass yield. It serves as a potential alternative to nonrenewable fossil fuels, thereby providing energy security sources [[62]21]. Switchgrass requires a minimal amount of water and nutrients and can grow on marginal croplands [[63]22]. Its rapid growth rate and broad adaptability contribute to a stable and high biomass supply. Switchgrass positively impacts the soil by improving soil quality, preventing erosion, and reducing soil nutrients [[64]23]. Switchgrass, like many other plants, is generally faced with extreme biotic and abiotic stresses. These stresses can be detrimental by causing retardation in plant growth, development, and even death [[65]24]. DT is a significant abiotic stress that limits switchgrass use as biofuel production. There is evidence of DT as an essential economic risk factor affecting biofuel production [[66]25]. Molecular mechanisms underlying DT responses in plants have been addressed in various articles [[67]26]. A previous report suggests DT could considerably reduce the yield and quality of biomass for biofuel production [[68]27]. The effect and response of switchgrass germplasms to DT stress have been evaluated in previous studies [[69]28–[70]30]. High temperatures in the Southern United States are projected to reduce switchgrass biomass in 2080–2090 [[71]22]. Similarly, various studies have reported the impact of high temperatures on switchgrass, emphasizing physiology, cell wall composition, biomass, and yield. A significant decrease in biomass yield was observed across various switchgrass genotypes due to the impact of high temperatures [[72]22, [73]31]. There is increasing research in switchgrass, and among the area of research is gene regulation. Transcriptome analysis has been used to determine genes associated with biomass production in switchgrass [[74]29]. The characterization of DT and HT responsive microRNAs has been recently reported [[75]18]. Besides, the role of microRNAs during DT and salt stress in switchgrass has been reported [[76]32]. Although switchgrass is an essential bioenergy crop, less information on the biology of switchgrass is available when imposed with abiotic stresses [[77]23]. The molecular mechanisms of the tolerance of switchgrass to hot and dry climates is not well studied [[78]18]. Therefore, understanding the effect of stress combinations in switchgrass will be important to reveal genes associated with important traits such as biomass and biofuel production in response to multiple environmental stresses. Additionally, breeding DT and HT resistant switchgrass cultivars will be an important milestone. Although several studies have reported switchgrass response to a single DT or HT stress, there are no reports as far as we know on the combination of DT and HT abiotic stresses in switchgrass, especially with prolonged exposure to DT and HT stresses. To better understand plant responses to the full complement of environmental stresses, it is important to compare data on single stresses with data on multiple stresses. It is also important to identify the early transcriptional response to DT and HT stress versus the prolonged exposure of switchgrass to these stresses. This will provide an idea of signaling cross-talk in systems biology [[79]33]. In this study, we used RNA-Seq approache to characterize and quantify gene expressions in response to DT and combined effects of DT and HT stresses in switchgrass. Results RNA-Seq data quality and summary A total of 6965 million paired-end reads were obtained from RNA-Seq samples. The number of reads in each sample was 129 million on average. Around 85% of the reads can be aligned to the reference genome. About 63% of reads were aligned to genic regions. To assess the similarities and differences among these samples, we performed a hierarchical cluster analysis of the RNA-Seq data (Fig. [80]1). We found that non-treated samples were grouped together except the 72 h DT treated samples. In the group of stress treated samples, DTHT samples were grouped together except 144 h DTHT sample, which clustered with the group of DT samples. Fig. 1. [81]Fig. 1 [82]Open in a new tab Hierarchical clustering analysis of Control, DT, and DTHT treated samples Analysis of DT and DTHT responsive genes in switchgrass From the analysis, many genes were identified in response to the DTHT compared to only DT stress. In total, 3912 out of 32,190 genes were identified as DT and 4635 as DTHT responsive genes. Among those, 1615 genes were shared between the DT and DTHT data sets, when DT samples were compared to plants exposed to combined DTHT stress. These commonly expressed genes likely play critical roles in DT and HT tolerance in switchgrass. Further analysis showed that 1432 out of 2282 of the up-regulated responsive genes were unique (Fig. [83]2A) and 1604 out of 2345 down-regulated genes were unique to DTHT (Fig. [84]2B). Similarly, for DT samples, 1307 out of 2157 up-regulated responsive genes were unique, while 1013 out of 1754 down-regulated genes were unique (Fig. [85]2A and B). Fig. 2. [86]Fig. 2 [87]Open in a new tab The number of common and specific up-regulated (A), and down-regulated (B) genes among switchgrass during DT and DTHT stress in the Venn diagram. The genes were significantly differentially expressed (DE) in more than one comparison of the time point, 0 h, 72 h, 96 h, 120 h, 144 h, and 168 h. DE genes for each comparison were quantified at log2 fold changes and P-value < 0.05 In our data, Pavir.6 [88]KG130600.v4.1 provided the best hit to Arabidopsis AT1G22360.1 (UDP-glucosyl transferase 85A2 (UGT85A2) and it is the only DT-responsive gene that showed both up and down-regulation between the time points after imposing DT treatment (Additional file [89]5, DT). This gene was significantly down-regulated at time points DT 96 h and DT 120 h after which its expression markedly up-regulated at 168 h. Through GO enrichment analysis (Fig. [90]3a, b, Additional file [91]6), we found that there were significantly enriched terms in all biological process, molecular function, and cellular component functional categories. In the biological process category, the enriched GO terms included photosynthesis, single-organism metabolic process, and metabolic process. GO enrichment analysis show that the GO terms; “response to stress” and “response to water”, with p-values (0.00042 and 0.00054, respectively) were smaller than 0.05 although the FDR values were above 0.05 (0.083 and 0.093, respectively). Eight out of 15,902 genes belonged to the GO term of response to water in the switchgrass genome whereas seven out of 3912 DT responsive genes also belonged to the GO term of response to water. In molecular function, some of the enriched terms were oxidoreductase activity, catalytic activity, and cofactor binding. In the cellular component category, the enriched terms were photosystem, photosynthetic membrane, and thylakoid part. We further performed KEGG enrichment analysis on the DT responsive genes. We found that these DEGs were enriched in the following KEGG pathways (Additional file [92]7): protein phosphatase 2C, glutaredoxin 3, homeobox−leucine zipper protein, jasmonate ZIM domain−containing protein, and solute carrier family, xyloglucan: xyloglucosyl transferase, HSP20 family protein, adenylate kinase, and UDP-glucuronate decarboxylase. Fig. 3. [93]Fig. 3 [94]Open in a new tab a. The Gene Ontology (GO) terms enriched by responsive genes to DT stress. The DEGs were annotated against the GO database. The GO terms are in the three GO domains (biological process, molecular function and cellular compartment). These terms were significantly enriched (p < 0.05) in combined DT and HT treated samples compared to control plants. The number of genes enriched in each term were plotted against the GO term. b. The Gene Ontology (GO) terms enriched by responsive genes to DTHT stress. The DEGs were annotated against the GO database. The GO terms are in the three GO domains ( biological process, molecular function, and cellular compartment). These terms were significantly enriched (p < 0.05) in combined DT and HT treated samples compared to control plants. The number of genes enriched in each term were plotted against the GO term Pavir.9NG755000.v4.1 which provided the best hit to (ATHCHIB, B-CHI, CHI-B, HCHIB, PR-3, PR3) is a basic chitinase gene was significantly down-regulated at 144/72 h and subsequently up-regulated after prolonged DT and HT stress at 168/96 h. Similarly, genes such as Pavir.5KG627200.v4.1, Pavir.2NG348700.v4.1 and Pavir.2NG348700.v4.1 with best hit to Arabidopsis genes encoding delta 1-pyrroline-5-carboxylate synthase 2 (AT3G55610.1), cytochrome P450, family 76, subfamily C (AT2G45550.1), polypeptide 4, and DUF1012 (AT5G43745.1) respectively were significantly down-regulated at 144/72 h (Additional file [95]5, DTHT). These genes at 168/96 h were significantly up-regulated after prolonged DT and HT stress, suggesting the possible role of these genes in protecting the plant during extreme environmental conditions. To study the functions of these responsive genes, GO enrichment analysis was performed. The main GO term from the enrichment analysis was the GO term (GO:0008152; metabolic process) which showed significant enrichment (FDR; 0.0014) (Fig. [96]3). None of the GO terms shows significant enrichment in combined DT and HT stress responsive genes, indicating that DTHT transcriptomic changes were not predictable from single stress treatments. In the category of biological process, there were 10 most enriched GO terms with P-value <= 0.05. These 10 GO terms were response to water, single-organism metabolic process, single-organism biosynthetic process, response to abiotic stimulus, organonitrogen compound metabolic process, photosynthesis, oxidation-reduction process, response to stress, nitrogen compound transport, and transmembrane transport respectively. We further performed KEGG enrichment analysis on the DTHT responsive genes (4635 genes). We found that these responsive genes were enriched in the following KEGG pathways (Additional File [97]7): adenylate kinase and protein phosphatase 2C. HT responsive genes in switchgrass The HT stress genes were deduced from the DEGs of DT and DTHT. In total, 2338 out of 32,190 genes were identified as HT responsive genes (Additional file [98]5). There were 1064 up-regulated genes and 1274 down-regulated genes. The functions of these responsive genes and GO annotation were presented (Additional file [99]6). In the category of biological process, these genes showed enrichment in the GO terms of organic cyclic compound catabolic process, organonitrogen compound catabolic process and heterocycle catabolic process, etc. In the category of molecular function, these genes showed enrichment in the GO terms of organic cyclic compound catabolic process, organonitrogen compound catabolic process and heterocycle catabolic process, etc. In the categories of cellular components, these genes showed enrichment in the GO terms of photosystem II oxygen-evolving complex, photosystem II, and thylakoid membrane. We also performed KEGG enrichment analysis on the HT specific responsive genes. We found that these responsive genes were enriched in the jasmonate ZIM domain-containing protein pathway (Additional File [100]7). Transcription Factors (TF) for DT, DTHT and HT responses The TFs identified from the analysis are shown in Table [101]1, and Additional file [102]8. These DT and DTHT responsive TFs belong to 45 different TF families. Out of 91,838 proteins on the switchgrass genome, 3948 were identified as transcription factors (TFs). A total of 1383 TFs were identified out of 32,190 genes that were used for identifying stress responsive genes. There were 209 genes identified as TFs out of 3912 DT responsive genes. Heat maps were generated to show expression patterns of these 209 genes in all the samples (Fig. [103]4A). Similarly, there were 220 genes identified as TFs out of 4635 DTHT responsive genes. A heat map was generated to show expression patterns of these 220 genes in all the samples (Fig. [104]4). A total of 106 genes out of the 2339 predicted HT responsive genes, were identified as TFs. Heat map was generated to show expression patterns of these 106 genes in all the samples (Fig. [105]4C). Table 1. Different families of TFs responsive to solely DT and combined DTHT stresses Transcription factor type DTvsCtrl DTHTvsCtrl DTHTvsDT bHLH 22 20 10 NAC 16 15 13 ERF 19 19 6 bZIP 17 17 5 MYB_related 10 17 10 MYB 12 15 7 WRKY 14 11 6 HD-ZIP 15 7 3 C3H 7 13 1 [106]Open in a new tab Fig. 4. [107]Fig. 4 [108]Open in a new tab Heat map with clusters based on FPKM values for A) DT vs Control, B) DTHT vs control and C) DTHT vs DT TFs. The Heat map shows a grouping of control samples and stress samples. Extended periods of DTHT to stress samples showed abundant up-regulated TFs (A and B) and down-regulated TFs (C) compared to their control samples. For example, there were more responsive TFs which were up-regulated at time 144/72 h compared to its control sample at Control 144/72 h (A) Pathway analysis of DT and HT responsive genes An overview of the secondary metabolism pathway is displayed in Fig. [109]5(A and B). We found a large number of plant secondary metabolites such as flavonoids, terpenes, and phenylpropanoids were down-regulated in DTHT vs control samples compared to DT vs control samples. Fig. 5. [110]Fig. 5 [111]Open in a new tab Metabolism overview in MapMan showing the DEGs between DT vs Control (A) and DTHT vs control (B) switchgrass samples. The log-fold ratio is indicated as a gradient with red color (down-regulated) and blue color (up-regulation) Co-expression network We performed weighted gene co-expression network analysis to identify genes involved in response to the DT and DTHT stresses. Most of co-expressed genes usually participate in the same biological processes [[112]34–[113]36]. In our co-expression analysis we identified 68 modules with distinct expression patterns (Additional file [114]11). To study whether the DEGs were enriched in some of the modules, Fisher’s exact test and multiple test correction (Benjamini-Hochberg method) were performed [[115]4] . Of the modules that have more than 100 genes, DT responsive genes were enriched in module 5, 7, 14, 17 and 25. DTHT responsive genes were enriched in module 1, 2, 3, 7 and 17. HT responsive genes were enriched in module 1, 2, 8, 9, 15, 16, 17 and 25. GO enrichment analysis of the genes of these modules were performed using agriGO. Results for GO enrichment are provided in (Additional file [116]12). Heat maps were generated for these 12 unique modules (Additional file [117]2). In module 7 and module 17, both DT responsive genes and DTHT responsive genes were enriched. In module 7 and module 17, genes were up-regulated after stress treatment. In module 7, the genes were enriched in GO terms of response to water, response to acid chemical, lipid biosynthetic process, and response to the oxygen-containing compound, or biological process. In module 17, the genes were enriched in GO terms of regulation of nucleic acid-templated transcription, regulation of RNA biosynthetic process, regulation of RNA metabolic process and regulation of transcription, DNA-templated, etc. for biological process. In module 1 (Fig. [118]6), most genes were up-regulated during the initial HT treatment at DTHT 96/24 h. Down-regulation of most of the genes in the same module occur and then up-regulated again at an extensive HT at DTHT168/96 h. Similarly, in module 8 which is enriched with HT responsive genes, showed upregulation of genes at the initial stage of imposing HT at DTHT96/24 h. In module 1, the genes were enriched in GO term biological processes such as translation, peptide biosynthetic process, amide biosynthetic process and peptide metabolic process. In module 8, the genes are enriched in GO terms including; multi-organism reproductive process, multi-multicellular organism process, cell recognition, and pollination for biological processes. Also, DTHT responsive genes were enriched in module 9 with most of the responsive genes recorded at time point DTHT96/24 h and DTHT120/48 h. A number of the genes recorded at DTHT96/24 h and DTHT120/48 h were enriched in different class of metabolic processes. Fig. 6. [119]Fig. 6 [120]Open in a new tab Heat map indicating genes enriched in module 1 from the WGCNA analysis. DTHT and HT responsive genes were enriched in module 1 DT and DTHT responsive genes in DroughtDB There were 386 genes from the switchgrass genome that have the best hits to Arabidopsis genes in the droughtDB [[121]37] Of these 386 genes, 172 were found in the 32,190 genes in this study. Detailed gene expression patterns of these 172 genes were shown in the heat map (Additional file [122]3). Out of these 172 genes, there were 35 DT responsive genes and 27 DTHT responsive genes in which 12 were common (Additional file [123]13). A list of the DT and DTHT genes have been indicated in Tables [124]2 and [125]3, respectively. The gene IDs, biological functions, the phenotype of mutants, references, tags of the