Quantification of chlorophylls and carotenoids
Different color performance of plant organs and tissues is mainly caused by the pigment difference [16]. Content of chlorophylls and carotenoids in four watermelon materials was measured. In accordance with the visual appearance, material WM102 which shows a dark green skin owns the highest content of chlorophylls (chlorophyll a and chlorophyll b) (Fig. 2). The other ones, ranging from high to low chlorophyll content are WT2 (light green rind), WT13 (light green-gray rind), and WT4 (yellow rind) in turn. Compared with the other three materials, the chlorophyll content in WT4 was significantly reduced. In contrast to chlorophylls, content of carotenoids in WT4 is considerably higher than that in the other three materials which show a similar carotenoids content (Fig. 2). These result indicates that different rind color appearance in watermelon was due to the varying content of chlorophylls and carotenoids.
Ultrastructural changes in the rind of watermelon
TEM of chemically fixed rind samples could help to reveal the ultrastructural differences. The most noticeable difference focuses on the chloroplast number and structure. Material WM102 with dark green rind contains the most chloroplasts with approximately eight chloroplasts in each field of a microscope (Fig. 3A), and chloroplast granas of WM102 are composed of a set of compact grana lamellae (Fig. 3E). In WT2 and WT13, the number of grana and grana lamellae decreased slightly, but the structure of chloroplasts was not significantly changed compared with that in WM102 (Fig. 3B, 3C, 3F, 3G). In the yellow-rind material WT4, the number of grana is significantly decreased, and the number of grana lamellae in each grana is also seriously decreased (Fig. 3D, 3H). Except for the number of grana and grana lamellae, the structure of grana also changed considerably. Because lacking grana lamellae, the size of the grana is observed to be smaller, and the shape appears to be a ball (Fig. 3H), which is obviously different from that in WM102, WT2, and WT13. The ultrastructural results illustrate that the difference in chloroplast number and structure may be the a reason for the different pigment content and rind color appearance.
De novo assembly of RNA-Seq reads
In the previous research, the gene responsible for yellow skin was delimited to the region of 1–91.42 kb on chromosome 4 of watermelon, but no candidate genes were found in this region [15]. We proposed that the target gene may not be assembled into the watermelon reference genome of cultivar 97103 [17] or cultivar Charleston Gray [18]. To provide possible reference information for the related gene study, the clean reads were de novo assembled and 542,314,676 clean reads were obtained after the 586,255,040 raw reads were screened (Additional file 1: Table S2). The 202,841 transcripts assembled with Trinity software [16] has a length of 325,191,335 bp (Additional file 1: Table S3) and 84,516 unigenes were obtained based on the assembled transcripts. Among the unigenes, 77.61% (65,595) are shorter than 1,099 bp, 20.90% (17,663) ranges from 1,099 to 4,999 bp, and only 1.49% (1,258) are longer than 4,000 bp (Additional file 1: Table S4). Compared with the 22,596 and 22,567 genes of the watermelon reference genome of cultivar 97103 and cultivar Charleston Gray, respectively, many more novel genes were obtained through de novo assembly. The newly assembled transcriptome would provide a reference resource for gene functional study in watermelon.
Functional annotation of the unigenes
After assembly, unigenes were aligned to the six functional databases (Additional file 2: Fig. S1; Additional file 1: Table S5). To further evaluate the BLAST result, the E-value and similarity distributions to the Nr database were calculated. Results revealed that about 39.05% of the sequences were matched with the database (E-value<10−45) (Fig. 4a) and 67.92% of the sequences had similarities over 80% (Fig. 4b). Species distribution analysis showed that Cucumis melo and Cucumis sativus are highly similar with BLASTx hit score of 24.66% and 21.17%, respectively (Fig. 4c).
In GO classification, 28,095 unigenes (33.12%) were assigned to 65 GO term subcategories. The molecular function process contained 15 (23.08%) subcategories, and cellular component (32.31%) and category biological process (44.62%) are the three mostly enriched GO terms (Additional file 3: Fig. S2). The Cellular Component category contains 21 subcategories, and the Cellular Component category includes the most enriched subcategories. The top five enriched GO terms were cell, cell part, membrane, organelle, and membrane part (Additional file 3: Fig. S2). In the Biological Process Process, 29 subcategories were found, and the most enriched term was metabolic process followed by cellular process (Additional file 1: Table S6, Additional file 3: Fig. S2). In the eukaryote-specific COG (KOG) classification, 40,522 unigenes were classified into 26 KOG categories (Additional file 4: Fig. S3). General function prediction accounted for the largest proportion (22.57%) followed by posttranslational modification, protein turnover, chaperones (9.80%) and signal transduction mechanisms (8.16%) (Additional file 4: Fig. S3; Additional file 1: Table S7). 21352 (25.23%) unigenes were mapped into 33 categories of pathway hierarchy level 2 and 5 categories of hierarchy level 1 (Additional file 5: Fig. S4). Among the 33 pathways, translation (10.84%) and carbohydrate metabolism (10.32%) were the most dominant groups followed by signal transduction (9.78%), overview (8.40%), energy metabolism (8.28%), and folding, sorting and degradation (7.52%) (Additional file 5: Fig. S4, Additional file 1: Table S8).
Identification of DEGs and validation of RNA-Seq Data
Gene expression levels of the samples were estimated with the FPKM value and the correlation coefficients of the twelve samples were also analyzed and are listed in Additional file 1: Table S9. As shown in Fig. 5, comparison of gene expression between WT13 and WM102 showed that 5,905 genes were significantly differentially expressed, with 3,690 up-regulated and 2,215 down-regulated. For WT13 vs WT2, there were 5,811 DEGs, with 3,101 being up-regulated and 2.710 down-regulated. For WT13 vs WT4, 6179 genes were significantly differentially expressed, with 3486 being up-regulated and 2693 down-regulated. In the comparison group of WT2 and WM102, 3,954 up-regulated and 2,627 down-regulated DEGs were identified. For WT2 vs WT4, 8149 genes were significantly differentially expressed, with 4313 up-regulated and 3836 down-regulated. In the comparison group of WT4 and WM102, 1961 u-pregulated and 1111 down-regulated DEGs were identified (Fig. 5). In all the comparison groups, only 62 DEGs in common were found (Fig. 6). A total of 14 unigenes were selected to validate the RNA-Seq data using qRT-PCR with three biological replicates. The qRT-PCR results showed that the gene expression patterns of the 14 unigenes had similar trends with the RNA-Seq data with a positive correlation coefficient of 1.124 (Additional file 6: Fig. S5), implying the accuracy of the data.
Expression of genes involved in chlorophyll metabolism
In higher plants, chlorophyll metabolism, including chlorophyll biosynthesis, chlorophyll cycle, and chlorophyll degradation, is catalyzed by a series of enzyme complexes. Changes in the expression patterns of any of these reactions may affect rind color appearance. In this experiment, chlorophyll and carotenoid content also varied in WM102, WT2, WT13, and WT4 in accordance with the different color appearances. Accordingly, the core genes that encoded the enzymes involved in the chlorophyll metabolic pathway were investigated. A total of 50 unigenes encoding 23 chlorophyll metabolism-related enzymes were identified (Fig. 7; Additional file 1: Table S10).
Among these unigenes, 9 enzymes (i.e., GSA, HEMC, HEMG, CHLH/D/I, CHLG, CAO, and RCCR) were found to belong to a single gene family, and the remaining ones belong to multigene families. As shown in Fig. 7, almost all the genes related to chlorophyll metabolism were highest expressed in WM102, which has the darkest rind and the highest chlorophyll content and number of chloroplasts. In WT4, expression of CHLD is considerably down-regulated compared with materials WM102, WT2, and WT13, with a similar expression pattern as other material(s). In WT13, the expression of the HEMD gene was down-regulated in the chlorophyll biosynthesis part but up-regulated in the chlorophyll degradation part of the SGR gene compared with materials WM102, WT2, and WT4. In WT2, the gene expression of HEME was up-regulated compared with the other three materials. Except for the above particular genes, the expression levels of all the related genes in these four materials appear to be similar (Fig. 7; Additional file 1: Table S10).
Expression of genes involved in carotenoid metabolism
Carotenoids can be biosynthesized by many living organisms, such as photosynthetic bacteria, cyanobacteria, algae, and higher plants. Carotenoids usually play vital roles in light harvesting, photoprotection, photomorphogenesis, nonphotochemical quenching, and lipid peroxidation [19, 20]. The various kinds and contents of carotenoids impart different colors to organisms, ranging from colorless to yellow, orange, and red. Because of the different coloration expression and carotenoids content in the four watermelon materials, the expression levels of genes involved in carotenoid metabolism were investigated.
A total of 80 unigenes encoding 18 carotenoid metabolism-related enzymes were identified (Fig. 8; Additional file 1: Table S11). Among these unigenes, only one gene (CCS) was found to belong to a single gene family, and the remaining genes were demonstrated to be from multigene families. Carotenoids are a diverse class of lipid-soluble isoprenoids built from the 5-carbon compound IPP [21, 22]. Biosynthesis of isoprenoid is limited by the first committed step catalyzed by DXS [23]. Six members were found in the DXS family. DXS1 and DXS5 were up-regulated in WT13 and WM102, respectively. Gene DXR catalyzes the first committed step of the MEP pathway for isoprenoid biosynthesis in plants [24]. The expression levels of DXR1 and DXR2 in WT13 are considerably higher than those of the other three materials. As the terminal gene of the MEP pathway, LytB can enhance the production of carotenoids [25]. Both LytB (LytB1 and LytB2) were down-regulated compared with the other three materials. Four members exist in the IPI family. A varied expression level was exhibited by materials WT4, WM102, WT2, and WT13 in genes IPI1 and IPI2. As a key enzyme in carotenoid biosynthesis, PSY belongs to the branching enzyme that determines the flux of carbon source toward carotenoids [26, 27]. The expression level of PSY1 in WT13 is up-regulated. Xanthophylls are produced from β-carotene with reactions of epoxidation-catalyzed ZEP [28]. In contrast to the similar expression level in WT13, the expression of the ZEP1 gene was significantly up-regulated (Fig. 8; Additional file 1: Table S11). Except for the mentioned genes with different expression levels above, similar expression levels were presented by most of the materials. Although the rind color of WT4 turns yellow and carotenoid content is higher than those of the other three green materials, no gene with significantly different expression levels was found.