Transcriptome Analysis of a Jujube (Ziziphus Jujuba Mill.) Cultivar Response to Heat Stress

Background: Heat stress (HS) is a common stress and inuences the growth and reproduction of plant species. We found and bred a putative heat-resistant jujube (Ziziphus jujuba Mill.) cultivar (JHR17) in previous study. Results: In the current study, we made the seedlings of ‘JHR17’ cultivar to be under HS (45°C) for 0, 1, 3, 5 and 7 days, respectively, and the leaf samples (HR0, HR1, HR3, HR5 and HR7) were collected accordingly. Fifteen cDNA libraries from ‘JHR17’ leaves were built with a transcriptome assay. The RNA sequencing (RNA-seq) and transcriptome comparisons were performed, and the results indicated that 1642, 4080, 5160 and 2119 differentially expressed genes (DEGs) were identied in HR1 vs. HR0, HR3 vs. HR0, HR5 vs. HR0 and HR7 vs. HR0, respectively. Gene Ontology (GO) analyses of the DEGs from these comparisons were implemented. Conclusion: It revealed that a series of biological processes, involved in stress response, photosynthesis and metabolism, were enriched successfully, suggesting that lowering or up-regulating these genes of processes might play important roles in response to HS. This study may contribute to understand the molecular mechanism of ‘JHR17’ cultivar response to HS, and be benecial for developing jujube cultivars to improve heat resistance. sepsis.

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 the highest production and good quality worldwide [13]. However, Xinjiang is particularly vulnerable to climate change, and has experienced signi cant climate warming in the past 40 years [14]. In recent years, the HS dramatically and repeatedly affected the production of quality of jujube fruits in Xinjiang. Therefore, searching and breeding heat-tolerant jujube cultivars might be one of feasible and important strategies to protect production and quality of jujube fruits.
Turpan, one of 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 putative heat-resistant jujube cultivar ('JHR17') in an orchard of Turpan by chance [16], and bred it in our laboratory successfully. In the current study, to investigate the transcriptomic change of 'JHR17' response to HS, we made its seedlings to be under 45 °C. In the 0, 1, 3, 5 and 7 days after HS-treatment, we checked the phenotypic and physiological features, and collected the samples for RNA-seq experiments. This study may provide a new insight into the transcriptional alterations in heat-resistant jujube cultivar responding to HS.

Phenotypic of jujube seedlings post HS
We found a putative heat-resistant jujube cultivar ('JHR17') in Turpan city in China, and inbred the cultivar in our laboratory. To obtain an overview of heat-tolerance phenotype of jujube cultivar, seedlings with 14 true leaves were subjected to treatment with HS (45 °C). 0 (control) 1, 3, 5 and 7 days, after the treatment, the leaves of all seedlings did not become withered (Fig. 1A), suggesting this jujube cultivar might be of heat-tolerance.
Stoma is an important channel for gas and water exchange between plants and the atmosphere, which can make adaptive adjustment under various stress conditions. The stomatal density and stomatal opening rate of leaves from each group were assessed. It indicated that the stomatal density and stomatal opening rate post heat treatment were signi cantly increased, and they showed a trend of rst increase and then decrease with the extension of heat treatment (Fig. 1B, Table 1). It suggested that the 'JHR17' could reduce the damage by passively changing stomatal density and opening rate. Through Illumina HiSeq X-ten platform, we generated over 0.402 billion pair-end reads, corresponding to an average of 26.8 million sequence reads per sample. Using HISAT2 [17]with default settings, about 68.9% clean reads were mapped to the jujube reference genome [11].
In order to understand the spatiotemporal expression patterns of all samples, principal component analysis (PCA) was performed. The three samples in the same time point could form independent clusters ( Fig. 2A). Moreover, Pearson correlation analysis for all pairs of RNA-seq samples was performed, demonstrating similar results (Fig. 2B), indicating that gene expression in three replications of every sample was homogeneous ( Fig. 2A). This suggests that the replicated samples produced data that are acceptable for further analyses.
Exploration of differentially expressed genes (DEGs) down-regulated gene were discovered. 1019 up-regulated and 1100 down-regulated genes were identi ed in HR7 vs. HR0. The numbers of the upregulated and downregulated DEGs were similar in each comparison, indicating that heat stress had effect on promoting and inhibiting the transcription of numerous genes. Moreover, there are 717 common DEGs among four comparisons above (Fig. 3B). In order to verify the validity of these results, ve genes from each list, with altered expression levels and FCs, were chosen as representatives to quantify their expression by RT-qPCR. Results shown in Fig. 5A con rm the robustness of global expression data obtained for both upregulated and downregulated genes (upper and lower panel, respectively).

