Comparative Transcriptome Analysis Reveals Gene Expression Differences in Eggplant (Solanum melongena L.) Fruits with Different Brightness

Fruit brightness is an important quality trait that affects the market value of eggplant. However, few studies have been conducted on eggplant brightness. In this study, we aimed to identify genes related to this trait in three varieties of eggplant with different fruit brightness between 14 and 22 days after pollination. Using RNA-Seq Gene Ontology and Kyoto Encyclopedia of Genes and Genomes enrichment analyses, we found that wax- and cutin-related pathways and differentially expressed genes displayed significant differences among different development stages and varieties. Scanning electron microscopy revealed that the wax layer was thinner in ‘30-1’ and ‘QPCQ’ than in ‘22-1’. Gas chromatography-mass spectrometry analysis revealed that wax content was significantly lower in ‘30-1’ than in ‘22-1’, which indicated that wax may be an important factor determining fruit brightness. We further identified and analyzed the KCS gene family, which encodes the rate-limiting enzyme of FA elongation in wax synthesis. The results provide an insight into the molecular mechanisms of fruit brightness in eggplants and further eggplant breeding programs.


Introduction
The eggplant (Solanum melongena L.) is a solanaceous crop that is grown worldwide. Additionally, its fruits are an important component of the human diet, providing vitamins and active phytochemical compounds. Fruit brightness is a highly valuable external quality trait that affects the market value of eggplants. Compared with dull fruit, the bright appearance of glossy fruit attracts consumers [1][2][3].
Previous studies on fruit brightness have mainly focused on other crops such as cucumber, tomato, Arabidopsis, and navel orange. Early genetic studies have shown that the dull fruit skin phenotype (D) is dominant over the glossy fruit skin phenotype (d) in cucumber [4,5]. Yuan et al. [6] found that D was mapped between two markers, ME23EM4 and CS15, by constructing a genetic linkage map using 224 recombinant inbred lines (RILs). Miao et al. [7] constructed a new genetic linkage map using 148 RILs, where the d gene was found to be mapped between two simple sequence repeat (SSR) markers, SSR06003 and SSR15818, on chromosome 5. Yang et al. [5] preliminarily mapped markers SCZ69 and SSR16203 by combining bulked segregant analysis and 11 polymorphic molecular markers on chromosome 5. For further high-resolution mapping of D/d genes by increasing the F2 population, they identified Csa016880 or Csa016887 as candidate gene D in cucumbers. In tomatoes, cutin deficiency has been reported to contribute to fruit glossiness [8][9][10]. However, Petit et al. [11] found that cuticle architecture was responsible for fruit glossiness by analyzing ethyl methanesulfonate mutants in tomatoes. In Arabidopsis, stem glossiness is The '22-1', '30-1', and 'QPCQ' are all breeding lines selected by the Shanghai Academy of Agricultural Sciences (SAAS), Shanghai, China. In brief, 'Guangdong green eggplant' and '765 long eggplant' were crossed. Their offspring were self-crossed for several generations to obtain '22-1' and '30-1' varieties. 'QPCQ' was local material collected from farmers. The three varieties were grown at the Horticultural Research Institute of SAAS. To ensure consistent fruit development, all fruits were artificially pollinated on the day of flowering and harvested after between 14 and 22 days. The peel of the fruits was cut off and frozen in liquid nitrogen. Three replicates of each variety were collected and stored at −80 • C until use.
Fruit peel brightness was measured using an intelligent gloss meter (3nh, Shenzhen, China) according to the operating manual. Briefly, the "power" button was pushed to turn on the instrument and "calibration" was selected from the main menu to enter the calibration interface. The instrument was placed in the calibration box for calibration. After completion of the calibration, the measurement "basic mode" was selected, the measurement hole was aligned with the fruits to be measured, and "measurement" was clicked to obtain the value. Fruits were measured in the field, and each variety was measured with three replications.

