RNA-Seq Provides New Insights into the Molecular Events Involved in “Ball-Skin versus Bladder Effect” on Fruit Cracking in Litchi

Fruit cracking is a disorder of fruit development in response to internal or external cues, which causes a loss in the economic value of fruit. Therefore, exploring the mechanism underlying fruit cracking is of great significance to increase the economic yield of fruit trees. However, the molecular mechanism underlying fruit cracking is still poorly understood. Litchi, as an important tropical and subtropical fruit crop, contributes significantly to the gross agricultural product in Southeast Asia. One important agricultural concern in the litchi industry is that some famous varieties with high economic value such as ‘Nuomici’ are susceptible to fruit cracking. Here, the cracking-susceptible cultivar ‘Nuomici’ and cracking-resistant cultivar ‘Huaizhi’ were selected, and the samples including pericarp and aril during fruit development and cracking were collected for RNA-Seq analysis. Based on weighted gene co-expression network analysis (WGCNA) and the “ball-skin versus bladder effect” theory (fruit cracking occurs upon the aril expanding pressure exceeds the pericarp strength), it was found that seven co-expression modules genes (1733 candidate genes) were closely associated with fruit cracking in ‘Nuomici’. Importantly, we propose that the low expression level of genes related to plant hormones (Auxin, Gibberellins, Ethylene), transcription factors, calcium transport and signaling, and lipid synthesis might decrease the mechanical strength of pericarp in ‘Nuomici’, while high expression level of genes associated with plant hormones (Auxin and abscisic acid), transcription factors, starch/sucrose metabolism, and sugar/water transport might increase the aril expanding pressure, thereby resulting in fruit cracking in ‘Nuomici’. In conclusion, our results provide comprehensive molecular events involved in the “ball-skin versus bladder effect” on fruit cracking in litchi.


Introduction
Fruit cracking is a physiological phenomenon where pericarp cracks due to the pressure from fruit internal growth exceeds the strength given by pericarp growth [1]. Fruit cracking is common in fruit crops such as apple [2], citrus [3], peach [4], and grape [5]. The economic value will be lost when fruit undergo cracking. Therefore, exploring the mechanism underlying fruit cracking has always been an important agricultural concern.
The occurrence of fruit cracking is closely related to internal or external environmental factors, of which calcium nutrition, temperature, and humidity are considered as the main physiological factors [6]. In general, calcium deficiency, low temperature, high air humidity, and soil drought will cause fruit cracking. Consistently, calcium spraying can enhance the pericarp strength via increasing the content of calcium in pericarp, thereby reducing the fruit cracking rate in pomegranate, kiwifruit, and litchi [7][8][9]. In pomegranate, if the fruit are exposed to drought stress during the early stage of development, the pulp will accumulate more dry matter and have lower osmotic potential, as a result, the pulp will expand significantly faster than pericarp upon fruit encountering heavy rain during fruit ripening, and fruit cracking occurs [10].
Fruit cracking is also closely associated with genetic factors as early as 1981, when Cuartero et al. pointed out that fruit cracking was regulated by a set of genes [11]. To date, it has been accepted that cell wall remodeling genes including expansin (EXP), endo-1,4-βglucanase (EG), polygalacturonase (PG), and xyloglucan endotrans-glycosylase (XET) play an important role in the control of fruit cracking since these genes can regulate the mechanical strength of the cell wall. In apple, the earlier expression of MdEXPA3 in pericarp by fruit bagging can induce the expansion of pericarp cells earlier, thereby decreasing the susceptibility of fruit cracking [12]. In sweet cherry, the expression level of genes related to cell wall remodeling and cuticular wax biosynthesis is higher in the pericarp of 'Kordia' than that in 'Bing', the difference between the expression of these genes may partially explain why 'Kordia' is more resistant to fruit cracking than 'Bing' [13]. Recently, RNA-Seq on both cracking-resistant and cracking-susceptible cultivars identified genes encoding XET1, XET2, glycerolipid triacylglycerol (TAG), and MADS box transcription factors (MDP) involved in fruit cracking in watermelon [14]. In addition, the decrease in the expression of PG and EXP via the RNA interference (RNAi) approach in tomato can promote generation of cuticular wax, protochitin, and cellulose in pericarp to increase its extensibility, therefore reducing the fruit cracking rate [15]. Taken together, though these findings reveal that genes related to cell wall remodeling, wax biosynthesis, and calcium signaling play an important role in fruit cracking, the molecular mechanism underlying this process is yet largely unknown.
Litchi (Litchi chinensis Sonn.), which originated in southern China, is an important tropical and subtropical fruit crop. Now, it is widely cultivated in over 20 countries such as India, Vietnam, Thailand, Madagascar, Australia, and South Africa. The area of industrial litchi cultivation across the world is about 81.9 million hectares, with a total of output of 4.3 million tons. However, fruit cracking has become a limiting factor in the sustainable development of some famous varieties such as 'Nuomici', which experiences a fruit cracking rate between 10-30% in general years, but in some special years such as when fruit encounter heavy rain during ripening, the fruit cracking rate is almost 90%. Litchi fruit growth pattern can be clearly divided into two stages. Stage I (last for about 50 days post anthesis), which accounts for two thirds of the whole fruit development period, is characterized by a slow growth of pericarp and seed coat, while stage II is characterized by a fast growth of aril and embryo [16]. It has been pointed out that the growth of aril begins when the pericarp grows almost to full size as the potential of aril growth is dependent on the space provided by the pericarp. This relationship between the growth patterns of pericarp and aril is called the "ball-skin versus bladder effect" theory [17]. According to this theory, a decrease in pericarp strength will provide a relatively smaller space for aril expansion, fruit cracking will occur when the aril expands, particularly under high temperature or high air humidity, and or a sudden heavy rain since these environmental conditions can facilitate the aril expanding [18]. It has been found that the expression level of two cell wall remodeling genes (LcXET1 and LcEG1) is higher in the pericarp of cracking-susceptible cultivar 'Nuomici' than that in the cracking-resistant cultivar 'Huaizhi', suggesting that these two genes might be involved in fruit cracking in litchi via regulating the pericarp extensibility [19,20]. Recently, two studies using transcriptome analysis in combination with metabolomic profiling found that the genes related to cell wall, water transport, calcium transport, and plant hormones were differentially expressed during the process of fruit cracking in the cracking-susceptible cultivar 'Baitangying' compared with the cracking-resistant cultivar 'Feizixiao', and suggested that the changes in the expression of these genes might lead to the decrease in the cell wall stability of pericarp, faster water transport, and the imbalance of endogenous hormones, thereby inducing fruit cracking in 'Baitangying' [21,22]. Though these findings provide new clues into the molecular mecha-nism underlying fruit cracking in litchi, two aspects should be given more attention: (i) the "ball-skin versus bladder effect" exists in litchi fruit cracking, which occurs as a result of interaction between pericarp strength and aril expanding pressure, thus it is necessary to explore the molecular events in both pericarp and aril; and (ii) fruit cracking in litchi, which is not a sudden phenomenon, is closely associated with fruit development at an early stage. Therefore, the differentially expressed genes (DEGs) during early fruit development (prior to fruit cracking) should be coordinated.
Weighted gene co-expression network analysis (WGCNA), which can connect gene expression patterns with phenotypes, has become a powerful strategy to identify key candidate genes that are correlated to specific phenotypic traits [23]. Thus, in this study, we selected the cracking-susceptible litchi cultivar 'Nuomici' and cracking-resistant cultivar 'Huaizhi' as materials, and samples including both pericarp and aril were collected from developing fruit and fruit that are under cracking for RNA sequencing. We then connected the differentially expressed genes with the daily cracking rate (DCR) by WGCNA to obtain the key modules and genes involved in fruit cracking. Finally, a preliminary model of molecular events in control of fruit cracking in 'Nuomici' litchi was discussed.

