Tolerance strategies revealed in tree peony (Paeonia suffruticosa; Paeoniaceae) ecotypes differentially adapted to desiccation

Premise of the Study Tree peony (Paeonia suffruticosa; Paeoniaceae) is well known for its ornamental value, edible oil, and medicinal properties. However, its growing area has been limited by drought that has been exacerbated by global climate change. Methods Gene expression profiles of a drought‐tolerant cultivar and a drought‐sensitive cultivar during dehydration and rehydration were investigated by transcriptome analysis. Expression patterns of unigenes related to drought and recovery response and unrelated to either cultivar were classified by hierarchical clustering and real‐time quantitative PCR (qPCR). Results A total of 81,725 unigenes with a mean length of 762 nucleotides that may play roles in drought response were identified. Unigenes were characterized as being involved in lipid transport metabolism, proline metabolism, and photosynthesis. In addition, plant hormone signaling pathway genes were also characterized as potentially being involved in drought response. Expression patterns of the 20 drought‐responsive unigenes verified by qPCR showed a differential expression pattern under either the drought or recovery treatment. Discussion This is the first report to identify and verify unigenes of tree peonies with differing water sensitivity during dehydration and rehydration. This study offers a valuable resource for candidate genes involved in drought and provides insight into the breeding of drought‐resistant tree peony cultivars.

Tree peony (Paeonia suffruticosa Andrews; Paeoniaceae) has been cultivated for more than 1600 years in China, and there are approximately 2000 tree peony cultivars worldwide (Wang, 1997). The first written description of this genus was in 200 BC as a medicinal plant. In the fifth century, with the selection of multiple flower shapes and colors, it became known as an ornamental plant (Li et al., 2011). Studies have shown that the seed oil of tree peony contains abundant unsaturated fatty acids that are beneficial to human health (Sarker et al., 1999;Su et al., 2016). Owing to its multiple uses, the peony genus has spread through Asia, to the Mediterranean, Caucasus Mountains, the mountainous regions of Europe, the United States, and Australia (Rogers, 1995).
The greatest number of cultivated varieties and the largest distribution area of tree peony are found in central China, which has long been characterized by an arid climate. Nevertheless, severe water deficiency stress can limit the cultivation area, lead to smaller leaves and flowers, inhibit the synthesis of organic substances and flower pigments, and reduce the ornamental value and seed yield of tree peony (Li et al., 2011). Recent studies, however, have mainly focused on oil extraction (Chen et al., 2016a, b;Cui et al., 2016;Han et al., 2016). In addition, efficient protocols for the micropropagation of tree peony and the effects of different medium compositions and exogenous hormones on the browning of tree peony callus in tissue culture were recently investigated (Wen et al., 2016;Zhou et al., 2016a). Surprisingly, however, the desiccation tolerance strategies in tree peony cultivars have not yet been investigated.
Transcriptome analysis that uses deep sequencing technology now permits large-scale gene expression detection in the absence of a reference genome. Although there have been several investigations of transcriptome sequencing of tree peony (Gai et al., 2012;Zhou et al., 2013;Zhang et al., 2014Zhang et al., , 2015Zhao et al., 2014;Barghini et al., 2015;Li et al., 2015;Shi et al., 2015;Gao et al., 2016;Wang et al., 2016b), studies of drought-responsive differential expression genes in tree peonies have not yet been reported in the literature. Two separate studies of reference gene selection in tree peony-one in plants with different flower colors and another during flower development-were recently reported Zhou et al., 2016b).
Screening of drought-tolerant tree peony cultivars revealed that 'Luo Yang Hong' (LYH) is tolerant to drought, whereas 'Wu Long Peng Sheng' (WLPS) is tolerant to flooding (Kong et al., 2011), making the two cultivars ideal study material for investigating mechanisms of drought response in plants. With a view toward improving plant structure, perfecting bloom quality, and mitigating damage from desiccation, this study used LYH and WLPS with their contrasting water sensitivity to characterize unigenes during dehydration and rehydration to explore the complex mechanisms of drought response networks.