RNA Extraction, Library Construction, and Illumina Sequencing
Total RNA was extracted using the MiniBEST Universal RNA Extraction Kit (TaKaRa, Tokyo, Japan) according to the manufacturer's instructions. Total RNA was quantified using a NanoDrop spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). The A260/280 ratios of all samples were above two. The 28S/18S ratio and RNA integrity number (RIN) values were determined using an Agilent 2100 system (Agilent, Santa Clara, CA, USA) [17]. The mRNA was purified from total RNA using poly T oligo-attached magnetic beads; an interrupting reagent was added to interrupt the mRNA and produce short fragments. The interrupted mRNA fragments were used as a template to synthesize onestrand cDNA with a six-base random primer. Subsequently, a two-strand synthesis reaction system was prepared to synthesize double-stranded cDNA. The purified double-stranded cDNA was repaired, poly A-tailed, and ligated with sequencing adapters. PCR amplification was performed after fragment size selection. The constructed libraries were qualified using an Agilent 2100 Bioanalyzer (Agilent) and sequenced using an Illumina HiSeq 2500. Sequencing generated 100-bp paired-end reads as raw reads. All of the generated raw sequencing reads were filtered to remove low-quality poly-N and low-quality reads using the software SOAPnuke (BGI). After filtering, the remaining reads were termed 'clean reads'. The data have submitted to NCBI SRA database with accession number PRJNA721241.