The molecular response to heat stress
To identify the pathways in which the DEGs were mainly involved, Gene Ontology (GO) enrichment analysis was conducted. It showed that 36, 58, 65 and 37 GO terms were identi ed in HR1 vs. HR0, HR3 vs. HR0, HR5 vs. HR0 and HR7 vs. HR0, respectively (P < 0.01).
Although the leaves of all seedlings did not become withered under high temperature, multiple DEGs were enriched in "response to stress" (GO: 0006950) and "response to heat" (GO: 0009408) terms for all four comparisons ( Fig. 4 and Table S2), indicating the 'JHR17' could be sensitive to HS normally and might establish a new steady-state balance of biological processes that enable the organism to function, survive at a higher temperature (45 °C) perfectly. We analyzed these DEGs enriched in "response to stress" and "response to heat" terms, it indicated that the expression level of multiple DEGs associated with response to heat stress were obviously up-regulated after HS (Fig. 4A).

Photosynthesis were affected by heat stress
Photosynthesis occurs in chloroplast, and is a sensitive biological process to high temperature in plants.
In order to explore how the HS affected the normal physiology of chloroplasts and photosynthesis, the DEGs were enriched in "chloroplast organization" (GO: 0009399), "chloroplast RNA processing" (GO: 0031425), "photosynthesis, light harvesting" (GO: 0009765) and "photosynthesis" (GO: 0015979) terms were checked and analyzed. To surprise, most of the DEGs enriched in the four terms above were upregulated by the HS. Moreover, this phenomenon was consistent with the results of the qRT-PCR experiments (Fig. 5). It suggested that the HS might not disrupt the physiology of chloroplasts and photosynthesis and promote the photosynthesis of 'JHR17'.
Metabolism were affected by heat stress HS always globally affects the metabolism of the plants. The "myo-inositol hexakisphosphate biosynthetic process" (GO: 0010264) was the only common term associated with metabolism identi ed in HR1 vs. HR0, HR3 vs. HR0, HR5 vs. HR0 and HR7 vs. HR0. However, multiple speci c terms associated with metabolism in four comparisons.

