Transcriptome Analysis of Jujube (Ziziphus jujuba Mill.) Response to Heat Stress

Heat stress (HS) is a common stress influencing the growth and reproduction of plant species. Jujube (Ziziphus jujuba Mill.) is an economically important tree with strong abiotic stress resistance, but the molecular mechanism of its response to HS remains elusive. In this study, we subjected seedlings of Z. jujuba cultivar “Hqing1-HR” to HS (45°C) for 0, 1, 3, 5, and 7 days, respectively, and collected the leaf samples (HR0, HR1, HR3, HR5, and HR7) accordingly. Fifteen cDNA libraries from leaves were constructed for transcriptomics assays. RNA sequencing and transcriptomics identified 1,642, 4,080, 5,160, and 2,119 differentially expressed genes (DEGs) in comparisons of HR1 vs. HR0, HR3 vs. HR0, HR5 vs. HR0, and HR7 vs. HR0, respectively. Gene ontology analyses of the DEGs from these comparisons revealed enrichment in a series of biological processes involved in stress responses, photosynthesis, and metabolism, suggesting that lowering or upregulating expression of these genes might play important roles in the response to HS. This study contributed to our understanding of the molecular mechanism of jujube response to HS and will be beneficial for developing jujube cultivars with improved heat resistance.


Introduction
Abiotic stresses, such as heat, drought, salinity, and cold, are major environmental constraints to crop production and food security all over the world [1,2]. In particular, with extreme weather and global warming, heat stress (HS) has received increasing concern and interest [3]. Although increasing temperatures are beneficial for crop production in some cooler regions of the world, their overall impact on global food security is negative [4]. HS harms cellular homeostasis and causes leaf etiolation, severe retardation in growth and development, increased risk of disease, and even death [5]. Temperature increases reduces global yields of major crops, such as wheat, rice, maize, and soybean [6], in addition to horticultural crops such as grapevine, almonds, apples, oranges, and avocados [7].
Jujube (Ziziphus jujuba Mill.) belongs to the Rhamnaceae family in the Rosales order [8] and is one of the oldest culti-vated horticultural crops with a long domestication history [9]. Jujube fruit is rich in vitamin C, phenolics, flavonoids, triterpenic acids, and polysaccharides and is widely consumed as a food or food additive [10]. It is now a major dry fruit crop with a cultivation area of 2 million ha, the main source of income for 20 million farmers, and a traditional herbal medicine for more than one billion people in Asia [11]. Jujube is well adapted to various biotic and abiotic stresses, especially drought and salinity, and is considered an ideal cash crop for arid and semiarid areas where common fruits and grain/oil crops do not grow well [8].
Xinjiang province of China is the core area of the arid region in Central Asia, which is one of the most arid regions in the world [12]. Jujubes are among the main agroeconomically important crops in Xinjiang, and those from this region have good quality and the highest production worldwide [13]. However, Xinjiang is particularly vulnerable to climate change and has experienced significant climate warming in the past 40 years [14]. In recent years, HS has dramatically and repeatedly affected the production of quality jujube fruits in Xinjiang. Therefore, identifying and breeding heat-tolerant jujube cultivars might be a feasible and important strategy for protecting the production and quality of jujube fruits.
Turpan, one of the northeastern cities of in Xinjiang province of China, has a unique temperate continental arid desert climate, with bright sunshine, high temperatures, and large day-night differences in temperature [15]. We found a jujube cultivar ("Hqing1-HR") in an orchard of Turpan by chance [16] and bred it in our laboratory successfully. To investigate the transcriptomic change in "Hqing1-HR" response to HS, we subjected its seedlings to 45°C stress. At 0, 1, 3, 5, and 7 days after HS-treatment, we assessed phenotypic and physiological features and collected samples for RNA sequencing experiments. This study may provide new insight into transcriptional alterations in heatresistant jujube cultivars responding to HS.

Plant Materials, Heat Treatment, and Sample Collection.
Green cuttings of "Hqing1-HR" were collected from the jujube orchard of Turpan in Xinjiang, China. The green cuttings were grown under greenhouse conditions with an automatic spray system (20~35°C with 90% humidity). When the cuttings had 7~9 true leaves, a total of 80 plants were transferred to a controlled growth chamber with a light/dark regime of 14/10 h at 30/20°C, 80% relative humidity, and light intensity of 300 μmol m −2 s −1 of photosynthetically active radiation.
Seedlings with 14 true leaves were cultured in the same chamber under the same conditions except with the temperature at 45°C. After 0 (control), 1, 3, 5, and 7 days of heat treatment, the 10th true leaves counting from bottom to top were collected from three different samples as biological repetitions. Leaf samples were immersed in liquid nitrogen and stored at -80°C for transcriptome sequencing.

