Weighted Gene Coexpression Network Analysis- Based Identi cation of Key Modules and Hub Genes that Regulate the Drought Stress Response in the Rice Drought-Sensitive Line PY6


 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. The dss-1 locus 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), 4 and 3 modules were strongly correlated with H2O2 and MDA, respectively. Likewise, GO analysis revealed that the photosynthesis-related GO terms were consistently overrepresented in H2O2-correlated modules. Functional annotation of the DEGs in the H2O2 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.We speculated that drought-induced photosynthetic inhibition leads to H2O2 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.

not only by infection from pathogens, such as viruses, bacteria, and fungi, but also in response to abiotic stresses, including drought and salinity [32]. All PR proteins in both monocots and dicotyledons can be classi ed into 17 families [33]. PR3 and PR4 proteins have plant chitinase activity, and the PR5 family of proteins includes permatins, osmotins, zeanatins, and thaumatin-like proteins [34]. PR10 genes belong to a multigene family in various plant species and exert RNase activity and ligand binding activity and are involved in posttranslational modi cation (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 [34].
Given the rapid development of next-generation deep sequencing technology, RNA-seq has been widely used to explore the transcriptomic pro ling 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 [35,36]. 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 [36]. Moreover, phytohormone signaling and Ca 2+ signaling have also been highlighted in transcript pro ling and have been revealed cross-talk among phytohormone signaling pathways [37,38]. 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 [39,40]. 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 [41,42]. 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 [43]. 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 [44][45][46][47][48][49].
Here, we identi ed 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 droughttolerant 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 pro ling 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.

Results
Identi cation of the dss-1 locus for drought-sensitive phenotypes in CSSLs In our previous study, 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 (Fig. 1A). The overlapping introgressed segments de ned by 167 molecular markers covered the whole genome of Lambayeque 1 (unpublished data). 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 identi ed as presenting a signi cant difference in survival rate (20.0%) compared with its recurrent parent PR403 (93.3%), whereas no signi cant 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 signi cant 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 F1 plants exhibited drought resistance, which was similar to that of the parental line PR403. Then, an F2 population was generated to determine the inherited features of dss-1. In a 151-individual F2 population, 111 members showed drought resistance, and 40 showed drought-sensitive traits (resistant:sensitive ≈ 3:1, χ 2 = 0.1788 < χ 2 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 2 O 2 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 [50] and coupled with the considerable potential for increased accumulation of superoxide and hydrogen peroxide in chloroplasts [51], 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 2 O 2 was detected in PY6 from points 40 to 15, whereas the enhanced content could be observed only at point 22 in PR403; therefore, the signi cantly increased accumulation of H 2 O 2 was detected in PY6 at sampling points 22 and 15 compared with that of their respective controls ( Fig. 2A). The marked accumulation of H 2 O 2 in the leaves of PY6 was con rmed by DAB staining (Fig. 2B). Consistently, the activity of POD, an enzyme for H 2 O 2 scavenging, was upregulated in response to the accumulation of H 2 O 2 (Fig. 2C). However, the activity of APX and CAT, another two major enzymes for H 2 O 2 scavenging, was not induced, and no signi cant 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 signi cant differences between the two genotypes with higher values in PY6. However, the content of SOA in the whole period of drought treatment was not signi cantly 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 signi cantly 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 2 O 2 , 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 [52]. 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 signi cant 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 signi cantly greater. However, PY6 and PR403 exhibited similar patterns in the ability to modulate their stomata in response to drought stress, and no signi cant difference in stomatal density was observed for either genotype (Fig. S2E). Therefore, these results suggested that severe dehydration in drought-treated PY6 might 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 pro ling of PY6 under drought stress
To further reveal the impact of dss-1 on drought stress-induced H 2 O 2 accumulation that may lead to the drought sensitivity of PY6, we investigated the transcriptional pro ling of PY6 and PR403 in response to drought stress through RNA-sEq. Principal component analysis (PCA) was initially performed based on all the transcriptional pro ling 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 identi ed 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- To verify the RNA-seq data quality, ve genes were 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  (Fig. 4C). 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 identi cation of drought-induced hub genes
To further investigate the impact of dss-1 on the regulatory network in response to drought stress and to identify the speci c 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 de ned as modules, and the genes in the same modules had high correlation coe cients. A total of 26 modules were identi ed 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 black, blue, and tan modules were negatively correlated with the majority of physiological traits, whereas the royalblue, brown, red, grey60, orange, and green modules 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 signi cant differences in the contents of H 2 O 2 and MDA were observed during drought stress, we focused on the correlations of H 2 O 2 and MDA with gene modules and identi ed hub genes that were strongly correlated with the two drought-related physiological indicators ( Fig. 2A Table S3-1~-7). Hereafter, we designated these hub genes as DEHGs.
Functional enrichment analysis of hub genes correlated with H 2 O 2 accumulation GO enrichment analysis was performed to functionally cluster the hub genes in the modules that were strongly correlated with H 2 O 2 . A total of 115 and 204 DEHGs were identi ed in the black and blue modules, respectively (Table S3- Table S5-1, -2). Interestingly, in agreement with previous reports [53], the expression of the majority of the DEHGs in the black and blue modules was downregulated without marked log 2 FC differences (|log 2 FC PY6 -log 2 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 grey60 module 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 2 FC PY6 -log 2 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-speci c DNA binding transcription factor activity (GO:0003700), and response to freezing (GO:0050826), were signi cantly consistently induced to a certain extent at point 22 and were downregulated at point 20 in PR403, whereas they were signi cantly 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 [54], 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 [55]. A different upregulation pattern of this gene was observed in PY6 and PR403 in response to drought stress ( Fig. 6C; Table 2).  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). Speci cally, OsWRKY70 (Os05g0474800) and OsWRKY76 (Os09g0417600) have been suggested to participate in ABA signaling and the biotic stress response [28,56].
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 downstream of the jasmonic acid pathway and is positively regulated by WRKY TFs in response to drought and high-salt stresses [32,57], 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 2 H 2 -type zinc nger protein, participates in ABA-OsMPK transduction, which results in the accumulation of H 2 O 2 [58,59]. 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 signi cantly 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.  Notes: a represents the samples collected at the 22% of soil moisture content; b represents the samples collected at the 20% of soil moisture content; c the value represents the log2 transformation of the fold change of the expression.
In total, 18 and 13 DEHGs were identi ed in the brown and royalblue modules, respectively (Fig. 7B, C; Table 4; Table S4). GO analysis revealed that the DEHGs in the brown module play different roles in various biological processes, including response to stress ( Table S4). In the royalblue module, only one GO term was signi cantly 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). Notes: a represents the samples collected at the 22% of soil moisture content; b represents the samples collected at the 20% of soil moisture content; c the value represents the log2 transformation of the fold change of the expression.

Discussion
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 pro les of PY6 in response to drought stress via RNA-seq and WGCNA.  [60]. The overaccumulation of ROS in cells causes severe oxidative damage to membranes (lipid peroxidation), proteins, RNA, and DNA molecules [61]. Higher plants have 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,62]. In this study, signi cant overaccumulation of H 2 O 2 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-  [63]. Thus, photosynthesis is the major biological process that occurs in chloroplasts. 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 [19]. According to the GO enrichment analysis, the drought-induced DEGs identi ed 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 black and blue modules 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 [53], the majority of hub genes clustering in the black and blue modules 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 [18].
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.  [65]. In this study, ve differentially expressed genes annotated as transcription factors, including 3 MYBs, 1 NAD, and 1 SHR, were highly clustered in the grey60 module, which was statistically correlated with the accumulation of H 2 O 2 during drought stress. Among the genes, Os01G0702700, Os05G0543600, and Os04G0508500 are MYB family transcription factors. They exhibited elevated expression in both genotypes in response to drought stress (Table 2) plant cells during drought stress, suggesting that the MYB TFs identi ed here were positively correlated with drought resistance (Fig. 2; Fig. 6C; Table 2). 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. Paradoxically, the overexpression of another laccase-encoding gene, OsChi1 (Os01g0827300), in Arabidopsis resulted in increased tolerance to drought and salinity stresses [66]. Thus, further study needs to be performed to explore the function of the laccases identi ed in the study and determine the discrepancy.
Moreover, it has been demonstrated that laccases are regulated by MYB TFs and are involved in lignin biosynthesis in secondary wall formation of Arabidopsis and Miscanthus [67,68]. Additionally, a rice laccase-encoding gene, OsChi1, has been reported to play an important role in the ROS signaling pathway [66], suggesting that, in the present study, MYB TF-regulated laccases with differential expression are involved in the improvement of drought stress tolerance (Fig. 8).
MDA-correlated hub genes that are differentially regulated in response to drought stress The majority of DEHGs clustering in the modules that were strongly correlated with MDA have been shown to be involved in plant defense signaling pathways. These included WRKY family TFs and PR family proteins (Table 3- (Table 3). Similarly, four members of the PR10 family encoding genes identi ed in the red module and three other families of PR proteins (PR3, PR4, and PR5) identi ed in the royalblue module exhibited similar altered expression patterns in response to drought stress. These 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 red module, encodes a group IIa WRKY TF in rice. It has been shown that overexpression of OsWRKY76 in rice plants suppresses a speci c 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 [28]. Moreover, OsWRKY13 was suggested to regulate the antagonistic cross-talk between drought and disease resistance pathways via the repression of SNAC1 and WRKY45-1 [27]. For the group of identi ed PR-encoding genes, OS12G0555000 (RSOsPR10), a gene encoding the rice root-speci c pathogenesisrelated protein PBZ1, was induced by drought and salinity stresses, as well as jasmonic acid; however, it was strongly inhibited by salicylic acid [32,57]. Moreover, the expression of RSOsPR10 was postulated to be repressed by SA-induced OSTGAs or OsWRKYs [32]. Thus, the WRKY TF-and PR family proteinencoded genes identi ed here may act as convergent nodes and may be responsible for the cross-talk between abiotic and biotic stresses. In addition, ZFP36 (OS03G0437200), a gene encoding a C 2 H 2 zinc nger 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 [58,59]. 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 After all of the data were combined 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 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 identi ed 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.

Conclusion
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 2 O 2 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. Speci cally, 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.

Methods
Plant materials, growth, and treatment conditions Rice seedlings were grown according to previously described methods [17], with minor modi cations. 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 [20].
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 overrepresented GO terms were visualized using ggplot2 and the GOplot R package [69]. A heatmap was generated by MEV v4.9 [70].

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 [41]. For the detection of superoxide, leaf pieces were vacuum-in ltrated 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 [17]. Ethics approval and consent to participate Not applicable.

Consent for publication
Not applicable.       Possible role of dss-1 in reprogramming the transcriptional pro les of rice in response to drought stress.
When rice is exposed to severe drought stress, the inhibition of photosynthesis leads to the overaccumulation of toxic ROS, including H2O2. ABA and H2O2 signals are integrated and perceived and trigger the activation of responsive signaling pathways. First, the transcription factors, MYBs, WRKYs, ZFP36, NACs, etc., was differentially reprogrammed. The downstream responsive genes, including a number of PR family proteins and laccases, are then regulated to respond to drought stress, and ameliorating the drought-induced damage of rice. The introgressed dss-1 may exert its effect on the hub genes and differentially alter their expression, therefore leading to the drought-resistant and -sensitive variation of rice.