RNA-Seq Data Analyses and DEG Identification
The clean reads were mapped to the eggplant reference genome (http://eggplant. kazusa.or.jp/ accessed on 27 June 2022) by hisat2 (hierarchical indexing for spliced alignment of transcripts) [18], with default parameters. The results of the comparison between the clean reads and the reference genome were stored in binary 'bam' files. The expression value from the fragments per kilobase per million mapped (FPKM) method was quantified using the cufflinks [19]. When calculating the difference in gene expression, htseq-count [20] was used to obtain the number of reads in each sample. The estimateSizeFactors function of the DESeq (2012) R package [21] was used to standardize the data, and the nbinomTest function was used to calculate the p-value and fold-change value of the different samples. p ≤ 0.05 and |log2 fold-change| ≥ 1 were set as thresholds for significantly differential expression.

Annotation and Classification of DEGs
Gene ontology (GO) annotation was performed using the Blast2GO program [22] and GO functional classification with a Pearson chi-squared test was performed using Web Gene Ontology Annotation Plot (WEGO) [23]. The DEGs were mapped to GO terms according to the analyses, and the number of DEGs in each term was calculated [17]. DEGs were mapped into the Kyoto Encyclopedia of Genes and Genomes (KEGG) metabolic pathway database and annotated using BLASTX. A hypergeometric test was applied to the results of GO and KEGG enrichment analyses to identify significant GO terms, enriched metabolic pathways, or signal transduction pathways in DEGs compared with the whole genome background. An adjusted p-value of ≤0.05 was used to define significantly enriched GO terms and KEGG pathways [24].

Validation of Gene Expression Profile by qRT-PCR
An Evo M-MLV RT kit (Accurate Biotechnology, Changsha, China) was used to synthesize cDNA from total RNA according to the manufacturer's instructions. qRT-PCR was performed on an ABI 7500 (Applied Biosystems, Waltham, MA, USA) using a 20-µL reaction mixture, containing 10 µL Hieff ® qPCR SYBR Green Master Mix (Low Rox Plus), 0.4 µL of each primer (0.2 µM final concentration), and 2 µL cDNA template, which was diluted four times with ddH 2 O. The PCR amplification procedure consisted of denaturation at 95 • C for 5 min, followed by 40 cycles of denaturation at 95 • C for 10 s, annealing, and extension at 60 • C for 30 s. Transcript levels were calculated using the 2-∆∆Ct method [25]. The SmEF1a gene (Sme2.5_01406.1_g00001.1) was used as an internal control. Three replicates were performed for the entire experiment. The primers used are listed in Table S1.

Cutin Monomer and Wax Analysis
Cuticular waxes of eggplant were extracted from 50 cm 2 peels. The peels were dipped into 15 mL of chloroform for 5 min. Subsequently, 100 µg docosane was added to the solution as internal standard. Extracts were dried under moderate nitrogen flux and lipids were derivatized 10 min at 100 • C by adding 100 pyridine and 100 µL of BSTFA (derivatization agent). After cooling down to room temperature, extracts were dried with nitrogen to a constant weight. Further, 500 µL of n-hexane was added to resuspend the sample and draw all samples through a 0.22-µm membrane. Wax was then analyzed and quantified using gas chromatography-mass spectrometry (GC-MS). Cutin monomer analyses were conducted as previously described [9].

Cryo-SEM
Fruit peels at different development stages were cut off and freeze-fixed using liquid nitrogen. Samples were then transferred to a vacuum-sputtering instrument for fracture under freezing conditions to expose their fresh fracture surface. Sublimation was carried out according to the condition of the sample to prevent the surface from being wrapped in ice. The sample preparation was completed after conductive spraying. The samples were placed on the cold stage of a scanning electron microscope (SEM, Hitachi S-4800, Tokyo, Japan) for observation and photo capture. The thickness of cuticle was measured from the SEM photos using Image-Pro plus 6.0 software (Media Cybernetics, Silver Springs, MD, USA).

Identification of KCS Family in Eggplant
An HMM [26] analysis and a simple modular architecture research tool (SMART) [27] were used to identify KCS protein families in the different varieties of eggplant. The HMM profile was downloaded from the Pfam database (http://pfam.xfam.org/ (accessed on 27 June 2022)) to obtain the Pfam serial number (PF08392.12). The hmmsearch function in the HMMER software was used to detect domains contained in the target protein sequence with an e-value ≤ 1 × 10 −3 . The results of HMMER sequence alignment were screened to remove protein sequences whose alignment length was less than 45% of the domain length of the HMM model while retaining the longest alternative splicing sequence. All non-redundant protein sequences were further analyzed using SMART (http://smart.embl-heidelberg.de/ (accessed on 27 June 2022)) for examination, and the same genes were confirmed as family members.

Phylogenetic Analysis
The identified 19 SmKCS protein sequences were aligned using ClustalX software with default parameters [28]. The phylogenetic tree was constructed using the neighbor-joining method in MEGA7.0 with 1000 bootstrap tests [29].

Structure and Motif Analysis
Information on gene location and exon-intron structures was acquired from the reference genome annotation files. The MEME program (http://alternate.meme-suite.org/ tools/meme (accessed on 27 June 2022)) was used to identify the conserved motifs of the KCS sequence, while the maximum motif search value was set at 15 with an optimum motif width of 10-100 amino acid residues.

Data Analysis
SPSS software (IBM, Armonk, NY, USA)was used for data analysis and Origin 9.0 (OriginLab, Northampton, MA, USA) was used for drawing. All data are presented as the mean ± standard deviation (SD) of three biological replicates. Tukey's at p = 0.01 or p = 0.05 was used to identify the significant differences.

Identification of Three Independent Eggplant Varieties with Different Fruit Brightness
To explore the differences in fruit brightness in S. melongena, we selected varieties with dull and glossy peels by visual evaluation and gloss meter measurements. The fruit brightness values of varieties '30-1' and 'QPCQ', both of which display a glossy peel, were 6.40 and 10.57 on the 14th day after pollination (DAP), respectively, which represented the fruit expansion period, and 12.73 and 13.23 on the 22nd DAP, respectively, which represented the fruit maturation period. The brightness value of variety '22-1' (dull peel) was 4.67 and 3.13 on the 14th and 22nd DAP, respectively ( Figure 1).

Identification of Three Independent Eggplant Varieties with Different Fruit Brightness
To explore the differences in fruit brightness in S. melongena, we selected varieties with dull and glossy peels by visual evaluation and gloss meter measurements. The fruit brightness values of varieties '30-1' and 'QPCQ', both of which display a glossy peel, were 6.40 and 10.57 on the 14th day after pollination (DAP), respectively, which represented the fruit expansion period, and 12.73 and 13.23 on the 22nd DAP, respectively, which represented the fruit maturation period. The brightness value of variety '22-1' (dull peel) was 4.67 and 3.13 on the 14th and 22nd DAP, respectively ( Figure 1).

Fruit Transcriptome Sequencing and Analyses
To further illustrate the molecular mechanisms of fruit brightness in the three eggplant varieties, RNA-Seq libraries were generated with fruit peels of 14  . The Q20 and Q30 of the reads (the proportion of the number of bases with a quality value greater than 20 and 30, respectively, compared to the total number of bases in the raw reads) in these data were over 97% and 92%, respectively (Table S2).
After removing the adapter sequences, poly-N, and low-quality reads from the raw data, the quality of the clean data was assessed. Approximately 48,891,631  and 46,712,450 (22-22) clean reads were generated for the '

Prediction of New Transcripts
New transcripts exhibited no annotation information in the reference genome. These new transcripts may have been new splicing subtypes of known genes or new transcripts of unknown genes. In this study, 30,726 new transcripts were detected, including 26,158 coding transcripts and 4568 noncoding transcripts (Table S4).

Single Nucleotide Polymorphisms and INDEL Detection
Single nucleotide polymorphisms (SNPs) are DNA sequence changes caused by single nucleotide changes, resulting in genome diversity among species or individuals. An average of 53,391 SNPs were found in the 22-1, 30-1, and QPCQ samples compared to that in the reference genome ( Table 1). The distribution of SNPs was mainly focused on the exons, introns, and intergenic areas ( Figure S1). An insertion-deletion (INDEL) is the insertion or deletion of small fragments of less than 50 bp in the sample relative to the reference genome. Similar to that of SNPs, the distribution of INDELs was also concentrated in exons, introns, and intergenic areas ( Figure S4). Table 1. SNP variant type summary. A-G: number of SNP for A-G variation. C-T: number of SNP for C-T variation. Transition: number of SNP for A-G and C-T variations. A-C: number of SNP for A-C variation. A-T: number of SNP for A-T variation. C-G: number of SNP for C-G variation. G-T: number of SNP for G-T variation. Transversion: number of SNP for A-C, A-T, C-G, and G-T variations.

Analyses of DEGs
The FPKM fragments method [30] was used to calculate the transcript abundance of genes among the different samples. DEGs were identified by setting a threshold of |log2 fold-change | > 1 and p < 0.05. As shown in Figure 2A Figure 2E-H). After comparing the DEGs of the different varieties, we found 37 genes that were differentially expressed in these four groups ( Figure 2I).

GO Functional Enrichment Analyses of DEGs
To further illustrate the potential function of the DEGs, their function nal classes were evaluated using GO enrichment analyses. GO is classified into the following three major functional categories: molecular function, cellular component, and biological process. As shown in Figures 3A and S3A-F, the DEGs were mainly clustered in the metabolic processes, cellular processes, single-organism processes, responses to stimuli, and developmental processes of the biological process category. Cell, cell part, organelle, and membrane dominated the cellular component category in all groups. The most common terms in the molecular function category were catalytic activity, binding, nucleic acid binding transcription factor activity, transporter activity, and molecular transducer activity. To further analyze the classification of upregulated and downregulated DEGs, we found that the most significantly different terms of the three GO categories in the seven groups were

GO Functional Enrichment Analyses of DEGs
To further illustrate the potential function of the DEGs, their function nal classes were evaluated using GO enrichment analyses. GO is classified into the following three major functional categories: molecular function, cellular component, and biological process. As shown in Figures 3A and S3A-F, the DEGs were mainly clustered in the metabolic processes, cellular processes, single-organism processes, responses to stimuli, and developmental processes of the biological process category. Cell, cell part, organelle, and membrane dominated the cellular component category in all groups. The most common terms in the molecular function category were catalytic activity, binding, nucleic acid binding transcription factor activity, transporter activity, and molecular transducer activity. To further analyze the classification of upregulated and downregulated DEGs, we found that the most significantly different terms of the three GO categories in the seven groups were consistent with the previous DEG classification ( Figures 3B and S3G

KEGG Pathway Database Enrichment Analyses of DEGs
According to the results of the DEG analysis, we performed KEGG pathway classification and enrichment analyses. KEGG pathways are divided into the following seven categories: cellular processes, environmental information processing, genetic information processing, human diseases (animals only), metabolism, organismal systems, and drug development. The DEGs in the seven groups were enriched in cellular processes, environmental information processing, genetic information processing, metabolism, and organismal systems. Among them, metabolism accounted for the largest portion (Figures 4A-C and 5A-D).

KEGG Pathway Database Enrichment Analyses of DEGs
According to the results of the DEG analysis, we performed KEGG pathway classification and enrichment analyses. KEGG pathways are divided into the following seven categories: cellular processes, environmental information processing, genetic information processing, human diseases (animals only), metabolism, organismal systems, and drug development. The DEGs in the seven groups were enriched in cellular processes, environmental information processing, genetic information processing, metabolism, and organismal systems. Among them, metabolism accounted for the largest portion (Figures 4A-C and 5A-D).

DEGs Related to Cutin and Wax in Eggplants of Different Brightness
After mapping the upregulated and downregulated DEGs in the KEGG database, we analyzed the top 30 pathways that were most significantly enriched. We found that several pathways related to the biosynthesis, metabolism, and transport of cutin, and wax accounted for a significant portion. In a previous study, cutin and wax affected fruit brightness [13]. We further analyzed the expression pattern of genes related to cutin and wax.

DEGs Related to Cutin and Wax in Eggplants of Different Brightness
After mapping the upregulated and downregulated DEGs in the KEGG database, we analyzed the top 30 pathways that were most significantly enriched. We found that several pathways related to the biosynthesis, metabolism, and transport of cutin, and wax accounted for a significant portion. In a previous study, cutin and wax affected fruit brightness [13]. We further analyzed the expression pattern of genes related to cutin and wax.  Figure 4D-F). Subsequently, we selected all pathways related to cutin and wax in these three groups. As shown in Table 2, 66, 38, and 39 DEGs involved in 'biosynthesis of unsaturated fatty acids', 'cutin, suberine, and wax biosynthesis', 'fatty acid biosynthesis', 'fatty acid degradation', 'fatty acid elongation', 'fatty acid metabolism', and 'ABC transporters' pathways were selected for further analysis in the three groups, respectively. After Venn plot analysis, 24 DEGs were common in the brightness of the three groups and displayed similar expression patterns ( Figure 6A). A total of sixteen DEGs were upregulated and eight DEGs were downregulated in the three groups (Table 3). Four genes encoding KCSs, the rate-limiting step enzymes in fatty acid elongation, were differentially expressed across the three groups. Sme2.5_22772.1_g00001.1, which encodes ECR, was downregulated in all three groups. Three genes belonging to the cytochrome P450 family were upregulated. Sme2.5_03548.1_g00005.1, which encodes fatty acyl-CoA reductase (FAR), was downregulated, as was Sme2.5_00227.1_g00015.1. Five DEGs belonging to the ABC transporter family were significantly differentially expressed among the three groups. Three of these genes were upregulated, whereas two were downregulated ( Figure 6C and Table 3).  Figure 4D-F). Subsequently, we selected all pathways related to cutin and in these three groups. As shown in Table 2, 66, 38, and 39 DEGs involved in 'biosynth of unsaturated fatty acids', 'cutin, suberine, and wax biosynthesis', 'fatty acid biosyn sis', 'fatty acid degradation', 'fatty acid elongation', 'fatty acid metabolism', and 'A transporters' pathways were selected for further analysis in the three groups, respectiv After Venn plot analysis, 24 DEGs were common in the brightness of the three groups displayed similar expression patterns ( Figure 6A). A total of sixteen DEGs were upre lated and eight DEGs were downregulated in the three groups (  Table 3).    When the KEGG pathway categories of group 30-14 were compared to that of 22-14, the 'ABC transporter' and 'fatty acid degradation' pathways showed significant differences.  Figure 6B). To further compare the differences among the three eggplant varieties, we created a Venn plot. Two genes (Sme2.5_07770.1_g00002.1 and Sme2.5_29857.1_g00001.1) were commonly downregulated in the four groups ( Figure 6B,C). Sme2.5_07770.1_g00002.1 encodes ABC transporter C family member 8, while Sme2.5_29857.1_g00001.1 encodes KCS.

