Study On The Mechanism of Color Formation In The Skin of Ophiopogon Japonicas Fruit Via Integrative Analysis of Metabolome And Transcriptome


 In this study, the metabolome and transcriptome profiles of ophiopogon japonicas fruit during five stages were analyzed to understand the mechanisms for the formation of color in the skin. Results showed that twenty eight types of flavonoids existed in the skin of the fruit, of which twenty five were anthocyanins and three were flavonols. Among the twenty five kinds of anthocyanins, Delphinidin 3-O-rutinoside, Delphinidin 3-O-glucoside, Petunidin 3-O-rutinoside, Cyanidin 3-O-galactoside and Delphinidin 3-O-galactoside were the major constitutes of the metabolites. Delphinidin 3-O-rutinoside and Delphinidin 3-O-glucoside were discovered to be responsible for the blue color of fruit due to their high correlation with genes in the connection network, especially for Delphinidin 3-O-rutinoside, which was the highest in the amount and accumulated rapidly to 67.6% in the last three stages. Analyzing the results of co-expression network between genes and anthocyanins by means of weighted gene co-expression network analysis (WGCNA), the biosynthesis genes of anthocyanin formed in the late stage were found to be F3’5’H, ANS and UFGT, and the regulatory gene was MYB, which determined the color change of O. japonicas fruit. The identification of key anthocyanins and dominant genes for blue skin of O. japonicas fruit has laid a theoretical basis for regulating the fruit pigmentation in ground cover plants.


Introduction
As a kind of avonoid 1 , anthocyanin consists of anthocyanidin backbone with sugar and acyl conjugates. More than 20 anthocyanidins have been discovered, but only six of them are prevalent in plants naming Cyanidin, Delphinidin, Pelargonidin, Malvidin, Peonidin and Petunidin. The color of anthocyanidin is dominated by the number of hydroxyl groups on the B-ring, and the more the hydroxyl groups, the bluer the color. In general, the blue color of plant tissue is related to the Delphinidin type anthocyanidin with two hydroxyl groups on the B-ring. Delphinidin type anthocyanins include Delphinidin and its derivatives Petunidin and Malvidin, of which the 3' and 5' positions of the B-ring are hydroxylated and/or methoxylated 2 . The accumulation of Delphinidin-type anthocyanins is the necessary but not the only requirement for developing blue/violet tissue. In some species, the interaction or participation of copigments and metal ions are also required 3 .
The biosynthesis pathway of anthocyanin is generally conserved in plants. The biosynthetic pathway for anthocyanin can be divided into two phases: phenylpropanoid and avonoid, both are catalyzed by a series of enzymes. The early avonoid biosynthesis genes (EBGs) are composed of phenylalanineammonia lyase (PAL), 4-coumaryl: CoA ligase (4CL), chalcone synthase (CHS), chalcone isomerase (CHI) and avanone 3-hydroxylase (F3H), which encodes common precursor for various of avonoids. The late anthocyanin biosynthesis genes (LBGs), such as avonoid-3′-hydroxylase (F3′H), avonoid-3′,5′hydroxylase (F3′5′H), Dihydro avonol 4-reductase (DFR), anthocyanidin synthase (ANS) and UDP-glucose: avonoid 3-O-glucosyltransferase (UFGT), on the other hand, are closely relevant to the accumulation of anthocyanins 4 . Additionally, the accumulation of anthocyanins is highly depended on transcription factors (TFs) which regulates the expression of structural genes, and MYB-bHLH-WD40 complex (MBW) is the most common one in TFs. MYB as the most important protein of MBW normally determines the activation or inhibition characteristics of MBW. The MYB gene can independently regulate the structural genes among the anthocyanin biosynthesis pathway [5][6][7] as well. Moreover, the expression of MYB gene in speci c tissues can distinguish target gene promoters, otherwise, unusual expression of MYB gene often leads to abnormal accumulation of anthocyanins in plant 8- 10 .
Integrating analysis of transcriptomics and metabolomics can give insight into transcriptional regulation and metabolic ux. To better know the mechanism of anthocyanin accumulation in fruits, the key genes are always screened and their correlation to the corresponding anthocyanin was disclosed 11  into several co-expression modules which accordingly helps identifying and obtaining hub genes from them 15 .
Ophiopogon japonicas belongs to Ophiopogoneae of Liliaceae family, and is widely used as both medicinal and ornamental plant. It is generally cultivated in the semi-shade of the forest margin or street lawn to cover the bare soil 16,17 . During the fruit season, O. japonicus fruit displays abundant and varying color of fruits from green to blue, which enriches the winter landscape, and accordingly, O. japonicus is also accepted as a superior fruit plant in garden plants. Tang studied the dynamic changes of total avonoids and metal elements in O. japonicus fruits, and found that the amount of Ca 2+ was high 18 . Cheng reported that the blue color of fruit extracts was greatly affected by pH value of the solution, which con rmed that anthocyanins were the main metabolites of O. japonicus fruits 19 . Anthocyanin dominant fruit were also found in other Ophiopogoneae species. Delphinin and Petunidin were obtained in blue fruit of O. jaburan, and the most effective co-pigment kaempferol was determined as well [20][21][22] . By now, four major anthocyanidins including two Delphinidins, one Petunidin and one Malividin have been identi ed in Liriope muscari black fruits 23 . However, few studies have been carried out on the gene regulation mechanism for the fruit color of Ophiopogoneae. It is of great signi cance to study the coloring mechanisms of O. japonicas fruit in order to elevating the ornamental performances of ground cover plants.
In this study, transcriptome and metabolome sequencing were performed on O. japonicas fruit sampling from ve growth stages. 25 anthocyanins and 3 avonols were obtained from specimens of the ve stages. By analyzing the results of transcriptome and RT-qPCR, the expression pro le of related genes for anthocyanin biosynthesis were identi ed. Based on dynamic changes of anthocyanins and genes, the crucial period for color changing was found. The key genes and anthocyanins which responsible for the blue skin of O. japonicas fruits were nally screened out through analyzing the correlation between genes and metabolites. This study not only revealed the mechanisms for anthocyanin accumulation and regulation in the skin of O. japonicas fruit from the aspects of transcriptomics and metabolomics, but also provided important insights in the molecular network of ornamental fruits, which founded molecular basis for the color cultivation of ornamental fruit.

