Melon Genome Regions Associated with TGR-1551-Derived Resistance to Cucurbit yellow stunting disorder virus

Cucurbit yellow stunting disorder virus (CYSDV) is one of the main limiting factors of melon cultivation worldwide. To date, no commercial melon cultivars resistant to CYSDV are available. The African accession TGR-1551 is resistant to CYSDV. Two major quantitative trait loci (QTLs) have been previously reported, both located near each other in chromosome 5. With the objective of further mapping the gene or genes responsible of the resistance, a recombinant inbred line (RIL) population derived from the cross between TGR-1551 and the susceptible cultivar ‘Bola de Oro’ was evaluated for resistance to CYSDV in five different assays and genotyped in a genotyping by sequencing (GBS) analysis. The major effect of one of the two QTLs located on chromosome 5 was confirmed in the multienvironment RIL assay and additionally verified through the analysis of three segregating BC1S1 populations derived from three resistant RILs. Furthermore, progeny test using the offspring of selected BC3 plants allowed the narrowing of the candidate interval to a 700 kb region. The SNP markers identified in this work will be useful in marker-assisted selection in the context of introgression of CYSDV resistance in elite cultivars.


Introduction
The Cucurbit yellow stunting disorder virus (CYSDV) is a Crinivirus of the family Closteroviridae [1,2]. The most important cucurbit crops are infected by this virus, and melon (Cucumis melo L.) is one of the species most severely affected.
CYSDV was first reported in the United Arab Emirates [3] and subsequently spread throughout the Middle East and the Mediterranean Basin [4][5][6][7][8][9][10], North and Central America [9,[11][12][13], and China [14]. Nowadays, this virus has become a significant threat throughout the cucurbit production regions in the southern United States, Mexico, and Central America [15] and has been reported as the most economically important of the viruses affecting cucurbit production in the southwestern United States [16]. The objective of this work was to better depict the genetic architecture of the CYSDV resistance derived from TGR-1551, with the aim of identifying molecular markers useful for MAS in breeding programs dedicated to the introgression of this resistance into elite cultivars. To attain this objective, several populations, a recombinant inbred line (RIL) population and derived advanced backcrosses, obtained from an original cross between TGR-1551 and the highly CYSDV-susceptible Spanish cultivar 'Bola de Oro' have been genotyped and evaluated for CYSDV resistance.

Multienvironment Phenotyping for Resistance to CYSDV of the RIL Population
The RIL population (F 7 /F 8 generations) derived from the cross TGR-1551 × 'Bola de Oro' (BO) was evaluated for resistance to CYSDV in five different years. Disease assessment was carried out by the evaluation of symptom development in all the assays. Additionally, the estimation of virus titer in plants by qRT-PCR was done in 2012 and 2013 (Supplementary Table S1).
All BO plants showed symptoms in all the assays, whereas TGR-1551 plants remained asymptomatic. Plants of the F 1 showed slight symptoms, later and less severe than those observed in BO plants. Viral accumulation was not detected in the resistant parent, TGR-1551. The virus titer in BO was either significantly higher than in the F 1 or similar, depending on the assay (Supplementary Table S1).
Correlations between symptom scores in the different years were moderate to high and highly significant in all cases (values between 0.367 and 0.839, p < 0.005 in all comparisons; Supplementary Table S2). Correlations between virus titer values measured by qRT-PCR in the 2012 and 2013 assays were high and highly significant in all comparisons (0.750-0.824, p < 0.0001; Table S2). In 2012 and 2013, it was possible to relate symptoms and virus titer. Correlation between symptom score and virus titer was high and highly significant for both years (0.513-0.776, p < 0.0001; Supplementary Table S2). The severity of symptoms does not always correlate with virus titer [32]. However, a significant correlation between symptom development and viral accumulation has been reported in melon viral diseases caused by viruses of different families, such as the geminivirus Tomato leaf curl New Delhi virus [33] or the potyvirus Watermelon mosaic virus [34]. Concretely, in the case of CYSDV, previous studies found a relationship between the time of symptoms appearance and viral accumulation, with delayed development being observed in plants with lower titers [35]. A positive correlation between CYSDV viral titer and whitefly transmission in cucurbit hosts has been reported [35]. Thus, the fact that lower symptoms found in TGR-1551-derived resistant plants correspond to lower viral titers is epidemiologically important in the reduction of inoculum sources in fields.
The effects of genotype (G), environment (E), and the interaction between them (GxE) were explored with data of the 46 coincident RILs included in the 2009, 2010, 2012, and 2013 assays. Genotype, environment, and the interaction between them each had a significant effect (p < 0.0001). Although differences among RILs were observed due to the effect of the interaction GxE, symptoms were significantly lower on average in 2009. On the contrary, symptom scores were significantly higher in 2012. B. tabaci confined in clip-cages was the inoculation method used in 2009, and plants were kept in pots in a glasshouse once inoculated. It seemed that the cultivation and growing systems led to lower symptom development. Massive inoculation was used in 2010, 2012, and 2013, and plants were grown in a plastic greenhouse. The highest temperatures were registered in 2012 (data not included), which could explain the higher symptom scores observed in this assay.
Resistance to CYSDV in TGR-1551 has been reported as associated with restriction of the virus movement and/or reduction in virus replication [23], which could be limiting the cascade of events associated with symptom induction in TGR-derived resistant plants [36]. Our results show that despite the high correlation between symptoms and viral accumulation, there are symptomless plants in which the virus is detected and, also, plants with symptoms of yellowing that could be attributed to CYSDV in which the virus cannot be detected in the moment of the analysis. The accumulation of the particles and restriction of the virus movement in TGR-1551 has previously been reported [23,27], which could explain the virus titer found in symptomless plants. On the other hand, as previously stated, yellowing symptoms caused by CYSDV are sometimes difficult to visually distinguish, not only from symptoms produced by other criniviruses, but also from nutritional disorders or phytotoxicity [5,28], which can explain the occurrence of symptomatic plants with a negative qRT-PCR. All these factors condition the evaluation of resistance against CYSDV and highlight the need to use efficient disease assessment methodologies in breeding programs. In this work, the availability of evaluation data (both symptoms and virus titer) for most of the RIL lines (several plants per RIL) allowed an accurate classification of each RIL as resistant or susceptible.

QTL Analysis With the RIL Population
A QTL analysis was performed using the GBS1 and phenotypic data of the RILs evaluated in each assay. The map used was the SNP map previously constructed with the GBS1 data from the RILs [37]. Chromosome 5 systematically appeared in all the assays as involved in resistance to CYSDV (Table 1). Two genomic regions were identified in this chromosome: one explained variation in symptom severity, while the other explained variation in virus titer detected by qRT-PCR. The QTL associated with symptom severity explained between 27% and 62% of the variance for this trait and, considering the intervals defined for each assay, spanned between 6,810,744 and 24,296,585 bp. For virus titer detected by qRT-PCR, the QTL identified explained between 49% and 53% of the variance, spanning between 24,791,006-27,121,114 bp (Table 1). Both intervals were near a region for which polymorphic markers were not available from the GBS1, of approximately 10 cM and 1.1 Mb (between 66.4 and 75.9 cM, from 25,229,866 to 26,193,386 bp) ( Figure 1). This would possibly interfere with the results of the QTL analyses.  The physical position (v3.6.1) is defined by the position of the markers flanking the QTL interval. 3 Closest marker to the LOD peak. Significance level in the Kruskal-Wallis test was 0.0001 for all the markers. 4 Mean of the genetic class TGR-1551 for the corresponding marker. 5 Mean of the genetic class 'Bola de Oro' for the corresponding marker. 6 Higher logarithm of the odds score. 7 Percentage of phenotypic variance explained by the QTL.
evaluation data gathered in all the assays). See Section 3 for details. 2 Interval position of the putative QTL on the genetic and physical maps according to an LOD drop of 2. The physical position (v3.6.1) is defined by the position of the markers flanking the QTL interval. 3 Closest marker to the LOD peak. Significance level in the Kruskal-Wallis test was 0.0001 for all the markers. 4 Mean of the genetic class TGR-1551 for the corresponding marker. 5 Mean of the genetic class 'Bola de Oro' for the corresponding marker. 6 Higher logarithm of the odds score. 7 Percentage of phenotypic variance explained by the QTL.

