Temporal transcriptomic differences between tolerant and susceptible genotypes contribute to rice drought tolerance

Drought-tolerance ensures a crop to maintain life activities and protect cell from damages under dehydration. It refers to diverse mechanisms temporally activated when the crop adapts to drought. However, knowledge about the temporal dynamics of rice transcriptome under drought is limited. Here, we investigated temporal transcriptomic dynamics in 12 rice genotypes, which varied in drought tolerance (DT), under a naturally occurred drought in fields. The tolerant genotypes possess less differentially expressed genes (DEGs) while they have higher proportions of upregulated DEGs. Tolerant and susceptible genotypes have great differences in temporally activated biological processes (BPs) during the drought period and at the recovery stage based on their DEGs. The DT-featured BPs, which are activated specially (e.g. raffinose, fucose, and trehalose metabolic processes, etc.) or earlier in the tolerant genotypes (e.g. protein and histone deacetylation, protein peptidyl-prolyl isomerization, transcriptional attenuation, ferric iron transport, etc.) shall contribute to DT. Meanwhile, the tolerant genotypes and the susceptible genotypes also present great differences in photosynthesis and cross-talks among phytohormones under drought. A certain transcriptomic tradeoff between DT and productivity is observed. Tolerant genotypes have a better balance between DT and productivity under drought by activating drought-responsive genes appropriately. Twenty hub genes in the gene coexpression network, which are correlated with DT but without potential penalties in productivity, are recommended as good candidates for DT. Findings of this study provide us informative cues about rice temporal transcriptomic dynamics under drought and strengthen our system-level understandings in rice DT.


Background
Water deficiency and drought is one of the most limiting and disastrous factors for plant adaptation and crop production. It is thus urgent for us to breed water-saving and drought-resistant crops [20,21]. Among various mechanisms of drought resistance, drought tolerance (DT) ensures a plant to maintain its normal life activities under drought and protect plant cell away from damages of dehydration [8]. It is compromised by diverse mechanisms (e.g. stomata conductance, osmotic adjustment, and protective molecules), which are temporally activated when a plant adapts to a naturally-occurred drought [8,29]. Briefly, the plant senses and transmits signals of water deficiency in first and then increases its osmolality to maintain turgor pressures by up-regulating various osmolytes [3,8]. Sooner or later, the plant represents various acclimation responses (e.g. closure of stomata, leaf rolling, decrease in photosynthesis, inhibitions of growth and development) to reduce water loss and consumption [6,11,30,43]. It needs to activate protective mechanisms, such as antioxidant enzymes and molecular chaperones, to protect cell from damages of dehydration at middle-later drought [31,44]. Finally, cell death occurs when the degree of dehydration exceeds its tolerance, representing as leaf senescence and/or fertile abortion [27,32]. The temporal physiological responses of a plant during drought should be determined by the regulation of relevant genes, which activates certain biological processes appropriately. Therefore, learning temporal transcriptomic dynamics of a plant adapting to a long-term drought can deepen our understanding in drought tolerance.
Rice is one of the most important cereal crops which provides > 50% stable food for global populations. The rice cultivation costs huge amount of water and is very sensitive to drought [2]. Water-saving and droughtresistance rice is thus required to ensure global food safety [21]. However, drought tolerance in rice is composed by thousands of genes with minor effects [10]. Hundreds of genes have been reported to be involved in plant responses to drought stress, forming a complicated gene network [10]. It is therefore the improvement of drought tolerance by manipulating on a single gene has limited success in the field, although many droughttolerant genes have been reported to have significant effects in the laboratory or small-scale field conditions. A systematic solution in gene engineering or breeding to improve drought tolerance in rice is required. The methodology of system biology, such as transcriptomics, is a powerful tool to explain complex traits. Many comparative transcriptomic studies on rice drought tolerance have been published in the last decade, contributing a lot in our understandings of drought tolerance at the systematic level [4, 5, 14, 18, 23-26, 34, 40, 46]. However, there are still several limitations remaining. First, the former studies preferred to use a single genotype [24,25,40] or a pair of genotypes with contrasting drought tolerance [4,5,23]. This may introduce genotype-specific bias. Second, drought tolerance in a crop requires stable productivity under drought rather than mere survival. In that case, many acclimation responses (e.g. stomata closure, leaf rolling, decrease in photosynthesis, etc.) inappropriately activated under normal conditions may lead to agronomic penalties. As a result, tradeoffs between DT and productivity have been frequently reported at the gene [15,45], genome [39], transcriptome [42], and population scales [22,41]. For the better understanding and utilization of drought-tolerant genes, such potential tradeoffs between drought tolerance and productivity in a crop should be seriously evaluated at the transcriptomic level. Finally, a comparative transcriptomic study can commonly obtain large number of differentially expressed genes (DEGs), most of which are byproducts accompanying with the plant's morphological and physiological changes. In fact, only a few candidate genes from transcriptomic studies have been functionally characterized, according to the large number of candidates proposed by these transcriptomic studies on drought tolerance. How to effectively isolate the cause gene of DT from the large number of drought-responsive genes (DRGs) is still a challenge for researchers.
In this study, twelve rice genotypes (six tolerant and six susceptible genotypes) were used to investigate their transcriptomic dynamics during a long-term progressive drought. The drought treatment is arranged in a modified field, where drought-avoidance by root systems is separated from drought tolerance by the design of shallow soil layers [23]. The involvement of enough genotypes could mitigate genotype-specific bias in generating drought-responsive genes (DRGs). Meanwhile, the temporal differences between tolerant and susceptible genotypes in the activated biological processes could provide valuable information to determine key biological processes and genes that contribute to drought tolerance. In addition, the weighted gene coexpression network analysis (WGCNA) are applied to extract modules associated with drought-tolerant traits and to identify candidate genes with systematic impacts on DT. We aimed to disclose the transcriptomic dynamics in rice adaptation to the long-term drought and explore the cause of drought tolerance at the transcriptomic level. It can provide novel insights into the system-level understanding of drought tolerance. caused dead leaves (ranged from 0.241 to 0.509) and significant reductions in biomass and yield among rice genotypes ( Table 1). The drought-tolerant index based on the relative biomass (DTIB) ranged from 0.261 to 1.058 (Table 1). We could then divide the twelve genotypes into two groups of contrasting DTs based on their DTIBs . Lower values of DTIB (0.352 ± 0.040) in the susceptible group mean these genotypes may have greater acclimation responses. In contrary, higher values of DTIB in the tolerant group (0.862 ± 0.080) indicated well maintained life activities under drought. We also measured the traits relevant to osmotic adjustment and anti-oxidization. The osmotic potentials in susceptible and tolerant groups were almost the same under well-watered conditions (W) ( Table 2). However, the osmotic potential of the tolerant group was significantly higher than that of the susceptible group during the drought period and recovery (R) stages ( Table 2). The content of H 2 O 2 in the tolerant group was significantly lower than that in the susceptible under both W and drought (D) conditions (Table 2). These results indicated osmotic adjustment and anti-oxidant capacities played essential roles in rice DT.

