Weighted gene coexpression network analysis-based identification of key modules and hub genes associated with drought sensitivity in rice

doi:10.21203/rs.3.rs-56040/v2

Download PDF

Research article

Weighted gene coexpression network analysis-based identification of key modules and hub genes associated with drought sensitivity in rice

https://doi.org/10.21203/rs.3.rs-56040/v2

This work is licensed under a CC BY 4.0 License

Journal Publication

published 20 Oct, 2020

Read the published version in BMC Plant Biology →

You are reading this older preprint version

Read the latest preprint version →

Background: Drought stress is an adverse factor with deleterious effects on several aspects of rice growth. However, the mechanism underlying drought resistance in rice remains unclear. To understand the molecular mechanism of the drought response in rice, drought-sensitive CSSL (Chromosome Single-substitution Segment Line) PY6 was used to map QTLs of sensitive phenotypes and to reveal the impact of the QTLs on transcriptional profiling.

Results: The QTL dss-1 was mapped onto the short arm of chromosome 1 of rice. According to transcriptomic analysis, the identified differentially expressed genes (DEGs) exhibited a downregulated pattern and were mainly enriched in photosynthesis-related GO terms, indicating that photosynthesis was greatly inhibited under drought. Further, according to weighted gene coexpression network analysis (WGCNA), specific gene modules (designating a group of genes with a similar expression pattern) were strongly correlated with H₂O₂ (4 modules) and MDA (3 modules), respectively. Likewise, GO analysis revealed that the photosynthesis-related GO terms were consistently overrepresented in H₂O₂-correlated modules. Functional annotation of the differentially expressed hub genes (DEHGs) in the H₂O₂ and MDA-correlated modules revealed cross-talk between abiotic and biotic stress responses for these genes, which were annotated as encoding WRKYs and PR family proteins, were notably differentially expressed between PY6 and PR403.

Conclusions: We speculated that drought-induced photosynthetic inhibition leads to H₂O₂ and MDA accumulation, which can then trigger the reprogramming of the rice transcriptome, including the hub genes involved in ROS scavenging, to prevent oxidative stress damage. Our results shed light on and provide deep insight into the drought resistance mechanism in rice.

Plant Molecular Biology and Genetics

Rice

Transcriptomic profiling

Drought-sensitive phenotype

WGCNA

photosynthesis inhibition

H2O2/MDA accumulation

DEGs encoding WRKYs/PR proteins

Abiotic stresses such as drought, chilling, heat, and salinity are widespread factors with deleterious effects on several aspects of plants, including changes in metabolism, growth, and development, in extreme cases leading to plant death. However, as sessile organisms, higher plant usually adopt the ‘overcome’ strategy upon encountering any extreme environmental stresses; this strategy, in contrast to that of animals, is preferentially to opt for the tolerance or avoidance of unfavorable circumstances. In general, various tolerance/avoidance mechanisms [1], such as regulation of endogenous hormone signals [2], increases in wax deposition in the cuticle [3], and stomatal regulation [4], are employed in response to abiotic stresses.

Rice (Oryza sativa L.) is a staple food crop worldwide. Drought has been one of the major abiotic constraints of rice production with increasing frequency of global water shortages. Two important strategies, dehydration tolerance and avoidance, are employed by plants to cope with drought stress. Dehydration avoidance depends on the development of a large and deep root system to uptake water from soil and reduce water loss via the modulation of stomatal opening. Dehydration tolerance mainly enhances drought resistance via complex mechanisms that occur during later stages of drought stress, such as antioxidation processes, ABA signal transduction, osmotic adjustment, and cell membrane protection [5, 6].

A large number of QTLs for root growth [7], osmotic pressure regulation [8], ROS regulation, and cuticular wax accumulation, which contribute to enhanced drought resistance, have been mapped to different locations of the rice genome via map-based cloning [9-11]. DRO1, one of the genes that has been well characterized, is involved in the regulation of rice root development under drought stress. DRO1 regulates cell elongation in the root tip, which causes downward bending of the roots to avoid water deficit by increasing deep rooting [11]. DST was identified as a negative regulator of H₂O₂ accumulation in guard cells and functions in the modulation of stomatal opening. Loss of DST function increases stomatal closure, resulting in enhanced drought tolerance in rice [10]. In addition, DS8 participates in the synthesis of plant cuticular wax and in the regulation of stomatal movement in response to drought stress [9]. Studies of these genes indicated that intricate mechanisms are employed by rice plants in the adaptation to drought stress. On the contrary, some drought stress-related QTLs were mapped through investigation of the drought sensitive traits. For example, 4 QTLs that contribute to increased transpiration rate under high vapour pressure deficit conditions were identified using F₇recombinant inbred lines of pearl millet [Pennisetum glaucum (L.) R. Br.]. All of these alleles were derived from drought-sensitive parent ICMB 841 [12]. Therefore, drought sensitivity is also an alternative trait for QTL identification in the study of drought tolerance mechanisms.

Chloroplasts are organelles that perform photosynthesis in plant cells. As a key process for the responsive regulation, when rice plants are subjected to abiotic stresses, photosynthetic activity is greatly inhibited [13]. The inhibition of photosynthesis reduces the utilization of absorbed light energy, resulting in the generation of reactive oxygen species (ROS), including singlet oxygen, H₂O₂, and superoxide anion in chloroplasts, which in turn can cause photodamage to photosystem proteins such as D1, a core component of photosystem II [14]. Plants have evolved a series of enzymatic and nonenzymatic antioxidant systems. The enzymatic system, which involves superoxide dismutase (SOD), catalase (CAT), and ascorbate peroxidase (APX), participates in scavenging ROS and preventing damage from oxidative stress [15]. SODs function in the dismutation of O²⁻ to H₂O₂ and O₂; the resultant H₂O₂ is then catalyzed into H₂O and O₂by the activity of CAT and APX [16]. The overaccumulation of ROS, causes severe injury to plant cells, even leading to programmed cell death. Therefore ROS homeostasis plays a key role in plant drought resistance regulation [18, 19]. It has been reported that the levels of antioxidant activity significantly enhanced with the increasing of POD, APX, CAT, SOD, and LOX enzymes in the drought tolerant genotypes of Persian walnut (Juglans regia L.) [20, 21]. Moreover, emerging evidence has indicated strong correlations among antioxidant mechanisms, photosynthesis regulation, and drought tolerance [18, 22, 23]. Specifically, ROS can act as secondary messengers that are involved in the stress signaling response, photosynthetic regulation, stomatal movement, and the development of higher plant [17]. In addition, ABA induces the accumulation of H₂O₂ via the activation of NADPH oxidase, which is activated by SnRK2 phosphorylation, eventually resulting in stomatal closure. However, on the other hand, the accumulation of ABA induces the expression of CAT, suggesting that a certain relationship exists between ABA signaling and H₂O₂balance [24-29].

Recent reports have suggested that WRKY transcription factors (TFs) are involved in the regulatory mechanism pathways of abiotic stresses, including those in response to drought and salinity [30, 31]. WRKY TFs constitute a plant-specific group of zinc finger transcription factors that are characterized by containing a conserved WRKY domain and that bind to a consensus cis-element W-box (TTGACT/C) [32]. The majority of the more than 100 WRKY members in rice are involved in the plant defense response in either a negative or positive manner [33-35]. Some WRKY members have important roles in abiotic stress responses; for example, OsWRKY13 has been suggested to regulate the antagonistic cross-talk between drought and disease resistance pathways via repression of SNAC1 and WRKY45-1 [30, 31]. Moreover, pathogenesis-related (PR) proteins have been well characterized as a group of proteins that are induced not only by infection from pathogens, such as viruses, bacteria, and fungi, but also in response to abiotic stresses, including drought and salinity [36]. All PR proteins in both monocots and dicotyledons can be classified into 17 families [37]. PR3 and PR4 proteins have plant chitinase activity, and the PR5 family of proteins includes permatins, osmotins, zeanatins, and thaumatin-like proteins [38]. PR10 genes belong to a multigene family in various plant species and exert RNase activity and ligand binding activity and are involved in posttranslational modification (phosphorylation) and phytohormone signaling in response to stress. The promoters of PR10 genes harbor cis-acting elements for WRKY, bZIP, ERF, and MYB transcription factors, suggesting PR10 genes have pivotal roles in the stress response pathways [38].

Given the rapid development of next-generation sequencing technology, RNA-seq has been widely used to explore the transcriptomic profiling of plants in response to challenges with adverse environmental conditions. Recent studies have revealed that antioxidant processes, the ABA pathway, and photosynthesis are the target biological processes for plants in response to environmental changes at the transcription level [39, 40]. Intriguingly, the photosynthesis rate was found to decline under drought stress, especially in drought-sensitive genotypes, suggesting that the enhancement of antioxidant capacity favors maintaining a higher photosynthetic activity and thereby improves the drought resistance of plants [40]. Moreover, phytohormone signaling and Ca²⁺ signaling have also been highlighted in transcript profiling and have been revealed cross-talk among phytohormone signaling pathways [41, 42]. In addition, a large number of transcription factor families, such as WRKY, bZIP, and MYB families, were found to be involved in ABA-mediated signal transduction under abiotic stress [43, 44]. However, the molecular mechanism of drought tolerance/sensitivity still needs to be further investigated.

Weighted gene coexpression network analysis (WGCNA) is a new technique that can be used to identify potential gene modules with the highest connectivity among genes, and these modules are correlated with certain phenotypes among the gene expression data [45, 46]. WGCNA involves the construction of a coexpression network that can indicate correlations among genes across samples. Based on the coexpression relationships, genes with similar expression are grouped into the same module, thus suggesting that genes in the same module may have similar functions or possibly have common biological regulatory roles [47]. This method has been successfully applied in numerous study cases that have been used to identify hub genes and to determine the relationships between gene expression data and relevant phenotypes [48-53].

Here, we identified a drought-sensitive rice CSSL PY6, which carries the substitution segment derived from the drought-sensitive variety Lambayeque 1, by backcrossing into the background of the drought-tolerant variety PR403. The dss-1 locus underlying the drought-sensitive phenotype of PY6 was harbored on the short arm of chromosome 1. To investigate the impact of dss-1 on the transcriptional profiling of PY6, RNA-seq was performed to identify genes that were differentially expressed under drought stress treatment. Further, WGCNA was applied to identify hub genes from the modules that were strongly correlated with drought-induced physiological indicators. By using this approach, we aimed to elucidate the dss-1-mediated potential regulatory network that is related to drought stress tolerance in rice.

Identification of the dss-1 locus for drought-sensitive phenotypes in CSSLs