Validation of RNA-seq by qRT-PCR
In order to verify the reliability of the transcriptome data, six DEGs were randomly selected for expression analysis by qRT-PCR experiments (Fig. 5). The expression patterns shown in the qRT-PCR results (Fig. 5B) were consistent with RNA-seq results (Fig. 5A), 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 [19], so the HS has becoming an agricultural problem in many areas in the world. HS generally impairs photosynthetic activity and the reduced water content caused by heat, and has negative effects on cell division and growth of crops. Thus, the HS is a major environmental factor limiting crop productivity and searching and breeding the heat-tolerant cultivars of crops is a suitable way to protect food production and ensure crop safety [20,21]. For example, the heatresistant cultivars were found in some major crops, including rice [22], maize [23], and wheat [24], but the heat-resistant cultivars of horticultural crops were seldom reported. In the current study, we found a putative heat-resistant jujube (Ziziphus jujuba Mill.) cultivar ('JHR17'), which can survive under serious HS (45 °C). To our knowledge, this is the rst report of heat-resistant cultivar in jujube species. Maybe, the 'JHR17' could be used to breed more good lines in the future.
Under high temperature conditions, plants exhibit short-term avoidance or acclimation mechanisms such as transpirational cooling, stomatal closure, and so on [25]. In our study, we found that there were no macroscopic phenotypic differences, such as wilting, leaf curl and yellowing in 'JHR17' jujube seedlings under different periods of high temperature stress (Fig. 1). However, the observations through scanning electron microscope showed that the stomatal density and opening rate of leaves were signi cantly affected by high temperature stress, which showed a trend of rapid increase and then slow decrease with the extension of high temperature stress time. Similar results has been reported in annual plants, such as soybean [26] and rice [27] yet. Stomatal development is very sensitive to environmental uctuations, such as temperature, osmotic stress, and carbon dioxide concentration [28]. Heat stress affects the expression of HSP90 [29,30], MUTE [31], SBH1 [32], AGL16 [33], et al., these genes are considered as regulators of stomatal differentiation by orchestrating the transcriptional network controlling the symmetric divisions.
The physiological effects of HS on plants have been extensively reported, but the understanding of underlying molecular mechanism remains limited. The expression levels of multiple genes would be affected by the HS, thus the RNA-seq analysis, which provides precise information on the transcriptomic changes that occur in response to abiotic stress including HS, is a suitable method. For examples, transcriptome pro ling of rice [22], barley [34], maize [35,36], Brachypodium distachyon [37] and Vitis vinifera (grape) [7,38] in response to HS has been performed, and some useful clues associated with molecular mechanism response to HS were obtained from these RNA-seq data. In jujube species, the RNA-seq experiments were also performed to explore the transcriptomic changes that occur in response to abiotic stress, including drought stress [39] and alkaline stress [40]. In the current study, transcriptomic analysis for jujube response to HS was rst carried out by RNA-seq experiments. It indicated that HS globally changed the expression levels of multiple genes, and 1642, 4080, 5160 and 2119 DEGs were found in HR1 vs. HR0, HR3 vs. HR0, HR5 vs. HR0 and HR7 vs. HR0, respectively. Moreover, functional analyses indicated that a huge number of DEGs were enriched in the terms associated with photosynthesis metabolism, suggesting that the 'JHR17' cultivar might be resistant to HS by upregulating or lowering the expression levels of these genes.

Conclusions
In this study, high temperature did no damage on the macroscopic phenotypic of 'JHR17'. However, the stomatal density and opening rate were signi cantly affected by high temperature stress. Furthermore, we conducted the transcriptome analysis of leaves and characterization of transcripts related to high temperature stress during the seedling stage in jujube using next generation sequencing approach. A total of 6606 DEGs were identi ed in 'JHR17' heat stress compared with the control treatment, and 1642, 4080, 5160 and 2119 DEGs were found in HR1 vs. HR0, HR3 vs. HR0, HR5 vs. HR0 and HR7 vs. HR0, respectively. 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 lowering or upregulation of genes in these processes may play an important role in response to heat stress. These results may contribute to understand the molecular mechanism of jujube response to high temperature stress.

Plant materials, heat treatment and sample collection
The green cuttings of the cultivar ('JHR17') were collected from the jujube orchard of Turpan in the Xinjiang of China. The green cuttings were grown in the greenhouse conditions with an automatic spray system (20 ~ 35℃ with 90% humidity).When the cutting seedlings had 7 ~ 9 true leaves, they were transferred to a controlled growth chamber with a light/dark regime of 14/10 h at 30/20℃, 80% relative humidity, and light intensity of 300 µmol m − 2 s − 1 of photosynthetically active radiation.
Seedlings with 14 true leaves were then cultured in the same chamber with the same conditions except with the temperature at 45℃. After 0 (control), 1, 3, 5, and 7d of heat treatment, the 10th true leaves were collected from all three different samples. The samples were immersed in liquid nitrogen and stored at -80℃ for subsequent experiments.