Phenotyping for Resistance to CYSDV and QTL Analysis in Three Selected BC1S1 Progenies Derived from Selected Resistant RILs
Three of the most resistant and vigorous RILs evaluated in most of the assays, RIL 278, RIL556, and RIL110, were selected to produce segregating generations ( Figure S1).
These RILs were backcrossed to the susceptible parent, BO, and one plant from each of the BC1 generations was selfed to generate the BC1S1 populations. The selection of the resistant RILs was done  is defined by the position of the markers flanking the QTL interval. 3 Closest marker to the LOD peak. Significance level in the Kruskal-Wallis test was 0.0001 for all the markers. 4 Higher logarithm of the odds score. 5 Percentage of phenotypic variance explained by the QTL.
As aforementioned, considering evaluation data gathered in all the assays, each RIL was classified as resistant or susceptible. The QTL obtained with this phenotypic value (response in Table 1) explained 44% of the variance, with the interval extending from 24,957,179 to 26,993,475 bp. This interval included the 10 cM region lacking markers ( Figure 1), which could account for the size of the interval.

Phenotyping for Resistance to CYSDV and QTL Analysis in Three Selected BC 1 S 1 Progenies Derived from Selected Resistant RILs
Three of the most resistant and vigorous RILs evaluated in most of the assays, RIL 278, RIL556, and RIL110, were selected to produce segregating generations ( Figure S1).
These RILs were backcrossed to the susceptible parent, BO, and one plant from each of the BC 1 generations was selfed to generate the BC 1 S 1 populations. The selection of the resistant RILs was done before the GBS1 was available. Genotyping by sequencing (GBS) results of the RILs confirmed later that RIL110 and RIL556 were both homozygous for the TGR allele for the whole chromosome 5 (Table S3). RIL278 was also TGR homozygous for almost the whole chromosome 5, with the exception of the region below 27,628,291, for which it was homozygous for BO alleles. BC 1 S 1 plants derived from RIL278 and RIL556 were phenotyped for resistance in 2016. Plants of the susceptible parent BO started showing symptoms from the first evaluation dates, which increased with time ( Figure 2). Viral accumulation was confirmed in all plants by qRT-PCR. Some plants of the resistant parent, TGR-1551, showed slight spotting, although it was later confirmed by qRT-PCR that this was not caused by the presence of the virus (Supplementary Table S1). Similarly, the yellowing observed on some evaluation dates in plants of the RILs 278 and 556 could not be attributed to the virus. Along with the parents, some plants of each BC 1 generation (selected RIL × BO) were evaluated. BC 1 plants showed milder symptoms than BO at all evaluation dates for both RILs. The first symptoms in the BC 1 derived from RIL556 appeared delayed with respect to BO. The virus was detected in all BC 1 plants but at lower levels than in BO. For subsequent analyses, the availability of new GBS data (GBS2) allowed the design of molecular markers in the 10 cM region lacking markers in GBS1 ( Figure 1) to better delimit the candidate interval. This new set of markers (CYSDV2) was used to genotype the 134 BC1S1556 plants; the evaluation period in this population extended to 8 weeks post infection (wpi), and virus titer was measured in all plants. To perform the QTL analysis of this population more accurately, a new map of the region was constructed with the genotyping data of the plants with both marker sets, CYSDV1 and CYSDV2, which covered 62.5 cM, corresponding to the region from 6,412,266 to 27,806,568 bp ( Figure 1). The new QTL analyses, using this map and the BC1S1556 population phenotyping results, confirmed the occurrence of a major QTL in this region of chromosome 5 ( Table 3). The intervals for the QTL considering the different evaluation dates were overlapping, and altogether covered the region between 25,619,503 and 26,688,074 bp. The closest marker to the LOD peak was cysdv63. Almost the same interval (25,619,503-26,629,653 bp) corresponded to the QTL obtained when assigning a phenotype to each RIL, considering the symptom evaluation and the virus titer (Table 3). In any case, the percentage of the variance explained for the different evaluation dates showed an increasing trend with time (from 33% to 56%). Evaluation 5 (8 wpi) was the date that explained the Genotypes are expressed as follows: TGR: TGR-1551; RIL: selected RIL, RIL278, 556, or 110; BC 1 : first backcross to 'Bola de Oro' of the selected RIL; BO: 'Bola de Oro'; BC 1 S 1 : selfing progeny from the BC 1 , classified according to the genotype for marker cysdv22 for RIL278 and marker cysdv63 for RILs 556 and 110, where BC 1 S 1 -B is homozygous for TGR-1551 allele, BC 1 S 1 -H is heterozygous, and BC 1 S 1 -A is homozygous for 'Bola de Oro' allele. Bars represent standard error. Different letters in the same evaluation date indicate significant differences among genotypes (p < 0.05, LSD test).
The first set of markers (derived from the RIL GBS) was designed for its use in the Agena Bioscience genotyping platform (CYSDV1; Supplementary Table S4). This set included markers covering the genomic regions of the two QTLs obtained from the RIL analysis. The QTL analysis was performed with symptom evaluation in the 133 BC 1 S 1 278 plants and the 134 BC 1 S 1 556 plants. (Table 2). A QTL was identified for both families in the region around markers cysdv22 (24,613,028 bp) and cysdv24 (26,766,636 bp), which explained 21% and 58% of the variation in symptom evaluation in both families, respectively ( Figure 1, Table 2). These results confirmed the QTL positions found in the RIL populations but supported the need to saturate the region between markers cysdv22 and cysdv24.
For subsequent analyses, the availability of new GBS data (GBS2) allowed the design of molecular markers in the 10 cM region lacking markers in GBS1 (Figure 1) to better delimit the candidate interval. This new set of markers (CYSDV2) was used to genotype the 134 BC 1 S 1 556 plants; the evaluation period in this population extended to 8 weeks post infection (wpi), and virus titer was measured in all plants. To perform the QTL analysis of this population more accurately, a new map of the region was constructed with the genotyping data of the plants with both marker sets, CYSDV1 and CYSDV2, which covered 62.5 cM, corresponding to the region from 6,412,266 to 27,806,568 bp ( Figure 1). The new QTL analyses, using this map and the BC 1 S 1 556 population phenotyping results, confirmed the occurrence of a major QTL in this region of chromosome 5 ( Table 3). The intervals for the QTL considering the different evaluation dates were overlapping, and altogether covered the region between 25,619,503 and 26,688,074 bp. The closest marker to the LOD peak was cysdv63. Almost the same interval (25,619,503-26,629,653 bp) corresponded to the QTL obtained when assigning a phenotype to each RIL, considering the symptom evaluation and the virus titer (Table 3). In any case, the percentage of the variance explained for the different evaluation dates showed an increasing trend with time (from 33% to 56%). Evaluation 5 (8 wpi) was the date that explained the highest percentage of the variation and had the highest LOD. The interval for the QTL in this evaluation date comprised the region between 25,982,529 and 26,629,653 bp, thus suggesting that this late phenotyping allowed a more accurate characterization of the population. A second QTL associated to the qRT-PCR was identified in the interval between markers cysdv17 and cysdv18 (22,651,076-24,296,585 bp), which explained 44% of the variation for the trait (Table 3).
To evaluate the effect of the symptom-related QTL in both progenies, BC 1 S 1 278 and BC 1 S 1 556 plants were grouped according to their genotype for the closest marker to the higher peak available in each case (cysdv22 for BC 1 S 1 278 and cysdv63 for BC 1 S 1 556). The average symptom score in each evaluation date was calculated for each genotypic class (Figure 2). In the case of BC 1 S 1 556, virus titer was also available. Symptom scores from the first evaluation dates were significantly higher in plants homozygous for the BO allele with respect to the rest of the genotypes in BC 1 S 1 . Symptoms in homozygous plants for the BO allele were slightly lower than in BO, corresponding to a slightly lower virus titer (evaluated in BC 1 S 1 556). Differences between heterozygous BC 1 S 1 plants and those homozygous for the TGR allele were significant for all evaluation dates in BC 1 S 1 278, but they were only significant at 8 wpi in BC 1 S 1 556. In the case of heterozygote BC 1 S 1 plants, the average symptom score was similar to those in the corresponding BC 1 , while virus titer (measured in the 556 progenies) was lower in the former. Average scores in plants homozygous for the TGR allele were higher than those observed in the resistant parent, TGR, and in the corresponding RILs (only in the final evaluation date in the case of BC 1 S 1 278), given that, as previously explained, the observed yellowing in the RILs was not attributable to the presence of virus.
To confirm the previous results, BC 1 S 1 progeny derived from RIL110 were phenotyped for resistance to CYSDV in 2017 (Supplementary Table S1). Percentage of plant infection did not reach 100% in the susceptible control BO at the end of the assay (8 wpi) and remained at levels similar to those obtained in the BC 1 S 1 278 assay, thus suggesting lower inoculum pressure with respect to the BC 1 S 1 556 assay (Figure 2). Similarly, the virus titer detected in the susceptible control in the BC 1 S 1 110 assay was significantly lower at 8 wpi than in the BC 1 S 1 556 assay. Resistant controls in the three assays displayed responses similar to those of TGR-1551, and the resistant RILs remained virus-free. On the contrary, symptom scores in BC 1 110 and BC 1 278 were higher than in BC 1 556, while virus titer was comparable in both assays.
The 146 BC 1 S 1 110 plants were classified according to their genotype for marker cysdv63 (Figure 2). Differences between the three genotypic classes were significant from the second evaluation date to the end of the assay. Symptoms in BC 1 S 1 plants homozygous for the BO allele were lower than in the assay with family 556, supporting the results obtained in the susceptible control BO. In any case, symptoms in BC 1 S 1 plants homozygous for the BO allele were slightly lower than in the susceptible parent, BO. Similarly, symptoms in BC 1 S 1 heterozygous for cysdv63 were lower than in the heterozygous BC 1 for most evaluation dates. Plants homozygous for the TGR allele showed slight symptoms, corresponding to the detection of virus titer, while TGR-1551 and RIL110 remained asymptomatic and virus-free.
Differences found among the three BC 1 S 1 populations (all segregating for the TGR-1551 candidate introgression in chromosome 5) could be explained by differences in the genetic background for the rest of the genome. Moreover, as previously stated, there is an important effect of environmental conditions on TGR-1551-derived resistance to CYSDV [27]. The most important differences in this work were found in heterozygous genotypes. These results agree with differences found in previous works, where the F 1 generation derived from the initial cross TGR-155 × BO behaved either as resistant [18] or susceptible [27]. The heterozygous plants in the three segregating BC 1 S 1 populations analyzed here showed (for the percentage of infection at the later evaluation date, 8 wpi) an intermediate response between TGR-1551 and BO. However, the evolution in these heterozygous plants from the segregating generations supported the previous findings that established later evaluation dates as more representative [27], as percentage of infection on the initial sampling dates did not differ from that of plants that were homozygous for the TGR-1551 allele in some of the assays.

