Comparative Transcriptome Analysis Identified Key Pathways and Genes Regulating Differentiated Stigma Color in Melon (Cucumis melo L.)

Stigma color is an important morphological trait in many flowering plants. Visual observations in different field experiments have shown that a green stigma in melons is more attractive to natural pollinators than a yellow one. In the current study, we evaluated the characterization of two contrasted melon lines (MR-1 with a green stigma and M4-7 with a yellow stigma). Endogenous quantification showed that the chlorophyll and carotenoid content in the MR-1 stigmas was higher compared to the M4-7 stigmas. The primary differences in the chloroplast ultrastructure at different developmental stages depicted that the stigmas of both melon lines were mainly enriched with granum, plastoglobulus, and starch grains. Further, comparative transcriptomic analysis was performed to identify the candidate pathways and genes regulating melon stigma color during key developmental stages (S1–S3). The obtained results indicated similar biological processes involved in the three stages, but major differences were observed in light reactions and chloroplast pathways. The weighted gene co-expression network analysis (WGCNA) of differentially expressed genes (DEGs) uncovered a “black” network module (655 out of 5302 genes), mainly corresponding to light reactions, light harvesting, the chlorophyll metabolic process, and the chlorophyll biosynthetic process, and exhibited a significant contribution to stigma color. Overall, the expression of five key genes of the chlorophyll synthesis pathway—CAO (MELO03C010624), CHLH (MELO03C007233), CRD (MELO03C026802), HEMA (MELO03C011113), POR (MELO03C016714)—were checked at different stages of stigma development in both melon lines using quantitative real time polymerase chain reaction (qRT-PCR). The results exhibited that the expression of these genes gradually increased during the stigma development of the MR-1 line but decreased in the M4-7 line at S2. In addition, the expression trends in different stages were the same as RNA-seq, indicating data accuracy. To sum up, our research reveals an in-depth molecular mechanism of stigma coloration and suggests that chlorophyll and related biological activity play an important role in differentiating melon stigma color.


Introduction
Flower pollination is a critical stage during reproduction in flowering plants, where the male gamete (pollen grains) from the anther comes into contact with the female gamete (stigma), and successful reproduction is highly dependent on efficient pollinators [1]. The pollinators are usually described either as an abiotic factor (wind) or a biotic factor (insect) [2]. The distinction in floral color patterns at the center of the flower directs the WGCNA presented the key pathway and DEGs involved in the formation of colors in melon stigmas. The relative expressions of candidate identified genes were checked using quantitative real time polymerase chain reaction (qRT-PCR) and validated with RNA-seq data. We believe that our research findings would enrich any further genetic and genomic resource studies of stigma color formation in numerous crop plants.

Pigment Content and the Structural Basis of the Stigma
Visual flower inspection showed that the stigma of MR-1 was green from S1 to S3; however, the stigma of M4-7 was green from S1 to S2 stages and gradually became yellow from S2 to S3 stages ( Figure 1). The endogenous chlorophyll content was checked in the stigmas of MR-1 and M4-7 during S3 stage (Figure 2A), which demonstrated a high chlorophyll content (0.29 ± 0.004 mg/g) in the MR-1 line compared to the M4-7 line (0.03 ± 0.001 mg/g). The carotenoid content of the MR-1 line remained higher (0.28 ± 0.01 mg/g) than in the M4-7 line (0.04 ± 0.004 mg/g) ( Figure 2B), but the carotenoid/chlorophyll ratio was much lower (0.07) in MR-1 than in M4-7 (2.33). This primary difference may have caused the yellow stigma of M4-7, even though it contained less carotenoid content. Endogenous chlorophyll and carotenoid content was quantified in the stigmas of both melon lines, and differentiation in the chloroplast ultrastructure was also observed. Moreover, comparative transcriptomes of the stigmas of both melon lines were performed, and WGCNA presented the key pathway and DEGs involved in the formation of colors in melon stigmas. The relative expressions of candidate identified genes were checked using quantitative real time polymerase chain reaction (qRT-PCR) and validated with RNA-seq data. We believe that our research findings would enrich any further genetic and genomic resource studies of stigma color formation in numerous crop plants.