Determination of phenotypic of jujube leaves
The leaves of the same part were rinsed and xed with FAA solution(70% ethanol)at 4℃. leaves were freeze-dried after dehydration by alcohol gradient series to critical drying point,then stuck the sample on table with conductive tapes, and coating the sample with a Pt lm using an ion sputtering instrument (Hitachi E-1045) and nally the SUPRA 55VP scanning electron microscope was used for observation of leaves at an accelerating voltage of 2.00 kV (Zeiss, German).

RNA Extraction, cDNA Library Construction and Illumina Sequencing
Total RNAs was extracted from 15 samples that represent for three biological replicates of JHR17 cultivar at ve treatment stages (0, 1, 3, 5, 7 days), using the RNAprep Pure Plant Kit (Tiangen, Beijing, China) according to the instructions of manufacturer. The extracted RNA was treated with DNase I (Promega, Madison, WI, USA) to remove DNA. RNA quality and quantity were assessed using Nanodrop 2000 Spectrophotometer (Thermo Fisher Scientic, Wilmington, USA) and an Agilent Bioanalyzer 2100 System (Agilent Technologies, CA, USA), respectively. The integrity of RNA was con rmed by 1% agarose gel electrophoresis.
RNA samples were pooled from equal amounts of RNA from three independent individuals, and then were used for libraries preparation and sequencing. The resulting libraries were sequenced using an Illumina HiSeq X-ten platform with paired-end 150 bp reads.
RNA-seq read processing and assembly Raw reads were generated by the Illumina HiSeq X-ten genome analyzer, and then were analyzed by Fast Q to assess the base quality and 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 [41] and then the remaining cleaned reads were mapped to the jujube reference genome sequences [11] using HISAT2 [17] with default settings.
Bioinformatic analysis FPKM (fragments per kilobase per million mapped reads), was used to evaluate the expression level of genes. To measure the FPKM value and screen out the differentially expressed genes (DEGs), we applied the software edge R [18]. The genes with RPKM < 0.1 in every sample were removed before analysis. To determine whether a gene was differentially expressed, we analyzed the results based on the fold change (FC) (FC ≥ 2 or ≤ 0.5) and false discovery rate (FDR) (FDR < 0.01).
To predict the gene function and calculate the functional category distribution frequency, Gene Ontology (GO) analyses were employed using DAVID bioinformatics resources [42]. Venn diagrams were built using an online available tool (http://bioinfogp.cnb.csic.es/ tools/venny/).

Validation of RNA-seq by qRT-PCR
In this study, to elucidate the validity of the RNA-seq data, quantitative real-time PCR (qRT-PCR) experiments were performed for some randomly selected DEGs. The designed primers were presented in Table 2. The same RNA samples for RNA-seq were used for qPCR. In each pooled sample, l µg of RNA was reversely transcribed using the Prime Script™ RT Reagent Kit (Takara, Dalian, China) according to the instructions of manufacturer. qPCR was performed on the Bio-Rad S1000 with Bestar SYBR Green RT-PCR Master Mix (DBI Bioscience, Shanghai, China). The PCR conditions were as follows: denaturing at 95 °C for 8 min, 38 cycles of denaturing at 95 °C for 15 s, annealing and extension at 60 °C for 1 min. Relative gene expression was calculated using the Livak and Schmittgen 2 −ΔΔCt method [43], normalized with the reference gene ZjH3 of jujube. PCR ampli cations were performed in triplicate for each sample.

Statistical analysis
All values of the all data were presented as mean ± standard deviation (SD). In order to determine the signi cance of differences between means, the Student's t-test (paired) was implemented, and a value P < 0.05 was considered to be statistically signi cant. Declarations opening rate of leaves from each group were assessed. It indicated that the stomatal density and stomatal opening rate post heat treatment were signi cantly increased, and they showed a trend of rst increase and then decrease with the extension of heat treatment