Comparative Transcriptomic Analyses of Different Jujube Cultivars Reveal the Co-Regulation of Multiple Pathways during Fruit Cracking

Fruit cracking is a common physiological disorder in many fruit species. Jujube (Ziziphus jujuba Mill.) is an economically valuable fruit in which fruit cracking seriously affects fruit yield and quality and causes significant economic losses. To elucidate cracking-related molecular mechanisms, the cracking-susceptible cultivars ‘Cuizaohong’ and ‘Jinsixiaozao’ and the cracking-resistant cultivar ‘Muzao’ were selected, and comparative transcriptome analyses of cracking and non-cracking ‘Cuizaohong’ (CC and NC), cracking and non-cracking ‘Jinsixiaozao’ (CJ and NJ), and non-cracking ‘Muzao’ (NM) were conducted. A total of 131 differentially expressed genes (DEGs) were common to the CC vs. NC and CJ vs. NJ comparisons. To avoid passive processes after fruit cracking, we also mainly focused on the 225 gradually downregulated DEGs in the CJ, NJ, and NM samples. The functional annotation of the candidate DEGs revealed that 61 genes related to calcium, the cell wall, the cuticle structure, hormone metabolism, starch/sucrose metabolism, transcription factors, and water transport were highly expressed in cracking fruits. We propose that expression-level changes in these genes might increase the turgor pressure and weaken mechanical properties, ultimately leading to jujube fruit cracking. These results may serve as a rich genetic resource for future investigations on fruit cracking mechanisms in jujube and in other fruit species.


Introduction
Chinese jujube (Ziziphus jujuba Mill.), which belongs to the family Rhamnaceae, is an economically, ecologically, and culturally valuable fruit crop species in China. Jujube fruits are a rich source of energy (i.e., calories) as well as nutritional and medicinal compounds, including carbohydrates, minerals, vitamin C, and cAMP [1]. Jujube is a major fruit tree species with considerable economic benefits for farmers in China. More than eight million tons of jujube fruits were produced in 2016, and jujube fruit cultivation serves as the primary source of income for nearly 20 million people in China [2]. However, fruit cracking has become a serious problem because it adversely affects the appearance of jujube fruits and increases the chances of microbial infections, which lead to significant losses in fruit yield and quality. Regarding jujube production, yield losses due to fruit cracking range from approximately 30% in normal years to more than 90% in severe years [3].
DEGs were identified. Functional categorization and co-expression analysis of cracking-100 related genes were conducted to reveal the various metabolic pathways that are involved 101 in fruit cracking. Our study provides candidate genes for genetically enhancing cracking 102 tolerance to increase fruit quality and yield and also lays the foundation for future studies 103 on the molecular mechanisms underlying fruit cracking in jujube and in other fruit spe-104 cies. 105 106 107 The three jujube cultivars 'Cuizaohong', 'Jinsixiaozao', and 'Muzao' that were used 108 in our experiments were provided by the National Foundation for Improved Cultivars of 109 Chinese Jujube, Cangzhou County, Hebei Province, China. The trees had been subjected 110 to the same orchard management strategy. A total of thirty uniform intact half-red fruits 111 of each jujube cultivar were collected with three biological replicates, in which ten indi-112 viduals per replicate were used for fruit cracking and fruit size measurements. The half-113 red fruit from the CC, NC, CJ, NJ, and NM cultivars were collected for RNA-seq, and the 114 cracking and non-cracking fruit samples from the same cultivar were collected from the 115 same tree for each replicate ( Figure 1A). The sampling data was September 10, 2017. The 116 collected fruit were divided into three biological replicates, with each containing at least 117 three fruits. Fruits without stones were frozen immediately in liquid nitrogen and were 118 stored at −80 °C for RNA extraction. 119 120 Figure 1. The phenotypes of the jujube at the half-red period (A) and evaluation of cracking re-121 sistance (B) and fruit size (C) of three jujube cultivars. CC: cracking 'Cuizaohong'; NC: non-cracking 122 'Cuizaohong'; CJ: cracking 'Jinsixiaozao'; NJ: non-cracking 'Jinsixiaozao'; NM: non-cracking 'Mu-123 zao', the same as below. Error bars represent standard errors from three biological replicates. Dif-124 ferent lowercase letters indicate a significant difference between genotypes (p < 0.05). 125 126 The fruit cracking (cracking rate and cracking index) and fruit size (fruit transverse 127 diameter, fruit longitudinal diameter, and fruit weight) of the three jujube cultivars were 128 measured for the three-year period of 2015-2017.Jujube fruit cracking was induced by 129 Figure 1. The phenotypes of the jujube at the half-red period (A) and evaluation of cracking resistance (B) and fruit size (C) of three jujube cultivars. CC: cracking 'Cuizaohong'; NC: non-cracking 'Cuizaohong'; CJ: cracking 'Jinsixiaozao'; NJ: non-cracking 'Jinsixiaozao'; NM: non-cracking 'Muzao', the same as below. Error bars represent standard errors from three biological replicates. Different lowercase letters indicate a significant difference between genotypes (p < 0.05).