Pigment Content and the Structural Basis of the Stigma
Visual flower inspection showed that the stigma of MR-1 was green from S1 to S3; however, the stigma of M4-7 was green from S1 to S2 stages and gradually became yellow from S2 to S3 stages ( Figure 1). The endogenous chlorophyll content was checked in the stigmas of MR-1 and M4-7 during S3 stage (Figure 2A), which demonstrated a high chlorophyll content (0.29 ± 0.004 mg/g) in the MR-1 line compared to the M4-7 line (0.03 ± 0.001 mg/g). The carotenoid content of the MR-1 line remained higher (0.28 ± 0.01 mg/g) than in the M4-7 line (0.04 ± 0.004 mg/g) ( Figure 2B), but the carotenoid/chlorophyll ratio was much lower (0.07) in MR-1 than in M4-7 (2.33). This primary difference may have caused the yellow stigma of M4-7, even though it contained less carotenoid content.  The internal ultrastructures of both melon stigmas were observed by transmission electron microscopy ( Figure 3). The results showed that the chloroplasts of both materials began to differentiate during S1 stage, but the granum in MR-1 was higher in quantity and was well-developed. At S2 stage, the amount of MR-1 granum remained significantly higher than M4-7; however, a small number of starch grains were observed in M4-7. The increase in granum during the two stages was consistent with the stigma of MR-1 and M4-7, exhibiting a green stigma.
At S3 stage, the MR-1 chloroplast structure exhibited the highest amount of granum, along with plastoglobulus. In contrast, the M4-7 chloroplast structure showed a breakdown of granum, with a high production of plastoglobulus and starch grains, and the color of the stigma completely changed to yellow. In general, microscopic observations were consistent with the changes in stigma color. Therefore, we suggest that differences in chloroplasts provide the structural basis for differences in stigma color. The internal ultrastructures of both melon stigmas were observed by transmission electron microscopy ( Figure 3). The results showed that the chloroplasts of both materials began to differentiate during S1 stage, but the granum in MR-1 was higher in quantity and was well-developed. At S2 stage, the amount of MR-1 granum remained significantly higher than M4-7; however, a small number of starch grains were observed in M4-7. The increase in granum during the two stages was consistent with the stigma of MR-1 and M4-7, exhibiting a green stigma.   The internal ultrastructures of both melon stigmas were observed by transmission electron microscopy ( Figure 3). The results showed that the chloroplasts of both materials began to differentiate during S1 stage, but the granum in MR-1 was higher in quantity and was well-developed. At S2 stage, the amount of MR-1 granum remained significantly higher than M4-7; however, a small number of starch grains were observed in M4-7. The increase in granum during the two stages was consistent with the stigma of MR-1 and M4-7, exhibiting a green stigma.

Overview of Transcriptomic Data
To investigate the major differences in the transcriptomic dynamics during the stigma developmental stages of the MR-1 and M4-7 lines, we performed principal component analysis (PCA) based on average FPKM values ( Figure 4). The stigma transcriptome of both melon lines exhibited similar function and activity at S1 stage (G1 and Y1) and S2 stage (G2 and Y2), but substantial differences at S3 stage (G3 and Y3) were subsequently observed. Mainly, the first two stages (S1 and S2) of both melon lines exhibited a closer relationship with each other, and G1G2 and Y1Y2 can be seen in one cluster, respectively. However, the G3 and Y3 of S3 stage were highly inconsistent with the first two stages (S1 and S2). The findings of PCA analysis showed a consistent trend with the color transformation process of the melon stigma. Therefore, we believe that differences in the transcriptional level at S3 determined the developmental specificities and stigma color of the contrasted melon lines (MR-1 and M4-7).

