‘Katy’ and ‘Kuijin’ apricot fruits were obtained from the Taidong experimental orchard of the Shandong Institute of Pomology (Taian, Shandong, China). Uniformly developing fruits at the same stage were collected from three disease- and insect-free per cultivar. Samples were collected during the young fruit period (42 days after full bloom; A), the color-changing period (52 days after full bloom; B), and the mature period (84 days after full bloom; F) ( Figure 1 ). Thus, six samples (KatyA, KatyB, KatyF, KuijinA, KuijinB, and KuijinF) with three biological replicates each (18 samples in total) were collected in 2018. The collected fruits were immediately frozen in liquid nitrogen and stored at −80°C until analyzed.

Freeze-dried pulp samples were crushed using zirconia beads and the MM 400 mixer mill (30 Hz for 1.5 min) (Retsch, Haan, Germany). The resulting powder was weighed and 100 mg was mixed with 1.0 mL 70% aqueous methanol and then incubated overnight at 4°C. After centrifuging the sample at 10,000 g for 10 min, thesupernatant was added to the CNWBOND Carbon-GCB SPE Cartridge (250 mg, 3 mL; ANPEL, Shanghai, China) and filtered (SCAA-104, 0.22 μm pores; ANPEL). The metabolites were quantified using a triple quadrupole mass spectrometry system in the multiple reaction monitoring mode and then identified using the MetWare database and other publicly available metabolite databases. The Analyst (v1.6.1) software (AB SCIEX, ON, Canada) was used to process the mass spectrometry data. After obtaining the mass spectra of the metabolites in different samples, the peak areas of all mass spectral peaks were integrated and the mass spectral peaks of the same metabolite in different samples were integrated and corrected (Fraga et al., 2010). Furthermore, a principal component analysis (PCA) was performed to analyze the metabolite data. A partial least squares discriminant analysis (PLS-DA) was completed to maximize the differences in metabolites between groups and identify the metabolites. The variable importance in projection (VIP) value of the OPLS-DA model was combined with the P-value or fold-change data in the univariate analysis to screen for differentially metabolites. Metabolites with a VIP ≥1 and a fold-change ≥2 or ≤0.5 were considered to be significantly different.