A set of 140 rice CSSLs was constructed by introgressing the genomic segments of the drought-sensitive variety Lambayeque 1 into the genome of the drought-resistant variety PR403 through backcrossing (unpublished data; Fig. 1A). To elucidate the molecular regulatory mechanism underlying the drought-sensitive phenotype of Lambayeque 1, a drought stress treatment was applied to identify the major QTLs for drought-sensitive traits. PY6, one of the CSSLs, was highlighted for its drought-sensitive phenotype. After 10 days of drought treatment, the treated CSSLs were rewatered and allowed to recover for 3 days. In the screening, PY6, a CSSL, was identified as presenting a significant difference in survival rate (20.0%) compared with its recurrent parent PR403 (93.3%), whereas no significant difference was observed in the survival rate of PY6 and its donor parent Lambayeque 1 (Fig. 1B-C). During the drought stress treatment, the relative water content of the seedlings declined in PY6 and PR403 with decreasing soil moisture. However, starting at 120 h, a significant difference was observed between PY6 and PR403, and the former exhibited more severe dehydration at the later stage of the treatment, suggesting that drought stress disrupted the balance of water uptake and evaporation in PY6 (Fig. 1D). PY6 carries a 10-Mb chromosome segment derived from Lambayeque 1 on the short arm of chromosome 1. These results suggested that the segment from Lambayeque 1 may harbor a major QTL locus for the drought-sensitive phenotype of PY6. We designated this QTL as dss-1.

To reveal the genetic basis underlying the drought-sensitive phenotype afforded by the dss-1 locus, PY6 was backcrossed with the recurrent parent PR403 to construct a mapping population. The F₁ plants exhibited drought resistance, which was similar to that of the parental line PR403. Then, an F₂ population was generated to determine the inherited features of dss-1. In a 151-individual F₂ population, 111 members showed drought resistance, and 40 showed drought-sensitive traits (resistant:sensitive≈3:1, χ² = 0.1788 < χ² _0.05 = 3.84), indicating that drought-sensitive trait is controlled by a single recessive gene. In turn, the recombinants were screened for genetic linkage analysis, and dss-1 was ultimately narrowed down to a 4.01-Mb interval (Fig. 1E).

Drought induced H₂O₂ and MDA accumulation in the leaves of PY6

Given that previous reports indicated that drought stress was known to inhibit photosynthetic activity in plants due to an imbalance between light capture and its utilization [54] and coupled with the considerable potential for increased accumulation of superoxide and hydrogen peroxide in chloroplasts [55], we examined the accumulation of ROS and the activity of its scavenging enzymes in the leaves of drought-treated PY6 plants. Leaf samples were harvested at 4 sampling points (Fig. S1). In contrast to PR403, the increased accumulation of H₂O₂was detected in PY6 from points 40 to 15, whereas the enhanced content could be observed only at point 22 in PR403; therefore, the significantly increased accumulation of H₂O₂ was detected in PY6 at sampling points 22 and 15 compared with that of their respective controls (Fig. 2A). The marked accumulation of H₂O₂ in the leaves of PY6 was confirmed by DAB staining (Fig. 2B). Consistently, the activity of POD, an enzyme for H₂O₂scavenging, was upregulated in response to the accumulation of H₂O₂ (Fig. 2C). However, the activity of APX and CAT, another two major enzymes for H₂O₂ scavenging, was not induced, and no significant difference was observed between either genotype during drought treatment (Fig. S2A, B). For the dynamic changes in superoxide activity, the content of SOA (superoxide anion) and the activity of its scavenging enzyme SOD at point 15 exhibited significant differences between the two genotypes with higher values in PY6. However, the content of SOA in the whole period of drought treatment was not significantly greater than that in the untreated control (point 40) (Fig. 2D, E; Fig. S2C). Moreover, the contents of MDA, a major reactive aldehyde resulting from the peroxidation of polyunsaturated fatty acid (PUFA) constituents of biological membranes and a compound that acts as an indicator of abiotic stress-induced biomembrane damage, were significantly higher in PY6 than in PR403 at points 22 and 15 (Fig. 2F).

Taken together, our results suggest that the drought-sensitive phenotype may be caused by the overaccumulation of ROS, especially H₂O₂, which could result in severe damage to plant cells and peroxidation of the cell membrane, thereby leading to the dehydration of rice plants under drought stress.

Accumulation of ABA and stomatal modulation during drought stress treatment

As mentioned above, PY6 exhibited more severe dehydration at the later stage of treatment compared with the early stage, suggesting that PY6 either had a higher transpiration rate or lower water uptake ability when subjected to drought stress. ABA is a key regulator of plant stress responses and regulates a series of physiological processes, including tolerance to abiotic stresses and stomatal closure [56]. To determine the role of ABA-mediated stomatal movement in drought-induced severe dehydration in PY6, the ABA content and stomatal aperture status in drought-treated samples were examined. ABA accumulated during drought stress treatment, although no significant difference between PY6 and PR403 was observed at any sampling point (Fig. S2D). Scanning electron microscopy was used to examine the stomata status of leaf samples. Compared with that in the untreated control samples, the percentage of completely closed stomata in the drought-treated samples was significantly greater. However, PY6 and PR403 exhibited similar patterns in the ability to modulate their stomata in response to drought stress, and no significant difference in stomatal density was observed for either genotype (Fig. S2E). Therefore, these results suggested that severe dehydration in drought-treated PY6 do not result from a disruption of ABA-mediated stomatal transpiration regulation but, rather, may be caused by the limitation of water uptake or by other water loss pathways.

Comparative transcriptional profiling of PY6 under drought stress

To further reveal the impact of dss-1 on drought stress-induced H₂O₂ accumulation that may lead to the drought sensitivity of PY6, we investigated the transcriptional profiling of PY6 and PR403 in response to drought stress through RNA-seq. Principal component analysis (PCA) was initially performed based on all the transcriptional profiling data to detect the variations among the samples under the drought treatment (Fig. 3A). PC1 explained 37.98% of the total variation. The samples of CK40 and CK22/20 of both PY6 and PR403 were clearly separated by PC1, indicating the effect of the development imposed on the variation, while PC1 also separated the drought-treated samples of point 22 from those of point 20, indicating an effect of drought treatment on the transcriptomes of the samples. PC2 explained 27.94% of the total variation and clearly separated drought-treated samples from the parallel untreated samples (Fig. 3A). Moreover, both PY6 and PR403 samples from point 20 exhibited marked divergence from other sampling points, indicating that dramatic transcriptomic reprogramming occurred at this sampling point.

Drought-induced differentially expressed genes (DEGs) were identified by comparison with untreated controls (Table S2). It is worth noting that the development-dependent DEGs were detected in the same sampling points of well-watered samples under the same criteria, and they were excluded from the DEG set to avoid misinterpretation of drought-induced transcriptomic responses (Table S2-1~-8). Finally, at the two sampling points of PR403, 3052 DEGs, including 1135 upregulated and 1917 downregulated DEGs, were detected at point 22 (Table S2-9), while 5821 DEGs were detected at point 20, including 2383 upregulated and 3438 downregulated DEGs (Table S2-10). In PY6, 3145 (1748 upregulated/1397 downregulated) and 5529 (2655 upregulated/2874 downregulated) DEGs were detected at points 22 and 20, respectively (Table 1; Table S2-11, -12). Among the resultant DEGs, Venn diagram analysis revealed that 1142 DEGs were common between PR403 and PY6 at point 22, while 3184 were common between PR403 and PY6 at point 20, which accounted for 38.4% and 48.7% of the total DEGs at the two sampling points, respectively (Fig. 3B-C; Table 1). These findings demonstrated that PY6 has a similar genetic background and a drought stress response as PR403.

To verify the RNA-seq data quality, five genes were randomly selected for qRT-PCR. Similar trends of expression changes were observed between the qRT-PCR and RNA-seq data, indicating that the transcriptomic data were reliable for further analysis (Fig. S3A-E; Table S1).

Gene ontology enrichment analysis

GO enrichment analysis was performed to cluster the identified DEGs into the functional categories of biological process (BP), cellular component (CC), and molecular function (MF) subgroups. The results showed that all the DEGs were enriched in the GO terms that were related to photosynthesis in the chloroplasts (Fig. 4A-D, P£0.01). For instance, in terms of the BP group, the GO terms involving photosynthesis (GO:0015979); photosynthesis, light harvesting (GO:0009765); photosynthesis, light harvesting in photosystem I (GO:0009768); photosynthesis, light reaction (GO:0019684); chlorophyll biosynthetic process (GO:0015995); and porphyrin-containing compound biosynthetic process (GO:0006779) were highly significantly overrepresented. The overrepresented GO terms, such as those involving chloroplast part (GO:0044434), the chloroplast envelope (GO:0009941), chloroplast thylakoid (GO:0009534), the chloroplast thylakoid membrane (GO:0009535), and the chloroplast thylakoid lumen (GO:0009543) in the CC category (Fig. S4A-D, P£0.01), and the overrepresented terms of the MF category, such as pigment binding (GO:0031409) and chlorophyll binding (GO:0016168) (Fig. S5A-D, P£0.01), supported the identified DEGs are related to photosynthetic modulation in response to drought stress. In addition, the superoxide metabolism-related GO terms superoxide metabolic process (GO:0006801) and regulation of superoxide metabolic process (GO:0090322) were highly enriched among the DEGs identified at point 20 of PR403, suggesting the involvement of ROS scavenging processes in the drought stress response (Fig. 4C, P£0.01). GO enrichment analysis suggested that photosynthetic modulation plays an important role in the regulatory network of rice plants in response to drought stress; therefore, chloroplasts are target organelles for the drought stress response.

Construction of the coexpression network and identification of drought-induced hub genes