Cuticle Condition of Eggplants of Different Brightness
Using the analyses mentioned above, we identified several pathways and genes related to cutin and wax that displayed significant differences among the different eggplant varieties. We hypothesized that the cuticle of the eggplant fruits may have contributed to the brightness. To verify our hypothesis, we used SEM to observe the fruit peels of the different brightness varieties of eggplant. As shown in Figure S5, the cuticle of '22-1' displayed relatively more irregular and circular-shaped, dome-like structures than those of 30-1 and QPCQ. After further analysis of the side structure of the peels using SEM, we found that the cuticular membrane of '22-1' was thicker than that of '30-1' and 'QPCQ' (Figure 8A-C). The average cuticle thickness of '22-1' reached 7.76 μm, whereas that of '30-1' and 'QPCQ' reached only 1.65 μm and 3.02 μm, respectively ( Figure 8D). We then

Identification of SmKCS Gene Family in Eggplant
As wax may contribute to fruit brightness and KCS is the rate-limiting enzyme of FA elongation, we further identified the entire SmKCS gene family in eggplants. A total of 21 candidate genes were identified using the blastp method in NCBI and the hidden Markov model (HMM) search program (PF08392.12). The length of SmKCS proteins ranged from 204 to 530 aa, with an average length of 426 aa and an average molecular weight of 47897.2 kda. The pI values ranged from 6.31 to 10 (Table 4).

