Comparative transcriptome analysis of Saposhnikovia divaricata to reveal drought and rehydration adaption strategies

Water scarcity has become one of the most prevalent environmental factors adversely affecting plant growth and development. Different species have developed multiple ways of drought resistance. Saposhnikovia divaricata is a commonly used traditional herb in East Asia. However, limited information is available on the drought response of this herb and further clarification of underlying molecular mechanism remains a challenge. In this study, a comparative transcriptome analysis was firstly conducted to identify the major pathways and candidate genes involved in the drought adaptive response of S. divaricata. The seedlings of S. divaricata were subjected to progressive drought by withholding water for 16 days followed by 8 days of rehydration. Transcriptome analysis identified a total of 89,784 annotated unigenes. The number of differentially expressed genes (DEGs) gradually increased with the deepening of drought and decreased after rehydration. Gene Ontology enrichment analysis and Kyoto Encyclopedia of Genes and Genomes pathway analysis suggested genes related to oxidoreductase activity, carbohydrate metabolism, plant hormone signaling pathway and secondary metabolism were important in drought response of S. divaricata. Specific genes involved in reactive oxygen species scavenging system (POD, Cu/Zn-SOD, APX), abscisic acid and jasmonic acid signaling pathway (PYL4, PP2Cs, JAR1, JAZ) and phenylpropanoid biosynthesis (4CL, CCR, CAD) underwent dynamic alterations under drought and rehydration. Finally, the expression pattern of 12 selected DEGs from the transcriptomic profiling was validated by real-time quantitative PCR. Our study laid a foundation for understanding the stress response of S. divaricata and other Apiaceae family plant at molecular level.


Introduction
Water scarcity has become one of the most recurrent and prevalent environmental factors adversely affecting plant growth and development. Plants respond to water deficit by combining a complex regulatory network that induces various morphological, physiological and metabolic changes [1]. With the progression of drought, accumulation of reactive oxygen species (ROS) causes cellular damage, thus affecting plant growth [2]. ROS scavenging system consists of antioxidant enzymes and non-enzymatic antioxidants coordinate together enabling plants to overcome oxidative damages [3]. In addition, plants respond to drought stress by adjusting certain endogenous hormone levels. Under drought stress, abscisic acid (ABA) is a fundamental regulator of plant morpho-physiological responses [4]. Exogenous application of jasmonic acid (JA) alleviated drought stress by improving antioxidative enzyme activities and modifying osmolyte concentration in different species [5][6][7]. On the other hand, plants subjected to water deficit conditions often produce a variety of secondary metabolites [8] including flavonoids, phenylpropanoids, alkaloids and terpenoids, which probably functioning as antioxidant, osmoprotectants and plant hormone precursors. In general, plants respond to drought stress through activation of osmotic adjustment and antioxidant system, fluctuation of plant hormone contents and variation of metabolic pathways. Both drought intensity and frequency may affect the plant response.
Saposhnikovia divaricata (Turcz.) Schischk is a perennial herb belonging to Apiaceae family and widely distributed in East Asia. It is known as "Fang Feng" in China, "Bouhu" in Japan and "Bangpung" in Korea [9]. As a commonly used traditional medicine, the dried root of S. divaricata have long been studied, exhibiting diverse pharmacological functions such as antipyretic, analgesic, anti-inflammatory and anti-arthritic effects [10]. At present, given the increasing demand and lack of wild resources, S. divaricata is cultivated on large scale in northeast China. Water management in cultivation has become a key problem. Our previous research revealed that water deficit affected the growth of Astragalus membranaceus [11,12] and Glycyrrhiza uralensis [13], meanwhile enhanced the drought related gene expression and metabolic regulation. In this study, a comparative transcriptomic analysis of S. divaricata under progressive drought and rehydration was firstly conducted to identify the major pathways and candidate genes involved in the adaptive response of S. divaricata.

