Chinese jujube (Ziziphus jujuba Mill.) is one of the global oldest and most extensively cultivated economic tree species (Qu and Wang, 1993), originating in China and exhibiting a cultivation history of over 3000 years. Its cultivation extends across 47 countries, spanning over 2 million hectares (Rashwan et al., 2020). They contain high hardiness and water use efficiency properties that favour cultivation in mountainous and desert areas, resulting in nutrient-rich properties and robust agronomic traits that have made them a globally popular health food (Muhammad et al., 2024). It has a unique flavour due to its amino acids, vitamins, triterpenoids and flavonoids (Zhang et al., 2021; Sapkota et al., 2023) and is widely used in the food processing industry and medical research due to its oxidative, anticancer, antimicrobial, anti-inflammatory, hepatoprotective, protective and gut-protective properties (Sobhani et al., 2020; Agrawal et al., 2023; Rahbardar et al., 2024).
Jujubes are among China’s highest-exported dried fruits, with the dried jujube market accounting for approximately 60% of the market share, while fresh jujubes represent only 10%. Unlike dried jujube fruits, fresh jujubes have limited storage resistance, and their taste, texture and aroma can change rapidly over a short period. These factors significantly impact fruit quality. Therefore, investigating the developmental maturation of fresh jujube fruits enhances our understanding of their quality. Previous physiological studies on fruit development have demonstrated that amino acids, fatty acids and reducing sugars contribute to various metabolic pathways and chemical reactions that influence fruit quality, including sweetness, acidity and other beneficial compounds (Song et al., 2019). Research on Z. jujuba cv. ‘Dongzao’ fruit development demonstrates an S-shaped growth curve for both the longitudinal and transverse diameters of the fruit, as well as single fruit quality. Sugar accumulation coincides with cell expansion, while the contents of acid, sugar and vitamin C show an initial increase followed by a decline, with malic acid and citric acid being the predominant organic acids (Li et al., 2021).
In recent years, reports on the transcriptome and genome data analysis of jujube have been published. Various scholars have conducted extensive research and exploration on key genes that influence fruit development and ripening, particularly focussing on cell wall metabolism during the ripening process of jujube fruits (Hou et al., 2018), the mechanism of accumulation of sugar and acid (Lu et al., 2022b), the gene network associated with ascorbic acid (AsA) biosynthesis and its regulatory mechanisms (Lu et al., 2022a), etc. which have gained a relevant study on the pivotal genes involved. Nodes with the highest degree in gene networks, known as ‘hubs’, are considered to have specific roles within the network (Sohrabi et al., 2023). In this study, we compared the transcriptomes of jujube fruit at various developmental stages and identified differentially expressed genes (DEGs), followed by functional enrichment analyses. The investigation of hub genes across different developmental stages provided a transcript-level reference, clarifying the dynamic developmental processes of Z. jujuba cv. ‘Dongzao’ fruits.
In this study, the jujube (Z. jujuba cv. ‘Dongzao’) served as the experimental material, sourced from a demonstration garden for fresh jujube production in Yangtake Township, Maigaiti, China (latitude: 38°25′-39°22′N, longitude: 77°28′-79°O5′E; altitude: 1178 m). Fruits were collected at four growth and developmental stages: fruit expanding (FE) stage, white-mature (WM) stage, early red (ER) stage and full red (FR) stage (Figure 1). After collection, the fruits were refrigerated and transported to the laboratory. Upon assessment of texture, sensory attributes and physiological characteristics, the peeled pulp was diced, thoroughly mixed and flash-frozen with liquid nitrogen, which was then preserved at -80°C until the extraction of total RNA and metabolite (Savoi et al., 2019).
Photographs of Z. jujuba cv. 'Dongzao' at four different stages of development (63 days, 78 days, 98 days, 108 days after flower blooming corresponding to the expansion stage, the FE stage, WM stage, ER stage, and FR stage, respectively). ER, early red; FE, fruit expanding; FR, full red; WM, white-mature.
The extraction of total RNA was performed from frozen jujube flesh by utilising TRIzol reagent (Thermo Fisher, 15596018), according to the instructions of the manufacturer. RNA-seq and bioinformatics analyses were conducted by LC-Bio Technology Co., Ltd. (Hangzhou, China). The quality and quantity of the RNA samples were assessed by employing a Bioanalyzer 2100 and an RNA 6000 Nano LabChip Kit (Agilent, CA, USA, 5067-1511). Library construction was realised using high-quality RNA samples with RIN values exceeding 7.0, with subsequent sequencing performed on the Illumina NovaseqTM 6000 platform.
High-quality clean data was gained by removing Reads containing joints and low-quality Reads (including Reads with a proportion of N > 10% and Reads with a number of bases of quality Q ≤ 10 accounting for >50% of the entire Read). Quality score and base contents of clean reads were counted. The HISAT2 (version 2.2.1) (Kim et al., 2019) software was utilised to compare clean reads to the reference genome (Z. jujuba var. ‘Spinosa’. GCA_020796205.1.genome. fa. Isolate: AT0. Assembly type: Haploid. Number of chromosomes: 12. Number of organelles: 2) to obtain information on the positioning of the clean reads on the reference genome (Kim et al., 2015). Subsequently, the alignment results of sequencing data were counted using the ‘BWA’ R package (version 1.9) (Li and Durbin, 2009). Variance testing of the comparison results was carried out via ‘GATK’ R package (version 4.3.0.0), containing single nucleotide polymorphism (SNP) and insertion-deletion (INDEL) (Li, 2014). String Tie software (version 2.2.1) was utilised to conduct Mapped Reads splicing and comparison with the original genome annotation information to discover new transcripts and unannotated novel genes of the species (Pertea et al., 2015). Sequence alignment of the discovered novel genes with NR, Swiss-Prot, gene ontology (GO), COG, KOG, Kyoto encyclopedia of genes and genomes (KEGG) databases was performed using DIAMOND software (version 2.0.15) (Buchfink et al., 2015). The amino acid sequence of the unannotated novel gene was then compared with the Pfam database to obtain the annotation information of the unannotated novel gene through HMMER software (version 3.3.2) (Ashburner et al., 2000). Finally, based on the reads count mapped to this gene and length of the gene, the fragments per kilobase of transcript per million fragments mapped (FPKM) of each gene was calculated.
The mRNA expression levels were compared via DESeq2 R package (version 1.34.0) in WM vs FE groups, ER vs WM groups, FR vs ER groups and FR vs FE groups (adj P<0.05 and |log2FC| >1) (Love et al., 2014). Furthermore, we performed the enrichment analysis of DEGs in different stages through ‘clusterProfiler’ R package (version 4.2.2) (adj P<0.05), including GO, KEGG and gene set enrichment analysis (GSEA) (Subramanian et al., 2005; Yu et al., 2012).
Initially, the correlations among DEGs within each group were analysed by applying the ‘stats’ R package (version 4.1.0). Genes with |r| >0.95 were selected for further analysis. Next, Cytoscape (version 3.8.2) was imported to construct a relational network diagram. Subsequently, hub genes for each group were identified using the degree algorithm in Cytohubba (Chin et al., 2014).
To further analyse the hub gene expression, their expression were compared across the four stages using the Kruskal–Wallis test function of the ‘rstatix’ R package (version 0.7.0).
The comparison efficiency of reads with the reference genome ranged from 79.83% to 81.82% for each sample (Supplementary Table S1). Analysis of single-base conversion types revealed relatively high frequencies of A → G, T → C, C → T and G → A mutations, with a consistent trend observed for other mutations (Figure 2A). A small proportion of genes exhibited a SNP density >1 ⋅kb−1 base, indicating a limited number of genes with SNP mutations (Figure 2B). The most frequent type of SNP variant was the downstream gene, followed by the upstream gene and intron mutations (Figure 2C). In terms of INDEL mutations, intron mutations were the most common, followed by downstream mutations (Figure 2D). Additionally, 3448 unannotated novel genes were identified. The logarithmic FPKM values for most of the 12 samples ranged from -2.5 to 2.5 (Figures 1E and 1F).
Global assessment of transcriptome data. (A) Sequencing of single-base conversion types. (B) Distribution density of SNP counts, which shows the quality of transcriptome libraries. (C) SNP annotation categorisation. (D) INDEL-annotated classification. INDEL, insertion-deletion; SNP, single nucleotide polymorphism.
In the WM vs FE groups, 1576 DEGs were identified (Figures 3A and 3B). Enrichment in the GO database showed that these DFGs are mainly concentrated in the pathways of biological processes such as stimulation, defence and immunity (Figure 3C). Specifically involving responses to biotic stimuli, defence responses and hydrolase activity related to hydrolysing O-glycosyl compounds (Figure 3C). Furthermore, the DEGs participated in 106 KEGG pathways, including the MAPK signalling pathway in plants, starch and sucrose metabolism and plant hormone signal transduction (Figure 3D). Additionally, GSEA results revealed enrichment in photosynthesis, nucleotide sugar and amino sugar metabolism (Supplementary Figure S1A).
Transcriptome analysis of fruit at the WM and FE stages. (A) Volcano plot showing the DEGs; the blue dots represent downregulated genes and the red dots represent upregulated genes. (B) Clustering heat maps of the top 10 DEGs, each row represents a gene and each column represents a sample, with the colour bands on the right indicating different gene expressions; the closer to yellow the higher the gene expression, the closer to green the lower the expression. (C) GO functional enrichment analysis of the genes in the FE and WM jujube groups. (D) KEGG functional enrichment analysis of the genes in the FE and WM jujube groups. DEGs, differentially expressed genes; FE, fruit expanding; GO, gene ontology; KEGG, Kyoto encyclopedia of genes and genomes; WM, white-mature.
Figures 4A and 4B illustrate 8561 DEGs observed between the ER and WM groups. These DEGs were enriched in cellular responses to various stimuli, including chemical and endogenous stimuli (Figure 4C). Enrichment analysis revealed 122 KEGG pathways, with alpha-linolenic acid metabolism being the only significant pathway (Figure 4D). Notably, highly expressed genes in ER vs WM were associated with tryptophan metabolism, whereas low-expression genes were linked to plantpathogen interactions (Supplementary Figure S1B).
Transcriptome analysis of fruit at the ER and WM stages. (A) Volcano plot demonstrating the DEGs; blue dots represent downregulated genes and the red dots represent upregulated genes. (B) Clustering heat maps of the top 10 DEGs, each row represents a gene and each column represents a sample, with the colour bands on the right indicating different gene expressions; the closer to the yellow the higher the gene expression, while the closer to the green the lower the expression. (C) GO functional enrichment analysis of genes in the WM and ER jujube groups. (D) KEGG functional enrichment analysis of genes in the WM and ER jujube groups. DEGs, differentially expressed genes; ER, early red; GO, gene ontology; KEGG, Kyoto encyclopedia of genes and genomes; WM, white-mature.
The 1073 DEGs were screened in FR vs ER groups (Figures 5A and 5B). In GO items, the DEGs were related to defence response, DNA-binding transcription factor activity and ADP binding (Figure 5C). In KEGG pathway, DEGs were markedly correlated with plant-pathogen interaction and alpha-linolenic acid metabolism (adj P<0.05) (Figure 5D).
Transcriptome analysis of fruit at the ER and FR stage. (A) Volcano plot showing the DEGs, the blue dots represent downregulated genes and the red dots represent upregulated genes. (B) Clustering heat maps of the top 10 DEGs, each row represents a gene and each column represents a sample, with the colour bands on the right indicating different gene expressions; the closer to the yellow the higher the gene expression, and the closer to the green the lower the expression. (C) GO functional enrichment analysis of genes in the ER and FR groups of jujuba. (D) KEGG functional enrichment analysis of genes in the ER and FR groups. DEGs, differentially expressed genes; ER, early red; FR, full red; GO, gene ontology; KEGG, Kyoto encyclopedia of genes and genomes.
Furthermore, 7884 DEGs were observed in the FR vs FE groups (Figures 6A and 6B). Within the GO category, enrichment was observed in secondary active transmembrane transporter activity (Figure 6C). DEGs were associated with 122 KEGG pathways, primarily involving the metabolism of various acids. However, none of these pathways was significant (adj P< 0.05). The GSEA suggested that highly expressed genes in FR vs FE were enriched in nucleotide excision repair, protein processing in the endoplasmic reticulum and other pathways (Supplementary Figure S1C). Conversely, low-expression genes were enriched in ribosome, plantpathogen interaction pathways and alpha-linolenic acid metabolism (Supplementary Figure S1C).
Transcriptome analysis of fruit at the FE and WM stages. (A) Volcano plot indicating the DEGs; blue dots represent downregulated genes and the red dots represent upregulated genes. (B) Clustering heat maps of the top 10 DEGs, each row represents a gene and each column represents a sample, with the colour bands on the right indicating different gene expressions; the closer to the yellow the higher the gene expression, and the closer to the green the lower the expression. (C) GO functional enrichment analysis of the genes in the FE and WM jujube groups. DEGs, differentially expressed genes; FE, fruit expanding; GO, gene ontology; WM, white-mature.
Five hub genes were identified in each group. As shown in Figures 7A–7D, the correlation between DEGs in each group was shown. In WM vs FE groups, the hub genes were LOC107414312, LOC107432086, LOC107409707, LOC107404967 and LOC107432771 (Figure 7A). Between ER and WM groups, the hub genes were LOC125421296, LOC107420835, LOC107425022, LOC107427546 and LOC107425432 (Figure 7B). In FR vs ER groups, hub genes were LOC107418664, LOC107424818, LOC107412925, LOC107416546 and LOC107415777 (Figure 7C). The hub genes were LOC107410556, LOC107417775, LOC107427813, LOC107414679 and LOC107433779 in FR vs FE groups (Figure 7D). In addition, we assessed the expression of hub genes in four stages of winter jujube. The results suggested that 4 hub genes were not markedly differential among four stages, namely LOC107409707, LOC107416546, LOC107415777 and LOC107414679 and the expression of the remaining hub genes was significantly different (Figure 8).
Gene interaction network diagram. The pink nodes represent non-hub genes obtained, the yellow nodes represent hub genes after correlation coefficient screening and black lines indicate the presence of correlation between genes. The black line indicates the correlation between genes: (A) WM vs FE; (B) ER vs WM; (C) FR vs ER; (D) FR vs FE. ER, early red; FE, fruit expanding; FR, full red; WM, white-mature.
Differential expression boxplot of the four hub genes. (A), LOC107409707; (B), LOC107416546; (C), LOC107415777; (D), LOC107414679.
Jujube (Z. jujuba cv. ‘Dongzao’) stands as a commercially significant fruit tree in China, renowned for its distinct flavour (Qu and Wang, 1993) and holds prominence in both dietary and traditional medicinal practices. Investigating the developmental and ripening processes of fruit plants aids in the understanding of fruit quality (García-Gómez et al., 2020; Nilo-Poyanco et al., 2021). Jujube undergoes diverse physiological and molecular transformations during maturation, necessitating a comprehensive study of its growth process. RNA-seq analysis serves as a vital tool for providing genetic insights into physiological and metabolic changes related to fruit development and ripening in fruit trees (Travisany et al., 2019; Nawaz et al., 2021). In this study, cDNA libraries were constructed and sequenced across four fruit stages: FE, WM, ER and FR, by applying an Illumina HiSeqTM6000 sequencing platform, leading to the discovery of 3448 unannotated novel genes.
Prior research employing RNA-Seq on Rubus idaeus cv ‘Nova’ and ‘Heritage’ highlighted differential expression of a larger number of DEGs associated with hormonal signalling and cell-wall disassembly between the flower and pink fruit stages of raspberries. And the same view was corroborated in our study, where jujube fruit-associated differential genes were shown to be significantly different at different developmental stages (Hyun et al., 2014; Travisany et al., 2019). In this study, the analysis of jujube fruits at various developmental stages revealed significant differences: 1576 DEGs were identified between the white-mature and FEs, 8561 DEGs between the early red and WMs, 1073 DEGs between the full red and ERs and 7884 DEGs between the full red and FEs. These findings indicated the crucial differences between libraries at different fruit developmental stages, facilitating the identification of more DEGs. Particularly noteworthy is the heightened gene expression and metabolic activity observed during the transition from white-mature to ERs, suggesting a period of increased fruit respiratory activity (Abbas and Saggar, 1989).
This study employed KEGG and GO analyses to determine the biological functions of DEGs. GO analysis revealed the metabolic enrichment disparities and associated changes in biological functions manifested by fruit gene expression (Tomczak et al., 2018). There were notable enrichment entries across the comparative periods examined in this study. Biological processes primarily focussed on abiotic stress responses. Additionally, the unmentioned entries indicated associations of these genes with the metabolism of various sugars and amino acids. The development of fruit plants is characterised by a series of internal and external morphological, physiological and transcriptomic changes that act in concert with stages of active cell division, expansion and response to the environment. Transcriptome differential genes showed significant enrichment in relation to abiotic stresses, suggesting a possible transition to enhanced defence between these periods (Ando et al., 2012). Concurrently, KEGG pathway analysis demonstrated DEGs in several pathways, notably enriched in various acid metabolisms, plant hormone signal transduction, starch and sucrose metabolism, alpha-linolenic acid metabolism and the MAPK signalling pathway in plants. These pathways were closely related to fruit development, suggesting a pivotal role for hormone- and sugar-acid-regulated genes in jujube fruit ripening (Wang et al., 2023).
In plants, α-linolenic acid is the main polyunsaturated fatty acid in cell membranes, which plays a role in maintaining membrane integrity (Hernández et al., 2016). Studies by previous scholars have shown that jujube has the highest number of aroma-active compounds associated with α-linolenic acid. This may be because these fatty acids have unsaturated double bonds and are more susceptible to the formation of aromatic compounds through oxidative reactions. (Sahaka et al., 2020; Gou et al., 2023). Fruit maturity is positively correlated with alpha-linolenic acid content, a phenomenon observed in Lycopersicon esculentum (Monaco et al., 2018) and Plukenetia volubilis (Goyal et al., 2022). The α-linolenic acid metabolic pathway has also been studied more, and some scholars believe that α-linolenic acid metabolism affects the fluidity of the cell membrane by adjusting the ratio of unsaturated phospholipids to cope with the damage caused by a variety of stress factors on the cell membrane function (Zhang et al., 2022; Wang et al., 2024). Some scholars also believe that α-linolenic acid metabolic pathway may be an important pathway in Osmanthus fragrans fruit development, in which the up-regulation of LOX, AOS and OPR genes may be the cause of α-linolenic acid accumulation (Ma et al., 2023). Our enrichment analysis revealed significant overexpression of α-linolenic acid metabolism across all fruit developmental stages, which suggests that this may also be an important pathway during fruit development of jujube.
The GSEA, a pivotal tool in genetic research, facilitates the identification of key biological pathways and processes linked to specific phenotypes (Powers et al., 2018). It has been widely employed in the study of plant fruit development processes (García-Gómez et al., 2020; Kim et al., 2021). Our GSEA findings revealed enrichment of photosynthesis, aminosugar and nucleotide sugar metabolism in the early stage, except for FR vs ER. Conversely, tryptophan metabolism and alpha-linolenic acid metabolism were more abundant in the late stage than in the early stage. These results partially aligned with those obtained from KEGG analyses, suggesting their significant roles in the fundamental development of jujube fruits.
Hub genes, which are highly connected nodes within a network, often affect the expressions of other genes, making them focal points in studies of critical biological processes (Haque et al., 2019; Sohrabi et al., 2023). In this study, 20 hub genes were identified through correlation coefficients among DEGs across various periods, employing the degree algorithm. Among these, 16 period-related hub genes exhibited significant differences (P< 0.05). For instance, the ZjPME15 gene (LOC107414312) presented marked variation between the FE and the WM. Pectin, a principal constituent of plant cell walls, was modified by pectin methyl esterase (PME), which affected cell wall elasticity and permeability, thereby crucially influencing plant growth and development. Notably, in the jujube variety ‘Dongzao’, protopectin content in fruit cell walls decreased from the white ripening stage to the full-red stage (Ng et al., 2013), suggesting a potential involvement of LOC107414312 in early fruit development through cell wall modification (Wang et al., 2021). Thus, we posited that LOC107414312 may be associated with jujube fruit cell wall metabolism. Furthermore, the involvement of certain hub genes in the developmental stages of jujube fruits has not been investigated. Various studies have highlighted their significant roles in regulating plant growth and development (Pecenková et al., 2011; Bavnhøj et al., 2021; Fujihara et al., 2021), with some being implicated in studies on plant stress tolerance (Coleman et al., 2021).
De novo assembly of transcripts and Illumina HiSeq™ 6000 sequencing, combined with KEGG term analysis and hub gene identification, provides valuable insights into candidate genes and DEGs involved in jujube fruit development and ripening. However, our study represents only an initial screening of hub genes across the four stages of jujube fruit growth. Further research is necessary to explore the precise regulatory mechanisms and functions of these genes. We are dedicated to investigating the roles of these genes in future studies.