Extensive transcriptome changes underlying the fruit skin color intensity variation in purple eggplant

Background: Fruit skin color intensity is one of the most important economic traits of purple eggplant. A wide diversity for fruit skin color intensity exists in purple eggplant and the accumulation of anthocyanins and chlorophylls of fruit skin mainly affected color intensity. However, limited information is available contributing to the molecular mechanisms underlying fruit skin color intensity variation in eggplant. Results: In the present study, variation of two purple eggplant advanced lines EP26 and EP28, with different fruit skin color intensity was investigated. Comparative transcriptome analysis of EP26 and EP28 identified a total of 2218 differential expressed genes (DEGs) at two different developmental stages. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis showed that these DEGs were mainly involved in flavonoid biosynthesis and photosynthesis. In addition, a total of 131 transcription factors including MYB , bHLH , WRKY , and NAC exhibited dynamic changes, which might be responsible for the variation of fruit pigments accumulation between EP26 and EP28. Conclusions: Taken together, these results expand our knowledge of molecular mechanisms underlying fruit skin color intensity variation in purple eggplant, which allowing for improvement of fruit coloration in eggplant breeding.

intensity exists within the group of dark purple colored cultivars, but intensive, uniform and lasting fruit color types are most preferred by consumers than the paler types, especially in China [3] .
The dark purple color of eggplant fruit skin has been contributed to the accumulation of both anthocyanins and chlorophylls [4,5] . Anthocyanins are an important class of flavonoids that represent a large group of plant secondary metabolites [6] . The anthocyanin biosynthesis pathway has been extensively studied. It is involved at least two types of genes: structural genes and regulatory genes which have been identified in many species [7][8][9][10] . Structural genes encode the enzymes directly participate in the helix-loop-helix (bHLH) transcription factors, and WD40 repeat proteins [11] . These factors form a ternary MYB-bHLH-WD40 (MBW) transcriptional complex, to activate the transcription of structural genes in the anthocyanin pathway. Many factors have been reported to regulate the anthocyanin biosynthesis, such as light, temperature, sucrose and plant hormones [12][13][14][15] .
Chlorophylls play an important role in solar light harvesting and energy transport to the reaction centers during photosynthesis [16] . Chlorophyll metabolic pathway can be divided into three phases: biosynthesis of chlorophyll a from glutamate, interconversion between chlorophyll a and b, and degradation of chlorophyll a into a non-fluorescent chlorophyll catabolite [17] . Almost all enzymes involved in the pathway are identified [18][19][20] . Several transcription factors have been reported as regulators of the chlorophyll biosynthesis pathway, such as LONG HYPOCOTYL5 (HY5), GOLDEN2-LIKE (GLK) [21][22][23] . However, limit information is available about the regulation of chlorophyll metabolism in eggplant.
In recent years, more researches have focused on the analysis of anthocyanin biosynthesis in eggplant. The structural genes and regulatory genes involved in the anthocyanin biosynthesis of eggplant fruit have been identified and characterized [10] . Jiang et al. (2016b) discovered that the anthocyanin biosynthesis of eggplant cultivar 'Lanshan Hexian' was regulated by light [24] . In addition, low temperature was reported to induce anthocyanin accumulation in eggplant [10] . However, the mechanism underlying fruit skin color intensity variation in eggplant remains poorly understood.
In the present study, comparative transcriptome analysis was conducted between two dark purple type materials of eggplant, EP26 and EP28, with different intensity of fruit skin coloration, resulting in a comprehensive data set for the identification of candidate genes responsible for the control of pigmentation intensity in dark purple eggplant. The results would provide better insights into the molecular basis underlying the intensity variation in dark purple eggplant.

Phenotypic characterization of EP26 and EP28
In the present study, two dark purple eggplant advanced lines EP26 and EP28, with different fruit peel color intensity were compared (Fig. 1a). Based on the lightness (L*) of fruit skin, the fruit color intensity of EP26 was deeper than that of EP28 both at DDA15 and DDA25 stage (Fig. 1b). The total content of anthocyanin and chlorophyll from the two different stages (DDA15 and DDA25) of EP26 and EP28 were measured, respectively. The content of anthocyanin of EP26 at DDA15 and DDA25 were higher than those of EP28 (Fig.   1c). However, the level of chlorophyll content of EP26 was lower than that of EP28 at both stages (Fig. 1d). In addition, a total of eight anthocyanin compounds in fruit skin of EP26 and EP28 at DDA25 were detected, among which delphinidin 3-O-rutinoside is the most abundant both in EP26 and EP28 ( Table 1).