To further investigate the impact of dss-1, a major QTL conferred PY6 drought sensitive phenotype, on the regulatory network in response to drought stress and to identify the specific genes that are strongly correlated with drought-induced physiological alterations in rice, we performed a WGCNA. After removing the genes with a low FPKM (FPKM<1), a total of 23178 genes were used to construct a scale-free coexpression network based on the soft-thresholding power of β =12 (Fig. S6). According to the WGCNA results, the clusters with highly interconnected genes were defined as modules, and the genes in the same modules had high correlation coefficients. A total of 26 modules (coded with different colors to indicate different modules) were identified via the Dynamic Tree Cut method (core parameter: MEDissThres = 0.25) (Fig. 5A). With respect to the correlations of physiological traits with overrepresented modules, intriguing results were observed in the correlations. In general, the modules black, blue, and tan were negatively correlated with the majority of physiological traits, whereas the modules royalblue, brown, red, grey60, orange, and green were positively correlated with the traits, indicating that the genes clustered in a module have similar altered expression patterns in response to drought stress (Fig. 5B). Considering that significant differences in the contents of H₂O₂ and MDA were observed during drought stress, we focused on the correlations of H₂O₂ and MDA with gene modules and identified hub genes that were strongly correlated with the two drought-related physiological indicators (Fig. 2A, B, F; Fig. 5B, C, D). With the cut-off threshold of GS (Gene Significance) > 0.4, in all 26 of the modules, black (741 genes; R²=-0.69, P=3.0 ´ 10^-5), blue (4810 genes; R²=-0.57, P=9.0 ´ 10^-⁴), grey60 (148 genes; R²=0.67, P=4.0 ´ 10^-5), and green (1099 genes; R²=0.61, P=3.0 ´ 10^-⁴) were considered the modules that were strongly correlated with H₂O₂, and red (913 genes; R²=0.63, P=2.0 ´ 10^-⁴), brown (3132 genes; R²=0.56, P=1.0 ´ 10^-³), and royalblue (45 genes; R²=0.68, P=4.0 × 10^-5) modules were considered strongly correlated with MDA (Fig. 5B, D). The genes in which the modules were clustered showed strong correlation with their respective modules (Fig. S7A, B). Furthermore, hub genes were identified based on the following cut-off thresholds: GS > 0.60 and MM (Module Membership) > 0.80 for the modules black, blue, green, and royalblue, and GS > 0.5 and MM> 0.8 for the modules grey60, red, and brown according to the threshold strengthening. After removing the genes with |log₂FC|<1, the expression of all the obtained hub genes was significantly differentially altered in PY6 and PR403 under drought stress and were highlighted in subsequent functional analyses (Fig. S8A-F; Table S3-1~-7). Hereafter, we designated these hub genes as DEHGs.

Functional enrichment analysis of hub genes correlated with H₂O₂ accumulation

GO enrichment analysis was performed to functionally cluster the hub genes in the modules that were strongly correlated with H₂O₂. A total of 115 and 204 DEHGs were identified in the modules black and blue, respectively (Table S3-1, -2). In line with the abovementioned results, these hub genes were highly enriched in the GO terms involving chloroplast (GO:0009507); the chloroplast envelope (GO:0009941); thylakoid (GO:0009579); chloroplast stroma (GO:0009570); chlorophyll binding (GO:0016168); photosynthesis (GO: GO:0015979); photosystem I (GO:0009522); photosynthesis, light harvesting (GO:0009765); and photosynthesis, light reaction (GO:0019684). These terms are related to photosynthesis (Qvalue £ 0.05) (Fig. S9A-B; Table S5-1, -2). Interestingly, in agreement with previous reports [57], the expression of the majority of the DEHGs in the modules black and blue was downregulated without marked log₂FC differences (|log₂FC_PY6-log₂FC_PR403|≥1, P<0.05) observed between PY6 and PR403 during drought treatment (Fig. 6A, B; Fig. S8A, B), implying that PY6 and PR403 have a common drought-responsive regulation and that their photosynthetic activities were inhibited in response to drought stress.

The module grey60 was enriched with 41 DEHGs (Table S3-3). The expression of most of these DEHGs was elevated under drought stress (Fig. S8C). Among these genes, eight showed differential expression with markedly different magnitudes at least at one sampling point between two genotypes, exhibiting a strong dynamic induction in response to drought treatment (|log₂FC_PY6-log₂FC_PR403|≥1, P<0.05) (Fig. 6C; Table 2). The expression of Os01G0842500 and Os03G0273200, both of which encodes proteins that are similar to laccase and both of which were enriched in GO terms involving the apoplast (GO:0048046)/lignin catabolic process (GO:0046274)/cytoplasmic membrane-bounded vesicle (GO:0016023)/oxidoreductase activity (GO:0016491), was consistently enhanced at two points of both PY6 and PR403 and showed higher magnitudes in PY6 (particularly at point 20) (Fig. 6C; Fig. S9C; Table 2; Table S5-3). Five genes (Os01G0702700, Os06G0131700, Os05G0543600, Os03G0433200, and Os04G0508500), which are annotated as transcription factors and mainly enriched in GO terms involving DNA binding (GO:0003677), sequence-specific DNA binding transcription factor activity (GO:0003700), and response to freezing (GO:0050826), were significantly consistently induced to a certain extent at point 22 and were downregulated at point 20 in PR403, whereas they were significantly induced only at point 20 in PY6, with the exception of Os03G0433200 (Fig. 6C; Table 2; Table S5-3). Os01G0702700, Os05G0543600, and Os04G0508500 are MYB family TFs, and Os06G0131700 is a NAC TF. A recent study suggested that Os06G0131700 (OsSWN1) is related to secondary cell wall formation [58], implying that secondary cell wall formation may play a certain role in drought stress resistance. In addition, Os08g0189200 is annotated as a gene encoding Germin-like protein 8-3 and was suggested to be involved in disease resistance [59]. A different upregulation pattern of this gene was observed in PY6 and PR403 in response to drought stress (Fig. 6C; Table 2).

A total of 96 hub genes were identified in the module green, including 36 downregulated and 60 upregulated genes in response to drought (Fig. S8D; Table S3-4). GO enrichment analysis revealed that the majority of the overrepresented GO terms were related to oxidation-reduction process (GO:0055114), response to salt stress (GO:0009651), carbohydrate metabolic process (GO:0005975), biosynthetic process (GO:0009058), and metabolic process (GO:0008152). These hub genes were sublocalized to distinct cellular components and have various molecular functions (Fig. S9D; Table S5-4). However, with the exception of iron-sulfur cluster binding (GO:0051536), no GO term was significantly overrepresented in our analysis (Table S5-4). It is worth noting that only 1 hub gene, OS04G0685300, which encodes a protein with harpin-induced 1 domain, showed significant differential expression with a greater upregulation magnitude at two points of PY6 (2.03/2.67) compared with those in PR403 (Table S3-4).

Taken together, the functional annotation of the hub genes identified in the H₂O₂-related modules suggested that both PY6 and PR403 experienced inhibition of drought-triggered photosynthetic activity, leading to the accumulation of H₂O₂in the rice plants. The hub genes with highly differential expression identified in the module grey60 may be responsible for the accumulation of H₂O₂ in rice plants.

Functional enrichment analysis of hub genes correlated with MDA accumulation

MDA is an important indicator of membrane lipid peroxidation. The modules red, brown, and royalblue were strongly correlated with MDA under drought treatment; there were 303, 304, and 13 hub genes identified in these modules, respectively (Table S3-5~-7). For the module red, GO analysis revealed that the hub genes are involved in biotic stress and defense responses (Fig. S10A; Table S5-5). Among these hub genes, 29 showed differential expression at least at one sampling point between two genotypes (|log2FC_PY6-log2FC_PR403|≥1, P<0.05) (Fig. 7A; Table 3). Strikingly, most of them consistently showed similar expression patterns in both genotypes, with elevated transcripts at point 22, after which the expression levels then declined to levels that were comparable to or even lower than those of their respective untreated controls (Fig. S8E; Table 3). GO analysis revealed that these genes were involved in biotic stress responses; for example, the GO terms involving defense response (GO:0006952), response to biotic stimulus (GO:0009607), response to stress (GO:0006950), and response to endogenous stimulus (GO:0009719) were overrepresented in the results (Fig. 7A; Table S5-5). Of these genes, 8 were WRKY family TFs (Fig. 7A; Table 3). Specifically, OsWRKY70 (Os05g0474800) and OsWRKY76 (Os09g0417600) have been suggested to participate in ABA signaling and the biotic stress response [31, 60]. OS01g0185900 and OS05g0322900 exhibited relatively higher upregulated expression at point 22, whereas two other genes, OS05g0571200 and OS11G0117600, had relatively higher expression at point 20 of PY6 than at PR403 (Fig. 7A; Table 3). Considering the accumulation of MDA in PY6 and PR403, these results implied that four genes may be positively correlated with the overaccumulation of MDA and lead to the death of PY6 under drought stress. Moreover, PR-10a (OS12g0555300), which functions at the downstream of jasmonic acid pathway and is positively regulated by WRKY TFs in response to drought and high-salt stresses [36, 61], and three other PR-10 protein family hub genes (Os12g0555200, Os12g0555500, and Os12g0555000), which were annotated as genes encoding probenazole-inducible PBZ1 proteins and are involved in disease resistance, exhibited a higher magnitude of upregulated expression at two sampling points of PY6 compared with PR403 (Fig. 7A; Table 3). BSR-d1/ZFP36 (Os03g0437200), which encodes a C₂H₂-type zinc finger protein, participates in ABA-OsMPK transduction, which results in the accumulation of H₂O₂[62, 63]. Similarly, this gene also had a higher upregulated expression in PY6 than in PR403 (Table 3). In addition, the expression of OS02g0121700 and OS02g0570400 was dramatically induced in PY6 during the whole period of drought treatment, whereas these genes were significantly induced only at point 22 of PR403, after which their expression levels decreased to those of the untreated control (Fig. 7A; Table 3). These genes encode a terpenoid synthase domain-containing protein and ent-kaurene synthase 1A, which are involved in gibberellin (GA) biosynthesis, suggesting a role of GA in the drought stress response in rice. OS04g0109100, a gene annotated as concanavalin A-like lectin/glucanase, had a similar expression pattern in both genotypes.

In total, 18 and 13 DEHGs were identified in the modules brown and royalblue, respectively (Fig. 7B, C; Table 4; Table S4). GO analysis revealed that the DEHGs in the module brown play different roles in various biological processes, including response to stress (GO:0006950), response to endogenous stimulus (GO:0009719), protein phosphorylation (GO:0006468), protein metabolic process (GO:0019538), lipid metabolic process (GO:0006629), carbohydrate metabolic process (GO:0005975), protein folding (GO:0005515), regulation of transcription (GO:0006355), and signal transduction (GO:0007165) (Fig. 7B; Table S5-6). Nearly all of these genes were induced at point 22 and were repressed at point 20 in response to drought in both genotypes. Although they have a similar expression pattern, a different extent of downregulation was observed in several hub genes; for example, OS11G0569300, OS01G0971800, OS04G0675400, and OS09G0243200 showed a lower extent of downregulation at point 20 of PY6, while OS07G0201500 and OS11G0209600 showed a lower extent of downregulation at point 20 of PR403 (Fig. 7B; Table S4). In the module royalblue, only one GO term was significantly overrepresented due to the relatively low number of DEHGs. The striking feature of the expression pattern of the hub genes was that their transcripts were dramatically induced by drought at point 22 of PY6, and their levels were also slightly reduced at the next sampling point. Interestingly, Os04g0494100 and Os11g0592200 were annotated as genes encoding PR3 and PR4 family proteins, which have plant chitinase activity. Os03g0663600 and Os12g0630200 are two genes encoding pathogenesis-related thaumatin-like proteins, which are members of the PR5 family, suggesting that these genes played roles in the cross-talk between drought and biotic stresses in this study (Fig. 7C; Table 4; Table S5-7).