Plant material treatment and sample collection
Four-year-old LYH and WLPS seedlings were cultured in pots using soil collected from the Luoyang National Peony Garden (Luoyang, China). The pots were buried deep in the ground from October until May to avoid freezing injury. The pots were then dug out and irrigated once every two days and cultured under natural conditions before they were used for the water deficiency treatments. Five individuals were used per treatment. The drought treatment (DR) was initiated by tap water irrigation until the soil moisture content reached 80%, after which plants were dehydrated naturally for seven days until the leaves had severely wilted. For the rehydration treatment (RE), the tree peonies were re-watered until the soil moisture content again reached 80%, after which they were cultured for one more day to let the leaf blades completely unfold. Tree peonies cultured in pots with a constant soil moisture content of 80% served as the control treatment (CK). The soil moisture content was measured by gravimetric methods (Bao, 2000). For all three treatments, the third and fourth leaves were sampled and immediately frozen in liquid nitrogen and stored at −80°C.

RNA extraction, cDNA library construction, and sequencing
The leaves sampled from the three different treatments (DR, RE, and CK) were assigned to six independent pools. Total RNA was extracted using the modified cetyltrimethylammonium bromide (CTAB) method (Gambino et al., 2008). The integrity of RNA was examined by an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, California, USA). An in-house library preparation method was used for mRNA sequencing, using fragment sizes of 200 bp. Libraries were quantified by an Agilent Bioanalyzer and qualified by the ABI StepOnePlus Real-Time PCR System (Thermo-Fisher Scientific, Waltham, Massachusetts, USA) during the quality control steps. Libraries were sequenced using the Illumina HiSeq 2000 (Illumina, Shenzen, Guangzhou, China) at the Beijing Genomics Institute (BGI). Each library was run on a separate lane of the HiSeq. The cDNA library was deposited in the National Center for Biotechnology Information (NCBI) Transcriptome Shotgun Assembly database (BioSample accession no. SRS1180651).