Materials And Methods
Plant materials and sampling O. japonica is widely planted in China, and accordingly, it is not restricted by any relevant legislation or guidelines. Fruits of O. japonicas were collected from a germplasm garden situated at Zhejiang A&F University, Hangzhou, China (N 29 • 56′ to 30 • 23′, E 118 • 51′ to 119 • 52′), and the collection of plant material was allowed by college for experimental research. The fruits were collected evenly from the scape of fteen to twenty plants located at the same area. Five growth stages for the collection were set based on the color and size of fruits. The rst stage (S1) was arranged on September of 2020, followed by the second (S2) on October, the third (S3) on November, the fourth (S4) on December and the fth (S5) on January, 2021. All the freshly picked fruits were frozen immediately in liquid nitrogen after sampling, and Determination on the Pigment content 200 mg fruit from the assigned growth stage was rstly ground into powder in 200 µl absolute ethyl alcohol (EA), then sealed and preserved in dark place. When the fruits turned colorless, the chlorophyll was considered to be dissolved in supernatant liquid. After ltrating the crushed fruits, the lter was collected and ready to be determined. The absorbance of extraction was measured at 470 nm, 649 nm, and 665 nm, respectively. The content of chlorophyll was computed according to the method described by Hu 43 . The total anthocyanin was extracted at 4℃ in dark by a mixed solution which was formulated as EA: H 2 O:HCl = 80:20:0.3 (calculated as volume, mL) 32 . The extracts were measured at the absorbance 530 nm and 600 nm. The relative content (amount) of total anthocyanin (TAC) was computed according to Liu's method 44 .

Sample of metabonomics preparation and disposition
The freeze-dried sample was crushed with zirconia beads using a mixer mill (MM 400, Retsch) at 30 Hz for 1.5 min. 50 mg powder was vortexed in 0.5 ml methanol/water/hydrochloric acid (799:200:1,V/V/V) for 10 min and ultrasonically extracted for 10 min. The extracts were then centrifuged at 12, 000 rpm under 4°C for 3 min to allow phase separation. The extraction was repeated twice and the total supernatants were ltrated and collected (PTFE, 0.22 µm; Anpel) for LC-MS/MS analysis.