Measurements of Fruit Cracking and Fruit Size
The fruit cracking (cracking rate and cracking index) and fruit size (fruit transverse diameter, fruit longitudinal diameter, and fruit weight) of the three jujube cultivars were measured for the three-year period of 2015-2017. Jujube fruit cracking was induced by means of a natural immersion method, and the cracking rate was measured at 48 h, as described by Du et al. [20]. The cracking index was calculated according to the grading system for jujube fruit cracking [21]. The transverse and longitudinal diameters (maximum value) of fruits were measured using a vernier caliper, and the fruit weight was measured with the help of a digital electric balance.

RNA Extraction and Sequencing
Total RNA was isolated using TRIzol (Invitrogen, Carlsbad, CA, USA) in accordance with the manufacturer's instructions. The RNA concentration was quantified using a NanoDrop 2000 (Thermo Fisher Scientific, Wilmington, DE, USA), and RNA integrity was evaluated using an RNA Nano 6000 Assay Kit on the Agilent Bioanalyzer 2100 system (Agilent Technologies, Palo Alto, CA, USA). A total amount of 1 µg RNA per sample was used as input material for the RNA sample preparations. Sequencing libraries were generated using the NEB Next Ultra small RNA Sample Library Prep Kit for Illumina (NEB, Ipswich, MA, USA) following the manufacturer's recommendations, and index codes were added to attribute sequences to each sample. The clustering of the index-coded samples was performed on a cBot Cluster Generation System using a TruSeq PE Cluster Kit v4-cBot-HS (Illumina, San Diego, CA, USA) in accordance with the manufacturer's instructions. After cluster generation, the library preparations were sequenced on an Illumina HiSeq Xten platform, and paired-end reads were generated. The isolation of the mRNA, fragment interruption, cDNA synthesis, adapter addition, PCR amplification, and RNA-seq were performed by Biomarker Technologies Co., Ltd. (Beijing, China).

Sequence Data Analysis
Raw data (raw reads) in the fastq format were first processed using in-house perl scripts. In this step, clean data (clean reads) were obtained by removing reads containing adapters and poly-Ns as well as low quality reads from the raw data. For the clean data, the percentages of the bases whose Phred values were more than 20 (Q20) and 30 (Q30) [22], the content of guanine and cytosine for the library in one tissue (GC content), and the frequency of reads with the same sequences (sequence duplication level) were calculated. The remaining high-quality clean sequencing reads were mapped onto the 'Dongzao' reference genome sequence (version 1.1) [23] using Tophat2 software [24].

