Two watermelon parent lines “LSW-177 (female as P₁, short-sized FP with less weight, thick diameter, and big cell structure in pith) and ZXG-1555 (male as P₂, lengthy FP with more weight, thin diameter, and small cell structure in pith)” were chosen as experimental materials on the basis of their primary phenotypes. The cytological and morphological differences of FP traits among both parent lines and their resultant genetic populations are shown respectively (Figure 1, Supplementary Figure 1).

A field experiment was conducted in the big plastic greenhouse (600 × 300 m²) at the XiangYang Experimental Base of Northeast Agricultural University (NEAU), Harbin (44°04′ N, E125°42′), China. In 2020, whole genome re-sequencing-based bio-informatic analysis were performed between two contrasted parents (P₁ and P₂) and novel SNP-derived CAPS markers were developed. These parents were also crossed to produce the F₁ progeny, and this F₁ was further self-pollinated to develop the respective mapping populations of 100-F₂ derived F₃ families.

In 2021 year, the seedlings of all genotypes were raised and 100-F₃ families (5 plants per family) were cultivated in a complete randomized design (CRD), along with P₁ (15 plants), P₂ (15 plants), and F₁ progeny (15 plants), respectively. The plants were cultivated at 65 cm apart with rows that were 45 cm apart, and standard cultivation practices were adopted to promote the healthy growth of plants and decrease the chances of major disease attacks, etc.

FPs of F₃ families were collected during the harvesting of a commercial fruit crop of watermelon at days after pollination (DAP) and maturity signs. The weight of FP (FPW) was freshly calculated using a digital weighing scale, and values were recorded in grams (g). Then, FPs were manually straightened and pedicel length (FPL) was measured. FP diameter (FPD) was checked with a phenotypic precision of thickness from the upper, middle, and lower positions. FPL and FPD were measured with a digital vernier caliper and data values were recorded in mm.