Plant material and drought treatment
Greenhouse pot experiment was conducted at a day/night temperature of 25 ± 5 o C/16 ± 3 o C with 50-70% relative humidity under 14/10 h light/dark cycle from March to July 2018. Seeds of S. divaricata were planted in 112 plastic pots (18 cm in diameter at the top and 25.5 cm deep). Three months later, seedlings were randomly divided into two groups making sure there are 6 seedlings in each pot. The control group (CK) was well watered while the treatment group was subjected to drought stress by withholding water for 16 days (DS) and then rewatered for 8 days (RW). On day 0, 4, 8, 12 and 16 of drought treatment and day 4 and 8 of subsequent rehydration, six pots of plants were randomly selected. The roots and fully expanded leaves were collected from each group for physiological measurements. The root samples of CK, DS8, DS12, DS16 and RW8 were collected and used for RNA sequencing. All the samples were immediately used for measurements or frozen in liquid nitrogen and stored at − 80 °C for further processing.

Physiological measurements
Soil water content was calculated as SWC (%) = (FW−DW)/ DW×100 [11]. Fresh soil was weighed as FW. Then, fresh soil was oven-dried at 105 o C to measure the dry weight as DW. The leaf relative water content was determined as LRWC (%) = (FW−DW)/(TW−DW)×100 [14]. After determining the fresh weight (FW), the petioles were submerged in water for 48 h to measure the saturated weight (TW) and then dried at 80 o C to determine the dry weight (DW). Fresh soil and leaves were collected from day 0, 4, 8, 12, 16 of drought as well as day 4 and 8 of subsequent rehydration. Each treatment was repeated for six times.
Superoxide dismutase (SOD) activity was determined as described by Dhindsa et al. [15]. Peroxidase (POD) activity was assayed based on the method of Hemeda and Klein [16]. Ascorbate peroxidase (APX) activity was determined according to the methods of Nakano and Asada [17]. Each treatment was repeated for three times.

RNA isolation, cDNA library construction and sequencing
Total RNA was extracted with TRIzol reagent (Invitrogen, CA, USA) following the manufacturer's protocol and genomic DNA was removed with DNase I (Takara, Dalian, China). 2100 Bioanalyzer (Agilent, CA, USA) and Nan-oDrop-2000 (Thermo Scientific, MA, USA) were used to ensure the samples were qualified for transcriptome sequencing. TruSeq™ RNA sample preparation kit (Illumina, CA, USA) was used for RNA-seq library construction. Illumina sequencing was conducted on an Illumina HiSeq™ 4000 platform. These experiments were completed by Shanghai Majorbio Bio-pharm Biotechnology Co. (http:// www. major bio. com, Shanghai, China).

