The experimental materials were collected from the National Jujube Germplasm Repository, Research Institute of Pomology, Shanxi Agricultural University (Taigu, Shanxi, China), the latitude is N37°20’, and the longitude is E112°29’. We selected two important species of the jujube genus, Ziziphus jujuba ‘Hupingzao’ (HPZ) and Ziziphus spinosa ‘Taigusuanzao’ (TGSZ). The age of the test trees was 30 years; the trees were healthy and showed normal growth under conventional water and fertilizer management conditions. According to the description of the maturity period in the jujube germplasm resources (Li, 2006), fruits were collected from two harvest stages: HPZ white ripening stage (HW) and red stage (HR) as well as TGSZ white ripening stage (SW) and red stage (SR) ( Figure 1 ). Each fruit was similar in size, without mechanical damage or disease. Each sample included fruits weighing 500 g. The selected fruits were immediately frozen in liquid nitrogen and stored at −80°C for metabolomics and transcriptomics analyses, with three biological replicates per group, that is, three trees with the same growth of HPZ and TGSZ were selected for sample collection, and each tree was a biological replicate.

The collected samples were placed in a lyophilizer (Scientz-100F), vacuum freeze–dried, and crushed using an MM 400 pulverizer (Retsch GmbH, Haan, Germany) for 1.5 min at 30 Hz. Next, 50 mg of powder was weighed and dissolved in 1.2 mL of 70% methanol solution for vortex extraction. This process was repeated 6 times for 30 s, followed by incubation for 30 min. The samples were then centrifuged at 12,000 × g for 3 min, after which the supernatant was passed through a 0.22-μm Millipore membrane filter and placed in a liquid sample bottle for ultra-performance liquid chromatography (UPLC) (ExionLC™ AD) and tandem mass spectrometry (MS/MS) analyses.

The data acquisition instrumentation system consisted of UPLC (ExionLC™ AD) and MS/MS (Applied Biosystems 4500 QTRAP).

The UPLC conditions were as follows: Agilent SB-C18 column (2.1 mm × 100 mm; 1.8 µm); mobile phase, 0.1% (v/v) formic acid water (solvent A) and acetonitrile with 0.1% (v/v) formic acid (solvent B); flow rate, 0.35 mL·min⁻¹, injection volume, 5 µL; and column temperature, 40°C. The gradient program was set as follows: 5% B at 0 min, with the proportion of B increasing linearly to 95% within 9 min and remaining constantfor 1 min, after which the proportion decreased to 5% between 10–11 min. The program was equilibrated at 14 min.

The MS conditions were as follows: electrospray ionization source temperature, 550°C; and ion spray voltage, 5500 V (positive mode)/−4500 V (negative mode).Ion source gas I, gas II, and curtain gas were set at 50, 60, and 25.0 psi, respectively, with high values ofthe parameters of collision-induced ionization. QQQ scans were used in multiple reaction monitoring (MRM) mode, and the collision gas (nitrogen) was set to medium. A specific set of MRM ion pairs was monitored during each stage, based on the metabolites eluted during the process.

Based on the Metware database, the material was qualitatively characterized according to the secondary spectrum information. Metabolite quantification was performed using MRM mode analysis of triple four-level rod MS (Wang A. M. et al., 2018). After obtaining the metabolite MS data in different samples, all chromatographic peaks were integrated and the peaks of the same metabolite in different samples were corrected.

To analyze the differences in alkaloids metabolites among the four samples, metabolites with a fold change of ≥2 or ≤0.5 were selected, and metabolites with VIP scores of ≥1 were simultaneously selected (Yuan et al., 2018).

RNA sequencing and assembly were performed by Beijing BioMarker Biological Technology Co., LTD (Beijing, China). Total RNA was extracted from fruits at HW, HR, SW, and SR stages. RNA concentration and quality were assessed using Agilent 2100 instrument. Then, cDNA libraries were constructed for Illumina sequencing. After removing adapters and low-quality reads, we aligned the clean reads to the reference genome and performed differential expression analysis using the DESeq2 method. The screening conditions for DEGs was log2 foldchange of ≥1 and FDR of <0.05.