ESI-Q TRAP-MS/MS
Linear ion trap and triple quadrupole (QQQ) scans were acquired by a triple quadrupole-linear ion trap mass spectrometer (Q TRAP). API 6500 Q TRAP UPLC/MS/MS System, equipped with an ESI Turbo Ion-Spray interface, was operated in positive and negative ion mode and controlled by Analyst 1.6.3 software (AB Sciex). The ESI source operation parameters were as follows: ion source with turbo spray, source temperature at 550℃; ion spray voltage (IS) at 5500 V (positive ion mode), curtain gas (CUR) at 35 psi, and the highly collision active dissociation (CAD). Each ion pair was scanned according to the optimized declustering potential (DP) and collision energy (CE). The metabolites eluted in each period were quanti ed using Multiple Reaction Monitoring (MRM).
Total RNA extraction and quality analysis Total RNA was extracted by RNA extraction kit (TIANGEN RNAprep FFPE Kit) according to the corresponding speci cations. The integrity of RNA was detected by Agilent 2100 Bioanalyzer (Agilent Technologies Inc., CA, USA) and agarose gel electrophoresis, while the purity and concentration of the RNA were analyzed by a NanoDrop™ One/OneC system (Thermo Fisher Scienti c, MA, USA).

Iso-Seq library construction and Data Analysis
Total RNA of the ve stages was maximized to acquire the full-length transcriptome. Iso-Seq library was prepared according to the Isoform Sequencing protocol (Iso-Seq) and using the Clontech SMARTer PCR cDNA Synthesis Kit protocol as described by Paci c Biosciences(www.pacb.com, PN 100-092-800-03). After PacBio system sequencing, 33.75 Gb subreads were obtained. SMRTlink 7.0 software was used to process these sequence data as follows. The circular consensus sequences (CCS) were generated from subreads of BAM les, and the resulted CCS contained full length non-chimeric (FLNC) sequences. The nal isoforms were gained from high-quality consensus FLNC reads after removing the redundant FLNC through Quiver and CD-HIT-v4.6.754, respectively.
Construction of Illumina RNA-seq library and estimation on expression level RNA sequences of fruits from ve growth stages, with three replicates for each stage were concluded, which means a total number of fteen libraries were established. Each sequencing library was prepared with 1 µg RNA and was generated from NEBNext® UltraTM RNA Library Prep Kit of Illumina® (NEB, USA).
The quantitation and quality of library was detected by Qubit2.0 Fluorometer and Agilent 2100 bioanalyzer, respectively. When passing the quality inspection, the library preparations were sequenced on an Illumina Novaseq platform and 150 bp paired-end reads were generated.
The expression level of each gene was determined by the fragments-per-kilobase of transcript-per-millionmapped reads (FPKM), which was calculated by featureCounts v1.5.0-p3. Differential expression analysis between two stages was performed using the DESeq2 R package (1.16.1). The obtained P-values were then adjusted according to the Benjamini and Hochberg's approach for controlling the false discovery rate. Genes with an adjusted P-value < 0.05 detected with DESeq2 were designated as differentially expressed genes (DEGs) (fold-change ≥ 2, false discovery rate < 0.01).

qRT-PCR analysis
The extraction of RNA was carried out according to the method mentioned above and PrimerScript™ RT Master Mix cDNA (TaKaRa, Japan) was used for cDNA synthesis. The primer pairs for each isoform/unigene were designed using primer 5. The obtained cDNA samples were analyzed by qRT-PCR in a 20 µl reacting solution formulated with TB Green™ Premix Ex Taq™ (TaKaRa, Japan), primer pairs and double distilled water. The qRT-PCR analysis was performed on a Light-Cycler 480II real-time PCR detection system (Roche, USA).
Statistical analysis SIMCA 14 was applied to conduct the Principal Component Analysis (PCA) and Variable Importance in Projection (VIP). One-way ANOVA was prformed by means of Duncan's multiple range test of the Statistical Product and Service Solutions program (version 19) (SPSS Inc.,Chicago, IL), and the graphs were plotted using Graphpad prism 8.0.

