Transcriptomics and Metabolomics Reveal the Critical Genes of Carotenoid Biosynthesis and Color Formation of Goji (Lycium barbarum L.) Fruit Ripening

Carotenoids in goji (Lycium barbarum L.) have excellent health benefits, but the underlying mechanism of carotenoid synthesis and color formation in goji fruit ripening is still unclear. The present study uses transcriptomics and metabolomics to investigate carotenoid biosynthesis and color formation differences in N1 (red fruit) and N1Y (yellow fruit) at three stages of ripening. Twenty-seven carotenoids were identified in N1 and N1Y fruits during the M1, M2, and M3 periods, with the M2 and M3 periods being critical for the difference in carotenoid and color between N1 and N1Y fruit. Weighted gene co-expression network analysis (WGCNA), gene trend analysis, and correlation analysis suggest that PSY1 and ZDS16 may be important players in the synthesis of carotenoids during goji fruit ripening. Meanwhile, 63 transcription factors (TFs) were identified related to goji fruit carotenoid biosynthesis. Among them, four TFs (CMB1-1, WRKY22-1, WRKY22-3, and RAP2-13-like) may have potential regulatory relationships with PSY1 and ZDS16. This work sheds light on the molecular network of carotenoid synthesis and explains the differences in carotenoid accumulation in different colored goji fruits.


Introduction
Carotenoids, synthesized in the plastids of photosynthetic and sink organs in plants, are essential molecules for photosynthesis and are responsible for natural colors, such as yellow, orange, and red [1,2]. Carotenoids are also vital for human health and nutrition because of their role as vitamin A precursors and antioxidants [3].

Carotenoid Differences in N1 and N1Y during Fruit Ripening
A previous study identified a bud mutant cultivar N1Y of N1, which has similar growth characteristics and fruit size to N1. However, the fruit color of N1Y is yellow, while N1 is red (Figure 1a). This could be due to a difference in carotenoid levels between N1 and N1Y during fruit ripening. We analyzed fresh fruit samples from the two cultivars at different stages of fruit ripening, namely M1, M2, and M3 (Figure 1a), to compare the carotenoid profile in N1 and N1Y. Total carotenoid content in N1 and N1Y continued to increase during fruit ripening, and N1Y fruit had significantly higher carotenoid levels than N1 during the M2 and M3 periods (Figure 1b). Based on the principal component analysis (PCA) of carotenoids, sampling of N1 and N1Y at different stages was readily repeatable and could be used for subsequent analysis (Figure 1c).
Our study identified 27 carotenoids in the fruits of N1 and N1Y (Figures 2 and S1 and Table S1). Phytoene, zeaxanthin, zeaxanthin dipalmitate, zeaxanthin palmitate, beta-cryptoflavin palmitate, and rubixanthin palmitate were among those whose levels rose during fruit ripening. In contrast, beta-carotene and xanthophyll levels decrease dramatically during fruit ripening ( Figure 2). During the M1 period, xanthophyll was the major carotenoid in N1 and N1Y fruits, accounting for 57.2% and 59.7% of the total carotenoid content, respectively. According to the PCA analysis for the M1 period, PC1 (95.27%) has a well explained variance, and xanthophyll was the main loading component Small letters indicate that N1 or N1Y have significant differences between M1, M2, and M3 (p < 0.05), as analyzed by Duncan's multiple tests; * means a significant difference between N1 and N1Y at M1, M2, and M3 periods (p < 0.05), as analyzed by Duncan's multiple tests; NS means no significant difference. (c) Principal component analysis for carotenoid content of N1 and N1Y fruits at the M1, M2, and M3 periods.
Our study identified 27 carotenoids in the fruits of N1 and N1Y (Figures 2 and S1 and  Table S1). Phytoene, zeaxanthin, zeaxanthin dipalmitate, zeaxanthin palmitate, beta-cryptoflavin palmitate, and rubixanthin palmitate were among those whose levels rose during fruit ripening. In contrast, beta-carotene and xanthophyll levels decrease dramatically during fruit ripening ( Figure 2). During the M1 period, xanthophyll was the major carotenoid in N1 and N1Y fruits, accounting for 57.2% and 59.7% of the total carotenoid con-  Thus, the M2 and M3 periods were critical for the carotenoid content difference between N1 and N1Y fruit. Furthermore, differences in the main carotenoids can explain the color variation between N1 and N1Y fruit at the M2 and M3 periods. Much of the The change in individual carotenoid content in N1 and N1Y fruits at different ripening stages. Small letters indicate that N1 or N1Y have significant differences between M1, M2, and M3 (p < 0.05), as analyzed by Duncan's multiple tests; * means a significant difference between N1 and N1Y at M1, M2, and M3 periods (p < 0.05), as analyzed by Duncan's multiple tests; NS means no significant difference. The green, orange, and red asterisks represent the primary carotenoid of N1 and N1Y fruits at the M1, M2, and M3 periods, respectively. Thus, the M2 and M3 periods were critical for the carotenoid content difference between N1 and N1Y fruit. Furthermore, differences in the main carotenoids can explain the color variation between N1 and N1Y fruit at the M2 and M3 periods. Much of the zeaxanthin in the N1 fruit was converted to zeaxanthin palmitate and zeaxanthin dipalmitate during the M2 and M3 periods, while the N1Y fruit continued to accumulate zeaxanthin. As a result, from M2 to M3, the N1 fruit gradually turned red, whereas the N1Y fruit turned yellow.

