The leaf color yellowing mutant material w-yl was obtained from the National Mid-term Genebank for Watermelon and Melon (Zhengzhou, China), the leaves in the whole growth period were yellow, including cotyledon and fruit. Normal green leaf material ZK was supplied by the Diploid Watermelon Genetics and Breeding Research Group of Zhengzhou Fruit Research Institute (ZZFRI) of Chinese Academy of Agricultural Sciences (CAAS). In this study, the mutant material was crossed with the ZK, and six generations were constructed: P₁ (the yellow parent w-yl), P₂ (the green leaf parent ZK), F₁ (orthogonal), BC₁P₁, BC₁P₂, and F₂. The materials were planted in a greenhouse at the Xinxiang Comprehensive Experimental Base of CAAS, with a row spacing of 1.5 m and a plant spacing of 0.4 m. The phenotype of leaf color was determined by visual observation.

Plant for chlorophyll were planted in an artificial climate chamber and treated with different environmental factors at three true-leaf stage: temperature 35°C/28°C, light intensity 30,000 Lx, namely HTHL(high temperature and high light); temperature 35°C/28°C, light intensity 12,000 Lx, namely as HTNL (high temperature and normal light); temperature 35°C/28°C, light intensity 5,000 Lx, marked as HTLL (high temperature and low light); temperature 28°C/25°C, light intensity 30,000 Lx, marked as NTHL (normal temperature and high light); temperature 28°C/25°C, light intensity 12,000 Lx, marked as NTNL (normal temperature and normal light); temperature 28°C/25°C, light intensity 5,000 Lx, marked as NTLL (normal temperature and low light); temperature 15°C/15°C, light intensity 30,000 Lx, marked as LTHL (low temperature and high light); temperature 15°C/15°C, light intensity 12,000 Lx, marked as LTNL (low temperature and normal light); 15°C/15°C, 5,000 Lx, labeled as LTLL (low temperature and low light). Light cycle was 16h/8h, humidity 80%. Each treatment set three replicates. Chlorophyll content was determined after 8 days of treatment.

Plant for chlorophyll precursors and transcriptome sequencing were grown in a smart greenhouse in ZZFRI of CAAS, the light cycle was 16h/8h, the temperature was 25°C/18°C and the light is natural light. The leaves were sampled after 8 days of treatment.

The third true leaf from five seedlings was sampled and mixed, weighed 0.1 g and put into a 15 mL centrifuge tube respectively, added 10 mL of 96% ethanol, and soaked in dark environment until the leaves turned completely white (Yuan et al., 2017). The absorbance A665, A649 and A470 at 665 nm, 649 nm and 470 nm were determined by UV spectrophotometer (UV-2600I, Shimadzu, Kyoto, Japan). The concentrations of chlorophyll a (chla), chlorophyll b (chlb), total chlorophyll (chla+b) and carotenoids were calculated using 96% ethanol as blank control. The equations are as following:

Chlorophyll content (mg·g^-1 = (C × V)/(W × 1000). C represents chlorophyll content, V represents the total volume of extract (mL), and W represents leaf mass (g).

The contents of main chlorophyll precursor in the process of chlorophyll synthesis were measured, among which δ-aminolevulinic acid (ALA) was determined according to the method of Dei (Dei, 2010) and the molar concentration of ALA was calculated with a molar extinction coefficient of 7.2 × 10⁴ mol^-1·cm^-1at 535 nm. Relative contents of protoporphyrin IX (protoIX), Mg-protoporphyrin IX (Mg-proto IX), and pchlide were determined according to the method of Rebeiz (Rebeizjames et al., 1975) and Lee (Lee et al., 1992). The relative mass molar concentration of Mg-Proto IX is presented as F 440 e x 595 : fluorescence emission intensity at 595 nm under 440 nm excitation light. Proto IX ( F 440 e x 633 ) ) = ( F   440 e x 633   − 0.25   ×   F   440 e x 622   − 0.24   ×   F   440 e x 640 ) / 0.95; P c h l i d e  ( 440 e x 640 ) =   ( F   440 e x 640   − 0.03   ×   F   440 e x 633 ) / 0.99.