Progeny Test Narrowing the Interval of the Major QTL and Confirmation of Its Effect in Advanced Backcross Selfing Populations
To complement the analysis of the RIL and BC 1 S 1 RIL derived populations, and to evaluate the effect of the introgression in plants with BO genetic background, 200 plants of the BC 3 (TGR-1551 × BO derived) population were genotyped with marker set CYSDV1. Fifty-two plants were selected for the presence of different introgressions in the candidate region of chromosome 5. These plants were additionally genotyped with a set of 124 markers evenly distributed throughout the genome. As an average, 82% of the genome in the 52 plants corresponded to regions homozygous for BO alleles. Fifteen BC3 plants were selected to obtain the BC 3 S 1 progenies. The selection was based on the presence of different introgressions in the candidate region of chromosome 5 (Table 4) and prioritized a high percentage of the BO genetic background for the rest of the genome.
The phenotyping/genotyping results of these offspring were compatible with the interval for the major QTL obtained in the analyses of the BC 1 S 1 556 plants (Table 4, Supplementary Table S1). Recombinants in the region allowed the delimitation of the candidate region. The interval obtained considering the different evaluation dates in BC 1 S 1 556 plants covered the region between 25,619,503 and 26,688,074 bp. Segregation among BC 3 S 1 plants derived from BC 3 166 and 198, and susceptibility of descendants from BC 3 105, 141, and 146, confirmed the location of the resistance gene between markers cysdv63 and cysdv65 (25,982,529-26,629,653 bp). All BC 3 progenies segregated according to the presence of the resistance gene(s) in the region between markers cysdv63 and cysdv65, producing susceptible selfing progenies (24, 28, 37, 78, 96, 105, 141, and 146), or segregating offspring (15, 19, 64, 95, 159, 166, and 198) (Table 4). The proposed interval explained the phenotype of the RILs; excluding recombinants in the region, there was cosegregation between the phenotype and the genotype in 87% of the homozygous lines for the rest of the RILs (Table S1).
The candidate interval for the major CYSDV resistance QTL obtained in this work (25,982,629,653 bp) contains 57 predicted genes, 51 of them annotated (Table S5). Several of the annotated genes in this region have resistance-related functions. In fact, the candidate interval overlaps with a 760 kb region with the highest concentration of resistance genes in the melon genome [38]. This region has been shown to be highly polymorphic at the intra-and interspecific levels, thus explaining differences in resistance found in different melon genotypes [38]. The interval proposed here contains three genes annotated as 'disease resistance protein' (MELO3C031332.2, MELO3C004320.2, and MELO3C031556.2), as well as related-resistance genes such as a receptor-like cytosolic serine/threonine-protein kinase (MELO3C004315.2) or two nucleotide-binding site-leucine-rich repeat proteins (MELO3C004319.2 and MELO3C031325.2). The most frequent class was the TMV resistance protein N-like (12 of the annotated genes in the candidate interval), a TIR-NBS-LRR gene having homology with resistance genes [38]. Moreover, this region includes the virus aphid transmission resistance gene (Vat) [38], carried by TGR-1551 [39]. The major QTL for resistance to powdery mildew caused by Podosphaera xanthii (Castagne) U. Braun & N. Shishkoff races 1, 2, and 5 derived from the resistance source used here, TGR-1551, has also previously been mapped to this region [40]. This will be an advantage in the breeding program for the introgression of TGR-1551 resistance to CYSDV, Aphis gosypii, or powdery mildew, given that the three resistances would be transferred simultaneously. In fact, the Yellow Canary breeding melon line 'Carmen', derived by backcrossing from the initial cross between TGR-1551 and 'Bola de Oro', confirmed this instance [31]. Although resistance to CYSDV was the only selection character used in the backcrossing program, the line 'Carmen' resulted in also being highly resistant to powdery mildew and Aphis gosypii. Furthermore, quality was assessed in 'Carmen'; despite carrying the introgression associated with the resistances, its commercial quality was confirmed. TGR-1551 is also resistant to Watermelon mosaic virus (WMV), and molecular markers tightly linked to the resistance QTL with a major effect on chromosome 11 have been developed [37]. Resistance to WMV can also be introgressed into elite cultivars in breeding programs with TGR-1551 as a donor parent. Apart from a major QTL of chromosome 11, a minor QTL also associated with resistance to WMV derived from this source is located in this same region of chromosome 5 [37]. Resistance to CYSDV in cucumber (Cucumis sativus L.) has been reported in accession PI 250147, characterized by the absence of symptoms [41]. One of the QTLs, which explained the highest percentage of this resistance, seemed to be linked in the repulsion phase to two loci conferring resistance to powdery mildew, derived from different sources. These resistances have been mapped to linkage group 5, in a region syntenic to melon chromosomes 9 and 10. A method to introgress resistance to CYSDV in melon chromosome 9 from an unknown source has been patented [42], but the region is not syntenic with the CYSDV-resistance region in cucumber.
The region of the melon genome interval containing the cluster of putative resistance genes has been studied in detail to analyze the genomic variability [43] in a study that confirmed the difficulty in obtaining an accurate and complete sequence by NGS due to the multiple, highly similar genes clustered in a relatively short region. This may have been the cause of the lack of markers in this region in the initial map developed from the GBS of the RIL populations. In any case, the availability of the subsequent GBS data allowed the identification of markers polymorphic between BO and TGR-1551 to design the marker set CYSDV2, which led to the narrowing of the candidate interval for the TGR-1551-derived resistance. Markers in both sets, CYSDV1 and CYSDV2, were initially implemented for their use in the Agena Bioscience genotyping platform. Some of the markers of both sets have been adapted for their analysis by high-resolution melting (HRM). Concretely, the HRM protocol for markers in the candidate regions has been set (Supplementary Table S4), which allows the efficient identification of the polymorphism between the BO and the TGR alleles by PCR-based markers. The availability of this type of marker, suitable for MAS in the context of breeding programs, is of special interest in the case of introgression of resistance to CYSDV. As previously stated, the yellowing symptoms caused by crinivirus are undistinguishable. Similar symptoms are caused by different criniviruses [28], and these symptoms are commonly confused with the effects of nutritional disorders or phytotoxicity [5]. Thus, the use of MAS is an effective method to circumvent these difficulties.
The analysis of different TGR-1551-derived generations allowed the narrowing of the candidate interval of resistance to CYSDV to a 700 kb region in chromosome 5. This resulted in the availability of molecular markers tightly linked to the resistance. The PCR-based markers developed here are an efficient resource for use in TGR-1551-derived resistance introgression in commercial melons. The breeding program for the introgression of TGR-1551-derived resistance to CYSDV, WMV, and powdery mildew in 'Bola de Oro' and 'Piel de Sapo' backgrounds has already been initiated. Future work will include further analysis in advanced backcrossed progenies with the purpose of better understanding the effect of the minor QTL obtained in some of the analyses in the region of markers cysdv14 and cysdv17. The linkage between this region and the region of the major QTL limits the availability of segregant generations appropriate to discriminate against the effect of this minor QTL. Moreover, RNAseq analysis will be carried out in order to identify differentially expressed genes in the candidate regions.