TRIzol reagent (Invitrogen, China) was used to extract total RNA from the pulp of ‘Kuijin’ and ‘Katy’ apricot fruits collected at three developmental stages (A, B, and F). The integrity of the extracted RNA and the presence of contaminating DNA were analyzed by agarose gel electrophoresis and the NanoPhotometer system (Implen GmbH, Munich, Germany). The 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) was used to verify the purity and quantity of the RNA. Triplicate samples of pulp at each developmental stage were used to construct 18 apricot transcriptome libraries, which were sequenced using the Illumina HiSeq™ 2000 system at the Beijing Genomics Institute to produce 150-bp paired-end reads. Raw sequencing reads were filtered by removing reads containing adapters or more than 10% unknown nucleotides (N) and low-quality reads [i.e., more than 50% low quality (Q <20) bases] using FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/ ). The clean reads were mapped to the apricot reference genome (https://www.ncbi.nlm.nih.gov/genome/?term=apricot ) using HISAT2 (v2.0.5) (Kim et al., 2013).

Gene expression levels were calculated on the basis of the fragments per kilobase of transcript per million fragments mapped (FPKM) value. Genes that were differentially expressed between two samples (KatyA vs KatyB, KatyA vs KatyF, KatyB vs KatyF, KuijinA vs KuijinB, KuijinA vs KuijinF, KuijinB vs KuijinF, KuijinA vs KatyA, KuijinB vs KatyB, and KuijinF vs KatyF) were analyzed using the DEGseq R package (Love et al., 2014). The non-standardized raw read count data for the genes were used, and the read count for the genes was determined using featureCounts (Li, 2014). After the differential gene expression analysis, the hypothesis test P-value was corrected according to the Benjamini–Hochberg method to obtain the false discovery rate (FDR). Genes with |log₂(fold-change)| ≥1 and FDR <0.05 were screened as DEGs, which were subsequently functionally annotated via Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG; https://www.genome.jp/kegg) pathway analyses (Yonekura et al., 2012).

Co-expression networks were constructed using the WGCNA (v1.29) package in R (Langfelder and Horvath, 2008). Among the 12593 genes, those with an average FPKM >1 (from three replicates) were used for the WGCNA. The modules were obtained using the automatic network construction function blockwise. The default settings were applied, with the exception of the soft power (18), the min module size (50), and the merge cut height (0.45). The eigengene value was calculated for each module and used to test the associations with the metabolism of flavor-related compounds.

Total RNA was extracted from the samples using the HiPure HP Plant RNA Mini Kit (Magen R4165-02) and then reverse transcribed to cDNA using the SUM Onetube RT Mixture III (gDNA removal) (SUMMER SUM7806a). On the basis of the gene names, the corresponding sequences were obtained and used to design primers. The relative expression levels of 12 candidate genes were determined by qPCR, which was completed using the CFX 96 Real-time PCR assay system (Bio-Rad, Hercules, CA, USA). The SYBR Green dye method was used, and PaACIIN was the internal reference gene. Seven structural genes and five transcription factor genes were candidate genes. The relative gene expression levels were calculated according to the 2^−ΔΔCt method. The qPCR program was as follows: 94°C for 10 min; 40 cycles of 94°C for 10 s, 60°C for 30 s, and 55°C for 30 s. The qPCR analysis was conducted using three biological replicates. Details regarding the qPCR primers are listed in Table S8 .

The SPSS 22.0 software (SPSS Inc., Chicago, IL, USA) was used for the one-way ANOVA with a Tukey test. The threshold for determining significant differences was P ≤ 0.05. Data are provided herein as the mean ± standard deviation.

To analyze the changes in the flavonoid content of apricot fruits at different developmental stages, the total flavonoid contents of ‘Katy’ and ‘Kuijin’ fruits collected in three periods were measured ( Figure 1 ). The flavonoid content decreased from 0.28 mg/g to 0.12 mg/g during the development of ‘Kuijin’ fruits, whereas it decreased from 0.23 mg/g to 0.05 mg/g during the development of ‘Katy’ fruits. The total flavonoid content was significantly higher in the ‘KuijinA’ samples than in the other samples ( Figure 2A ). The changes in the flavonoid metabolites in ‘KuijinA’ were investigated further ( Table S1 ). The PCA results indicated that the biological replicates had similar metabolite profiles, with KuijinA significantly separated from KuijinB, KuijinF, KatyA, KatyB, and KatyF along the principal component 1 (PC1) axis and the PC2 axis ( Figure 2B , Figure S1 ). A total of 572 metabolites were identified in the ‘Katy’ and ‘Kuijin’ fruits at three developmental stages ( Table S2 ). These metabolites, which were divided into 30 groups, included 111 flavonoids (12 catechin derivatives, 16 flavanones, 25 flavones, 9 flavone C-glycosides, 29 flavonols, 1 flavonolignan, 7 isoflavones, 7 anthocyanins, and 5 proanthocyanidins) ( Figure 2C ). Upon comparing the classification of flavonoid compounds, it was found that the levels of flavanol, flavanone, and flavonec-glycosides were significantly higher during the ‘KuijinA’ compared to other stages. Therefore, a further comparison was made to analyze the types and contents of flavonoid compounds during each stage.

After comparing the flavonoid compounds, it was observed that 16, 38, 39, 52, 42, 30, 32, 36, and 50 compounds were differentially accumulated in the comparisons of KatyA vs KatyB, KatyA vs KatyF, KatyB vs KatyF, KuijinA vs KuijinB, KuijinA vs KuijinF, KuijinB vs KuijinF, KuijinA vs KatyA, KuijinB vs KatyB, and KuijinF vs KatyF, respectively as shown in ( Figure 3A ). The result indicates that the group with the highest number of accumulated flavonoid compounds is KuijinA vs KuijinB, followed by KuijinF vs KatyF.

The flavonol species were the most abundant among differential flavonoid species in the nine groups ( Figure 3A ). There were 20 differentially accumulated flavonols in the KuijinA vs KuijinB and KuijinF vs KatyF comparisons, which was more than in the other comparisons ( Figure 3A ). The differentially accumulated flavonols among the samples were identified by analyzing the relative flavonol contents. The KatyA vs KatyB, KatyA vs KatyF, KatyB vs KatyF, KuijinA vs KuijinB, KuijinA vs KuijinF, KuijinB vs KuijinF, KuijinA vs KatyA, KuijinB vs KatyB, and KuijinF vs KatyF comparisons revealed 1, 12, 11, 13, 11, 2, 14, 6, and 19 flavonols with decreasing contents, respectively ( Figure 3B , Table S3 ). Although there are more different accumulate flavonols between KuijinA vs KatyA and KuijinF vs KatyF, ‘KatyA’ shows a greater decrease in flavonol content compared to its increase, while ‘KuijinF’ has more flavonols decreasing than increasing compared to’KuijinA’, indicating that ‘KuijinA’ has a higher flavonol content than both ‘KatyA’ and ‘KuijinF’.

Moreover, the following 10 flavonols accumulated more in the ‘KuijinA’ than in the other samples: quercetin 7-O-rutinoside (pmb0711), quercetin O-acetylhexoside (pmb3026), quercetin 3-O-rutinoside (pme0202), quercetin 3-alpha-L-arabinofuranoside (pme0361), kaempferol 3-O-rutinoside (pme0369), kaempferol 3-O-robinobioside (pme1606), quercetin 4′-O-glucoside (pme3129), quercetin 3-O-glucoside (pme3211), kaempferol 3-O-galactoside(pme3268), and quercetin 7-O-β-D-glucuronide (pme3442)( Figure 3C ). These might be the key metabolites influencing the total flavonol content in the ‘KuijinA’.

The flavonoid content of ‘Kuijin’ in its young fruit stage is higher than other stages, with the greatest variety of flavanols. This is mainly due to the presence of 10 types of flavanols. ( Figure 3A ). The relative contents of 10 flavonols were significantly higher in KuijinA than in the other samples ( Figure 3C ), including quercetin 3-O-rutinoside (rutin), kaempferol 3-O-rutinoside (nicotiflorin), and kaempferol 3-O-robinobioside (biorobin) ( Figure 3C ). Accordingly, these flavonols in apricot fruits are glycosylated at 3-O.

The transcriptomes of the ‘Katy’ and ‘Kuijin’ fruits in three developmental stages were analyzed. The transcriptome sequencing data for 18 apricot samples, the sample correlation heat map ( Figure S2 ), and the PCA plot ( Figure S3 ) revealed biological replicates of the same species are grouped together. A total of 138.84 Gb clean data were obtained, with 108 Gb high-quality reads. The proportion of Q20 bases exceeded 96%, whereas the proportion of Q30 bases was between 89.15% and 91.57%. The GC content was greater than 45% ( Table S4 ).

Different numbers of DEGs were detected in nine pairs of comparisons, including 990 for KatyA vs KatyB, 8,004 for KatyA vs KatyF, 5557 for KatyB vs KatyF, 3952 for KuijinA vs KuijinB, 9023 for KuijinA vs KuijinF, 6543 for KuijinB vs KuijinF, 4005 for KuijinA vs KatyA, 1105 for KuijinB vs KatyB, and 2299 for KuijinF vs KatyF ( Table S5 ).

The grouping with significant differences in flavonol content consists of three pairs: KuijinA vs KuijinF, KuijinA vs KuijinB, and KuijinA vs KatyA ( Figure 2C ). KEGG enrichment analysis revealed that phenylpropanoid biosynthesis, flavonoid biosynthesis, and flavonol biosynthesis were enriched in three comparison groups: KuijinA vs KuijinF, KuijinA vs KuijinB, and KuijinA vs KatyA( Figures 4A–C ). There were 166, 100, and 94 differentially expressed genes (DEGs) related to flavonols in the three comparison groups, respectively ( Figure 4D ). A Venn diagram analysis was performed on the differentially expressed genes of the three groups, which revealed that there were 53 commonly DEGs ( Figure 4D ). Among 53 DEGs, three structural genes strongly correlated to the biosynthesis of flavonols were identified (Pearson correlation coefficients > 0.8, p value< 0.05), including PARG09190, PARG15135, PARG17939. According to the analysis of gene function and RNA-Seq, three structural genes were positively correlated to the flavonols content ( Table S6 ).

To gain further insight into the regulation of flavonols biosynthesis, we carried out WGCNA to investigate the co-expression gene modules and the critical modules involved in flavonols biosynthesis. A total of 8 co-expression modules were identified according to their expression patterns ( Figure 5 ). The correlation between the gene matrix of different modules and the flavonols content was analyzed, and the correlation and corresponding e-value were presented in a digital form in the grid where each module and trait intersect ( Figure 5 ). According to the ‘module character’ correlation analysis, the turquoise module showed a significant positive correlation with 10 flavonols. There were 4897 genes in this module.

Among the 4897 genes in the turquoise module, 66 structural genes were identified involved in phenylpropanoid biosynthesis, flavonoid biosynthesis, and flavonol biosynthesis pathways. To generate the regulatory network associated with flavonols biosynthesis, we constructed the flavonol metabolic pathway in Figure 6A , and examined the structural genes involved in flavonols biosynthesis. Interestingly, 66 DEGs including PARG09190, PARG15135, PARG17939, were recognized both in the biosynthesis of flavonols in the turquoise module, whose expression was highly correlated with the 10 flavonols content ( Figure 6B ), suggested that these DEGs correspond to the regulation of flavonols biosynthesis.

According to weight value, the top ten transcription factors associated with each of the three structural genes were selected from 4897 genes in the turquoise modules. A total of 28 transcription factors were screened ( Figure 6C ). The structural genes and transcription factors were organized into a connection network using Cytoscape software. In the turquoise module, the top 10 transcription factors (TFs) associated with PARG09190 were PARG08234 (AP2/ERF-ERF), PARG23153 (C2C2-LSD), PARG11889 (Tify), PARG08766 (AP2/ERF-ERF), PARG07854 (MYB), PARG10875 (TCP), PARG21139 (MYB), PARG21235 (AP2/ERF-ERF), PARG27864 (AP2/ERF-ERF) and PARG28076 (C2C2-YABBY)( Figure 7A ). The top 10 TFs and transcriptional regulatory factor (TR) associated with PARG15135 were PARG26144 (bHLH), PARG09370 (C3H), PARG27864 (AP2/ERF-ERF), PARG10875 (TCP), PARG27049 (B3), PARG08611 (HB-HD-ZIP), PARG18460 (TRAF), PARG26295 (Others), PARG18567 (bHLH), PARG13130 (MYB)( Figure 7B ). The top 10 TFs associated with PARG17939 were PARG02410 (E2F-DP), PARG02071 (GRAS), PARG01683 (WRKY), PARG01970 (AP2/ERF-ERF), PARG02947 (bHLH), PARG02441 (SBP), PARG00953 (HB-HD-ZIP), PARG01507 (C2C2-YABBY), PARG02307 (WRKY), PARG00355 (bZIP)( Figure 7C ). Interestingly, two common transcription factors, including PARG27864 and PARG10875 associated with PARG09190 and PARG15135 in the module, whose expression was highly correlated with the 10 flavonols content ( Figure 7D ).

To identify potential genes related to the metabolism of 10 flavonol compounds, candidate genes and regulators were analyzed using Cytoscape software. Three structural genes were identified: CYP98A2 (PARG09190), 4CL (PARG15135), and FLS (PARG17939). Based on the weight value, 28 transcription factors are associated with 3 structural genes. Two of the transcription factors are not only associated with PARG09190 but also with PARG15135, indicating their critical importance in the flavonols biosynthesis.

To validate the RNA-seq (i.e., transcriptome sequencing) data and the contribution of specific genes to flavonol synthesis, degradation, and transport, we selected 12 genes for a qPCR analysis. The qPCR data for seven structural genes and five transcription factor genes were consistent with the FPKM values derived from the RNA-seq data ( Figure S4 ).

Metabolites are the final products of cellular regulatory processes, and their levels vary in response to genetic and environmental changes (Fiehn, 2002; García-Gómez et al., 2021). Flavonols are major flavonoids with biological functions (Barbehenn and Constabel, 2011). For example, they participate in auxin transport and responses to oxidative stress (Agati et al., 2012; Brunetti et al., 2013). Previous studies on apricot fruits primarily focused on fruit flavors (e.g., sweetness due to sugars and acidity) and flesh colors (Zhang et al., 2019; Zhou et al., 2020).

In the current study, we observed that for ‘Katy’ and ‘Kuijin’, flavonol glycosides were more abundant in young fruits than in mature fruits ( Figure 2A ). We applied a comprehensive targeted metabolomics approach to analyze the secondary metabolites of ‘Katy’ and ‘Kuijin’ fruits during three developmental stages. A total of 572 metabolites were detected, and the changes in their contents in developing fruits were explored ( Figure 2C ).

Our analyses showed that the flavonoid content was highest in KuijinA, with flavonol glycosides accounting for a considerable proportion of the accumulated flavonoids ( Figure 2C ). According to the molecular structure of the flavonols, eight glycosides were detected, including three monoglycosides (glucose and galactoside), four polyglycosides (rhamnoside, rutinoside, neohesperidoside, and robinobioside), and one glucuronide ( Figure 3C ). These results suggested that the flavonols were mainly glycosylated at 3-O, 5-O, 7-O, and the C-position, which is consistent with the findings of earlier studies (Barreca et al., 2011; Xiao et al., 2014).

Glycosidases catalyze the production of flavonol glycosides, including monoglycosides (e.g., glucosides, galactosides, and rhamnosides), diglycosides [e.g., rutoside (6-O-α-L-rhamnose-D-glucose) and neohesperidin (2-O-α-L-rhamnose-D-glucose)], and polyglycosides (Xiao, 2017; Peterson and Dwyer, 1998).

The mechanisms mediating flavonoid synthesis have been elucidated in many plant species. There has recently been significant progress in the research on the flavonoid synthesis in fruit crops, including grape(Boss et al., 1996), strawberry(Almeida et al, 2007), apple(Takos et al, 2006; Vimolmangkang et al., 2013), and pear(Fischer et al, 2007). The key enzymes involved in flavonol synthesis include flavonol synthase (FLS)(Martens et al., 2002), which belongs to the 2-oxoglutarate-dependent dioxygenase superfamily (Park et al., 2019; García-Gómez et al, 2021). Transgenic tobacco (Nicotiana tabacum) plants expressing OsFLS produce pale pink or white flowers; their petals contain significantly more and less kaempferol-3-O-rutinoside and anthocyanins, respectively, than the wild-type petals (Park et al., 2019). Considering the differences between KuijinA and the other examined samples in terms of their flavonol contents, three DEGs may be critical for flavonol accumulation of which PaFLS (PARG17939), CYP98A2(PARG09190), Pa4CL(PARG15135) may be essential for flavonol production.

In Arabidopsis thaliana, UGT73C6 catalyzes the transfer of glucose from UDP-glucose to the 7-OH position of kaempferol-3-O-rhamnoside and quercetin-3-O-rhamnoside(Jones et al., 2003). In the current study, PaUGT73C (PARG22962) expression was highly correlated with the kaempferol 3-O-rutinoside (nicotiflorin) content. An earlier study confirmed that AtUGT74F1 glycosylates quercetin at the 4′- and 7-hydrogen sites (Lim et al., 2004). In our study, PaUGT74F2 (PARG29757) expression was correlated with the quercetin 7-O-rutinoside content.

The flavonoid pathway is regulated by a highly conserved MYB–bHLH–WD transcriptional complex (Stracke et al., 2001; Wang et al., 2011; Schaart et al., 2013; Xu et al., 2014; Xu and Huang, 2014). In grapes, VvMYBF1 regulates VvFLS1 expression, thereby modulating flavonol biosynthesis (Czemmel et al., 2009). Moreover, AtMYB11, AtMYB12, and AtMYB111 promote flavonol synthesis in A. thaliana (Stracke et al., 2007), while PbMYB12b positively regulates flavonol production in ‘Red Zaosu’ pear fruit (Zhai et al., 2019).

VvMYB86 may negatively regulate anthocyanin synthesis, or positively regulate flavonol synthesis (Sun, 2020). PARG02716 and AtMYB86 belong to the SG13 subfamily. Vimolmangkang et al (2013) overexpression of apple SG4 MYB protein was found MdMYB3 promotes the synthesis of anthocyanin and flavonol in tobacco. PARG14380 and AtMYB6 belong to the SG4 subfamily. Additionally, the AtMYB11, AtMYB12, and AtMYB111 genes reported by Stracke et al. also belong to the SG7 subfamily (Stracke et al., 2007). The three MYB TFs (PARG07517, PARG07518, and PARG27682) were putatived to have similar expression. AtMYB12 was originally identified as a flavonol-specific transcriptional activator in Arabidopsis thaliana, and this has been confirmed by ectopic expression in tobacco. AtMYB12 is able to induce the expression of additional target genes in tobacco, leading to the accumulation of very high levels of flavonols (Zhai et al., 2019).

Bovy et al (2002) expressed the maize transcription factor genes LC and C1 (MYC-type) in the fruits of genetically modified tomato plants to increase flavonoid levels; these fruits accumulated relatively large amounts of kaempferol (flavonol) in their flesh. In our study, the transcription factor genes associated with flavonol accumulation included 12 belonging to the bHLH family, of which the highest kwith value is bHLH106 (PARG26144) ( Table S7 ).

PARG27864 was included in the turquoise module following the WGCNA ( Table S7 ). The AP2/ERF transcription factors reportedly regulate flavonol synthesis in apple (Bovy et al., 2002). AtERF11 - Facilitates flavonol synthesis in celery by directly binding to the PAL gene promoter(Xu et al., 2014; Xu and Huang, 2014). Pti4 - Promotes the expression of several key genes in the flavonol synthesis pathway in tomato (Li et al., 2006). ERF98 - Affects fruit quality and alkaloid content in pomegranate by regulating the expression of genes involved in flavonol synthesis pathway (Zhang et al., 2018).