Dynamic of Fruit Cracking Rate between 'Nuomici' and 'Huaizhi'
As shown in Figure 1, fruit cracking occurred beginning at 55 days after anthesis in 'Nuomici (NMC)' and the daily cracking rate (DCR) was increased thereafter, while the 'Huaizhi (HZ)' fruit did not crack throughout the whole fruit growth and development period. The DCR of 'Nuomici' was between 55 and 62 days, and 62 and 69 days after anthesis was 5.01% and 5.77%, respectively ( Figure 1B). As a result, the cumulative fruit cracking rate of 'Nuomici' was as high as 43.12% ( Figure 1C). 'Baitangying' [21,22]. Though these findings provide new clues into the m anism underlying fruit cracking in litchi, two aspects should be given mo the "ball-skin versus bladder effect" exists in litchi fruit cracking, which oc of interaction between pericarp strength and aril expanding pressure, thus to explore the molecular events in both pericarp and aril; and (ii) fruit cra which is not a sudden phenomenon, is closely associated with fruit deve early stage. Therefore, the differentially expressed genes (DEGs) during ea opment (prior to fruit cracking) should be coordinated.
Weighted gene co-expression network analysis (WGCNA), which can expression patterns with phenotypes, has become a powerful strategy to id didate genes that are correlated to specific phenotypic traits [23]. Thus, in selected the cracking-susceptible litchi cultivar 'Nuomici' and cracking-re 'Huaizhi' as materials, and samples including both pericarp and aril were developing fruit and fruit that are under cracking for RNA sequencing. nected the differentially expressed genes with the daily cracking rate (DCR to obtain the key modules and genes involved in fruit cracking. Finally, model of molecular events in control of fruit cracking in 'Nuomici' litchi w

Dynamic of Fruit Cracking Rate between 'Nuomici' and 'Huaizhi'
As shown in Figure 1, fruit cracking occurred beginning at 55 days a 'Nuomici (NMC)' and the daily cracking rate (DCR) was increased therea 'Huaizhi (HZ)' fruit did not crack throughout the whole fruit growth and period. The DCR of 'Nuomici' was between 55 and 62 days, and 62 and anthesis was 5.01% and 5.77%, respectively ( Figure 1B). As a result, the cu cracking rate of 'Nuomici' was as high as 43.12% ( Figure 1C).

Transcriptome Sequencing, Mapping, and Annotation
RNA was extracted from pericarp and aril, respectively, and quality was analyzed using an Agilent 2100 bioanalyzer. Subsequently, sixty RNA samples were sequenced with an Illumina HiSeq 2500. After stringent quality assessment and data filtering, 646 million reads were produced by 150 bp paired end sequencing ( Table 1). The reads with base quality greater than Q30 (Q ≥ 30) and no ambiguous "N" were defined as high-quality. Approximately 78% of the reads (502 million mapped reads) were successfully aligned to the litchi genome using hisat software (Table 1). For annotation, the gene expressed in pericarp and aril were compared against four databases (NR: RefSeq non-redundant proteins, SWiss-Prot: Swiss-Prot Protein Sequence Database, COG: cluster of orthologous group, and KEGG: Kyoto Encyclopedia of Genes and Genomes) using BLASTX. A large number of litchi genes including 27,621 genes in Nr, 20,815 genes in Swiss-Prot, 10,507 genes in KEGG, and 9912 in COG were annotated ( Figure 2), of which 5204 genes could be shared by these four databases, and most genes could be annotated at whatever two databases. RNA was extracted from pericarp and aril, respectively, and quality was analyzed using an Agilent 2100 bioanalyzer. Subsequently, sixty RNA samples were sequenced with an Illumina HiSeq 2500. After stringent quality assessment and data filtering, 646 million reads were produced by 150 bp paired end sequencing ( Table 1). The reads with base quality greater than Q30 (Q ≥ 30) and no ambiguous "N" were defined as high-quality. Approximately 78% of the reads (502 million mapped reads) were successfully aligned to the litchi genome using hisat software (Table 1). Q30 percent: The percentages of the bases whose Phred values were more than 30; GC content: The content of guanine and cytosine for the library in one tissue.
For annotation, the gene expressed in pericarp and aril were compared against four databases (NR: RefSeq non-redundant proteins, SWiss-Prot :Swiss-Prot Protein Sequence Database, COG: cluster of orthologous group, and KEGG: Kyoto Encyclopedia of Genes and Genomes) using BLASTX. A large number of litchi genes including 27,621 genes in Nr, 20,815 genes in Swiss-Prot, 10,507 genes in KEGG, and 9912 in COG were annotated ( Figure 2), of which 5204 genes could be shared by these four databases, and most genes could be annotated at whatever two databases.   First, differentially expressed genes (DEGs) in pericarp or aril were screened at the same time point using FPKM (Fragments Per Kilobase Million) values according to differential expression greater than two-fold (False discovery rate, FDR < 0.01, log 2 > 1 or < −1). Then, the DEGs in pericarp and aril at all time points were integrated and defined as P1 and A1, respectively. Results showed a total of 5624 genes in P1 and 4399 genes in A1. Then, DEGs in P1 and A1 were divided into 12 and six independent modules according to WGCNA, respectively ( Figure 3A,B). In pericarp, two module genes (turquoise and blue modules) were found to be significantly negative correlated to DCR, while two module genes (brown and green modules) were significantly positively correlated to DCR. In aril, two module genes (green and turquoise modules) were found to be significantly positively correlated to DCR, and one module (blue module) was significantly negatively correlated to DCR ( Figure 3C,D). Here, the DEGs that were significantly correlated to DCR in pericarp (3080 genes) and aril (2402 genes) were defined as PM and AM, respectively.

Screening of Genes Related to Fruit Cracking by Weighted Gene Co-Expression Network Analysis (WGCNA)
First, differentially expressed genes (DEGs) in pericarp or aril were screened at the same time point using FPKM (Fragments Per Kilobase Million) values according to differential expression greater than two-fold (False discovery rate, FDR< 0.01, log2 >1 or < −1). Then, the DEGs in pericarp and aril at all time points were integrated and defined as P1 and A1, respectively. Results showed a total of 5624 genes in P1 and 4399 genes in A1. Then, DEGs in P1 and A1 were divided into 12 and six independent modules according to WGCNA, respectively ( Figure 3A,B). In pericarp, two module genes (turquoise and blue modules) were found to be significantly negative correlated to DCR, while two module genes (brown and green modules) were significantly positively correlated to DCR. In aril, two module genes (green and turquoise modules) were found to be significantly positively correlated to DCR, and one module (blue module) was significantly negatively correlated to DCR ( Figure 3C,D). Here, the DEGs that were significantly correlated to DCR in pericarp (3080 genes) and aril (2402 genes) were defined as PM and AM, respectively. To further screen DEGs that are more closely related to fruit cracking, we first compared the FPKM values between 'Nuomici' and 'Huaizhi' in pericarp and aril at 62 d and 69 d after anthesis since the fruit of 'Nuomici' were undergoing cracking during this period. We named the DEGs in pericarp (2480 genes) and aril (3158 genes) at these two time points as P2 and A2, respectively. We then compared the DEGs between P2 and PM, A2 and AM. We found that 639 DEGs were shared between P2 and PM and 1097 DEGs were common between A2 and AM ( Figure 4). Thus, these DEGs were used for further analysis. Genes with similar expression patterns were clustered. Each branch in the (A) and (B) represents one gene, and the color below each branch represents the co-expression module. Twelve and six modules in (C) and (D) are grouped by genes with the same color. The asterisk (*) indicates that the module has a significant correlation coefficient (p < 0.05) with DCR. The number above (or below) the asterisk represents the gene numbers in this module.
To further screen DEGs that are more closely related to fruit cracking, we first compared the FPKM values between 'Nuomici' and 'Huaizhi' in pericarp and aril at 62 d and 69 d after anthesis since the fruit of 'Nuomici' were undergoing cracking during this period. We named the DEGs in pericarp (2480 genes) and aril (3158 genes) at these two time points as P2 and A2, respectively. We then compared the DEGs between P2 and PM, A2 and AM. We found that 639 DEGs were shared between P2 and PM and 1097 DEGs were common between A2 and AM ( Figure 4). Thus, these DEGs were used for further analysis.

Clustering of Differentially Expressed Genes
According to the "ball-skin versus bladder effect" theory, we hypothesized that the

Clustering of Differentially Expressed Genes
According to the "ball-skin versus bladder effect" theory, we hypothesized that the genes involved in fruit cracking in litchi might be closely related to pericarp extensibility and rapid aril expanding, so we performed cluster analysis of the DEGs in pericarp and aril via SOTA (Self-Organizing Tree Algorithm) clustering. As shown in Figure 5, the DEGs (639 genes) in pericarp were divided into seven clusters (A1-A7), of which 241 DEGs were grouped into A1 cluster, which accounts for the largest proportion, and their expression level was higher in 'Huaizhi' than that in 'Nuomici' during the whole fruit development period. A total of 180 DEGs were grouped into cluster A5 and their expression level was higher in 'Huaizhi' than that in 'Nuomici' beginning at 40 d post anthesis. In contrast, the expression level of genes in the A7 cluster was higher in 'Nuomici' than that in 'Huaizhi' beginning at 40 d post anthesis. In addition, the DEGs in cluster A4 and A6 might be less closely related to the pericarp extensibility when compared with that in other clusters since the genes in cluster A4 and A6 were expressed differentially at only one time point. Thus, DEGs (583 genes) in clusters except A4 and A6 were selected for further analysis. The DEGs (1097 genes) in aril were divided into four clusters (B1-B4). The expression of genes in cluster B1 did not change in 'Nuomici', but showed a higher expression level in 'Huaizhi'. The DEGs in cluster B2 displayed a decreased expression trend in both 'Nuomici' and 'Huaizhi', but their transcripts were higher in 'Huaizhi' than 'Nuomici'. The DEGs in both cluster B3 and B4 did not alter in 'Huaizhi', but showed a higher expression level in 'Nuomici'. The difference is that the DEGs of cluster B3 in 'Nuomici' displayed a constant decline expression trend, while the expression pattern of DEGs in cluster B4 in 'Nuomici' was first increased and then remained unchanged.

Functional Annotation of Differentially Expressed Genes
To explore the functions of DEGs in the regulation of fruit cracking in litchi, we performed functional annotation of candidate DEGs in pericarp (583 genes) and aril (1097 genes). As shown in Figure 6, both DEGs in pericarp and aril were divided into  16 groups including calcium, cell wall, development, DNA/RNA, energy/chloroplast/ mitochondria, hormone metabolism, lipid metabolism, starch/sucrose metabolism, oxidation/reduction, protein/ribosome, secondary metabolism, signaling, stress, transcription factor, transporter, and others/unknown. Next, to better explore the function of these DEGs in the regulation of fruit cracking, the patterns of their expression level in 'Nuomici' compared with that in 'Huaizhi' was displayed by heatmaps. As shown in Figure 7, in pericarp, there were 14 DEGs related to calcium, all of them were significantly downregulated in 'Nuomici' involving calcium transport, binding, and signaling. Eight DEGs were found to be related to cell wall metabolism, of which four genes were significantly upregulated in 'Nuomici', involving cell wall degradation and modification. Interestingly, a total of 30 DEGs were related to hormones. Among them, 26 were significantly downregulated in 'Nuomici', of which eight were associated with auxin synthesis, transportation, and signal transduction. Fourteen DEGs were related to lipid metabolism. Among them, 12 DEGs were significantly downregulated in 'Nuomici' involving wax synthesis. Thirty-three DEGs were related to transcription factors, of which 29 DEGs were significantly downregulated in 'Nuomici' including the WRKY, bZIP, bHLH, and MYB transcription factors. In aril (Figure 8), there were 35 DEGs related to hormones, of which 24 were significantly upregulated in 'Nuomici', 10 were involved in auxin synthesis, transport, and signal transduction. Ten DEGs were related to starch/sucrose metabolism, of which nine were significantly upregulated in 'Nuomici' involving sucrose, β-glucoside, and trehalose transport. There were 65 DEGs that were transcription factors including WRKY, bHLH, DOF, NAC, and MYB. Forty-six DEGs were related to transport including 35 upregulated genes and nine downregulated genes in 'Nuomici' involving water, carbohydrate, and mineral elements and the transport of various small molecule compounds. Taken together, these results indicate that most of the DEGs involved in calcium, hormone metabolism, and wax metabolism were significantly downregulated in 'Nuomici' pericarp, while genes related to cell wall remodeling were significantly upregulated in 'Nuomici' pericarp. Additionally, most of the DEGs involved in hormone metabolism, starch/sucrose metabolism, and transport pathway were significantly upregulated in 'Nuomici' aril.
including calcium, cell wall, development, DNA/RNA, energy/chloroplast/mitocho hormone metabolism, lipid metabolism, starch/sucrose metabolism, oxidation/redu protein/ribosome, secondary metabolism, signaling, stress, transcription factor, porter, and others/unknown. Next, to better explore the function of these DEGs regulation of fruit cracking, the patterns of their expression level in 'Nuomici' com with that in 'Huaizhi' was displayed by heatmaps. As shown in Figure 7, in pericarp were 14 DEGs related to calcium, all of them were significantly downregula 'Nuomici' involving calcium transport, binding, and signaling. Eight DEGs were fo be related to cell wall metabolism, of which four genes were significantly upregula 'Nuomici', involving cell wall degradation and modification. Interestingly, a tota DEGs were related to hormones. Among them, 26 were significantly downregula 'Nuomici', of which eight were associated with auxin synthesis, transportation, and transduction. Fourteen DEGs were related to lipid metabolism. Among them, 12 were significantly downregulated in 'Nuomici' involving wax synthesis. Thirty DEGs were related to transcription factors, of which 29 DEGs were significantly regulated in 'Nuomici' including the WRKY, bZIP, bHLH, and MYB transcription f In aril (Figure 8), there were 35 DEGs related to hormones, of which 24 were signifi upregulated in 'Nuomici', 10 were involved in auxin synthesis, transport, and transduction. Ten DEGs were related to starch/sucrose metabolism, of which nine significantly upregulated in 'Nuomici' involving sucrose, β-glucoside, and tre transport. There were 65 DEGs that were transcription factors including WRKY, DOF, NAC, and MYB. Forty-six DEGs were related to transport including 35 upreg genes and nine downregulated genes in 'Nuomici' involving water, carbohydrat mineral elements and the transport of various small molecule compounds. Tak gether, these results indicate that most of the DEGs involved in calcium, hormone m olism, and wax metabolism were significantly downregulated in 'Nuomici' pe while genes related to cell wall remodeling were significantly upregulated in 'Nu pericarp. Additionally, most of the DEGs involved in hormone metabolism, star crose metabolism, and transport pathway were significantly upregulated in 'Nu aril.

Validation of RNA-Seq Data
In order to verify the RNA-Seq data, 18 DEGs in pericarp and 22 DEGs in aril w randomly selected for qRT-PCR (Quantitative Real-time Polymerase Chain Reacti analysis. As shown in Figure 9, the qRT-PCR results showed that the expression patte

Validation of RNA-Seq Data
In order to verify the RNA-Seq data, 18 DEGs in pericarp and 22 DEGs in aril were randomly selected for qRT-PCR (Quantitative Real-time Polymerase Chain Reaction) analysis.
As shown in Figure 9, the qRT-PCR results showed that the expression patterns of 12 genes in pericarp were consistent with RNA-Seq results including PAE1(Pectinesterase inhibitors, Lc. 10 Next, the correlation coefficients between RNA-Seq and qRT-PCR were calculated ( Figure 10). We found that the pairwise correlation coefficient reached a significant level higher than 0.6, indicating that expression data obtained from RNA-Seq could reflect the real transcript abundance.

Discussion
Litchi fruit cracking is a physical process of mechanical breaking. An antagonism exists between the internal aril and external pericarp, called the "ball-skin versus bladder effect". Fruit cracking will occur upon the pressure from the internal aril expanding, which exceeds the strength from the pericarp growth.
The pericarp strength is dependent on the pericarp development status and structure. Endogenous hormones play an important role in pericarp development. Previous studies have shown that the application of NAA (Naphthylacetic acid) can reduce the fruit cracking rate in several fruit trees including litchi, citrus, sweet cherry, and apple, which may be due to the promotion of pericarp growth by NAA [24][25][26][27]. Recently, based on the analysis of metabolism in pericarp, Wang et al. pointed out that IAA deficiency in pericarp might be one main factor inducing fruit cracking in 'Baitangying' litchi [22] were downregulated in the pericarp of cracking-susceptible cultivar 'Nuomici' when compared with that in cracking-resistant cultivar 'Huaizhi', which might lead to the retardation of pericarp growth and development in 'Nuomici'. In addition, four genes that encode the rate-limiting enzymes for ethylene biosynthesis (Lc.4.493: 1-aminocyclopropane -1 -carboxlic acid synthas, Lc.9.493: l-aminocycloProPane-1-carboxylic acid oxidas, Lc.9.503: l-aminocycloPro-Pane-1-carboxylic acid oxidas, Lc.9.506: l-aminocycloProPane-1-carboxylic acid oxidas) and two gibberellin-regulated genes (Lc.1.532 and Lc.1.534) were found to be downregulated in 'Nuomici' pericarp. It has been reported that ethephon and gibberellin can reduce the cracking rate of litchi [28,29]. We thus speculated that both ethylene and gibberellin might coordinate with auxin to positively regulate the development of litchi pericarp. In one hand, calcium, as an important component of cell wall to improve its mechanical strength, plays an important role in preventing pericarp cracking [30]. On one hand, calcium, as an important development and stress response [31]. Calcium sensors can transfer intracellular Ca 2+ signals into cellular physiological responses by binding Ca 2+ to their own EF hand binding domain. Calmodulin like protein (CML) is one of the important calcium receptors Figure 10. Coefficient analysis between gene expression patterns between qRT-PCR and RNA-Seq data. Extremely significant difference levels (p < 0.01) are indicated with two asterisks (**).

Discussion
Litchi fruit cracking is a physical process of mechanical breaking. An antagonism exists between the internal aril and external pericarp, called the "ball-skin versus bladder effect". Fruit cracking will occur upon the pressure from the internal aril expanding, which exceeds the strength from the pericarp growth.
The pericarp strength is dependent on the pericarp development status and structure. Endogenous hormones play an important role in pericarp development. Previous studies have shown that the application of NAA (Naphthylacetic acid) can reduce the fruit cracking rate in several fruit trees including litchi, citrus, sweet cherry, and apple, which may be due to the promotion of pericarp growth by NAA [24][25][26][27]. Recently, based on the analysis of metabolism in pericarp, Wang et al. pointed out that IAA deficiency in pericarp might be one main factor inducing fruit cracking in 'Baitangying' litchi [22]. In this study, eight auxin-related DEGs (Lc.2.1082: AUX/IAA protein, Lc.3.1076: indole-3-acetate O-methyltransferase, Lc.3.1569: Auxin transporter-like protein, Lc.4.1204: AUX/IAA protein, Lc.7.1212: Gretchen Hagen3, Lc.8.1468: Auxin-induced protein, Lc.14.1070: ARF, Lc.14.1626: ARF) were downregulated in the pericarp of cracking-susceptible cultivar 'Nuomici' when compared with that in cracking-resistant cultivar 'Huaizhi', which might lead to the retardation of pericarp growth and development in 'Nuomici'. In addition, four genes that encode the rate-limiting enzymes for ethylene biosynthesis (Lc.4.493: 1-aminocyclopropane -1 -carboxlic acid synthas, Lc.9.493: l-aminocycloProPane-1-carboxylic acid oxidas, Lc.9.503: l-aminocycloProPane-1-carboxylic acid oxidas, Lc.9.506: l-aminocycloProPane-1-carboxylic acid oxidas) and two gibberellin-regulated genes (Lc.1.532 and Lc.1.534) were found to be downregulated in 'Nuomici' pericarp. It has been reported that ethephon and gibberellin can reduce the cracking rate of litchi [28,29]. We thus speculated that both ethylene and gibberellin might coordinate with auxin to positively regulate the development of litchi pericarp. In one hand, calcium, as an important component of cell wall to improve its mechanical strength, plays an important role in preventing pericarp cracking [30]. On one hand, calcium, as an important development and stress response [31]. Calcium sensors can transfer intracellular Ca 2+ signals into cellular physiological responses by binding Ca 2+ to their own EF hand binding domain. Calmodulin like protein (CML) is one of the important calcium receptors [32,33]. It has been shown that spraying calcium compounds can reduce the fruit cracking rate of sweet cherry and grape [34,35]. In litchi, one study also showed that the application of calcium chloride or calcium nitrate at the early stage of fruit development can increase the calcium concentration in litchi pericarp and reduce the fruit cracking rate [9]. In our study, we found that 14 DEGs in pericarp were related to calcium and all of them were downregulated in 'Nuomici' including six calmodulin like proteins (Lc.0.2152, Lc.0.3086, Lc.7.1332, Lc.13.587, Lc.13.888, Lc.14.824). We thus speculated that the lower expression level of calcium sensor genes might be due to lower calcium content in 'Nuomici' pericarp. The extensibility of pericarp is also closely related to wax deposition on the surface of pericarp [36]. C18 fatty acyl-ACPs, which can elongate fatty acid elongase (FAE) complexes, play important roles in the regulation of the biosynthesis of wax compounds [37]. Here, we found that the expression of two acyl carrier protein genes (Lc.8.658, Lc.8.2082) was significantly lower in 'Nuomici' pericarp than that in 'Huaizhi', implying that the wax component of the cell wall improved its mechanical strength and plays an important role in preventing pericarp cracking [30]. On the other hand, Ca 2+ , which can also act as a second messenger to regulate plant deposition in 'Nuomici' pericarp, was less than that in 'Huaizhi'. The increased expression of cell wall remodeling genes such as PME, PG, and EXP can disassemble the cell wall and reduce the mechanical strength of the pericarp [38,39]. In this study, it was found that the expression level of one expansin (Lc.1.2435) and one β-d-xylanase gene (Lc.0.4352) was higher in 'Nuomici' pericarp than that in 'Huaizhi', suggesting that the strength of 'Nuomici' pericarp was reduced, which is consistent with previous findings that the content of cellulose, hemicellulose, and water-insoluble pectin in the pericarp of 'Huaizhi' was higher than that in 'Nuomici' [40]. In addition, we also found that most of the DEGs related to development were downregulated in 'Nuomici' pericarp. These genes were reported to be involved in cell division and differentiation (Figure 7). Mutation of ACTIN-RELATED PROTEINS 2/3 in Arabidopsis results in misdirected expansion of various cell types [41]. Arabidopsis wat1 mutants have a defect in cell elongation and display a dwarf phenotype [42]. Interestingly, these development-related genes were downregulated in 'Nuomici' pericarp, indicating that the development of 'Nuomici' pericarp may be hindered. In support of this hypothesis, Wang et al. found that the pericarp cells of the cracking-susceptible cultivar 'Hongpinuo' litchi stopped cell division earlier than that in the cracking-resistant cultivar 'Baipinuo' litchi [43].
Litchi fruit cracking is not only caused by the gradual decrease in pericarp strength, but also by the sudden increase in the pressure resulting from the aril expanding [18]. The aril pressure comes from the aril expanding, which is closely related to its osmotic potential. It has been suggested that a decrease in osmotic potential in litchi aril is associated with the degradation of macromolecules and the accumulation of assimilates in aril, and the acceleration of solute transport from other tissues into the aril [44]. In this study, we found that genes encoding one glycosyltransferase family protein (Lc.12.1389), four β-glucosidase (Lc.0.11, Lc.0.3431, Lc.0.157, Lc.0.3975), and one glycogen synthase kinase (Lc.0.1659) showed higher transcript abundance in 'Nuomici' than that in 'Huaizhi', implying that the accumulation of these osmotic regulation substances might be higher in the 'Nuomici' aril than that in the 'Huaizhi' aril. In addition, one potassium transporter gene (Lc.9.1824), one GDPmannose transporter gene (Lc.11.1938), four sugar transporter genes (Lc.8.604, Lc.10.1667, Lc.0.1911, Lc.12.1579), and one organic cation/carnitine transporter gene (Lc.3.650) also showed a higher expression level in 'Nuomici'. Previous studies have shown that the potassium transporter, GDP-mannose transporter, and organic cation/carnitine transporter (OCT) are involved in the transport of key osmotic regulation substances such as potassium, mannose, sucrose, and carnitine, respectively [45][46][47][48][49][50][51][52][53]. Taken together, these findings suggest that the osmotic potential might be lower in the 'Nuomici' aril than that in the 'Huaizhi' aril. Aquaporin (PIP) is a protein located on the cell membrane (inner membrane protein) and forms a "pore" to control water transport through the cell membrane [54]. During strawberry fruit ripening, overexpression of PIP can result in a fast water absorption in fruit by promoting water permeability [55]. Our RNA-Seq and qRT-PCR analysis showed that in one PIP gene (Lc.2.583), accumulation gradually increased and was higher in 'Nuomici' during fruit development, and showed a maximum expression difference between the 'Nuomici' aril and 'Huaizhi' aril during fruit cracking. In addition, two other PIP genes (Lc.0.806, Lc.2.3184) also had a higher expression level in the 'Nuomici' aril during fruit cracking. Based on these findings, we speculated that higher expression of these three PIP genes in 'Nuomici' might lead to stronger water absorption potential, so fruit cracking will be aggravated when fruit encounter a heavy rain [56]. The accumulation of auxin (IAA) and abscisic acid (ABA) in aril is necessary for the growth and development of the aril, but a higher content of IAA or ABA can induce fruit cracking in litchi [57]. ILR1 could release free IAA by cleaving IAA conjugates [58]. In this study, we found that the transcript abundance of one IAA-amino acid hydrolase (ILR, Lc.0.190) was higher in the 'Nuomici' aril during fruit cracking than that in the 'Huaizhi' aril, suggesting that the level of free auxin might be higher in the 'Nuomici' aril than that in 'Huaizhi'. Consistently, the transcript abundance of genes related to auxin signaling including Auxin response protein (Lc.4.1203, Lc.9.967, Lc.4.1203), Auxin-induced protein (Lc.14.358), Small Auxin Up-Regulated genes (Lc.10.123), ARF (Lc.10.2141), and Auxin-regulated protein (Lc.8.1150) was higher in the 'Nuomici' aril than that in 'Huaizhi'. Zeaxanthin epoxidase (ZEP) can produce the precursor of ABA biosynthesis via catalyzing carotenoids in plant [59]. In this study, we found that one ZEP (Lc.0.938) was expressed higher in the 'Nuomici' aril than in the 'Huaizhi' aril, implying that the ABA content might be higher in the 'Nuomici' aril that in the 'Huaizhi' aril.
Transcription factors (TFs) play critical roles in the regulation of plant growth and development through modulating target genes transcription that are involved in specific plant life. To date, there is no TF that has been proven to be involved in fruit cracking directly. Here, we found that genes encoding several families of TFs such as WRKY (Lc.12.1749, Lc.5.1288, Lc.2.1691), MYB (Lc.9.2085, Lc.1.2165, Lc.14.257), and LOB (Lc.7.1393, Lc.1.1523, Lc.7.1160) were differentially expressed during the fruit development and cracking between 'Nuomici' and 'Huaizhi'. It is reported that overexpression of grape vvWRKY11 in Arabidopsis can strengthen the cold resistance by adjusting the osmotic potential [60]. In strawberry, upregulation of FaMYB44.2 can promote the accumulation of sugar and acid in fruit [61]. In banana, MaLBDs are involved in regulating fruit ripening through transcriptional activation of the EXPANSIN expression [62]. However, how these TFs are involved in fruit cracking in litchi needs to be further elucidated.
Based on the expression patterns of these candidate genes and in combination with the "ball-skin versus bladder" theory, we proposed a model of molecular mechanism underlying the fruit cracking in 'Nuomici' (Figure 11). During the fruit growth, the transcript level of genes related to growth-promoted hormones IAA, GA and ETH (Gretchen was high in aril, which could enhance the osmotic potential. It is of note that the transcripts of PIP were also high in aril, which could promote the water transport into aril when fruit encounter a heavy rain, thereby accelerating the fruit cracking rate. Taken together, we propose that these DEGs abovementioned could weaken the pericarp strength, but enhance the aril expanding pressure, thereby making the 'Nuomici' fruit susceptible to cracking.  10.166710. , sugar transporter: Lc.0.1911 was high in aril, which could enhance the osmotic potential. It is of note that the transcripts of PIP were also high in aril, which could promote the water transport into aril when fruit encounter a heavy rain, thereby accelerating the fruit cracking rate. Taken together, we propose that these DEGs above-mentioned could weaken the pericarp strength, but enhance the aril expanding pressure, thereby making the 'Nuomici' fruit susceptible to cracking. Figure 11. A preliminary gene network involved in litchi fruit cracking. In pericarp, the expression level of genes related to hormones IAA, GA, and ETH is low, which might directly repress the pericarp growth or indirectly downregulate specific transcription factors to repress the expression of growth-promoted genes, as a consequence, the pericarp might provide a limited space for aril expanding. In addition, the transcripts of genes related to calcium transport and signaling, and wax synthesis is also low in pericarp, while the transcript level of genes related to cell wall remodeling is high in pericarp. As a result, the pericarp of 'Nuomici' might be both developmental retardation and weak mechanical strength. During the aril growth, the transcripts of genes related to IAA and ABA is high, which might directly promote aril expanding or indirectly upregulate transcription factors such as WRKY, bHLH, DOF, NAC, and MYB to increase growth-promoted gene expression. Additionally, the transcript level of genes involved in starch/sucrose metabolism and transport is high, which could cause high osmotic potential. Notably, the expression level of one PIP is also high in aril, which together imply that the fruit cracking could be accelerated when fruit encounter a heavy rain since the aril might possess a strong capacity of water absorption to enhance the aril expanding pressure. Collectively, these DEGs might weaken the pericarp strength, but enhance the aril expanding pressure, thereby making the 'Nuomici' fruit susceptible to cracking.

Sample Collection and Data Statistics
Three 28 year-old 'Huaizhi' litchi trees and three 28 year-old 'Nuomici' litchi trees Figure 11. A preliminary gene network involved in litchi fruit cracking. In pericarp, the expression level of genes related to hormones IAA, GA, and ETH is low, which might directly repress the pericarp growth or indirectly downregulate specific transcription factors to repress the expression of growth-promoted genes, as a consequence, the pericarp might provide a limited space for aril expanding. In addition, the transcripts of genes related to calcium transport and signaling, and wax synthesis is also low in pericarp, while the transcript level of genes related to cell wall remodeling is high in pericarp. As a result, the pericarp of 'Nuomici' might be both developmental retardation and weak mechanical strength. During the aril growth, the transcripts of genes related to IAA and ABA is high, which might directly promote aril expanding or indirectly upregulate transcription factors such as WRKY, bHLH, DOF, NAC, and MYB to increase growth-promoted gene expression. Additionally, the transcript level of genes involved in starch/sucrose metabolism and transport is high, which could cause high osmotic potential. Notably, the expression level of one PIP is also high in aril, which together imply that the fruit cracking could be accelerated when fruit encounter a heavy rain since the aril might possess a strong capacity of water absorption to enhance the aril expanding pressure. Collectively, these DEGs might weaken the pericarp strength, but enhance the aril expanding pressure, thereby making the 'Nuomici' fruit susceptible to cracking.

Sample Collection and Data Statistics
Three 28 year-old 'Huaizhi' litchi trees and three 28 year-old 'Nuomici' litchi trees with similar initial fruit bearing were randomly selected at Xili Orchard in 2018 (Shengzhen, China). At 13 days after anthesis, twenty fruit bearing shoots located in different directions from each tree were tagged. The tagged shoots were used to monitor the fruit cracking rate and the remaining shoots were used for sample collection. Cracked fruit on tagged shoots were removed after each statistic. Pericarp samples were collected at 13,26,40,47,55,62, and 69 days after anthesis, while aril samples were harvested beginning at 55, 62, 69 days after anthesis. All samples were immediately frozen in liquid nitrogen and stored at −80 • C for future analysis. Each tree was treated as a biological replicate.

cDNA Library Construction and Illumina Sequencing
A total of 100 mg of frozen samples of pericarp or aril were ground into fine powders in liquid nitrogen and total RNA was extracted using RNAprep Pure Plant Kit (TIANGEN, Beijing, China) according to the manufacturer's instructions. Each RNA sample was subjected to DNase digestion (TaKaRa, Dalian, China) to remove any remaining DNA. The RNA content was quantified by spectrophotometry (BioPhotometer plus, Eppendoff, Germany) and checked on 1.2% denaturing agarose gels. The RIN (RNA integrity number) values (>8.0) of these samples were assessed using an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). The construction of the libraries and the RNA-Seq was performed by the Biomarker Biotechnology Corporation (Beijing, China). All procedures for cDNA library construction were performed via the standard Illumina sample preparation protocol. Sequencing of the purified libraries was carried out on an IlluminaHiseqTM2500 (Illumina, San Diego, CA, USA).

Sequence Data Analysis
Clean data were obtained by removing reads containing adapter, poly-N, and low quality reads from raw data. The Q30 and GC content of the clean data were calculated and evaluated. Clean paired end reads were mapped to the reference genome using Hisat software [63]. Read counts and normalized gene expression levels of FPKM values were obtained using Stringtie [64]. Then, all of the genes were searched against the NR, Swiss-Prot, COG, and KEGG annotation.

Screening of Differentially Expressed Genes (DEGs) Related to Fruit Cracking
First, co-expression network analysis was used to screen gene sets related to fruit cracking in litchi. The FPKM values between 'Nuomici' and 'Huaizhi' in pericarp or aril at the same time point were compared and DEGs were extracted according to differential expression greater than two-fold (FDR < 0.01, log 2 > 1 or < −1). After that, the DEGs in pericarp at each time point were combined and defined as P1, and DEGs in aril were defined as A1. Then, the genes in P1 and A1 were divided into different modules using WGCNA with a default parameter setting, respectively. Finally, the WGCNA package was used to calculate the correlation and significance among different gene modules and DCR. The DEGs in modules of pericarp and aril that had significant correlation with DCR were defined as PM and AM, respectively. Meanwhile, the DEGs between 'Nuomici' and 'Huaizhi' in pericarp and aril during fruit cracking were also screened and recorded as P2 and A2, respectively. Next, the DEGs shared by AM and A2, PM and P2 were identified as candidate genes that are involved in the fruit cracking of 'Nuomici'.
In addition, the SOTA cluster analysis was performed on the DEGs shared by AM and A2, PM and P2 using MeV software (http://mev.tm4.org), TreeView (http://jtreeview. sourceforge.net), which was used to visualize the clustering of genes [65,66].

Quantitative RT-PCR (qRT-PCR) Validation
RNA was reverse-transcribed into cDNA by TransScript ® One-Step gDNA Removal and cDNA Synthesis SuperMix (TransGen, Beijing, China) with gDNA Remover. Q-PCR was performed on a LightCycler 480II (Roche, Basel, Switzerland) using the SYBR Green Realtime PCR Master Mix (TOYOBO, Osaka, Japan) as the readout. The running parameters of the PCR amplifications included the following conditions: 95 • C for denaturation at 30 s, followed by 40 cycles of denaturation at 95 • C for 5 s, annealing at 60 • C for 30 s, and extension at 72 • C for 30 s. Then, PCR products were analyzed by the melting curve to confirm the specificity of amplification. Primers were designed using Primer5 and are listed in Table S1. In this study, LcEF1α-normalized expressions are presented in Figure 9, which were reproducible after LcActin normalization [67]. Relative gene expression levels were calculated using the 2 − CT method with three independent biological replicates [68]. Primers used here are listed in Table S1 and the qPCR melting curves are shown in Figure S1.

Data Deposition
All the raw read sequences were deposited in the NCBI (National Center for Biotechnology Information) sequence read archive under BioProject accession PRJNA681070.