Carotenoid Biosynthetic Genes in Cabbage: Genome-Wide Identification, Evolution, and Expression Analysis

Carotenoids are natural functional pigments produced by plants and microorganisms and play essential roles in human health. Cabbage (Brassica oleracea L. var. capitata L.) is an economically important vegetable in terms of production and consumption. It is highly nutritious and contains β-carotene, lutein, and other antioxidant carotenoids. Here, we systematically analyzed carotenoid biosynthetic genes (CBGs) on the whole genome to understand the carotenoid biosynthetic pathway in cabbage. In total, 62 CBGs were identified in the cabbage genome, which are orthologs of 47 CBGs in Arabidopsis thaliana. Out of the 62 CBGs, 46 genes in cabbage were mapped to nine chromosomes. Evolutionary analysis of carotenoid biosynthetic orthologous gene pairs among B. oleracea, B. rapa, and A. thaliana revealed that orthologous genes of B. oleracea underwent a negative selection similar to that of B. rapa. Expression analysis of the CBGs showed functional differentiation of orthologous gene copies in B. oleracea and B. rapa. Exogenous phytohormone treatment suggested that ETH, ABA, and MeJA can promote some important CBGs expression in cabbage. Phylogenetic analysis showed that BoPSYs exhibit high conservatism. Subcellular localization analysis indicated that BoPSYs are located in the chloroplast. This study is the first to study carotenoid biosynthesis genes in cabbage and provides a basis for further research on carotenoid metabolic mechanisms in cabbage.


Introduction
Cabbage (Brassica oleracea L. var. capitata L., 2n = 18) is a member of the family Cruciferae and one of the most economically important vegetable crops cultivated worldwide [1]. Cabbage is rich in dietary fiber, vitamin B1 (VB1), vitamin B (VB2), calcium, and iron. Most importantly, cabbage also contains β-carotene, lutein, and other antioxidant carotenoids. It is one of the best vegetables recommended by the World Health Organization [2,3].
Carotenoids are naturally occurring pigments mainly found in plants, algae, and photosynthetic bacteria. Carotenoid pigments are mostly C40 lipophilic isoprenoids, comprising eight isoprene units joined in a head-to-tail fashion, and the central unit has a reverse connection [4]. The number of conjugated double bonds determines the color of carotenoids; the more the number of conjugated double bonds, the deeper the red color [5]. So far, over 800 carotenoids have been found in nature [6]. Carotenoids are divided into two main groups based on whether they contain oxygen or not: carotenes and xanthophylls. Carotenes are composed of carbon and hydrogen atoms and include phytoene, lycopene,