Identification of differentially expressed genes (DEGs) between EP26 and EP28
In order to systematically study the phenotypic differences between EP26 and EP28, the transcriptomes of four eggplant fruit peel samples with three repeats (A: EP26 at DAA15, B: EP26 at DAA25, C: EP28 at DAA15, D: EP28 at DAA25) were obtained using RNA-Seq. In total, twelve libraries (two materials × two developmental stages × three biological replicates) were constructed. Differences in gene expression in the fruit peel of EP26 and EP28 at two developmental stages were assessed (A vs C and B vs D). DEGs were identified with the expression fold (log 2 Ratio ≥ 1) and false discovery rate ≤ 0.01 as the thresholds.
A total of 1456 DEGs (913 up-regulated and 543 down-regulated) and 1163 DEGs (821 upregulated and 342 down-regulated), were detected in A vs C and B vs D, respectively ( Fig.   2a). Overall, a total of 2218 DEGs were found between EP26 and EP28 at the two stages.
Of all these DEGs, 401 genes with 231 up-regulated and 111 down-regulated were shared at both stages (Fig. 2b).

Gene ontology and KEGG pathway analyses of differentially expressed genes
The functions of the DEGs detected in the two developmental stages between EP26 and EP28 were classified according to the GO database respectively. Go terms were divided into three functional categories: biological process, cellular component, and molecular function. In the biological process category, the DEGs of the two stages were most enriched in the oxidation-reduction process (GO: 0055114). In the molecular functional category, the highly enriched terms among three stages were ATP binding (GO: 0005524) and structural constituent of ribosome (GO: 0003735). In the cellular component category, the GO terms cell periphery (GO: 0071944) were most enriched in A vs C comparison, while the GO terms including photosystem I (GO: 0009522), photosystem I reaction center (GO: 0009538), chloroplast thylakoid membrane (GO: 0009535) and photosystem II (GO: 0009523) were highly enriched terms in B vs D comparison (Table S1).
KEGG pathway enrichment analyses were performed according to the DEGs in the two stages respectively, to further explore the biological functions. As shown in Figure 3

Genes involved in anthocyanin biosynthesis
The anthocyanin biosynthetic pathway is an extension of the general flavonoid pathway [6] .
In the present study, structural genes participated in each step of the anthocyanin biosynthesis pathway were identified. Compared to EP28, the predicted encoding genes of

Genes involved in chlorophyll metabolism
In current study, DEGs related to chlorophyll biosynthesis were identified based on KEGG pathway annotations. Compared to EP28, only four genes involved in chlorophyll biosynthesis including HEMA (encoding glutamyl-tRNA reductase), CHLM (encoding Mgprotoporphyrin IX methyltransferase), CRD (encoding Mg-protoporphyrin IX monomethylester cyclase), and PORA (encoding protochlorophyllide oxidoreductase A) were up-regulated in EP26 at DDA25 stage. However, none of DEGs involved in chlorophyll biosynthesis were identified between EP26 and EP28 at DDA15 stage (Table 2). Moreover, compared to EP28, the expression levels of DEGs related to photosynthesis were mostly up-regulated in EP26 at DDA25 based on KEGG enrichment. Particularly, a total of 11 unigenes encoding chlorophyll a-b binding protein, which transfer excitation energy to the reaction centers in photosystem I and photosystem II, showed increased expression levels in EP26 at DDA25 ( Table 2).

Identification of differentially expressed transcription factors
Transcription factors play important roles in regulating gene expression in anthocyanin synthesis. In the present study, a total of 131 TFs were found to be differentially expressed in EP26 at least one pairwise combination (A vs C and B vs D) (Table S2). These TFs could be classified into 29 families, among which bHLH constituted the largest group (20 genes), followed by MYB (16 genes), HD-ZIP (9 genes), ERF (8 genes), NAC (8 genes), and WRKY (7 genes) (Fig. 5a).
Among the 20 DEGs encoding bHLH TFs, only one DEG Sme2.5_03973.1_g0005 was shared by both stages, which was down-regulated in EP26 compared with EP28. The expression level of Sme2.5_00592.1_g00005, which is homologus to AtTT8 increased in the comparison of EP26 and EP28 at DDA25. Among the 16 unigenes encoding MYB TFs, two unigenes Sme2.5_05099.1_g00002 and Sme2.5_00492.1_g00008 were shared by both A vs C and B vs D comparison. Sme2.5_05099.1_g00002 was up-regulated and phylogenetic analysis showed that this unigene was homologous to Arabdopsis MYB113, which has been reported to be involved in the control of anthocyanin biosynthesis.