Phylogenetic Analysis of SmKCS Proteins
To study the relationships between eggplant SmKCS proteins, we constructed a phylogenetic tree with the SmKCS protein sequences using MEGA7.0. As shown in Figure 9, SmKCS proteins were categorized into four subgroups. Among the 21 SmKCS proteins, 5 belong to subgroup I, 3 belong to subgroup II, 6 belong to subgroup III, and 7 belong to subgroup IV. The four SmKCS genes we mentioned above displayed significant differences in RNA-Seq. Among them, Sme2.5_00826.1_g00007.1 and Sme2.5_00014.1_g00009.1 are classified in group I. Sme2.5_01347.1_g00005.1 and Sme2.5_29857.1_g00001.1 are classified in group III and group IV, respectively.

Structure and Conserved Motif Analyses of the SmKCS Gene Family
To examine the gene structure diversity of SmKCS genes in eggplants, we analyzed their exon-intron organization. The four SmKCS genes contained no introns, five SmKCS genes contained two introns, three SmKCS genes contained three introns, and one gene

Structure and Conserved Motif Analyses of the SmKCS Gene Family
To examine the gene structure diversity of SmKCS genes in eggplants, we analyzed their exon-intron organization. The four SmKCS genes contained no introns, five SmKCS genes contained two introns, three SmKCS genes contained three introns, and one gene contained four introns. The remaining eight genes each contained one intron. Sme2.5_00826.1_g00007.1 and Sme2.5_00014.1_g00009.1 contained one intron, while Sme2.5_01347.1_g00005.1 and Sme2.5_29857.1_g00001.1 contained two introns (Table S6 and Figure 10).
To further understand the structure of KCS genes in eggplants, conserved mo were analyzed using the MEME tool. As shown in Figure 10 and  (Figure 10).  To further understand the structure of KCS genes in eggplants, conserved motifs were analyzed using the MEME tool. As shown in Figure 10 and Table S7,