Transcriptome Profiles of N1 and N1Y Fruits in Different Ripening Periods
The gene expression profiles of N1 and N1Y fruits at the M1, M2, and M3 periods were examined using RNA sequencing to identify the essential genes for carotenoid biosynthesis at various ripening stages. We obtained an average of 7 G clean bases per sample, with 93.38% of bases scoring Q30 (Table S3). Over 92% of the reads mapped to the reference genome of goji (Table S4) Table S5). N1Y-M1 vs. N1Y-M2, N1Y-M1 vs. N1Y-M3, and N1Y-M2 vs. N1Y-M3 shared 943 common DEGs, which were primarily involved in alpha-linolenic acid metabolism, photosynthesis, carotenoid biosynthesis, and starch and sucrose metabolism (Figure 4b). The difference between N1 and N1Y at the M2 and M3 periods was significantly responsible for the change in carotenoid content and the resulting fruit color. Therefore, the DEGs of N1-M1 vs. N1-M2, N1-M1 vs. N1-M3, N1Y-M1 vs. N1Y-M2, and N1Y-M1 vs. N1Y-M3 were analyzed. Finally, 5614 DEGs were found in signaling and cellular processes, carotenoid biosynthesis, and carbohydrate, starch, sucrose, and amino acid metabolism ( Figure 4c). We selected 13 genes (Table S13) to validate our RNA-seq data. qRT-PCR was used to examine the expression levels of these genes in N1 and N1Y during M1, M2, and M3 ( Figure S2). The results revealed that the RNA data were valid and trustworthy, as the expression patterns of these genes matched the RNA-seq results.

Expression Profiles of Genes of the Carotenoid Biosynthetic Pathway during Goji Fruit Ripening
To further investigate the critical genes of carotenoid biosynthesis during goji fruit ripening, we analyzed the expression of the carotenoid biosynthesis genes in N1 and N1Y fruits ( Figure 5). Three phytoene synthase genes (PSY1, PSY2-1, and PSY2-2) involved in

Identification of WGCNA Modules of Carotenoid Biosynthesis in Goji Fruit
A weighted gene co-expression network analysis (WGCNA) was performed for all N1 and N1Y fruit gene expressions to identify critical genes involved in carotenoid biosynthesis, and 17 WGCNA modules were identified (Figure 7a, Table S8). We discovered that the main carotenoids of N1 and N1Y fruits at the M1, M2, and M3 stages were zeaxanthin, phytoene, xanthophyll, zeaxanthin dipalmitate, zeaxanthin palmitate, beta-cryptoflavin palmitate, and rubixanthin palmitate. Therefore, we primarily sought modules highly related to these substances. According to the WGCNA, the yellow (0.94, p < 0.05), brown (0.70, p < 0.05), and tan modules (0.65, p < 0.05) had a positive correlation with total carotenoid content, while the pink (−0.9, p < 0.05) and light cyan modules (−0.64, p < 0.05) had a negative correlation. The cyan PSY1 and ZDS16 were found to be essential genes for goji fruit carotenoid biosynthesis during fruit ripening. The brown module consisted of PSY1 and ZDS16 (Table S8). Moreover, total carotenoid, zeaxanthin, xanthophyll, beta-carotene, zeaxanthin dipalmitate, neoxanthin, violaxanthin, and lutein palmitate was also strongly correlated with the brown module ( Figure 7a). Therefore, the genes of the brown module were worthy of further analysis and study. Meanwhile, we also used the Mfuzz software to examine the expression trend of genes in goji fruit at various stages of ripening. The findings revealed that differential genes could be divided into six clusters (Figure 7b, Table S9). Among them, PSY1 and ZDS16 were classified in cluster 2.