Validation of differential gene expressed genes by qRT-PCR
Based on the RNA-Seq results, gene expression of thirteen selected DEGs was further validated using qRT-PCR with gene-specific primers (Table S3). These genes included anthocyanin biosynthetic genes and regulatory genes, as well as genes involved in photosynthesis. A linear regression analysis showed an overall correlation coefficient of R 2 =0.7501, indicating the results of qRT-PCR were consistent with those of the transcriptome analysis ( Figure 6).

Discussion
Fruit peel color intensity is a key trait for dark purple type of eggplant. Despite the importance of eggplant worldwide as a food crop and considerable interest in eggplant fruit secondary metabolites for their nutritive value, limited research has been conducted to understand the molecular mechanism underlying the fruit skin color intensity variation [4] . In the present study, EP26 exhibits good and stable performance on fruit skin color intensity compared to EP28. Also, comparative transcriptome analysis identified a set of key genes for variation of fruit skin color intensity between EP26 and EP28. The analysis of identified DEGs, together with phenotypic characterization of EP26 and EP28 will provide insights into the molecular basis underlying the intensity variation in fruit skin pigmentation of eggplant.

Fruit skin color intensity variation between EP26 and EP28
Although the fruit skin color intensity of both two lines was deeper at DDA25 than those at DDA15, the color intensity of EP26 was deeper than that of EP28 at each stage. The anthocyanin contents of EP26 were also higher than those of EP28 at both stages.
However, the chlorophyll content of EP26 was lower than that of EP28 at DDA15 and DDA25, respectively, which suggests that the anthocyanin content variation rather than the chlorophyll content variation contributes to fruit skin color intensity variation at the physiological level. Moreover, eight kinds of anthocyanins were found both in EP26 and EP28. It has been reported that delphinidin 3-(p-coumaroylrutinoside) -5-glucoside is the major anthocyanin in Japanese type eggplant and delphinidin 3-rutinoside is the major anthocyanin in non-Japanese type eggplant [29,30] . In the fruit peel of EP26 and EP28, which both are non-Japanese type eggplant, the major anthocyanin identified is delphinidin 3-rutinoside, which is consistent with previous studies.

Structural genes involved in anthocyanin biosynthesis
Most of the structural genes involved in anthocyanin biosynthesis have been isolated and cloned in eggplant [10] . These genes have been divided into early biosynthesis genes including CHS, CHI, and F3H , while the late biosynthesis genes including F3′H, F3′5′H, DFR, ANS, and UFGT [6] . Previous study suggested that higher expression of CHS, DFR, and ANS in eggplant fruit with violet vs. white fruit color result in high vs. very low anthocyanin concentration [4] . The expression levels of SmCHS, SmCHI, SmF3′5′H, SmDFR, and SmANS, except for SmF3H, showed a positive correlation with anthocyanin accumulation in eggplant fruit peel of cv. Lanshan Hexian [10] . In addition, DFR, F3H, and UFGT were considered as key enzymes in eggplant fruit anthocyanin biosynthesis [31] . In our study, the up-regulation of almost all of the structural biosynthetic genes, from CHS to UGFT in EP26 compared to EP28 at DDA25 suggests that fundamental transcriptional regulation of anthocyanin biosynthetic pathways could be a major factor in EP26 for deeper fruit peel color intensity. Particularly, the expression levels of CHI (unigene Sme2.5_00188.1_g00020), F3H (unigene Sme2.5_00015.1_g00020), and UFGT (unigene Sme2.5_00228.1_g00013) increased in parallel with the accumulation of anthocyanins, indicating these three genes might be the key structural genes underlying fruit peel color intensity variation between EP26 and EP28.

