Analysis of drought and heat stress response genes in rice using co-expression network and differentially expressed gene analyses

Data collection

Multiple gene expression profiling datasets, including high-throughput sequencing (Illumina HiSeq 2000/Illumina HiSeq 4000/Illumina NovaSeq 6000) and array datasets (Affymetrix Rice Genome Array Platforms), were sought and retrieved from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). These high-highthroughput sequencing datasets included GSE221542, GSE168650 (Kan et al., 2022), and GSE159816 (Zu et al., 2021). These array datasets included GSE136746 (Ps et al., 2017), GSE41648 (Sharma et al., 2021), GSE14275 (Hu, Hu & Han, 2009), GSE93917 (Wang et al., 2020), and GSE83378 (Wei et al., 2017) (Table 1). Gene symbols for these GEO datasets were annotated using the National Center for Biotechnology Information (NCBI), Rice Annotation Project database (RAP-db) (https://rapdb.dna.affrc.go.jp/), and the Rice Genome Annotation Project (http://rice.uga.edu/index.shtml). The data were processed using R (version 4.2.3) and RStudio (version 2023.03.0) software. GSE221542 contains 15 samples, including three water levels and two heat levels, each with three replicates.

Weighted gene co-expression network analysis of drought/heat response genes

Counts per million were computed to standardize the sequencing depth of RNA-seq data using the R package “edgeR” (Robinson, McCarthy & Smyth, 2010). Using the weighted gene co-expression network analysis (WGCNA) (Langfelder & Horvath, 2008) package in R (version 4.2.3), aco-expression network was constructed using the following steps. First, the average expression of each gene under different levels of drought or heat stress was calculated, and genes that did not exhibit any changes in expression were filtered out. Second, normalization of gene expression levels to a range of 0–1 was followed by the calculation of Pearson’s correlation coefficients, which is used to measure the similarity of co-expression between genes. Third, to ensure a scale-free network distribution, an appropriate beta value was selected for the adjacency matrix weights to construct a topological overlap matrix for module clustering and segmentation. Finally, to select modules related to drought or heat responses, the relationship between each network module and the sample phenotype was analyzed.

Gene Ontology (GO) terms were used to enrich selected genes (Tian et al., 2017). The analysis results were presented using the R package “clusterProfiler” for visualization (Yu et al., 2012). Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment (Kanehisa & Goto, 2000) analysis was also performed using the R package “clusterProfiler” (Yu et al., 2012). Using the CytoHubba (Chin et al., 2014) plugin of Cytoscape (3.9.1), based on the shortest paths, every gene of the key module was scored using the Maximal Clique Centrality (MCC) method, and the top 20 hub genes were selected.

Differentially expressed gene analysis with DESeq2 and GO enrichment in R

Differentially expressed gene (DEG) analysis was performed using the R package DESeq2 (Love, Huber & Anders, 2014). Raw count data from the RNA-seq experiments were imported into R, and genes with low expression were filtered using the “filterByExpr” function. Next, the “DESeqDataSetFromMatrix” function was used to create a DESeq2 object, which was then used to estimate size factors and dispersions using the “estimateSizeFactors” and “estimateDispersions” functions, respectively. A false discovery rate cutoff of 0.05 was applied to identify genes that were significantly differentially expressed, based on an absolute log2 fold change ≥ 1 and an adjusted p-value ≤ 0.05. All data analyses were performed using R software (version 4.2.3).

GO enrichment analysis was also performed to analyze DEGs using a previously described approach (Tian et al., 2017; Yu et al., 2012).

Intersection of hub genes and DEGs for candidate key genes

The top 20 hub genes from the filtered key modules were compared with the DEGs obtained from the filtering process. Based on their intersection, the candidate key genes along with their log2 fold change values were obtained. The Rice Gene Index (RGI) (https://riceome.hzau.edu.cn/) was used to determine the gene ID corresponding to the rice gene (Yu et al., 2023).

We searched for datasets on drought or heat treatments in the GEO database (Table 1). Count data were processed using the same method as above but not filtered for log2 fold change ≥ 1 and p-value ≤ 0.05. For array data, online GEO2R analysis was performed, and a matrix table containing the log2 fold change, p-value, and adjusted p-value data was downloaded.

Using the “pheatmap” package in R (version 4.2.3), the log2 fold change calculated from the different array or count data treatments was clustered and plotted. Key genes with high and stable expression levels were selected for further experiments.

Plant materials

The model rice variety Oryza sativa Nipponbare was subjected to appropriate environmental conditions, drought stress, and heat stimulation, as well as RNA extraction for quantitative real-time polymerase chain reaction (qRT-PCR).

Nipponbare rice seeds were germinated for 3 days at 30 °C. cultivated in Yoshida Rice Medium (Coolaber, Beijing, China) for 10 days (12 h light at 30 °C and 12 h dark at 27 °C every day). The control group was directly sampled. In the drought treatment group, the samples were grown for 10 days in the medium containing 2 g/L mannitol. In the heat treatment group, on the 10th day, rice plants grownwhich in normal medium werewas exposed to a temperature of 40 °C for 1 h and then sampled. The stress modeling method reference was based on the GSE221542 dataset.

RNA extraction

Whole shoot tissues (100 mg) from different treatment groups were weighed and placed in a grinding tube containing steel beads. The grinding tubes were immersed in liquid nitrogen for 10–20 min. Finally, the samples were freeze-ground at −20 °C for 120 s and returned to liquid nitrogen for storage. RNA extraction was performed using the FastPure Universal Plant Total RNA Isolation Kit (Vazyme, Nanjing, China), and the extracted total RNA was stored at −80 °C.

Quantitative real-time PCR

cDNA was synthesized using the Revert Aid First Strand cDNA Synthesis Kit (Thermo Scientific, Watham, MA, USA). qRT-PCR analysis was performed using a LightCycler 96 (Roche, Basel, Switzerland). eEF1 was used as the reference gene (Ambavaram & Pereira, 2011). Gene sequences were searched using Phytozome (https://phytozome-next.jgi.doe.gov/), and qRT-PCR primer sequences were designed using the primer blast tool of NCBI (https://blast.ncbi.nlm.nih.gov/Blast.cgi). The primers used in this study are listed in Table 2. See Supplementary 1 for the MIQE checklist (Bustin et al., 2009).

The relative expression level of target genes was calculated based on the 2^−ΔΔCt method for normalization (Livak & Schmittgen, 2001). The normalized qRT-PCR data were analyzed using the t-test to determine statistically significant differences in gene expression between the control and experimental groups (Wilson & Worcester, 1942). Statistical significance was set at p < 0.05.

Construction of co-expression network

The workflow followed in this study is depicted in Fig. 1.

WGCNA was applied to analyze the GSE221542 dataset, with a scale-free topology model fitting degree of 0.8 and a soft threshold of 30 selected for network construction (Figs. 2A–2B). A hierarchical clustering process was used to create a tree-like structure representing genes. Subsequently, gene modules were determined using the dynamic cutting method, followed by the calculation of the eigenvector value of each module. Similar modules were then merged to identify distinct modules, which were assigned different colors for better visualization (Fig. 2C).

Co-expression network module analysis

Six modules, namely, black (1,550 genes), green (1,646 genes), dark orange (3,658 genes), dark magenta (535 genes), royal blue (9,432 genes), and gray (234 genes), were obtained. The modules showed either positive or negative correlation with drought or heat stress, and the genes within these modules were either upregulated or downregulated, suggesting that the genes respond differently under different stress conditions. The green module with heat and the dark magenta module with drought had the highest positive correlation coefficients (0.98 and 0.71, respectively) (Fig. 3A). According to the scatter plots, the genes in the green module were highly correlated with heat stimulation, whereas the genes in the dark magenta module showed a weak association with drought stress (Figs. 3B–3E). The other modules showed low correlation with heat or drought stress (Fig. S1).

Cytoscape software was used to process the dark magenta and green modules separately and visualize the co-expression network obtained from WGCNA. Genes in the network were scored using the maximal clique centrality (MCC) method, and the top 20 hub genes with the highest correlations with other genes were selected (Figs. 3F–3G). These genes were located at the most central positions in the co-expression network, and they may play a central regulatory role in drought and heat stress.

In addition, KEGG enrichment analysis showed that genes in the green module are involved in processes such as protein synthesis in the endoplasmic reticulum and RNA splicing, whereas genes in the dark magenta module are involved in essential processes such as carbon metabolism, synthesis of amino acids and coenzyme factors, and glycerolipid metabolism. In addition, both gene modules are involved in carbon fixation in photosynthetic organisms (Fig. 4A). The results of the GO enrichment analysis showed that genes in the green module were related to biological processes such as cellular response to stimuli, phosphorylation, and signaling, whereas genes in the dark magenta module were involved in phosphorus metabolism and phosphate-containing compound metabolic processes. Although genes of both modules are expressed in the cytoplasm and vesicles, genes of the green module are expressed in membrane-bound organelles only. In the green module, the expressed proteins exhibited transferase activity and nucleic acid binding, whereas those in the dark magenta module exhibited catalytic activity and metal ion binding (Fig. 4B).

Construction of the gene co-expression network narrowed the range of candidate genes, and the 20 hub genes obtained by screening made the follow-up study more convenient.

DEG analysis

Differential gene analysis was performed on two modules of the GSE221542 dataset: severe drought stress and control and prolonged heat stress and control (Fig. 5A). The stress group with severe drought and long-term heat shock was selected for analysis (Control: GSM6883305-7, Drought: GSM6883299-301, Heat: GSM6883311-3). We found 484 and 1,559 had increased and decreased expression levels, respectively, in the drought stress group, whereas 1,876 and 3,158 DEGs were upregulated and downregulated, respectively, in the heat stress group.

GO analysis indicated that under drought or heat stress, macromolecule metabolism was the most altered biological process, and hormonal responses and other response activities were also altered. The proteins expressed by the two groups of DEGs were mainly binding proteins and mostly located in membrane-bound organelles and vesicles. However, the DEGs in the drought stress group were mostly related to ATP function, whereas those in the heat stress group were involved in nucleic acid and DNA-related functions (Fig. 5B).

The GO enrichment analysis results for the green module and heat stress–induced DEGs shared many similarities. In terms of biological processes, both involved cellular response activities. Regarding molecular functions and cell components, both involve a large proportion of proteins with nucleic acid binding that commonly act on membrane-bound organelles or vesicles. Therefore, the prediction of heat stress-related genes in the green module was expected to be more accurate (Figs. 4B and 5B). DEG analysis in the dataset is an important basis for the subsequent screening of key genes.

Further screening of hub genes

The top 20 hub genes in the dark magenta and green modules intersected with the DEGs under drought and heat stress, respectively, resulting in the selection of two candidate key genes associated with drought stress and 20 genes associated with heat stress (Fig. 5C).

Green module: The RNA-seq dataset used was GSE168650 (Kan et al., 2022), which contained RNA-seq data for two different genotypes of rice subjected to heat treatment and their corresponding controls. The data type was the RAW count. DEGs were analyzed using the same method without setting a threshold filter to identify key genes and their relative expression levels. GEO2R was used to analyze the expression levels of key candidate genes in multiple array datasets, including GSE136746 (Ps et al., 2017), GSE41648 (Sharma et al., 2021), and GSE14275 (Hu, Hu & Han, 2009). Cluster analysis was performed separately according to the original data types. Genes with high expression levels and consistent expression levels among different samples were selected from the heatmaps of both RNA-seq and array data (Figs. 6A–6B, Fig. S1A). Five key genes with the best overall performance were selected: OsDjC53, MBF1C, BAG6, HSP23.2, and HSP21.9.

Dark magenta module: the dataset GSE159816 (Zu et al., 2021) was downloaded, which contained two lines of rice subjected to drought treatment and their corresponding controls. We also analyzed the expression levels of key candidate genes in the GSE93917 (Wang et al., 2020) and GSE83378 (Wei et al., 2017) array datasets. The results showed that the expression levels of UGT83A1 and OsCPn60a1 did not show the same trend in multiple datasets; however, they were classified as key genes for further confirmation (Figs. 6C–6D, Fig. S2B).

Verification of key genes

qRT-PCR was used to verify changes in the expression levels of key genes in rice subjected to drought or heat stress conditions (Supplementary 2). The results showed that OsDjC53, MBF1C, BAG6, HSP23.2, and HSP21.9 were significantly overexpressed in rice under heat stress conditions (Fig. 7A), whereas the expression levels of UGT83A1 and OsCPn60a1 significantly decreased in rice under drought stress conditions (Fig. 7B). In summary, the five candidate genes in the green module may be the key genes associated with the heat stress response in rice.