The freshly harvested FPs were sampled, dissected in longitudinal form, and 1 cm paraffin-embedded sections were prepared. All samples were then fixed in an equal fixative solution of 3.5% formaldehyde, 5% glacial acetic acid, and 50% ethanol for complete embedding. A marked series of ethanol (70, 85, 95, and 100%) was used for the ultra-dehydration of fixed samples followed by the xylene and ethanol ratio. A total of 10-μm microtomes were prepared, stained by using Toluidine Blue, and examined at 40 × magnification (0.1 mm scale) of optical microscope as mentioned (Yang et al., 2021). The total cell numbers of each pedicel sample were photographed and visually checked by using the digital scale bar of image processing software “ImageJ” (https://imagej.nih.gov/ij/).

The leaf material (0.2 g) was taken from two-week-old seedlings of all cultivated genotypes, snap-frozen with cryogenic liquid nitrogen, and preserved at an ultra-high-temperature of −80°C. A good quality gDNA was extracted using the cetyl trimethyl ammonium bromide (CTAB) protocol with slight modifications (Allen et al., 2006). The quantified gDNA was fragmented with Bioruptor (Thermo Fisher Scientific, USA) and 200–300 bp fragments were pooled-up. A sequencing library of each parent line gDNA was then prepared as follows: (1) the gDNA fragments were subjected to end-repairing and A-tailing; (2) the resulting fragments were ligated with bubble adapters that contained a sequential barcode, and then amplified with polymerase chain reaction (PCR). Next-generation sequencing (NGS, High-throughput Illumina Sequencing^TM 2500) based whole genome re-sequencing (>25× genome depth coverage) was performed at Beijing Genomics Institutes (BGI) company, Shenzhen, PR China.

The total re-sequencing data of each watermelon line were quality-checked using the pre-processing FASTX-Toolkit v0.0.13. We first filtered the raw/low-quality adapter sequences, then went through a series of data processing for quality trimming, and obtained the valid data of clean read bases by using the filter parameters of SOAPnuke software (developed by BGI) as follows: (1) the reads that match > 25% of adapter sequences were deleted for filtered adapter, (2) the reads that match more than > 50% of bases having a quality value of <20 were removed, and (3) then reads with > 3% N were removed. The resultant clean end bases of paired-end sequencing of both parent lines were aligned to the de novo assembled reference genome (v2) of watermelon (97103) using the Burrows-Wheeler Aligner (BWA)-MEM algorithm (https://sourceforge.net/projects/bio-bwa/files/) (Li and Durbin, 2009). The Binary Alignment/Map (BAM) files were similarly used for fixing the filtered alignment of re-sequenced bases, and marked duplicate bases were discarded using the Picard tools (v1.106), and the finally sorted BAM sequences were used for subsequential analysis.

Whole-genome SNP variants were called by using the Soaps v.1.03 (Tian et al., 2009), the Sequence Alignment/Map format of SAMtools software package v1.14 (http://www.htslib.org/download/) (Li et al., 2009), and Genome Analysis Toolkit (GATK) v4.2.0.0 (https://www.broadinstitute.org/gatk/) (McKenna et al., 2010). This program allows the local alignment and calibration to generate the SNP variants per sample in VarScan2 file format. The output file format consisted of indexed information related to the positions of variants along with coverage of quality mapping and in-depth read statistics (Ruffalo et al., 2011).

All the filtered SNP variant call format (VCF) were annotated based on their physical positions, and categorized according to their genetic variation using SnpEff v4.3 (Cingolani et al., 2012). The binary (.bin) database of reference annotated watermelon genome files (.gff3) and genomic sequences were downloaded and aligned for annotation of genetic variant effects by their regions (intergenic and exonic) and impact (high, low, moderate, and modifier), and similarly classified by their functional classes (missense, non-sense, and silent), respectively. The web-based generated HTML and VCF files were used for identifying the genetic locations of introns, exons, start-stop codons, upstream-downstream regions, splice regions, and 5' to 3' end UTR regions.

Whole-genome SNP-CAPS loci pairs were identified using the following steps as mentioned in our earlier study (Amanullah et al., 2020). Firstly, putative SNPs sequences were detected from the mapped re-sequencing reads by setting the default protocol of BWA and SAM tools, respectively. Then, flanking sequences with >20 quality mapping scores and 500 base pairs (bp) were pair-wise filtered at each SNP locus. The frequency of finally retrieved SNP-CAPS loci pairs were analyzed by the self-compiled perl script and SNP2CAPS tool, with an input parameter setting and required cutting sites of different endonucleases (Thiel et al., 2004), and the further missing CAPS loci were also retrieved through the same procedure by using SnapGene viewer (v6.0). Finally, a total of 40–45 suitable CAPS loci pairs were used for primer designing and oligo synthesizing.

The physical positions of whole-genome SNP-CAPS loci pairs were cross-checked and PCR primers (forward and reverse sequences) were developed using the latest primer designing tool “Primer Premier (v6.20)” as mentioned (Amanullah et al., 2020). The default preferred parameters were as followed; optimal length of 18–22 bp, primary melting temperature (T_m) of 50–60°C and 3°C gradients in secondary annealing (T_a) temperature, GC contents of 40–60%, and amplification length of >500 bp. The exported primers were then marked with a unique ID and abbreviated as “chromosome number, respective name of cutting site endonuclease, and a number of markers on each chromosome”, respectively.

The oligo-synthesized primers were preliminary screened using the optimized polymerase chain reaction (PCR) in a 20 μL mixture as mentioned (Amanullah et al., 2021). Then, the exact polymorphism was cleaved by performing the enzyme digestion reaction in a total of 15 μL mixture (5 PCR product and 10 ul of digestion mixture), incubated at the respective endonuclease temperature (37°C) for a total 3-h duration. The finally digested products were visually tested on 0.8% agarose gel electrophoresis at 150-volts for 35 min, gel images were captured, and polymorphism level was validated among three testing genotypes (P₁, P₂, and F₁), respectively. The finally screened polymorphic markers were effectively genotyped with newly developed F₂ population, and homozygous and heterozygous allelic bands were coded for mapping of respective QTLs.

A genetic linkage map was developed and putative QTLs were detected by using the default parameter settings of the integrated IciMapping (v4.2) tool as followed (Meng et al., 2015). The coded genotypic data of allelic bands of 100-F₂ mapping population was grouped and anchored over the whole genome chromosomes according to their exact physical positions and maximum likelihood means, using the “rippling” method. The default Chi-squared method was used to determine the segregation ratio of genotypic markers and distorted markers were eliminated at p-value of > 0.001. The default kosambi mapping function and recombination frequency value (0.35 rf) were used to estimate the genetic distances, intervals, and positions in centimorgans (cM), respectively.

Quantitative trait loci were detected by using the inclusive composite interval mapping with additive effects (ICIM-ADD) approach using the sliding parameter (1 cm) and permutation testing (1000×) at p-value of 0.05 (Amanullah et al., 2021). The significant QTLs were defined as above the threshold level of the default set cut-off logarithm of odd (LOD) score (2.50) and genome-wide type I error at α = 0.05. QTLs with >10% phenotypic variance explained (PVE) and high LODs were named major-effect QTLs, and QTLs with <10% PVE and low LODs were named minor-effect QTLs. Then, all the mapped QTLs were abbreviated as follows; trait name, chromosome number, and QTL number, respectively.

The detected physical positions of flanking CAPS marker sequences were pair-wise aligned over the genomic sequences of re-sequencing data of both parent materials and de novo reference genome database of watermelon (97103, v2) (http://cucurbitgenomics.org/), using the web-based application of Integrative Genomics Viewer (IGV, v2.12.2). Then, total candidate genes underlying the adjacent genomic regions were predicted according to their functional annotation of GO database and major SNP mutations within gene coding region were identified, respectively.

The collected quantitative dataset of phenotypic traits was computed on a Microsoft Excel worksheet (v2016) and statistically analyzed, e.g., descriptive statistics of means, standard deviations (SDs), and frequency distributions, by using the default quantitative methodology of statistics software suit (IBM-SPSS, v26.0). The Pearson's correlation coefficients of the observed dataset were analyzed and graphically represented by using the R program (v4.1.2) with its RStudio interface and functional packages “ggplot2” and “corrplot”, respectively.

The mapping statistics of whole-genome re-sequencing reads of two contrasted parent lines of watermelon is shown respectively (Supplementary Table S1). A total of 98.93% (416 Mb) of de novo assembled reference genome was covered by re-sequencing reads of LSW-177 (P₁). In sum, 13,347,1440 reads were QC-passed and properly paired, and 271,440 reads were as supplementary reads, from which 13,204,4185 reads were efficiently mapped, 12 million (12,200,6330) reads (about 91.60%) were properly paired (map Q≧5), as followed by other types of sequencing reads.

A total of 98.94% (416.50 Mb) of the reference genome was covered by the re-sequencing reads of ZXG-1555 (P₂). In sum, 69,708,578 were QC-passed reads and paired-end sequencing reads, and 1,831,128 were supplementary reads. About 68,971,776 reads (< than total) were successfully mapped, and 60,649,000 reads (about 89.35%) were properly paired (map Q≧5), as followed by another type of re-sequencing reads. The whole-genome re-sequenced reads of both parent lines suggested a very close relationship with the reference genome (V2) of watermelon (97103). Overall, the re-sequencing of both parent lines exhibited that P₁ parent reads were mapped higher than P₂ parent reads, and these reads were further used for bioinformatics analysis and molecular experimentation.

Whole-genome detected SNP variants were pre-filtered and then characterized respectively (Figure 2, Supplementary Table S2), by using the bioinformatics analysis and aligning the mapped re-sequenced reads of two contrasted parent lines over the chromosomal genetic length of reference genome of watermelon (97103, v2). A total of 363,083,641 bp of genetic length was covered by whole-genome chromosomes (Cla97Chr01 to Cla97Chr11), among which Cla97Chr02 exhibited the maximum genomic coverage of 37.92 Mb (37,915,939 bp) length, Cla97Chr04 exhibited a minimum coverage of 27.11 Mb (27,110,815 bp) length, and the remaining chromosomes exhibited the average length of genomic coverage.

Further, SNP variants density was checked within 1-Mb window-sized chromosomal regions of whole-genome (Cla97Chr01 to Cla97Chr11). A total of 49,89,869 SNP variants with an average rate of 72 variants/kb were observed, and these variants were randomly localized over the whole genome chromosomal length. The majority of SNP variants (532,615) were noticed within the genetic length (35.93 Mb) of Cla97Chr09 with an average of 70 variants/kb. The least number of SNP variants (326,114) were exhibited within the genetic length (29.50 Mb) of Cla97Chr06 with an average rate of 90 variants/kb. Regarding the total average variant rate per kb, a minimum average of 68 variants/kb was observed within the genetic length of Cla97Chr08 and a maximum average of 72 variants/kb was observed in Cla97Chr06.

Whole-genome detected SNP variants were well-characterized into transversion-type and transition-type variants as shown in Figure 3A. A total of 99.96% transversions-type (1,688,505 counts) and 99.97% transitions-type (3,301,319 counts) variants were detected across the whole genome re-sequenced reads, respectively. Among the transversion variants, the maximum DNA-base mutations were noticed for 26.24% of AC (209,682 counts) → CA (233,530 counts), 26.40% of TG (211,563 counts) → GT (234,335), and 30.85% of AT (260,317 counts) → TA (260,656 counts), and minimum substitution mutations were noticed for 16.48% of CG (139,583 counts) → GC (138,884 counts), respectively. Among the transition variants, 50% of AG (756,313 counts) → GA (895,322 counts) substitution mutations were observed, which was lower as compared to the 49.97% of CT (849,919 counts) → TC (754,765 counts). However, the ratio of total transitions-type to total transitions-type (TSs/TVs) was 1.95%.

The annotation of SNP variants was performed to assess the proportions of potential genetic variants by their possible function, impact, and region, along with counts and percentages (%) across the re-sequenced watermelon parent lines (Figure 3B). The highest rate of genetic variants (4,186,928 counts) was identified in 48.50% of intergenic_region, 1,946,227 variants were identified in 22.54% of the upstream_region, 1,696,837 variants were identified in the 19.65% of the downstream_region, and 661,236 variants were identified in 7.6% of the intron_region. However, the lowest variants (131,565 counts) were present in 1.52% of the intron_region, respectively.

The SNP variant effects were organized into 4 types, e.g., “high_effect type”, “moderate_effect type”, “low_effect type”, and “modifier_type”. A total of 98.36% of SNP variants were identified with modifier impact in the intron or exon region, 1,538 variants (0.02%) were with high_effects, 80,132 variants (0.93%) were with low_effects, and 60,034 variants (0.70%) were with low_effects. A total of three SNP variant types were identified according to their functional categories. A total of 45.53% of missense variants (60,467 counts) were identified, which switched one amino acid into chain form by connecting it with others. There was a total of 0.53% non-sense SNP variants (704 counts), which caused a shorter chain due to the early stop codon. However, 53.94% of silent mutations (caused by a single-base substitution mutation) were detected with 71,649 counts, where the encoding function of proteins was never altered. The missense to silence class ratio was about 0.84.

Whole-genome re-sequencing reads (>25× depth coverage) of two-parent lines were aligned over the reference genome (v2) of watermelon (97103), respectively. A total of 1,82,949 CAPS loci pairs were excavated by filtering the existence and non-existence sites of suitable endonucleases in each SNP location (Supplementary Table S1). A minimum of 11,688 CAPS loci pairs were positioned in Cla97Chr08 and a maximum of 20,830 CAPS loci pairs were noticed in Cla97Chr09. Overall, 5-endonucleases (Msp I, Hind II, EcoR I, BsaH I, Pst I, BamH I) sites of major CAPS loci were perfectly matched and these were selected for the design of whole-genome markers according to experimental conceptualization.

In sum, 476 pairs of randomly selected CAPS sequences were used for primer designing and further identification of codominant polymorphisms among the DNA of P₁, P₂, and F₁ progenies. The amplified fragmented PCR products and enzyme digested products were cleaved at 0.8% agarose gel electrophoresis. A total of 221 pairs of CAPS primers exhibited the properly digested products within exact base pair (bp) lanes of polymorphic bands and validated the 46.42% cleaved codominant polymorphisms. The remaining 255 CAPS pairs produced the same band patterns (did not show polymorphism) or failed to show proper enzyme digestion, so we excluded those primer pairs. Finally, the enzyme digested products and polymorphic bands of P₁, P₂, and F₁ genotypes were visually observed in different bp lanes (Figure 4). We noticed that all the F₂ genotyped polymorphic markers were located on each chromosome by their proper physical positions (bp), fragment product lengths, and endonuclease sites (Figure 5).

A genetic linkage map was developed by using the coded genotypic data of 221 pairs of novel codominant CAPS markers within 100-F₂ mapping population developed from two parent materials of watermelon (Figure 5, Table 1). The genotypic percent for allele frequency was observed as AA: 24.7%, BB: 24.4%, and AB: 50.7%, respectively. The orders of codominant markers were randomly distributed across the whole-genome linkage map based on their physical positions (bp) and screened fragment lengths. Out of 221 pairs, a total of 200 CAPS markers (90.49%) exhibited the exact ratio of Mendel's segregation (1:2:1) at the p-value of > 0.05, except for the 21 pairs of CAPS markers (9.05%), which indicated the distorted segregation ratio.

The final genotyped markers were exactly positioned with goodness-of-fit over the whole genome chromosomes (Cla97Chr01 to Cla97Chr11). The developed linkage map covered a total genetic length of 2,630.49 cM with an average of 11.90 cM. An average genetic length ranged from a minimum of 157.65 cM (Cla97Chr03) to a maximum of 323.17 cM (Cla97Chr08), respectively. The marker numbers were also varied and ranged from 14 (Cla97Chr03) to 27 (Cla97Chr04) with 209.33 cM length and 7.75 average; however, most of the chromosomes contained a maximum and a normal number of markers (>20 markers). The proper genetic positions of the markers facilitated the possible development of a quality genetic linkage map and fine QTL mapping.

The collinearity analysis was also performed between whole-genome genotypic CAPS markers information and the developed physical and genetic map (Supplementary Figure 2). The scatter plot of the collinearity test also depicted positive and linear regression (r2 = 0.92 to 0.98) values, because the straight curve illustrated an appropriate, increased, and fitting trend of the marker's physical positions (bp) and genetic positions (cM) over the whole-genome chromosomes. Three chromosomes (Cla97Chr04, Cla97Chr06, and Cla97Chr08) contained a few markers with partial segregation that were somewhat deviated from their best fit line, and an inversion in the collinearity curve also represented the extended genetic gaps between the flanking positions of adjacent markers. We suggest that further re-sequencing data based on high-density markers should overcome this problem.

The quantitative dataset of FP traits was collected from two contrasted parent lines (LSW-177 and ZXG-155), F₁ progeny, and respective F_2:3 mapping population. The statistical data analysis depicted normal frequency distribution, transgressive segregation, and highly significant correlation coefficients between the observed traits as shown respectively (Figure 6, Supplementary Table S3).

We checked the longitudinal sections of paraffin-embedded FP samples at 40× magnification (0.1 mm scale) (Figure 1B). There were significant differences in the centered pith cell numbers of FPs of P₁ (LSW-177), P₂ (ZXG-1555), and F₁ samples, and calculated cell numbers of P₁ were 36, P₂ were 92, and F₁ were 90, respectively. The optical cell size was also different between the observed pedicel samples, implying that the longer pedicel phenotype of P₂ parent and F₁ progeny is mainly responsible for the greater pedicel weight variation in F_2:3 mapping population.

The FPL (mm) values of both parent lines and their F₁ exhibited different means due to the major variability in pedicel lengths (long and short pedicel). The values of P₁, P₂, and F₁ progeny were 61.31 ± 1.08 mm, 119.90 ± 2.20 mm, and 100.75 ± 2.00 mm, respectively. The F₁ pedicel length showed higher mean values (mm) than the P₁ parent values but very close to the mean values of P₂ parent. In the F_2:3 population, the FPL mean value was 16.52 ± 3.91 mm, and the range varied from 9 mm to 52 mm, respectively. Further, a standard distribution frequency and a transgressive segregation ratio were noticed in phenotypic analysis of pedicel traits (Figure 6A-1). The observed phenotypic characteristic showed a complete type of dominance effect, inherited from the FP length trait of P₂ parent.

The FPW (g) values of both parent lines and their F₁ also exhibited distinct means (Figure 6A-2) due to the significant variability in pedicel weight morphology (more and less weight). For P₁, P₂, and F₁ progeny, the averaged values were 2.20 ± 0.01 g, 3.70 ± 0.02 g, and 3.60 ± 0.02 g, separately. The values of FPW of F₁ progeny were also higher than female P₁ parent values and tightly related with the P₂ parent values. The overall mean of FP weight of the F₂ population was 2.11 ± 0.94 (g), and variability ranged from the lowest value of 1 g to the highest value of 5 g, with a uniform frequency distribution. However, the male parent line (P₂) seemed to show a major inheritance effect with more FP weight.

The phenotypic values of FP diameter (FPD) of contrasted watermelon parent lines were also distinctive. For P₁, P₂, and F₁ progeny, the mean values were 6.13 ± 0.06 mm, 5.28 ± 0.05 mm, and 5.15 ± 0.05 mm, respectively. The values of FPD of F₁ progeny depicted the close genetic inheritance of pedicel diameter received by the P₂ parent (with less diameter but a connection with more pedicel weight). In the F₂ mapping population, the mean value of FPD was 5.36 ± 0.93 mm and range were from 3 mm to 8 mm, respectively (Figure 6A-3). A uniform distribution of FP weight frequency and variability was in line with differential QTL effects on different chromosomes. Overall, a transgressive frequency was observed for the collected datasets of F_2:3 mapping population, suggesting the important quantitative nature of FP traits.

The correlation matrix indicated a good connection among the FP-associated traits (Figure 6B). A positive and highly significant correlation coefficients (0.75^***, 0.57^***, 0.22^*) were observed among the FP traits (FPL, FPW, FPD) at both levels of ^***p- value of < 0.01 and ^*p-value of < 0.05, respectively. It was noticed that a close relationship between the entire variability of FP traits was present in the F_2:3 mapping population, which suggested a clear quantitative genetic nature, inherited from both parents and F₁ progeny. The well-connection between observed variables also proposes that elongated FPs with a moderate diameter have more pedicel weight than those of short-length pedicels with less weight.

Quantitative trait loci were mapped by combining the coded genotypic data of codominant markers and the means of quantitative datasets of FPL, FPW, and FPD of 100-F_2:3 mapping individuals. Overall, 6 QTLs (5 major-effect and 1 minor-effect) were mapped. Among them, 1 major-effect QTL of FPL, 1 major-effect QTL of FPW, and 4 QTLs (3 major-effect and 1 minor-effect) of FPD were localized between the flanking genetic regions of codominant CAPS markers, respectively. The details of genetic effects of mapped QTLs, e.g., detected LOD scores, additive (Add) effects, dominance (Dom) effects, phenotypic variance explained (PVE%), flanking markers (left-right), their physical positions (bp), and genetic positions (cM) are shown respectively (Figures 5, 7, Table 2).

For the FPL, a single major-effect QTL (qFPL4.1) was mapped near the bottom end of chromosome (Cla97Chr04). The LOD score of 2.52, Add effect of 13.45, and Dom effect of −8.25 were observed. This QTL explained a total of 10.46% PVE for individual trait effect, and detected positive Add effect significantly described the main genetic contribution of elongated FPs in the derived mapping population. Regarding the genetic location (cM), qFPL4.1 was positioned at 191 cM between the flanking regions of two left-right codominant CAPS markers “4Ms37 and 4Hd38”. These two markers were positioned at 187.50 cM and 195.50 cM, and a short genetic interval of 8 cM was noticed between those genetic regions. The physical position (24392485-25067584 bp) of two flanking markers exhibited a total of 675.10 kb genetic interval for detected QTL (qFPL4.1), respectively.

Interestingly, a single major-effect QTL (qFPW4.1) was found to be co-localized with a major-effect QTL of FPL (qFPL4.1) near the bottom end of chromosome (Cla97Chr04), by sharing a common genetic region. For the qFPW4.1, LOD score of 2.94, Add effect of 0.43, and Dom effect of −0.48 were observed, and this QTL explained a total of 13.60% PVE with individual trait effect, and highly positive Add effect significantly described its major contribution for variation in FPs with more weight in derived mapping population. Regarding the genetic location (cM), qFPW4.1 existed at 191 cM between the flanking sections of two left-right codominant CAPS markers “4Ms37-4Hd38” positioned at 186.50 cM and 195.50 cM, and a short genetic interval of 7.5 cM was noticed between them. The same physical position (24392485-25067584 bp) and genetic interval of 675.10 kb were noticed for both QTLs (qFPL4.1 and qFPW4.1) between the same flanking markers. Further, these two QTLs mutually explained a total of 24.06% PVE and mutually indicated the significant contribution of genetic inheritance of lengthy and weighted FPs in the derived F_2:3 mapping population, which might be mainly transferred by P₂ parent.

A total of 4 QTLs (3 major-effects and 1 minor-effect) of FPD (qFPD2.1, qFPD2.2, qFPD8.1, and qFPD10.1) were mapped on three different chromosomes (Cla97Chr02, Cla97Chr08, and Cla97Chr10), showing an unlinked and diverged multi-genic architecture in the watermelon genome, respectively. A major-effect QTL (qFPD2.2) was identified just below the flanking position of nearly positioned minor-effect QTL (qFPL2.1) at the bottom region of chromosome (Cla97Chr02), and detected with an LOD score of 3.54, Add effect of 0.33, Dom effect of 0.63, and 12.76% of PVE. This QTL (qFPL2.2) was genetically positioned at 315 cM between two flanking left-right codominant CAPS markers “2Ms48-2Ms49” pinpointed at 304.50 cM and 323 cM position, with a genomic interval of 18 cM, respectively. The physical position (36931835-37879344 bp) between two markers exhibited the 947.51 kb genetic distance for qFPD2.2. Another single minor-effect QTL (qFPD2.1) was similarly dentified at the bottom region of chromosome (Cla97Chr02), and detected with an LOD score of 2.84, Add effect of 0.33, Dom effect of 0.38, and 8.64% of PVE. This QTL (qFPL2.1) was positioned at 254 cM between two flanking left-right codominant CAPS markers “2Ms46-2Hd47” pinpointed at 238.50 cM and 280 cM, having a genomic interval of 42 cM, respectively. The physical position (35039739-35988127 bp) of these two markers exhibited 948.39 kb genetic interval for qFPD2.1.

Another major-effect QTL (qFPD8.1) was discovered separately at the end of the chromosomal region of Cla97Chr08, and detected with an LOD score of 3.31, Add effect of −0.35, Dom effect of −0.56, and 10.99% of PVE. This QTL was positioned at 202 cM between left-right flanking codominant CAPS markers “8Ms37-8Ms40” pinpointed at 193.50 cM and 211.50 cM, with genomic interval of 18 cM, respectively. The physical position (25356537-26107916 bp) of these two markers exhibited a total of 751.38 kb genetic distance for qFPD8.1. On Cla97Chr10, a separate major-effect QTL (qFPD10.1) was also identified with a detected LOD score of 3.63, Add effect of.54, Dom effect of −0.11, and 13.13% of PVE. Regarding the genetic position, qFPL10.1 was positioned at 132 cM between two left-right flanking codominant CAPS markers “10Ba27-10Ba28” positioned at 116.50 cM and 145.50 cM, exhibiting a genetic interval of 29 cM between them, respectively. The physical position (22813961-23673195 bp) of these two markers revealed total 859.24 kb genetic distance for qFPD10.1. However, a total of 45.52% of mutually explained PVE was noticed for all mapped QTLs of FPD. One minor-effect and two major-effect QTLs (qFPD2.1, qFPD2.2, and qFPD10.1) demonstrated the positive Add effects, predicted by their main genetic inheritance for increased pedicel diameter; however, another major-effect QTL (qFPD8.1) explained the negative additive effect by showing the inheritance of a pedicel with less diameter and weight.

The pair-wise alignment of physical sequences of delimited flanking QTLs and reference genome database of watermelon (97103, v2) revealed a total of 454 putative genes related with FPL, FPW, and FPD were predicted, and presented with different generic gradients of circle diagram (Figure 8), and detailed functional GO annotation of defined genes and detected SNPs in coding region is also mentioned respectively (Supplementary Tables S4–S8).

For the FPD, a total of 386 genes were identified in the 4 QTL regions, e.g., 116 genes in the 948.39 kb region of qFPD2.1, 125 genes in the 947.51 kb region of qFPD2.2, 98 genes in the 751.38 kb region of qFPD8.1, and 47 genes were identified in the 859.24 kb region of qFPD10.1), respectively. The genes annotation analysis described that most of the genes have known functions but very few of them exhibited unknown functions. All the annotated genes were also categorized in the molecular function, cellular components, and biological processes, respectively. However, some genes exhibited more explicit description, e.g., Cla97C02G047510 (Auxin-responsive family protein), Cla97C02G048250 (ABSCISIC ACID-INSENSITIVE 5-like protein 1), Cla97C02G048310 [Late embryogenesis abundant (LEA) hydroxyproline-rich glycoprotein family], Cla97C02G048480 [Late embryogenesis abundant (LEA) hydroxyproline-rich glycoprotein family], Cla97C02G050460 and Cla97C02G050490 [Late embryogenesis abundant (LEA) hydroxyproline-rich glycoprotein family], Cla97C02G050510 (GDSL esterase/lipase family), Cla97C08G157940 (Cytokinin riboside 5'-monophosphate phosphoribohydrolase), Cla97C08G158260 (Cytokinin dehydrogenase 5), Cla97C10G193930 (Late embryogenesis abundant protein 3L-1), Cla97C10G193940 (Late embryogenesis abundant protein 5), Cla97C10G194030 (Ethylene responsive transcription factor 2b), and Cla97C10G194120 (Long cell-linked locus protein).

For the FPL and width, a common genetic position was shared by two co-QTLs (qFPL4.1 and qFPW4.1) which exhibited a total of 68 genes in 675.10 kb region, respectively. Some genes similarly exhibited more explicit descriptions, e.g., Cla97C04G076840 (germin-like protein subfamily 3 member 2), Cla97C04G077140 and Cla97C04G077150 (IAA-amino acid hydrolase ILR1), and Cla97C04G077450 (Auxin efflux carrier component). Overall, non-synonymous SNPs (changing the codon to one that codes for a different amino acid) were found in the majority of annotated genes; however, a greater number of gene SNPs were noticed in the CDS coding regions of ZXG-1555 (P₂) as compared to LSW-1555 (P₁).