Regulatory genes involved in anthocyanin biosynthesis
It has been revealed that R2R3-MYB and bHLH transcription factors, together with WDR proteins play an important role in regulating the transcription of structural genes of anthocyanin biosynthesis [2,4,6,12] . In eggplant, SmMYB1 and SmMyb C was found to have higher transcript levels in violet eggplant compared to white eggplant fruits [4] . In this study, a total of 16 MYB TFs with different expression patterns was identified between EP26 and EP28. One unigene Sme2.5_05099.1_g00002 was homologus to AtMYB113 which was involved in regulation of anthocyanin biosynthesis and its transcript level increased in accordance with anthocyanin accumulation, which suggests its possible involvement in regulating the structural genes of anthocyanin biosynthesis in EP26. It has been proved that overexpression of Sme2.5_05099.1_g00002 activated abundant anthocyanin accumulation in the regenerating shoots of purple eggplant cultivar 'Zichang' [26] . The expression pattern of another MYB TF Sme00021.1_g00013 showed a negative correlation with the anthocyanin concentration, implying a possible role in repressing anthocyanin synthesis in EP26.
Among the 20 DEGs of bHLH TFs identified in our study, one up-regulated bHLH, eight down-regulated and seven up-regulated, four down-regulated bHLHs were found between EP26 and EP28 at DDA15 and DDA25, respectively. However, only one down-regulated bHLH Sme2.5_03973.1_g0005 was shared by both stages. Furthermore, the expression level of Sme00592.1_g00005, an orthologue of TT8, was up-regulated in EP26 at DDA25.
TT8 was first reported to positively regulate the expression of anthocyanin biosynthetic genes in Arabidopisis, and was potentially involved in the regulation of anthocyanin biosynthetic genes in pak choi [32] and peanut [33] . These finding suggest that the bHLH family is involved in the regulation of anthocyanins synthesis by different expression pattern.
Transcription factors of WRKY, NAC, bZIP families have also been reported to play important roles in regulating anthocyanin biosynthesis. In apple, overexpression of MdWRKY11 promoted the expression of F3H, FLS, DFR, ANS, and UFGT and increased the accumulation of flavonoids and anthocyanin [34] . NAC TF was found to be highly upregulated in blood-fleshed peaches compared to non-red-fleshed peaches [35] . SmHY5 (bZIP family) was identified to function as anthocyanin regulation factor in response to light signal in eggplant [24] . These transcription factors are indirectly regulate anthocyanin biosynthesis via MBW complex. The details of transcription factors in our study were summarized in Table S2. Some TFs belong to the same family showed different expression patterns, indicating a complex transcriptional regulatory network.

DEGs involved in chlorophyll metabolism
In the present study, only four genes including HEMA, CHLM , CRD, and PORA involved in chlorophyll (Chl) biosynthesis were found down-regulated in EP28 compared to EP26 at DDA25. The HEMA gene catalyzes the initial enzymatic step of tetrapyrrole biosynthesis, which eventually leads to chlorophyll production [36] . However, the expression pattern of these genes was not associated with Chl content between EP26 and EP28. Notably, a total of 11 unigenes encoding chlorophyll a/b-binding protein was down-regulated in EP28 compared to EP26 at DDA25. Previous study showed that coordinated synthesis of chlorophyll and the chlorophyll a/b-binding proteins is critical to the development of functional light-harvesting complexes [37] . Considering anthocyanin biosynthesis often affected by light [6] , the extensive alteration in the expression of genes involved in photosynthesis possibly reduced the photosynthetic capacity, which in turn affects the biosynthesis of anthocyanin in EP28.

Conclusions
In this study, we investigated variation in fruit skin color intensity between EP26 and EP28. Higher anthocyanin contents and lower chlorophyll contents were observed in EP26 with deeper fruit skin color intensity. Comparative transcriptome analysis identified a total of 2218 differential expressed genes between EP26 and EP28 at two developmental stages. KEGG analysis showed that these DEGs were mainly involved in flavonoid biosynthesis and photosynthesis. Enhanced anthocyanin biosynthesis resulted in higher anthocyanin contents seems to be the major factor account for the deeper color intensity of fruit skin in EP26. Our findings provide insights into the molecular mechanism underlying the variation of fruit skin color intensity of purple eggplant.

Plant materials
Two dark purple eggplant advanced lines EP26 and EP28 were used in this study. Both of the two materials were obtained from Jiangsu Academy of Agricultural Science. EP26 and EP28 were both Chinese type with long fruit and purple calyx. The color intensity of EP26 fruit is deeper and more stable than that of EP28. Well-developed plants of EP26 and EP28 were grown in a greenhouse at Jiangsu Academy of Agricultural Sciences. Fruits were collected at two different developmental stages: young fruit, 15 days after anthesis (DAA15) and commercial maturity, 25 days after anthesis (DAA25). Three fruits from each plant and each stage were harvested. The fruits of each stage were divided into 3 groups, one for transcriptome sequencing, one for anthocyanin and chlorophyll content measurement, as well as HPLC analysis, and the other for quantitative RT-PCR analysis.
The fruit peels were collected and immediately frozen in liquid nitrogen and stored at -80 ℃ until use. Three independent biological replications were performed for each stage mentioned above.