Identification of Key Transcription Factors Associated with Carotenoid Biosynthesis Pathway
Many studies have shown that transcription factors (TFs) play vital regulatory functions in the biosynthesis of carotenoids in plants [1,4]. The genes of the brown module were worthy of further analysis and study. As a result, we conclude that some TFs in the brown module may be involved in the regulation of carotenoid synthesis. We identified 63 TFs from the brown module (Table S10) that can be classified into two groups based on their expression levels ( Figure 8). Since PSY1 and ZDS16 were more important for carotenoid biosynthesis in goji fruit during ripening, we propose that TFs meet the following three criteria to be considered as potential regulators of carotenoid synthesis. First, the correlation with PSY1 and ZDS16 must be more than 0.9 or less than −0.9 (Table S11). There should also be more than a 2-fold increase or decrease in the expression in M2 or M3 compared to M1 (Figure 8). Lastly, the network weight values with PSY1 and ZDS16 must be greater than 0.5 (Table S12). We finally selected three TFs (CMB1-1, WRKY22-1, and WRKY22-3) that positively regulate the carotenoid synthesis and one TF (RAP2-13-like) that negatively regulates the carotenoid synthesis ( Figure 9). There should also be more than a 2-fold increase or decrease in the expression in M2 or M3 compared to M1 (Figure 8). Lastly, the network weight values with PSY1 and ZDS16 must be greater than 0.5 (Table S12). We finally selected three TFs (CMB1-1, WRKY22-1, and WRKY22-3) that positively regulate the carotenoid synthesis and one TF (RAP2-13like) that negatively regulates the carotenoid synthesis ( Figure 9).  ER REVIEW 15 of 21 Figure 9. Patterns of transcription factor regulation of carotenoid biosynthesis in goji fruit. The black background indicates that the substance can be detected in N1 and N1Y fruits. The font color represents the substance color. Small letters indicate that N1 or N1Y have significant differences between M1, M2, and M3 (p < 0.05), as analyzed by Duncan's multiple tests; * means a significant difference between N1 and N1Y at M1, M2, and M3 periods (p < 0.05), as analyzed by Duncan's multiple tests; NS means no significant difference.

Discussion
As a traditional Chinese medicinal and nutritional plant, goji (Lycium barbarum L.) has attracted much interest from scientists worldwide, especially regarding the extraction of goji carotenoids and their role in human health [18]. Previous studies discovered 18 carotenoids in goji fruit extracts [7]. Using liquid chromatography-tandem mass spectrometry (LC-MS/MS), the current study identified 27 carotenoids in N1 and N1Y fruits during fruit ripening (Figures 2 and S1). Peng, Ma, Li, Leung, Jiang and Zhao [9] indicated that the carotenoid content in goji berries increases as the fruit ripens, and zeaxanthin dipalmitate, which makes up 31-56% of all the carotenoids in fruits, is the primary carotenoid. In the present study, we further analyzed the carotenoids content in N1 and N1Y in the three stages of fruit ripening from M1 to M3 (Figure 1b). In N1 and N1Y fruits in Figure 9. Patterns of transcription factor regulation of carotenoid biosynthesis in goji fruit. The black background indicates that the substance can be detected in N1 and N1Y fruits. The font color represents the substance color. Small letters indicate that N1 or N1Y have significant differences between M1, M2, and M3 (p < 0.05), as analyzed by Duncan's multiple tests; * means a significant difference between N1 and N1Y at M1, M2, and M3 periods (p < 0.05), as analyzed by Duncan's multiple tests; NS means no significant difference.