In this study, CSSL PY6 was used to characterize a QTL locus, dss-1, for its drought-sensitive phenotype and to investigate the impact of dss-1 on the reprogramming of transcriptional profiles of PY6 in response to drought stress via RNA-seq and WGCNA.

Drought-induced overaccumulation of H₂O₂ may lead to the sensitive phenotype of dss-1

Reactive oxygen species (ROS) (e.g., O^2.-, H₂O₂, OH˙, ¹O₂) are unavoidable toxic byproducts of plants generated in response to abiotic/biotic stresses [64]. The overaccumulation of ROS in cells causes severe oxidative damage to membranes (lipid peroxidation), proteins, RNA, and DNA molecules [65]. Higher plant has thus evolved dedicated scavenging pathways to protect themselves from ROS toxicity, including pathways involving detoxifying enzymes such as CAT, SOD, POD, and APX, as well as the antioxidant ascorbate-glutathione (GSH) cycle [16, 66]. In this study, significant overaccumulation of H₂O₂ was detected in PY6 compared with PR403 during drought treatment (Fig. 2A, B). Regarding the drought-sensitive phenotype of PY6, which was afforded by the introgression of dss-1 into the background of PR403, our results suggested that important relationships were evident among dss-1, H₂O₂ accumulation, and the sensitive phenotype. That is, we speculated that the introgression of dss-1 might directly or indirectly lead to the overaccumulation of H₂O₂ in rice cells. The overaccumulation of H₂O₂ causes oxidative damage to the cell membrane (lipid peroxidation), which was reflected by the measurements of MDA and others, subsequently leading to the death of rice plants under severe drought stress conditions.

Photosynthesis regulation, a common regulatory process for rice in response to drought stress

Chloroplasts are important organelles for the fixation of light energy necessary for the biological activity of higher plant and all other life forms in our biosphere [67]. Beyond their normal role in photosynthesis, increasing evidence suggests that chloroplasts are targeted organelles that are involved in environmental stress responses by positively or passively tuning photosynthetic activity to adapt to environmental changes [22]. According to the GO enrichment analysis, the drought-induced DEGs identified either in PR403 or PY6 were mainly enriched in photosynthesis-related GO terms (Fig. 4; Fig. S4; Fig. S5). Additionally, the analysis of differentially expressed hub genes that clustered in the modules black and blue showed that a similar set of GO terms that are related to the photosynthesis process in chloroplasts were overrepresented (Fig. S9A, B). In terms of the altered expression during drought treatment, in line with the results of previous studies [57], the majority of hub genes clustering in the modules black and blue consistently declined to a similar extent in both PY6 and PR403 (Fig. S8A, B), implying that the photosynthetic activity of both genotypes was dramatically inhibited in response to drought stress. The inhibition of photosynthesis reduces the utilization of absorbed light energy; excess light energy results in the generation of toxic ROS [14]. Therefore, photosynthetic activity regulation could be a common regulatory process that is shared among diverse varieties for the adaptation of rice in response to environmental stresses. The inhibition of photosynthesis in both genotypes of this study might be correlated to the accumulation of H₂O₂ in rice cells. The differential accumulation of H₂O₂ in the leaves of PY6 and PR403 results from the effects of the introgression of dss-1 into PR403 on the different expression changes of the hub genes identified in the related modules. Furthermore, previous studies considered that the improvement of photosynthetic capability was important for drought resistance of plants [40]. However, in this study, no significant expression difference was observed in the photosynthesis-related DEGs identified in either PR403 or PY6, suggesting that the dss-1-induced drought-sensitive phenotype was not caused by the severe inhibition of photosynthesis but, rather, by overaccumulation of H₂O₂.

ROS-correlated hub genes that were differentially regulated in response to drought stress

Transcription factors are key players in the regulatory networks of plants in response to unfavorable stresses. The MYB family is one of the largest transcription factor families in plants; its members have a conserved MYB-binding domain, and the MYB family TFs have been shown to be induced by H₂O₂ in soybean and play essential roles in response to abiotic stress [68, 69]. In this study, five DEGs annotated as transcription factors, including 3 MYBs, 1 NAD, and 1 SHR, were highly clustered in the module grey60, which was statistically correlated with the accumulation of H₂O₂ during drought stress. Among the genes, Os01G0702700, Os05G0543600, and Os04G0508500 are MYB family transcription factors and showed a higher magnitude of upregulation at point 22 of PR403. Considering the relatively low accumulation of H₂O₂ and MDA in PR403, MYB family TFs may be induced by drought stress or other signaling messengers and may trigger the activation of downstream components of drought-responsive signaling pathways. This action may also include the upregulation of stress-responsive genes involved in increasing the content of osmo-protectants and the activities of POD, CAT, and SOD, as well as maintaining a sublethal content of H₂O₂ in plant cells during drought stress, suggesting that the MYB TFs identified here were positively correlated with drought resistance (Fig. 2; Fig. 6C; Table 2) [70]. Moreover, WRKY TFs, a family of transcription factors involved in plant defense responses [71], were found to be strongly correlated with MDA and 8 differentially expressed WRKY TFs identified in the module red showed a higher degree of induction at point 22 and a lower level at point 20 in PY6 compared with PR403 (Table 3).

PR proteins are another group of proteins that is involved in plant defense responses [36]. In this study, four PR10 family encoding genes identified in the module red and three other families of PR proteins (PR3, PR4, and PR5) identified in the module royalblue exhibited similar altered expression patterns in response to drought stress. Among the identified PR-encoding genes, OS12G0555000 (RSOsPR10), a gene encoding the rice root-specific pathogenesis-related protein PBZ1, was previously reported to be induced by drought and salinity stresses, as well as jasmonic acid; however, it was strongly inhibited by salicylic acid [36, 61]. Besides, the expression of RSOsPR10 was postulated to be repressed by SA-induced OSTGAs or OsWRKYs [36]. Our results suggested that cross-talk occurs between abiotic and biotic stress-responsive pathways during drought stress. Recent reports have revealed that some common components are shared by the signaling pathways involved in abiotic and biotic stress responses: the phytohormones abscisic acid, salicylic acid, and jasmonic acid; cis-acting regulatory elements; and components of protein kinase cascades [2]. Plants can rapidly respond to environmental changes by the synergistic or antagonistic regulation of cross-talk among pathways via their convergent nodes [2]. OsWRKY76 (OS09G0417600), a hub gene of the module red, encodes a group IIa WRKY TF in rice. It has been shown that overexpression of OsWRKY76 in rice plants suppresses a specific set of PR genes and increases the expression of abiotic stress-associated genes, which resulted in increased susceptibility to Magnaporthe oryzae (M. oryzae) and improved tolerance to chilling stress [31]. In addition, OsWRKY13 was suggested to regulate the antagonistic cross-talk between drought and disease resistance pathways via the repression of SNAC1 and WRKY45-1 [30]. Thus, the WRKY TF- and PR family protein-encoded genes identified here may act as convergent nodes and may be responsible for the cross-talk between abiotic and biotic stresses.

In addition, strikingly, Os01G0842500 and Os03G0273200, both of which were annotated as putative laccase-encoding genes, consistently exhibited elevated expression in both genotypes in response to drought. However, PY6 had a relatively higher upregulated expression at point 20, implying that the overaccumulation of the transcripts of the genes for laccase may be correlated with the drought-sensitive phenotype of PY6. It has been demonstrated that a rice laccase-encoding gene, OsChi1, has been reported to play an important role in the ROS signaling pathway [74], suggesting that, in the present study, MYB TF-regulated laccases with differential expression are involved in the improvement of drought stress tolerance (Fig. 8). However, paradoxically, the overexpression of another laccase-encoding gene, OsChi1 (Os01g0827300), in Arabidopsis resulted in increased tolerance to drought and salinity stresses [74]. Thus, further study needs to be performed to explore the function of the laccases identified in the study and determine the discrepancy. Moreover, ZFP36 (OS03G0437200), a gene encoding a C₂H₂ zinc finger protein (Brs-d1), is a key factor for regulating rice cell redox status and is involved in ABA-induced upregulation of the expression and activity of SOD and APX. Overexpression of ZFP36 enhanced tolerance to drought and oxidative stresses but reduced resistance to rice blast [62, 63]. However, dramatically increased expression of ZFP36 was exhibited in the drought-sensitive genotype PY6, suggesting that this gene may play a synergistic role in the regulatory network with WRKY TFs and PRs in rice plants in response to drought stress.

Impact of dss-1 on the reprogramming of hub genes in response to drought stress

Taken together, a working model was proposed here to depict the possible dss-1-mediated mechanism of rice plants in response to drought stress at the transcriptional level (Fig. 8). When rice is exposed to severe drought stress, stress-induced dehydration of rice plants and stomatal closure of leaf guard cells result in the inhibition of photosynthesis activity and the enhancement of ABA biosynthesis. The inhibition of photosynthesis reduces the utilization of absorbed light energy and results in the generation of toxic ROS, including H₂O₂. Overaccumulation of ROS can cause oxidative stress damage to plant cells, leading to the drought-sensitive phenotype of rice. Drought stress signals, ABA accumulation, and H₂O₂ signals are integrated and perceived by rice cells and trigger the activation of responsive signaling pathways. First, the expression of a series of TFs, including MYBs, WRKYs, ZFP36, NACs, etc., was differentially reprogrammed, and these TFs induced cross-talk between drought and biotic stress responses. Induction of the TFs then regulates the expression of downstream responsive genes, including a number of PR family proteins and laccases, to mobilize multiple aspects of in vivo resources to respond to drought stress, and eventually reducing the content of ROS in plant cells and ameliorating the drought-induced sensitive phenotype of rice. The introgressed dss-1 may exert its effect on the hub genes that were identified in our study to differentially alter their expression, therefore leading to the differential accumulation of ROS between drought-resistant and drought-sensitive varieties. The differential activation of ROS-related responsive pathways could cause the differences between the drought resistance/sensitive phenotypes of rice varieties.

In this study, we used a drought sensitive plant material to fine map QTL dss-1 and investigate its effect on the transcriptomic profiling of rice. In subsequent, we will identify a drought sensitive QTL allele from Lambayeque1 by positional cloning; the drought tolerant allele of dss-1 from PR403 can also be isolated as well. They will then be studied for the elucidation of drought tolerant mechanisms and drought tolerant allele DSS-1 will be applied for rice breeding.