Identification and Analysis of Orthologs between B. oleracea and A. thaliana
To identify the CBGs and protein sequences of cabbage, we performed BLAST search on the cabbage genome and proteome database using A. thaliana carotenoid biosynthetic genes and protein sequences, with a cutoff E-value ≤ 10 −10 and coverage ≥ 0.75. Syntenic orthologous genes between A. thaliana and B. oleracea were identified based on the collinearity of flanking genes sequence similarity (E-value ≤ 10 −20 ). The specific distribution of B. oleracea carotenoid biosynthetic genes (BoCBGs) on the chromosome was analyzed using MapChart software (http://mapinspect.software.informer.com/, accessed on 24 April 2021).

Non-Synonymous/Synonymous Substitution (Ka/Ks) Ratios of Gene Pairs among B. oleracea, A. thaliana, and B. rapa
PAML software [43] was used to calculate the ratio of Ka/Ks rates of all orthologous gene pairs to estimate the selection mode of CBGs between B. oleracea, A. thaliana, and B. rapa. Ka/Ks ratio less than 1 indicates pure selection, ratio equal to 1 represents neutral selection, and ratio greater than 1 indicates positive selection.

Expression Analysis of Carotenoid Biosynthetic Genes in B. oleracea
The expression pattern of BoCBGs was investigated using RNA-Seq data [44]. Seven tissues, including root, leaf, stem, flower, bud, callus, and silique of B. oleracea accession '02-12' were used for RNA extraction. Total RNA was isolated from the cabbage tissues using TIANGEN RNAprep Pure Plant Kit and then reverse-transcribed to cDNA using PrimeScript™ RT reagent kit (TaKaRa, Kyoto, Japan), according to the manufacturer's instructions.

Exogenous Phytohormone Spraying of Cabbage Leaves
The leaves of cabbage were sprayed with different phytohormones to explore their effect on CBGs expression. Four exogenous phytohormones: ETH, ABA, SA, and MeJA at concentrations of 100 mg/L, 50 mg/L, 100 µM, and 100 µM, respectively, were sprayed. The leaves were sampled at 2, 4, 6, 12, and 24 h after spraying. The control was sprayed with double-distilled water alone. The relative expressions of BoPSY.1, BoPDS.1, BoZDS, and BoLYC genes were determined using qRT-PCR. The primers for qRT-PCR are listed in Table S1.

Phylogenetic Analysis
PSY protein sequences from 24 plant species were analyzed. MEGA6.0 was used to construct the phylogenetic tree, and Poisson correction model was used for distance computation. Node support was assessed by 1000 bootstrap replicates.

Identification of Carotenoid Biosynthetic Genes in B. oleracea
Genes associated with the carotenoid synthesis pathway in A. thaliana were analyzed using KEGG pathway database and TAIR database. A total of 47 genes were implicated in carotenoid biosynthesis in A. thaliana, of which 21 genes were shown to participate in the MEP pathway and 26 gene-encoding carotenoid biosynthetic enzymes (Table 1). A total of 62 CBGs were identified in cabbage, and 9 AtCBGs (GGPS4, GGPS7, GGPS6, GGPS8, GGPS9, GGPS10, GGPS11, GGPS12, and LUT1) showed no cabbage orthologs. Among the 62 BoCBGs, 58 were syntenic orthologs of the AtCBGs, and 4 BoCBGs had no syntenic relationships ( Figure 1).   In this study, BoCBGs were divided into three sub-genomes (LF, MF1, and MF2). Among the 58 syntenic homologous genes, there are 27, 17, and 14 genes in LF, MF1, and MF2, respectively (Table 1). In addition, a collinearity relationship was observed between B. oleracea and B. rapa. In total, 43 pairs of carotenoid biosynthetic genes were mapped on chromosomes using Circos software ( Figure 1).

Genomic Distribution on Chromosomes
Based on the genome annotation file, 62 BoCBGs were mapped to nine chromosomes. Among them, 46 were located on the chromosomes and 16 on the scaffold. The distribution of the 46 genes mapped on the chromosomes is shown in Figure 2, with 6, 3, 10, 4, 7, 2, 5, 6, and 3 genes located on chromosome C01-C09, respectively.

Comparative Evolutionary Analyses of Orthologous Gene Pairs for Carotenoid Biosynthetic Genes
Ka/Ks value is the ratio of non-synonymous mutation (Ka) to synonymous mutation (Ks) of genes. In this study, PAML software [43] was used to calculate the Ka/Ks of homologous gene pairs to explore the evolutionary selection mode of CBGs in A. thaliana, B. oleracea and B. rapa. Among the BoCBGs or BrCBGs, there were three copies of AtPSY gene, and two copies of AtPDS3, AtLUT2, and AtCHY1 genes. Eventually, 20 orthologous gene pairs were selected from A. thaliana, B. oleracea, and B. rapa ( Figure 3). Relative to A. thaliana, the Ka/Ks values of carotenoid homologous genes in B. oleracea and B. rapa were roughly the same, although the Ka/Ks ratios of BoGGPS2, BoCHY1.1, BoCHY1.2, and BoCHY2.1 in B. oleracea genome were slightly higher than that of B. rapa genome. Also, BoLUT2.1 and BoLUT5.1 were lower than those in B. rapa. The Ka/Ks values shown in Figure 3 are all less than 1, indicating that these genes are subjected to purify selection in evolution.

Genomic Distribution on Chromosomes
Based on the genome annotation file, 62 BoCBGs were mapped to nine chromosomes. Among them, 46 were located on the chromosomes and 16 on the scaffold. The distribution of the 46 genes mapped on the chromosomes is shown in Figure 2, with 6, 3, 10, 4, 7, 2, 5, 6, and 3 genes located on chromosome C01-C09, respectively.

Comparative Evolutionary Analyses of Orthologous Gene Pairs for Carotenoid Biosynthetic Genes
Ka/Ks value is the ratio of non-synonymous mutation (Ka) to synonymous mutation (Ks) of genes. In this study, PAML software [43] was used to calculate the Ka/Ks of homologous gene pairs to explore the evolutionary selection mode of CBGs in A. thaliana, B. oleracea and B. rapa. Among the BoCBGs or BrCBGs, there were three copies of AtPSY gene, and two copies of AtPDS3, AtLUT2, and AtCHY1 genes. Eventually, 20 orthologous gene pairs were selected from A. thaliana, B. oleracea, and B. rapa ( Figure 3). Relative to A. thaliana, the Ka/Ks values of carotenoid homologous genes in B. oleracea and B. rapa were roughly the same, although the Ka/Ks ratios of BoGGPS2, BoCHY1.1, BoCHY1.2, and BoCHY2.1 in B. oleracea genome were slightly higher than that of B. rapa genome. Also, BoLUT2.1 and BoLUT5.1 were lower than those in B. rapa. The Ka/Ks values shown in Figure 3 are all less than 1, indicating that these genes are subjected to purify selection in evolution.

Expression Analysis of Orthologous for Carotenoid Biosynthetic Genes among B. oleracea, B. rapa and A. thaliana
In this study, 20 CBGs in A. thaliana were studied to explore the differential expression of homologous genes in B. oleracea and B. rapa genomes.     The tissue expression of important genes in the carotenoid pathway was analyzed to further understand the expression of CBGs in B. oleracea ( Figure 5). The results showed that BoGGPS1 was expressed in all organs, including roots, leaves, stems, flowers, buds, silique, and callus, while BoGGPS2 was only expressed in flower buds. The expression of BoPSY.1 was highest in callus and flower buds, followed by flowers, stems, leaves, and silique. Almost no expression of BoPSY.1 was detected in the roots. Meanwhile, the expression of BoPSY.2 was highest in flowers, followed by buds and callus, and almost no expression was detected in other tissues. BoPSY.3 was expressed in all tissues, with the highest expression detected in the stem. On the other hand, BoPDS3.1 and BoPDS3.2 were expressed in all tissues, however, BoPDS3.1 was mainly expressed in callus, while BoPDS3.2 was mainly expressed in the flowers. The expression of BoZDS was highest in the roots, while that of BoLYC was highest in the flowers. Gene expression analysis revealed an extensive variance between paralogs of each CBG.

BoPSY.1, BoPDS.1, BoZDS, and BoLYC Respond to Exogenous Phytohormone Treatments
Specific concentrations of ETH, ABA, SA, and MeJA were sprayed on cabbage leaves at the seedling stage to explore their effects on CBGs (BoPSY.1, BoPDS.1, BoZDS, and BoLYC). ETH at 100 mg/L promoted the expression of the four genes ( Figure 6). Specifically, the expression of BoPSY.1 increased at 2 h after treatment and was the highest at 12 h. Meanwhile, the expressions of BoPDS.1, BoZDS, and BoLYC were highest at 4 and

Evolution of PSY Genes
The protein sequences of PSY in 24 species were analyzed using MEGA 6.0 software to investigate the evolutionary relationship ( Figure 7 and Table S2). As shown in Figure 7, the phylogenetic tree was consistent with the taxonomic conclusion, indicating that PSY has a high conservatism, which provides a theoretical basis for plant genetic molecular evolution. The PSY proteins of 24 species were divided into two categories. Haematococcus clustered into one category alone. The protein cluster of the other plant PSY belonged to angiosperms and was further divided into five subgroups: rice, maize, and sorghum were grouped under Commelinidae; alfalfa, carrot, strawberry, and apple belonged to Rosidae; lilium brownii and narcissus belonged to the Liliidae; Dilleniidae was divided into two groups, one group includes Arabidopsis and cabbage (Cruciferae), papaya (Carica Linn), pumpkin, and watermelon (Cucurbitaceae), the other group includes persimmon and kiwifruit clustered with Asteridae; Asteridae includes tomato, chilli, medlar, sweet potato, cape jasmine, Osmanthus, and agastache. Notably, BoPSY.1 and BoPSY.3 clustered on one branch and were more closely related to AtPSY than BoPSY.2 (Figure 7). Motif-based sequence analysis showed that all PSY proteins derived from angiosperms contain six conserved motifs, but the Haematococcus PSY only contains four motifs, lacking motif3 and motif6.
in cabbage. After spraying 50 mg/L ABA, the expression of BoPSY.1 increased gradually and peaked after 24 h, while the expression of BoZDS increased first and then decreased. SA at 100 μmol/L promoted the expression of BoPSY.1, but had minimal effect on the other three genes. MeJA at 100 μmol/L increased the expressions of BoPSY.1 and BoPDS.1 to a certain extent, which then decreased and peaked after 24 h. Meanwhile, the expression of BoZDS gene peaked after 2 h of MeJA treatment. Collectively, these results indicate that ETH, ABA, and MeJA can promote CBGs expression in cabbage. The Y-axis and X-axis represent the relative expression level and the time course of exogenous phytohormone treatment, respectively. Leaves were sampled at 0, 2, 4, 6, 12, and 24 h after ETH, ABA, SA, and MeJA spraying. Data represent the mean ± SD of three technical repetitions.

Evolution of PSY Genes
The protein sequences of PSY in 24 species were analyzed using MEGA 6.0 software to investigate the evolutionary relationship ( Figure 7 and Table S2). As shown in Figure  7, the phylogenetic tree was consistent with the taxonomic conclusion, indicating that PSY has a high conservatism, which provides a theoretical basis for plant genetic molecular evolution. The PSY proteins of 24 species were divided into two categories. Haematococcus clustered into one category alone. The protein cluster of the other plant PSY belonged to angiosperms and was further divided into five subgroups: rice, maize, and sorghum were grouped under Commelinidae; alfalfa, carrot, strawberry, and apple belonged to Rosidae; lilium brownii and narcissus belonged to the Liliidae; Dilleniidae was divided into two groups, one group includes Arabidopsis and cabbage (Cruciferae), papaya (Carica Linn), pumpkin, and watermelon (Cucurbitaceae), the other group includes persimmon and kiwifruit clustered with Asteridae; Asteridae includes tomato, chilli, medlar, sweet potato, cape jasmine, Osmanthus, and agastache. Notably, BoPSY.1 and BoPSY.3 clustered on one branch and were more closely related to AtPSY than BoPSY.2 ( Figure 7). Motif-based sequence analysis showed that all PSY proteins derived from angiosperms contain six conserved motifs, but the Haematococcus PSY only contains four motifs, lacking motif3 and motif6.  The results of Plant-PLoc subcellular localization prediction indicated that BoPSY proteins are localized in the chloroplast. Subcellular location of BoPSYs in tobacco (Nicotiana benthamiana L.) was examined using Confocal Laser Scanning Microscopy (FV3000). All three fluorescence signals of BoPSY-GFP fusion vector (BoPSY.1-GFP, BoPSY.2-GFP, and BoPSY.3-GFP) were observed in the chloroplast, consistent with the prediction of subcellular localization (Figure 8).

Characterization of Carotenoid Biosynthetic Genes in B. oleracea
Whole genome duplication provides rich genetic material for the expansion of gene families or the evolution of new genes in plants [46,47]. A whole-genome triplication (WGT) event occurred in B. oleracea after its divergence from A. thaliana [44]. There are 47 genes associated with carotenoid synthesis in Arabidopsis. In this study, 62 CBGs were identified in cabbage, among which 58 were syntenic orthologs of Arabidopsis, and only four had no syntenic relationships, suggesting that they could have evolved from some duplicate genes of Arabidopsis. However, nine AtCBGs (GGPS4, GGPS7, GGPS6, GGPS8, GGPS9, GGPS10, GGPS11, GGPS12, and LUT1) showed no B. oleracea orthologs, indicating that they might have been lost during evolution of B. oleracea. Notably, B. oleracea and A. thaliana belong to the Brassicaceae family. Phylogenetic analysis revealed that the triplicated B. oleracea genome segments diverged from a common ancestor soon after the divergence of the A. thaliana and Brassica lineages [48,49]. The multiple copies of the BoCBGs that are syntenic to genes in A. thaliana were generated from the WGT. The fragment with the highest gene densities was in the sub-genome LF, the fragment with moderate gene densities was in the sub-genome MF1, and the fragment with the least

Characterization of Carotenoid Biosynthetic Genes in B. oleracea
Whole genome duplication provides rich genetic material for the expansion of gene families or the evolution of new genes in plants [46,47]. A whole-genome triplication (WGT) event occurred in B. oleracea after its divergence from A. thaliana [44]. There are 47 genes associated with carotenoid synthesis in Arabidopsis. In this study, 62 CBGs were identified in cabbage, among which 58 were syntenic orthologs of Arabidopsis, and only four had no syntenic relationships, suggesting that they could have evolved from some duplicate genes of Arabidopsis. However, nine AtCBGs (GGPS4, GGPS7, GGPS6, GGPS8, GGPS9, GGPS10, GGPS11, GGPS12, and LUT1) showed no B. oleracea orthologs, indicating that they might have been lost during evolution of B. oleracea. Notably, B. oleracea and A. thaliana belong to the Brassicaceae family. Phylogenetic analysis revealed that the triplicated B. oleracea genome segments diverged from a common ancestor soon after the divergence of the A.
thaliana and Brassica lineages [48,49]. The multiple copies of the BoCBGs that are syntenic to genes in A. thaliana were generated from the WGT. The fragment with the highest gene densities was in the sub-genome LF, the fragment with moderate gene densities was in the sub-genome MF1, and the fragment with the least genes was in the sub-genome MF2. There were 27, 17, and 14 genes located in LF, MF1, and MF2, respectively. The distribution of CBGs was in line with the gene separation status in the whole genomic sequences [44].
Further, we compared the non-synonymous/synonymous (Ka/Ks) ratio of orthologous gene pairs between A. thaliana, B. oleracea, and B. rapa to analyze the evolutionary relationship of orthologous gene pairs for CBGs among the three species (Figure 3). The Ka/Ks ratio is a measure of selective pressures acting on genes. There were no significant differences in orthologous gene pairs between A. thaliana-B. oleracea and A. thaliana-B. rapa lineages. Therefore, we speculated that BoCBGs and BrCBGs could have undergone similar negative selection. However, comparing the expression differences between orthologous gene pairs in CBG revealed functional variation of a few gene pairs in B. oleracea and B. rapa.

Exogenous Phytohormones Regulate the Expression of Carotenoid Biosynthetic Genes
Regulation of carotenoid biosynthesis in plants is a complex process that is regulated by multilevel factors. Exogenous phytohormones affect the accumulation of carotenoids. Ethylene (ET) plays a vital role in plant ripening and carotenoid accumulation. Thus, inhibiting ET synthesis in fruits can hinder fruit ripening and lycopene accumulation [50]. Abscisic acid (ABA) is one of the most important plant phytohormones. Its synthetic precursors include neoxanthin and violaxanthin. These precursors regulate ABA synthesis, which in turn mediate the synthesis of carotenoids in plants [51,52]. Salicylic acid (SA) and Methyl Jasmonate (MeJA) also play an essential role in plant growth and stress resistance.
In the present study, different concentrations of exogenous phytohormones (ETH, ABA, SA, and MeJA) were sprayed on cabbage leaves to explore their effects on CBGs expression. The expression of four CBGs (BoPSY.1, BoPDS.1, BoZDS, and BoLYC) was analyzed. The result showed that ETH, ABA, and MeJA could promote the expression of CBGs in cabbage at specific concentrations. SA at 100 µmol/L did not significantly affect the expression of the four genes, consistent with the report of Gao, who concluded that the carotenoid synthesis pathway is insensitive to SA [53]. These results will guide the adoption of exogenous phytohormone spraying to increase carotenoid contents in cabbage.

Phytoene Synthase (PSY) Is a Key Enzyme in the Carotenoid Biosynthetic Pathway
PSY is the first special enzyme in the carotenoid biosynthesis pathway, catalyzing the condensation of two GGPP molecules into phytoene [54]. Overexpression of PSY increased carotenoid content and substantially improved β-carotene synthesis in canola seeds, cassava roots, potato tubers, and endosperms [55][56][57]. In tomato, PSY1 is highly expressed in fruits, PSY2 is essential for carotenoid synthesis in leaf tissues, and PSY3 potentially functions in roots [32][33][34][35].
In this study, BoPSY.1 and BoPSY.3 clustered on one branch, with BoPSY.2 in a separate branch. However, the three BoPSYs clustered on a specific branch in Arabidopsis, where a single PSY gene was shown to regulate phytoene synthesis in all tissues [21]. Three genes encoding PSY enzymes were identified in B. oleracea and B. rapa after WGT event [58]. In this study, three BoPSYs genes in cabbage were all localized in the chloroplast, which is similar to maize PSY1 [59]. However, examination of tissue expression of BoPSYs showed that they are differentially expressed in cabbage tissues ( Figure 5). Thus, we speculate that the BoPSYs might function in different cabbage tissues, like the PSYs expression patterns in tomato [32][33][34][35]; however, further studies are needed to verify this hypothesis.

Conclusions
In total, 62 carotenoid biosynthetic genes were identified in cabbage. These genes in B. oleracea species underwent similar negative selection with those in B. rapa species. However, the differential expression pattern of the duplicated CBGs occurred after polyploidization.
The expression of BoCBGs can be regulated by phytohormones. ETH, ABA, and MeJA can promote the expression of BoPSY.1, BoPDS.1, BoZDS, and BoLYC, which provides a basis for spraying phytohormones to improve the carotenoid content of cabbage. As the key enzyme in the carotenoid biosynthetic pathway, three BoPSY proteins were all localized in chloroplasts and exhibited different expression patterns in various cabbage tissues. This study uncovered carotenoid metabolic mechanisms in B. oleracea, proving a basis for developing cabbage varieties with high carotenoid content via genetic engineering.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/genes12122027/s1, Figure S1: Carotenoid synthesis pathway; Table S1: The primer sequences information used in this study; Table S2  Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: Data will be available on reasonable request.