Phenotypic
Determination of Jujube Leaves. Samples from the same part of each leaf were rinsed and fixed with FAA solution (70% ethanol) at 4°C. Leaves were freezedried after dehydration using an alcohol gradient series to the critical drying point and then stuck to the table using conductive tape. Samples were coated with a Pt film using an ion sputtering instrument (Hitachi E-1045), and a SUPRA 55VP scanning electron microscope (Zeiss, German) was used for observation of leaves at an accelerating voltage of 2.00 kV.
2.3. RNA Extraction, cDNA Library Construction, and Illumina Sequencing. Total RNAs were extracted from 15 samples representing three biological replicates of "Hqing1-HR" at five treatment stages (0, 1, 3, 5, and 7 days of HS) using an RNAprep Pure Plant Kit (Tiangen, Beijing, China) according to the instructions of the manufacturer. Extracted RNA was treated with DNase I (Promega, Madison, WI, USA) to remove DNA. RNA quality and quantity were assessed using a Nanodrop 2000 Spectrophotometer (Thermo Fisher Scientific, Wilmington, DE, USA) and an Agilent Bioanalyzer 2100 System (Agilent Technologies, Santa Clara, CA, USA), respectively. The integrity of RNA was confirmed by 1% (w/v) agarose gel electrophoresis.
RNA samples containing equal amounts of RNA were pooled from three independent individuals and then used for library preparation and sequencing. The resulting libraries were sequenced using an Illumina HiSeq X-ten platform with paired-end 150 bp reads.
2.4. RNA-seq Read Processing and Assembly. Raw reads were generated by the Illumina HiSeq X-ten genome analyzer and were analyzed using Fast Q to assess the base quality. Reads were cleaned by removing adaptor sequences, low-quality sequences including empty reads, and sequences containing <10% bases with a Phred quality score < 20. The transcriptome was assembled using StingTie V1.3.1 [17], and the remaining cleaned reads were mapped to the jujube reference genome sequences [11] using HISAT2 [18] with default settings.