Relationship between the color development and the content variation of chlorophyll and anthocyanin
The growth of O. japonicas fruits was divided into ve stages according to the color of skin and size of fruit. Figure. 1a-c presented the growth of fruit, dimension of long axis (length) and the hundred-grain weight (HGW). During S1 and S2, green was the dominant color of the fruit and it became darker from S1 to S2. S3 was the transition period for color development and when entering the fourth stage, blue turned to be the prevailing color. From S4 to S5, the blue turned deeper. The average length of fruit on S1was as low as 5.83 ± 0.20 mm, but increased continuously by 1.11-fold, 1.48-fold and 1.09-fold from S1 to S2, S2 to S3 and S3 to S4, respectively, then leveled off from S4 to S5. The variation of HGW had similar tendency as the length. It grew from 3.38 ± 0.23 g on the rst stage to 11.07 ± 0.60 g on the fth stage, the weight gain of which increased signi cantly from stage 1 to stage 4 (P < 0.05) and slowed down from stage 4 to stage 5.
The content of chlorophyll and anthocyanin accumulated in the skin of O. japonicas fruits was displayed in Fig. 1d. The amount of chlorophyll climbed rapidly from 14.38 ± 1.62µg/g FW to 25.50 ± 2.84 µg/g FW during the rst two stages and declined quickly to the lowest from S2 to S4. On the contrary, the content of total anthocyanin was near zero during the rst two stages, but soared in the following three stages, reaching the highest in S4, followed by a slight drop on S5. Interestingly, the relative content of chlorophyll was equivalent to that of anthocyanin around stage three. When correlating the proportion of chlorophyll and anthocyanin to the color of the fruit, it is understandable that the fruit developed from green dominant color to blue prevailing color, with a transition period appeared around stage three. Hereby, the phenotypic changes of O. japonicas fruits skin were highly consistent with the content of anthocyanin and chlorophyll.

The Metabolites pro les of O. japonicas Fruits skin
The fruits of O. japonicus were collected from the ve growth stages and the composition of the skin was analyzed by means of metabolome. 28 types of avonoids were obtained, of which 25 was classi ed as anthocyanin compounds and 3 was avonols. Among the 25 anthocyanin compounds, 6 was the main and most common anthocyanins, which proved that O. japonicas fruits contained complete biosynthesis pathway for the anthocyanin. Further study on the accumulation of the 6 types of anthocyanin compounds and avonol compounds during every growth stage was conducted and the results were presented in Fig. 1-e. It is evident that Delphinidin derivatives were the major metabolites, which accounted near 80% in the last three stages. Principal component analysis (PCA) was adopted to investigate the metabolites of the ve growth stages and the results were shown in Fig. 2. Samples from S1 and S2 were obviously clustered, which con rmed that there was no signi cant change in the type and amount of metabolites during the rst two stages. But the following three stages could be de ned clearly through PC1 explaining 66.4% and PC2 explaining 12.6%. Based on the analyses above, S3 was considered to be the key period for color variation of fruits.

Differentially Accumulated Metabolites (DAMs) Analysis of O. japonicas Fruit skin
To better understand the major metabolites affecting the coloration of O. japonicas Fruit, differentially accumulated metabolites (DAMs) analysis was applied under the condition of VIP ≥ 1, p < 0.05, Log2FC ≥ 1 or Log2FC≤-1. A total number of 21 avonoids were found which included 19 anthocyanins and 2 avonols. The heat maps (Fig. 3. a) showed that most anthocyanins increased sharply in the amount during S3 and leveled off during S5, of which Delphinidin 3-O-rutinoside, Delphinidin 3-O-glucoside, Petunidin 3-O-rutinoside, Cyanidin 3-O-galactoside and Delphinidin 3-O-galactoside varied dramatically and became the major metabolites at the ripening stage, accordingly, they were considered to be the key metabolites for the transition of fruit color from green to blue. In addition, among the ve anthocyanins, Delphinidin 3-O-rutinoside was the most abundant one, which accounted for 69.2% of the DAMs in S4, 95.57 times higher than that at S1.
The expression pro les of avonols, on the other hand, were similar to anthocyanins in the rst four stages, but decreased signi cantly at S5, which could be the key metabolites to distinguish the last two stages. The content of Quercetin 3-O-glucoside and Dihydromyricetin declined rapidly at S5, 7.05-fold and 2.39-fold lower, respectively, in comparison with S4.