De novo transcriptome assembly and functional gene annotation
Trinity software (https:// github. com/ trini tyrna seq/ trini tyrna seq/ wiki) was used for clean read assembly [18]. Using BLASTx tool with an E-value less than 10 −5 , the sequences of all assembled unigenes were functionally annotated to Clusters of Orthologous Groups (COG) database, Gene Ontology (GO) database, Kyoto Encyclopedia of Genes and Genomes (KEGG) database, Swiss-Prot database, Protein family (Pfam) database and non-redundant (NR) database.

Screening and analysis of differentially expressed genes
The abundance of all genes was normalized and calculated using uniquely mapped reads by the FPKM method (fragments per kilobase per million reads) [19]. A set of differentially expressed genes (DEGs) between the treatment group and the control group were gained by differential expression analysis using DESeq2 [20].

3
The screening criteria of DEGs were corrected with P value < 0.05 and |log 2 fold change| ≥1. Using the DEGs detected at different stages of drought stress and rehydration, GO and KEGG enrichment was implemented by Goatools and KOBAS platform [21,22].

Validation of transcriptome sequencing using real-time quantitative PCR
The expression of 12 DEGs was validated by real-time quantitative PCR (RT-qPCR) analysis. First-strand cDNA was synthesized by TransScript® Uni All-in-one First-Strand cDNA Synthesis SuperMix for qPCR kit (Transgen Biotech, Beijing, China). The primers were designed using the Primer Express 3.0.1 and listed in Supplementary Table S1. RT-qPCR was completed by TransStart® Tip Green qPCR SuperMix Kit (Transgen Biotech, Beijing, China) and run on ABI7500 Real-Time PCR System (Thermo Fisher Scientific, USA). EF 1-α was used as internal control [23]. The thermal cycling conditions were 30 s at 95 °C, followed by 40 cycles of 5 s at 95 °C, 34 s at 60 °C. Dissociation curve analysis was performed to determine the target specificity. The relative gene expression was calculated by 2 -ΔΔCt method and three biological replicates were conducted for all reactions.

Effects of drought and rehydration on relative water content
At 8 days after drought treatment, the SWC had no significant change. SWC dropped sharply to 13.2% and 9.5% on 12 and 16 days after treatment, respectively. After rewatering, SWC in the treatment was restored to the control level (Fig. 1a). LRWC maintained at a higher level from 0 to 12 days. LRWC decreased significantly from 90.5 to 42.6% on day 16 after drought, which was 53% lower than that of the control. Similarly, the LRWC reverted to the control level after rehydration (Fig. 1a).

Effects of drought and rehydration on the antioxidant enzyme activities
The activities of SOD, POD and APX in the root increased as the drought intensified. SOD and APX activities showed a similar trend and reached the highest level after 16 days of drought which increased by 57% and 111%, respectively ( Fig. 1b and d). In contrast, POD activity in the roots reached its maximum value on day 12. After rewatering, the activities of all three enzymes decreased to the normal level ( Fig. 1c).

Sequencing and de novo transcriptome assembly
The root samples of S. divaricata from CK, DS8, DS12, DS16 and RW8 were collected and used for RNA sequencing. A total of 116.66 GB clean data were obtained from 15 samples after removal of sequence adapter-containing reads, poly-N-containing reads and other low quality reads. On average, clean data was greater than 7.07 GB for each sample and Q30 value was more than 92% (Supplementary  Table S2). The filtered clean reads were assembled with Trinity program producing 167,447 unigenes with mean size of 807.49 bp and N50 length of 1253 bp (Supplementary Table S3). The mapped ratio was 76.80-79.08% and the unigenes ranged from 201 to 14,945 bp in length. These data suggested that the sequencing results were reliable to be used for subsequent analysis.  . 2a). Blasting into NCBI_nr database revealed that species most contributed to the annotation was Daucus carota (57.24%), followed by Quercus suber (4.39%) and Helianthus annuus (3.01%) (Fig. 2b).

Functional annotation and differential expression analysis
To identify DEGs in response to progressive drought and rehydration, differential expression analysis was conducted between the group of DS8, DS12, DS16, RW8 and CK ( Fig. 2c and d). Generally, the number of DEGs gradually increased with the intensification of drought and decreased after rehydration in comparison with the control. Venn diagrams showed that 145 DEGs were co-expressed in all four groups, with the highest overlap between DS12 and DS16, indicating similar expression patterns in DS12 and DS16 samples. Downregulated genes played a dominant role in moderate and severe treatment, while most of the DEGs upregulated during mild drought and after rewatering (Fig. 2d). For example, downregulated genes increased from 1831 to 5771 in DS12 and DS16, and then dropped to 676 in RW8.

GO and KEGG pathway enrichment analysis of differentially expressed genes
The up and downregulated DEGs of S. divaricata roots were further subjected to GO enrichment analysis to identify the major biological processes involved following the criteria of false discovery rate (FDR) < 0.01 (Table 1). Compared with the control, 117 upregulated DEGs at DS8 were correlated to oxidoreductase activity. Polysaccharide metabolism and O-glycosyl hydrolysis were enriched in the upregulated DEGs at DS8. However, the DEGs clustered to these GO terms were all downregulated at DS12 and DS16. Interestingly, zinc ion binding and transition metal ion binding function were enriched in the upregulated DEGs at DS16.    After rehydration, DNA-binding, protein kinase activity, carbohydrate metabolism and adenyl nucleotide binding function were enriched in a considerable number of upregulated DEGs suggesting the gradual completion of damage repair and cell recovery. KEGG pathway analysis of DEGs showed that pathways involved in drought response can be mainly divided into eight physiological categories, such as signal transduction, metabolism of terpenoids and steroids, lipid metabolism, biosynthesis of other secondary metabolites, carbohydrate metabolism, energy metabolism, replication and repair and environmental adaption (Fig. 3). Most of the pathways were enriched in upregulated DEGs at DS8 and RW8, whereas enriched in downregulated DEGs at DS12 and DS16. The highest enrichment factors among all treatment groups were plant hormone related signaling pathway. Linolenic acid metabolism, phenylpropanoid biosynthesis and glutathione metabolism pathway in the root were mostly enriched in upregulated DEGs, indicating the activation of these processes under mild stress. Starch and sucrose metabolism was enriched in upregulated DEGs at DS16 indicating carbohydrates were indispensable for plant to maintain survival under critical conditions.