Bioinformatic
Analysis. FPKM (fragments per kilobase per million mapped reads) was used to evaluate the expression level of genes. The software edge R [19] was used to measure the FPKM values and identify differentially expressed genes (DEGs). Genes with RPKM < 0:1 in every sample were removed before analysis. DEGs were determined on the basis of fold change (FC) (FC ≥ 2 or ≤ 0:5) and false discovery rate (FDR) (FDR < 0:01).
To predict gene function and calculate the distribution frequency of functional categories, Gene Ontology (GO) analyses were employed using DAVID bioinformatics resources [20]. Venn diagrams were generated using a tool available online (http://bioinfogp.cnb.csic.es/tools/venny/).

2.6.
Validation of RNA-seq Using Quantitative Reverse-Transcription PCR (qRT-PCR). To verify the validity of the RNA-seq data, qRT-PCR experiments were performed for randomly selected DEGs. Primers are presented in Table 1. The same RNA samples used for RNA-seq were used for qRT-PCR. In each pooled sample, l μg of RNA was reverse transcribed using a Prime Script™ RT Reagent Kit (Takara, Dalian, China) according to the instructions of manufacturer. The qPCR was performed on a Bio-Rad S1000 with Bestar SYBR Green RT-PCR Master Mix (DBI Bioscience, Shanghai, China). PCR conditions were as follows: denaturing at 95°C for 8 min, 38 cycles of denaturing at 95°C for 15 s, and annealing and extension at 60°C for 1 min. Relative gene expression was calculated using the Livak and Schmittgen 2 −ΔΔCt method [21], normalized with the reference gene ZjH3 of jujube. PCR amplifications were performed in triplicate for each sample.

Statistical Analysis.
All values of all data are presented as mean ± standard deviation (SD). To determine the significance of differences between means, Student's t-test (paired) was implemented, and a value P < 0:05 was considered statistically significant.

Phenotype of Jujube Seedlings Post HS.
To obtain an overview of the heat-tolerance phenotype of "Hqing1-HR," seedlings with 14 true leaves were subjected to treatment with HS (45°C). At 0 (control) 1, 3, 5, and 7 days after the treatment, none of the seedlings displayed withered leaves ( Figure 1(a)), suggesting this jujube cultivar might be of heat tolerant.
Stoma are important channels for gas and water exchange between plants and the atmosphere and can make adaptive adjustments under various stress conditions. We assessed the stomatal density and stomatal opening rate of leaves from each group. Stomatal density and stomatal opening rate were significantly increased post heat treatment, and they showed a trend of first increasing and then decreasing with the extension of heat treatment ( Figure 1 Table 2). This suggested that the "Hqing1-HR" could reduce the damage by passively changing stomatal density and opening rate.
Through the Illumina HiSeq X-ten platform, we generated over 0.402 billion paired-end reads, corresponding to an average of 26.8 million sequence reads per sample. Using HISAT2 [18] with default settings, approximately 68.9% clean reads were mapped to the jujube reference genome [11].
To understand the spatiotemporal expression patterns of all samples, we performed principal component analysis (PCA). The three samples collected at each time point formed independent clusters (Figure 2(a)). Moreover, Pearson's correlation analysis for all pairs of RNA-seq samples was performed, demonstrating similar results (Figure 2(b)) and indicating that gene expression in the three replications of every sample was homogeneous (Figure 2(a)). This suggested that the replicated samples produced data acceptable for further analyses.
In HR1 vs. HR0, we identified 902 upregulated genes and 740 downregulated genes. There were 1,850 upregulated and 2,230 downregulated genes in HR3 vs. HR0. In HR5 vs. HR0, 2,167 upregulated and 2,993 downregulated genes were discovered. Meanwhile, 1,019 upregulated and 1,100 downregulated genes were identified in HR7 vs. HR0. The numbers of the upregulated and downregulated DEGs were similar in each comparison, indicating that heat stress promoted and inhibited the transcription of numerous genes. Moreover, there were 717 common DEGs among the four comparisons ( Figure 3(b)).
Although leaves of none of the seedlings become withered under high temperature, multiple DEGs were associated with "response to stress" (GO: 0006950) and "response to heat" (GO: 0009408) terms in all four comparisons ( Figure 4 and Tables S5-S8), indicating that "Hqing1-HR" might be sensitive to HS under normal conditions but establish a new steady-state balance of biological processes enabling the organism to function and survive well at higher temperatures (45°C). Analysis of DEGs enriched in "response to stress" and "response to heat" terms indicated that the expression levels of multiple DEGs associated with the response to HS were clearly upregulated after HS. Previous study indicated that cells can activate an ancient signaling pathway leading to the transient expression of HSPs in response to heat stress [22]. Indeed, it was found that multiple DEGs encode the heat shock proteins (HSPs), including HSP17 (gene12298), HSP18 (gene3931 and gene3933), HSP21 (gene2239), HSP 22 (gene7584), HSP 23 (gene6955), HSP26 (gene22046), HSP70 (gene2890 and gene4042), HSP83 (gene21597), HSP90 (gene9467), HSF30 (gene12331), and HSC-2 (gene22447) in the current study.
To explore how HS affected the normal physiology of chloroplasts and photosynthesis, we analyzed the DEGs associated with "chloroplast organization" (GO: 0009399), "chloroplast RNA processing" (GO: 0031425), "photosynthesis, light harvesting" (GO: 0009765), and "photosynthesis" (GO: 0015979) terms. To our surprise, most of the DEGs enriched in these four terms above were upregulated by HS. This suggested that HS might not disrupt the physiology of chloroplasts or photosynthesis and instead might promote the photosynthesis in "Hqing1-HR." 3.6. Metabolisms Is Affected by Heat Stress. HS always affects the global metabolism of plants. "Myo-inositol hexakisphosphate biosynthetic process" (GO: 0010264) was the only common term associated with metabolism identified in HR1 vs. HR0, HR3 vs. HR0, HR5 vs. HR0, and HR7 vs.
HR0. However, we identified multiple specific terms associated with metabolism in the four comparisons ( Figure 4 and Tables S5-S8).

Validation of RNA-seq by qRT-PCR.
To verify the reliability of our transcriptome data, six DEGs were randomly selected for expression analysis using qRT-PCR experiments ( Figure 5). The expression patterns shown by the qRT-PCR results ( Figure 5(b)) were consistent with the RNA-seq results ( Figure 5(a)), with PCCs higher than 0.9.

Discussion
The global air temperature is predicted to rise by 0.2°C per decade, which will lead to temperatures 1.8-4.0°C higher than the current level by 2100 [23]. HS is therefore becoming an agricultural problem in many areas in the world. HS generally impairs photosynthetic activity and reduces water content caused by heat and has negative effects on cell division and growth of crops. Thus, HS is a major environmental factor limiting crop productivity, and identifying and breeding the heat-tolerant cultivars of crops are a valuable way to protect food production and ensure crop safety [24,25]. For example, the heat-tolerant cultivars have been identified in some major crops, including rice [26], maize [27], and wheat [28]; however, heat-tolerant cultivars of horticultural crops are seldom reported. In the current study, we found a putative heat-tolerant jujube (Ziziphus jujuba Mill.) cultivar ("Hqing1-HR") that can survive under serious HS (45°C). To our knowledge, this is the first report of a heattolerant cultivar of jujube. "Hqing1-HR" could be used to breed more heat-tolerant lines in the future.
Under high-temperature conditions, plants exhibit shortterm avoidance or acclimation mechanisms such as transpirational cooling and stomatal closure [29]. We observed no macroscopic phenotypic differences, such as wilting, leaf curl, or yellowing, in "Hqing1-HR" jujube seedlings under different durations of high-temperature stress ( Figure 1). However, scanning electron microscopy of leaves revealed that stomatal density and opening rate of leaves were significantly affected by high-temperature stress, showing a trend of rapid increase and then slow decrease with the extension of high temperature stress duration. Similar results have been reported in annual plants, such as soybean [30] and rice [31]. Stomatal development is very sensitive to fluctuations in environmental conditions such as temperature, osmotic stress, and carbon dioxide concentration [32]. HS affects the expression of HSP90 [33,34], MUTE [35], SBH1 [36], AGL16 [37], and other genes that are considered regulators of stomatal differentiation by orchestrating the transcriptional network controlling symmetric divisions. In the current study, multiple upregulated DEGs have the capacities to produce HSPs, including HSP17, HSP18, HSP21, HSP 22, HSP 23, HSP26, HSP70, HSP83, HSP90, HSF30, and HSC-2, suggesting HSPs might be molecular chaperones which prevent the formation of nonspecific and harmful protein aggregates and assist proteins in the acquisition of their native structures.     International Journal of Genomics The physiological effects of HS on plants have been extensively reported, but our understanding of the underlying molecular mechanisms remains limited. Expression levels of multiple genes are affected by HS; thus, RNA-seq analysis, which provides precise information on the transcriptomic changes occurring in response to abiotic stress, including HS, is a suitable method for elucidating these mechanisms. For example, transcriptome profiling of rice [26], barley [38], maize [39,40], Brachypodium distachyon [41], and Vitis vinifera (grape) [7,42] in response to HS has generated useful clues associated with the molecular mechanism of the response to HS. In jujube species, RNAseq experiments have also been performed to explore the transcriptomic changes that occur in response to abiotic stress, including drought stress [43] and alkalinity stress [44]. In the current study, we performed transcriptomic analysis for the jujube response to HS using RNA-seq experiments. These indicated that HS changes global expression levels of multiple genes, and we found 1,642, 4,080, 5,160, and 2,119 DEGs in HR1 vs. HR0, HR3 vs. HR0, HR5 vs. HR0, and HR7 vs. HR0 comparisons, respectively. Moreover, functional analyses indicated that a considerable number of DEGs were enriched in terms associated with photosynthesis metabolism, suggesting that "Hqing1-HR" might be tolerant to HS by upregulating or lowering the expression levels of these genes.

Conclusions
In this study, high temperature did not damage to the macroscopic phenotype of "Hqing1-HR." However, stomatal density and opening rate were significantly affected by high-temperature stress. We conducted the transcriptome analysis of leaves and characterization of transcripts related to high-temperature stress during the seedling stage in jujube using a next-generation sequencing approach. A total of 6,606 DEGs were identified in "Hqing1-HR" under heat stress compared with the control treatment, and 1,642, 4,080, 5,160, and 2,119 DEGs were found in HR1 vs. HR0, HR3 vs. HR0, HR5 vs. HR0, and HR7 vs. HR0, respectively. The GO enrichment analysis showed that a series of biological processes related to stress response, photosynthesis, and metabolism were enriched during high-temperature stress, suggesting that down-or upregulation of genes in these processes may play an important role in the response to HS. These results contributed to our understanding of the molecular mechanism of jujube responses to high-temperature stress.

Data Availability
The RNA-seq data has been deposited in NCBI Gene Expression Omnibus (GEO) under accession code GSE136047. Table S1. The DEGs were detected in the comparisons of HR1 vs. HR0. Table S2. The DEGs were detected in the comparisons of HR3 vs. HR0. Table S3. The DEGs were detected in the comparisons of HR5 vs. HR0. Table S4. The DEGs were detected in the comparisons of HR7 vs. HR0. Table  S5. Functional analysis of the DEGs from HR1 vs. HR0. Table S6. Functional analysis of the DEGs from HR3 vs. HR0. Table S7. Functional analysis of the DEGs from HR5 vs. HR0. Table S8. Functional analysis of the DEGs from HR7 vs. HR0. (Supplementary Materials)