Overview of Transcriptomic Data
To investigate the major differences in the transcriptomic dynamics during the stigma developmental stages of the MR-1 and M4-7 lines, we performed principal component analysis (PCA) based on average FPKM values ( Figure 4). The stigma transcriptome of both melon lines exhibited similar function and activity at S1 stage (G1 and Y1) and S2 stage (G2 and Y2), but substantial differences at S3 stage (G3 and Y3) were subsequently observed. Mainly, the first two stages (S1 and S2) of both melon lines exhibited a closer relationship with each other, and G1G2 and Y1Y2 can be seen in one cluster, respectively. However, the G3 and Y3 of S3 stage were highly inconsistent with the first two stages (S1 and S2). The findings of PCA analysis showed a consistent trend with the color transformation process of the melon stigma. Therefore, we believe that differences in the transcriptional level at S3 determined the developmental specificities and stigma color of the contrasted melon lines (MR-1 and M4-7).

Identification of DEGs during Stigma Development Stages
To identify the DEGs of the melon stigma color, we compared the FPKM values of each gene in both lines, MR-1 and M4-7, during different developmental stages. DEGs were retained with a fold change of >2, and false discovery rate (FDR) correction was set at p < 0.01. The comparative transcriptome data were analyzed at different developmental stages, where both lines showed a consistent trend during the first two stages (S1 and S2) and a significant increase in DEGs was observed at S3 ( Figure S1). We also compared the transcriptome data of the same stage for both varieties and observed that the number of

Identification of DEGs during Stigma Development Stages
To identify the DEGs of the melon stigma color, we compared the FPKM values of each gene in both lines, MR-1 and M4-7, during different developmental stages. DEGs were retained with a fold change of >2, and false discovery rate (FDR) correction was set at p < 0.01. The comparative transcriptome data were analyzed at different developmental stages, where both lines showed a consistent trend during the first two stages (S1 and S2) and a significant increase in DEGs was observed at S3 ( Figure S1). We also compared the transcriptome data of the same stage for both varieties and observed that the number of DEGs was similar at S1 and S2 stages but doubled at S3 stage (G1 vs. Y1: 2523; G2 vs. Y2: 2614; G3 vs. Y3: 5312) ( Figure S2). However, the total DEGs of S3 stage seemed more likely to be associated with melon stigma color.
The regulatory network involved in plant development is very complex and involves various aspects of biological processes. To further illustrate the key biological functions of the transcriptome at different stages of MR-1 and M4-7, we conducted gene ontology (GO) enrichment analysis of all DEGs during all stages. It was noticed that the identified DEGs of both melon lines were mainly enriched in many of the same terms, e.g., plant organ morphogenesis (GO term: 1905392), response to auxin (GO term: 0009733), and regulation of hormone levels (GO term: 0010817) (Table 1). Overall, the biological processes involved in melon stigma development were consistent, despite the color differences. It is noteworthy that chlorophyll and the related biological pathways (photosynthesis and chloroplast localization) were enriched in S2 and S3 only. Therefore, we visualized the expression profiles of DEGs between MR-1 and M4-7 at S3 using MapMan ( Figure 5A). The highest number of DEGs in the light reaction pathway and further analysis revealed that the transcriptional activity of the MR-1 line was higher than that of the M4-7 line, in terms of the Calvin cycle and the light reaction ( Figure 5B). enrichment analysis of all DEGs during all stages. It was noticed that the identified DEGs of both melon lines were mainly enriched in many of the same terms, e.g., plant organ morphogenesis (GO term: 1905392), response to auxin (GO term: 0009733), and regulation of hormone levels (GO term: 0010817) (Table 1). Overall, the biological processes involved in melon stigma development were consistent, despite the color differences. It is noteworthy that chlorophyll and the related biological pathways (photosynthesis and chloroplast localization) were enriched in S2 and S3 only. Therefore, we visualized the expression profiles of DEGs between MR-1 and M4-7 at S3 using MapMan ( Figure 5A). The highest number of DEGs in the light reaction pathway and further analysis revealed that the transcriptional activity of the MR-1 line was higher than that of the M4-7 line, in terms of the Calvin cycle and the light reaction ( Figure 5B).

Figure 5.
MapMan overview of the metabolic pathways (A) and photosynthesis (B) with the differential expression profile in MR-1 as compared to M4-7 during S3. Differentially expressed genes (DEGs) between MR-1 and M4-7 at S3 were loaded into MapMan to generate the overview. On the log2 scale, dark blue and dark red represent higher and lower expressions in MR-1 compared to M4-7, respectively.