qPCR Verification
In order to validate the accuracy of RNA-seq data, we chose 15 genes related to wax biosynthesis, metabolism, and transportation, which were mentioned above in the 22 vs. 14 groups for qRT-PCR analysis (Table 3 and Figure 11). In total, 2 genes were downregulated and the other 13 genes were upregulated in the 22 vs. 14 groups. These results of the qRT-PCR analysis were similar to those of the RNA-seq data.

qPCR Verification
In order to validate the accuracy of RNA-seq data, we chose 15 genes related to wax biosynthesis, metabolism, and transportation, which were mentioned above in the 22 vs. 14 groups for qRT-PCR analysis (Table 3 and Figure 11). In total, 2 genes were downregulated and the other 13 genes were upregulated in the 22 vs. 14 groups. These results of the qRT-PCR analysis were similar to those of the RNA-seq data.

Discussion
Fruit brightness is an important commodity trait that affects consumer consumption. Research on fruit brightness has mainly focused on plants such as cucumbers, tomatoes, Arabidopsis, and navel oranges, but seldom on eggplants. We selected three varieties of eggplant with different brightness levels to investigate the dynamics and differential expression of genes using RNA-Seq. Across the 18 samples, the Q20 and Q30 of the reads Figure 11. Verification of the expression of selected DEGs by qRT-PCR.