Data were analyzed using Microsoft Excel v. 2007 (Microsoft Corporation, Redmond, WA, USA). PCA was performed using the Statistical Analysis System v. 9.2 (SAS Institute, Cary, NC, USA). Tbtools (http://doi.org/10.1101/289660) and OriginPro9.0 (OriginLab Corporation, Northampton, MA, USA) were used to construct the figures. The pearson method was used for association analysis of differential alkaloid metabolites and differentially expressed genes, and the values were normalized with log2, with a threshold of 0.8 for association analysis and a p-value threshold of 0.05.

In this study, we analyzed two important species of the jujube genus, Ziziphus jujuba ‘Hupingzao’ (HPZ) and Ziziphus spinosa ‘Taigusuanzao’ (TGSZ). A total of 44 alkaloid metabolites were detected in the following 4 samples: HW, HR, SW, and SR samples, as shown in the Venn diagram ( Figure 2A ; Supplement 1 ). Among these 44 metabolites, 38 were detected in all 4 samples. Three metabolites including N-p-coumaroylagmatine, frangufoline, and scutianineC were detected only in sour jujube fruits at two different developmental stages. Three metabolites (N-acetylputrescine, agmatine, and tryptamine) were not detected in HR samples but were detected in the other three samples. Then, we classified the identified alkaloid metabolites ( Figure 2B ), based on structural differences between metabolites, 44 metabolites could be divided into 5 classes, namely, alkaloids (n = 31), aporphine alkaloids (n = 6), isoquinoline alkaloids (n = 3), phenolamine (n = 2), and plumerane (n = 2).

A comparative group analysis of SW vs. SR, HW vs. HR, HW vs. SW, and SW vs. SR in four samples (HW, HR, SW, and SR) is shown in Figure 3 . A total of 34 DAMs were obtained from the 4 comparison groups ( Supplement 2 ). In the SW vs. SR comparison group, 17 DAMs were identified; compared with SW, 15 DAMs were down regulated and only 2 DAMs (stepharine and shikonin) were up regulated in SR. In the HW vs. HR comparison group, we identified 13 DAMs; compared with HW, 10 DAMs were down regulated and only 3 DAMs (stepharine, isoboldine, and shikonin) were up regulated in HR. Thus, in both Chinese jujube and sour jujube, DAMs were more accumulated in the white ripening stage than in red stage. These metabolites mainly included N-feruloylagmatine, N-acetylputrescine, jubanineA, frangulanine, cocamidopropylbetaine, choline, and agmatine. Conversely, the accumulation of stepharine and shikonin in Chinese jujube and sour jujube was higher in the red stage than in the white stage. In the HW vs. SW comparison group, 21 DAMs were identified; compared with HW, 16 DAMs were up regulated and 5 were down regulated in SW. In the HR vs. SR comparison group, we identified 24 DAMs; compared with HR, 14 DAMs were up regulated and 10 DAMs were down regulated in SR. This also indicates that the accumulation of DAMs in sour jujube was higher than that in Chinese jujube; moreover, the accumulation was higher in the white ripening stage than in the red stage. These metabolites mainly include 3-hydroxypropylpalmitate glc-glucosamine, acetylcholine, Fer-agmatine, frangufoline, jubanineA, N-feruloylagmatine, N-p-coumaroyl agmatine, scutianineC, serotonin, and tryptamine. However, the accumulation of DAMs such as 6-deoxyfagomine, lotusine, and nuciferine was higher in Chinese jujube than in sour jujube.

By further analyzing the characteristics of DAMs in the different groups ( Figure 4 ), we found that there were one to six metabolites that differed in one or more comparison groups; in particular, two metabolites, jubanine A and N-feruloylagmatine, differed in all comparison groups. In the SW vs. SR, HW vs. SW HR vs. SR comparison groups, there were one, two, and six DAMs with single differences, respectively ( Table 1 ), but there were no single differential metabolites in the HW vs. HR comparison group.

Principal component analysis (PCA) was performed using 34 DAMs; two principal components, PC1 and PC2, explained 93.26% of the total variation ( Figure 5A ). The first principal component (PC 1) explained 74.60% of the total variation, and large positive values were associated with the relative content of N-acetylputrescine and agmatine, suggesting that these two metabolites contributed significantly to PC1.The second principal component (PC2) contributed 18.66% of the total variation, and large positive values were associated with the relative content of N-nornuciforine and yuzirine, suggesting that these two metabolites contributed significantly to PC 2. The white ripening stage and red stage of Chinese jujube and sour jujube were separated in the PCA scatterplot, mainly due to their differences in PC1 and PC2.

In scatter plots derived from the PCA of the DAMs based on the four samples, PC1 and PC2 together explained 95.74% of the total variation ( Figure 5B ). PC1 and PC2 explained 77.25% and 18.49% of the total variation respectively. Moreover, HW and SW positively contributed to PC1, HR positively contributed to PC2. The PC1 of N-nornuciforine, yuzirine, N-acetylputrescine, (-)-caaverine, and agmatine were higher and indicated higher contents of these five alkaloid metabolites in HW and SW. The PC2 of N-nornuciforine and yuzirine was higher and indicated higher content of these two metabolites in HR.

To further investigate the change characteristics of DAMs during different harvest stages in Chinese jujube and sour jujube, the relative content of alkaloid metabolites in all samples was standardized by z-value according to screening criteria, and analyzed by a k-average algorithm. The DAMs showed six change trends in the four samples (class 1–6; Figure 6 ). The number of metabolites classified in classes 1 and 5 were 10 and 5, respectively and showed the highest content in SW. Among them, jubanine A and N-feruloylagmatinein in class 1 showed evident differences among the four samples. N-acetylputrescineand agmatine in class 5 showed distinct differences in PCA. There were five DAMs in class 2, which showed the highest in HR, among them N-nornuciforine and yuzirine with obvious differences were classified class 2. Further, there were five metabolites in class 3, with descending order of concentration of HW > HR > SW > SR, indicating higher content of these metabolites in Chinese jujube than in sour jujube. Two metabolites, stepharine and shikonin, were classified as class 4, and were present at high levels in the red stage.

Next, to investigate the metabolism of alkaloids, we performed KEGG pathway enrichment analysis and functional annotation of DAMs ( Figure 7 ). The DAMs in the 4 comparison groups were enriched in 12 metabolic pathways, including glycine, serine, and threonine metabolism (ko00260); arginine and proline metabolism (ko00330); tryptophan metabolism (ko00380); glycerophospho lipid metabolism (ko00564); one carbon pool by folate (ko00670); nicotinate and nicotinamide metabolism (ko00760); indole alkaloid biosynthesis (ko00901); amino acyl-tRNA biosynthesis (ko00970); metabolic pathways (ko01100); biosynthesis of secondary metabolites (ko01110); carbon metabolism (ko01200); and ABC transporters (ko02010). All comparison groups had DAMs enriched in the ko00330 and ko01100 metabolic pathways, whereas those enriched in the ko00760 metabolic pathway only existed in the HR vs. SR comparison group. DAMs enriched in ko00670, ko00970, and ko01200 metabolic pathways were only present in the SW vs. SR comparison group.

According to the transcriptomics data and KEGG enrichment characteristics of DAMs, 11 potential pathways involved in alkaloid synthesis and metabolism were screened out, and we analyzed the genes involved in their metabolic pathways ( Figure 8 ). The 11 pathways were as follows: glycine, serine, and threonine metabolism (ko00260); arginine and proline metabolism (ko00330); tryptophan metabolism (ko00380); glycerophospho lipid metabolism (ko00564); one carbon pool by folate (ko00670); nicotinate and nicotinamide metabolism (ko00760), isoquinoline alkaloid biosynthesis (ko00950); tropane, piperidine, and pyridine alkaloid biosynthesis (ko00960); amino acyl-tRNA biosynthesis (ko00970); carbon metabolism (ko01200); and ABC transporters metabolism (ko02010). Overall, there were between one and 224 genes involved in one or more metabolic pathways. Most genes (n = 285) were enriched in ko01200 pathway, of which 224 genes were only involved in this pathway. Further, 98, 33, and 21 genes were enriched in ko00564, ko02010, and ko00760 pathways, respectively, and these genes were only involved in these enriched metabolic pathways.

We analyzed the genes enriched in the metabolic pathway of alkaloids and screened the DEGs in four samples. A total of 259 DEGs were identified ( Figure 9A ), among which 251 were expressed in all samples; no DEGs were expressed in only one sample. Gene 30628 was expressed in SW and SR samples, indicating that it was only expressed in sour jujube; conversely, Ziziphus_jujuba_newGene_5605 was expressed in HW and HR, indicating that it was only expressed in Chinese jujube. Two genes, gene32374 and gene28293, were expressed in HW and SW samples, indicating that they were only expressed in the white ripening stage.

By analyzing the DEGs in different comparison groups, we found that one to 98 genes were different in one or more comparison groups ( Figure 9B ). Two genes, gene19198 and gene15104, were differentially expressed in all groups. Some DEGs only existed in a single comparison group. The highest number of DEGs was observed in HW vs. SW group (n = 52), whereas the lowest number of DEGs was noted in HR vs. SR group (n = 5). Further, the number of DEGs in the two comparison groups, SW vs. SR (32) and HR vs. SR (43), was relatively small ( Figure 9C ), whereas the other two comparison groups had a large number of DEGs, with 177 and 188 DEGs in HW vs. HR and HW vs. SW groups, respectively. Simultaneously, we analyzed the up regulation and down regulation of DEGs in each comparison group. The number of up regulated and down regulated genes showed a small difference in SW vs. SR and HW vs. HR groups, whereas it displayed a large difference in the HR vs. SR and HW vs. SW groups.

To analyze the expression trends of DEGs in the different samples, we performed k-average algorithm analysis. The DEGs showed six change trends in the four samples class1–6 ( Figure 10 ). The number of genes in each class was 25, 27, 28, 46, 36, and 97, respectively. The expression of DEGs in class1 was highest in SW, and in class3 was highest in HR. The expression levels of DEGs in HR, SW, and SR in classes 5 and 6 were similar. However, the expression of DEGs in class5 was lowest in HW, and in class6 was highest in HW. The expression levels of class2 genes were highest in HW and HR and lowest in SW and SR. The genes in class4 had the lowest expression in HW and the highest expression in SR.

We conducted correlation analysis of the important screened alkaloid metabolites and DEGs, revealing that each metabolite had multiple candidate genes that might be involved in regulating its metabolism ( Figure 11 ). In total, 5 genes were associated with the change in shikonin (pmf0557) metabolism, and 36 genes were associated with the change in stepharine (P2126) metabolism, all of which were involved in the metabolism of a single substance. Six genes were found to be involved in the synthesis and metabolism of N-acetylputrescine (pme2693) and agmatine (pmb0501), three genes were involved in the synthesis and metabolism of jubanine A (N5961) and N-feruloylagmatine (pmb0496), and five genes were involved in the synthesis and metabolism of N-nornuciforine (P3309) and yuzirine (Wmzzp003295).

Finally, we constructed an expression heat map of DEGs related to important DAMs ( Supplement 3 ; Figure 12 ). Based on the differences in their expression levels among different groups, 11 candidate genes, including gene10804, gene12087, gene16783, gene18355, gene21879, gene22155, gene30225, gene30346, gene32374, gene5151, and Ziziphus_jujuba_newGene_1074, were screened out. According to the gene annotation information, the 11 candidate genes were annotated as 1 UGT74B1, 1 frmA, 1 pckA, 1 mmsA, 1 gcvH, 1 TR1, 2 cysK, 1 pfkA, 1 ALDH, and 1 accC gene, respectively. Simultaneously, KEGG pathway annotation demonstrated that these genes were involved in the pathways related to alkaloid metabolism.