In this work, we selected two Cucurbita moschata accessions (PI 604506 and PI 381814), previously reported (Sáez et al., 2016) as symptomless or with slight symptoms after whitefly and sap inoculation with ToLCNDV, to study the genetic control of the resistance. PI 604506 is the improved pumpkin cultivar Large Cheese from the United States and PI 381814, an Indian landrace. The Chinese C. moschata accession PI 419083 was used as susceptible control. Seeds of the three accessions were first provided by the United States Department of Agriculture-National Plant Germplasm System (USDA-NPGS) genebank, then fixed by selfing and multiplied by the cucurbits breeding group at the Institute for the Conservation and Breeding of Agricultural Biodiversity (COMAV), and stored at the COMAV genebank.

To generate the viral inoculum source, susceptible zucchini plants were agroinfiltrated with an infectious clone based on the Spanish isolate of ToLCNDV (99% nucleotide identity with the sequences of the A and B viral genomic particles: KF749224 and KF749225 (Juárez et al., 2014), following the procedure described in Sáez et al. (2016).

The tissue of symptomatic leaves from 15 days post-ToLCNDV agroinoculation plants was crushed in a mortar together with inoculation buffer [50 mM potassium phosphate (pH 8.0), 1% polyvinylpyrrolidone 10, 1% polyethylene glycol 6000, 10 mM 2-mercaptoethanol and 1% activated charcoal] in a 1:4 (w/v) proportion (López et al., 2015). The homogenate was used to mechanically inoculate all plants at the stage of one true leaf, dusting on the true leaf and on one cotyledon with Carborundum 600 mesh and scratching with a cotton swab dipped in the blend. Inoculated plants were grown in a climatic chamber, and disease progression was monitored. Symptomless plants 15 days after mechanical inoculation (dpi) were reinoculated to avoid escaping to the infection.

Ten seeds of each C. moschata accession were disinfected and germinated as described by Sáez et al. (2016). Seedlings were transplanted to pots and grown in climatic chamber under controlled conditions (photoperiod of 16 h day at 25°C and 8 h night at 18°C and 70% of relative humidity). Subsequently, plants were moved to a greenhouse and crossed to obtain three F₁ progenies: F₁ PI 419083 × PI 604506, F₁ PI 419083 × PI 381814, and F₁ PI 604506 × PI 381814. Eight plants of each parent and the corresponding hybrids were mechanically inoculated with ToLCNDV as described above and phenotyped according to symptomatology and viral accumulation as described below.

Eight additional plants of the C. moschata parents were cultivated in a greenhouse along with eight plants of the F₁ progenies. To generate segregating populations, F₁ plants were selfed to obtain F₂ progenies and backcrossed to plants of PI 604506, PI 381814, and PI 419083 to generate the BC1_{PI 604506}, BC1_{PI 381814}, and BC1_{PI 419083} populations, respectively. All these segregating populations were screened against ToLCNDV with the same inoculation and phenotyping methodology, using three plants of each C. moschata accession as controls. F₂ and BC1 derived from F₁ PI 419083 × PI 381814 were obtained later because of the influence of the local climate conditions in PI 381814 vegetative growth, causing late-flowering and slow development of fruits. Hence, we studied first the genetic control of the resistance to ToLCNDV in the segregating populations derived from PI 604506, and results were validated in F₂ and BC1 coming from F₁ PI 419083 × PI 381814.

Symptomatology was evaluated in all plants at 15 and 30 dpi using the visual scale described by López et al. (2015). Symptoms score ranged from 0 (absence of symptoms) to 4 (highly severe symptoms), classifying as resistant those plants with symptoms scored 0 or 1 and as susceptible those with symptoms scored from 2 to 4. The goodness-of-fit between the expected and observed segregation ratios resistant/susceptible plants was analyzed by chi-squared (χ²) test (p < 0.05) in the F₂ and BC1-segregating populations.

The relative ToLCNDV accumulation in each plant was determined at 30 dpi by quantitative PCR (qPCR). Total DNA from apical leaves was extracted using the cetyltrimethyl ammonium bromide (CTAB) method (Doyle and Doyle, 1990) and quantified using a NanoDrop 1000 spectrophotometer (Thermo Scientific, Waltham, MA, United States). DNA was diluted with sterile-deionized water to a final concentration of 5 ng μl^–1. Three biological replicates were done for each parental genotype, and all plants of the assay were analyzed in three technical replicates using a LightCycler^® 480 System (Roche). In each qPCR reaction, 15 ng of genomic DNA were used as templates, in a final volume of 15 μl. We used 7.5 μl of 2 × iTaq^TM universal SYBR^® Green Supermix (BIO-RAD) and 1.5 μl (100 nM) of each primer and 1.5 μl of H₂O. Primers ToLCNDVF1 (5′-AATGCCGACT ACACCAAGCAT-3′, positions 1145–1169) and ToLCNDVR1 (5′-GGATCGAGCAGAGAGTGGCG-3′, positions 1399–1418) were used for the amplification of a 273-bp fragment of viral DNA-A. The single-copy gene CpACS2 was amplified in all samples as internal control using the primers CpACS2F (5′-ACT CGATCAACTTCGAGCAAA-3′) and CpACS2R (5′-GCCTA TCCAAAGACCTCGGCCTTCCC-3′). Both ToLCNDVF1/R1 and CpACS2 primers were used in previous works by Sáez et al. (2016). Cycling conditions consisted of incubation at 95°C for 5 min, 45 cycles of 95°C for 5 s and 60°C for 30 s. Relative ToLCNDV levels were calculated using the 2^–ΔCt expression, a variation of the Livak method (Livak and Schmittgen, 2001; Bio-Rad Laboratories, 2006), where ratio (target/reference) = 2^–ΔCt = 2^{–[Ct (viral target) – Ct (reference gen)]}.

PI 604506 and PI 419083 accessions were included in an RNAseq analysis, performed in the frame of a de novo assembly of the zucchini genome project (Montero-Pau et al., 2018), and their transcriptome sequences were used to generate the single nucleotide polymorphism (SNP) panel used here. SNPs were selected by aligning each sequence to the version 1 of the C. moschata cv. Rifu genome (Sun et al., 2017), available at the Cucurbit Genomics Database¹. We used the Bowtie2 tool with the very-sensitive-local argument. Variant calling was performed using Freebayes version 1.0.2 (Garrison and Marth, 2012), excluding alignments from the analysis if they had a mapping quality < 40, alleles with quality under 20, and filtering SNPs with minimum count of 10. A set of 137 SNPs evenly distributed throughout the C. moschata genome (Supplementary Table 1) were selected and used to genotype PI 604506, PI 419083, their derived F₁, and 134 plants of the corresponding F₂ population.

All plants were genotyped using the Agena Bioscience iPLEX^® Gold MassARRAY (Agena Biosciences) system at the Epigenetic and Genotyping unit of the University of Valencia [Unitat Central d’Investigació en Medicina (UCIM), Faculty of Medicine, Spain]. Total DNA was extracted from the tissue of young leaves, using the protocol described above, and quantified and adjusted to 15 ng μl^–1. F₂ genotyping results were run in MAPMAKER 3.0 (Lander et al., 1987; Lincoln et al., 1992) with the Kosambi map function, obtaining the genetic position of each marker.

To identify markers linked to the resistance to ToLCNDV derived from the PI 604506 accession, a QTL analysis was performed using symptoms at 15 and 30 dpi and ToLCNDV relative accumulation at 30 dpi as quantitative traits, and a qualitative score of resistance (0 susceptible phenotype and 1 resistant phenotype) assigned to each plant according to symptoms and viral accumulation. We used the Kolmogorov–Smirnov test to check the normality assumption of traits distribution. Since the traits were not normally distributed, Kruskal–Wallis non-parametric test was used for QTL detection using the MapQTL version 4.1 software (Van Ooijen, 2009), considering as significant associations those with p < 0.05. Since 2^(–ΔCt) values have a skewed distribution, we used the original ΔCt data for QTL analysis. The binary qualitative trait of resistance was also analyzed by logistic regression model, with a significance level of α = 0.05.

In addition, a composite interval mapping approach (CIM, Zeng, 1994) was applied in Qgene 4.0 (Joehanes and Nelson, 2008), using the genetic map previously generated with this F₂. The logarithm of the odds ratio (LOD) threshold was calculated using a 1,000 permutations test per trait, for p < 0.05. The percentage of phenotypic variance explained (R²), the additive and dominance effects, degree of dominance, and the interval position of the QTL in accordance with a 2-U LOD drop was estimated for the highest significant peak LOD. Loci identified with both methods (Kruskal–Wallis and CIM) were considered true QTLs of putative interest.

The previous analysis allowed detecting a major QTL responsible for the resistance in chromosome 8. To confirm this QTL in additional C. moschata-segregating populations and to introgress the candidate region in chromosome 8 of C. moschata in the zucchini (C. pepo) background (the cucurbit crop more severely affected by ToLCNDV), a new set of 19 SNPs of the chromosome 8 candidate region was implemented in a new Agena Bioscience platform. These new SNPs were selected to be useful for both purposes. The transcriptomic sequences of PI 604506, PI 381814, and PI 419083 (obtained in the RNAseq analysis by Montero-Pau et al., 2018) were aligned to the C. pepo genome (Zucchini accession MU-CU-16), available at the Cucurbit Genomics Database², using Bowtie2. Integrative Genomics Viewer (IGV) (Thorvaldsdóttir et al., 2013) was used to detect variations between sequences, and those polymorphic SNPs between resistant (PI 604506 and PI 381814) and susceptible (PI 419083 and MU-CU-16) genotypes were selected. This Agena platform was employed to genotype a subset of 131 plants of the previously genotyped F₂ (PI 419083 × PI 604506), 121 of F₂ (PI 419083 × PI 381814), 31 BC1_{PI 604506}, and 73 of BC1_{PI 381814}.

For further saturation of the candidate region, five additional SNPs, not integrated in the new Agena Bioscience set, were designed with the same requirements and used to genotype the F₂ (PI 419083 × PI 604506) population by high-resolution melting (HRM) (Vossen et al., 2009). PRIMER3 software (Untergasser et al., 2012) was employed to design the oligonucleotides for the HRM analysis. The genomic positions of all these new SNPs (Agena Bioscience platform and HRM markers) and their flanking sequences are shown in Supplementary Table 2.

A new map of chromosome 8 was constructed with 24 SNP markers (3 and 16 SNPs from the first and second Agena platforms, respectively, and 5 HRM), using genotyping results of F₂ (PI 419083 × PI 604506). MAPMAKER 3.0 (Lander et al., 1987; Lincoln et al., 1992) software and the Kosambi map function were employed to generate the new map. The genetic distances of the new map were used in a second QTL analysis, with the F₂ (PI 419083 × PI 604506) population, following the same procedure described above. Means of symptom scores at 30 dpi of plants from F₂ (PI 419083 × PI 381814), BC1_{PI 604506}, and BC1_{PI 381814} populations classified according to the marker classes (a, b, and h for F₂ and h and a for BC1) were analyzed by ANOVA and Bonferroni multiple range tests using STATGRAPHIC Centurion XVI.I statistic software, to evaluate differences between means, considering statistically significant differences when p ≤ 0.01.

An interspecific cross between the ToLCNDV-susceptible C. pepo accession MU-CU-16 (Sáez et al., 2016) and the resistant C. moschata accession PI 604506 provided five F₁ seeds that were germinated as described above. Four seedlings were moved to a greenhouse and selfed to obtain F₂ (MU-CU-16 × PI 604506) generation. The remaining F₁ seedling and 176 plants of F₂ (MU-CU-16 × PI 604506) were screened by mechanical inoculation of ToLCNDV. Symptoms and viral titers were determined by the same procedure described above.

This Cucurbita-interspecific F₂ population was genotyped with the new Agena Bioscience platform and the five HRM SNP markers of chromosome 8. The genotyping results were used to construct a new genetic map of chromosome 8 and to perform an additional QTL analysis as described above.

To obtain a more detailed view of the underlying genomic variation in the candidate region, both C. moschata resistant and susceptible parents (PI 604506 and PI 419083) were fully sequenced. Raw reads are deposited in the National Center for Biotechnology Information (NCBI) under BioProject PRJNA604046. Genomic DNA was obtained from fresh tissue using CTAB extraction, and a pair-end library (2 × 150 bp) was built for each accession. Libraries were sequenced as part of an Illumina HiSeq 2000 lane by Polar Genomics (Ithaca, NY, United States). Reads were cleaned using the ngs_crumbs software³ to eliminate adapters, low-quality bases (Phred quality < 25 in a 5-bp window), reads shorter than 50 bp, and duplicated sequences. Clean reads were mapped against the reference C. moschata genome using bwa-mem (Li, 2013), and variant calling was performed using Freebayes version 1.1.0 (Garrison and Marth, 2012) after filtering reads with a mapping quality cutoff MAPQ lower than 57. To study the potential effect of the genetic changes, SNPs were annotated based on its predicted effect on the gene using SNPEff v4.3 (Cingolani et al., 2012). Differences in sequencing genome coverage between both accessions were studied to explore possible genomic deletions. Read coverage along the candidate region was calculated using samtools v.1.9 (Li et al., 2009), and we checked if coverage deviated from the 99% confidence interval of the observed coverage for each accession assuming a log-normal distribution. Confidence interval for the log-normal distribution was calculated using the function elnorm of R package “EnvStats” (Millard, 2013). In addition to that, the structural variant caller Manta v.1.6 (Chen et al., 2016) was used to check for differential large insertion/deletions. Identification of putative paralogs of the genes in the candidate region was done with OrthoMCL (Li et al., 2003).

Identification of syntenic regions between C. moschata and C. pepo and C. melo was done by nucleotide basic local alignment search tool (BLAST) of each gene within the candidate region of C. moschata against the other two genomes. BLAST hits were filtered using an E value cutoff of 10^–20 and a minimum overlap between sequences of 70%. For C. pepo, to inspect for possible insertion/deletions, a dot plot comparing chromosome 17 region of C. pepo and chromosome 8 of C. moschata was built based on the alignment of both sequences using LAST (Kielbasa et al., 2011). For C. melo, the module of Tripal “SyntenyViewer,” available in cucurbitgenomics.com, was used to visualize the synteny.

New C. moschata and C. pepo genome assemblies have become recently available^{4 ,5} (online availability since November 2019), but after finishing the analysis that we showed here. Our results were checked through alignments with the new assemblies to avoid misinformation.

In addition to the analysis of the genomic sequences, SNPs discovered using the available RNAseq data (Montero-Pau et al., 2018) from the three C. moschata accessions used as parentals in the previous crosses and six additional C. moschata from different origins that exhibited susceptibility to ToLCNDV in previous works (López et al., 2015; Sáez et al., 2016) were also annotated using SNPEff.

The inoculation assay showed that the two C. moschata accessions resistant to ToLCNDV, PI 604506 and PI 381814, remained totally symptomless or with only slight symptoms (score from 0 to 1) at 30 dpi, contrasting with the severe mosaic developed in the susceptible control (score 4), PI 419083 (Figure 1). F₁ plants of the two susceptible × resistant crosses were highly susceptible, displaying a similar symptomatology as PI 419083 at 30 dpi. Conversely, F₁ progeny derived from the cross between the two C. moschata-resistant accessions remained symptomless throughout the essay (Figure 1).

Strong correlation between symptom severity and viral titers was observed (r² = 0.73, p = 0.030) after measuring relative ToLCNDV accumulation by qPCR. In accordance with their resistant behavior, PI 604506, PI 381814, and the F₁ (PI 604506 × PI 381814) had viral titers, on average, 7.8 × 10³ times lower than those of the susceptible control PI 419083 and the two F₁ derived from it (Figure 1).

The fact that F₁ progenies derived from the two susceptible × resistant crosses were susceptible, while the F₁ derived from the resistant × resistant cross was resistant, suggests a recessive genetic control of the resistance in both accessions, controlled by common genes. A further analysis of the genetic control of the resistance was performed in F₂ and BC1 populations.

F₂ and BC1_{PI 604506} populations, derived from the F₁ PI 419083 × PI 604506, segregated for symptoms severity. Table 1 shows resistant:susceptible plants segregation, according to symptomatology at 30 dpi. At the end of the assay, 38 plants of F₂ remained symptomless (score 0), and 96 showed severe symptomatology (scores 2–4). The X² test indicated that this segregation fitted to a 1:3 (resistant:susceptible) ratio expected for a single recessive gene for resistance (p = 0.43) (Table 1). To further characterize the response to ToLCNDV, virus accumulation was estimated in the segregating population F₂ (PI 419083 × PI 604506) by qPCR (Figure 2). On average, viral titer strongly correlated to symptoms severity following an exponential model (r² = 0.82, p = 0.035). All plants developing mosaic, deformation, or short internodes had high viral titers, whereas in the symptomless plants, ToLCNDV accumulation was detected at very low concentrations. On average, the viral accumulation (2^–ΔCt) in susceptible plants was 2.2 × 10³ times higher than in resistant plants. Since viral titer is in concordance with symptoms development, symptom scores were used to phenotype the response to ToLCNDV in plants of the remaining F₂ and BC1 populations. In BC1_{PI 604506}, 33 plants were resistant (score 0) and 26 were susceptible (scores 2–4). This segregation also fitted to a 1:1 ratio expected for a single recessive gene (p = 0.44) (Table 1). In accordance with the occurrence of a single recessive gene controlling the resistance, all plants of the BC1_{PI 419083} generation had severe symptoms at the end of the assay.

Symptom segregation ratios observed in the F₂ (PI 419083 × PI 381814) and BC1_{PI 381814} populations also fitted to one recessive gene for resistance null hypothesis in X² test (Table 1). Forty and 43 plants of F₂ and BC1_{PI 381814}, respectively, remained symptomless (score 0), and 81 and 30 plants showed severe symptoms (scores 2–4), with p = 0.047 and p = 0.16, in both respective populations (Table 1).

In accordance with the F₁ results, the 160 plants of the F₂ derived from the resistant × resistant cross PI 604506 × PI 381814 remained totally symptomless along all the assay.

The F₂ (PI 419083 × PI 604506) population was genotyped with the 137 SNPs markers evenly distributed throughout C. moschata genome and used to construct a linkage map that included 20 linkage groups and spanned a total of 2,681.5 cM of genetic distance, with an average of 22.92 cM between markers (Supplementary Table 1). The linkage map was used to identify QTLs involved in ToLCNDV resistance in C. moschata, based on genotyping and phenotyping results (symptoms scores at 15 and 30 dpi, virus titer at 30 dpi, and the qualitative resistance score) of F₂ (PI 419083 × PI 604506) population. QTL analysis, performed using non-parametric Kruskal–Wallis test (KW) followed by CIM, resulted in the detection of a major QTL in chromosome 8 (Table 2), validated by logistic regression of the qualitative trait of resistance (data not shown). Four QTLs, all located in almost the same genetic position, showed significant association with all the traits evaluated, explaining a proportion between 29.0 and 45.0% of the observed phenotypic variance. All QTLs (ToLCNDVCm_Sy15-8, ToLCNDVCm_Sy30-8, ToLCNDVCm_VT30-8, and ToLCNDVCm_Re-8) were located close to D133 (physical position, 1,366,729 bp), with LOD peaks between 10.06 and 17.31.

To validate the major QTL identified in the previous analysis and to increase marker density in the candidate region, F₂ (PI 419083 × PI 604506) population was genotyped with the new Agena Bioscience-HRM SNPs set of chromosome 8. Twenty-one out of the 24 new markers (Supplementary Table 2) were polymorphic in this population, despite all of them were selected in silico as SNP variants between both parents using the IGV software. Genotyping results were employed to generate a new linkage map in this region (Supplementary Table 2), covering 72.5 cM, with an average distance between consecutive markers of 3.15 cM, and two clusters of linked markers at 0 and 11.7 cM genetic positions. The QTL analysis was performed using the new map and the new genotyping results (using one selected marker of each of the two clusters of completely linked SNPs) (Table 3). ToLCNDVCm_Sy15-8 QTL, associated to the variation of symptoms at 15 dpi, was identified again with both non-parametric Kruskal–Wallis and CIM analysis, near D133 (located at 18.8 cM in this new map) and with similar explained variance, LOD peak, additive, and dominant effects. However, ToLCNDVCm_Sy30-8, ToLCNDVCm_VT30-8, and ToLCNDVCm_Re-8 QTLs, corresponding to traits measured at the end of the assay (30 dpi), when differences are clearer between resistant and susceptible plants, were closely linked to a new marker with both analysis methods. The closest markers (those linked at 11.7 cM DPM37, DMP39, DMP11, DMP10, DMP42, DMP43, DMP44, DMP41, and snp_8202510 markers) are included in the interval position of the same QTLs identified with Kruskal–Wallis and CIM (between 8 and 14 cM) (Figure 3A) and validated with logistic regression of the qualitative trait of resistance, according to their physical and genetic position. The interval of the four QTLs was flanked by DPM34 and D133 markers, with physical positions 561,788 and 1,366,729 bp, respectively. After this further QTL analysis of chromosome 8, the proportion of explained variance was increased, with percentages between 29.5 and 66.0% of R².

The Agena Bioscience SNPs panel of chromosome 8 was used to genotype the BC1_{PI 604506} derived from the PI 419083 × PI 604506 cross, and the F₂ and BC1_{PI 381814} populations derived from the PI 419083 × PI 381814 cross. Mean of symptoms scores at 30 dpi were calculated for each genotypic class of selected SNPs located within the defined QTL interval (DMP35 and DMP39) and compared in Figures 3B–D. The lowest level of symptoms was observed when plants in the three populations had the PI 604506 or PI 381814 homozygous genotype (b), in DPM35 or DPM39 indistinctly. Plants heterozygous (h) or PI 419083 homozygous (a) in both markers displayed significantly more severe symptomatology.

Consistently with the results obtained in F₁ from susceptible × resistant C. moschata crosses, severe symptoms were developed by F₁ C. pepo MU-CU-16 × C. moschata PI 604506 plants at 15 and 30 dpi (Figure 4). This result supports that resistance in PI 604506 has a recessive genetic control. F₂ (MU-CU-16 × PI 604506) plants segregated for symptomatology and viral accumulation. Symptoms, including upward and downward curling and severe mosaic of young leaves, short internodes, and bad distorted development, were observed in 124 and 151 F₂ (MU-CU-16 × PI 604506) plants at 15 and 30 dpi, respectively. The number of resistant plants decreased from 52 to 25 between 15 and 30 dpi. Nine plants had bad development or died in the course of the infection. On average, virus titers determined by qPCR at 30 dpi were in concordance with symptoms development, with mean of relative viral accumulation expressed as 2^(–ΔCt) of 1.04 ± 0.31 and 49,571.67 ± 9,670.31 in resistant and susceptible plants, respectively. The observed segregation proportion was adjusted to the expected ratio resistant/susceptible plants, in case of one recessive gene responsible on the genetic control of resistance to ToLCNDV at 15 dpi (X² = 2.1894, p = 0.14), but not at 30 dpi (X² = 10.312, p = 0.0014).

The genetic map of chromosome 8 generated with the genotyping results of the Agena Bioscience-HRM SNPs in the F₂ (MU-CU-16 × PI 604506) gave a total genetic length of 21.4 cM, with an average genetic distance between successive markers of 0.98 cM (Supplementary Table 2).

The QTL analysis performed in this population show that the QTLs identified in the C. moschata populations were stable in the cross with the C. pepo accession MU-CU-16. ToLCNDVCm_Sy15-8, ToLCNDVCm_Sy30-8, ToLCNDVCm_ VT30-8, and ToLCNDVCm_Re-8 were located in the same region that in C. moschata (Table 4), physically mapped in chromosome 17. The highest R² value (65%) was explained by the ToLCNDVCm_Sy15-8, associated to DMP39 as the nearest marker to the peak LOD. R² values were lower in QTLs related to advanced stages of the ToLCNDV infection, mainly in the viral titer at 30 dpi (ΔCt) trait. In these cases, the nearest markers to the LOD peaks were DMP39 and snp_7926165 in ToLCNDVCm_Sy30-8 (Kruskal–Wallis and CIM tests, respectively), snp_7926165 in ToLCNDVCm_VT30-8 and DMP35, and snp_7926165 in ToLCNDVCm_Re-8 (Kruskal–Wallis and CIM tests, respectively). Logistic regression validate the occurrence of ToLCNDVCm_Re-8 QTL. According to the 2-LOD drop confidence intervals, the position interval where the four QTLs are comapping in chromosome 17 of C. pepo genome (v.4.1) is delimited between DMP34 (7,658,175 bp) and DMP41 (8,165,929).

After both QTL analysis of chromosome 8, a consensus candidate region considered as responsible for ToLCNDV resistance in C. moschata, was established between DMP34 (561,788) and snp_8202510 (1,116,660).

The alignment between the reference assemblies of C. moschata and C. pepo used for mapping purposes in the current paper⁶ and the new assemblies available in November 2019⁷ showed no significant effect on the QTL region studied here (Supplementary Figure 1). Consequently, we keep working with the previous reference versions of both genomes.

A total of 53.2 and 31.5 million genomic clean reads were obtained from PI 604506 and PI 419083, respectively, and approximately more than 97% of them mapped against the C. moschata v.1 reference genome. No large structural variants were found between both accessions, and the read genome coverage was similar among them (Figure 5), which indicates that there are no deletions causing the observed phenotype. Some genomic positions show significant deviations for the expected coverage in both accessions (Figure 5), which could indicate some assembly errors on the reference genome.

After filtering for mapping quality, 28.2 and 18.6 million reads were kept. A total of 1,220,940 SNPs were found to be variable between both parental accessions, and 2,748 were located in the candidate region in chromosome 8. Out of them, nine SNPs had a predicted high impact (either a frameshift or missense variant, a stop codon gain/loss, or a splice site variant) and located within six genes (Supplementary Table 3). Two of these markers are located in the same genes where SNPs used in mapping (snp_7926165 and DMP44) were detected to be linked to ToLCNDV resistance [CmoCh08G001470 encoding a BZIP transcription factor bZIP80 (835,327 to 841,749 bp) and CmoCh08G001770 encoding an unknown protein (1,047,526–1,051,835 bp)]. The remaining seven SNPs with predicted high impact were located in three additional genes of this interval [CmoCh08G001130 encoding a Ribosome inactivating protein (583,200–588,238 bp), CmoCh08G001780 encoding a putative transmembrane protein (1,051,479–1,053,847 bp) and CmoCh08G001880 coding a IQ-DOMAIN 14-like protein (1,097,864–1,102,974)]. In addition, some other SNPs with low, moderate, or unknown modifying effect are placed in genes related to plant virus resistance (Supplementary Table 3).

In addition to the genomic SNPs, the transcriptomic sequences of the three parentals and the six additional susceptible C. moschata accessions provided 731 SNPs in the candidate region, 94 of them were fixed for different alleles in the PI 604506-resistant accession and in the seven susceptible accessions (Supplementary Table 3). PI 381814 transcriptomic sequence had a low coverage in the candidate region, and it was not possible to identify common polymorphisms between the two resistant accessions, PI 604506 and PI 381814. Three SNPs were detected with high predicted effect, all of them were common to those found in the genomic sequences analyzed and were located in three genes (Supplementary Table 3) (CmoCh08G001130 encoding a ribosome-inactivating protein, CmoCh08G001470 encoding a BZIP transcription factor bZIP80 and, CmoCh08G001770 encoding an unknown protein).

The structure of the candidate region was studied in more detail. A whole genome duplication likely occurred in the species that originated the Cucurbita genus (Montero-Pau et al., 2018). In fact, the search for putative paralogs of the genes in the chromosome 8 region indicated that 68 out of the 86 genes in the chromosome 8 candidate region could be assigned to an orthogroup, and 58 of them presented at least one paralog gene. These paralog genes are widespread along the genome (Figure 6), although it seems that there is a conserved duplicated region of chromosome 8 on chromosome 17. Interestingly, some genes of the candidate region have been identified as single copy in chromosome 8 (Supplementary Table 4), without paralog genes in other chromosomes, which is consistent with a major QTL responsible of ToLCNDV resistance.

We also studied the synteny of this region with the susceptible C. pepo, which is phylogenetically closely related to C. moschata. BLAST alignment showed synteny between chromosome 8 region and chromosome 17 from 7,658,023 to 8,205,474 bp of C. pepo (see Figure 7 and Supplementary Table 4). Gene order and orientation is preserved for most genes, but there is one region showing INDELs. Interestingly, the region with a major insertion in C. pepo, from 8,108,962 to 8,113,419 bp, is the region in which the MAD-box transcription factor CmoCh08G001760 maps. This region correspond to position 1,024,011 bp of C. moschata, located between the 5′ untranslated region (UTR) and the first exonic region of this gene. Specific analysis of this C. pepo insertion sequence allowed to detect a long terminal repeats (LTR) retrotransposon of Ty1-copia Retrofit/Ale kind, of 3,692 bp length located from 8,109,186 to 8,113,548 bp. This transposable sequence was previously annotated using the annotation procedure for repetitive sequences described in Montero-Pau et al. (2018). Supplementary Table 6 shows the annotation results and the fasta sequence of the region. Although this insertion is absent in both resistant (PI 604506) and susceptible (PI 419083) C. moschata accessions (Figure 5), many polymorphic SNPs between them are located in this gene, including 5′ UTR and 3′ UTR variants (1,023,872 and 1,047,775 bp), and a missense variant with moderate effect (1,043,369).

BLAST search of the C. moschata QTL region against C. melo found several syntenic regions. In the case of C. melo, highly significant alignments were obtained against chromosome 11 where a major QTL associated with resistance to ToLCNDV is located (Sáez et al., 2017). Results show inversions in SNPs positions between both species, with at least two points of inversion events and loss of information regions (Supplementary Table 2). This syntenic relationship was confirmed with the information displayed by the SyntenyViewer of cucurbitgenomics.org tool. Using chromosome 8 of C. moschata as query genome and location, circular representation showed regions of synteny with eight chromosomes of melon DHL92 (v3.6.1), including the candidate region of chromosome 11 (Figure 8A). Figures 8B,C show the syntenic blocks where ToLCNDV resistance-linked QTLs are located (coded as cmomedB906 and cmomedB910 in the database), the genomic position covered, and the graphic synteny relationship in both blocks. Furthermore, statistical significance of synteny between homologous genes in the candidate region of C. moschata and C. melo is presented in Supplementary Table 5. Seventeen genes are shared by both candidate regions, including the MAD-box transcription factor CmoCh08G001760 and the transmembrane protein CmoCh08G001780 where INDELs or high-effect SNPs have been identified.