De novo assembly and protein-coding region prediction
Reads with adapters, unknown nucleotides larger than 5%, and lowquality reads (bases quality ≤10) were discarded. Only reads longer than 90 bp were used for assembly. Reads from all treatments and/ or cultivars were assembled together by Trinity 3.4 (Grabherr et al., 2011; open source code publicly available at http://TrinityRNASeq. sourceforge.net). Because a reference genome is not available for tree peony, reads were mapped to the assembled unigene set.
Unigenes were first aligned by BLASTX (E-value <0.00001) to protein databases in the following order of NCBI's nonredundant protein database (nr), Swiss-Prot, Kyoto Encyclopedia of Genes and Genomes (KEGG), and NCBI's Clusters of Orthologous Groups (COG) database. Unigenes aligned to a higher-priority database were not aligned to lower-priority databases. The best alignment results were used to decide sequence direction of unigenes. When a unigene was not aligned with any of the above databases, ESTScan was used to determine its sequence direction (Iseli et al., 1999).
Proteins with the highest ranks in the BLAST results were used to decide the coding region sequences of unigenes, then the coding region sequences were translated into amino sequences with the standard codon table. Unigenes that could not be aligned to any database were scanned by ESTScan (Iseli et al., 1999), producing nucleotide sequence (5′-3′) direction and amino sequence of the predicted coding region.

Real-time quantitative PCR (qPCR) verification analysis
Twenty dehydrin-related unigenes were selected for the assessment of expression profiles. Total RNA was converted into singlestranded cDNA using an M-MLV reverse transcriptase (Promega Corporation, Madison, Wisconsin, USA). The ABI StepOnePlus Real-Time PCR System (Thermo-Fisher Scientific) was utilized to perform the expression profile verification. The reaction was carried out as described in our previous publication (Pang et al., 2015), using three biological replicates per sample. The relative gene expression levels of the selected unigenes were normalized to 18S rRNA and calculated using the 2-ΔΔCt method (Livak and Schmittgen, 2001). The primer sequences are shown in Table 1.

Sequencing output statistics, assembly metrics, and proteincoding region classification
The total number of clean nucleotides generated from the six libraries of the two tree peony cultivars (LYH and WLPS) by three treatments (CK, DR, and RE) exceeded 4.6 Gbp. The number of clean reads of LYH-CK, LYH-DR, LYH-RE, WLPS-CK, WLPS-DR, and WLPS-RE were 54, 52, 52, 52, 51, and 52 million, respectively. The Q20 values exceeded 97%, and the GC contents were approximately 46% of all samples, which meant that the sequencing data were robust for further analysis. After assembling all sequences from all samples, 81,725 unigenes were obtained. Their aggregate length was 62,310,011 nucleotides, with a mean length of 762 nucleotides.
A total of 41,808 protein-coding regions from the unigenes were translated into amino sequences. Detailed information of their length distributions is given in Appendix 1. The COG analysis showed that 14,768 unigenes were assigned to 25 classifications. The largest category was 'General function prediction only'; 'Transcription, ' 'Replication, recombination, and repair, ' 'Post transcriptional modification, protein turnover, chaperones, ' and 'Signal transduction merchanism' were comparatively high; and 'RNA processing and modification' had the smallest number of responding unigenes (Fig. 1).

Functional annotation of the unigenes
A total of 43,977 unigenes were successfully allocated to the three main gene ontology categories: biological process, molecular function, and cellular component (Fig. 2). The biological process category contained 22 classes subsumed under five larger groups: cellular process (18,310), metabolic process (17,809), singleorganism process (12,415), stimulus (8451), and biological regulation (6850). The cellular components category consisted of 17 classes, dominated by cell (21,823), cell part (21,822), and organelle  Significantly enriched gene ontology terms with a larger cluster frequency (i.e., ≥9%) were identified. For the biological process category, the oxygen-containing compound and oxidation-reduction process were more common not only among treatments (CK, DR, and RE) within a cultivar but also between the cultivars (LYH and WLPS) for a given treatment. The unigenes responded to stimuli, including abiotic, endogenous, biotic, and chemical stimuli that were identified extensively in LYH and WLPS. Unigenes' response to stress was detected only in the CK vs. RE treatments of LYH (Appendix 2). For the molecular function category, unigenes were identified in DR vs. RE of LYH, while oxidoreductase and catalytic activity were largely detected both between cultivars and among the treatments. None of the other candidate molecular functions were identified, except the above two functions between LYH and WLPS (Appendix 3). For the cellular component category, the membrane, cell periphery, plasma membrane, and extracellular region were extensively detected in the LYH and WLPS treatments separately. Within the same treatment, however, none of the unigenes responded to the cellular component between the cultivars (LYH and WLPS) apart from two exceptions: the membrane detected in DR and the extracellular region identified in RE (Appendix 4).  genotype yet responding to dehydration and rehydration were detected. Hierarchical clustering of LYH-CK, LYH-DR, LYH-RE, WLPS-CK, WLPS-DR, and WLPS-RE indicated that CK, DR, and RE were clustered together, revealing dehydration-and rehydration-responsive unigenes (Fig. 4). We note that our differential expression results are preliminary, as our experiments lack replication. Follow-up experiments will be needed to fully validate our observations.

Pathway analysis using the KEGG database
BLAST analysis of 81,725 unigenes against the KEGG database was performed to analyze gene products during metabolism processes and related gene functions in cellular processes. A total of 23,518 unigenes were involved in 128 KEGG pathways. More than one-fifth (21.77%) were related to metabolic pathways, whereas 11.82% were related to the biosynthesis of secondary metabolites, 5.51% were related to plant-pathogen interaction, and 4.69% were related to plant hormone signal transduction ( Table 2). The KEGG pathways of unigenes with annotation for the two cultivars (LYH vs. WLPS) within the same treatment and among treatments (CK, DR, and RE) within the same cultivar were also analyzed. The metabolic pathways, biosynthesis of secondary metabolites, plantpathogen interactions, and plant hormone signal transduction had similar percentages in all pathways identified. Abscisic acid, jasmonic acid, ethylene, brassinosteroids, salicylic acid, gibberellins, cytokinin, and auxin signaling pathways were all detected in this study. It is interesting to note that when analyzing the pathway between LYH and WLPS, the plant hormone signal transduction pathway was not detected in CK and RE (Table 3). LYH and WLPS exhibited higher plant-pathogen interaction after dehydration than the control, which might be caused by plant interaction with a pathogen such as arbuscular mycorrhizal fungi. However, it was unclear why a plant-pathogen interaction apparently went undetected after rehydration.

Unigene validation by qPCR
To validate the expression profiling of the dehydrin-responsive unigenes, 20 genes predicted to participate in the dehydrin response pathway were selected (1) to determine their relative expression in the dehydration (DR), rehydration (RE), and nontreatment of tree peony (CK) and (2) to validate the transcriptome sequencing results. Abundance of the target genes was normalized relative to the abundance of 18S RNA; the Ct values (i.e., the number of cycles corresponding to the inflection point from baseline to exponential growth) of 18S rRNA for all samples ranged from 24.0 to 26.0. The results of the qPCR verification showed a differential expression pattern under both the DR and RE treatments. Dehydrin Xero 2-like was significantly upregulated after dehydration but then downregulated during rehydration of tree peony seedlings (Table 4).

DISCUSSION
Drought stress is one of the main abiotic stresses, and it may alter plant growth, metabolism, and yield (Ajithkumar and Panneerselvam, 2014). In tree peonies cultivated in central and northwestern China, water deficiency is a common problem. This drought stress limits the growth of leaves and flowers, inhibits the synthesis of organic compounds and anthocyanin, and reduces seed yield (Li et al., 2011). Plants that receive drought signals initiate a range of physiological, morphological, and biochemical defense responses at both the cellular and molecular level (Verslues et al., 2006). An overexpression of genes in response to drought stress could alleviate drought-induced damage while promoting plant growth. Drought tolerance strategies, as revealed by transcriptome sequencing in poplar (Barghini et al., 2015) and sorghum (Fracasso et al., 2016), have uncovered a number of drought-responsive genes. To diminish the damage caused by water deficiency in tree peony, two cultivars with contrasting water sensitivity were selected for unigene characterization to investigate the molecular mechanisms driving their drought response.
Plants can adapt to desiccation stresses and stay alive by alternating the accumulation of osmolytes (Parida et al., 2007). Proline, one of the most important osmolytes, is quickly accumulated and involved in the plant response to dehydration to maintain a cellular balance of water content and turgor potential (Vendruscolo et al., 2007). A previous study showed that proline accumulates after dehydration and then decreases to the initial level after rehydration in both LYH and WLPS tree peony cultivars (Li et al., 2014). Proline variation during dehydration and rehydration corresponded to transcriptome analysis and was consistent with that reported in previous publications (Gechev et al., 2012;Hossain et al., 2016). In the present study, numerous unigenes related to the proline metabolism process, including proline-rich proteins (Unigene8963_All, Unigene11447_All), proline transporters (Unigene11945_All), and an osmotic precursor (CL5318.Contig1_All) were identified in LYH and WLPS under the dehydration and rehydration treatments. One of the most harmful effects of drought stress is increased production of reactive oxygen species (ROS) (Miller et al., 2010;Bartwal et al., 2016). When plants suffer drought stress, active oxygen metabolism is strengthened, which generates a large amount of O 2 − , H 2 O 2 , and OH − (Bian and Jiang, 2009). However, plants possess an evolved antioxidant defense system that enables them to maintain ROS at a low quantity to protect cells from excessive and permanent oxidative damage. Therefore, the ability of plants to clear the active oxygen can reflect their ability to tolerate stress (Hossain et al., 2016). Induction of the antioxidant enzymes' expression level makes the antioxidant system an efficient mechanism to control ROS accumulation, both temporally and spatially (Hossain et al., 2016). These redox-sensitive proteins may be oxidized by ROS directly or indirectly via non-enzymatic compounds, such as glutathione (GSH), which are major players in redox signaling when antioxidants are involved (Shao et al., 2008;Klumpen et al., 2016;Thangamani et al., 2016). The unigene corresponding to glutamate (CL7346.Contig2_All) displayed a different responsive pattern in tree peony when exposed to dehydration and rehydration.
Large pigment-protein complexes are the most significant factors that function during photosynthesis (Qin et al., 2015). Lightharvesting complex stress-related proteins catalyze excess energy dissipation in photosystems I (PSI) and II (PSII) (Pinnola et al., 2015). The unigenes involved in photosynthesis include 14 subunit complexes in PSI and another 14 subunit complexes in PSII, which were examined in the present study. Specifically, the chlorophyll a/b-binding protein that encodes the light-harvesting complex (LHC), which captures sunlight and transfers the excitation energy to the core in higher plants, was obtained. Four unique subunits (PsaG, PsaH, PsaN, and PsaO) of the PSI-LHCI super complex in higher plants were detected in the present study. As a member of the ROS, H 2 O 2 and its production occur in the chloroplasts, peroxisomes, and mitochondria of plants. Drought-treated plants had a significantly increased ROS content and diminished operating and  maximum efficiencies of PSII photochemistry (Ryan et al., 2014). Reduced photosynthetic pigment contents resulting from drought stress might decrease ROS formation by regulating chlorophyll synthesis and other components of the photosynthetic machinery. ROS generation is considered to be closely related to lipid peroxidation under drought stress (Ryan et al., 2014;Uzilday et al., 2014). For example, lipid peroxidation analysis showed that transgenic Agrostis stolonifera L. root exhibited less cellular damage when compared with the wild type under drought stress conditions (Xu et al., 2016). The alleviation of adverse effects of drought stress is partially attributable to an increased antioxidant ability and decreased lipid peroxidation induced by early ROS accumulation (Xing et al., 2016). Tobacco plants treated with low and moderate levels of riboflavin accumulated higher levels of ROS and lipid peroxide with enhanced drought tolerance (Deng et al., 2014). In the present study, 666 unigenes involved in lipid transport and metabolism were identified according to the COG classification. This clearly illustrates the close relationship between lipid peroxidation and the drought stress response in plants.
Calcium mobilization is one of the downstream events modulated by H 2 O 2 (Neill et al., 2002). The calcium ion (Ca 2+ ) functions as a secondary messenger in modulating diverse physiological processes that are important for stress adaptation in plants. Both Ca 2+ and Ca 2+ /calmodulin (CaM)-binding protein and transcription factors have been identified, and their functional analysis suggests that they play key roles in plant stress signaling pathways (Reddy et al., 2011). Previous studies have indicated that drought stress activates ABA-dependent and ABA-independent gene expression (Yoshida et al., 2014). The cis-regulatory element ABA-responsive element (ABRE) (CACGTG [T/C/G]) and their coupling elements ([C/A]ACGCG[T/C/G]) in the upstream region were observed in the upregulated genes (Kaplan et al., 2006). Hence, it was concluded that for some specific Ca 2+ transients, ABREs function as Ca 2+ -responsive cis-regulatory elements (Reddy et al., 2011).
Heat stress can trigger the higher expression of heat-shock proteins (HSPs), which might coordinate with other stress-response mechanisms to mitigate cellular damage and re-establish cellular homeostasis (Wang et al., 2004). Copper applied to tree peony revealed an increase in dehydration-responsive element-binding (DREB) protein . In the present study, we identified one class II HSP isoform 1 (CL7474.Contig3_All) and one heavy metal-associated isoprenylated plant protein (HIPP) (Unigene32639_All), both of which were unrelated to genotype but responsive to dehydration and rehydration. Regulation of HSP and HIPP by dehydration and rehydration in tree peony illustrates the synergistic interaction of drought with other stress-response mechanisms to alleviate cellular damage and re-establish cellular homeostasis.

CONCLUSIONS
Transcriptome profiling analysis demonstrated unigene response to dehydration and rehydration in tree peony, namely MYB, AP2/ ERF, NAC, bHLH, RING-H, HSP, and HIPP. These newly identified unigenes will increase our understanding of drought stress-responsive mechanisms, and they may be quite useful as novel genes for the molecular breeding of tree peony to improve its drought tolerance. Further research is necessary to reveal and understand how antioxidant enzymes interact with key hormones in the signaling responses of plants to drought stress.

ACKNOWLEDGMENTS
This work was supported by the National Natural Science Foundation of Henan Province (grant no. 162300410079), the Education Department Project of Henan Province (grant no. 17A180003), and the Innovative Research Team in Henan University of Science and Technology (grant no. 2015TTD003).

DATA ACCESSIBILITY
The cDNA library was deposited in the National Center for Biotechnology Information (NCBI) Transcriptome Shotgun Assembly database (BioSample accession no. SRS1180651).