Differential Expression Analysis
The jujube gene expression levels were estimated using fragments per kilobase of transcript per million fragments mapped (FPKM) [25]. Differential expression analyses of two samples were performed using edgeR software [26]. A p-value < 0.05 and an absolute value of log 2 fold change (FC) ≥ 1 were used as thresholds to determine the significant differences in gene expression levels.
The DEGs were clustered using the Short Time-series Expression Miner (STEM) software (http://www.sb.cs.cmu.edu/stem). The DEGs with similar expression patterns were included in the same profile. The profiles with p < 0.05 were identified as significantly enriched modules.

Functional Annotation, Enrichment Analysis and Transcription Factor Identification
Gene function was annotated based on the following databases: the NCBI nonredundant protein database (Nr), the Clusters of Orthologous Groups database (COG), Swiss Institute of Bioinformatics databases (Swiss-Prot), and the evolutionary genealogy of genes: Nonsupervised Orthologous Groups database (eggNOG) using BLAST (E value ≤ 1 × 10 −5 ). Gene Ontology functional enrichment and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were carried out using the GOseq R packages [27] and KOBAS software [28]. DEGs were significantly enriched in GO terms, and KEGG pathways were determined at p-values ≤ 0.05. Software iTAK was used to identify transcription factors (TFs) [29].

Co-Expression Analysis
The co-expression between DEGs was analyzed using SPSS (version 19) by calculating the Spearman's correlation coefficients of their FPKM values. Then, the co-expression of these DEGs was visualized using Cytoscape [30].

Validation by qRT-PCR Analysis
Total RNA was isolated from fruit tissue as described above. The first-strand cDNA was synthesized from 500 ng total RNA using a FastQuant RT Kit (with gDNase) (Tiangen, Beijing, China) in accordance with the manufacturer's instructions. The specific primers were designed using Beacon Designer 7.9 for qRT-PCR analysis. The ZjACT gene was used as an internal control (housekeeping gene) [31]. All of the primer sequences that were used for qRT-PCR are provided in Table S1. qRT-PCR was carried out on an ABI 7500 Fast Real-Time PCR System (Applied Biosystems, USA) using SYBR green (Aidlab Biotechnology, Beijing, China) in a 20 µL reaction mixture containing 1 µL of diluted cDNA, 400 nM of each primer, and 10 µL 2 × SYBR Green PCR Master Mix. The PCR reaction protocol was 94 • C for 3 min, 40 cycles of 94 • C for 20 s, 55 • C for 20 s, and 72 • C for 30 s. Relative gene expression levels were determined using the 2 −∆∆Ct method.

Statistical Analysis
SPSS 19.0 and EXCEL 2013 were used for statistical analysis. Analysis of statistical significance was performed according to one-way analysis of variance (ANOVA) with Duncan. Each experiment was repeated three times independently, and p < 0.05 was considered to represent significance.

Analysis of Fruit Characteristics in Three Jujube Cultivars
The cracking resistance of three jujube cultivars were measured for three continuous years. Two cultivars, 'Cuizaohong' and 'Jinsixiaozao', exhibited a high cracking rate and cracking index, but no 'Muzao' fruit were cracked ( Figure 1B). Based on the cracking rate and index, 'Cuizaohong' and 'Jinsixiaozao' were cracking-susceptible cultivars, and 'Muzao' was a strong cracking-resistant cultivar. In terms of fruit size, there was no significant difference in terms of the fruit transverse diameter, fruit longitudinal diameter, and fruit weight between 'Jinsixiaozao' and 'Muzao', which demonstrated contrasting crack resistance ( Figure 1C), which indicated that these two cultivars are ideal materials for investigating the fruit cracking mechanisms in jujube [19].

Transcriptome Sequencing Results
Three independent replicates of each fruit sample (CC, NC, CJ, NJ, and NM) were subjected to RNA-seq analyses to identify the potential molecular mechanisms that were responsible for jujube fruit cracking. A total of 26.36, 26.32, 29.15, 27.45, and 28.83 million raw reads were obtained from the CC, NC, CJ, NJ, and NM transcriptome libraries, respectively (Table 1). After quality control, more than 84.21% clean reads (the percentages of Q30 and GC being 86.48-91.50% and 44.26-46.12%, respectively) were obtained and were used for downstream analyses (Table S2). The alignment results showed that most of the clean reads (70.30-73.59%) from the 15 libraries could be mapped to the reference genome (Table S3). On average, approximately 32.13 (62.23%), 32.59 (59.93%), and 35.34 (62.54%) million reads were uniquely mapped on the reference genome using Tophat2 for 'Cuizaohong', 'Jinsixiaozao', and 'Muzao', respectively. From the saturation curve, even the genes with low expression levels (FPKM < 1) were saturated ( Figure S1). This indicated that the number of obtained reads was sufficient to cover most of the expressed genes. The gene coverage analysis showed no bias in the curves. The results showed that the reads were evenly distributed among the genes in the reference genome ( Figure S2). Overall, the data indicated that our Illumina sequencing was of high quality and could be used for further analyses.  Figure 2A). A Venn diagram indicated the DEG distribution among the three comparisons ( Figure 2B). A total of 21 downregulated genes overlapped in CC vs. NC, CJ vs. NJ, and NJ vs. NM. An analysis of the DEGs was performed between the two cracking-susceptible jujube cultivars ('Cuizaohong' and 'Jinsixiaozao'), and 131 DEGs showed similar expression patterns in the cracking and non-cracking fruits, reflecting the common cracking-related genes in both cultivars. When these two cultivars were exposed to fruit cracking, the number of downregulated genes (125 DEGs) was greater than the number of upregulated genes (six DEGs), and 'Cuizaohong' had more DEGs than 'Jinsixiaozao', indicating more complex cracking response pathways in 'Cuizaohong'. The gene coverage analysis showed no bias in the curves. The results showed that the 226 reads were evenly distributed among the genes in the reference genome ( Figure S2). Over-227 all, the data indicated that our Illumina sequencing was of high quality and could be used 228 for further analyses. 229  NM. An analysis of the DEGs was performed 239 between the two cracking-susceptible jujube cultivars ('Cuizaohong' and 'Jinsixiaozao'), 240 and 131 DEGs showed similar expression patterns in the cracking and non-cracking fruits, 241 reflecting the common cracking-related genes in both cultivars. When these two cultivars 242 were exposed to fruit cracking, the number of downregulated genes (125 DEGs) was 243 greater than the number of upregulated genes (six DEGs), and 'Cuizaohong' had more 244 DEGs than 'Jinsixiaozao', indicating more complex cracking response pathways in 245 'Cuizaohong'. In 'Cuizaohong', the log2FC of the DEGs ranged from −11.27 to 2.59, while in 249 'Jinsixiaozao', this parameter fluctuated between −4.68 and 4.38 (Table S4). In 250 'Cuizaohong', gene5475 was the most downregulated gene (log2FC = −11.27). It was anno-251 tated as encoding an unknown protein. The gene30138 showed the greatest increase in 252 expression (log2FC = 2.59) and was also annotated as encoding an unknown protein. In 253 'Jinsixiaozao', gene5490 exhibited the greatest expression level (log2FC = −4.68) and en-254 coded linoleate 13S-lipoxygenase 2-1, chloroplastic. The expression levels of gene28904 255 displayed the most dramatic repression (log2FC = 4.38) and was also annotated as 256 In 'Cuizaohong', the log 2 FC of the DEGs ranged from −11.27 to 2.59, while in 'Jinsixiaozao', this parameter fluctuated between −4.68 and 4.38 (Table S4). In 'Cuizaohong', gene5475 was the most downregulated gene (log 2 FC = −11.27). It was annotated as encoding an unknown protein. The gene30138 showed the greatest increase in expression (log 2 FC = 2.59) and was also annotated as encoding an unknown protein. In 'Jinsixiaozao', gene5490 exhibited the greatest expression level (log 2 FC = −4.68) and encoded linoleate 13S-lipoxygenase 2-1, chloroplastic. The expression levels of gene28904 displayed the most dramatic repression (log 2 FC = 4.38) and was also annotated as encoding an unknown protein. Because the functions of the most up-and downregulated genes in each library are unknown, further analyses to identify their functions are needed.