Gene expression patterns of avonoid biosynthesis
Based on KEGG pathway, 293 genes regulating anthocyanin and avonol biosynthesis were identi ed, and 181 of them was differentially expressed ( Table 1). As shown in Fig. 3. b, a part of EBGs became low in expression during the fruit ripening stages, which mainly involved in eight PALs, two C4Hs, two CHSs, ve CHIs and one F3H. In comparison, the LBGs including DFRs, ANSs, UFGTs and OMTs, shared the same expression pro le with anthocyanins except one DFR and ve OMTs. Owing to the high expression, LBGs were considered to be positively correlated with the accumulation of anthocyanin. In addition to the biosynthesis genes, the transcription factors (TFs) also affected the accumulation of anthocyanin. Therefore, a total of 817 TFs regulatory genes assigned to 67 gene families were identi ed, of which 351 belonged to DEGs. The major TFs of MYBs and bHLHs, which control the biosynthesis of anthocyanins, were focused on and their expression levels were analyzed in detail. The expression level of 7 MYBs was found to be consistent with the accumulation of anthocyanins, and three of them including maidong_transcript_11129, maidong_transcript_7157 and maidong_transcript_6578 were highly expressed in S3. Nevertheless, bHLH gene (maidong_transcript_6829), the only differentially expressed, after reaching the highest expression level in S2 decayed gradually in the subsequent periods (Fig. 3. c). : Shared biosynthesis genes of anthocyanin and avonol. : Speci c biosynthesis genes of anthocyanin. : Speci c biosynthesis genes of avonol. : Two of common TFs that regulate avonoid synthesis.
Obviously, the FLSs and FG3s genes, which are vital for avonol biosynthesis expressed strongly at S3 and S4, and correspondingly the content of Quercetin 3-O-glucoside and Dihydromyricetin was high in the fruit. The expression level at S5 was signi cantly lower than in S4, which were consistent with the variation of avonols.
Identifying anthocyanin-related DEGs by Co-expression network analysis 4004 DEGs (with mean FPKM > 45) were selected in this study to perform WGCNA. According to the standard WGCNA networks, the soft power was set as 13, min Module Size as 30 and the merge Cut Height as 0.15. The DEGs were clustered into 11 tree branches, and each one represented a module displaying different color (Fig. 4., Fig. 5.). The highly related DEGs in the same module were co-expressed according to the standard WGCNA, but the color of Cyanidin 3-O-glucoside corresponding to the modules were different from those of the other DAMs, which suggested that the gene regulatory network of Cyanidin 3-O-glucoside was independent. Further analysis revealed that six modules were highly positively correlated to DAMs (P < 0.05), while three modules were highly negatively correlated to DAMs (P < 0.05, Table 2). Based on the analyses above, the centrality weights of modules nodes were obtained by means of CentiScaPe (http://apps.cytoscape.org/apps/centiscape#cy-app-releases-tab), that was, by calculating the geometric means of degree centrality, closeness centrality and betweenness centrality with the edge weights between 0.2 and 0.8, the Hub genes were then screened from the highest centrality weight. Among the positive modules, hub genes that were predicted to encode the synthesis enzymes of anthocyanins were more abundant in purple module (Fig. 6.), including ve ANSs (maidong_transcript_18463, maidong_transcript_14175, maidong_transcript_12015, maidong_transcript_20037 and maidong_transcript_22160), F3'5'H (maidong_transcript_2898) and UFGT (maidong_transcript_9469). The maidong_transcript_6056 (unknown protein) identi ed as hub gene in purple module was highly correlated to the anthocyanin synthesis genes including CHS, three ANSs, UFGT and F3'5'H, which suggested that maidong_transcript_6056 gene played an important role in the synthesis of anthocyanin. Except for purple module, four hub genes including C4H (maidong_transcript_3366), two OMTs (maidong_transcript_18253 and maidong_transcript_9185)and GST maidong_transcript_22603 sourced from pink module, greenyellow module and midnightblue module in group , were also related to the synthetic pathways of anthocyanin (Table 2). In addition to the genes described above, two iso avone 2'-hydroxylase (I2'H) genes (maidong_transcript_3111 and maidong_transcript_3164) ascribed to midnightblue module were found to be homologous to soybean gmI2'H 24 , but their expression increased continuously during the ve stages, indicating that the synthesis of iso avone was slightly different from that of anthocyanin. Thirty hub genes were screened from the three modules of group , in which the DEGs negatively responded to DAMs (Table 2), and 9 of hub genes were annotated in KEGG. Among the nine genes, LHCB1 (light-harvesting Complex II chlorophyll A/B binding protein 1), which belonged to the red module and encoded light-harvesting Complex protein, participated in the photosynthesis. The expression of LHCB1 that directly affected the content of chlorophyll 25 declined continuously during the ve stages, and the expression level dropped 35.01 times at S3 in comparison with S1. Maidong_transcript_13873, as the hub gene of cyan module, was annotated as caltractin protein, one of the Ca-binding protein, which displayed a continuous drop in gene expression from S2 to S5. However, hub genes of turquoise module which negatively respond to DAMs were not annotated in KEGG.
Relationship between gene expression and the amount of anthocyanin in O. japonicas fruit skin In order to investigate the gene regulatory network of anthocyanins biosynthesis in O. japonicas fruits, ve anthocyanins and 181 DEGs (Table 1) pro les were selected and the Pearson's correlation coe cients were calculated. Correlations with a coe cient of R 2 > 0.8 were adopted to visualize a network, which was displayed in Fig. 7. There were 5 anthocyanins and 89 DEGs in the diagram, including 75 structural genes and 14 regulatory genes. The highest centralities assigned to the module's nodes were obtained according to the method described in co-network. DEGs and DAMs included in the six nodes were two UFGTs (maidong_transcript_3906 and maidong_transcript_13974), two regulatory genes MYBs (maidong_transcript_7157 and maidong_transcript_6578) and two anthocyanins (Delphinidin 3-Orutinoside and Delphinidin 3-O-glucoside), respectively. Among the six highly correlated DEGs and DAMs, the correlation between two MYBs and UFGTs (maidong_transcript_3906) reached 0.8841 and 0.9090, respectively. In addition, the expression levels of the two MYBs were closely related to the accumulation of Delphinidin-types anthocyanin. When analyzing their pro le in the fruits, the variation of these candidate nodes was basically similar, and their expression level or content were signi cantly increased in S3, which suggested that these DEGs and DAMs were responsible for color formation of O. japonicas fruit.

Validation via qRT-PCR
Based on the results of WGCNA and the correlation between DAMs and DEGs, qRT-PCR was conducted on 11 hub genes screened from purple module (Fig. 6.) and co-network (Fig. 7.). Taking RNA-binding protein Musashi (MSI) and ubiquitin-fold modi er conjugating enzyme 1 (UFC1) as the reference genes, the 11 hub genes for the O. japonicas fruits at ve stages were veri ed and presented in Fig. 8. a-k. To validate the accuracy of RNA-seq analysis, the coe cient of Pearson's correlation was computed. Except for one MYB gene (maidong_transcript_6578), the expression patterns of RNA-seq and RT-qPCR were highly consistent, especially for ANSs, and the R 2 ranged from 0.7543 to 0.9824 with p < 0.05.
The expression level of 11 hub genes was low or scarce during the rst two periods. When entering S3, the expression level climbed gradually, and the accumulation of anthocyanin promoted the transition of O. japonicas fruit from green to blue. This gene expression reached the peak at S4, and the content of anthocyanin also attained the maximum. Subsequently in S5, the expression level of genes dropped rapidly except for one MYB (maidong_transcript_6578), and the accumulation of anthocyanin ceased. Based on the analyses above, MYB, ANS, UFGT and F3'5'H were considered to be the key types of genes for anthocyanin synthesis.

Discussion
The color of plant organs is generally controlled by the type and amount of anthocyanins. Delphinin was reported to be the most common anthocyanin in blue organs of plant 26 . Studies on the mature fruit of O. japonicas also revealed that Delphinin pigments were dominant, followed by Petunidins, which agreed with the blue color of its fruit and con rmed that large proportion of Delphinin was the requirement for blue organ 27 . When the content of Petunidin and Malvidin increases, the carboxylation and/or methoxylation reactions happened on the 3 ' and 5 ' positions of the B-ring of Delphinin reduced the reactivity of hydroxyl groups 28 , consequently, the color shifts red slightly 29 . Alternatively, when the content of Cyanidin is equivalent to that of Delphinin, the petal or fruit often shows purple phenotypes 30,31 . The synthesis of Anthocyanin is regulated by a variety of structural genes, and it is di cult to identify key gene only through analyzing the expression pro le of DEGs. Therefore, WGCNA, as a mean of dividing coexpression DEGs into many modules, was introduced in this experiment. Hub genes were identi ed as the highest correlation in each module and could represent the expression features of the module 35 . In positive modules that correlated with anthocyanin (P < 0.05), the most LBGs that identi ed as hub genes enriched in purple module, including F3'5'H, three ANSs, and UFGT. RT-qPCR validation con rmed that dynamic changing between relative quantitation and FKPM of these genes were similar. Based on the analyses above, F3'5'H, ANS and UFGT were considered as the key genes for the accumulation of Delphinidin, which was proved in other blue/violet plants. F3'5'H was located at the rst node of anthocyanin biosynthesis pathway, and when it competed with F3'H for dihydrokaempferol, a metabolic intermediate, Delphinidin and Cyanidin will be generated. The F3'5'H gene is so important that plant will fail to produce blue/purple owers without it 36 . The expression of F3'5'H was enhanced from S3 laid foundation for the accumulation of Delphinidin in the skin of O. japonicas fruits. However, only high expression of F3'5'H was not enough. ANS and UFGT as LBGs were also required for the synthesis of anthocyanin, and the expression level of which determined the color of organs. A total of 181 DEGs related to anthocyanin biosynthesis were identi ed by KEGG pathway, the LBGs of which expressed increasingly during S3 except one DFR and ve OMTs. Combined with the analysis of WGCNA, ANS, UFGT and F3'5'H were con rmed as the hub genes, the expression levels of which were highly consistent with the variation of Delphinin content. It was con rmed that LBGs regulated the accumulation of anthocyanin in O. japonicas fruits, which was similar with previous studies. For example, Chen reported that the color of eggplant petal was controlled by a single dominant gene called Flower Anthocyanidin Synthase (FAS), which was encoded into the ANS protein. A single base-pair deletion of FAS at the 438th position will lead to transcription halt and the failure of anthocyanin accumulation, as a result, the plant produced mutant with white-ower 37 . Zhang con rmed that ANS was a key gene in Hosta plantaginea purple petals 38 . Under the catalyzation of ANS, anthocyanidins will produce exposed hydroxyl groups, which could be modi ed via glycosylation and the formation of intramolecular bonds made them stable under acidic conditions. UFGT, as a modi er gene, also signi cantly affects the accumulation of blue/violet anthocyanins. It is reported that the mechanism for exogenous calcium signal to induce the accumulation of the anthocyanin in Vitis vinifera was the interaction between calmodulin and UFGT 39 . Noda succeeded in cultivating blue chrysanthemum by introducing UGFT, whereas non-glycogenated transgenic chrysanthemums showed purple 40 . These results suggested that the high expression of F3'5'H, ANS and UFGT genes was essential for the accumulation of Delphinin.
In addition to LBGs, TFs also regulate the accumulation of Anthocyanin, of which MYB is the most common and critical one. In this research, correlation analyses disclosed that MYB gene located at the center of the regulatory network (Fig. 7.) and was co-expressed with the key gene UFGT (maidong_transcript_3906). Moreover, the expression pro le of MYB gene was closely related to the content of Delphinidin 3-O-rutinoside and Delphinidin 3-O-glucoside. These results were consistent with the previous studies. Masumi reported that MYB19LONG in Lilium brownii was the essential gene involving in producing anthocyanin pigmentation, which was stimulated by different vascular bundle signals to develop anthocyanin pigmentation related brushmarks or raised spots 41 . Experiments of SSMYB1 conducted by overexpression in Arabidopsis thaliana and VIGS silence in Sapium sebiferum leaves con rmed that the accumulation of anthocyanins in S. sebiferum leaves was depended on SSMYB1 42 . However, the variation of MYB gene expression always causes abnormal anthocyanin accumulation in fruit (e.g., Citrus grandis 8 , 'Fuji' apple (Ralls Janet × Red Delicious) 9 and Elaeis guineensis 10 ). Therefore, it was speculated that MYB protein activated the expression of LBGs, leading to the accumulation of Delphinidin in O. japonicas fruits.

Conclusions
In summary, the blue skin of O. japonicas fruits was originated from Delphinidin, and of which, Delphinidin 3-O-glucoside and Delphinidin 3-O-rutinoside were the key metabolic substances for the color transformation. In addition, the transcription factor MYB gene and downstream genes LBGs including F3'5'H, ANS and UFGT were the key genes for anthocyanin biosynthesis. These results will not only serve the purpose for elevating fruit pigmentation of O. japonicas and the other proximal species, but also provide important ways to the identi cation of related genes and proteins involved in blue fruit development.      Relationships between module eigengenes (ME, rows) and traits as DAMs (columns).

Figure 6
Correlation networks in the purple module.

Figure 7
The correlation networks of 5 major anthocyanins and 89 DEGs, including 75 structural genes and 14 regulatory genes.