Discussion
As a traditional Chinese medicinal and nutritional plant, goji (Lycium barbarum L.) has attracted much interest from scientists worldwide, especially regarding the extraction of goji carotenoids and their role in human health [18]. Previous studies discovered 18 carotenoids in goji fruit extracts [7]. Using liquid chromatography-tandem mass spectrometry (LC-MS/MS), the current study identified 27 carotenoids in N1 and N1Y fruits during fruit ripening (Figures 2 and S1). Peng, Ma, Li, Leung, Jiang and Zhao [9] indicated that the carotenoid content in goji berries increases as the fruit ripens, and zeaxanthin dipalmitate, which makes up 31-56% of all the carotenoids in fruits, is the primary carotenoid. In the present study, we further analyzed the carotenoids content in N1 and N1Y in the three stages of fruit ripening from M1 to M3 (Figure 1b). In N1 and N1Y fruits in the M1 period, xanthophyll was the main carotenoid, accounting for 57.2% and 59.7% of the total carotenoid content, respectively (Figure 2). Zeaxanthin and zeaxanthin dipalmitate, which together accounted for 56.2% and 77.3% of the carotenoid content, respectively, were the primary carotenoids in N1 and N1Y fruits during the M2 period (Figure 2). At the M3 period, the main carotenoid components of N1 fruits were phytoene, zeaxanthin palmitate, and zeaxanthin dipalmitate, which made up 59.6% of the total carotenoid content. However, the primary carotenoids in N1Y fruit were zeaxanthin, phytoene, and zeaxanthin palmitate, which accounted for 87.9% of the carotenoid content ( Figure 2). Liu, Zeng, Sun, Wu, Hu, Shen and Wang [17] compared the carotenoid accumulation in red fruit (Lycium barbarum L.) and black fruit (Lycium ruthenicum Murr.). They found that the principal carotenoids in red and black fruit at the green fruit stage (S1) are xanthophyll, violaxanthin, and β-carotene, but, as the fruit ripens, these levels decrease. In addition, the zeaxanthin content of red fruit significantly increases from color break (S2) to the ripe fruit stage (S4) but is undetectable in the black fruit. The present study further explores the carotenoid changes in goji fruit during ripening. As the N1 and N1Y fruits ripened, we observed an increase in phytoene, zeaxanthin, zeaxanthin dipalmitate, zeaxanthin palmitate, beta-cryptoflavin palmitate, and rubixanthin palmitate. Moreover, the amounts of beta-carotene and xanthophyll in N1 and N1Y fruit considerably reduced from M1 to M3 (Figure 2).
Fruit color, an essential appearance and quality character, is determined by the buildup of different pigment compounds. The most significant among them is carotenoid accumulation, which allows plants to take on yellow, orange, and red colors [10,13,[19][20][21][22]. Most goji cultivars have red fruit color, while a few are black, yellow, purple, and white ( Figure S3). So far, the cause of the color variation between different goji cultivars is still unknown. Some studies suggest that the color difference between red and black fruits is due to the accumulation of anthocyanins in black fruits and carotenoids in red fruits [2,23]. In the present study, the M2 and M3 periods were critical for the difference in carotenoid content and fruit color between N1 and N1Y fruits (Figure 1a,b). During the M2 and M3 periods, zeaxanthin in the N1 fruits was largely converted to zeaxanthin palmitate and zeaxanthin dipalmitate, while the N1Y fruits continued to accumulate zeaxanthin. Thus, the color of the N1 fruit gradually turned red, while the N1Y fruit color appeared yellow from the M2 to M3 period (Figure 1c).
The genes involved in carotenoid metabolism have been identified and investigated in several plants [2,24]. The genes involved in plant carotenoid biosynthesis usually have multiple copies that exhibit tissue-specific expression, such as PSY1, PSY2, and PSY3, which are specifically expressed in the fruits, leaves, and roots of tomatoes and citrus, respectively [25,26]. All tissues express PSY-A, whereas PSY-B is found only in the leaves and roots of watermelons [19]. Most citrus species have two PDS genes and more than three ZDS genes [27]. At least two LCYBs and two LCYEs are present in tomatoes, papaya, and citrus. The LCYB1 gene of tomato is highly expressed in vegetative tissues, while the LCYB2 gene is significantly expressed in fruits and flowers [28][29][30][31]. However, few reports about goji carotenoid metabolism exist. Zhao, et al. [32] cloned the LcPSY, LcPDS, LcZDS, LcLCYB, LcLCYE, LcCHXE, LcZEP, and LcNCED from Lycium chinense, as those genes were constitutively expressed in the leaves, flowers, and fruits. Liu, Zeng, Sun, Wu, Hu, Shen and Wang [17] found that similar enzyme activities of PSY1, CYC-B, and CRTR-B2 were observed in red fruits (Lycium barbarum L.) and black fruits (Lycium ruthenicum Murr.), suggesting that the undetectable carotenoid content in black fruits was not due to inactivation of carotenoid biosynthetic enzymes. In addition, expression of carotenoid cleavage dioxygenase 4(CCD4) was higher in black fruit than in red fruit.
In the current study, we used RNA sequencing to evaluate the gene expression profiles of N1 and N1Y fruits at M1, M2, and M3 periods (Figure 3a). A total of 35,615 genes were expressed in N1 and N1Y fruits during fruit ripening. Based on the gene expression levels, we conclude that PSY1, PDS1, Z-ISO, ZDS16, crtISO4, LCYE, LCYB, CYP97A2, CYP97A3, CrtR2, VDE2, ZEP3, and NSY may play vital roles in the synthesis of carotenoids during goji fruit ripening (Figure 6a, Table S6). Furthermore, PSY1 and ZDS16 were correlated more strongly with the total carotenoid content of goji fruit than PDS1, Z-ISO, crtISO4, LCYE, LCYB, CYP97A2, CYP97A3, CrtR2, VDE2, ZEP3, and NSY (Figure 6b, Table S7). However, during fruit ripening, there was no significant difference in the expression of PSY1 and ZDS16 between N1 and N1Y (Figure 6a). Thus, variations in the expression levels of PSY1 and ZDS16 may not be the cause of the variations in carotenoids and color between N1 and N1Y fruits. Since it is unclear how zeaxanthin converts to zeaxanthin palmitate and zeaxanthin dipalmitate in goji fruit, the mechanism needs further investigation.
The regulation of plant carotenoids is challenging due to changes in concentration and composition during development. Numerous studies have demonstrated that TFs regulate carotenoid accumulation by influencing plant growth [2]. In tomato fruit, the MADS-box transcription factor RIN influences carotenoid accumulation by directly upregulating PSY1, ZISO, ZDS, and CRTISO [33,34]. Zhang, et al. [35] reported that SICMB1 regulates ethylene production and carotenoid accumulation during tomato fruit ripening via interactions with SIMADS-RIN, SIMADS1, SIAP2a, and TAGL1. In dark-grown Arabidopsis, phytochromeinteracting factor 1 (PIF1) can bind the PSY promoter and inhibit carotenoid accumulation [36]. Similarly, ethylene-responsive transcription factor AtPAP2.2 can bind to the promoters of PDS and PYS and influence carotenoid biosynthesis [11]. Yuan, et al. [37] indicated that SIWRKY35 could positively regulate the synthesis of tomato carotenoids by directly inducing SIDXS1 to rewire metabolism toward the MEP pathway. Hu, et al. [38] hypothesized that bHLH, MYB, NAC, MADS, and WRKY are involved in the biosynthesis of carotenoids and anthocyanins in Coffea arabica. There are few studies on the TFs responsible for carotenoid production in goji fruits. Yin, et al. [39] reported that two R2R3-MYBs and two BBXs may play key roles in the carotenoid biosynthesis of goji. In the present study, we analyzed the gene expression of N1 and N1Y fruit during fruit ripening by WGCNA, and seventeen WGCNA modules were identified (Figure 7a, Table S8). Among them, brown modules include PSY1 and ZDS16 (Table S8), and strongly correlate with total carotenoid, zeaxanthin, phytoene, xanthophyll, zeaxanthin dipalmitate, zeaxanthin palmitate, beta-cryptoflavin palmitate, and rubixanthin palmitate (Figure 7a). Thus, we conclude that many TFs in the brown module regulate carotenoid synthesis in goji fruit. We identified 63 TFs from the brown module (Table S10) that may be divided into 2 groups based on expression levels ( Figure 8). According to the TF expression, network weight values, and correlation with PSY1 and ZDS16, we conclude that one MADS-box TF (CMB1-1) and two WRKY TFs (WRKY22-1 and WRKY22-3) may favorably influence the synthesis of carotenoids by regulating PSY1 and ZDS16. One ethylene-responsive TF (RAP2-13-like) may negatively control the synthesis of carotenoids by PSY1 and ZDS16 (Figure 9).

Transcriptomic Analysis
The transcriptomic analysis method used in this study was modified from a previous study [23]. In brief, the clean reads were mapped to the genome of goji [40] using HISAT2 [41]. Differential expression was analyzed using DESeq 2 [42], and genes that met the criteria |log2 fold change| ≥ 1 and false discovery rate (FDR) < 0.05 were classified as differentially expressed. The raw RNA-seq data are freely available from the National Center for Biotechnology Information (accession: PRJNA946044).

Weighted Gene Co-Expression Network Analysis (WGCNA)
R package WGCNA [43] was used for identifying highly co-expressed gene modules with carotenoid content. WGCNA network construction and module detection were performed using an unsigned topological overlap matrix (TOM), a power β of 7, a deep split of 2, a minimum module size of 30, and a branch merge cut height of 0.25. The module eigengene value was calculated and used to evaluate the association of modules with carotenoid content in the eighteen samples.

Gene Expression Trend Analysis
R package Muffz (10.18129/B9.bioc.Mfuzz) was used to analyze gene expression trends in all samples. The c parameter was set to 6, and M-estimate was used for m.

Conclusions
In the present study, the difference in carotenoid biosynthesis and color formation in N1 (red fruit) and N1Y (yellow fruit) during fruit ripening were revealed by transcriptomics and carotenoid profiling. Finally, 27 carotenoids were identified in N1 and N1Y fruits at the M1, M2, and M3 periods. Among them, M2 and M3 were the critical periods for carotenoid and color differences between N1 and N1Y fruit. The zeaxanthin in the N1 fruits was significantly converted to zeaxanthin palmitate and zeaxanthin dipalmitate from M2 to M3, while the N1Y fruits continued accumulating zeaxanthin. As a result, N1 fruits progressively take on a reddish hue, while the color of N1Y fruits appears yellow from M2 to M3. However, the mechanism of conversion of zeaxanthin into zeaxanthin palmitate and zeaxanthin dipalmitate in goji fruits necessitates further investigation. Our findings from WGCNA analysis, gene trend analysis, and correlation analysis suggest that PSY1 and ZDS1 may be crucial genes in the synthesis of carotenoids during goji fruit ripening. Additionally, there are four TFs (CMB1-1, WRKY22-1, WRKY22-3, and RAP2-13-like) that may be involved in the carotenoid biosynthesis of goji fruit by regulating PSY1 and ZDS1. This work provides a new understanding of the molecular network of carotenoid synthesis and explains the variations in carotenoid accumulation in different colored goji fruits.
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/plants12152791/s1, Figure S1: The carotenoid standards and the raw data for quantification including calibration curves; Figure S2: The q-PCR validation of RNA-seq data; Figure S3: Fruit color of different wolfberry varieties; Table S1: Carotenoids dected in N1 and N1Y fruit at M1, M2, and M3; Table S2: The principle component analysis (PCA) of carotenoid, the samples from N1 and N1Y at M1, M2, and M3; Table S3: The RNA-sequencing data quality; Table S4: The genomic alignment of RNA-sequencing data; Table S5: List of DEGs between N1 and N1Y at different period; Table S6: Expression of genes related to carotenoid synthesis; Table S7: Carotenoids and genes correlation coefficient; Table S8: Genes list of 17 modules; Table S9: Gene expression trend analysis; Table S10: List of candidate transcription factors in the brown module; Table S11: Transcription factors and genes correlation coefficient; Table S12: The weight vlue between candidate transcription factors and genes in the brown module; Table S13: Primers for qRT-PCR in this study.

Data Availability Statement:
The raw RNA-seq data are freely available from the National Center for Biotechnology In-formation (accession: PRJNA946044).