Transcriptomic Analysis Reveals Regulatory Networks for Osmotic Water Stress and Rewatering Response in the Leaves of Ginkgo biloba

: To elucidate the transcriptomic regulation mechanisms that underlie the response of Ginkgo biloba to dehydration and rehydration, we used ginkgo saplings exposed to osmotically driven water stress and subsequent rewatering. When compared with a control group, 137, 1453, 1148, and 679 genes were differentially expressed in ginkgo leaves responding to 2, 6, 12, and 24 h of water deﬁcit, and 796 and 1530 genes were differentially expressed responding to 24 and 48 h of rewatering. Upregulated genes participated in the biosynthesis of abscisic acid, eliminating reactive oxygen species (ROS), and biosynthesis of ﬂavonoids and bilobalide, and downregulated genes were involved in water transport and cell wall enlargement in water stress-treated ginkgo leaves. Under rehydration conditions, the genes associated with water transport and cell wall enlargement were upregulated, and the genes that participated in eliminating ROS and the biosynthesis of ﬂavonoids and bilobalide were downregulated in the leaves of G. biloba . Furthermore, the weighted gene coexpression networks were established and correlated with distinct water stress and rewatering time-point samples. Hub genes that act as key players in the networks were identiﬁed. Overall, these results indicate that the gene coexpression networks play essential roles in the transcriptional reconﬁguration of ginkgo leaves in response to water stress and rewatering.