Fruit skin color intensity measurement
The skin color intensity of fruits of EP26 and EP28 of each stage was analyzed according to Gao et al. (2016) [25] . The L* (lightness), a* (redness and greenness), and b* (yellowness and blueness) were measured using a hand-held spectrophotometer (CR-400, Minolta, Japan). The L* is an indicator of fruit skin color intensity, as lower lightness generally means deeper color. The measurement was performed with three fruits of each stage of EP26 and EP28. Three independent biological replications were performed.

Total anthocyanin and chlorophyll analysis
The total anthocyanins of EP26 and EP28 at two different stages (DAA15 and DAA25) were measured using spectrophotometric differential pH method according to   [26] . The extraction of chlorophyll of EP26 and EP28 at two different stages was performed according to   [27] . The analysis mentioned above was performed in three independent biological replicates.

ESI-Q TRAP-MS/MS analysis of anthocyanins
The sample preparation, anthocyanin identification and quantification were performed at Wuhan MetWare Biotechnology Co., Ltd. (www. metware.cn) following the standard procedures. LIT and triple quadrupole (QQQ) scans were acquired on a triple quadrupolelinear ion trap mass spectrometer (Q TRAP), API 6500 Q TRAP LC/MS/MS System, equipped with an ESI Turbo Ion-Spray interface, operating in a positive ion mode and controlled by Analyst 1.6 software (AB Sciex). The ESI source operation parameters were as follows: ion source, turbo spray, source temperature 500°C, ion spray voltage (IS) 5500 V, ion source gas I (GSI), gas II(GSII), curtain gas (CUR) were set at 55, 60, and 25.0 psi, respectively.
The collision gas (CAD) was high. Instrument tuning and mass calibration were performed with 10 and 100 μmol/L polypropylene glycol solutions in QQQ and LIT modes, respectively. QQQ scans were acquired as MRM experiments with collision gas (nitrogen) set to 5 psi. DP and CE for individual MRM transitions were done with further DP and CE optimization. A specific set of MRM transitions were monitored for each period according to the metabolites eluted within this period.

RNA extraction, library preparation and Transcriptome sequencing
The total RNA extracted from fruit peels of EP26 and EP28 was used to construct RNA-Seq libraries, and every sample was analyzed in three biological replicates. Total RNA of each sample was isolated using the RNeasy Plant Mini Kit (Qiagen, Hilden, Germany). RNA quantity, integrity and purity were analyzed using Nanodrop and Agilent Bioanalyzer 2100 (Agilent Technologies, CA, USA), respectively. The libraries were sequenced and analyzed by Biomarker Corporation (Beijing, China) on an Illumina HiSeq 4000 platform.

Bioinformatics Analysis
Filtered clean reads were mapped to an eggplant reference genome (http://eggplant.kazusa.or.jp). FPKM (fragments per kb per million fragments) was used to calculate unigene expression levels. The FDR (False Discovery Rate) control method was applied to determine the threshold of p value. An FDR < 0.01 and the absolute value of log2Ratio ≥ 1 were used as the threshold for the judgment of the significance of the gene expression differences. The different expressed genes were then subjected to GO and KEGG Ontology enrichment analysis.

Quantitative Real-Time PCR Analysis
A total of thirteen differentially expressed unigenes were selected for quantitative realtime PCR (qRT-PCR) analysis. First-strand cDNA was generated from the total RNA isolated from fruit peels of EP26 and EP28 at DDA15 and DDA25, respectively, using the ReverTran Ace qPCR RT kit (Toyobo). Specific primers were designed using Primer 3

Availability of data and materials
The datasets used and analysed during the current study are available from the corresponding author on reasonable request.      Quantitative real-time PCR analysis of the thirteen DEGs between EP26 and EP28

Consent for publication
at two different stages.

Supplementary Files
This is a list of supplementary files associated with the primary manuscript. Click to download. Table S2.xlsx   Table S3.xlsx Table S1.xlsx