Identification of WGCNA Modules Associated with Photosynthetic Reactions
The co-expression networks were constructed on the basis of pairwise correlations of gene expression across all samples. We identified 31 modules by weighted gene co-expression network analysis (WGCNA) as shown in the constructed dendrogram ( Figure 6A). The analysis of module-trait relationships revealed a "black" network module (655 out of 5302 genes), which was highly correlated with a green stigma ( Figure 6B). Gene ontology (GO) enrichment of genes in the black module significantly depicted photosynthesis, light harvesting, a generation of precursor metabolites, and a chlorophyll biosynthetic process (Figure 7).

Figure 5.
MapMan overview of the metabolic pathways (A) and photosynthesis (B) with the differential expression profile in MR-1 as compared to M4-7 during S3. Differentially expressed genes (DEGs) between MR-1 and M4-7 at S3 were loaded into MapMan to generate the overview. On the log2 scale, dark blue and dark red represent higher and lower expressions in MR-1 compared to M4-7, respectively.

Identification of WGCNA Modules Associated with Photosynthetic Reactions
The co-expression networks were constructed on the basis of pairwise correlations of gene expression across all samples. We identified 31 modules by weighted gene co-expression network analysis (WGCNA) as shown in the constructed dendrogram ( Figure 6A). The analysis of module-trait relationships revealed a "black" network module (655 out of 5302 genes), which was highly correlated with a green stigma ( Figure 6B). Gene ontology (GO) enrichment of genes in the black module significantly depicted photosynthesis, light harvesting, a generation of precursor metabolites, and a chlorophyll biosynthetic process (Figure 7).   WGCNA can also be employed to construct the gene networks, in which each node represents a hub gene, and connecting lines between hub genes represent the co-expression of correlations. Here, we only show the top ten hub genes for 'radiality' values (Figure 8). Mainly, four of these hub genes (MELO03C011657, MELO03C014648, MELO03C021788, and MELO03C006053) are annotated as chlorophyll-related (Table 2).

Gene_ID
Gene Annotation MELO3C014477 Cyclic nucleotide-gated ion channel 1-like isoform X3 MELO3C006751 DNA-binding protein MELO3C012015 WAT1-related protein MELO3C011657 Tetrapyrrole-binding protein, chloroplastic MELO3C023604 Glutathione transport system permease gsiD MELO3C014648 Plastid lipid-associated protein MELO3C021788 Photosystem I reaction center subunit III MELO3C006053 Protein CURVATURE THYLAKOID 1A, chloroplastic MELO3C009659 WD-repeat protein, putative WGCNA can also be employed to construct the gene networks, in which each node represents a hub gene, and connecting lines between hub genes represent the co-expression of correlations. Here, we only show the top ten hub genes for 'radiality' values ( Figure 8). Mainly, four of these hub genes (MELO03C011657, MELO03C014648, MELO03C021788, and MELO03C006053) are annotated as chlorophyll-related (Table 2). Table 2. The gene annotation of 10 hub genes.

Gene_ID
Gene Annotation

MELO3C014477
Cyclic nucleotide-gated ion channel 1-like isoform X3 MELO3C006751 DNA-binding protein MELO3C012015 WAT1-related protein MELO3C011657 Tetrapyrrole-binding protein, chloroplastic MELO3C023604 Glutathione transport system permease gsiD MELO3C014648 Plastid lipid-associated protein MELO3C021788 Photosystem I reaction center subunit III MELO3C006053 Protein CURVATURE THYLAKOID 1A, chloroplastic MELO3C009659 WD-repeat protein, putative MELO3C015140 Cytochrome P450 Therefore, we further determined the expression of five key genes of the chlorophyll synthesis pathway in MR-1 and M4-7 at different stages of stigma development by performing qRT-PCR. The results are consistent with our prediction that the expression of these genes gradually increases during the stigma development of MR-1 and decreases in M4-7 at S2 (Figure 9). In addition, the expression trends in different stages were similar to RNA-seq, indicating the accuracy of the RNA data. Therefore, we further determined the expression of five key genes of the chlorophyll synthesis pathway in MR-1 and M4-7 at different stages of stigma development by performing qRT-PCR. The results are consistent with our prediction that the expression of these genes gradually increases during the stigma development of MR-1 and decreases in M4-7 at S2 (Figure 9). In addition, the expression trends in different stages were similar to RNA-seq, indicating the accuracy of the RNA data.