Taken together, our results suggest that chloroplasts are the targeted organelles for PR403 and PY6 in response to drought stress regulation. Functional annotation analysis of the hub genes in the modules that were strongly correlated with H₂O₂ and MDA accumulation during drought treatment provided us with a global picture of the impact of dss-1 on the reprogramming of the regulatory mechanism(s) at the transcriptional level. Specifically, the differential expression patterns of MYBs, WRKYs, PR protein encoding genes, ZFP36, etc., in both genotypes implied their roles in the variation in ROS and MDA accumulation, thereby resulting in the sensitive phenotype of rice in response to drought stress.

Plant materials, growth, and treatment conditions

PY6 was screened in a set of CSSLs that were constructed by introgression of the genomic fragments of donor parent Lambayeque 1 into recurrent parent PR403. Rice seedlings were grown according to previously described methods [13], with minor modifications. In brief, rice plants were grown in equal amounts of paddy soil. Additionally, an equal amount of water was applied routinely. The drought treatment was performed according to previously described methods [23]. In the 17-day of drought treatment, soil moisture was monitored regularly to determine the stages of drought stress using an EM-WSYP soil moisture detector (Hengmei Company, China). The penultimate leaves of 3-week-old soil-cultured rice seedlings (three biological replicates) grown in the greenhouse (35 ℃ day/30 ℃ night) were collected at four soil moisture points, 40 (40% soil moisture, representing well-watered conditions), 22 (22% soil moisture, drought treatment for 96 h after 40% soil moisture conditions), 20 (20% soil moisture, drought treatment for 120 h), and 15 (15% soil moisture, drought treatment for 168 h) (Fig. S1), and flash frozen in liquid nitrogen for RNA-seq and physiological index measurements. The samples harvested at point 15 were only for physiological measurement analyses. Drought stress treatment for hydroponically cultivated seedlings in 96-well plates was performed in rice nutrient solution supplemented with 20% (w/v) PEG 6000. The survival rate was determined after 5 days of treatment and rewatering for 3 days.

RNA-seq library construction and transcriptomic data processing