Discussion
Fruit brightness is an important commodity trait that affects consumer consumption. Research on fruit brightness has mainly focused on plants such as cucumbers, tomatoes, Arabidopsis, and navel oranges, but seldom on eggplants. We selected three varieties of eggplant with different brightness levels to investigate the dynamics and differential expression of genes using RNA-Seq. Across the 18 samples, the Q20 and Q30 of the reads were higher than 97% and 92%, respectively (Table S2), reflecting high sequencing reliability. The average mapping ratio of the samples to the reference genome was 82% (Table S3) and was comparable among the different samples.
To further illustrate the potential function of the DEGs, we conducted GO and KEGG pathway enrichment analyses. After analyzing the most significantly enriched 30 pathways in KEGG pathway analysis, we found that 'biosynthesis of unsaturated fatty acids', 'fatty acid biosynthesis', 'fatty acid elongation', 'fatty acid metabolism', 'fatty acid degradation', 'cutin, suberine, and wax biosynthesis', and 'ABC transporter', which are related to biosynthesis and metabolism, and 'transport of cutin and wax', represented a significant portion of the pathways ( Figures 4D-F and 5E-H). This result was consistent with that of previous studies. One study on tomato mutants showed that glossiness and increased stiffness in fruit peels are associated with cutin deficiency [31]. Petit et al. [11] found that a more complex cuticle architecture could be responsible for fruit brightness. In Arabidopsis, wax load or an alteration in wax compounds contributes to stem glossiness [12]. The SEM images obtained in this study ( Figures S5 and 8A-C) indicate that '22-1' has a thicker and more irregular cuticle, which may have contributed to the dull appearance of '22-1' fruits.
We identified pathways related to cutin and wax in the three 22 vs. 14 groups. In total, 24 DEGs were selected as candidate genes using a Venn plot ( Figure 6A). The selected genes encode proteins including LACS, KCS, ECR, WSD, and FAR. Two genes, which encode KCS and ABCC proteins, were downregulated at different fruit development stages among the different eggplant varieties. ECR, which constitutes one of the FAE complexes, regulates the fourth step of the reduction reaction in FA elongation. TSC13/YDL015 was the first ECR-type gene isolated from yeast and is related to the synthesis of sphingolipids and the elongation of VLCFA [32]. The CER10 gene, which encodes the ECR protein, displays a homologous relationship with TSC13 in A. thaliana. Loss of function of the CER10 gene can lead to a decrease in leaf cuticle wax, seed triacylglycerol, and sphingolipids [32]. Sme2.5_22772.1_g00001.1, which encodes the ECR protein, was downregulated in the 22 DAP vs. 14 DAP groups in our study ( Figure 6C and Table 3). In contrast, Sme2.5_02292.1_g00008.1, which encodes a homologous protein of at5g57800, was upregulated in the 22 vs. 14 groups. At5g57800 is also referred to as CER3 and plays an important role in the alkane synthesis pathway [33]. Overexpression of the CER1 gene increases the content of very long-chain alkanes and reduces the permeability of the cuticle, thereby reducing the non-stomatal loss of water and enhancing drought tolerance in Arabidopsis [34]. CER1 interacts with the WAX2, CER3, and CYTB5 genes. The interaction of these genes in yeast can convert very long-chain acyl-CoAs into very long-chain alkanes through redox reactions [35]. KCS catalyzes the polymerization of malonyl-CoA and long-chain acyl-CoA. Malonyl-CoAs and long-chain acyl-CoAs exhibit strict substrate specificity and their types determine the rate of the cyclic reaction and the acyl chain length of the final acyl-CoA products [36]. In Arabidopsis, two groups of KCS genes have been identified; the FAE1 KCS gene group, of which a total of 21 have been discovered so far, and the ELO KCS gene group. Currently, four of these genes have been annotated, but their functions have not yet been studied. Among them, KCS1, KCS2, KCS6, KCS9, and KCS20 have been reported to regulate wax biosynthesis [36][37][38][39][40][41][42]. In our study, the expression of four KCS genes was significantly different between the three 22 vs. 14 groups ( Figure 6C and Table 3). The expression of Sme2.5_29857.1_g00001.1, which encodes KCS6, also displayed significantly differential expression in groups 30-14 vs. 22  . The ABC transporter located on the plasma membrane achieves the transportation of wax components across the plasma membrane. Five DEGs (three genes being upregu-lated and two genes being downregulated) belonging to the ABC transporter family were identified in three, 22 vs. 14 groups ( Figure 6C and Table 3). Sme2.5_07770.1_g00002.1, which encodes ABC transporter C family member 8, was downregulated at different fruit development stages among the three eggplant varieties. Based on the above analyses, we suggest that '30-1' and 'QPCQ' may affect wax synthesis by affecting the elongation of fatty acids and the export of wax components. The reduced wax content subsequently affects the cuticle structure of the peel, thereby generating a glossy appearance of the fruit (Figure 12). Studies in recent years have found that TFs play an important role in regulating growth and development in plants. Mutation of the Glossy 1 gene in maize causes cuticle wax on the leaves to nearly disappear. Its homologous genes cer1 and wax2 mutants in Arabidopsis also displayed cuticle wax defects on the stem [36,43]. In Arabidopsis, MAH1 was identified to catalyze the conversion of alkanes into secondary alcohols and ketones [44]. WIN1/SHN1, a TF related to wax biosynthesis, directly or indirectly activates the transcription of genes encoding enzymes in cutin and wax biosynthesis. Overexpression of WIN1 increases wax content and improves drought tolerance in Arabidopsis [45]. Arabidopsis TFs MYB16 and MYB106 participate in cuticle development [46]. MYB96 regulates the expression of KCS1, KCS2, KCS6, WSD1, CER2, and KCR1 genes [47]. MYB30, MYB41, and CFL1 act as positive regulators in cutin and wax biosynthesis and cuticle development [48][49][50]. In our present research, 41 TFs displayed significant differences in fruit development. Among them, seven genes belonged to the MYB family ( Figure 7C). After protein interaction analysis using the string website, we found that MYB30 (Sme2.5_00097.1_g00008.1) may participate in cuticle formation by interacting with KCS6 (Sme2.5_29857.1_g00001.1) and CER4 (Sme2.5_03548.1_g00005.1) ( Figure 7B). Sme2.5_00097.1_g00008.1 may influence wax formation by interaction with KCS6 and CER4, and this needs further verification in the future.