Differences in Chloroplast Development of the Contrasted Melon Lines (MR-1 and M4-7)
Color phenotypes are often considered a communication signal in plants whose genetic diversity is driven by the direct natural selection of biochemical (pigmentation) and

Differences in Chloroplast Development of the Contrasted Melon Lines (MR-1 and M4-7)
Color phenotypes are often considered a communication signal in plants whose genetic diversity is driven by the direct natural selection of biochemical (pigmentation) and structural (cell shape) characteristics, resulting in color. At the physiological level, pigment content and chloroplast development were associated with coloration in higher plants [10]. The central accumulation of photosynthesis takes place in the chloroplast, which is composed of a chloroplast and a thylakoid membrane matrix (platform for photosynthesis's light reactions) [26].
Herein, transmission electron microscopic observations revealed a large difference in the process of chloroplast development between different stages of contrasted melon lines (MR-1 and M4-7). Some studies have reported that differences in chloroplast development affect chlorophyll content and thus affect leaf coloration [27][28][29][30]. It was reported that a tea cultivar with a yellow leaf exhibits a lack of chlorophyll content due to abnormal chloroplast structures under high light intensity; however, the DEGs were supposed to be predominantly involved in chloroplast development and photosynthetic pigment synthesis [31]. RNA-seq analysis revealed that the chloroplast ultrastructure and multiple pigment biosynthesis have a significant impact on leaf color formation in mutants of Anthurium andraeanum 'Sonate' [10], and nitrogen deficiency also triggers irregular chloroplast arrangement in rice leaves due to the gradual development of unbalanced thylakoids [32]. In our study, the number of cystoids in MR-1 chloroplasts, compared to M4-7 chloroplasts, was found to be excessive from S1 to S3; however, the cystoid membrane was clearer, and granules were more uniform. Carotenoid accumulation in the cystoid membrane was higher in MR-1 than in M4-7 [33,34]. Carotenoids are potent antioxidants that act as an effective barrier to high-energy blue light [35,36], so the higher carotenoid content of the MR-1 stigma may help its female flowers to be able to withstand stronger light. Similarly, in the rice mutant (CH1), abnormal chloroplast development along with thinner grana was found to be associated with a reduction in pigmentation, an elevated Chl a/b, and reduced photosynthetic capacity [37], but further experiments are needed to demonstrate this assumption.

A Complex Regulatory Network Exists for Melon Stigma Color Development
Transcriptome analysis has evidently shown changes in gene expression, which helped to gain insight into biological processes and revealed spatiotemporal expression patterns across the genome [30]. In this study, the transcriptomes of two melon lines (MR-1 with a green stigma and M4-7 with a yellow stigma) were characterized to identify the candidate DEGs involved in three stages (S1, S2, and S3). It was found that the number of DEGs of S1 and S2 was lower than that of S3 ( Figures S1 and S2), and PCA analysis similarly showed a consistent trend with the color transformation process of contrasted melon stigma ( Figure 4). Therefore, it is evidently assumed that the difference in the transcriptional level at S3 may determine the developmental specificities and stigma color of the MR-1 and M4-7 lines.
It has been stated that many metabolites with complex regulatory networks play simultaneous roles in the process of plant growth and development [38]. For instance, phenylpropanoid and lignin pathways are significantly involved in the transition from pit-hardening to the cell enlargement stage and subsequently to the ripening stage in peaches [39]. In this study, GO enrichment analysis presented several important biological processes, cellular and metabolic processes, DNA binding, oxidoreductase activity, transcription regulation, and nucleus and cell wall activities associated with melon stigma development during S1 and S2 stages. These results suggest that the biological processes involved in the development of melon stigma were irrespective of the stigma color. Similar phenomena were observed in peach and olive flesh, where carotenoids accumulated without any noticeable changes in the early stages but with a significant accumulation during the ripening period [40,41].

An Important Biological Process Played a Key Role in Chlorophyll Synthesis Variation and the Formation of Different Colors in Melon Stigmas
The WGCNA of 5302 DEGs enabled us to identify the "black" module that includes 655 genes, which were highly correlated with stigma color. According to the GO enrichment analysis of the "black" module, the critical stage (S3) of color presentation showed a maximum number of expressed genes related to Photosystem I, Photosystem II, chloroplast, photosynthesis, and light reactions. This illustrated the key role of chlorophyll in the color presentation of melon stigmas. Thus, the WGCNA analysis package [42,43] has been extensively used for similar analyses in many studies [44][45][46][47], and it was also meaningful in the biological sense of our study.
Chlorophyll and carotenoids are the main pigments contributing to significant color formation [48]. It was reported that different compositions and different concentrations of chlorophyll and carotenoids result in a variety of colors, e.g., red, orange, yellow, green, light green, and white, in the flesh of Cucurbit fruits [49,50]. In this study, we further examined the endogenous chlorophyll and carotenoid content of two contrasted melon lines and observed the metabolic characteristics of the varying stigma coloration. The endogenous chlorophyll content of the green stigma of MR-1 was much higher compared to that of the yellow stigma of M4-7, and this result is consistent with the characteristics of the transcriptome data. In addition, the carotenoid content of the green stigma of MR-1 was higher than that of M4-7, but the overall carotenoid content of both materials was small, and we speculated a limited role in stigma coloration. There are other studies with similar results, which indicate that chlorophyll-deficient mutant tomatoes turned yellow due to abnormal chloroplast development [13]; however, yellow naringenin chalcone and carotenoids clearly exist in the peels of dark-green fruits [51].
In this study, four of the 10 candidate genes were obtained by WGCNA analysis, and they were directly related to the chlorophyll content. The differential roles of these genes have been reported in a few published studies. The MELO3C006053 gene encodes curvature thylakoid 1 (CURT1), which is a major contributor for the internal shaping of the chloroplast [52]. In Arabidopsis thaliana, CURT1 proteins are present in etioplasts, and a loss of CURT1 results in the looser packing of the PLB paracrystalline lattice, which leads to a faster disassembly of PLBs, reduces chlorophyll synthesis, and triggers the accumulation of LHCII [53]. The MELO3C021788 gene was annotated as the Photosystem I (PSI) reaction center, which is a large super-complex protein that catalyzes the light-dependent oxidation of plastocyanin and reduces the ferredoxin; however, its biogenesis process requires a specific regulatory and quality control network [54]. The MELO3C014648 gene encodes a plastid lipid-associated protein (PAP), which belongs to the fibrillin family related to chromoplast fibrils, thylakoids, and photosynthetic antenna complex [55]. This PAP protein has been reported to perform an important role in the development of chloroplasts in watermelon and is similarly involved in the chlorophyll metabolic pathway, regulating the pale green flesh color due to the maximum endogenous chlorophyll accumulation [56,57]. The MELO3C011657 encodes the tetrapyrrole-binding protein, which has been significantly demonstrated in the Arabidopsis papp5 mutant, and it has been shown that the tetrapyrrolebinding protein is mainly associated with PAPP5 in complex components involved in the downstream of the plastid signaling pathway of tetrapyrrole Mg-ProtoIX/Mg-ProtoIX-ME [58]. The complexity of these components can act as a negative regulator of PhANG expression during chloroplast development, which ultimately affects the chloroplast morphology associated with plant flower color morphology [59]. Most importantly, we believe that our current research findings demonstrate comprehensive bio-information about the candidate pathways and associated genes involved in the determination of melon stigma color.

Plant Materials and Sampling Stages
The seeds of contrasted melon lines, MR-1 with a green stigma and M4-7 with a yellow stigma, were obtained from the Laboratory of Molecular Genetics and Breeding in Watermelon and Melon, at Northeast Agricultural University, Harbin, China. Thirty plants of each melon line were cultivated in a plastic greenhouse at Xiangyang Agricultural Experiment Station. The plant developmental stages were regularly checked, and stigmas of both melon lines (MR-1 and M4-7) were collected during the three key developmental stages, including the first day (the S1 stage, where stigmas of both melon lines were green), the seventh day (the S2 stage, where the stigmas of MR-1 were green, but the stigmas of M4-7 ranged from green to partially yellow), and the twelfth day (the S3 stage, where the stigmas of MR-1 were still green, but the stigmas of M4-7 were completely yellow), respectively ( Figure 1). The required amount of samples was collected and used for RNA sequencing analysis and other subsequent experimentation.

Metabolite Detection and Transmission Electron Microscopy
The collected stigmas of both melon lines were subjected for total carotenoid profiling using the Plant Carotenoid ELISA Kit (Guiechem) and a microplate reader. A sufficient quantity of stigmas (0.1 g) was taken from the S3 stage (with completely differentiated stigma colors) of both melon lines, and samples were then immersed in 0.9 mL of lysate and homogenized with a bullet blender tissue homogenizer for 5 min. The final solution was later centrifuged at 5000× rpm for 5 min, and the upper supernatant was removed for subsequent testing. The absorption rate of the final extracts was measured at a 450 nm wavelength using a microplate reader, and carotenoids were determined.
For endogenous chlorophyll quantification, 0.1 g of stigmas were cut from both melon lines, immersed in 1 mL of acetone (80%, v/v), and homogenized using a bullet blender tissue homogenizer for 5 min. The absorption of chlorophyll extracts was measured at 663.6 and 646.6 nm wavelengths using a Hitachi U2800 microplate reader, and chlorophyll was determined according to a previously reported method [60]. Similarly, the collected stigmas of three key stages were subsequently checked to observe the major differentiation in chloroplast ultrastructure with a transmission electron microscope, according to a previously reported method [61].

RNA Sequencing and Transcriptomic Analysis
The stigmas of both melon lines were collected at all three stages. Samples were quickly frozen in cryogenic liquid nitrogen and stored at an ultra-low temperature of −80 • C until further experimentation. The RNA-seq libraries were produced with 4 µg of total isolated RNA following the instructions of the Illumina TruSeq RNA sample preparation kit (FC-122-1001). After amplification and purification, the libraries with an average size of 300 bp were checked on a 2% low-range ultra-agarose gel (BIO-RAD). RNA quality (RIN > 8) and library size were assayed on a 2100 Bioanalyzer (Agilent Technologies). The RNA-seq library construction, Illumina sequencing, and read mapping were performed at Novogene (Beijing, China, https://www.novogene.com/, accessed on 18 May 2022).
The quality of the raw sequenced reads from all samples was checked using FastQC (v0.11.2). All clean-end reads obtained from each sample were then mapped onto the melon reference genome (DHL92, v3.6.1) using the default parameter setting of TopHat (v2.0.11). The unique mapped reads of each specific transcript were counted using HTSeq (v0.6.1). The raw RNA-seq data were uploaded to the Sequence Read Archive (SRA) database (PRJNA839169, 18 independent libraries including biological replicates in this bioproject) for access by the scientific community (https://www.ncbi.nlm.nih.gov/bioproject/?term= PRJNA839169, accessed on 18 May 2022).
The gene abundances were calculated and normalized to RPKM (reads per kb per million reads). The edgeR package (http://www.rproject.org/, accessed on 18 May 2022) was used to identify the differentially expressed genes (DEGs) across the contrasted groups with fold changes ≥2 and a false discovery rate (FDR) <0.05. Pearson's correlation coefficients between the independent biological replicates for each sample were calculated and demonstrated using the R package (heatmap). DEGs were then subjected to an enrichment analysis of GO (gene ontology) pathways, gene numbers were calculated for every term, and significantly enriched GO terms in DEGs compared to the genome background were defined by a hypergeometric test. The calculated p-value was subject to FDR correction, taking an FDR of ≤0.05 as a threshold.
The co-expression networks were constructed using the weighted gene co-expression network analysis (WGCNA, v1.47) package in the R language programming tool (v4.2.0) [46]. After filtering the genes, gene expression values were imported into the WGCNA to construct co-expression modules. Expression correlation coefficients of the remaining genes were then calculated to search a suitable soft threshold for building gene networks using a scale-free topology model [62]. To identify biologically significant modules, module engines were then used to calculate the correlation coefficients. The intramodular connectivity (function soft connectivity) of each gene was calculated, and the top 1 or 5% of genes with the highest connectivity tended to be hub genes. The networks were visualized using Cytoscape_3.3.0. For identification of the genes in each module, GO enrichment pathway analyses were conducted to analyze the biological functions of the modules, and a q-value of <0.05 was used as the threshold after correction.
MapMan figures were obtained by running the Mercator tool (http://mapman.gabipd. org/web/guest/mercator, accessed on 18 May 2022) with default parameters to assign MapMan bins to melon transcripts [63]. The Log2 fold changes as obtained from the DESeq output were used as the MapMan input to represent expression changes. The Bioconductor package Pathview version 1.6.0 [64] and Pathview parameters were set as default ones, and the limit parameters were set as follows: limit = list (gene = 5, cpd = 1). As per Pathview default settings, log2-fc values for boxes representing more than one gene were summed.
We further explored the key hub genes that had important biological functions in each module. First, the original WGCNA output of this study was divided into different modules, which are described in Section 2.4, and modules with significantly high correlation values and where p < 0.05 were identified. Subsequently, we defined key hub genes as highly connected and differentially expressed in each significant "black" network module. We then analyzed the hub genes for each heat-map module and visualized them in Cytoscape.

qRT-PCR Analysis
First strand complementary DNA (cDNA) was synthesized from 100 ng of total isolated RNA using the TRUEscript 1st Stand cDNA Synthesis Kit (Aidlab, Beijing, China) [65]. The gene primers were designed using Primer Premier 6.0 software, and relative expression levels of 5 key genes of the chlorophyll synthesis pathway were analyzed, and transcriptome data between MR-1 and M4-7 was validated using qRT-PCR. The assays were performed in 10 uL of reaction mixture of the SYBR Green I Master Mix, including 20 ng of cDNA and 300 nM of each gene primer. Three biological replicates of each stigma tissue sample and at least three technical replicates of each biological replicate were used for the subsequent expression analysis of five key genes. The transcript level of each key gene was normalized for each sample using the most suitable internal control gene (Actin, MELO3C023264), and fold changes in relative gene expression was calculated using the 2 −∆∆CT (−delta delta CT) method.

Conclusions
In this study, comparative transcriptome analysis of the green and yellow stigmas of two different melon lines (MR-1 and M4-7) was performed. The primary differentiation in the chloroplast development of the contrasted stigmas provide morphological support for endogenous pigment content variations. The WGCNA analysis identified a "black" module that was found to be highly associated with stigma color and enriched in photosynthesis, light reactions, light harvesting, the chlorophyll metabolic process, and the chlorophyll biosynthetic process. Meanwhile, the relative expression of five key genes involved in the chlorophyll synthesis pathway was consistently higher in MR-1 (green stigma) than in M4-7 (yellow stigma). These findings evidently suggest that the chlorophyll and its related biological process played an important role in the formation of different colors in melon stigmas. Thus, we believe that our research findings will deliver valuable insights on the in-depth molecular mechanisms of stigma color formation in melons.

Acknowledgments:
We are grateful for the valuable assistance provided by research group members of our laboratory.

Conflicts of Interest:
The authors have declared no competing interest. The funders had no role in the design of the study, in the collection, analysis, or interpretation of data, in the writing of the manuscript, or in the decision to publish the results.

FPKM
Fragments per kilobase of exon model per million mapped fragments FDR False discovery rate DEGs Differentially expressed genes GO Gene ontology WGCNA Co-expression network analysis PSI Photosystem I CUTRT1 Curvature thylakoid 1