DEGs involved in antioxidant system, plant hormone signal transduction and secondary metabolism
To further explore the molecular mechanism, a series of candidate genes were screened, including DEGs relevant to antioxidant system and plant hormone signal transduction, as well as DEGs concerning secondary metabolism.
27 genes annotated to antioxidant enzyme system were selected, of which 20 genes were related to POD (Fig. 4a). Most genes encoding POD were upregulated at DS8 and downregulated at DS12 and DS16 with the opposite for SOD, APX and glutathione peroxidase, which were mainly induced by moderate and severe drought. Glutathione S-transferase is the main detoxifying system for plants under extreme conditions. In this study, two glutathione S-transferase genes (DN103932_c0_g1, DN95127_c1_g5) were significantly upregulated under severe drought and the other two genes (DN91685_c0_g1, DN105946_c0_g2) mainly responded to mild stress and rehydration.
ABA participates in various physiological responses including drought induced stomatal closure and seedling growth inhibition [4]. ABA receptor genes (PYL4) first upregulated at DS8 and then downregulated at DS16, while its co-receptor protein phosphatase 2 C genes (PP2C) showed a reverse trend which in turn released protein kinase  Plant hormone signal transduction  11  0  0  18  10  42  97  11   8  7  3  0  7  4  1  0  0  0  t  n  a  l  p  -y  a  w  h  t  a  p  g  n  i  l  a  n  g  i  s  K  P  A  M   0  3  1  5  0  0  0  0  0  s  i  s  e  h  t  n  y  s  o  i  b  d  i  o  n  e  p  r  e  t  o  n  SnRK2s from inhibition (Fig. 4b). Jasmonic acid signaling pathway genes were mostly induced by drought suggesting their involvement in drought response of S. divaricata. Jasmonic acid-amido synthetase genes JAR1 (DN108879_c2_ g1, DN108238_c1_g1) significantly upregulated with the intensification of drought catalyzing the addition of L-isoleucine to jasmonate (Fig. 4b). Additionally, the expression of genes encoding jasmonate ZIM-domain protein (JAZ) and transcriptional activator (MYC2) also varied within a wide limit under drought and rehydration.
It has been reported that S. divaricata root contained high content of phenylpropanoids [8]. A total of 103 DEGs were enriched in phenylpropanoid biosynthetic pathway according to KEGG analysis. PAL, 4CL, CCR and CAD gene families were upregulated with the increase of drought. HCT, COMT and F5H mainly downregulated under moderate and severe drought and returned to the control level after rehydration (Fig. 5). Peroxidases were found to participate in the oxidative polymerization of phenylpropanoids [24]. Most peroxidase genes were induced by mild stress while a few were induced by moderate and severe drought.

Validation of drought-responsive genes by RT-qPCR
To ensure the reliability of sequencing data, the expression patterns of 12 genes identified as associated with water deficit responses were evaluated by RT-qPCR (Supplementary  Table S1). These genes demonstrate different functions and participate in various pathways, including 4 genes involved in plant hormone signal transduction, 4 genes involved in antioxidant system and 4 genes involved in secondary metabolism. As shown in Supplementary Fig. S1, the RT-qPCR results of 12 genes were highly consistent with the results of RNA-seq, which confirmed the reliability of the RNA-sequencing. Fig. 4 Response of selected genes involved in antioxidant system, plant hormone signal transduction and terpene metabolic processes during drought and rehydration of Saposhnikovia divaricata. a selected antioxidant enzyme related genes; b plant hormone signal transduction related genes