Plant Material
The following populations obtained from an original cross between the resistant melon line TGR-1551 and the susceptible Spanish melon cultivar 'Bola de Oro' (BO) have been evaluated against CYSDV (Supplementary Figure S1): 1.
Three BC 1 S 1 progenies derived from crosses between three resistant RILs and BO; 3.
Fifteen BC 3 S 1 progenies derived from 15 BC 3 (× BO) selected for their genotype in the candidate region. . Not all RILs were evaluated in all the assays due to seed availability. The different backcross and selfing populations were phenotyped in additional screening assays according to the methodology described below.

Inoculation Method
Healthy whitefly colonies (Bemisia tabaci) reared on plants of the susceptible Spanish melon accession 'ANC-57' were used in all the inoculations and assays. Controlled CYSDV-infected plants of this accession showing clear and typical symptoms of CYSDV were used as inoculum source of the virus in all the inoculations.
RIL assays conducted in 2009 and 2011 were carried out by using B. tabaci confined in clip-cages following the methodology developed by [4], in which 60 whiteflies were fed on the inoculum source for 48 h (acquisition time) and then transferred to young plants (1-2 true leaves) to feed for 48 h (transmission time). Once inoculated, plants were sprayed with imidacloprid to kill the whiteflies and transplanted to pots in the glasshouse. Assays conducted in 2010, 2012, and 2013 were carried out by massive inoculation. Plants of each RIL at the stage of 2-3 true leaves were transplanted to a plastic greenhouse where they were randomly distributed interspersed with infected BO plants. One week after transplanting, viruliferous whiteflies were released in the plot to enhance virus transmission. In this case, any spraying against whiteflies was done over the course of the experiment.
The clip-cage method was used in the inoculation of subsequent generations.