Introduction
Drought is one of the most critical environmental constraints that limit trees' growth and productivity [1]. Current climate models predict significant increases in the frequency, duration, and intensity of droughts because of global warming in the future [2]. Moreover, recent high-profile evidence has documented that extreme droughts can kill millions of trees within a short timeframe [3]. Therefore, it is urgent to better understand the molecular mechanisms underpinning the response of trees to drought.
The impact of drought on molecular processes has been well studied in herbaceous plants and model woody species such as poplar [4,5]. In herbaceous plants, response to water deficit begins with signal perception, in which histidine kinase 1 (HK1) can function as a drought sensor in Arabidopsis [6]. As soon as an extracellular dehydration stimulus is perceived, second messenger molecules such as ROS are immediately generated [7]. ROS not only serves as an alarm signal but also can unrestrictedly oxidize many cellular components and even lead to hypersensitive cell death [8]. Plant cells scavenge excess ROS by evoking at least one antioxidant defense system such as glutathione S-transferase (GST) [9]. Plants also induce the level of heat shock proteins (HSPs) and dehydrins (members of the late embryogenesis abundant proteins' LEAs family) to stabilize membranes and proteins in the process of denaturation [8]. Indeed, several genes encoding HSPs, LEAs, and enzymes involved in eliminating ROS are upregulated by dehydration in various plants [5,10]. It is noteworthy that phytohormones such as abscisic acid (ABA) can also contribute to the orchestration of responses to water deficit [11]. Endogenous ABA rapidly accumulates upon dehydration, stimulating a signaling cascade that leads to downstream responses such as stomatal closure [12]. In leaves, ABA-induced stomatal closure is efficient in reducing water loss and plant growth because of a reduction in transpiration and inhibition of CO 2 uptake for photosynthesis [13].
Plants can also regulate the concentrations and types of secondary metabolites to overcome drought stress [14]. For instance, the concentrations of flavonoids have been found to be stimulated in response to water deficit in the leaves of willow and sea buckthorn [15,16]. A decrease in terpenoids has been documented in the leaves of Eucalyptus when dehydrated [17], but opposite results were found in Quercus ilex and Pinus halepensis under aridity stress [18]. These results suggest that the effect of drought on the amounts of secondary metabolites is complex and differ between tree species. Thus, more studies need to be performed to explore the impact of water deficit on secondary metabolites in various woody plants.
In natural conditions, drought stress is often intermittent, and the growth and productivity of plants are associated with the recovery ability after rehydration [19]. Therefore, the rewatering process is a crucial phenomenon to be considered [20]. Rewatering can offset physiological damage caused by drought stress to a certain extent in herbaceous plants [21]: for instance, the concentrations of malondialdehyde (MDA), H 2 O 2 , superoxide anion, and superoxide radicals induced by drought recover during subsequent rehydration in the leaves of tobacco (Nicotiana tabacum), alfalfa (Medicago truncatula), and purslane (Portulaca oleracea) [21][22][23].
Substantial progress has been made in understanding the response of model trees to dehydration and rewatering at the transcriptomic levels [24,25]. For instance, two recent studies have explored transcriptional reconfiguration of poplar in response to water deficit and rehydration [24,25]. One study has examined global transcriptome in the shoot apical meristem of poplar during drought and rewatering treatment [25]. The droughtresponsive transcriptomic network in the leaves of Populus nigra has been reported to be correlated with response to osmotic stress, response to stimulus, and response to ABA, and rehydration-responsive transcriptomic network to be associated with the biosynthesis of tetrapyrroles and porphyrin, and cellular component biogenesis [24]. Currently, little is known about the physiological response and transcriptional reconfiguration in the leaves of non-model trees in response to drought and rewatering.
Ginkgo biloba is a Chinese native tree species, and has been introduced to many countries owing to its vital ecological and economic value [26,27]. For instance, ginkgo leaves are widely used as a Chinese herb and feed additives due to their valuable flavonoids and terpene lactones content [28,29]. Additionally, G. biloba reveals a high level of ambient adaptability and therefore possesses a very long lifespan [30]. Nevertheless, drought studies in ginkgo have mainly concentrated on physiological changes such as the antioxidant system, secondary metabolites, and phytohormones [31,32]. Indeed, only one study has reported gene expression in the leaves of G. biloba in response to water deficit [32]. Thus, it is critical for us to illustrate the transcriptomic regulation in the leaves of G. biloba in response to drought and rewatering. Here, we exposed ginkgo saplings to polyethylene glycol (PEG)-6000-induced osmotic water stress and rewatering. The objectives of this study were to characterize genome-wide transcriptional reconfiguration stirring in ginkgo leaves in response to water stress and rewatering, as well as to identify key genes in these responses.

Plant Materials and PEG Treatment
One-year-old saplings of G. biloba were cultivated hydroponically with Hoagland solution in a plant growth chamber (day/night temperate: 25/20 • C; light per day, 16 h; light intensity, 200 µmol m −2 s −1 ; relative humidity: 60%) for two weeks. The cultivation solution for these saplings was refreshed every other day. Subsequently, 42 healthy plants that exhibited similar growth performance were selected for water stress treatment. PEG-6000 (20%, w/v) solution, which can quickly induce water stress, was applied in the cultivation solution to simulate drought stress. Mature leaves (leaf plastochron index = 9-12) were collected chronologically at 0 h (control, CK), 2 h (water stress 1, D1), 6 h (D2), 12 h (D3), and 24 h (D4), respectively. After 24 h osmotically driven water stress, the hydroponics solution was renewed to simulate rewatering. Then, the mature leaves were collected at 24 h (Rewatering 1, R1) and 48 h (R2).
Six plants were assigned to each sampling point. To avoid injury effects, sampling was explored from distinct sets of plants. The harvested leaves from each plant were wrapped in tin foil and frozen immediately in liquid nitrogen. The frozen samples were milled into fine powder and stored at −80 • C until their further analysis. Within the same sampling point, equal amounts of powder from each of two plants were combined, thus yielding three pooled samples for every sampling point.

RNA Extraction, cDNA Library Preparation, and Sequencing
Total RNA was isolated using a Mini BEST Plant RNA Extraction Kit (Takara, Dalian, China) following the manufacturer's protocols. To remove genomic DNA contamination, total RNA samples were treated with Turbo DNA-free Kit (AM1907, Ambion, CA, USA). The RNA quality was checked using a Bioanalyzer 2100 (Agilent, Santa Clara, CA, USA). The RNA quantity was determined using an RNA 6000 Nano Chip (5067-1511, Agilent, CA, USA). The sequencing libraries were established using a TruSeq™ RNA Sample Preparation v2 protocol (Illumina, San Diego, CA, USA). In brief, poly A-containing mRNA was purified from total RNA, followed by fragmentation and the generation of doublestranded cDNA. After end repair, adapter ligation, and PCR amplification, the libraries were sequenced on a HiSeq 4000 instrument (Illumina), obtaining pair-end 150 bp reads.
To obtain high quality reads, raw reads were processed by removing low-quality reads and trimming adaptor sequences by using in-house Perl scripts. The sequencing data were deposited in the Genome Sequence Archive (GSA, https://ngdc.cncb.ac.cn/gsa/, accessed on 1 December 2021) (accession number CRA004752). Three cDNA libraries per sampling point were generated and sequenced in this study.

Sequence Analysis
The high-quality reads were mapped to the G. biloba reference genome (http://gigadb. org/dataset/100613, accessed on 1 December 2021) using the spliced aligner software (Hisat2, version 2.2.4) [33] with default parameter. The read counts of genes were quantified by using htseq-count software as suggested previously [34]. The differential expression of genes was calculated using the read counts with the DESeq2 package [35]. Additionally, the fragments per kilobase of exon model per million mapped reads (FPKM) values of genes were calculated and used for the weighted gene coexpression network analysis (WGCNA). Resulting p-values of statistics were corrected to the false discovery rate (FDR). Significantly differentially expressed genes (DEGs) were determined with the threshold of FDR < 0.05 and |log 2 (fold change)| ≥ 1.

Annotation, Functional Categorization, Gene Ontology, and Kyoto Encyclopedia of Genes and Genomes Analysis
To annotate the significantly DEGs, the coding sequences of DEGs were retrieved from the G. biloba database. The closest homolog of a G. biloba gene in the genome database of A. thaliana was determined by translated nucleotide BLAST (BLASTX) of the coding sequence hit against the protein sequence database (version Araport11) of A. thaliana. Afterwards, the gene identifiers of A. thaliana were submitted to the MapMan software package for functional analysis as suggested by Luo et al. [36].
Gene ontology (GO) enrichment analysis of significantly DEGs was performed using Blast2GO (version 4.0) [37]. Metascape was used to test the Kyoto encyclopedia of genes and genomes (KEGG) pathway enrichment analysis of significantly DEGs [38]. The p-value threshold for significance of GO and KEGG terms was 0.05.

Weighted Gene Coexpression Network Analysis
The WGCNA of significantly DEGs was performed using the WGCNA package in R as suggested previously [39,40]. Briefly, the significantly DEGs with FPKM greater than 5 were imported into WGCNA to establish coexpression modules using the blockwiseModules function with the power set as 10, mergeCutHeight set as 0.2, and minModuleSize set as 30.
Hub genes of modules that were highly associated with water stress or re-water treatments were determined based on the values of module eigengene (K ME ≥ 0.9). The coexpression networks were visualized using Cytoscape (version 3.7.1).

Validation of RNA Sequencing Data
Total RNA was isolated and purified as mentioned above for RNA sequencing. Nine genes were randomly selected for validation using real-time quantitative PCR (RT-qPCR) as proposed by Lu et al. [41]. Specific primers for selected genes are presented in Table S1. Glyceraldehyde-3-phosphate dehydrogenase (GAPDH) was used as internal control (Table S1).

Transcriptome Changes in Response to Osmotically Driven Water Stress and Rewatering
In total, 21 cDNA libraries from the control group (0 h, CK), water stress (2 h, D1; 6 h, D2; 12 h, D3; 24 h, D4), and rehydration groups (24 h, R1; 48 h, R2) were constructed. After filtering and sequence trimming, nearly 45.8 million to 62.7 million clean bases were obtained per transcriptome library (Table S2). The Q30 was more than 94.8%, and the GC content was between 44.7% and 46.8% (Table S2). More than 88.2% of the high-quality reads could be mapped to the G. biloba genome (Table S2). These results suggested that our RNA-seq data were reliable.
To interpret the major biological processes (BPs) induced by the water stress and rewatering, we conducted GO enrichment analysis at each of the treatment time-points (D1-D4, as well as R1 and R2) ( Figure 2). We found that eight BPs were significantly enriched throughout the studied period of osmotic water stress and rewatering treatment ( Figure 2). These eight common ontologies were: response to a stimulus, response to a chemical, response to stress, response to a temperature stimulus, response to an organic substance, response to a hormone, response to an endogenous stimulus, and response to an abiotic stimulus ( Figure 2). Half of the enriched ontologies belonged to response to a stimulus, implying that this category dominated the common set. Moreover, this dominant category contained not only genes involved in the response to an external environmental stimulus (e.g., response to an abiotic stimulus) but also genes related to the response to an endogenous stimulus, which suggests that the signaling of water stress and rewatering and the plant endogenous development programs have a significant cross-talk. We also found that five GO terms were enriched in the period between D2 and D4, as well as R1 and R2 ( Figure 2). These five common GO terms were: the oxidation-reduction process, the metabolic process, signal transduction, flavonoid biosynthesis, and metabolic processes. Consistently, the KEGG pathway enrichment analysis revealed that flavonoid biosynthesis and monolignol biosynthesis pathways were enriched in D2, D3, D4, R1, and R2 ( Figure S2). These results indicated that flavonoid biosynthesis plays a pivotal role in the response to water stress and rehydration in ginkgo leaves. To interpret the major biological processes (BPs) induced by the water stress and rewatering, we conducted GO enrichment analysis at each of the treatment time-points (D1-D4, as well as R1 and R2) ( Figure 2). We found that eight BPs were significantly enriched throughout the studied period of osmotic water stress and rewatering treatment ( Figure  2). These eight common ontologies were: response to a stimulus, response to a chemical, response to stress, response to a temperature stimulus, response to an organic substance, response to a hormone, response to an endogenous stimulus, and response to an abiotic stimulus ( Figure 2). Half of the enriched ontologies belonged to response to a stimulus, implying that this category dominated the common set. Moreover, this dominant category contained not only genes involved in the response to an external environmental stimulus (e.g., response to an abiotic stimulus) but also genes related to the response to an endogenous stimulus, which suggests that the signaling of water stress and rewatering and the plant endogenous development programs have a significant cross-talk. We also found that five GO terms were enriched in the period between D2 and D4, as well as R1 and R2 ( Figure 2). These five common GO terms were: the oxidation-reduction process, the metabolic process, signal transduction, flavonoid biosynthesis, and metabolic processes. Consistently, the KEGG pathway enrichment analysis revealed that flavonoid biosynthesis and monolignol biosynthesis pathways were enriched in D2, D3, D4, R1, and R2 ( Figure  S2). These results indicated that flavonoid biosynthesis plays a pivotal role in the response to water stress and rehydration in ginkgo leaves.  Table S3.
To further decipher the functions of the water stress-responsive and/or rewateringresponsive DEGs, DEGs that had homologs in Arabidopsis were selected, and these corresponding unique Arabidopsis gene identifiers (AGIs) were assigned to functional categories using Mapman (Table S3). Several water stress-responsive and/or rewatering-responsive genes were associated with the Mapman categories of hormone metabolism, stress, misc, development, transport, and cell wall (Table S3).
In the ginkgo leaves under osmotic water stress treatment, in the hormone metabolism category, the transcript level of the ginkgo homolog to Nine-cis-epoxycarotenoid dioxygenase 9 (NCED9) encoding an NCED, which catalyzes the first step of ABA biosynthesis [42], was elevated in D1 and D3 when compared with that in the control group (Table S3). The transcript level of Cytochrome P450, family 707, subfamily A, polypeptide 1 (CYP707A1) encoding a protein with ABA 8 -hydroxylase activity that is involved in ABA catabolism [43] increased in D2 vs. CK of ginkgo leaves (Table S3). The transcript level of a ginkgo homolog to GA requiring 3 (GA3), which is involved in gibberellins (GAs) biosynthesis [44], was reduced in D1, D2, and D3, whereas the transcript level of Gibberellin 2-oxidase 6 (GA2OX6), encoding a gibberellin 2-oxidase that converts bioactive gibberellins into inactive forms [45], increased in D2 when compared with that in CK of ginkgo leaves (Table S3). The mRNA level of a ginkgo gene homologous to Arabidopsis lipoxygenase 1 (LOX1), encoding a lipoxygenase participates in jasmonate biosynthesis [46], was stimulated in D3 vs. CK, and the transcript level of LOX3 was elevated in D2 vs. CK (Table S3). The transcript level of Ethylene responsive element binding factor 2 (ERF2), which encodes members of the ERF/AP2 TF family and the protein of which contains an AP2 domain [47], was reduced in D2, D3, and D4 of ginkgo leaves compared to the control group (Table S3).
In ginkgo leaves under osmotically driven water stress treatment, in the stress category, the expression levels of Heat shock protein 70 (HSP70) and HSP23.5 were upregulated in D2, and the transcript level of HSP90.5 was induced in D3, while HSP17.4, HSP17.6II, HSP18.2, HSP21, HSP22.0, HSP81.1, and HSP101 were elevated in D4 when compared with those in the CK leaves of G. biloba (Table S3). The mRNA levels of genes encoding heat shock factors (e.g., HSF3 and HSFA4A) that control resistance to environmental stress [48] were upregulated in D2 when compared with those in CK of ginkgo leaves (Table S3). Additionally, the transcript level of a ginkgo homolog to Arabidopsis Dehydration-responsive protein-related was induced in D2, and the mRNA level of Early responsive to dehydration 4 (ERD4) increased in D3 of ginkgo leaves compared to those in CK (Table S3). In the misc category, the mRNA level of a ginkgo gene encoding Glutathione s-transferase tau 16 (GSTU16), which is involved in scavenging ROS [49], was stimulated in D2 vs. CK, and the mRNA levels of GSTU20 and GSTU28 increased in D3 vs. CK (Table S3). Additionally, the expression level of GSTU19 was elevated by both D2 and D3 in the leaves of G. biloba (Table S3). In the development category, the expression level of ginkgo homolog to LEA (late embryogenesis abundant) protein-related increased in D4 compared to CK of ginkgo leaves (Table S3). To further decipher the functions of the water stress-responsive and/or rewateringresponsive DEGs, DEGs that had homologs in Arabidopsis were selected, and these corresponding unique Arabidopsis gene identifiers (AGIs) were assigned to functional catego- In the ginkgo leaves under water stress treatment, in the transport category, a ginkgo homolog to Arabidopsis tonoplast intrinsic protein 2;1 (AtTIP2;1) encoding one of the TIPs related to water transport across the tonoplast [50] was downregulated in D1, D2, and D3 when compared with that in CK of ginkgo leaves (Table S3). Similarly, the transcript level of TIP1;1 decreased in D2 and D3, and the mRNA level of TIP4;1 was reduced in D4 when compared with those in CK leaves of G. biloba (Table S3). The expression level of Plasma membrane intrinsic protein 2;8 (PIP2;8), which is involved in water transport across the plasma membrane [51], decreased in D2, D3, and D4 compared to that in CK of ginkgo leaves (Table S3). The mRNA level of PIP2;8 was also validated by RT-qPCR ( Figure S1, Table S1). In the cell wall category, a ginkgo homolog to Arabidopsis expansin A1 (EXPA1), which causes loosening and extension of plant cell walls [52], was downregulated in D2 and D4 compared to the CK (Table S3). The gene encoding the Expansin-like B3 precursor (EXLB3) was also downregulated in D2 and D3 compared to that in CK of ginkgo leaves (Table S3). Additionally, 12 h and 24 h dehydration led to downregulation of Xyloglucan endotransglucosylase/hydrolase 9 (XTH9), Xyloglucan endotransglycosylase 6 (XTR6), and Endoxyloglucan transferase A4 (EXGT-A4) in the leaves of G. biloba (Table S3). Additionally, the mRNA levels of Gb_12742 and Gb_12743 encoding Abscisic stress-ripening protein 2 (ASR2), as a candidate for drought tolerance in Oryza sativa [53], were stimulated in D2 compared to CK (Table S3).
In ginkgo leaves under the rehydration treatment, in the hormone metabolism category, the mRNA level of NCED4 increased in R1 and R2, and the transcript levels of NCED9 and CYP707A1 were elevated in R2 compared to the control group (Table S3). The expression of Gibberellin 20-oxidase 2 (GA20OX2), which is involved in GA biosynthesis [54], increased in R1 and R2, whereas the transcript level of GA2OX6 was reduced in R1 compared to those in CK (Table S3). In the stress category, the expression levels of three genes that encode HSP70, HSP70-1, and HSP90.7 were reduced in both R1 and R2, and the transcript levels of HSP18.2 and HSP21 decreased in R2 when compared with those in CK (Table S3). In the misc category, when compared with the control group, the transcript levels of Glutathione S-transferase phi 9 (GSTF9), GSTU18, GSTU19, and GSTU28 decreased at the 24 h and 48 h rewatering time points in the leaves of G. biloba (Table S3). In the development category, a ginkgo gene homologous to LEA protein-related was downexpressed in R1 and R2, and the mRNA level of LEA5 was downregulated in R2 when compared with those in CK (Table S3). The transcript levels of two genes (Gb_07622 and Gb_07623) encoding Dehydrin, known as a group of plant LEA proteins [55], were reduced after 24 h and 48 h rehydration (Table S3). In the transport category, a gene that encodes TIP4;1 was upregulated in R2 compared to CK (Table S3). In the cell wall category, the transcript level of XTR6 was reduced in R1 and R2 compared to the CK (Table S3). In contrast, when compared to CK, the mRNA level of XTR2 was induced in 24 h rehydration, and the expression of EXGT-A2, EXGT-A4, and EXLA1 increased after 48 h of rehydration (Table S3). Additionally, the expression level of ASR2 was reduced in R2 when compared with that in the CK leaves of G. biloba (Table S3).

WGCNA Network in Ginkgo Leaves under Osmotic Water Stress and Rewatering
To detect specific genes that play pivotal roles in response to osmotically driven water stress and/or rehydration in ginkgo leaves, WGCNA was performed using the DEGs under each comparison (Figure 3). A total of 17 coexpression gene modules were obtained ( Figure 3A). Moreover, the DEGs in each of these modules varied from 23 to 873 ( Figure 3A). These modules are related to diverse samples ( Figure 3B). In the royal blue module, 42 genes were highly associated with D1 ( Figure 3B, Table S4). In the brown module, 873 genes were highly associated with D2 and D3 ( Figure 3B, Table S5). In the light yellow module, 73 genes were identified as specific to the 24 h dehydration ( Figure 3B, Table S6). In the purple module, 135 genes were tightly associated with R1 ( Figure 3B, Table S7). In the blue module, 698 genes positively correlated with R2 ( Figure 3B, Table S8). For each sample-correlated module, a gene coexpression network was constructed (Figure 4 and S3). Subsequently, GO enrichment analysis was conducted for these five modules ( Figure 5). The term of response to stress was significantly enriched in all these modules ( Figure 5). The ontology of the response to an abiotic stimulus was significantly enriched in the brown, light yellow, and blue modules ( Figure 5). Moreover, the terms of the response to an endogenous stimulus and response to a hormone were enriched in the royal blue, brown, and blue modules ( Figure 5). The ontology of signal transduction was enriched in the royal blue and brown modules, and the GO category of positive regulation of the developmental process was enriched in the blue module ( Figure 5). For each sample-correlated module, a gene coexpression network was constructed (Figures 4 and S3). Subsequently, GO enrichment analysis was conducted for these five modules ( Figure 5). The term of response to stress was significantly enriched in all these modules ( Figure 5). The ontology of the response to an abiotic stimulus was significantly enriched in the brown, light yellow, and blue modules ( Figure 5). Moreover, the terms of the response to an endogenous stimulus and response to a hormone were enriched in the royal blue, brown, and blue modules ( Figure 5). The ontology of signal transduction was enriched in the royal blue and brown modules, and the GO category of positive regulation of the developmental process was enriched in the blue module ( Figure 5).  To further determine the key genes in the royal blue, brown, light yellow, purple, and blue modules, the hub genes that show most connections in each coexpression network as indicated by their high KME values were identified (Figures 4 and S3). In particular, in the coexpression network of the royal blue module, hub genes included Related to AP2.4-A (RAP2.4A) and RAP2.4B ( Figure 4A), which encode transcription factors belonging to the abiotic stress-associated DREB A-6 clade [56]. In the coexpression network of the brown module, Heat shock factor 3 (HSF3), Mitochondrial GRPE 2 (MGE2), which controls the nucleotide-dependent binding of mitochondrial HSP70 [57], Glycine-rich RNAbinding protein 3 (GR-RBP3), involved in RNA transcription or processing during environmental stress [58], and Leaf rust 10 disease-resistance locus receptor-like protein kinase-like 1.2 (LRK10L1.2), which is positively involved in ABA signaling [59], were identified as hub genes ( Figure S3A). In the coexpression of the light yellow module, three genes that encode the HSPs (HSP101, HSP17.8-CI, and HSP17.4-CIII) were recognized as hub genes ( Figure 4B). In the coexpression network of the purple module, hub genes included Disease resistance protein UNI (UNI), which participates in disease resistance through the SA signaling pathway [60], and Constitutive expression of PR genes 5 (CPR5), which modulates SA To further determine the key genes in the royal blue, brown, light yellow, purple, and blue modules, the hub genes that show most connections in each coexpression network as indicated by their high K ME values were identified (Figures 4 and S3). In particular, in the coexpression network of the royal blue module, hub genes included Related to AP2.4-A (RAP2.4A) and RAP2.4B (Figure 4A), which encode transcription factors belonging to the abiotic stress-associated DREB A-6 clade [56]. In the coexpression network of the brown module, Heat shock factor 3 (HSF3), Mitochondrial GRPE 2 (MGE2), which controls the nucleotide-dependent binding of mitochondrial HSP70 [57], Glycine-rich RNA-binding protein 3 (GR-RBP3), involved in RNA transcription or processing during environmental stress [58], and Leaf rust 10 disease-resistance locus receptor-like protein kinase-like 1.2 (LRK10L1.2), which is positively involved in ABA signaling [59], were identified as hub genes ( Figure S3A). In the coexpression of the light yellow module, three genes that encode the HSPs (HSP101, HSP17.8-CI, and HSP17.4-CIII) were recognized as hub genes ( Figure 4B). In the coexpression network of the purple module, hub genes included Disease resistance protein UNI (UNI), which participates in disease resistance through the SA signaling pathway [60], and Constitutive expression of PR genes 5 (CPR5), which modulates SA and the unfolded protein response to balances between plant growth and stress responses [61] ( Figure 4C). In the coexpression network of the blue module, two genes (GA20OX2 and Multiprotein bridging factor 1C, MBF1C), which belong to the hormone metabolism category, and BCL-2-associated athanogene 6 (BAG6), which regulates programmed cell death and stress responses [62], were regarded as hub genes ( Figure S3B). and the unfolded protein response to balances between plant growth and stress responses [61] ( Figure 4C). In the coexpression network of the blue module, two genes (GA20OX2 and Multiprotein bridging factor 1C, MBF1C), which belong to the hormone metabolism category, and BCL-2-associated athanogene 6 (BAG6), which regulates programmed cell death and stress responses [62], were regarded as hub genes ( Figure S3B).

Analysis of Genes Associated with Flavonoid and Bilobalide Biosynthesis in Ginkgo Leaves under Water Stress and Rewatering
The expression of genes that participate in flavonoid and bilobalide biosynthesis was also evaluated ( Figure 6). Most genes that take part in flavonoid biosynthesis were upregulated by osmotically driven water stress, but were downregulated by rewatering in the leaves of G. biloba (Figure 6A,C). For instance, the mRNA levels of Chalcone synthase (CHS) and Flavonol synthase (FLS) were elevated after 2 h, 6 h, 12 h, and 24 h drought, whereas the mRNA levels of these genes decreased after 12 h and 24 h rehydration ( Figure 6A,C). The transcript levels of CHS, Cinnamoyl-CoA reductase 1 (CCR1), and Flavanone 3-hydroxylase (F3H), involved in flavonoid biosynthesis, were also validated by RT-qPCR ( Figure  S1 and Table S1). Similarly, we also found that several genes such as Hydro xymethylglutaryl-CoA synthase 1 (HMGS1), Hydroxy methylglutaryl-CoA reductase 1 (HMGR1), and HMGR2, which participate in bilobalide biosynthesis [63], were upregulated under osmotic water stress, but downregulated under rewatering conditions ( Figure 6B,C). These results indicated that the biosynthesis of flavonoid and bilobalide might be stimulated by water stress but will probably be suppressed by rehydration.

Analysis of Genes Associated with Flavonoid and Bilobalide Biosynthesis in Ginkgo Leaves under Water Stress and Rewatering
The expression of genes that participate in flavonoid and bilobalide biosynthesis was also evaluated ( Figure 6). Most genes that take part in flavonoid biosynthesis were upregulated by osmotically driven water stress, but were downregulated by rewatering in the leaves of G. biloba (Figure 6A,C). For instance, the mRNA levels of Chalcone synthase (CHS) and Flavonol synthase (FLS) were elevated after 2 h, 6 h, 12 h, and 24 h drought, whereas the mRNA levels of these genes decreased after 12 h and 24 h rehydration ( Figure 6A,C). The transcript levels of CHS, Cinnamoyl-CoA reductase 1 (CCR1), and Flavanone 3-hydroxylase (F3H), involved in flavonoid biosynthesis, were also validated by RT-qPCR ( Figure S1 and Table S1). Similarly, we also found that several genes such as Hydro xymethylglutaryl-CoA synthase 1 (HMGS1), Hydroxy methylglutaryl-CoA reductase 1 (HMGR1), and HMGR2, which participate in bilobalide biosynthesis [63], were upregulated under osmotic water stress, but downregulated under rewatering conditions ( Figure 6B,C). These results indicated that the biosynthesis of flavonoid and bilobalide might be stimulated by water stress but will probably be suppressed by rehydration.

DEGs Related to Stress, Development, Transport, Cell Wall, and Hormone Metabolism Might Play Pivotal Roles in Response to Simulated Water Deficit Stress and Rewatering in the Leaves of G. biloba
Drought, one of the most essential stresses, has a negative impact on plant growth and development [10]. Plenty of previous studies have revealed that drought stress brings Drought, one of the most essential stresses, has a negative impact on plant growth and development [10]. Plenty of previous studies have revealed that drought stress brings about excessive ROS in plants [11]. Indeed, elevated concentration of H 2 O 2 under 24 h of PEG6000-simulated drought stress has been reported in ginkgo leaves [32]. Due to their high toxicity, ROS can damage proteins, DNA, and lipids and even trigger cell death [8]. Plants have developed efficient antioxidant defense systems to scavenge ROS [64]. For instance, transcriptional upregulation of APXs, GPXs, and GSTs that encode antioxidant system components was observed in the leaves of J. curcas and poplar under water deprivation [5,10,65]. Consistently, the mRNA levels of GSTU18/19/20/28 and GSTF9/10 encoding glutathione transferases increased in water stress-exposed ginkgo leaves. Moreover, HSPs and LEAs are usually induced to help plants deal with water-deficit stress [8]. HSPs have diverse cellular functions such as detoxifying ROS, protein folding, as well as maintaining homeostasis against drought stress [66]. The precise functions of LEAs are still unknown, but they are supposed to participate in protective processes as hydration buffers, antioxidants, and both enzyme and membrane stabilizers [67]. Previous studies found that HSPs and LEAs were upregulated by drought in the leaves of diverse plants [5,10,68,69]. Consistent with these results, the mRNA levels of HSP17.4/18.2/21/22.0/81.1/101, LEA5, and HSF3/4A that might function as a positive transcriptional regulator of HSPs [48] increased in water stress-treated ginkgo leaves.
Aquaporins' cellular membranes are responsible for water transport and are essential for sustaining water balance in plants response to drought stress [70]. PIPs and TIPs were two main classes of water-transport membrane proteins that adjust cell turgor and maintain the water status in plants [71]. PIPs are located in the plasma membrane, and participate in water transport across it. TIPs encode tonoplast intrinsic proteins related to water transport across the tonoplast. Reduced transcript levels of PIPs and TIPs was reported in poplar leaves and maize ears [5,72]. Consistently, in this study, water stress led to downregulation of PIP2;8, and TIP1;1/2;1/4;1 in the leaves of G. biloba. In addition, aquaporin can act synergistically with XTR to regulate cell growth [73]. In this study, XTR6 was downregulated in the water stress-treated ginkgo leaves. We also found that XTH9, EXGT-A4, EXPA1/13, and EXLA1 involved in turgor-driven cell wall enlargement [74,75] were down-expressed in the leaves of G. biloba when exposed to water deficit. The suppressed cell wall enlargement can result in smaller cells, thus reducing plant leaf expansion. Therefore, these results suggested that transcriptional downregulation of PIPs, TIPs, XTR6, XTH9, EXGT-A4, EXPA1/13, and EXLA1 might play key roles in maintaining water balance and reducing cell growth in the leaves of water-deprived G. biloba.
Plant hormones have a great impact on stress signal transduction and adjusting the physiological responses to dehydration [76]. Variations in the transcript levels of genes involved in phytohormone biosynthesis often occur in plants in response to drought stress [1]. For instance, increased transcript levels of NCEDs, CYP707A1 and LOXs that encode crucial enzymes in the biosynthesis of ABA and JA were observed in the leaves of Arabidopsis, Jatropha curcas, and poplar [5,10,77]. Consistently, our study also found water stress-induced mRNA levels of NCED9, CYP707A1, LOX1, and LOX3 in ginkgo leaves. In actuality, the foliar concentration of ABA increased significantly by osmotically driven water stress in ginkgo leaves [32]. ABA and JA have been regarded as involving in stress-associated stomatal closure [1,78]. Therefore, upregulating NCED9, CYP707A1, LOX1, and LOX3 might improve drought tolerance by regulating stomatal closure in ginkgo leaves. Additionally, GA1 and GA3 involved in GAs biosynthesis were downregulated in the leaves of G. biloba exposed to water deficit. Similar results were also observed in water stress-treated maize and poplar [79]. Plants always reduced GAs concentrations to cope with water deprivation by stimulating osmotic adjustment [80]. Thus, reduced mRNA levels of GA1 and GA3 may promote water-deficit tolerance by maintaining the leaf turgor in ginkgo leaves.
Rapid and efficient recovery after removing drought is of significance for plant adaption and tolerance. For instance, increases in mRNA levels of genes associated with cell division and expansion as well as transport and reduced expression of genes that participated in ABA biosynthesis and antioxidants were found in Medicago truncatula, J. curcas, and Sorghum bicolor [10,81,82]. In line with these results, we found upregulation of GA20OX2, TIP4;1, XTR2, EXGT-A2, EXGT-A4 and EXLA1, which participate in the biosynthesis of GA, water transport, and cell wall enlargement, as well as downregulation of HSP18.2/21/70/90.7, GSTU18/19/28, GSTF9, and LEA5, which are related to antioxidant systems, in the leaves of G. biloba when treated with rewatering. These results revealed that the cell expansion inhibition and oxidative stress induced by water stress might be efficiently alleviated by rehydration in ginkgo leaves. Noticeably, the transcript levels of NCED4 and NCED9, which participate in ABA biosynthesis, increased under rehydration in ginkgo leaves. The higher mRNA levels of NCED4 and NCED9 suggest that water stress-induced enhancement of ABA biosynthesis is also active during the watered recovery stage in the leaves of G. biloba. In other words, ginkgo exhibited a memory for improving ABA biosynthesis even in the post-water stress stage. Memory genes with function in the ABA-regulated pathway after dehydration stress were also identified in other plant species [83,84], and they are effective means for plants to promote their resistance to environmental stress [83,84].

Gene Coexpression Networks Probably Contribute to the Response to Osmotic Water Stress and Rewatering in Ginkgo Leaves
Differentially expressed genes often co-ordinate in plants' response to drought and rewatering [85]. Hence, previous studies have dissected gene coexpression networks in plants' response to water deficit and rehydration [85][86][87]. A coexpression network underlying changes in the photosynthetic index and antioxidant enzymes has been constructed in the leaves of Pennisetum giganteum under drought and rehydration [19]. The coexpression networks of DEGs underlying the physiological response to water deprivation and rewatering have also been dissected in the root and leaves of Trachycarpus fortune [85]. In addition to these studies, we have identified several modules of coexpressed genes that positively correlate with the drought and rewatering treatment in ginkgo leaves. Noticeably, the royal blue, brown, light yellow, purple, and blue modules are highly associated with D1, D2, D3, D4, R1, and R2, respectively. This result suggested that the coexpressed genes in these five modules exert a crucial function in responding to osmotic water stress and rehydration in ginkgo leaves.
Furthermore, the hub genes (such as RAP2.4A and RAP2.4B from the royal blue module; HSF3, MGE2, GR-RBP3, and LRK10L1.2 from the brown module; HSP101, HSP17.8-CI, and HSP17.4-CIII from the light yellow module; UNI and CPR5 from the purple module; GA20OX2, MBF1C, and BAG6 from the blue module) in these five modules were identified. We propose that these hub genes act as key players in response to water stress and rewatering in the leaves of G. biloba. For instance, RAP2.4 can bind to both the ethyleneresponsive GCC-box and the dehydration-responsive element [56]. Overexpression of RAP2.4 brings about higher tolerance to water-deficit stress in A. thaliana [56]. HSF3 can act as a key positive regulator in HSPs' expression for heat and drought stress responses in plants [8,48]. LRK10L1.2 belongs to receptor kinase proteins (RLKs) that can sense environmental signals [32]. RLKs act as key genes responding to drought in ginkgo leaves, which has been described elsewhere [32]. Additionally, LRK10L1.2 positively participates in ABA signaling, and the lrk10l1.2 mutant exhibits lower tolerance to dehydration stress than wild-type A. thaliana [59]. Overall, these results suggested that gene coexpression networks represent important components of the coordinated response to water stress and rehydration in the leaves of G. biloba.

The Expression of Genes That Participated in the Biosynthesis of Flavonoids and Bilobalide Are Induced by Water Stress but Reduced by Rehydration in Ginkgo Leaves
Flavonoids and terpene lactones can function in eliminating ROS, blocking oxidation, and preserving plant growth during abiotic stress [1,88]. In addition, ginkgo leaves are famous for containing valuable flavonoids and terpene lactones [89,90]. Thus, it is important to understand the expression of genes involved in the biosynthesis of flavonoids and lactones terpene under water stress and rewatering conditions. Previous studies have shown that changes in transcript levels of genes related to flavonoids and lactones terpene biosynthesis were induced by water deprivation in Cistus creticus, Gossypium hirsutum, and Vitis vinifera [91][92][93]. Similarly, our study also found that CHS, FLS, HMGS1, and HMGR1/2, involved in flavonoids and bilobalide (as a kind of terpene lactones), were upregulated in water stress-treated ginkgo leaves. However, opposite results were observed in rewatering conditions. These results suggested that flavonoids and bilobalide might be elevated by water deficit but are probably decreased by rehydration. Moreover, these results imply that to obtain high-quality ginkgo leaves with more flavonoids and bilobalide, proper drought without rewatering could be applied before harvest.

Conclusions
In conclusion, transcriptome changes in response to water stress and rewatering were analyzed in ginkgo leaves. These DEGs, potentially involved in water stress and rewatering responses, include genes related to the biosynthesis of ABA, GA, and JA; ROS scavenging; the biosynthesis of flavonoids and bilobalide; water transport; and cell wall enlargement. Furthermore, weighted gene coexpression networks were constructed, and hub genes such as HSF3, LRK10L1.2, and HSPs, which act as key players in such networks, were identified. This study helps understand the molecular basis of water stress and rewatering responses in ginkgo leaves and supplies hints for evaluating crucial candidate genes for improving trees' drought tolerance.