Conclusions
This study selected three eggplant varieties, '22-1', '30-1', and 'QPCQ', with different fruit brightness levels for comparative transcriptomic analysis. After DEG enrichment analyses of GO and KEGG, we found that many pathways and genes related to cuticle cutin and wax content significantly differed between varieties. In total, 24 DEGs encoding proteins, such as KCS, ECR, FAR, and ABC transporters, were selected as candidate genes Studies in recent years have found that TFs play an important role in regulating growth and development in plants. Mutation of the Glossy 1 gene in maize causes cuticle wax on the leaves to nearly disappear. Its homologous genes cer1 and wax2 mutants in Arabidopsis also displayed cuticle wax defects on the stem [36,43]. In Arabidopsis, MAH1 was identified to catalyze the conversion of alkanes into secondary alcohols and ketones [44]. WIN1/SHN1, a TF related to wax biosynthesis, directly or indirectly activates the transcription of genes encoding enzymes in cutin and wax biosynthesis. Overexpression of WIN1 increases wax content and improves drought tolerance in Arabidopsis [45]. Arabidopsis TFs MYB16 and MYB106 participate in cuticle development [46]. MYB96 regulates the expression of KCS1, KCS2, KCS6, WSD1, CER2, and KCR1 genes [47]. MYB30, MYB41, and CFL1 act as positive regulators in cutin and wax biosynthesis and cuticle development [48][49][50]. In our present research, 41 TFs displayed significant differences in fruit development. Among them, seven genes belonged to the MYB family ( Figure 7C). After protein interaction analysis using the string website, we found that MYB30 (Sme2.5_00097.1_g00008.1) may participate in cuticle formation by interacting with KCS6 (Sme2.5_29857.1_g00001.1) and CER4 (Sme2.5_03548.1_g00005.1) ( Figure 7B). Sme2.5_00097.1_g00008.1 may influence wax formation by interaction with KCS6 and CER4, and this needs further verification in the future.

Conclusions
This study selected three eggplant varieties, '22-1', '30-1', and 'QPCQ', with different fruit brightness levels for comparative transcriptomic analysis. After DEG enrichment analyses of GO and KEGG, we found that many pathways and genes related to cuticle cutin and wax content significantly differed between varieties. In total, 24 DEGs encoding proteins, such as KCS, ECR, FAR, and ABC transporters, were selected as candidate genes regulating the wax synthesis and transport in the 22 Table S1: The primers used in qRT-PCR; Table S2: Summary of transcriptome sequencing data; Table S3: Summary of reads mapping to the reference genome; Table S4: New transcript type statistics; Table S5: Cutin and wax content in varieties '22-1' and '30-1'; Table S6: The intron-exon structure of eggplant SmKCS genes; Table S7: Motifs in SmKCS.