Disease Assessment
Virus symptoms were visually assessed once a week over a one-month period ( [29]. During the Spring 2016 and Spring 2017 assays, three BC 1 S 1 progenies derived from crosses between three resistant RILs (278, 556, and 110) and BO were evaluated. Virus inoculations were carried out by using B. tabaci confined in clip-cages as explained above. Once inoculated, plants were transplanted to the greenhouse. Evaluation of CYSDV was based on the percentage of leaves showing typical CYSDV symptoms of infection (number of leaves affected/total leaf number). qRT-PCR analyses to estimate the virus titer in the plants were carried out every week (2016) or every two weeks (2017).
Additionally, in 2019, BC 3 S 1 progenies (12 to 20 plants each) derived from 15 BC 3 plants selected according to their genotype in the candidate region of chromosome 5 were inoculated with CYSDV using the clip-cage method. Then, plants were transferred to the greenhouse and evaluated by symptoms and qRT-PCR as described for the RIL110 BC 1 S 1 population.

Statistical Analyses
Analysis of variance (ANOVA) using data from four of the RIL evaluations was carried out. Correlations between pairs of environments were estimated by using the Pearson correlation coefficient. Comparisons of pairs of means in different assays were performed using the LSD test with a probability level of p < 0.05. Statistical analyses were conducted using Statgraphics Centurion XVI.I software (StatPoint Technologies, Inc., Warrenton, VA, USA).

Markers and Genotyping Methods
Total DNA was extracted from young leaves following the method described by [44] with minor modifications [45]. DNA concentration was measured using spectrophotometry in a Nanodrop ND-1000 Spectrophotometer v.3.5 (LabTech International, Heathfield, UK) and adjusted to 10 ng/µL for the genotyping analysis.
The SNPs used in this study were designed from data originating from different GBS analysis. The whole RIL population (148 RILs), both parents (BO and TGR-1551) and their F 1 , constituted the plant materials for the first GBS (GBS1 assay). A total of 16 SNPs, located on the candidate region on chromosome 5, were selected from this GBS and constituted the first set of markers (CYSDV1 ;  Supplementary Table S4). Preferably, markers were selected that corresponded to those included in the map used in the analysis. When necessary to evenly cover the region, markers were chosen among the markers discarded for the map because of segregation distortion (only two markers). An additional GBS experiment, GBS2, conducted to perform genetic diversity studies (including BO and TGR-1551, among many other genotypes), was also used as a source for 24 extra SNPs in order to increase the resolution of the candidate region (marker set CYSDV2; Supplementary  Table S4). All the SNPs were implemented for their use in the Agena Bioscience platform (https: //agenabio.com/products/massarray-system/); use of this platform was carried out in the Epigenetic and Genotyping unit of the University of Valencia, Unitat Central d´Investigació en Medicina (UCIM), Spain. Marker set CYSDV1 was used to genotype the BC 1 S 1 , BC 3 , and BC 3 S 1 generations. Generations BC 1 S 1 556 and BC 3 S 1 were also genotyped with marker set CYSDV2. BC 3 plants were genotyped with an existing panel of 124 SNPs evenly distributed throughout the genome, also implemented in the Agena Bioscience platform. This SNP set had been previously validated in populations derived from ibericus × acidulus melon crosses [45][46][47]. Some of the markers in CYSDV1 and CYSDV2 sets have been adapted to a PCR-based protocol for their analysis by high-resolution melting (HRM) (Supplementary Table S4).

QTL Analyses
QTL analysis was performed using the genotyping results of GBS1 and RIL CYSDV phenotypes. Phenotypic data used were the percentage of infection at the end of each assay and qRT-PCR values, when available. Moreover, a final phenotypic value was assigned to each RIL, considering the five evaluations, so that each RIL was classified as resistant or susceptible. First, the Kruskal-Wallis nonparametric test was used for QTL detection (MapQTL version 4.1 software) [48]. In addition, a composite interval mapping (CIM) approach was performed [49], using a window size of 15 cM and five cofactors (Windows QTL Cartographer v.2.5-009) [50]. For qualitative traits (response, where plants were classified as either susceptible or resistant), qGene v.4.4.0 was the software used [51]. The LOD threshold was determined by a permutation test (1000 cycles). Loci detected by both Kruskal-Wallis and CIM methods were considered robust QTLs. QTL interval was defined as a 2-LOD drop from the peak. The phenotypic effect, expressed as the percentage of phenotypic variance explained, R 2 , and the additive (when possible) and dominance effects were estimated for each QTL. The map used in the QTLs analyses for the RIL population was that constructed with the GBS data of the RILs [37].
Additional QTL analyses were performed with BC 1 S 1 progenies from RILs 278 and 556. BC 1 S 1 progenies from RIL 278 were genotyped with a set of SNP markers selected from the GBS of the RILs (CYSDV1), covering the candidate region in chromosome 5. A new set of markers (CYSDV2) that saturated the candidate region was designed; it was used, together with CYSDV1, to genotype the BC 1 S 1 progenies from RIL 556. A new map of the region was constructed for each of the progenies, and these maps were used for the respective QTL analyses. The software used was MAPMAKER 3.0 [52]. The map was generated using the Kosambi map function.

Conclusions
No commercial varieties resistant to CYSDV are presently available. The African accession TGR-1551 has been reported as resistant to Bemisia tabaci and to the virus itself. The work reported has allowed the narrowing of the candidate interval for the major QTL associated with resistance to CYSDV derived from TGR-1551 to a region of approximately 700 kb. The SNP markers provided here are useful in marker-assisted selection (MAS) in breeding programs aimed at the introgression of CYSDV resistance. CYSDV is exclusively transmitted by its insect vector, Bemisia tabaci. Moreover, yellowing symptoms caused by CYSDV are difficult to distinguish from symptoms caused by other viruses, and yellowing provoked by crinivirus is easily confused with the effects of nutritional disorders, phytotoxicity, or other causes. Thus, the availability of markers for MAS in the context of melon breeding programs is essential for accelerating the introgression of resistance to CYSDV into elite cultivars.
Supplementary Materials: Supplementary materials can be found at http://www.mdpi.com/1422-0067/21/17/5970/ s1. Figure S1. Generations used and experiments carried out. Table S1. Symptom evaluation and qRT-PCR result after inoculation with CYSDV of the RIL population and advanced derived backcrosses. Table S2. Correlations between parameters used for disease assessment (symptom scores and qRT-PCR Ct values) in the RILs population in the different assays. Table S3. Genotype of RILs 110, 278 and 556 for the markers in GBS 1. Table S4. Marker sets (CYSDV1 and CYSDV2) used in genotyping of progenies. Table S5. Gene loci in the region of the candidate interval for the major QTL for resistance to CYSDV derived from TGR-1551.