A transcriptome library was generated by a NEBNext^® Ultra^TM RNA Library Prep Kit for Illumina^® (NEB, USA), following the manufacturer’s instructions. After digestion, purification, and PCR amplification, the cDNA fragments (250~300 bp) were purified by AMPure XP beads (Beckman, USA). cDNA sequencing was performed on an Illumina 2500 platform, and paired-end reads were generated. After removing low-quality reads from the data, the read depth of 6 Gb-clean reads that is sufficient for optimal transcriptome coverage universally were obtained and mapped to the reference genome (the Rice Annotation Project Database, https://rapdb.dna.affrc.go.jp/index.html). DEGs between each treatment point and untreated controls were detected using DESeq, with |log2FC| ≥1 (FC, the abbreviation for fold change) and padj £0.05 used as thresholds to identify drought-induced DEGs. Additionally, the development-dependent DEGs were detected in the same sampling points of well-watered samples with the same criteria. Because the development-dependent DEGs could be caused by the development of rice plants, they were excluded from the final dataset.

GO enrichment analysis was performed via PlantRegMap (http://plantregmap.cbi.pku.edu.cn/). The overrepresented GO terms were visualized using ggplot2 and the GOplot R package [75]. A heatmap was generated by MEV v4.9 [76].

Coexpression network analysis

The WGCNA R package was used to construct a coexpression network and analyze the correlation of modules and physiological data according to a previously described method [45]. Raw transcriptomic datasets of all samples were filtered to remove all genes with an FPKM (Fragments Per Kilobase Million) <1, even for a single replicate of any sampling point of the samples. Clean data were used to generate a coexpression network. First, the soft threshold power was estimated by the pickSoftThreshold function, which provides the appropriate soft-thresholding power for network construction. Second, a correlation matrix was constructed based on the obtained soft threshold power, and afterward, a topological overlap matrix (TOM) was calculated from the transformed correlation matrix. Finally, the genes were grouped based on the topological overlap dissimilarity (1-TOM) by average hierarchical clustering using the hclust function. Gene modules were then identified using a dynamic tree cutoff algorithm (minimum cluster size of 30, merging threshold function of 0.25). Module membership (MM) was calculated based on Pearson correlations between the expression level and the module eigengenes to identify hub genes within the modules. A relatively high MM indicates that those genes have a relatively high connectivity within the module.

To correlate the physiological data with the network, module eigengenes were correlated with each physiological data point. Gene significance (GS) was used to correlate the physiological data with the expression data of individual genes.

Analysis of the movement of the stomatal aperture

Four-week-old seedlings were selected for scanning electron microscopy (SEM) analysis. The middle part of penultimate leaves was cut into 0.5-cm pieces, which were subsequently flash frozen in liquid nitrogen to maintain stomatal morphology. An FEI Quattro S scanning electron microscope (Thermo Fisher Scientific, USA) with a freezing stage was used to examine the leaf stomatal apertures. The opening status and density of stomata were analyzed with ImageJ software. Four random fields were selected for counting the number of stomata in each replicate.

Quantitative real-time reverse transcription PCR analysis

Total RNA was isolated from leaf samples (three biological replicates per sample). cDNA synthesis was performed by reverse transcription (RT) with a Thermo Scientific Revert Aid First Strand cDNA Synthesis Kit (Cat# K1622) according to the manufacturer’s protocol. The sequences of genes were downloaded from RAP-DB (https://rapdb.dna.affrc.go.jp/index.html). The sequences of the primers of the relevant genes were obtained from qPrimerDB (https://biodb.swu.edu.cn/qprimerdb/) (Table S1) [77]. qPCR was performed on a Roche LightCycler 480 Real-Time PCR System in 10 μL reactions with a SYBR green PCR Master Mix Kit (BIO-RAD, USA), following the manufacturer’s protocol. The relative expression of each gene was calculated according to the 2^-^DD^CT method [78]. The actin gene (Os11g0163100) was used as an internal reference for analysis.

Determination of the relative water content of the leaves

Leaf samples were collected at each sampling point and weighed to determine the fresh weight (FW). The samples were then oven-dried to a constant dry weight (DW). The relative water content was subsequently calculated according to the following equation: RWC (%) = [(FW– DW) / FW] × 100%.

Examination of ABA, ROS and enzyme scavenging activity

Fresh leaf samples (0.15 g) of PR403 and PY6 were collected and ground into a fine powder in liquid nitrogen, after which they were resuspended in 1.35 mL of PBS buffer (10 mM, pH 7.2). After centrifugation at 5000 rpm for 5 min, the supernatant was collected for subsequent assays. The contents and activity of H₂O₂, SOA, ABA, SOD, CAT, POD, and APX were measured using an enzyme-linked immunosorbent assay (ELISA) kit (MSKBIO, China) following the manufacturer’s instructions.

3,3′-Diaminobenzidine (DAB) staining was performed to detect the accumulation of H₂O₂ in the leaves. The leaves of 3-week-old seedlings were collected for staining. The samples were subsequently vacuum-infiltrated for 10 min and then incubated in DAB solution (1 mg/mL in water, pH 5.8) for 12 h in the dark. For the detection of superoxide, leaf pieces were vacuum-infiltrated as described above in 0.1% nitro blue tetrazolium (NBT) together with 0.1% Triton X-100 in 10 mM phosphate buffer (pH 7.2) for 10 min. The samples were then boiled in absolute ethanol to remove chlorophyll [13].

Availability of data and materials

All data generated or analysed during this study are included in this published article [and its supplementary information files].

Acknowledgments

We thank Lijun Luo, a senior researcher at Shanghai Academy of Agricultural Sciences, for providing the parental rice seeds. We thank the reviewers and editors for the careful reading of our manuscript and constructive comments, and we appreciate the technical assistance concerning the scientific instrument platform of the State Key Laboratory for Conservation and Utilization of Subtropical Agro-bioresources (Guangxi University) and are grateful for the staff performing the SEM analysis in this study.

Funding

This research was supported by grants from the National Natural Science Foundation of China (CN) (31671646), the Guangxi Hundred-Talents Program (2015), the Guangxi innovation-driven development special funding project (Guike-AA17204070), and the State Key Laboratory for Conservation and Utilization of Subtropical Agro-bioresources (SKLCUSA-a201907, -a201918, and -a201801). The funding bodies had no role in the design of the study and collection, analysis, and interpretation of data and in writing the manuscript.

Author contributions

JL and RL conceived and designed the experiments. BY, JBL, DW, YL, and WC performed all the lab experiments. BY, JBL, and DW conducted the field trials and collected samples for RNA-seq. JL, BY, and WC processed the data. JL drafted the manuscript. JL, RL, BY, and SW revised the manuscript. All authors approved the final version of the manuscript.

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

CSSL: Chromosome single-substitution segment line; DEGs: Differentially expressed genes; MDA: Malondialdehyde; WGCNA: Weighted gene coexpression network analysis; ROS: Reactive oxygen species; SOD: Superoxide dismutase; CAT: Catalase; APX: Ascorbate peroxidase; TFs: Transcription factors; PR: Pathogenesis-related; POD: Peroxidase; PUFA: Polyunsaturated fatty acid; PCA: Principal component analysis; BP: Biological process; CC: Cellular component; MF: Molecular function; DEHGs: Differentially expressed hub genes; TOM: Topological overlap matrix; MM: Module membership; GS: Gene significance; SEM: Scanning electron microscopy; FW: Fresh weight; DW: Dry weight; DAB: 3,3′-Diaminobenzidine; NBT : Nitro blue tetrazolium

de Zelicourt A, Colcombet J, Hirt H: The Role of MAPK Modules and ABA during Abiotic Stress Signaling. Trends Plant Sci 2016, 21(8):677-685.
Verma V, Ravindran P, Kumar PP: Plant hormone-mediated regulation of stress responses. BMC Plant Biol 2016, 16:86.
Zhu X, Xiong L: Putative megaenzyme DWA1 plays essential roles in drought resistance by regulating stress-induced wax deposition in rice. Proc Natl Acad Sci U S A 2013, 110(44):17790-17795.
Huber AE, Melcher PJ, Pineros MA, Setter TL, Bauerle TL: Signal coordination before, during and after stomatal closure in response to drought stress. New Phytol 2019, 224:675-688.
Hu H, Xiong L: Genetic engineering and breeding of drought-resistant crops. Annu Rev Plant Biol 2014, 65(1):715-741.
Luo LJ: Breeding for water-saving and drought-resistance rice (WDR) in China. J Exp Bot 2010, 61(13):3509-3517.
Lou Q, Chen L, Mei H, Wei H, Feng F, Wang P, Xia H, Li T, Luo L: Quantitative trait locus mapping of deep rooting by linkage and association analysis in rice. J Exp Bot 2015, 66(15):4749-4757.
Lilley JM, Ludlow MM, McCouch SR, O'Toole JC: Locating QTL for osmotic adjustment and dehydration tolerance in rice. J Exp Bot 1996, 47:1427-1436.
Huang L, Chen L, Wang L, Yang Y, Rao Y, Ren D, Dai L, Gao Y, Zou W, Lu X et al: A Nck-associated protein 1-like protein affects drought sensitivity by its involvement in leaf epidermal development and stomatal closure in rice. Plant J 2019, 98(5):884-897.
Huang XY, Chao DY, Gao JP, Zhu MZ, Shi M, Lin HX: A previously unknown zinc finger protein, DST, regulates drought and salt tolerance in rice via stomatal aperture control. Genes Dev 2009, 23(15):1805-1817.
Uga Y, Sugimoto K, Ogawa S, Rane J, Ishitani M, Hara N, Kitomi Y, Inukai Y, Ono K, Kanno N et al: Control of root system architecture by DEEPER ROOTING 1 increases rice yield under drought conditions. Nat Genet 2013, 45(9):1097-1102.
Aparna K, Nepolean T, Srivastsava RK, Kholová J, Rajaram V, Kumar S, Rekha B, Senthilvel S, Hash CT, Vadez V: Quantitative trait loci associated with constitutive traits control water use in pearl millet [Pennisetum glaucum (L.) R. Br.]. Plant Biol 2015, 17(5):1073-1084.
Cen W, Liu J, Lu S, Jia P, Yu K, Han Y, Li R, Luo J: Comparative proteomic analysis of QTL CTS-12 derived from wild rice (Oryza rufipogon Griff.), in the regulation of cold acclimation and de-acclimation of rice (Oryza sativa L.) in response to severe chilling stress. BMC Plant Biol 2018, 18(1).
Kojima K, Oshita M, Nanjo Y, Kasai K, Tozawa Y, Hayashi H, Nishiyama Y: Oxidation of elongation factor G inhibits the synthesis of the D1 protein of photosystem II. Mol Microbiol 2007, 65(4):936-947.
Xu PL, Guo YK, Bai JG, Shang L, Wang XJ: Effects of long-term chilling on ultrastructure and antioxidant activity in leaves of two cucumber cultivars under low light. Physiol Plant 2008, 132(4):467-478.
Mori IC, Schroeder JI: Reactive oxygen species activation of plant Ca2+ channels. A signaling mechanism in polar growth, hormone transduction, stress signaling, and hypothetically mechanotransduction. Plant Physiol 2004, 135(2):702-708.
Liu J, Zhang C, Wei C, Liu X, Wang M, Yu F, Xie Q, Tu J: The RING Finger Ubiquitin E3 Ligase OsHTAS Enhances Heat Tolerance by Promoting H2O2-Induced Stomatal Closure in Rice. Plant Physiol 2016, 170(1):429-443.
Ramachandra Reddy A, Chaitanya KV, Vivekanandan M: Drought-induced responses of photosynthesis and antioxidant metabolism in higher plants. J Plant Physiol 2004, 161(11):1189-1202.
Wang Z, Wang F, Hong Y, Huang J, Shi H, Zhu JK: Two Chloroplast Proteins Suppress Drought Resistance by Affecting ROS Production in Guard Cells. Plant Physiol 2016, 172(4):2491-2503.
Lotfi N, Soleimani A, Vahdati K, Çakmakçı R: Comprehensive biochemical insights into the seed germination of walnut under drought stress. Scientia Horticulturae 2019, 250:329-343.
Lotfi N, Vahdati K, Hassani D, Kholdebarin B, Amiri R: Peroxidase, Guaiacol Peroxidase and Ascorbate Peroxidase Activity Accumulation in Leaves and Roots of Walnut Trees in Response to Drought Stress. Acta Horticulturae 2010(861):309-316.
Gan P, Liu F, Li R, Wang S, Luo J: Chloroplasts- Beyond Energy Capture and Carbon Fixation: Tuning of Photosynthesis in Response to Chilling Stress. Int J Mol Sci 2019, 20(20):5046.
Zong W, Tang N, Yang J, Peng L, Ma S, Xu Y, Li G, Xiong L: Feedback Regulation of ABA Signaling and Biosynthesis by a bZIP Transcription Factor Targets Drought-Resistance-Related Genes. Plant Physiol 2016, 171(4):2810-2825.
Bailey-Serres J, Parker JE, Ainsworth EA, Oldroyd GED, Schroeder JI: Genetic strategies for improving crop yields. Nature 2019, 575(7781):109-118.
Liu C, Mao B, Ou S, Wang W, Liu L, Wu Y, Chu C, Wang X: OsbZIP71, a bZIP transcription factor, confers salinity and drought tolerance in rice. Plant Mol Biol 2014, 84(1-2):19-36.
Okamoto M, Peterson FC, Defries A, Park SY, Endo A, Nambara E, Volkman BF, Cutler SR: Activation of dimeric ABA receptors elicits guard cell closure, ABA-regulated gene expression, and drought tolerance. Proc Natl Acad Sci U S A 2013, 110(29):12132-12137.
Shim JS, Oh N, Chung PJ, Kim YS, Choi YD, Kim JK: Overexpression of OsNAC14 Improves Drought Tolerance in Rice. Front Plant Sci 2018, 9:310.
Tang N, Ma S, Zong W, Yang N, Lv Y, Yan C, Guo Z, Li J, Li X, Xiang Y et al: MODD Mediates Deactivation and Degradation of OsbZIP46 to Negatively Regulate ABA Signaling and Drought Resistance in Rice. Plant Cell 2016, 28(9):2161-2177.
Wang P, Zhao Y, Li Z, Hsu CC, Liu X, Fu L, Hou YJ, Du Y, Xie S, Zhang C et al: Reciprocal Regulation of the TOR Kinase and ABA Receptor Balances Plant Growth and Stress Response. Mol Cell 2018, 69(1):100-112 e106.
Xiao J, Cheng H, Li X, Xiao J, Xu C, Wang S: Rice WRKY13 regulates cross talk between abiotic and biotic stress signaling pathways by selective binding to different cis-elements. Plant Physiol 2013, 163(4):1868-1882.
Yokotani N, Sato Y, Tanabe S, Chujo T, Shimizu T, Okada K, Yamane H, Shimono M, Sugano S, Takatsuji H et al: WRKY76 is a rice transcriptional repressor playing opposite roles in blast disease resistance and cold stress tolerance. J Exp Bot 2013, 64(16):5085-5097.
Ross CA, Liu Y, Shen QJ: The WRKY Gene Family in Rice (Oryza sativa). J Iintegr Plant Biol 2007, 49(6):827-842.
Ross CA, Liu Y, Shen QJ: The WRKY Gene Family in Rice (Oryza sativa). Journal of Integrative Plant Biology 2007, 49(6):827-842.
Group RWW: Nomenclature report on rice WRKY's - Conflict regarding gene names and its solution. Rice (N Y) 2012, 5(1):3.
Pandey SP, Somssich IE: The role of WRKY transcription factors in plant immunity. Plant Physiol 2009, 150(4):1648-1655.
Takeuchi K, Gyohda A, Tominaga M, Kawakatsu M, Hatakeyama A, Ishii N, Shimaya K, Nishimura T, Riemann M, Nick P et al: RSOsPR10 expression in response to environmental stresses is regulated antagonistically by jasmonate/ethylene and salicylic acid signaling pathways in rice roots. Plant Cell Physiol 2011, 52(9):1686-1696.
Liu J-J, Ekramoddoullah AKM: The family 10 of plant pathogenesis-related proteins: Their structure, regulation, and function in response to biotic and abiotic stresses. Physiol Mol Plant P 2006, 68(1-3):3-13.
Agarwal P, Agarwal PK: Pathogenesis related-10 proteins are small, structurally similar but with diverse role in stress signaling. Mol Biol Rep 2014, 41(2):599-611.
Varoquaux N, Cole B, Gao C, Pierroz G, Baker CR, Patel D, Madera M, Jeffers T, Hollingsworth J, Sievert J et al: Transcriptomic analysis of field-droughted sorghum from seedling to maturity reveals biotic and metabolic responses. Proc Natl Acad Sci U S A 2019, 116( 52):27124-27132.
Zhang Z-F, Li Y-Y, Xiao B-Z: Comparative transcriptome analysis highlights the crucial roles of photosynthetic system in drought stress adaptation in upland rice. Scientific Reports 2016, 6(1):19349.
Huang L, Zhang F, Zhang F, Wang W, Zhou Y, Fu B, Li Z: Comparative transcriptome sequencing of tolerant rice introgression line and its parents in response to drought stress. BMC Genomics 2014, 15(1):1026.
Kim TW, Youn JH, Park TK, Kim EJ, Park CH, Wang Z, Kim SK, Kim TW: OST1 Activation by the Brassinosteroid-Regulated Kinase CDG1-Like 1 in Stomatal Closure. Plant Cell 2018, 30(8):1848‐1863.
Borah P, Sharma E, Kaur A, Chandel G, Mohapatra T, Kapoor S, Khurana JP: Analysis of drought-responsive signalling network in two contrasting rice cultivars using transcriptome-based approach. Sci Rep 2017, 7:42131.
Okay S, Derelli E, Unver T: Transcriptome-wide identification of bread wheat WRKY transcription factors in response to drought stress. Mol Genet Genomics 2014, 289(5):765-781.
Langfelder P, Horvath S: WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics 2008, 9:559.
Zhang B, Horvath S: A general framework for weighted gene co-expression network analysis. Stat Appl Genet Mol Biol 2005, 4(1):17.
Zhou XG, Huang XL, Liang SY, Tang SM, Wu SK, Huang TT, Mo ZN, Wang QY: Identifying miRNA and gene modules of colon cancer associated with pathological stage by weighted gene co-expression network analysis. Onco Targets Ther 2018, 11:2815-2830.
Clarke C, Madden SF, Doolan P, Aherne ST, Joyce H, O'Driscoll L, Gallagher WM, Hennessy BT, Moriarty M, Crown J et al: Correlating transcriptional networks to breast cancer survival: a large-scale coexpression analysis. Carcinogenesis 2013, 34(10):2300-2308.
Deng T, Liang A, Liang S, Ma X, Lu X, Duan A, Pang C, Hua G, Liu S, Campanile G et al: Integrative Analysis of Transcriptome and GWAS Data to Identify the Hub Genes Associated With Milk Yield Trait in Buffalo. Front Genet 2019, 10:36.
Du J, Wang S, He C, Zhou B, Ruan YL, Shou H: Identification of regulatory networks and hub genes controlling soybean seed set and size using RNA sequencing analysis. J Exp Bot 2017, 68(8):1955-1972.
Greenham K, Mockler TC, Weinig C, Guadagno CR, Ewers BE, Gehan MA, McClung CR: Temporal network analysis identifies early physiological and transcriptomic indicators of mild drought in Brassica rapa. eLIFE 2017, 6:e29655.
Shaik R, Ramakrishna W: Genes and co-expression modules common to drought and bacterial stress responses in Arabidopsis and rice. PLoS One 2013, 8(10):e77261.
Zhang X, Feng H, Li Z, Li D, Liu S, Huang H, Li M: Application of weighted gene co-expression network analysis to identify key modules and hub genes in oral squamous cell carcinoma tumorigenesis. Onco Targets Ther 2018, 11:6001-6021.
FOYER CH, NOCTOR G: Oxygen processing in photosynthesis: regulation and signalling. New Phytol 2000, 146(3):359-388.
Robinson JM, Bunce JA: Influence of Drought‐Induced Water Stress on Soybean and Spinach Leaf Ascorbate‐Dehydroascorbate Level and Redox Status. Int J Plant Sci 2000, 161(2):271-279.
Sah SK, Reddy KR, Li J: Abscisic Acid and Abiotic Stress Tolerance in Crop Plants. Front Plant Sci 2016, 7:571.
Ma X, Xia H, Liu Y, Wei H, Zheng X, Song C, Chen L, Liu H, Luo L: Transcriptomic and Metabolomic Studies Disclose Key Metabolism Pathways Contributing to Well-maintained Photosynthesis under the Drought and the Consequent Drought-Tolerance in Rice. Front Plant Sci 2016, 7(1886):1886.
Chai M, Bellizzi M, Wan C, Cui Z, Li Y, Wang G-L: The NAC transcription factor OsSWN1 regulates secondary cell wall development in Oryza sativa. J Plant Biol 2015, 58(1):44-51.
Davidson RM, Manosalva PM, Snelling J, Bruce M, Leung H, Leach JE: Rice Germin-Like Proteins: Allelic Diversity and Relationships to Early Stress Responses. Rice 2010, 3(1):43-55.
Zhang L, Gu L, Ringler P, Smith S, Rushton PJ, Shen QJ: Three WRKY transcription factors additively repress abscisic acid and gibberellin signaling in aleurone cells. Plant Sci 2015, 236:214-222.
Takeuchi K, Hasegawa H, Gyohda A, Komatsu S, Okamoto T, Okada K, Terakawa T, Koshiba T: Overexpression of RSOsPR10, a root-specific rice PR10 gene, confers tolerance against drought stress in rice and drought and salt stresses in bentgrass. Plant Cell, Tissue and Organ Culture (PCTOC) 2016, 127(1):35-46.
Li W, Zhu Z, Chern M, Yin J, Yang C, Ran L, Cheng M, He M, Wang K, Wang J et al: A Natural Allele of a Transcription Factor in Rice Confers Broad-Spectrum Blast Resistance. Cell 2017, 170(1):114-126 e115.
Zhang H, Liu Y, Wen F, Yao D, Wang L, Guo J, Ni L, Zhang A, Tan M, Jiang M: A novel rice C2H2-type zinc finger protein, ZFP36, is a key player involved in abscisic acid-induced antioxidant defence and oxidative stress tolerance in rice. J Exp Bot 2014, 65(20):5795-5809.
Choudhury FK, Rivero RM, Blumwald E, Mittler R: Reactive oxygen species, abiotic stress and stress combination. Plant J 2017, 90(5):856-867.
Mittler R: Oxidative stress, antioxidants and stress tolerance. TRENDSin Plant Science 2002, 7(9):405-410.
Mittler R: Abiotic stress, the field environment and stress combination. Trends Plant Sci 2006, 11(1):15-19.
Cai Z, Jia P, Zhang J, Gan P, Shao Q, Jin G, Wang L, Jin J, Yang J, Luo J: Genetic analysis and fine mapping of a qualitative trait locus wpb1 for albino panicle branches in rice. PLoS One 2019, 14(9):e0223228.
Li C, Ng CKY, Fan L-M: MYB transcription factors, active players in abiotic stress signaling. Environ Exp Bot 2015, 114:80-91.
Wang N, Zhang W, Qin M, Li S, Qiao M, Liu Z, Xiang F: Drought Tolerance Conferred in Soybean (Glycine max. L) by GmMYB84, a Novel R2R3-MYB Transcription Factor. Plant Cell Physiol 2017, 58(10):1764-1776.
Naser L, Kourosh V, Bahman K, Reza A: Soluble sugars and proline accumulation play a role as effective indices for drought tolerance screening in Persian walnut (Juglans regiaL.) during germination. Fruits 2010, 65(2):97-112.
Bagnaresi P, Biselli C, Orrù L, Urso S, Crispino L, Abbruscato P, Piffanelli P, Lupotto E, Cattivelli L, Valè G: Comparative transcriptome profiling of the early response to Magnaporthe oryzae in durable resistant vs susceptible rice (Oryza sativa L.) genotypes. PLoS One 2012, 7(12):e51609.
He F, Machemer-Noonan K, Golfier P, Unda F, Dechert J, Zhang W, Hoffmann N, Samuels L, Mansfield SD, Rausch T et al: The in vivo impact of MsLAC1, a Miscanthus laccase isoform, on lignification and lignin composition contrasts with its in vitro substrate preference. BMC Plant Biol 2019, 19(1):552.
Zhou J, Lee C, Zhong R, Ye ZH: MYB58 and MYB63 are transcriptional activators of the lignin biosynthetic pathway during secondary cell wall formation in Arabidopsis. Plant Cell 2009, 21(1):248-266.
Cho HY, Lee C, Hwang SG, Park YC, Lim HL, Jang CS: Overexpression of the OsChI1 gene, encoding a putative laccase precursor, increases tolerance to drought and salinity stress in transgenic Arabidopsis. Gene 2014, 552(1):98-105.
Walter W, Sanchez-Cabo F, Ricote M: GOplot: an R package for visually combining expression data with functional analysis. Bioinformatics 2015, 31(17):2912-2914.
Saeed AI, Sharov V, White J, Li J, Liang W, Bhagabati N, Braisted J, Klapa M, Currier T, Thiagarajan M et al: TM4: A Free, Open-Source System for Microarray Data Management and Analysis. BioTechniques 2003, 34(2):374-378.
Lu K, Li T, He J, Chang W, Zhang R, Liu M, Yu M, Fan Y, Ma J, Sun W et al: qPrimerDB: a thermodynamics-based gene-specific qPCR primer database for 147 organisms. Nucleic Acids Res 2018, 46(D1):D1229-D1236.
Livak KJ, Schmittgen TD: Analysis of relative gene expression data using real-time quantitative PCR and the 2(-△△ C(T)) Method. Methods 2001, 25(4):402-408.

Table 1 DEGs detected in PR403 and PY6 under drought stress treatment

Genotypes	Category of DEG	Sampling point
Genotypes	Category of DEG	22	20
	Up-regulated	1135	2383
PR403	Down-regulated	1917	3438
	Total	3052	5821
	Up-regulated	1748	2655
PY6	Down-regulated	1397	2874
	Total	3145	5529
Common	Up-regulated	459	1202
	Down-regulated	683	1982
	Opposite-regulated	477	532
	Total	1619	3716

Table 2 Differentially expressed hub genes identified in the module grey60

Gene_id	PR403_22^a	PR403_20^b	PY6_22	PY6_20	Description	GO enrichment
Os01G0842500	4.69^c	4.40	5.85	8.49	Similar to Laccase (EC 1.10.3.2).	GO:0005507, GO:0009809, GO:0016722, GO:0046274, GO:0048046, GO:0052716, GO:0055114
Os03G0273200	4.65	3.70	6.79	9.00	Similar to Laccase (EC 1.10.3.2).	GO:0005507, GO:0009809, GO:0016722, GO:0046274, GO:0048046, GO:0052716, GO:0055114
Os01G0702700	2.20	0.81	0.41	1.16	Similar to Transcription factor MYB86 (Myb-related protein 86) (AtMYB86) (Myb homolog 4) (AtMyb4).	GO:0003677
Os05G0543600	2.56	1.51	0.79	2.28	Similar to Myb-related protein.	GO:0003677
Os04G0508500	3.65	3.34	0.50	1.18	MYB family transcription factor, putative, expressed.	GO:0003677
Os06G0131700	2.38	1.80	0.98	2.16	Similar to NAM protein.	GO:0003677, GO:0005634, GO:0006355
Os08G0189200	7.90	7.40	3.48	6.64	Germin-like protein 8-3, Disease resistance	GO:0005618, GO:0030145, GO:0033609, GO:0045735, GO:0046564, GO:0048046
Os03G0433200	3.04	2.62	1.94	2.63	SHR transcription factor, Regulation of the number of cortex cell layers in the root, Coordination of stomatal patterning	GO:0003700, GO:0005634, GO:0006355, GO:0008356, GO:0009956, GO:0043565, GO:0045930, GO:0048366, GO:0055072

Notes: ^a represents the samples collected at the 22% of soil moisture content; ^brepresents the samples collected at the 20% of soil moisture content; ^c the value represents the log2 transformation of the fold change of the expression.

Table 3 Differentially expressed hub genes identified in the module red

Gene_id	PR403 _22^a	PR403 _20^b	PY6 _22		PY6 _20	Description	GO enrichment
OS01G0730700	1.71^c	-0.16	3.00	0.15		WRKY transcription factor 14 (WRKY14).	GO:0003700, GO:0006355, GO:0043565
OS05G0474800	3.36	0.38	4.46	0.36		WRKY transcription factor 70.	GO:0003700, GO:0006355, GO:0043565
OS09G0417600	2.76	-1.65	3.82	0.17		WRKY transcription factor 76.	GO:0003700, GO:0006355, GO:0009617, GO:0009620, GO:0009751, GO:0010200, GO:0043565
OS08G0386200	1.95	-0.28	3.46	0.99		WRKY transcription factor 69.	GO:0003700, GO:0006355, GO:0043565
OS01G0185900	0.19	-3.15	2.43	-2.67		Similar to WRKY 1 (Fragment).	GO:0003700, GO:0006355, GO:0043565
OS05G0322900	1.25	-1.98	3.04	-0.17		Similar to WRKY transcription factor 45.	GO:0003700, GO:0009862, GO:0009864, GO:0010200, GO:0043565, GO:0045892, GO:0050832, GO:1900056
OS05G0571200	3.60	-0.37	5.70	2.53		Similar to WRKY transcription factor 19.	GO:0003700, GO:0006355, GO:0043565
OS11G0117600	3.59	-0.18	5.29	3.60		Similar to WRKY transcription factor 50 (Fragment).	GO:0003700, GO:0006355, GO:0043565
OS03G0437200	1.00	-1.14	5.03	0.67		C2H2-type zinc finger protein, Abscisic acid-induced antioxidant defence, Water stress and oxidative stress tolerance	GO:0003700, GO:0005634, GO:0006355, GO:0043565, GO:0044212, GO:0046872
OS10G0377300	5.42	0.66	3.21	-4.34		Similar to Homeobox-leucine zipper protein HOX8.	GO:0003700, GO:0005634, GO:0006355, GO:0043565
OS12G0555200	3.14	1.10	6.30	3.50		Similar to Probenazole-inducible protein PBZ1.	GO:0004864, GO:0004872, GO:0005634, GO:0005737, GO:0006952, GO:0009607, GO:0009738, GO:0010427, GO:0043086, GO:0080163
OS12G0555500	2.03	0.42	4.62	2.85		Pathogen resistance protein PBZ1,C 17 kDa RNase, Disease resistance	GO:0004864, GO:0004872, GO:0005634, GO:0005737, GO:0006952, GO:0009607, GO:0009738, GO:0010427, GO:0043086, GO:0080163
OS12G0555000	3.75	0.52	6.81	4.27		Similar to Probenazole-inducible protein PBZ1.	GO:0004864, GO:0004872, GO:0005634, GO:0005737, GO:0009738, GO:0010427, GO:0043086, GO:0080163
OS12G0555300	3.35	0.35	5.23	2.70		Similar to Pathogenesis-related protein PR-10a.	GO:0006952, GO:0009607
OS02G0282000	3.52	1.30	4.64	2.65		NB-ARC domain containing protein.	GO:0043531
OS01G0134700	2.02	-0.96	3.78	0.29		Calmodulin binding protein-like family protein.	GO:0002229, GO:0003700, GO:0005516, GO:0005634, GO:0006355, GO:0010112, GO:0010224, GO:0043565, GO:0071219, GO:0080142, GO:1900426
OS01G0508100	2.48	-1.35	4.13	0.89		Ferritin/ribonucleotide reductase-like family protein.	GO:0010112
OS02G0121700	3.68	-0.23	4.63	2.82		Terpenoid synthase domain containing protein.	GO:0000287, GO:0009507, GO:0010333, GO:0016114, GO:0034007, GO:0042742, GO:0043693

Table 3 Differentially expressed hub genes identified in module red (continued)

Gene_id	PR403 _22^a		PR403 _20^b		PY6 _22		PY6 _20	Description	GO enrichment
OS02G0570400	3.55^c	-0.11		5.06		2.08		Similar to Ent-kaurene synthase 1A.	GO:0000287, GO:0006952, GO:0008152, GO:0010333, GO:0034277
OS03G0103300	3.42	-1.76		4.33		1.68		Hybrid proline- or glycine-rich protein, Control of low-temperature germinability, Pre-harvest sprouting resistance	GO:0005783, GO:0009506, GO:0009707, GO:0050832
OS04G0226600	1.63	-1.76		3.56		-0.63		Similar to OSIGBa0135L04.3 protein.	GO:0004674, GO:0005524, GO:0005886, GO:0006468, GO:0006952, GO:0009506
OS05G0368000	4.03	-1.95		5.31		0.50		Non-protein coding transcript.	GO:0010112
OS07G0537900	2.97	-0.27		3.49		0.77		Similar to SRK3 gene.	GO:0004674, GO:0005524, GO:0005886, GO:0006468, GO:0006952, GO:0009506, GO:0016021
OS06G0210400	1.63	-2.81		1.79		-1.61		Concanavalin A-like lectin/glucanase, subgroup domain containing protein.	GO:0004674, GO:0005524, GO:0006468, GO:0016021, GO:0030246
OS04G0109100	5.74	-1.52		8.78		2.92		Concanavalin A-like lectin/glucanase domain containing protein.	GO:0004674, GO:0005524, GO:0006468, GO:0030246
OS04G0370900	3.91	-1.01		5.03		0.77		Similar to H0607F01.5 protein.	GO:0004672, GO:0005509, GO:0005515, GO:0005524, GO:0006468, GO:0016021, GO:0030247
OS04G0371000	4.22	-1.80		4.98		0.05		Four-helical cytokine family protein.	GO:0030247
OS08G0514000	4.02	0.95		6.61		1.49		Protein kinase, catalytic domain containing protein.	GO:0004672, GO:0005524, GO:0006468, GO:0016021, GO:0030246
OS03G0802500	1.90	-1.45		3.13		0.42		ATPase, AAA-type, C core domain containing protein.	GO:0005524, GO:0016021

Table 4 Differentially expressed hub genes identified in the module royalblue

Gene_id	PR403 _22^a	PR403 _20^b	PY6 _22	PY6 _20	Description	GO enrichment
Os01g0298400	1.91^c	-0.31	3.15	0.57	Myb transcription factor domain containing protein.	GO:0003677
Os01g0392600	-0.81	-2.14	2.14	-2.02	Conserved hypothetical protein.
Os03g0663600	1.87	0.68	4.31	2.03	Similar to Pathogenesis-related thaumatin-like protein.
Os12g0630200	0.44	-2.47	2.76	-0.64	Thaumatin%2C pathogenesis-related family protein.
Os11g0592200	0.33	-2.22	2.14	-1.30	Similar to Chitin-binding allergen Bra r 2 (Fragments).	GO:0042742, GO:0050832
Os04g0494100	1.63	0.65	4.83	1.19	Similar to Chitinase.	GO:0000272, GO:0004568, GO:0005618, GO:0006032, GO:0006952, GO:0008061, GO:0009611, GO:0009617, GO:0010262, GO:0016998
Os06g0487620	1.97	1.70	2.31	0.98	Conserved hypothetical protein.
Os06g0493100	5.99	0.26	5.93	1.68	Conserved hypothetical protein.	GO:0016021
Os06g0586300	-0.01	1.13	5.08	3.48	Hypothetical protein.
Os06g0588450	-1.00	-1.82	4.64	2.74	Hypothetical conserved gene.
Os06g0591200	4.28	-2.17	8.51	4.08	Conserved hypothetical protein.

Fig.S1.tif
Fig. S1 Diagram showing the sampling points and the corresponding phenotypes of the plant materials.
Fig.S21.tif
Fig. S2 Determination of APX activity, endogenous ABA level, and stomatal status of PR403 and PY6 under drought treatment. (A) APX activity. (B) CAT activity. (C) NBT staining of leaf samples. (D) Endogenous ABA content. (E) Leaf stomatal opening status during drought treatment. The top panel shows three levels of stomatal aperture: completely open, partially open and completely closed. The bottom panel shows the percentages of three levels of stomatal opening in PY6 and PR403 (n = 100 stomata for PR403 and PY6). The different letters at the top of each column in (A), (B), (D), and (E) indicate statistically significant differences based on ANOVA with Tukey’s HSD test (P < 0.05). Scale bars = 0.5 cm in (B) and = 5 µm in (D).
Fig.S3.tif
Fig. S3 Validation of the RNA-seq data by qRT-PCR. The ACTIN gene (Os11g0163100) was used as an endogenous reference for qPCR. Os01g0164600 (A), Os01g0289600 (B), Os02g0115700 (C), Os03g0319400 (D), and Os04g0610400 (E) were selected for qPCR. The sequences of the primers used are shown in Table S1.
Fig.S4.tif
Fig. S4 Cellular component GO terms significantly overrepresented in the DEG GO enrichment analysis. (A) DEG-enriched GO terms at sampling point 22 of PR403. (B) DEG-enriched GO terms at sampling point 22 of PY6. (C) DEG-enriched GO terms at sampling point 20 of PR403. (D) DEG-enriched GO terms at sampling point 20 of PY6. Bubble size is proportional to the number of each GO-term, and the color represents the -log10 (Qvalue).
Fig.S5.tif
Fig. S5 Molecular function GO terms significantly overrepresented in the DEG GO enrichment analysis. (A) DEG-enriched GO terms at sampling point 22 of PR403. (B) DEG-enriched GO terms at sampling point 22 of PY6. (C) DEG-enriched GO terms at sampling point 20 of PR403. (D) DEG-enriched GO terms at sampling point 20 of PY6. Bubble size is proportional to the number of each GO-term, and the color represents the -log10 (Qvalue).
Fig.S6.tif
Fig. S6 Soft threshold power estimation in the WGCNA.
Fig.S7.tif
Fig. S7 Scatter plot of module eigengenes in the modules significantly correlated with H2O2 and MDA. (A) H2O2 accumulation correlated with the modules black, blue, green, and grey60. (B) MDA accumulation correlated with the modules red, brown, and royalblue.
Fig.S8.tif
Fig. S8 Heatmap showing that the expression patterns of the hub genes in the modules correlated with H2O2 and MDA accumulation. (A) Module Black. (B) Module Blue. (C) Module Grey60. (D) Module Green. (E) Module Red. (F) Module Brown. All the data used in the analysis were subjected to log2 transformation.
Fig.S9.tif
Fig. S9 GO terms overrepresented in the GO enrichment analysis of hub genes in the modules correlated with H2O2 accumulation. The GO terms of biological processes, cellular components, and molecular functions were overrepresented in the modules (A) black, (B) blue, (C) grey60, and (D) green, respectively. Bubble size is proportional to the number of each GO-term, and the color represents the -log10 (Qvalue).
Fig.S10.tif
Fig. S10 GO terms overrepresented in the GO enrichment analysis of hub genes in the modules correlated with MDA accumulation. The GO terms of biological processes, cellular components, and molecular functions were overrepresented in the modules (A) red, (B) brown, and (C) royalblue, respectively. Bubble size is proportional to the number of each GO-term, and the color represents the -log10 (Qvalue).
TableS1.docx
TableS2.xlsx
TableS3.xlsx
TableS4.rtf
TableS5.xlsx

Download PDF

Journal Publication

published 20 Oct, 2020

Read the published version in BMC Plant Biology →

Editor assigned by journal
23 Sep, 2020
Submission checks completed at journal
22 Sep, 2020
Editor invited by journal
22 Sep, 2020

You are reading this older preprint version

Read the latest preprint version →

Weighted gene coexpression network analysis-based identification of key modules and hub genes associated with drought sensitivity in rice

Status:

Journal Publication

Version 2

Abstract

Figures

Background

Results

Discussion

Conclusion

Methods

Declarations

Abbreviations

References

Tables

Supplementary Files

Status:

Journal Publication

Version 2