Leaves of five plants were selected as a sample from w-yl and ZK, with three biological replicates respectively. Total RNA was extracted using RNeasy Plant Mini Kit (Beijing Tiagen), following the manufacturer’s instructions. Then RNA was reversely transcribed to cDNA, and the cDNA fragments were segmented by PCR. Finally, the double-stranded PCR product is thermally denatured to form single-stranded circular DNA, which is then formatted into a final library. The cDNA library was sequenced by BGISEQ-500 system (BGI-Shenzhen, China) with reads of 100bp in length.

The sequencing data were screened to obtain Clean reads, which were then mapped into the ‘97103’ watermelon genome (http://cucurbitgenomics.org/organism/21) using Bowtie2. Gene expression levels were calculated using FPKM (million fragments per kilobase). Based on KEGG (http://www.genome.jp/kegg/) and GO (http://www.geneontology.org/) database for gene annotation and function assignment. Differentially expressed genes (DEGs) were set as gene fold change ≥2.00 and false discovery rate ≤0.001. Through GO enrichment and KEGG enrichment pathways, the significantly enriched metabolic pathways were screened and compared with the whole genome background. Functional classification of DEGs was performed according to GO and KEGG annotation results and official classification, and FDR ≤ 0.01 was set as significant enrichment.

Total RNA was extracted by plant RNA kit (Huayue Yang Biotechnology Co., LTD.). A total of 1.0 μg of RNA was used for cDNA synthesis using the PrimeScript RT kit and gDNA Eraser (TaKaRa) according to the manufacturer’s protocol. Primers were designed using NCBI online tools (https://www.ncbi.nlm.nih.gov/tools/primer-blast/), and synthesized by Sangon Biotech (Shanghai, China). All the primer sequences were shown in Supplementary Table 1 . Quantitative real-time PCR reaction procedure and system were as described previously (Yuan et al., 2022a). All primers are shown in Supplementary Table S1 . The 2^−ΔΔCt method was used to calculate relative gene expression values (Kenneth and Thomas, 2002).

Leaf DNA of 30 individual plants with yellowed and green extreme phenotypes in F2 population were selected for the construction of two extreme sequencing mixed pools, and parental DNA was used to construct the parental pools for sequencing analysis. The depth of parental sequencing was 20×, and the depth of extremely mixed-pool sequencing was 30×. Sequencing was performed by Biomarker Technologies Co, LTD (Beijing, China) using Illumina HiSeq2000. The sequencing read length was 150 bp.

Raw reads were filtered to remove reads containing adapter, and reads containing >5% N and low-quality reads (the number of bases with quality value Q ≤ 10 accounted for more than 50% of the whole read) were used to obtain clean reads for subsequent analysis. Clean reads were mapped to the ‘97103’ watermelon genome (http://cucurbitgenomics.org/organism/21) using BWA software. Then GATK (4.0.4.0) and SNPeff (4.3) were used to annotate the mutation sites, and single nucleotide polymorphisms (SNPs) and insertion-deletion polymorphisms (InDels) were identified.

The SNP-index algorithm was used to establish the target region to find the significant difference in genotype frequency between the pool, and Δ(SNP-index) was used for statistics. In this project, the DISTANCE method was used to fit the ΔSNP-index, and then the region above the threshold was selected as the region related to the trait according to the association threshold. The stronger association between SNP and trait, the closer Δ(SNP-index) to 1.

Using the cDNA of green leaf ZK as template, specific primers were designed to amplify the CDS regions of Cla97C02G036010, Cla97C02G036030, Cla97C02G036040, Cla97C02G036050 and Cla97C02G036060, and the primers were shown in Supplementary Table 1 . BamHI (GGATCC) restriction sites were added to both ends of the primers and inserted into the cucumber green mottle mosaic virus (CGMMV) gene silencing vector PV190 by homologous recombination to construct virus-induced gene silencing (VIGS) vector. The dual vector was transformed into Agrobacterium tumefaciens GV3101.

Induction and inoculation of A. tumefaciens according to Liu (Liu, 2019) When watermelon seedlings were at cotyledon stage, the induced A. tumefaciens was injected from the back of watermelon cotyledon with 1 mL syringe. The blank control (Blank, B), water control (Water, W), medium control (YT medium, Y), blank vector control (PV190, P) and PDS gene positive control (PDS) were set up respectively. Three biological replicates were set up for each treatment. Two weeks after injection, leaf phenotype was observed, and samples were collected for ultrastructural analysis, chlorophyll content measurement and gene expression analysis.

The above-mentioned leaves with phenotype after A. tumefaciens were used as materials, fixed with 4% glutaraldehyde (configured with pH 7.2 phosphate buffer) overnight at 4°C, rinsed with phosphate buffer three times, fixed with 1% osmium tetroxide for 1 h, rinsed with phosphate buffer three times, dehydrated with 30%, 50%, 70%, 80%, 95%, 100% ethanol and acetone step by step for 5 min, and finally embedded with resin. After sectioning, they were stained with 2% uranyl acetate saturated alcohol and lead citrate for 15 min, and the chloroplast ultrastructure was observed under transmission electron microscope (HT7700, Hitachi, Japan).

All data graphs were analyzed by Office 2016 software. Differences were analyzed by SPSS 18.0 software, and one-way ANOVA was used for statistical analysis, p< 0.05 (n = 3) was considered significant difference.

The leaves of w-yl showed yellow throughout the whole growth period ( Figure 1A ), and the color did not change with environmental changes, such as temperature and light intensity ( Figure 1B ). Under different temperature and light intensity, there were no significant differences in the contents of chla, chlb, chla+b and carotenoids.

In addition, phenotypic data showed that all F₁ plants appeared green leaves, indicating that the yellow mutation was recessive. For F₂ plants, among the 237 progeny in the summer of 2018, 178 plants had green leaves and 59 plants had yellow leaves; among the 993 progeny in spring of 2019, 730 had green leaves and 263 had yellow leaves ( Table 1 ). The χ² test of green and yellow leaves in the two seasons showed that the separation pattern was consistent with the Mendelian separation ratio of 3:1 (χ² > χ² _0.05 = 3.841). Furthermore, for the backcross progeny BC₁P₁ and BC₁P₂, the yellow leaf plants were 29 and 0 respectively, indicating that the yellowing mutation of watermelon leaves conformed to the genetic pattern controlled by a single recessive nuclear gene, and green leaves were dominant to yellowing.

Previous studies had demonstrated that there are significant differences in chlorophyll content and photosynthetic indicators (Ren et al., 2019). To further validate the difference, the chlorophyll precursors, including ALA, protoIX, Mg-ProtoIX and pchlide were analyzed (Figure 2). The results showed that the contents of four indexes detected in the w-yl were significantly lower than those in ZK, which explained the low chlorophyll content to a certain extent.

A total of 6 samples were measured by RNA-seq, including 3 samples for w-yl and 3 samples for ZK, yielding an average of 6.06 Gb of data per sample. The average rate of genome alignment was 89.40%, and the average rate of gene set alignment was 65.81% ( Figure 3A ). For the comparison group w-yl-vs-ZK, a total of 19,261 genes were detected, and there were 18,323 shared genes, including 2,471 DEGs ( Figure 3A ), with 848 up-regulated DEGs and 1893 down-regulated DEGs ( Figure 3B ).

GO enrichment and KEGG pathway enrichment analyses were carried out to better understand the function of DEGs. For GO enrichment, 17 of the top 20 selected GO terms were related to photosynthesis process, including 4 terms related to photosystem, such as photosystem (GO:0009521), photosystem I (GO:0009522), Photosystem II (GO:0009523), light harvesting in photosystem I (GO:0009768); 4 terms involved in photosynthesis, such as photosynthesis (GO:0015979), photosynthetic membrane (GO:0034357), photosynthesis—light reaction (GO:0019684) and photosynthesis—light harvesting (GO:0009765); 6 terms involved in thylakoid, such as thylakoid (GO:0009579), thylakoid membrane (GO:0042651), chloroplast thylakoid (GO:0009534), chloroplast thylakoid membrane (GO:0009535), plastid thylakoid (GO:0031976) and plastid thylakoid membrane (GO:0055035); 3 terms relate to pigments, such as tetrapyrrole binding (GO:0046906), chlorophyll binding (GO:0016168) and protein-chromophore linkage (GO:0018298) ( Figure 3C ; Supplementary Table S2 ). These results showed that w-yl and ZK had significant differences in photosynthesis.

For KEGG pathway enrichment, the two most significant enrichment pathways of the top 20 pathways were photosynthesis—antenna proteins and plant—pathogen interaction ( Figure 3D ; Supplementary Table S3 ). Due to the importance of antenna protein for photosynthesis ( Figure 4A ), we focused on the analysis of antenna protein-related DEGs, and completely screened 16 DEGs that encoded antenna protein ( Figure 4B ). LHCI and LHCII, as important components of photosystem I complex and photosystem II complex, are composed of four and six small components, respectively. For LHCI, the number of DEGs that encoded LHCI Chl a/b binding protein 1 (Lhca1), Lhca2, Lhca3 and Lhca4 was 1, 2, 1 and 2, respectively. For LHCII, the number of DEGs that encoded LHCII Chl a/b binding protein 1 (Lhcb1), Lhcb2, Lhcb3, Lhcb4, Lhcb5 and Lhcb6 was 4, 1, 1, 2, 1 and 1, respectively. The expression levels of all the 16 DEGs in ZK were significantly higher than those of w-yl, and the fold change was between 2.6 and 14.0 ( Figure 4B ).

To verify the accuracy of RNA-seq data, 12 DEGs (Cla97C02G035950, Cla97C02G035960, Cla97C02G035980, Cla97C02G036070, Cla97C02G036090, Cla97C02G036110, Cla97C02G036130, Cla97C02G036140, Cla97C02G036150, Cla97C02G036160, Cla97C02G036190 and Cla97C02G036200) of the 29 genes in the interval were selected to conduct qPCR ( Figure 5 ). The results showed that expression patterns of 12 DEGs were highly consistent with those of genes in RNA-seq data, which demonstrated that the RNA-seq data are reliable.

In order to quickly identify the key candidate genes related to leaf color in the F₂ population, 30 green and 30 yellow leaf progeny were selected and sequenced on the Illumina platform. A total of 51.0 Gb clean bases were generated with an average depth of about 26.5×. Finally, we identified 266,255 SNPs between w-yl and ZK, and 83,373 SNPs between the F₂ pools. According to the SNP-index values of w-yl and ZK, the Δ(SNP-index) value of approximately 7.42 Mb genome region (11,540,000-18,960,000) on chromosome 2 was greater than the threshold ( Figure 6A ). These results indicated that this region might contain the key gene of watermelon leaf yellow traits.

In order to further locate the candidate genes for yellowing leaf, the chromosome region of the variation between w-yl and ZK were analyzed. A total of 12 pairs of InDel molecular markers ( Supplementary Table S4 ) were developed for the candidate region based on 233 F₂ populations. The results verified that these genes were in the range of 11.54 Mb—18.96 Mb on chromosome 2. Subsequently, based on the determination of leaf color phenotype data and individual exchange genotype, 12 recombinant individuals were further screened using 1883 F₂ populations. Finally, it is found that the candidate interval corresponds to the 2.217 Mb region of InD14,179,011—InD16,396,362 ( Figure 6B ). There were 29 genes in this region and annotated them according to the watermelon reference genome ( Figure 3C ; Table 2 ). Notably, compared with ZK, Cla97C02G036010, Cla97C02G036020, Cla97C02G036030, Cla97C02G036040, and Cla97C02G036050 were completely absent in w-yl, and Cla97C02G036060 had partial base deletion, suggesting that they were the key genes determining w-yl leaf color mutation ( Supplementary Figure S1 ).

In addition, the results of agarose gel electrophoresis and qPCR showed that Cla97C02G036010, Cla97C02G036020, Cla97C02G036030, Cla97C02G036040, Cla97C02G036050 and Cla97C02G036060 could not be amplified in w-yl, as well as Cla97C02G036020 also had no target product in ZK ( Figure 7 ). The RNA-seq results also showed the same results ( Figure 6C ). These results further proved the importance of Cla97C02G036010, Cla97C02G036030, Cla97C02G036040, Cla97C02G036050 and Cla97C02G036060 in leaf yellowing.

In order to verify the gene function of the candidate genes, cucumber mosaic virus-mediated VIGS vector was used to perform gene silencing assay on ZK leaves. The results showed that at 16 days after inoculation (DAI), the plants inoculated with water ( Figure 8B ), medium ( Figure 8C ) and blank vector ( Figure 8D ) showed no significant difference in phenotype compared with the blank control ( Figure 8A ), while the positive control plants inoculated with PDS gene showed virus symptoms at DAI16, with severe true leaf pucking and chlorosis ( Figure 8E ). Watermelon plants silencing Cla97C02G036010 ( Figure 8F ) and Cla97C02G036030 ( Figure 8G ) showed symptoms of disease at DAI17, and their true leaves were slightly wrinkled and mottled greenish yellow. Watermelon plants silencing Cla97C02G036040 ( Figure 8H ), Cla97C02G036050 ( Figure 8I ) and Cla97C02G036060 ( Figure 8J ) showed obvious virus symptoms at DAI13, with obvious true leaf wrinkling and large area mottled yellow.

Then the expression levels of the silenced genes were detected, when compared with the control group (B, W, Y, P and PDS), their expression levels were significantly reduced. Among them, the expression of Cla97C02G030640 decreased most sharply, which were 2.9%, 2.9%, 2.8%, 3.2% and 27% of the control group, respectively ( Figure 9 ). Besides, the results of chlorophyll content of leaves with phenotype showed that there was no significant difference in chla, chlb and chla+b content among groups B, W, Y and P, while the chlorophyll contents of silenced PDS group was significantly lower than that of the former four groups. For the five silenced genes, the contents of chla, chlb and chla+b were significantly lower than those of B, W, Y and P control groups. Among them, the contents of chla, chlb and chla+b in silenced genes Cla97C02G030640 and Cla97C02G030660 were the lowest and had not significantly different from those in silenced PDS group ( Figure 10 ).

Furthermore, the ultrastructure of chloroplast was analyzed to analyze the reasons for these phenomenon. As a result, the chloroplast ultrastructure of the silenced gene Cla97C02G030610 ( Figure 11B ), Cla97C02G030630 ( Figure 11C ) and Cla97C02G030650 ( Figure 11D ) did not change compared with the blank control ( Figure 11A ), and all contained normal grana lamella (GL) and plastid globule (PL). However, the ultrastructure of silenced gene Cla97C02G030640 and Cla97C02G03060 was significantly changed, and the chloroplast structure may be damaged (red dotted circle area). For silenced gene Cla97C02G030660 ( Figure 11F ), there was no PL, and PG stratification was not obvious, appearing in a fuzzy state. Especially for the silenced gene Cla97C02G030640 ( Figure 11D ), there was no PL and no obvious PG, speculating that PG was degraded.

Taken together, these results indicated that candidate genes play an important role in causing leaf yellowing, especially Cla97C02G030640.