Summary of differentially expressed genes (DEGs) detected in twelve rice genotypes at six time points
The averaged total base for each sample was 4.48*10 9 bases with the uniquely mapped rate of 87.0%. The Q30 rate was as high as 92.6% (Additional file 2: Table S1). There were 49,470 genes annotated to the reference genome MSU7.0 totally detected to be expressed in this study (FPKM> 0.001), ranging from 33,351 to 42,130 among 12 rice genotypes (Additional file 2: Table S2). Among these expressed genes, 14,354 genes were determined as DEGs at least once among all rice genotypes throughout six experimental time points (D 1 to D 5 , and R), covering 99 known drought-tolerant genes (Additional file 1: Figure S2a, Additional file 2: Table S3).
For each rice genotype, the number of DEGs detected during drought period (D 1 to D 5 ) (defined as drought-responsive genes, DRGs) ranged from 3432 to 6660 (Additional file 1: Figure S3a). The averaged ratio of common DRGs shared by any two rice genotypes was 34.9 ± 6.0% (mean ± SD), suggesting a great genotype-specific manner in rice response to drought. Interestingly, the tolerant genotypes and the Table 1 Agronomic traits (mean ± standard deviation) and drought-tolerance index (DTI) evaluated on twelve rice genotypes under well-watered (W) and drought (D) conditions. "*", "**", and "***" indicate significant differences between W and D at p < 0.05, p < 0.01, and p < 0.001 via independent t test. DTIB or DTIY is calculated as P d /P W *(P d /P ad ). P d is referred as biomass/yield in D, P w is referred as biomass/yield in W, and P ad is referred as averaged biomass/yield of total rice genotypes in D  Figure S3a), indicating distinguished transcriptomic bases between the two groups. Meanwhile, the number of DRGs detected in each genotype was negatively correlated with its DTIB (Fig. 1a). It was also noteworthy that the DTIB was positively correlated with the proportion of upregulated DRGs among rice genotypes (Fig. 1b), indicating positive roles played by upregulated DRGs in DT. Based on the GO enrichment, upregulated and downregulated DRGs were relevant to different classifications of biological processes (Additional file 2: Table S4). For example, upregulated DRGs tended to be relevant to "carbohydrate metabolic process", "response to abiotic stimulus", "cellular component organization and biogenesis", and "cell homeostasis", while downregulated DRGs tended to be relevant to "response to external/ extracellular stimulus", "nucleotide and nucleic acid metabolic process", "cell death", "photosynthesis" and "generation of precursor metabolites and energy" (Additional file 2: Table S4).
To avoid positive false of our transcriptomic data, only the gene represented as a DEG at least in two genotypes at a time point were chosen for further analysis, resulting in 8270 DEGs in the final dataset (Additional file 1: Figure S2b). The excluded 6084 genes (accounting for 42.3%) conferred 13 function-characterized droughttolerant genes (Additional file 1: Figure S2b

Cluster analysis of rice genotypes under W or D conditions based on expression of total genes and DRGs
Cluster based on the FPKMs of total expressed genes could not separate samples from D and W/ or D and recovery (R) at all six time points (Additional file 1: Figure  S5), while cluster based on the FPKMs of DEGs could majorly separate D samples from W/R samples along with the progressive drought (Additional file 1: Figure  S6). Moreover, the tolerant genotypes and the susceptible genotypes could be majorly separated by their FPKMs of DRGs under drought (Additional file 1: Figure  S6). It was still necessary to point out that the susceptible and tolerant genotypes could be well separated by the cluster analysis based on their SNPs called from total transcripts (Additional file 1: Figure S7a) or called from transcripts of DEGs (Additional file 1: Figure S7b). These results allowed us to combine transcriptomic data of tolerant or susceptible genotypes together to form a tolerant or a susceptible group.

Differences in DEGs detected in tolerant and susceptible genotypes at six time points
There were 1444-3893 DRGs in the susceptible group during the drought period and 1358-2410 recovery related genes (RRGs) at the recovery stage (Additional file 1: Figure S3b, Fig. 1c). The percentage of upregulated DRGs was gradually decreased along with the progressive drought (Fig. 1d). As expected, the tolerant group contains more upregulated DRGs. Based on the presence (1) and absence (0) of DEGs in the tolerant group and the susceptible group, they could not be clearly separated at the early drought (D1 and D2), while they were separated at the later drought stages (D3, D4, and D5) along the first and second coordinates (Additional file 1: Figure S8a). However, they were not well separated at the R stage based on the presence and absence of RRGs, There were 3355 common DRGs shared by the tolerant group and the susceptible group during drought period, accounting for 57.8% of total DRGs. There were 2476 susceptible-specific and 1171 tolerant-specific DRGs (Additional file 1: Figure S9). The tolerant-specific DRGs tended to be relevant to catabolic process, carbohydrate metabolic process, and cellular component organization and biogenesis (Additional file 1: Figure S10). The susceptible-specific DRGs tended to be relevant with lipid metabolic process, protein metabolic process, nucleotide and nucleic acid metabolic process, and protein modification process (Additional file 1: Figure S10). Meanwhile, we also detected 605 DRGs having significant differences in fold changes between the tolerant group and the susceptible group (defined as the tolerant group and the susceptible group different DRGs, TS-different DRGs). In most cases (84.6%, 918 in 1085), higher absolute log 2 FC values of these TS-different DRGs were observed in the tolerant group. They were relevant to GO classifications of response to abiotic stimulus, response to stress, carbohydrate metabolic process, catabolic process, and secondary metabolic process (Additional file 1: Figure S10). At recovery, we detected 589 susceptible-specific and 920 tolerant-specific RRGs (Additional file 1: Figure S9).
Temporal regulation of DEGs in the tolerant group and the susceptible group during drought and at recovery The time-series regulation of DRGs could form 16 modes (clusters) respectively in the tolerant (Additional  Table S5). The distribution of osmoticcorrelated or H 2 O 2 -correlated DRGs in the 16 modes was not random (Additional file 1: Figure S11, S12). Osmoticcorrelated DRGs in the susceptible group majorly distributed in mode-1, 2, 5, 7, and 13, while these of the tolerant group majorly distributed in mode-1, 2, 3, 6, 7, 8, and 11. H 2 O 2 -correlated DRGs in the susceptible group majorly distributed in mode-1 and 3, while these of the tolerant group majorly distributed in mode-1 and 2.
The regulation of RRGs during the recovery process (W 5 -D 5 -R) formed eight modes (Additional file 1: Figure S13a). There were 1463 RRGs representing similar modes of timeseries regulation between susceptible and tolerant groups, accounting for 66.9% total common RRGs (Additional file 1: Figure S13b). RRGs belonging to the same mode had similar biological functions as revealed by the PCoA based on enriched GOBPs (Additional file 1: Figure S13c). We thought RRGs in cluster-7 play key roles in recovery as they were specifically activated (upregulated) during the recovery process. RRGs of cluster-7 had great differences in enriched GO biological processes between the tolerant group and the susceptible group (Additional file 2: Table S6, Additional file 1: Figure S13c), providing additional evidence for its important roles in drought recovery.

Temporal differences in biological processes between the tolerant group and the susceptible group
There were 356 GOBPs totally enriched (FDR < 0.05) based on DEGs in the tolerant group and the susceptible group at six time points (Additional file 1: Figure S14). During drought period, 78, 80, 66, 29, 56, and 46 GOBPs belonged to Type I, Type II, Type III, Type IV, Type V, and Type VI, respectively (Additional file 1: Figure S14). GOBPs of Type I were commonly enriched by tolerant and susceptible genotypes, suggesting their universal roles in responses to drought. They were majorly relevant to transport, carbohydrate metabolic process, secondary metabolic process, response to endogenous stimulus, response to abiotic stimulus, and response to biotic stimulus. GOBPs of Type II and Type IV could be considered as the DT-featured BPs. They tended to be relevant to nucleotide and nucleic acid metabolic process ( Table  S7). In contrast, GOBPs of Type III and V could be considered as susceptibility-featured BPs. They tended to be relevant to lipid metabolic process (GO:0006629), secondary metabolic process (GO:0019748), response to external stimulus (GO:0009605), response to endogenous stimulus (GO:0009719), and response to biotic stimulus (GO:0009607) (Additional file 2: Table S7). Finally, there were 179 GOBPs enriched at the recovery stage by the tolerant and the susceptible groups, ten among which were recovery-specific, referring to growth (GO: 0040007), reproductive process (GO:0022414), metal ion homeostasis (GO:0055065), developmental process involved in reproduction (GO:0003006), etc. (Additional file 1: Figure S14).
Based on the enriched GOBPs, the tolerant group and the susceptible group could be separated along with the first coordinate in the PCoA, particularly at the later drought (D 3 -D 5 ) (Additional file 1: Figure  S8b). It indicated the tolerant group and the susceptible group had temporal differences in the activated biological processes in responses to drought, which varied at six experimental time points. The temporal differences in activated biological processes could be revealed by the correspondence analysis (Fig. 2). Early differences (D 1 and D 2 ) in biological processes between the tolerant group and the susceptible group mainly referred to GO classifications of nucleic acid metabolic process (GO:0006139), carbohydrate metabolic process (GO:0005975), generation of precursor metabolites and energy (GO:0006091), tropism (GO: 0009606), response to external stimulus (GO: 0009605), response to biotic stimulus (GO:0009607), cell homeostasis (GO:0019725), and biosynthetic process (GO:0009058). Differences at middle drought (D 3 and D 4 ) mainly referred to protein metabolic process (GO:0019538), lipid metabolic process (GO: 0006629), secondary metabolic process (GO:0019748), catabolic process (GO:0009056), response to endogenous stimulus (GO:0009719), and transport (GO: 0006810). Later differences (D 5 ) mainly referred to cellular component organization and biogenesis (GO: 0016043) and cell death (GO:0008219) (Fig. 2). It was noteworthy that regulation of programmed cell death (GO:0043067) and negative regulation of cell death (GO:0060548) were uniquely enriched in the susceptible group in the middle and later drought (Additional file 1: Figure S14). Differences in GOBPs at the recovery stage mainly referred to growth (GO: 0040007), multicellular organismal development (GO: 0007275), signal transduction (GO:0007165), and secondary metabolic process (GO:0019748) (Fig. 2, Additional file 1: Figure S14).

Cross-talks among phytohormones in tolerant and susceptible genotypes
The proportion of DRGs in total genes in the metabolic pathway of eight phytohormones ranged from 34.5 to 52.4% (Table 3). The occupation of group-distinguished genes (susceptible-specific, tolerant-specific, and TS different) ranged from 22.2 to 60.0% (Table 3). The regulation of ABA-relevant genes (80.0%, 8 out of 10) were mainly upregulated during drought period, while GA-(77.3%, 17 out of 22), auxin-(37.1%, 13 out of 35), and ethylene-relevant (48.1%, 26 out of 54) genes were mainly downregulated (Additional file 1: Figure S24). The regulation of other phytohormone-relevant genes varied among time-points and groups (Additional file 1: Figure  S24). We can find that phytohormones relevant to growth and development (auxin, BR, and GA) have higher proportion of group-distinguished genes ( Table 3). In contrast, the widely acknowledged stress-response phytohormone, ABA, had the lowest proportion of group-distinguished genes (20.0%, 2 out of 10) (Table 3). This result indicated that differences in the drought acclimation should contribute to their differences in rice DT. In addition, we observed a better coordination among JA, SA, BR, auxin, cytokinin, and ethylene in the tolerant genotypes, reflected by higher averaged PCC values among genes relevant to these phytohormones (Additional file 1: Figure S25).

Correlations of DRGs with measured DT traits
There were 183 DRGs negatively correlated (Pearson correlation coefficient > 0.4, p < 0.001) with H 2 O 2 content, while 383 DRGs positively correlated with it. Meanwhile, there were 459 DRGs negatively correlated with osmotic potential while 1292 DRGs positively correlated with it. The overlap between H 2 O 2 -correlated and osmolalitycorrelated DRGs was very rare (Additional file 1: Figure  S26a). As expect, absolute values of Log 2 (fold change) of shared H 2 O 2 -correlated and osmolality-correlated DRGs were significantly higher in the tolerant group than those in the susceptible group (Additional file 1: Figure S27). Twenty-three GOBPs were enriched by H 2 O 2 -correlated DRGs (Additional file 2: Table S8), which mainly referred to GO classifications as metabolic process, cellular process, and transport (Additional file 2: Table S9). Among these enriched GOBPs, energy coupled proton transport (down electrochemical gradient) and ATP synthesis coupled proton transport were tolerant-specific (Type II) (Additional file 2: Table S8). Meanwhile, 136 GOBPs were enriched by osmolality-correlated DRGs (Additional file 2: Table S8), which mainly referred to GO classifications of carbohydrate metabolic process, biosynthetic process, nucleotide and nucleic acid metabolic process, etc. (Additional file 2: Table S9). Among these enriched GOBPs, 13 and 14 GOBPs were of type II and Type IV, respectively (Additional file 2: Table S8).
There were 1463 DRGs significantly correlated with the drought-tolerant index (DTI) (|PCC| > 0.576, p < 0.05) while there were 1403 DRGs significantly correlated with the biomass in the well-watered condition. Among 590 DRGs both correlated with DTI and biomass, a large proportion (499, 81.1%) located at the region of III and VII (Fig. 3), indicating their potential contradictory effects on DT and growth. This result indicated a gene, which is beneficial to DT, may potentially bring negative effects to growth. Interestingly, the tolerant- (Fig. 3b) and susceptible-specific (Fig. 3c) DRGs had opposite patterns in regulating DT and growth (Fig. 3b, c). In addition, a considerable proportion of H 2 O 2 -and osmolality-correlated genes were also correlated with biomass-W (Additional file 1: Figure S26b). This indicated DT mechanisms of ROS scavenging and osmotic adjustment can also have impacts on rice growth.

Coexpression network and candidate hub genes for drought-tolerance and drought-recovery
Based on the result of WGCNA, nine biologically significant modules were generated when the cutoff threshold of edges was set at 0.16, which contained 52 to 3330 nodes ( Table 4). The eigengene of the turquoise module was negatively correlated with the osmotic potential, while the eigengenes of the brown and purple modules were positively correlated with the osmotic potential (Table 4). This result indicated potential associations of the three modules with osmotic adjustment. The eigengenes of the tan, turquoise, yellow, and green modules  (Table 4). This result indicated the potential associations of above modules with ROS scavenging. Meanwhile, each module was involved in many DT-featured biological processes based on the result of GO enrichment (Additional file 2: Table S10). For instance, many carbohydrate transporters, such as D-ribose transport (GO:0015752), pentose transport (GO:0015750), D-xylose transport (GO: 0015753), sorbitol transport (GO:0015795), mannitol transport (GO:0015797), and dicarboxylic acid transport (GO:0006835), were enriched by genes in the brown module. The regulation of these transporters shall contribute to the osmotic adjustment (Additional file 2: Table S10). Meanwhile, the involvement in NADP biosynthetic process (GO:0006741) of the yellow module could be the potential explanation for its correlation with the H 2 O 2 content (Additional file 2: Table S10). It is also noteworthy that the upregulated ABA relevant genes were mainly distributed in the blue and brown module (Additional file 2: Table S11), indicating that the association of the two modules with DT following the ABA-dependent manner. Other six phytohormones (ET, SA, auxin, JA, KT, and GA) were mainly distributed in the turquoise module (Additional file 2: Table S11), whose nodes were mainly downregulated (Table 4). It indicated the turquoise module played a role in the acclimation response. The hub genes in each module played crucial roles in regulating DT and/or acclimation responses. Twenty hub nodes were therefore recommended as good candidates for DT improvement. Most of the candidates were correlated with traits of DT and involved in at least one DT-featured biological process (Additional file 2: Table S12). Meanwhile, they had minor penalties on rice growth as they were not oppositely correlated with DTIB and biomass in the wellwatered field (Additional file 2: Table S12). For drought-recovery, we suggested thirteen candidates, which were involved in DT-specific biological processes at the recovery (Additional file 2: Table S13). They belonged to cluster-4 and cluster-7 in the tolerant group (Additional file 2: Table  S13). The candidates had relative higher expressions in callus based on data collected from the Rice Genome Annotation Project (http://rice.plantbiology.msu.edu/index. shtml) (Additional file 2: Table S13).

Tolerant and susceptible rice genotypes have distinguished transcriptomic features
The transcriptomic studies of stress tolerance via comparison between single tolerant and single susceptible rice genotype always recorded great transcriptomic differences. In this study, the mean ratio of common DRGs shared between any two rice genotypes is about one-third, which is consistent with previous studies (33.4% in Moro versus IR64, [5]; 34.5% in N22 versus IR64, [34]; 35.1% in DD versus IR20, [4]). This means a lot of different DRGs between genotypes follow the genotype-specific pattern. To avoid positive and negative false by the single genotype, we determined DRGs from 12 rice genotypes by its presence frequencies. Interestingly, the tolerant genotypes and the susceptible genotypes could be well separated by the presence and absence of total DRGs. Meanwhile, they could be distinguished by the expressions of DRGs during drought. These results indicate the transcriptomic difference between the tolerant group and the susceptible group should contribute to DT.

Temporal differences in biological processes between susceptible and tolerant genotypes contributions to rice drought tolerance
Plant adapts to the naturally occurred drought progressively by various stepwise morphological and physiological responses [8,29]. Although there have been numerous Table 4 Description of nine modules in the weighted gene coexpression network. PCC is the abbreviation of Pearson's correlation coefficient. *. **, and *** indicate the eigengene is significantly correlated with the trait at p < 0.05, p < 0.01, and p < 0.001, respectively genome-wide transcriptomic studies on rice drought tolerance [4,14,23,25,34,45,46], little is known about rice temporal transcriptomic dynamics in responses to a longterm progressive drought in field conditions. At the gene scale, less DRGs were detected in tolerant genotypes or at early stages of drought, indicating a stable transcriptome under drought is associated with DT. Meanwhile, upregulated DRGs should play an important role in DT, as tolerant genotypes have more upregulated DRGs. This is consistent with empirical selections of upregulated DEGs as candidate genes in many previous transcriptomic studies [18,33].
At the scale of biological process, the tolerant genotypes and the susceptible genotypes share many biological processes (53.4%) during the drought period (Type I), reflecting the common acclimation responses or tolerant mechanisms. Meanwhile, we also detect great temporal differences in activated biological processes between tolerant and susceptible genotypes (Type IIT ype VI), which could help us to separate the tolerant mechanism from the passive acclimation response. Logically, these biological processes uniquely activated in tolerant genotypes (defined as Type II) shall contribute to DT. For example, many biological processes of carbohydrate metabolism (fucose, raffinose, and trehalose) are uniquely enriched in middle-latter drought by DRGs of tolerant genotypes (Fig. 4). They have been proven to be important in DT due to their roles in cellwall synthesis, membrane protection, and photoprotection [9,16,28,37]. Similar biological processes, including rhamnogalacturonan metabolic process, tryptophan and aspartate transports, etc. shall also play a role in DT at different periods (Fig. 4). Besides, if a biological process is activated earlier (e.g. hexose biosynthetic process, gluconeogenesis, etc.) and/or lasted longer (e.g. regulation of protein deacetylation, regulation of histone deacetylation, transcriptional attenuation, transcription antitermination, ferric iron import, ferric iron transport, etc.) during drought in the tolerant group, it shall be another cause of the better DT. For example, histone deacetylases play an essential role in stress responses and several genes regulating histone deacetylation (histone deacetylase and acetyltransferases) have been reported to be associated with DT [7,19,47]. However, roles of histone deacetylation played in rice DT have not been fully understood. In this study, we recorded three histone deacetylase genes (LOC_Os02g12380, LOC_ Os05g36920, and LOC_Os05g36930) as DRGs. It could be valuable for testing their functions in DT. In contrast, if a biological process is initialed earlier in the susceptible group and/or then shared by the tolerant group (defined as Type V), it is very likely a drought-degree dependent process. Interestingly, protein ubiquitination, always associated with protein degradation, belongs to this type. Some other biological processes in Type III and V (e.g. regulation of cell death) reflect earlier and severer occurrences of negative impacts by drought on the susceptible genotypes. There are also many susceptibility-featured biological processes at various time points, particularly in metabolic processes of lipid and secondary metabolites (Fig. 4). Genes relevant to these biological processes could be used as molecular markers for drought susceptibility. Finally, the drought recovery also contributes to DT. It is not a simple reverse process of drought adaption as the involvement of recovery-specific DEGs and biological processes. We consider that the quick recovery in photosynthesis, growth, and development makes significant senses in drought recovery as they represent great differences between the tolerant group and the susceptible group (Fig. 4).
In a word, learning temporal transcriptomic dynamics during long-term drought can strengthen our understanding of drought tolerance and mine valuable drought-tolerant genes that refers to key biological processes. Temporal activations of genes and biological processes can further guide applications of drought-tolerant genes by temporally regulating their expressions via advanced gene engineering.

Good balance between tolerance and productivity in drought-tolerant genotypes
Drought tolerance in crops requires relatively higher productivity under drought. However, wild plants evolve various acclimation mechanisms (e.g. stomata closure, stomata closure, leaf rolling, decrease in photosynthesis, etc.) to ensure good survival by inhibitions of growth and development [6,11,30]. As a result, a gene of DT may possess negative impacts (penalties) on growth and development under normal conditions, which can limit its application in breeding. The existence of penalty in growth has been reported in many previous functional studies of drought-tolerant genes, such as OsIAA6 [15], OsABF1 [45], and OsCYP96B4 [35]. It means once we enhance the drought tolerance by over-expressing or inhibiting a drought-tolerant gene, it may cause unwanted yield penalty in normal conditions. We shall thus fully consider the pleiotropic effect from some droughttolerant genes when utilizing it. In this study, we found hundreds of DRGs have contradictory correlations with productivity potential and DT, which confirms the tradeoff between DT and productivity and provides valuable cues to select propriate candidate genes in breeding.
Another valuable observation is that the regulation of a DRG does not always tend to increase drought tolerance, particularly when it has contradictory effects between DT and productivity. This phenomenon is very common as reported for many drought-tolerant genes, such as SNAC2 [12], OsIAA6 [15], OsABF1 [45] etc. Interestingly, the tolerant genotypes and the susceptible genotypes are seemed to have different strategies in drought adaptation. The tolerant genotype rebuilds the balance between DT and growth much better than the susceptible genotype, by activating tolerant-specific DRGs. This leads to its stable productivity (biomass) recorded under drought. In contrast, the regulation of susceptible-specific DRGs tend to increase the biomass under normal conditions rather than enhance DT. The imbalance between DT and growth may accelerate water loss and over-reproduction of ROS, finally leading to the earlier cell death in susceptible genotypes and poorer productivity.

Conclusion
The tolerant genotypes and the susceptible genotypes have distinguished genetic and transcriptomic bases, which underlies their different morphological and physiological responses to drought. The tolerant genotypes have their featured biological processes, which are temporally-activated through the progressive drought. Meanwhile, the tolerant and the susceptible genotypes represent some differences in their photosynthesis and cross-talks among phytohormones under drought. A considerable transcriptomic tradeoff was detected as Fig. 4 Description of drought tolerance-and susceptibility-featured morphological responses, transcriptomic dynamics, and biological processes during drought and at the recovery stage. * indicates significant difference between tolerant and susceptible genotypes; † indicates different temporal patterns between tolerant and susceptible genotypes. C: carbohydrate metabolic process; G: general biological process; L: lipid metabolic process; N: nucleic acid metabolic process; P: protein metabolic process, S: secondary metabolic process; T: transport hundreds of DRGs are oppositely correlated with DT and productivity. The tolerant genotype has a better balance between DT and productivity under drought by regulating the tolerant-specific DRGs. All these transcriptomic differences disclosed in this study contribute to rice drought tolerance. This knowledge on rice temporal transcriptomic dynamics under drought strengthens our system-level understandings in rice drought tolerance and provides us informative cues in finding and utilizing DT genes.

Plant materials
Twelve rice genotypes were selected to investigate their transcriptomic dynamics in responses to a long-term progressive drought, as well as their morphological and physiological performances. These plant materials were collected from International Rice Research Institute and conserved in the seed bank of Shanghai Agrobiological Gene Center (http://seed.sagc.org.cn/). The twelve rice genotypes have great differences in drought tolerance, which makes them a good system to study transcriptomic bases of drought tolerance. Meanwhile, they are important breeding materials for water-saving and drought-resistance rice (WDR) [20,21]. Once we learn the transcriptomic features for drought tolerance in one or more genotypes, we can apply them in WDR breeding as the donors of drought tolerance.

Field experiment design
Rice seedlings of twelve genotypes were grown in the modified field in the drought-resistance facility at Baihe Experimental Station (Shanghai, China) in 2014. The canopy can be closed on rainy days to keep rainfall out of the experimental field and to enable continuous drought. Different rice genotypes may have varied root lengths, making their capacities of accessing water at depth different. To avoid any influences caused by drought-avoidance, we imitated the design of experiments in column limited pots, making the depth of the soil layer in our experimental field only 30 cm. With this design, the developments of roots to the depth were neutralized among these rice genotypes and therefore their differences of capacities to access water at depth could be largely mitigated. Plants were grown in a plot of 10 rows× 10 hills. Three replicates for each plot were designed in both drought-treated (D) and well-watered conditions (W). The interval between hills was 18 cm. The field arrangement was followed with the single factor randomized block design. We made a minor adjustment in their dates of germination and transplanting to ensure they can meet the drought-stress before heading (Additional file 2: Table S6). We started the drought treatment on 16th, July and continued the artificial drought for 38 days. The drought-stressed field was rewatered on the afternoon of 22th, August. We measured the soil water content once every 3-5 days in the drought-stressed field to monitor and categorize the degree of drought stress at the depth of 30 cm at twelve sites distributing evenly. By monitoring the soil-water content, we can determine the time points of sampling. The averaged coefficient of variance (C.V.) of soil water content was 6.96%, reflecting homogenous levels of soil-water content across our field (Additional file 1: Figure S1).

Measurements and analyses of morphological and physiological traits
The osmotic potential was measured to reflect the capacity of osmotic adjustment for each rice genotype via the Vapro™ vapor pressure osmometer (Wescor Model 5600). We measured the content of H 2 O 2 using the assay kits (product#A064, Nanjing Jiancheng Bioengineering Institute, Jiangsu, China) to reflect the redox status of a rice genotype during the drought. They were measured from leaf samples of three replicates (plots). Plant height and number of tillers were measured from five individuals in each plot. These physiological and morphological traits were measured at six time points at different drought stage indicted by the soil-water content and/or developmental stages by the observation (Additional file 1: Figure S1): 24th of July (D1, later tillering), 29th of July (D2, booting), 5th of August (D3, booting), 11th of August (D4, early flowering), 22nd of August (D5, later flowering), and 23rd of August (recovery stage, R). Meanwhile, three leaf samples for each rice genotype were mixed harvested from each plot at these six time points. They were immediately put in liquid nitrogen for RNA-sequencing. The proportion of dead leaves was observed and estimated on the 23rd, August from five individual plants per plot. We measured the seed-set rating, biomass, and grain yield for each rice genotype after the harvest from all three plots with eight plants per plot. As the seed production of a rice genotype under drought was greatly influenced by its heading date, the relative biomass is more suitable for estimating its drought-tolerance. The drought-tolerant index based on biomass (DTIB) or yield (DTIY) for each rice genotype was calculated as: (biomass or yield in D/ that in W) * (biomass or yield in D/ averaged biomass or yield of all genotypes in D). Twelve genotypes could be divided into a tolerant group and a susceptible group by their DTIB. We applied the independent t-test or oneway ANOVA (SNK method) to detect any significant differences in the measured morphological and physiological traits between the tolerant group and the susceptible group in well-watered and drought-treated fields at each time point.

Procedures of RNA-sequencing
Leaf samples of rice genotypes S3, S11, S17, S26, and S28 from three replicates were mixed for RNA extraction and sequencing. Other rice genotypes (S6, S9, S12, S14, S18, S24, and S31) had two or three biological replicates for RNA-sequencing as leaf samples from two or three replicates were sent for RNA extraction independently. We extracted the total RNA using the PureLink® Plant RNA Reagent (Life Technologies). We used the qualified RNA samples for library construction following the specifications of the TruSeq® RNA Sample Preparation v2 Guide (Illumina) and conducted the RNA-sequencing with Illumina Hiseq 2500 in Shanghai Majorbio Biopharm Technology Co., Ltd. (Shanghai, China). We used SeqPrep to strip adaptors and/or merging paired reads with overlap into single reads (https://github.com/jstjohn/SeqPrep) and used Sickle to remove low-quality reads (https://github. com/najoshi/sickle). We then assembled the clean data using the software Cufflinks and mapped them to the reference genome (Nipponbare, msu7.0) with mitochondrial and chloroplast genomes (http://rice.plantbiology.msu. edu/) via Tophat with no more than two base mismatches allowed in the alignment [36]. The general information of the RNA-sequencing data was provided in Table S1. The raw data for each sample had been submitted to the NCBI Sequence Read Archive (SRA) under the genotype number PRJNA609211.

Data analysis
Determination of the differentially expressed gene (DEG) for each rice genotype between well-watered (W) and droughtstressed (D) treatments We determined the gene expression levels with the Fragment Per Kilobase of exon per Million fragments mapped (FPKM) method via the widely applied software Cuffdiff [36]. The quantification of gene expression by RNA-sequencing was well validated in our previous study using same samples of two rice genotypes [23]. Since rice genotypes S6, S9, S12, S14, S18, S24, and S31 had biological replicates, we determined their DEGs via a false discovery rate (FDR) < 0.05 and a logarithm twofold change |log2FC| ≥ 1. Given the mixed nature of the cDNA library of S3, S11, S17, S26, and S28, we thus determined its DEGs with a p-value < 0.05 and |log2FC| ≥ 1. Additionally, rice plant commonly delayed its development when it exposed to drought (Additional file 2: Table S15). Therefore, a part of differences in terms of transcripts between samples from the well-watered and drought-stressed conditions should be a result of delayed development. Development-dependent genes, which had differences between any two adjacent time points under the well-watered condition, were excluded from the final database for each rice genotype. The DEG detected during drought period was defined as the drought-responsive gene (DRG), while the DEG detected at recovery was defined as the recovery related gene (RRG).
A cluster analysis was conducted to determine relationships among twelve rice genotypes based on their activated DEGs (if a gene was determined as a DEG in a rice genotype, it was scored as "1"; if not, it was scored as "0"). Moreover, as DEGs with low frequencies might be occasionally detected and had low probabilities to be associated with drought-tolerance, we thus exclude DEGs determined only once (~44.5%) among 12 rice genotypes at a time point. Cluster analyses were conducted to see whether rice genotypes within a group (tolerant or susceptible) had common transcriptomic features, using FPKMs of total detected genes and DEGs at six time points.
DEGs having differences between the tolerant group and the susceptible group should contribute to droughttolerance. We thus defined the group-specific DEG which was detected only in the tolerant or susceptible group. We also defined the group-different DEG, which possessed significant difference in the fold change between the tolerant group and the susceptible group by independent t test at a time point. Analyses of Gene Ontology (GO) enrichment was applied for groupspecific and group-different DEGs by. Among three categories of GO terms (cell component, molecular function, and biological process), we paid attentions to GO terms of biological processes (GOBP). We conducted analyses of GO enrichment using the software GOATOOLS (https://github.com/tanghaibao/GOatools) [17].

Temporal gene regulations and relevant biological processes in the tolerant group and the susceptible group
Time-series analyses on gene regulations via hierarchical clustering (Euclidean method) were conducted in the tolerant group and the susceptible group, respectively. The drought period (from D1 to D5) and recovery process (W5-D5-R) were separately analyzed. When conducting time-series analyses during drought period, mean fold changes of DEGs within a group (tolerant or susceptible group) were used. Additionally, correlation analysis of time-series regulation for each DEG was conducted between tolerant or susceptible groups. When conducting time-series analyses at recovery stage, eight regulation modes (clusters) were defined by regulations (up-regulation, unchanged, or down-regulation) from "W5 to D5" and "D5 to R". Analyses of GO enrichment (FDR < 0.05) were conducted for each mode of the tolerant group and the susceptible group during the recovery process, respectively. PCoA was also conducted for each mode based on its enriched GOBPs.
To describe temporal biological processes in responses to the long-term drought, analyses of GO enrichment (FDR < 0.05) were conducted based on DEGs at each time points in the tolerant group and the susceptible group, respectively. PCoA was conducted to describe similarities of transcriptome features between the tolerant group and the susceptible group at six time points using their corresponding DEGs and enriched GOBPs. The enriched GOBPs were further functionally categorized in to dozens of classifications by a web tools named "GO Terms Classifications Counter" (http:// www.animalgenome.org/cgi-bin/util/gotreei) [13]. We applied correspondence analysis to investigate differences of GO classifications between the tolerant group and the susceptible group in responses to the drought at various time points.
We further investigated differences of GOBPs between the tolerant group and the susceptible group in six important classifications of biological processes (p < 0.05 as confidence interval): carbohydrate metabolic process (GO:0005975), protein metabolic process (GO:0006464), lipid metabolic process (GO:0006629), secondary metabolic process (GO:0019748), nucleic acid metabolic process (GO:0006139, short for nucleobase, nucleoside, nucleotide and nucleic acid metabolic process), and transport (GO:0006810). Based on temporal patterns of GO enrichments during drought (D1 to D5), terms of GOBP could divided in to six modes: (I) having the same temporal pattern between the tolerant group and the susceptible group; (II) DT-specific (uniquely enriched in tolerant group); (III) susceptibility-specific (uniquely enriched in susceptible group); (IV) initiated earlier and/ or lasted longer in the tolerant group; (V) initiated earlier and/or lasted longer in the susceptible group; (VI) having other different temporal patterns.
As regulation of photosynthesis is a primary acclimation response to drought stress in plant [11] and has been reported to be associated with drought-tolerance in crops [23,46], we investigated temporal expressions of photosynthesis relevant DEGs (genes annotated to GO:0009767, GO:0009772, GO:0015979, and GO:0022900) in the tolerant group and the susceptible group. Cluster analyses were conducted to see whether rice genotypes within a group (tolerant or susceptible) had common transcriptomic features in terms of photosynthesis, using their FPKMs.

The weighted gene co-expression network analysis (WGCN A)
Co-expression gene network modules (highly coexpressed clusters of genes) were inferred by WGCNA. The automatic one-step network construction and module detection method with default settings were used, which include an unsigned type of topological overlap matrix (TOM). To gain biologically significant modules in a weighted gene co-expression network, the cutoff threshold of edges was determined by the method described previously, at where the network density displays a minimal value [1]. Edges above the cutoff indicated significant relationships of co-expression and were thus determined as significant edges. Only nodes with significant edges were retained in the gene co-expression network. Their values of module eigengenes were calculated and used to test their associations with morphological, physiological, and agronomic traits by correlation analyses. A gene with the top 15% connectivity in a module was determined as a hub gene. Analyses of GO enrichment (FDR < 0.05) were also conducted for each module to investigate its potential biological functions.

Correlation analyses of between fold changes/expressions of DRGs and drought-tolerance/agronomic traits
We conducted a correlation analysis between averaged fold changes (mean fold changes from D 1 to D 5 ) of a DRG and drought-tolerances (DTIB) among twelve rice genotypes (|PCC| > 0.576 and p value < 0.05 for significance). As biomass in the well-watered field represented the growth and productivity, we also conducted correlation analysis between the biomass and averaged gene expression (mean of FPKMs from W1 to W5) in the well-watered field (|PCC| > 0.576 and p value < 0.05 for significance) among twelve rice genotypes. Contradictories of coefficients (positively correlated with DTIB but negatively correlated with biomass; or reverse) from two correlation analyses above could reflect tradeoffs between drought-tolerance and growth at the scale of gene expression. We also conducted correlation analyses between gene expression levels and osmotic potential /H 2 O 2 content using all samples during drought.

Supplementary Information
The online version contains supplementary material available at https://doi. org/10.1186/s12864-020-07193-7.  Figure S14. A heatmap describing temporal differences of enriched GO biological processes (Bonferroni corrected p < 0.05) between the tolerant and the susceptible groups at six time points (D1-D5 and R). Supplementary Figure S15. A heatmap describing temporal differences of enriched GO biological processes (p < 0.05) in the GO classification of carbohydrate metabolic process (GO:0005975) between the tolerant and the susceptible groups at six time points (D1-D5 and R). Supplementary Figure S16. A heatmap describing temporal differences of enriched GO biological processes (p < 0.05) in the GO classification of lipid metabolic process (GO:0006629) between the tolerant and the susceptible groups at six time points (D1-D5 and R). Supplementary Figure S17. A heatmap describing temporal differences of enriched GO biological processes (p < 0.05) in the GO classification of protein metabolic process (GO:0019538) between the tolerant and the susceptible groups at six time points (D1-D5 and R). Supplementary Figure S18. A heatmap describing temporal differences of enriched GO biological processes (p < 0.05) in the GO classification of nucleobase, nucleoside, nucleotide, and nucleic acid metabolic process (GO:0006139) between the tolerant and the susceptible groups at six time points (D1-D5 and R). Supplementary Figure S19. A heatmap describing temporal differences of enriched GO biological processes (p < 0.05) in the GO classification of transporter (GO:0006810) between the tolerant and the susceptible groups at six time points (D1-D5 and R). Supplementary Figure S20. A heatmap describing temporal differences of enriched GO biological processes (p < 0.05) in the GO classification of secondary metabolic process (GO:0019748) between the tolerant and the susceptible groups at six time points (D1-D5 and R). Supplementary Figure S21. Cluster analyses of rice genotypes in the drought field (D) at time points D1 to D5 (a-e) and recovery stage (R) (f) based on expressions of photosynthesis-relevant genes. Supplementary Figure  S22. Cluster analyses of rice genotypes in the well-watered field (W) at time points W1 to W5 (a-e) based on expressions of photosynthesisrelevant genes. Supplementary Figure S23. A heatmap of regulations of photosynthesis-relevant DRGs during drought period (D1-D5) and at recovery (R) by mean log 2 (fold changes) in tolerant (T) and susceptible (S) groups. Supplementary Figure S24. A heatmap of regulations of phytohormone-relevant DRGs during drought period (D1-D5) and at recovery (R) by mean log 2 (fold changes) in tolerant (T) and susceptible (S) groups. Supplementary Figure S25. The matrix of mean Pearson's correlation coefficients (PCCs) calculated from DRGs among eight phytohormones. The above matrix is of tolerant genotypes while the below one is of susceptible genotypes. Supplementary Figure S26. Venn diagram of (a) osmolality-and H 2 O 2 -correlated drought-responsive genes and (b) DTIB-and biomass-correlated drought-responsive genes. Supplementary Figure S27. Comparisons of absolute values of Log 2 (Fold change) of osmolality-(a) and H 2 O 2 -correlated (b) drought-responsive genes (DRGs) between the tolerant and the susceptible groups at six time points. Bars indicate standard errors. "*", "**", and "***" indicate significant differences at p < 0.05, p < 0.01, and p < 0.001 by independent t test between the tolerant and the susceptible groups.
Additional file 2: Supplementary Table S1. Basic information of RNAseq for all samples. Supplementary Table S2. Total expressed genes (FPKM> 0.001) detected in twelve rice genotypes. Supplementary Table  S3. Function-studied drought-tolerant genes, their correlation with drought-tolerance index by biomass (DTIB) and biomass measured in well-watered fields. Supplementary Table S4. Number of enriched GOBPs by upregulated and downregulated DRGs in various classifications of biological processes. Number in shade indicates the GO classification preferred by upregulated or downregulated DRGs. Supplementary  Table S5. Correlations of fold changes from common droughtresponsive genes between the tolerant and the susceptible groups during drought. Supplementary Table S6. Number of biological processes in various GO classifications enriched by recovery related genes in different clusters. Supplementary Table S7. Distributions of GO biological processes of different types in various GO classifications. Percentage in shade indicates the GO classification preferred by the type of GO biological processes. Supplementary Table S8. Enrichments of GO biological processes by osmolality-and H2O2-correlated drought responsive genes. Supplementary Table S9. Distribution of biological processes enriched by osmolality-and H2O2-correlated drought responsive genes in various GO classification. Supplementary Table S10. The enrichment of DT-featured biological processes in Gene Ontology for nine modules. Supplementary Table S11. Distribution of genes relevant to photosynthesis and phytohormones in the modules. Supplementary Table S12. Information of candidate genes of drought-tolerance. Supplementary Table S13. Information of candidate genes for drought-recovery.