The red-skinned female parent is Kaorino originated from Japan (https://www.hinshu2.maff.go.jp/vips/cmm/apCMM112.aspx?TOUROKU_NO=19529&LANGUAGE=Japanese). The white-skinned male parent 2012-W02 is a selectively bred seedling derived from the Japanese strawberry variety IROHA-001 (https://www.hinshu2.maff.go.jp/vips/cmm/apCMM112.aspx?TOUROKU_NO=24429&LANGUAGE=Japanese). The pink-skinned accession Fenyu No.1 is the F1 offspring of Kaorino and 2012-W02. All strawberry plants were propagated by strawberry stolons. The stolon seedlings were planted in the greenhouse of Zhijiang Base of Hangzhou Academy of Agricultural Sciences in September 2022, and flowering occurred in mid-October 2022.

Fruit samples were collected from three different plants at 15 days after flowering (white fruit stage, S1), 20 days after flowering (turning stage, S2) and 28 days after flowering (ripening stage, S3). Three biological replicates were performed for each sample. The collected samples were then rapidly frozen in liquid nitrogen and stored at -80°C. These samples were prepared for subsequent metabolomic and transcriptomic analysis.

Metabolite profiling was performed using a widely untargeted metabolome method. The freeze-dried fruit skin from each sample was crushed into powder using a mixer mill (MM 400, Retsch). The lyophilized powder (100 mg) was weighed and extracted with 1.2 mL of 70% methanol. The samples were centrifuged at 12,000 rpm for 10 minutes. After filtering the supernatant, the filtrate was analyzed by UPLC-ESI-MS/MS system (UPLC: Nexera X2 system, Shimadzu, Kyoto, Japan; MS, 4500 Q TRAP, Applied Biosystems, Waltham, MA, USA).

Metabolites were identified according to the m/z values, secondary fragments, and isotopic peaks. Principal component analysis (PCA) was performed to classify and discriminate between samples using the R package pcaMethods. Orthogonal partial least-squares-discriminant analysis (OPLS-DA) (Galindo-Prieto et al., 2015) was performed using the R software package tools (Thevenot, 2016) to calculate the variable influence on projection (VIP) of metabolites. Metabolite levels between different comparison groups were analyzed using paired t-tests and their fold changes were calculated. The p-value of the t-tests was adjusted for false discovery rate (FDR). The differentially regulated metabolites (DRMs) were determined based on the thresholds of VIP>1, log2(fold change)>1 and p ≤ 0.05. The DRMs were then mapped to the KEGG, HMDB and LipidMaps databases for functional and biological process annotation.

The fruit skins (1 cm wide and 0.2 cm thick along the fruit skin) in each sample were obtained. A total of twenty-seven samples (three plants, three fruit development stages and three replicates of each sample) were prepared for RNA extraction. Total RNA was extracted using a PureLink Plant RNA Kit (Invitrogen, Carlsbad, CA, USA), according to the manufacturer’s instructions. An Agilent 2100 Bioanalyzer was used to determine the quality of the RNA, after which the mRNA was purified using poly-T oligo. The library was constructed using the NEBNext Ultra RNA Library Pre Kit (NEB, USA). After quality assessment, the library preparations were sequenced on an Illumina HiSeq 2500 platform, and 150 bp paired-end reads were generated. Clean reads were obtained from the original data by removing reads with adapters, poly-N sequences, and low-quality reads. The clean reads were then mapped to the reference genome (Fragaria x ananassa Royal Royce Genome v1.0 Assembly & Annotation, https://www.rosaceae.org/Analysis/12335030) using HISAT2 software (Kim et al., 2015).

Gene expression levels were analyzed using the fragments per kilobase per million reads (FPKM) method and calculated using Cuffquant and Cuffnorm (v2.2.1) (Trapnell et al., 2010). DEGSeq2 was used to identify DEGs according to the criteria of a log 2-fold change>1 and q ≤ 0.01. Functional annotation of the transcripts was performed using the phyper package in R software. Gene function was annotated based on the following databases: Nr (NCBI non-redundant protein sequences); Nt (NCBI non-redundant nucleotide sequences); Pfam (Protein family); KOG/COG (Clusters of Orthologous Groups of proteins); Swiss-Prot (A manually annotated and reviewed protein sequence database); KO (KEGG Ortholog database); GO (Gene Ontology). Gene Ontology (GO) enrichment analysis of the differentially expressed genes (DEGs) was implemented using the GO seq R packages based on the Wallenius non-central hyper-geometric distribution (Young et al., 2010), which can adjust for gene length bias in DEGs. The KOBAS software (Mao et al., 2005) was used to test the statistical enrichment of DEGs in KEGG pathways.

After obtaining the DEGs and DRMs, we selected the pathways commonly enriched in both DEGs and DRMs for further analysis. The corresponding transcripts and metabolites were filtered based on the annotation information from the KEGG database. Then, the corresponding transcripts and metabolites were mapped to the relevant KEGG pathways.

To evaluate the accuracy of the transcriptome data, we selected eight DEGs, Fxa1Ag101086, Fxa7Cg103039, Fxa7Dg100111, Fxa7Dg102232, Fxa7Dg102859, Fxa4Cg202053, Fxa1Bg201911 and Fxa1Cg102117 for expression validation. qRT-PCR was performed on a Roche-Light Cycler^® 96 system (Roche Diagnostics, Pleasanton, CA, USA). For each reverse transcription reaction, a total volume of 20 µL was prepared containing 10 µL of 2× qPCRBIO SyGreen Mix (PCR Biosystems, London, UK), 1 µL of each of the forward and reverse primers (10 pmoles) and 60 ng/µL of cDNA as template. qRT-PCR was carried out with denaturation at 95°C for 10 min and 45 cycles of amplification with denaturation at 95°C (20 s), annealing at 55°C (20 s) and elongation at 72°C (25 s). Three technical replicates of each of the three biological replicates were tested. The relative expression levels of each gene were analyzed via the 2^–ΔΔCt method (Livak and Schmittgen, 2001). FaACTIN was used as the reference gene for expression analysis (Amil-Ruiz et al., 2013). The primers were designed by using Primer Premier 5 (Premier Biosoft, CA, USA). All primers used for qRT-PCR are listed in Supplementary Table S1 .

Distinct differences in fruit skin color were observed among the red skin accession Kaorino, the white skin accession 2012-W02, and the pink skin accession Fenyu NO.1 at three fruit development stages, especially at the S3 stage ( Figure 1A ). At the S1 stage, the 2012-W02 had light green skin color, while Kaorino and Fenyu NO.1 had white skin color. As the fruit developed, the skin color of 2012-W02 gradually changed to white by the S3 stage. In contrast, the skin color of the Fenyu NO.1 and Kaorino fruits changed from white to pink and finally to red. The skin color of the Kaorino fruits was significantly darker than that of Fenyu NO.1 fruits and appeared bright red at the S3 stage of fruit development ( Figure 1A ).

To further investigate the differences in metabolite accumulation among these three accessions, metabolomic profiles were obtained for each accession at the three different developmental stages. A total of 1,057 metabolites were obtained using a metabolome detection system ( Supplementary Dataset S1 ). The PCA analysis revealed a clear trend of isolation among the three accessions at different developmental stages based on the detected metabolite variables ( Figure 1B ). The first principal component (PCo1, 33.35%) effectively distinguished between the S1, S2 and S3 stages, and the second principal component (PCo2, 13.43%) separated the 2012-W02, Kaorino and Fenyu NO.1 accessions ( Figure 1B ). According to the enrichment of these metabolites in the KEGG database, the top 20 most enriched pathways are shown in Figure 1C . The majority of the identified metabolites were found to be involved in the pathway of purine metabolism, followed by the pathways of pyrimidine metabolism, biosynthesis of phenylpropanoids, biosynthesis of plant secondary metabolites, biosynthesis of various plant secondary metabolites and biosynthesis of plant hormones ( Figure 1C ). In addition, HMDB and LipidMaps annotations were performed for these metabolites ( Supplementary Figure S1 ), and the anthocyanin-related metabolites were then identified by integrating the results of KEGG, HMDB, and LipidMaps annotations. A total of 69 anthocyanin-related metabolites were identified, including one metabolite involved in anthocyanin biosynthesis, 8 metabolites involved in flavone and flavonol biosynthesis, 17 metabolites involved in flavonoid biosynthesis, 29 metabolites involved in phenylpropanoid biosynthesis, and 14 metabolites involved in isoflavonoid biosynthesis ( Figure 2 ; Supplementary Dataset S2 ). All these anthocyanin-related metabolites accumulate differently in three cultivars at different growth stages. For these phenylpropanoid biosynthesis related metabolites, half of the metabolites were accumulated more at S3 stage in three accessions ( Figure 2 ). More p-Coumaroyl quinic acid was found in red skin line Kaorino than in white skin line 2012-W02 at S2 stage, and the same pathway related metabolite, 5-Hydroxyconiferyl alcohol, was also more accumulated in red skin line Kaorino than the other two lines at S1 stage ( Figure 2 ). Likewise, the other anthocyanin-related metabolites involved in flavonoid biosynthesis, flavone and flavonol biosynthesis and isoflavonoid biosynthesis also showed the more accumulated at S3 stage in three accessions. At the S2 developmental stage, the metabolite related to isoflavonoid biosynthesis, Coumestrol, was most abundant in red skin line Kaorino and least abundant in white skin line 2012-W02. Moreover, the metabolites, Peonidin 3-O-glucoside directly participated the pathway of anthocyanin biosynthesis was detected and most accumulated in pink skin line Fenyu NO.2.

A total of 514 DRMs were filtered from the 1,057 metabolites with a threshold of VIP>1, log 2-fold change>1 and p ≤ 0.05. There were 34 DRMs in the 2012-W02 -vs.- Kaorino comparison group that were detected in three stages simultaneously and the number of DRMs in S1, S2, and S3 were 184, 182, and 187 respectively ( Figure 3A ). The 15 out of 34 DRMs were more accumulated in Kaorino ( Supplementary Dataset S3 ). For the 2012-W02 -vs.- Fenyu NO.1 comparison group, a total of 112, 114, and 134 DRMs were detected in three stages, respectively ( Figure 3A ). The 24 DRMs commonly present in three stages including 9 DRMs up-regulated in Fenyu NO.1 ( Supplementary Dataset S3 ). In addition, a total of 112, 117, and 133 DRMs were identified in the Fenyu NO.1 -vs.- Kaorino comparison group at the three stages, respectively, and 14 DRMs were commonly present at the stages of S1, S2, S3 ( Figure 3A ). Most of these 14 DRMs were more abundant in Fenyu NO.1 than in Kaorino ( Supplementary Dataset S3 ). More DRMs were identified in the 2012-W02-vs.- Kaorino comparison group than in the other two comparison groups ( Figure 3A ). In the three comparison groups, more DRMs were identified in the S3 stage than in the S1 and S2 stages ( Figure 3A ).

The DRMs annotated in both the KEGG and HMDB databases were filtered, resulting in a total of 151 DRMs obtained. Based on the HMDB annotation results, the top ten categories of these 151 DRMs included organooxygen compounds, carboxylic acids and derivatives, fatty acylspyrimidine nucleotides and prenol lipids ( Figure 3B ). The changes in DRM contents were significantly influenced by the different fruit skin colors ( Figure 3C ; Supplementary Dataset S4 ). The DRMs belonging to flavonoids (syringetin, baicalin, isoswertisin 2’’-rhamnoside, chrysoeriol and luteolin 7-O-beta-D-glucoside) showed higher accumulation levels in white skin accession 2012-W02 and pink skin accession Fenyu NO.1 compared to red skin accession Kaorino ( Table 1 ; Figure 2 ). The levels of two DRMs belonging to isoflavonoids, wedelolactone and daidzin, gradually increase with the depening of skin color, with the lowest levels found in the white skin accession 2012-W02 and highest in the red skin accession Kaorino ( Table 1 ).

The KEGG pathway enrichment for DRMs was further conducted to investigate the key pathways involved in skin coloration. The results showed that the DRMs were significantly enriched in glyoxylate and dicarboxylate metabolism, amino acid metabolism, biosynthesis of phenylpropanoids, flavone and flavonol biosynthesis and isoflavonoid biosynthesis ( Supplementary Figure S2 ). Overlap analysis of the DRMs from the nine compared groups revealed that 35, 18 and 22 DRMs overlapped at the S1, S2 and S3 developmental stages, respectively ( Figure 3D ). For the comparison groups of 2012-W02 -vs.- Kaorino and 2012-W02 -vs.- Fenyu NO.1, there were 12 DRMs that occurred commonly at three stages ( Supplementary Dataset S5 ). Among these commonly identified DRMs, the content of syringetin decreased as the skin color deepened, with the highest in the white skin accession 2012-W02 and the lowest in the red skin accession Kaorino ( Table 1 ). In addition, the DRMs of 2,7,4’-trihydroxyisoflavanone were also commonly identified at three developmental stages at the compared group of 2012-W02 -vs.- Fenyu NO.1, and they were downregulated in the Fenyu NO.1 ( Table 1 ).

Transcriptome sequencing of the skins of the three accessions at different develop stages yied a range of 38,297,972 to 53,015,656 clean reads, with Q30 values ranging from 93.13% to 94.47% ( Supplementary Table S2 ). The mapping results showed that the number of mapped reads ranged from 36,360,914 to 50,224,230, accounting for 92.65% to 95.66% of the total number of clean reads ( Supplementary Table S2 ). This indicated that the transcriptome data were of high quality and suitable for further analysis.

A total of 19,918 DEGs were identified among these nine compared groups (2012-W02-vs.- Kaorino, 2012-W02-vs.- Fenyu NO.1 and Fenyu NO.1 -vs.- Kaorino at S1, S2 and S3 development stage, respectively) ( Table 2 ). At three developmental stages, the number of DEGs including up-regulated and down-regulated DEGs was shown in Table 2 . The results showed that there were more DEGs in the 2012-W02-vs.- Kaorino (number of DEG at three development stage: 8,903, 9,641 and 10,578) than in the 2012-W02-vs.- Fenyu NO.1 (number of DEGs at three development stage: 5,222, 5,097 and 6,308) and Fenyu NO.1 -vs.- Kaorino (number of DEGs at three development stage: 4593, 4,859 and 4,415), and the number of DEGs at the S2 and S3 stages was greater than that at the S1 stage in the three compared groups. The 591 DEGs were detected spontaneously in three comparison groups at S1 stage, while the 679 and 869 DEGs were detected at S1 and S2 respectively ( Figure 4A ). In addition, two comparison groups, 2012-W02-vs.- Kaorino and 2012-W02-vs.- Fenyu NO.1, had more co-occurring DEGs at each stage than the other two comparison groups ( Figure 4A ). A total of 1,294 DEGs were detected to co-occur at three stages in the 2012-W02-vs.- Kaorino and 2012-W02-vs.- Fenyu NO.1 comparison groups ( Supplementary Dataset S6 ). These 1,294 DEGs were identified between white skin accession 2012-W02 and red and pink skin accession Kaorino, Fenyu NO.1, and we believe that these genes are closely related to the development of strawberry fruit skin color.

To elucidate the biological functions of these DEGs, GO and KEGG enrichment analyses were conducted. The plant pathways of starch and sucrose metabolism, DNA replication, terpenoid backbone biosynthesis and circadian rhythm were significantly enriched among the DEGs in the W02-vs.- Kaorino comparison group at all three stages ( Figure 4B ). For the comparison groups of the Fenyu NO.1 -vs.- Kaorino, starch and sucrose metabolism, flavonoid biosynthesis and ABC transporters were significantly enriched in DEGs at all three stages ( Figure 4B ). In addition to the starch and sucrose metabolism pathway, the pathways of flavonoid biosynthesis; glyoxylate and dicarboxylate metabolism; valine, leucine and isoleucine biosynthesis; monobactam biosynthesis; as well as glycine, serine and threonine metabolism were consistently enriched in DEGs at three stages in the 2012-W02-vs.- Fenyu NO.1 comparison group ( Figure 4B ). GO enrichment analysis of the DEGs in each comparison group revealed that membrane and membrane part were the most significantly enriched terms in the cell component subcategory. Catalytic activity, binding, and transporter activity were the most enriched terms by DEGs in the molecular function subcategory, while metabolic process was the most enriched term in the biological processes subcategory ( Figure 4C ).

To identify the genes related to the fruit skin color in strawberry plants, the DEGs commonly identified in the Fenyu NO.1-vs.- Kaorino, 2012-W02-vs.- Fenyu NO.1 and 2012-W02-vs.- Kaorino comparison groups at the three fruit development stages were analyzed. A total of 19 DEGs were identified and shown to be involved in pathways related to anthocyanin biosynthesis ( Table 3 ; Supplementary Dataset S7 ). Among them, the expression levels of two anthocyanidin 3-O-glucosyltransferase (BZ1)-encoding genes, Fxa7Dg102859 and Fxa7Ag203300, which are involved in anthocyanin biosynthesis (ko00942), were found to be highest in red skin accession Kaorino at all three developmental stages compared to the other two accessions. Conversely, these genes showed the lowest expression levels in the white skin accession 2012-W02 ( Table 3 ; Supplementary Dataset S7 ). Similar to the expression tendencies of these two genes, the genes enriched in the Flavonoid biosynthesis pathway (ko00941), including the four flavanone 3-hydroxylase (F3H)-encoding genes (Fxa1Bg201013, Fxa1Ag101086, Fxa1Cg101029, Fxa1Dg200969), four chalcone synthase-encoding genes (Fxa7Ag200133, Fxa7Bg200132, Fxa7Cg100092, Fxa7Dg100111), five chalcone-flavanone isomerase-encoding genes (Fxa7Bg202436, Fxa7Dg102232, Fxa7Ag202526, Fxa7Cg102357, Fxa7Dg101819), and two dihy-droflavonol 4-reductase (DFR)-encoding genes (Fxa2Ag103849, Fxa2Bg203601), were also exhibited highest expression levels in Kaorino and lowest in 2012-W02 at all three fruit development stages. In addition, one 4-coumarate-CoA ligase (4CL)-encoding gene (Fxa7Cg103039) involved in the pathway of phenylpropanoid biosynthesis (ko00940) was also most expressed in Kaorino and lowest expressed in 2012-W02 at all three fruit development stages ( Table 3 ; Supplementary Dataset S7 ). In contrast to these genes, the vinorine synthase (shikimate O-hydroxycinnamoyltransferase, HCT)-encoding gene (Fxa4Cg202053), annotated in pathways of phenylpropanoid biosynthesis (ko00940) and flavonoid biosynthesis (ko00941), showed lowest expression level in Kaorino and highest in 2012-W02 ( Table 3 ; Supplementary Dataset S7 ).

The skin color of strawberry fruits is influenced not only by genes that directly catalyzing anthocyanin biosynthesis but also by TFs. A total of 559 differentially expressed TFs were identified in strawberries of different skin colors at different stages, including members of the MYB, bHLH, and WRKY families ( Supplementary Dataset S8 ). To investigate the TFs essential for anthocyanin biosynthesis in strawberry skin color formation, we focus on the four TF families, MYB, bHLH, WD40 and WRKY. A total of 18 TFs were identified and differentially expressed among three accession samples at different developmental stages ( Figure 5 ; Table 3 ). For the MYB family of TFs, the gene Fxa6Cg101113, which encodes an AS1-like protein was up-regulated in white skin accession 2012-W02 and down-regulated in red skin accession Kaorino and pink skin accession Fenyu NO.1 at three stages ( Figure 5 ; Table 3 ). In addition, the R2R3-MYB member MYB10-1 (Fxa1Bg201911) was significantly downregulated in white skin accession 2012-W02 ( Figure 5 ; Table 3 ). It showed lower expression levels at the S1 and S2 stages compared to Fenyu NO.1, and lower expression at the S2 stage when compared to Kaorino ( Table 3 ; Supplementary Dataset S7 ). Seven bHLH family members were found to be differentially expressed between the red skin accession Kaorino and the white skin accession 2012-W02 at all three developmenqtal stages ( Figure 5 ; Table 3 ). Among them, three bHLHs (Fxa2Bg201632, Fxa2Dg202804 and Fxa6Dg101934) were upregulated in Kaorino compared to 2012-W02, while four TFs (Fxa3Bg203691, Fxa3Cg103622, Fxa4Bg103194 and Fxa6Ag102265) were downregulated in Kaorino compared to 2012-W02. The Fxa2Bg201632 gene, which encodes bHLH62, was upregulated in the pink skin accession Fenyu NO.1 at three developmental stages compared with the white skin accession 2012-W02 ( Figure 5 ; Table 3 ). In addition, the Fxa6Ag102265 gene encoding bHLH130 was downregulated in the Kaorino at three developmental stages compared with Fenyu NO.1 ( Table 3 ; Supplementary Dataset S7 ). For the WD40 family of TFs, the eight genes were expressed at lower levels in 2012-W02 than in Kaorino and Fenyu NO.1, three of them (Fxa1Bg202086, Fxa1Cg102222 and Fxa1Dg202079) were downregulated in 2012-W02 at S3 stage as compared to Kaorino and Fenyu NO.1 ( Table 3 ; Supplementary Dataset S7 ). In addition, the gene Fxa1Cg102117, annotated as WRKY13, was found to be differentially expressed in all comparison groups. Its expression decreased with increasing skin color, with the highest level observed in 2012-W02 and the lowest level in Kaorino ( Table 3 ; Supplementary Dataset S7 ).

Numerous studies have shown that strawberry fruit color is mainly controlled by anthocyanin biosynthesis. In this study, pathways related to anthocyanin biosynthesis, including phenylpropanoid biosynthesis (ko00940) (https://www.kegg.jp/pathway/map00940), isoflavonoid biosynthesis (ko00943) (https://www.kegg.jp/pathway/map00943) and flavone and flavonol biosynthesis (ko00944) (https://www.kegg.jp/pathway/map00944) were commonly enriched in DRMs and DEGs.

In the phenylpropanoid biosynthesis pathway (ko00940), four metabolites, p-coumaroyl quinic acid, 5-hydroxyconiferyl alcohol, coniferin and sinapyl alcohol, and 21 genes, PAL (Fxa6Dg101303, Fxa6Cg101365 and Fxa6Ag101540), 4CL (Fxa7Cg103039), HCT (Fxa4Cg202053), CAD (Fxa2Dg203268), CCR (Fxa7Ag200183), F5H (Fxa1Ag101125), REF1 (Fxa1Ag100449), UGT72E (Fxa1Bg201201, Fxa1Dg201152), CYP73A (Fxa3Bg203470, Fxa4Bg100519, Fxa3Ag103799, Fxa3Dg203252, Fxa3Cg103403), CYP98A (Fxa1Dg202656) and COMT (Fxa7Bg203060, Fxa7Dg102804, Fxa7Cg102967 and Fxa7Ag203226), were differentially accumulated or expressed in the three accessions with different skin color. The levels of p-coumaroyl quinic acid and 5-hydroxyconiferyl alcohol were found to be higher in the red skin accession Kaorino compared to 2012-W02 and Fenyu NO.1. Consistently, the expression level of F5H, an upstream gene involved in 5-hydroxyconi coumaroyl feryl alcohol biosynthesis, was higher in the red skin accession Kaorino than in 2012-W02 and Fenyu NO.1. In addition, the metabolites of coniferin and sinapyl alcohol were up-regulated in the white skin accession 2012-W02 compared to red skin accession Kaorino and pink skin accession Fenyu NO.1 ( Figure 6 ; Supplementary Dataset S9 ). Correspondingly, five upstream genes associated with these two metabolites, including four COMT genes and the CAD gene, were also found to be upregulated in 2012-W02 when compared to Kaorino and Fenyu NO.1. The PAL genes showed consistent upregulation in red skin accession Kaorino at all the three developmental stages compared to Fenyu NO.1 and 2012-W02. Additionally, these genes were also found to be upregulated in pink skin accession Fenyu NO.1 at the S2 and S3 stages when compared with 2012-W02. Similarly, gene 4CL showed the highest expression level in red skin accession Kaorino, while the lowest expression level in 2012-W02 across all three developmental stages ( Figure 6 ; Supplementary Dataset S9 ). CCR, CPY73A, and UGT72E genes were also highly expressed in red skin accession Kaorino. The expression of the REF1 gene was significantly higher in Kaorino and Fenyu NO.1 compared to 2012-W02, with the highest level observed in Fenyu NO.1. In contrast, the expression levels of HCT and CYP98A were decreased as the skin color deepened in all three developmental stages, which exhibited a highest in white skin accession 2012-W02.

In the isoflavonoid biosynthesis pathway (ko00943), four DRMs, coumestrol, 2’-Hydroxydaidzein, 2,7,4’-Trihydroxyisoflavanone and daidzin, and six DEGs, I2’H (Fxa2Dg202884), CYP93B2_16 (Fxa5Ag203364), HIDH (Fxa4Ag101504, Fxa2Ag101983), IF7MAT (Fxa7Ag200152), and VR (Fxa7Cg101124), were identified ( Figure 6 ; Supplementary Dataset S9 ). Among them, the content of coumestrol was upregulated in the red skin accession Kaorino compared to Fenyu NO.1 at both S2 and S3 stages and compared to 2012-W02 at S2 stage. Conversely, 2,7,4’-trihydroxyisoflavanone was downregulated in the red skin accession Kaorino compared to 2012-W02 at all three stages. The metabolites 2’-hydroxydaidzein and daidzin were more abundant in the pink-skinned accession Fenyu NO.1 than in Kaorino and 2012-W02. Similarly, the expression level of the CYP93B2_16 was also highest in the pink-skinned accession Fenyu NO.1 across the three accessions ( Figure 6 ; Supplementary Dataset S9 ). Interestingly, the expression levels of all the DEGs in this pathway were found to be up-regulated in the white skin accession 2012-W02 compared to red and pink skin accessions ( Figure 6 ; Supplementary Dataset S9 ).

A total of five DRMs, syringetin, 3,7,4’-Tri-O-methylquercetin, chrysoeriol, acacetin, luteolin 7-O-beta-D-glucoside, and nine DEGs, CYP75B1 (Fxa5Dg201116), FG3 (Fxa5Dg201948), UGT73C6 (Fxa6Bg101474) and IF7MAT (Fxa6Bg104065, Fxa6Bg104066, Fxa6Ag104642, Fxa7Ag200152, Fxa6Cg103966, and Fxa6Ag104436), were identified in the flavone and flavonol biosynthesis pathway (ko00944) ( Figure 6 ; Supplementary Dataset S9 ). All the five metabolites showed the highest accumulation in the white skin accession 2012-W02 and the lowest accumulation in the red skin accession Kaorino). Interestingly, the upstream gene of these DRMs, CYP75B1, which encodes flavonoid 3’-monooxygenase, was also upregulated in 2012-W02 compared to Kaorino and Fenyu NO.1 ( Figure 6 ; Supplementary Dataset S9 ). Similarly, the FG3 gene, which encodes flavonol-3-O-glucoside/galactoside glucosyltransferase, showed the highest expression in white skin accession 2012-W02 compared to the other two accessions ( Figure 6 ; Supplementary Dataset S9 ). The gene UGT73C6, which encodes flavonol-3-O-L-rhamnoside-7-O-glucosyltransferase, showed the highest expression level in pink skin accession Fenyu NO.1 at the S3 stage compared to 2012-W02 and Kaorino. The six IF7MAT genes, encoding isoflavone 7-O-glucoside-6’’-O-malonyltransferase showed higher expression levels in the red skin accession Kaorino than in pink skin accession Fenyu NO.1 and white skin accession 2012-W02.

To validate the RNA-Seq results, eight DEGs associated with anthocyanin biosynthesis were selected for qRT-PCR validation. The transcriptional levels of these genes determined by qRT-PCR analysis showed a similar trend to that observed in the RNA-Seq data ( Figure 7 ), indicating the high reliability of the RNA-Seq data.