Discussion
Drought tolerance is a complex process of polygenic regulation in plants. In this study, a comparative transcriptome analysis was conducted to investigate the molecular mechanism underlying stress response of S. divaricata, a species of Apiaceae family without reference genome [25].
The activity of antioxidant systems has been used as one of the effective indicators of plant drought resistance [2]. In this study, the activities of SOD, POD and APX increased as the drought intensified suggesting that these enzymes work together to remove ROS which is consistent with our previous research [26] and similar results can also be found in other species of Apiaceae family such as Bupleurum chinense [27], Ferula assafoetida [28] and celery (Apium graveolens L.) [29]. Furthermore, we identified and selected 27 DEGs associated with ROS scavenging system. Among them, genes encoding POD were found to be the most abundant. According to KEGG pathway enrichment analysis, POD also participated in the phenylpropanoid metabolic pathway, catalyzing the final stage of lignin synthesis [30]. Most POD genes were upregulated after eight days of drought, while Cu/Zn-SOD and APX genes were mainly induced by moderate drought. It suggests that the enzymes play major roles in scavenging ROS are different at different stages of drought. In accordance with these observations, physiological response of S. divaricata root showed that POD activity reached its maximum value at DS12, while SOD and APX activities peaked at DS16.
Under limited water conditions, the increase of endogenous ABA level transcriptionally regulates downstream genes protecting plants from dehydration damage [4]. Signaling of the ABA is mediated by PYR/PYL/RCAR receptors, which bind to and inhibit PP2Cs releasing protein kinases SnRK2s from inhibition [31]. Overexpressing PtPYRL1 and PtPYRL5 in transgenic poplar plants significantly enhanced the osmotic stress resistance by positively regulating ABA signaling [32]. In our study, genes encoding PYLs significantly upregulated at DS8 and RW8, while genes encoding PP2Cs showed a reverse trend which upregulated at DS12 and DS16. SnRK2 genes were both up and downregulated under drought and rehydration. These results indicated that ABA signaling pathway was closely related to the drought accommodation of S. divaricata. JA and its derivatives are  significant in plant adaptive response to abiotic stress [33]. As for JA signaling, JAR1 genes were upregulated with the progression of drought suggesting the gradual accumulation of jasmonoyl-L-isoleucine. The expression of JAZ, a repressor of JA, was significantly downregulated at DS12 and DS16. Transcriptional activator MYC2 genes were more highly induced at the early stage. These findings were confirmed by previous studies which have shown that ABA and JA signaling pathways were both involved in priminginduced stress tolerance [34]. Instead of acting independently, it is more likely that plant hormones coordinate together contributing to better accommodation to environmental stress, which needs to be further examined. Plant secondary metabolism is the result of long-term evolutionary interaction between plant and environment and plays a critical role in plant environmental adaption, especially for perennial herbs. Drought could improve the accumulation of secondary metabolites in plants [7]. The main chemical components of S. divaricata root are phenylpropanoids, terpenoids and flavonoids [10]. KEGG pathway analysis showed a total of 103, 43 and 21 DEGs were enriched in phenylpropanoids, terpenoids and flavonoids biosynthesis, respectively. In S. divaricata root under drought, phenylpropanoids biosynthesis was enriched at DS8 and RW8 in the upregulated DEGs, whereas downregulation played a dominant role at DS12 and DS16. PAL and 4CL are two key genes in the early stages of lignin and flavonoid biosynthesis. CCR and CAD consecutively convert acyl-CoA to corresponding alcohols. In this study, two PAL, three 4CL, two CCR , and three CAD genes were found and most of them showed a drought induced pattern, upregulating with the deepening of drought. Plenty of peroxidase genes likely generate different transcripts for oxidative polymerization of different monolignol units. Most peroxidase genes were induced by mild stress and the other few were induced by moderate and severe drought. Taken together, it can be inferred that phenylpropanoids took part in the drought defense of S. divaricata root.

Conclusion
We reported the characterization of drought and rehydration responsive genes of S. divaricata based on comparative transcriptomic profiling. A total of 89,784 unigenes were annotated to at least one database. The number of DEGs gradually increased with the intensification of drought and decreased after rehydration. S. divaricata responded to drought and rehydration through a diversified and integrated network consist of genes relevant to ROS scavenging, plant hormone signaling pathway and phenylpropanoid biosynthesis. These results would provide useful information and serve as a transcriptomic database for stress resistance of S. divaricata and other Apiaceae family plant.