Functional Annotation of DEGs
To determine the functions of the DEGs, their sequences were used as queries to screen six databases. In total, 710 (98.47%), 444 (97.58%), and 1621 (96.26%) DEGs in the CC vs.  (Table S5). A GO analysis was performed to classify the DEGs according to their functions. Among the DEGs, 434 (37.94%), 246 (52.09%), and 930 (41.50%) in the CC vs. NC, CJ vs. NJ, and NJ vs. NM comparisons, respectively, were assigned GO terms (Table S5). The most common biological process GO term among the DEGs that were detected in the three comparisons was the 'oxidation-reduction process'. In the molecular function category, 'lupeol synthase activity' and 'hydroxymethylglutaryl-CoA reductase (NADPH) activity' were significantly enriched among the DEGs that were detected in the comparisons of cracking and non-cracking fruits (CC vs. NC and CJ vs. NJ). However, regarding the NJ vs. NM comparison, the most enriched molecular function GO term was 'oxidoreductase activity'. For all of the comparisons, the significantly enriched cellular component GO terms were 'cell wall' and 'plant-type cell wall' (Figure 3A-C).
encoding an unknown protein. Because the functions of the most up-and downregulated 257 genes in each library are unknown, further analyses to identify their functions are needed. 258

259
To determine the functions of the DEGs, their sequences were used as queries to 260 screen six databases. In total, 710 (98.47%), 444 (97.58%), and 1,621 (96.26%) DEGs in the 261 CC vs. NC, CJ vs. NJ, and NJ vs. NM comparisons, respectively, matched sequences in at 262 least one database (Table S5). A GO analysis was performed to classify the DEGs accord-263 ing to their functions. Among the DEGs, 434 (37.94%), 246 (52.09%), and 930 (41.50%) in 264 the CC vs. NC, CJ vs. NJ, and NJ vs. NM comparisons, respectively, were assigned GO 265 terms (Table S5). The most common biological process GO term among the DEGs that 266 were detected in the three comparisons was the 'oxidation-reduction process'. In the mo-267 lecular function category, 'lupeol synthase activity' and 'hydroxymethylglutaryl-CoA re-268 ductase (NADPH) activity' were significantly enriched among the DEGs that were de-269 tected in the comparisons of cracking and non-cracking fruits (CC vs. NC and CJ vs. NJ). 270 However, regarding the NJ vs. NM comparison, the most enriched molecular function GO 271 term was 'oxidoreductase activity'. For all of the comparisons, the significantly enriched 272 cellular component GO terms were 'cell wall' and 'plant-type cell wall' (Figure 3A-C). To identify the biological pathways that are activated in the selected jujube fruits, we 283 mapped the annotated sequences to the reference pathways in the KEGG database. In this 284 analysis, 80, 62, and 104 significantly enriched functional categories were identified for 285 the DEGs in the CC vs. NC, CJ vs. NJ, and NJ vs. NM comparisons, respectively (Table 286  S6). Two pathways ('Amino sugar and nucleotide sugar metabolism' and 287 To identify the biological pathways that are activated in the selected jujube fruits, we mapped the annotated sequences to the reference pathways in the KEGG database. In this analysis, 80, 62, and 104 significantly enriched functional categories were identified for the DEGs in the CC vs. NC, CJ vs. NJ, and NJ vs. NM comparisons, respectively (Table S6). Two pathways ('Amino sugar and nucleotide sugar metabolism' and 'phenylpropanoid biosynthesis') were common among the top 20 KEGG pathways for the DEGs detected in the three comparisons ( Figure 3D-F). For the two cracking-susceptible jujube cultivars, 'flavonoid biosynthesis' was significantly enriched among the DEGs that were revealed by the comparison of the cracking and non-cracking fruits (CC vs. NC and CJ vs. NJ, Figure 3D,E). However, in the NJ vs. NM comparison, most of the DEGs were associated with 'starch and sucrose metabolism' (Figure 3F).
DEGs detected in the three comparisons ( Figure 3D-F). For the two cracking-susceptible 289 jujube cultivars, 'flavonoid biosynthesis' was significantly enriched among the DEGs that 290 were revealed by the comparison of the cracking and non-cracking fruits (CC vs. NC and 291 CJ vs. NJ, Figure 3D,E). However, in the NJ vs. NM comparison, most of the DEGs were 292 associated with 'starch and sucrose metabolism' (Figure 3F). 293

Identification of Differentially Expressed Transcription Factors
To explore the roles of TFs in fruit cracking, we analyzed the differential expression of the TFs in CC vs. NC, CJ vs. NJ, and NJ vs. NM. A total of 124 TFs from 26 families, such as AP2/ERF (18), MYB (13), WRKY (10), and bZIP (eight), were differentially expressed in all three groups ( Figure S3). Among these TFs, 24 were in CC vs. NC, 44 were in CJ vs. NJ, and 79 were in NJ vs. NM ( Figure 6). The most abundant TFs were AP2/ERF family members in each comparison.

Co-Expression Analysis of Cracking-Related Genes
To explore the DEG functions related to the regulation of jujube fruit cracking, we functionally annotated candidate DEGs, including 131 genes common to the CC vs. NC and CJ vs. NJ comparisons as well as 225 gradually downregulated genes in the CJ, NJ, and NM samples. The DEGs in the Venn diagram and trend analyses were divided into seven groups (Ca, cell wall, cuticle structure, hormone metabolism, starch/sucrose metabolism, TF, and water transport) ( Figure 7A,B). Regarding the CC, NC, CJ, and NJ samples, the expression levels of 33 cracking-related DEGs (one water transport gene, eight cell wall metabolism genes, three starch and sucrose metabolism genes, five cuticle structure genes, two Ca transport genes, two ABA metabolism genes, three indoleacetic acid (IAA) metabolism genes, two jasmonic acid (JA) metabolism genes, and seven TFs) were upregulated in the cracking fruits of the two cracking-susceptible jujube cultivars ( Figure 7A). In the CJ, NJ, and NM samples, the expression levels of 34 cracking-related DEGs (one water transport gene, four cell wall metabolism genes, one cuticle structure gene, four Ca transport genes, five ABA metabolism genes, two IAA metabolism genes, two GA metabolism genes, and fifteen TFs) were gradually downregulated ( Figure 7B). Additionally, one water transport gene (gene6278), one ABA metabolism gene (gene30271), one IAA metabolism gene (gene15774), and three TFs (gene6845, gene12280, and gene30640) were common to the two analyses.

349
To explore the DEG functions related to the regulation of jujube fruit cracking, we 350 functionally annotated candidate DEGs, including 131 genes common to the CC vs. NC 351 and CJ vs. NJ comparisons as well as 225 gradually downregulated genes in the CJ, NJ, 352 and NM samples. The DEGs in the Venn diagram and trend analyses were divided into 353 seven groups (Ca, cell wall, cuticle structure, hormone metabolism, starch/sucrose metab-354 olism, TF, and water transport) ( Figure 7A,B). Regarding the CC, NC, CJ, and NJ samples, 355 the expression levels of 33 cracking-related DEGs (one water transport gene, eight cell 356 wall metabolism genes, three starch and sucrose metabolism genes, five cuticle structure 357 genes, two Ca transport genes, two ABA metabolism genes, three indoleacetic acid (IAA) 358 metabolism genes, two jasmonic acid (JA) metabolism genes, and seven TFs) were upreg-359 ulated in the cracking fruits of the two cracking-susceptible jujube cultivars ( Figure 7A). 360 In the CJ, NJ, and NM samples, the expression levels of 34 cracking-related DEGs (one 361 water transport gene, four cell wall metabolism genes, one cuticle structure gene, four Ca 362 transport genes, five ABA metabolism genes, two IAA metabolism genes, two GA metab-363 olism genes, and fifteen TFs) were gradually downregulated ( Figure 7B). Additionally, 364 one water transport gene (gene6278), one ABA metabolism gene (gene30271), one IAA me-365 tabolism gene (gene15774), and three TFs (gene6845, gene12280, and gene30640) were com-366 mon to the two analyses. The expression levels of the genes related to Ca, cell wall, cuticle structure, hormone metabolism, starch/sucrose metabolism, TF, and water transport in jujube underwent a network analysis to examine the relationship between these processes. The expression levels of one water transport gene, nine cell wall metabolism genes, three starch and sucrose metabolism genes, five cuticle structure genes, six Ca transport genes, fourteen hormone metabolism genes, and seventeen TFs were highly correlated ( Figure 7B). The cell wall metabolism-related, hormone metabolism-related, and TF expression levels were positively correlated with the expression of the analyzed cracking-related genes. Moreover, the expression levels of the water transport-related and cuticle structure-related genes were highly positively correlated with the expression of the starch/sucrose metabolism-related and Ca-related genes, respectively.

qRT-PCR Validation
To validate the transcriptome data, we selected 11 cracking-related genes and evaluated their expression profiles using qRT-PCR. The gene-specific primers used in the qRT-PCR analysis are listed in Table S1. Figure 8 illustrates the expression levels of the DEGs as determined by qRT-PCR and RNA-seq. The qRT-PCR results indicated that the expression profiles of the 11 DEGs were consistent with the corresponding FPKM values derived from the RNA-seq analysis. These genes exhibited a similar expression pattern when both methods were used. The similar results that were obtained from the qRT-PCR and RNA-seq analyses indicated that the RNA-seq data were reproducible and reliable. The expression levels of the genes related to Ca, cell wall, cuticle structure, hormone 376 metabolism, starch/sucrose metabolism, TF, and water transport in jujube underwent a 377 network analysis to examine the relationship between these processes. The expression 378 levels of one water transport gene, nine cell wall metabolism genes, three starch and su-379 crose metabolism genes, five cuticle structure genes, six Ca transport genes, fourteen hor-380 mone metabolism genes, and seventeen TFs were highly correlated ( Figure 7B). The cell 381 wall metabolism-related, hormone metabolism-related, and TF expression levels were 382 positively correlated with the expression of the analyzed cracking-related genes. Moreo-383 ver, the expression levels of the water transport-related and cuticle structure-related genes 384 were highly positively correlated with the expression of the starch/sucrose metabolism-385 related and Ca-related genes, respectively. 386

387
To validate the transcriptome data, we selected 11 cracking-related genes and evalu-388 ated their expression profiles using qRT-PCR. The gene-specific primers used in the qRT-389 PCR analysis are listed in Table S1. Figure 8 illustrates the expression levels of the DEGs 390 Figure 7. Co-expression of cracking-related genes. Heatmaps of the expression of candidate genes related to fruit cracking in two cracking-susceptible cultivars (A) and in cracking-susceptible and -resistant cultivars (B). (C) Network representation of cracking-related genes (orange circles: water transport, green circles: cell wall, blue circles: starch and sucrose metabolism, pink circles: hormone metabolism, green triangles: cuticle structure, blue triangles: calcium, pink triangles: transcription factor). Nodes represent genes, and edges represent relationships between any two genes. Edges represent positive correlations as determined by a Spearman correlation coefficient > 0.7. as determined by qRT-PCR and RNA-seq. The qRT-PCR results indicated that the expres-391 sion profiles of the 11 DEGs were consistent with the corresponding FPKM values derived 392 from the RNA-seq analysis. These genes exhibited a similar expression pattern when both 393 methods were used. The similar results that were obtained from the qRT-PCR and RNA-394 seq analyses indicated that the RNA-seq data were reproducible and reliable.

399
Fruit cracking is a serious disorder that negatively affects fruit marketability. Our 400 limited knowledge of the molecular characteristics of jujube has made it difficult to rec-401 ommend appropriate methods to prevent fruit cracking. In the present study, we hypoth-402 esized that jujube fruit cracking is due to an imbalance in the turgor pressure and exocarp 403 mechanical strength, which is regulated by a set of genes ( Figure S4). 404 Fruit cracking is the result of excessive water uptake, which causes the exocarp to 405 burst [4]. In the current study, one DEG (gene6278) involved in water transport was more 406 highly expressed in the cracking fruits than it was in the non-cracking fruits, and its ex-407 pression was gradually downregulated in the CJ, NJ, and NM samples ( Figure 7A,B). This 408 DEG encodes the aquaporin TIP1-3, which facilitates the transport of water and small 409 neutral solutes across cell membranes [32]. In lychee, an aquaporin gene is commonly 410 more highly expressed in cracking fruits than it is in non-cracking fruits [8]. We speculated 411 that TIP1-3 may be the key gene that controls jujube fruit cracking by regulating water 412 transport. 413 The osmotic potential of a fruit is another factor that is correlated with cracking. 414 Sugar is an important component of fruits because of its role in fruit growth and 415

Discussion
Fruit cracking is a serious disorder that negatively affects fruit marketability. Our limited knowledge of the molecular characteristics of jujube has made it difficult to recommend appropriate methods to prevent fruit cracking. In the present study, we hypothesized that jujube fruit cracking is due to an imbalance in the turgor pressure and exocarp mechanical strength, which is regulated by a set of genes ( Figure S4).
Fruit cracking is the result of excessive water uptake, which causes the exocarp to burst [4]. In the current study, one DEG (gene6278) involved in water transport was more highly expressed in the cracking fruits than it was in the non-cracking fruits, and its expression was gradually downregulated in the CJ, NJ, and NM samples ( Figure 7A,B). This DEG encodes the aquaporin TIP1-3, which facilitates the transport of water and small neutral solutes across cell membranes [32]. In lychee, an aquaporin gene is commonly more highly expressed in cracking fruits than it is in non-cracking fruits [8]. We speculated that TIP1-3 may be the key gene that controls jujube fruit cracking by regulating water transport.
The osmotic potential of a fruit is another factor that is correlated with cracking. Sugar is an important component of fruits because of its role in fruit growth and development. An increase in the total soluble sugar content in maturing fruits leads to a decrease in the tissue osmotic potential, which results in increased water absorption and growth as well as a greater likelihood of cracking [33]. Li et al. [34] reported that the total soluble sugar content is higher in cracking-susceptible cultivars than it is in cracking-resistant cultivars. Regarding the comparisons of the DEGs between cracking and non-cracking fruits, 'starch and sucrose metabolism' was revealed as an enriched KEGG pathway among the common DEGs. Additionally, the expression of starch/sucrose metabolism-related genes (gene21709, gene21813, and gene32506) was upregulated in cracking fruits ( Figure 7A). This finding is consistent with the results of an earlier study by Wang et al. [35], in which starch/sucrose metabolism-related genes were highly expressed in the cracking fruits of cracking-susceptible cultivars. In the present study, the expression of starch/sucrose metabolism-related genes was positively associated with the expression of water transportrelated genes. The transcripts for the genes involved in starch/sucrose metabolism were highly abundant in the cracking fruits of the cracking-susceptible cultivar, which may result in a decrease in the tissue osmotic potential. Notably, the TIP transcript levels were also high in the cracking fruits, which may lead to an increase in the amount of water entering fruits during heavy rainfall, thereby accelerating fruit cracking.
The mineral content is correlated with fruit cracking, and a preharvest spray application of mineral salts can restrict the uptake of water by fruits [33,36]. Calcium is critical for maintaining the structural integrity and firmness of fruit cell walls because it decreases cell wall permeability, leading to restricted water uptake [37]. The GO term 'Ca ion transport' was significantly enriched among the DEGs that were detected in the CC vs. NC comparison ( Figure 4). Two Ca transport-related genes (gene16635 and gene17677) were more highly expressed in the cracking fruits than they were in the non-cracking fruits. Moreover, four Ca transport-related genes (gene15072, gene28419, gene28420, and gene8932) had gradually downregulated expression levels in the CJ, NJ, and NM samples ( Figure 7B). These DEGs encode the Ca-binding protein CML (CML) or Ca-transporting ATPase (ACA), which are crucial for plant growth and development as well as for responses to various external stimuli [38,39]. The CML and ACA expression levels are reportedly downregulated in the cracked lychee pericarp [8]. According to Yang et al. [40], in cracking 'Lingwu Long-jujube' fruits, Ca 2+ precipitates are distributed at low levels in the exocarp cells, which is in contrast to the relatively large amounts of Ca 2+ precipitates in the mesocarp cell walls, cell wall edges, and intercellular space. The discrepancy between our results and those of an earlier investigation on lychee may be associated with the differences in the examined plant materials. In our study, the exocarp and mesocarp were sampled together for sequencing analyses. Accordingly, future studies will need to analyze the exocarp and mesocarp separately.
The exocarp developmental stages and mechanical properties, which are closely related to cell wall metabolism, affect the occurrence of fruit cracking [41]. In our study, several cell wall metabolism-related GO terms, including 'cell wall' and 'plant-type cell wall', were assigned to the DEGs that were detected in three comparisons. The Venn diagram and trend analyses confirmed that the metabolic activities of the cell wall were enriched during fruit cracking. According to the DEG annotations, 12 cell wall metabolismrelated genes might be involved in jujube fruit cracking, of which two XTH genes (gene14330 and gene19193), two expansin genes (EXP, gene21836 and gene8126), and one PME gene (gene28523) were more highly expressed in the cracking fruits than in they were in the non-cracking fruits ( Figure 7A,B). Similar results were reported for atemoya, in which the expression levels of most PME, EXP, and XTH genes are higher in cracking fruits than they are in non-cracking fruits [10]. The suppressed expression of a polygalacturonase (SlPG) and SlEXP1 in ripening tomato fruits decreases cell wall disassembly, which leads to increased cracking resistance [42]. Thus, these cell wall-related genes contribute to cell wall disassembly and decrease the elasticity of the exocarp.
The structure and function of the cuticle, which is an important part of the fruit exocarp, have important effects on the mechanisms underlying fruit cracking [43]. The plant cuticle comprises a complex mixture of aliphatic compounds, which are composed of lignin and highly unsaturated FAs [44]. Lignin is an abundant phenylpropanoid polymer [45], and α-linolenic acid is a major unsaturated FA [46]. In the present study, 'phenylpropanoid biosynthesis' was one of the top 20 KEGG pathways among the DEGs detected in the CC vs. NC, CJ vs. NJ, and NJ vs. NM comparisons (Figure 3). Regarding the common DEGs in two comparisons of cracking and non-cracking fruits and the gradually downregulated DEGs in jujube cultivars that are susceptible and resistant to fruit cracking, the functional annotations revealed that these DEGs are related to unsaturated FA oxidation and many processes affecting the cuticle structure ( Figure 4). The high expression levels of a lipoxygenase (LOX) in the lychee cultivar 'Baitangying' suggest that the cuticle of this cracking-susceptible cultivar might have relatively few unsaturated FAs [9]. In our study, two LOX genes (gene4863 and gene5490) were highly expressed in cracking jujube fruits ( Figure 7A), implying that they may positively affect fruit cracking-related networks. Hence, we speculate that the susceptibility of jujube cultivars to fruit cracking is related to the changes in the oxidation of unsaturated FAs and the cuticle structure.
Hormonal homeostasis is an important factor that influences fruit cracking [32]. The accumulation of IAA and ABA in the aril is necessary for the growth and development of the aril, but high IAA or ABA contents can induce fruit cracking in lychee [47]. Auxin regulates cell size, stimulates carbohydrate translocation to the fruit, and increases cell wall elasticity [48]. Applications of ABA increase cracking because this treatment increases water movement into and promotes the enlargement of the fruit [42,49]. In this study, the KEGG functional analysis of Profile 7 (gradually downregulated profile) revealed a significant enrichment in the ABA-related 'carotenoid biosynthesis' pathways ( Figure 5B). An analysis of the data generated in the current study identified one IAA metabolism-related gene (gene15774) and one ABA metabolism-related gene (gene30271) with important roles during the fruit cracking process and in cultivars that differ in terms of fruit cracking resistance ( Figure 7A,B). The gene15774 encodes an auxin-binding protein that probably mediates the rapid effects of IAA on cell elongation [50]. The protein encoded by NCED1 (gene30271) catalyzes key ABA biosynthesis steps, resulting in an increase in fruit firmness and an extended shelf-life [51]. Abscisic acid affects cell wall catabolism in ripening tomato fruits by downregulating the expression of major catabolic genes [51]. Our co-expression analysis detected a strong correlation between hormone metabolism and other cracking-related metabolic activities. Consistent with the findings of previous studies on atemoya [10] and lychee [12,47], we observed that the upregulated expression of ABA-related and IAArelated DEGs in cracking jujube fruits can lead to increased ABA and IAA levels in fruits, which might promote fruit cracking directly or indirectly by upregulating the expression of specific TFs, such as MYB, WRKY, and NAC, to increase growth-related gene expression.
Transcription factors play critical roles in responding to various types of abiotic stress, including water stress [52]. In the current study, we determined that the genes encoding TFs from several families, including the WRKY (gene30640, gene1197, gene12181, gene17919, gene21416, and gene22590), MYB (gene12280, gene24181, and gene7685), and ERF (gene6845, gene21778, gene32981, gene21158, and gene9811) families were more highly expressed in cracking fruits than they were in non-cracking fruits. Moreover, the expression levels of gene6845 (ERF), gene12280 (MYB), and gene30640 (WRKY) were gradually downregulated in the CJ, NJ, and NM samples ( Figure 7A,B). Liao et al. [53] identified ClERF4 as a candidate gene that was responsible for rind hardness, which is associated with watermelon fruit cracking resistance. The overexpression of the grape VvWRKY11 gene in Arabidopsis thaliana can adjust the osmotic potential to increase the resistance to cold conditions [54]. In another study involving A. thaliana, the expression of the sheepgrass LcMYB4 gene positively regulated the accumulation of soluble sugars in transgenic plants [55]. Transcription factors regulate plant growth and development by modulating the transcription of specific target genes [12]. In our study, we revealed a positive correlation between TFs and fruit crackingrelated processes. However, precisely how these TFs are involved in jujube fruit cracking remains to be determined.

Conclusions
We conducted comprehensive RNA-seq analyses of five jujube fruit types, namely CC, NC, CJ, NJ, and NM. A total of 131 DEGs were common to the CC vs. NC and CJ vs. NJ comparisons, and 225 DEGs were gradually downregulated in the CJ, NJ, and NM samples. Functional enrichment analyses suggested that DEGs were related primarily to cell wall metabolism, cuticle structure, and hormone metabolism, which highlight the importance of the metabolic pathway that is involved in fruit cracking. Among the two categories (131 and 225) of DEGs, 33 and 34 genes were closely associated with fruit cracking, respectively. Six genes (TIP1-3, NCED1, ABP19A, ERF2, MYB, and WRKY40) were common to the Venn diagram and trend analyses. The upregulation of the genes associated with starch/sucrose metabolism and water transport might increase the expanding pressure, whereas the genes related to cell wall metabolism and cuticle structure might decrease the exocarp mechanical strength, thereby making jujube fruit susceptible to cracking. The TFs and genes associated with hormones, such as auxin and ABA, play critical roles in the regulation of fruit cracking via the above-named pathways. Our findings may help to elucidate the molecular mechanisms underlying fruit cracking.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/genes13010105/s1, Figure S1: The saturation curve of different ripening stages. Figure S2: The Coverage of sequencing. Figure S3: Distribution of differentially expressed transcription factors in gene families. Figure S4: Hypothesized mechanisms of fruit cracking in jujube. Table S1: Primer used for qRT-PCR. Table S2: Summary of sequence assembly after Illumina sequencing. Table S3: Clean reads mapped to the reference genome. Table S4: Complete lists of differentially expressed genes in three groups based on pairwise comparisons. Table S5: Annotation summaries of DEGs.