Genome-Wide Association Analysis of Resistance to Pseudoperonospora cubensis in Citron Watermelon

Cultivated watermelon (Citrullus lanatus) is an important vegetable crop to growers and consumers around the world (Food and Agriculture Organization Corporate Statistical Database; FAOSTAT; https://www.fao.org/faostat/en). The crop is an excellent source of moisture (hydration) and sugars as well as beneficial bioactive compounds like citrulline, lycopene, and beta carotene (Nagal et al. 2012; Rimando and Perkins-Veazie 2005). Historically, breeding for resistance to biotic stress in C. lanatus has been slow because of the narrow genetic bottleneck exacerbated by improvement sweeps for yield and fruit quality traits (Levi et al. 2001). Cultivated watermelon is related to other members of the genus Citrullus including citron watermelon (C. amarus), which can be intercrossed with successful pollinations (Chomicki and Renner 2015; Jarret et al. 2017). The C. amarus (citron watermelon) has native variability for resistance to several phytopathogens and pests that affect C. lanatus and is thus an attractive germplasm to improve the sweet-fleshed watermelon (Ali et al. 2012; Levi et al. 2013; Wechter et al. 2012).

Cucurbit downy mildew (CDM), caused by Pseudoperonospora cubensis ([Berkeley and M.A. Curtis] Rostovzev) is a devastating foliar disease that attacks many species within the Cucurbitaceae family (Lebeda and Cohen 2011; Ojiambo et al. 2015). Symptoms associated with CDM include chlorosis, necrosis, and stunting (Lebeda and Cohen 2011). P. cubensis survives on a wide host range of cultivated and wild cucurbits in both temperate and tropical areas around the world (Lebeda 1992; Palti and Cohen 1980). It is an obligate biotrophic oomycete that cannot overwinter in areas with a killing frost (Holmes et al. 2015). In the United States, the pathogen overwinters in warmer regions like the Gulf of Mexico, southern Florida, and Texas (Ojiambo et al. 2015). Annually, in the spring, the sporangia are transported by air currents through the southern to northern states of the United States, infecting crops along the dispersal path (Nusbaum 1948; Ojiambo et al. 2015).

P. cubensis is a geographically and genetically diverse oomycete phytopathogen that can be classified into mating types, clades, pathotypes, and races (Lebeda and Widrlechner 2003; Runge et al. 2011; Thomas et al. 1987; Wallace et al. 2020). There are two mating types (A1 and A2) within P. cubensis needed for sexual reproduction and oospore production (Cohen and Rubin 2012). Mating type A1 preferentially infects members of the Cucumis genus while mating type A2 has a greater affinity for watermelon and butternut squash hosts (Lebeda et al. 2014; Thomas et al. 2017). This diversity and host specialization exhibited by P. cubensis makes cucurbit resistance breeding elusive and is partly responsible for the recent CDM epidemics around world (Cohen et al. 2015; Holmes et al. 2015).

Management of CDM can be achieved through the following methods: Use of downy mildew biosurveillance like spore traps to track sporangia movement and forecast systems to predict times of earliest risk to cucurbit production (Ojiambo et al. 2011; Ojiambo and Holmes 2011; Rahman et al. 2017). Early planting to escape disease outbreaks during the growing season is also feasible, especially when there are cucurbit varieties that can tolerate cold temperatures during the planting season. Fungicide application during the growing season can control CDM disease epidemics and reduce crop loss (Blum et al. 2011; Ojiambo et al. 2010). Growing tolerant and resistant varieties is a more promising option to control CDM epidemics because of the efficient yield protection and less environmental concern with fungicide residues (Cespedes-Sanchez et al. 2015). Therefore, it is vital to actively screen cucurbit germplasm materials for sources of resistance to multiple isolates and races of CDM to avoid breakdown in host disease resistance (Cohen et al. 2015; Holmes et al. 2015).

Recent advances in cucurbit genomics, especially sequencing of sweet fleshed and citron watermelons to develop chromosome-scale genomes, coupled with whole-genome resequencing of most cultivated and wild members of the genus Citrullus, has made genome-wide association study (GWAS) an attractive approach to dissect the genetic architecture of numerous traits in watermelon (Guo et al. 2019; Wu et al. 2019; http://cucurbitgenomics.org/). These genome analysis efforts have accelerated development of millions of single nucleotide polymorphic (SNP) markers that can be deployed in genetic studies such as GWAS. GWAS analysis experiments rely on linkage disequilibrium patterns and historical meiotic events that have accumulated during breeding, selection, and evolution within germplasm diversity panels to detect genomic regions underlying important traits in plants (Ingvarsson and Street 2011; Zhu et al. 2008). In this study, we utilized a densely genotyped U.S. Department of Agriculture’s National Plant Germplasm System collection of citron watermelon to determine the genetic control of CDM resistance (mating type A2 P. cubensis isolate) using GWAS.

Materials and Methods

Plant germplasm and growth conditions.

A total of 122 C. amarus genotypes were evaluated for resistance to CDM infection (Supplementary Table S1). The genotypes were originally received from the U.S. Department of Agriculture’s National Plant Germplasm System in Griffin, GA. Individual genotypes were self-pollinated to reduce heterozygosity that is common in cross-pollinating species. The majority of the genotypes (110) were collected/received from Africa; two were from Asia, while Europe and North America each had five genotypes (Supplementary Table S1). Seeds for each genotype were sown in 50-cell trays (The HC Companies Inc., Twinsburg, OH) filled with Metro-Mix 360 (Sun Gro Horticulture, Agawam, MA). To ensure optimum growth, the seedlings were watered as needed and a balanced water-soluble (NPK: 20-20-20) fertilizer was applied at a rate of 5 g liter⁻¹ (Scotts, Marysville, OH). Average temperature in the seedling growth room was 25.5°C (ranging from 18 to 32.9°C) and 65% relative humidity. Seedlings were grown under 13 h of light provided by balanced LED light fixtures (G.E. Arize Lynk lights 3 Red:1 Blue; HortAmericas, Bedford, TX). The experiment was conducted as a randomized complete block design with two replications over two screening tests.

CDM pathogen maintenance, inoculation, and disease resistance assessment.

A virulent isolate of P. cubensis (WAL-01) was collected in Charleston, South Carolina from symptomatic leaves of a butternut squash cultivar Waltham in 2018 (Toporek et al. 2021). The WAL-01 isolate was maintained through periodic spray inoculations of 2-week-old ‘Waltham’ seedlings grown in 50-cell trays (The HC Companies Inc.) as described in Toporek et al. (2021). To prepare for inoculation of the C. amarus collection, sporulating WAL-01 inoculated ‘Waltham’ leaves were placed in a 1-liter beaker filled with 300 ml of sterile deionized water and shaken vigorously to dislodge the sporangia. The sporangial solution was poured through double-layered cheesecloth into a sterile beaker to remove leaf tissues from the filtrate. A hemocytometer was used to determine the initial concentration of the inoculum, which was then adjusted to a final concentration of 2 × 10⁴ sporangia per ml using sterile deionized water. The 3-week-old C. amarus seedlings were inoculated using a Paasche H-series airbrush sprayer (Paasche Airbrush Co., Chicago, IL) with a no. 3 spray tip at 25 PSI. The inoculum was sprayed to incipient runoff on the adaxial side of previously marked leaves for all genotypes. Inoculated seedlings were immediately moved into a humidity chamber and kept 24 h at 100% relative humidity, 26°C, and ambient light. Seedling trays were then moved to the growth chamber at 25.5°C, 65% relative humidity and 13 h of light until disease rating. Photographs of inoculated leaves were taken for all genotypes at 7 days postinoculation. The photographs were analyzed for percentage leaf area infected (LAI) including both chlorotic and necrotic regions using the image analysis software ASSESS v.2.0 (American Phytopathological Society Press, St. Paul, MN). The ASSESS disease rating for LAI ranged from 0% (healthy green leaf) to 100% (fully chlorotic and/or necrotic leaf) as shown in Figure 1. To process the photographs of inoculated leaves in ASSESS, we used the automatic panel and the “L*a*b*” color space with the rest of the parameters set to default settings.

Download as PowerPoint

DNA extraction.

Seeds of 122 C. amarus accessions were obtained from the U.S. Department of Agriculture’s National Plant Germplasm System to assemble the diversity panel. A single plant (S₀) for each accession was grown to maturity under standard greenhouse conditions and self-pollinated by hand to reduce heterozygosity in the mapping panel. This process was repeated to produce S₂ seed for phenotyping. Genomic DNA was extracted from young leaf tissue of the S₁ plants (a single plant per accession) using DNeasy plant mini kits (QIAGEN, Hilden, Germany). DNA was quantified with a Qubit fluorometer (Thermo Fisher Scientific, Waltham, MA). Genomic DNA (100 ng/μl) was sent to the Roy J. Carver Biotechnology Center at the University of Illinois at Urbana-Champaign for whole-genome resequencing.

Whole genome resequencing.

Shotgun genomic libraries were prepared with the KAPA Hyper Library Construction Kit (Roche, Basel, Switzerland). Libraries were quantified with qPCR and sequenced on one lane of a NovaSeq 6000 (Illumina, San Diego, CA). A NovaSeq S2 Reagent Kit (Illumina) was used for cluster generation and sequencing with 151 cycles from each end of the fragments. The resulting paired-end reads were 150 nucleotides in length. The software bcl2fastq v.2.20 (Illumina) demultiplexed samples, generated “fastq” files, and trimmed adapters from the 3′ end of the reads.

Genotyping.

Perl scripts used to remove duplicated reads were obtained from https://github.com/Sunhh/NGS_data_processing/blob/master/drop_dup_both_end.pl). Low-quality reads were filtered out with the following settings in the tool trimmomatic v.0.38 (http://www.usadellab.org/cms/?page=trimmomatic; Bolger et al. 2014): SLIDINGWINDOW:4:20 LEADING:3 TRAILING:3 HEADCROP:10 MINLEN:40. The tool Burrow-Wheeler Aligner (BWA v.0.7.17; https://github.com/lh3/bwa; Li and Durbin 2009) was used to align the remaining high-quality reads to the C. amarus reference genome of USVL246-FR2 downloaded from http://cucurbitgenomics.org/ftp/genome/watermelon/USVL246/. Reads were tagged if they originated from a single fragment and assigned to a read group with the program Picard v.2.18.7 (http://broadinstitute.github.io/picard). A reference sequence dictionary was created with Picard. The software SAMtools v.0.1.8 (https://sourceforge.net/projects/samtools/files/samtools/0.1.8/; Li et al. 2009) was used to index the reference genome. The GATK Best Practices workflow (Broad Institute) was used for variant discovery and initial quality control filtering (DePristo et al. 2011; Van der Auwera et al. 2013). All steps of the workflow were completed with the software GATK v.3.6 (https://gatk.broadinstitute.org/hc/en-us; McKenna et al. 2010). SNPs with a minor allele frequency <0.05, >10% missing data, >2 alleles, and read depth ± 1 standard deviation from the overall mean depth were removed from the SNP dataset with the program Vcftools v.0.1.15 (Danecek et al. 2011). Based on the 10% missing data filtering cutoff above, no accessions were excluded. The mean read depth across sites was 11.8 with standard deviation = 10.6; SNP markers with a read depth <1.2 or >22.4 were excluded.

Assessment of kinship, population genetic structure, and linkage disequilibrium.

Cryptic relatedness (kinship) among all genotypes was determined by computing a kinship coefficient matrix using centered identity-by-state method (Endelman and Jannink 2012) implemented in the program TASSEL v.5.2 (https://www.maizegenetics.net/tassel; Bradbury et al. 2007) using whole genome SNP markers. Population genetic structure was examined using neighbor-joining clustering and principal component analysis (PCA) algorithms (Price et al. 2006) implemented in TASSEL v.5.2 (Bradbury et al. 2007). The neighbor-joining dendrogram was visualized with the program FigTree (http://tree.bio.ed.ac.uk/software/figtree/). Principal components (PCs) were plotted in two dimensions and visualized using the software platform R (R Core Team 2018) to reveal population stratification. The generated PCs were used as model covariates to account for population structure in GWAS. Five PCs were used in GWAS to account for population structure based on calculation of the cumulative eigenvalue contributions and scree plot analysis. Linkage disequilibrium (LD) across the 122 citron watermelon accessions for every chromosome was assessed as the squared correlation coefficient (r²) between pairs of SNP markers and implemented in the software TASSEL v.5.2 (Bradbury et al. 2007). For every chromosome, LD decay graphs were plotted with physical distance (bp) versus r² for all marker pairs using nonlinear regression and visualized in R (R Core Team 2018; Remington et al. 2001). The LD decay rate was measured as the chromosomal distance where the r² fell to half its maximum value (Huang et al. 2010).

Phenotypic data analysis.

Analysis of variance (ANOVA) for LAI by P. cubensis within the C. amarus genotypes was performed using the “PROC MIXED” procedure in the program SAS v.9.4 (SAS Institute 2013) following the statistical model

where Y_ijk is the LAI value of the ith genotype in the kth replication of the jth screening test, µ is the grand mean, G_i is the effect of the ith genotype, T_j is the effect of the jth screening test, GT_ij is the interaction term between the ith genotype and the jth screening test, rep(T)_jk is the effect of kth replication within the jth screening test, and ε_ijk is the error term assumed to be normally distributed with mean = 0. Normality tests were conducted on the combined residuals of the LAI data using the “PROC UNIVARIATE” procedure in SAS v.9.4 that revealed nonnormal distribution (Shapiro−Wilk test; P = 0.0001). The LAI raw data were converted to best linear unbiased predictors (BLUPs) using the lme4 package in R based on all the components of the statistical model above with all factors set as “random” (Bates et al. 2015). The BLUPs for each genotype were then used as the trait input file for TASSEL v.5.2 to perform GWAS.

The variance components for computing broad-sense heritability (H²) estimates on an entry mean basis were generated with the “PROC VARCOMP” procedure in SAS v.9.4 using the restricted maximum-likelihood estimation method (SAS Institute 2013). Broad sense heritability across tests on an entry mean basis was computed as described in Holland et al. (2003) and shown in the following equation:

where Var(G) is genotypic variance, Var(GT) is the genotype by test variance, and Var(Error) is the residual/experimental error variance. The denominators t and r represent the number of screening tests and replications, respectively.

Genome-wide association analysis.

Genome-wide association analysis was performed using two models, i.e., the general linear model (GLM) and the mixed linear model (MLM) in TASSEL v.5.2 (Bradbury et al. 2007). The GLM methodology includes SNP markers and population structure covariates from PCA in the model. The MLM methodology incorporates SNP markers and the population structure covariates from PCA as fixed effects, and the cryptic relatedness (kinship) matrix was included as a random effect (Yu et al. 2006). The mathematical notation for GLM was:

while for MLM it was:

where Y is the vector of BLUPs for LAI, X is the vector of SNP effects, P is the vector of the population structure effects from PCA, K is the vector of the kinship effects, and ε is the vector of the residuals assumed to be normally distributed with mean = 0. To estimate the amount of phenotypic variability explained (R²) by each of the significant SNP markers, we used the R² values generated in the GLM and MLM statistics output files from TASSEL v.5.2 (Bradbury et al. 2007). The R² values were multiplied by 100 and expressed as a percentage. To determine effects of the biallelic significant SNPs across the C. amarus genotypes, alleles of each peak SNP and their effect on response to CDM was compared using a two-tailed distribution t test and visualized with boxplots generated in R (R Core Team 2018). To correct for multiple marker-trait association testing, we used a false discovery rate ≤ 0.05 (Benjamini and Hochberg 1995). To visualize marker-trait associations, Manhattan and Quantile−Quantile plots were generated in R (Turner 2014).

Candidate gene identification.

To detect putative candidate genes, immediate genomic regions upstream and downstream of the peak SNPs were browsed over an ±80-Kb window in the annotated C. amarus (USVL246-FR2) reference genome (http://cucurbitgenomics.org/ftp/genome/watermelon/USVL246/). Plausible candidate genes included those whose description(s) and gene ontology functions were relevant to molecular plant disease resistance. Additionally, a literature search was conducted to find orthologs and functional gene data in other crop systems similar to the candidate genes uncovered in this research.

Results

Descriptive summary statistics, broad sense heritability, and ANOVA.

Phenotypic variability for LAI averaged over two screening tests ranged from 1.14 to 57.76% (Fig. 2). Broad sense heritability for LAI was moderately high (0.55). The statistical distribution of LAI was left-skewed toward more genotypes resistant to CDM (Fig. 2). Genotype and screening test had significant effects (P = 0.0001) on the observed variability for LAI (Table 1). Genotype and test explained 40.2 and 5% of the observed variability in LAI, respectively (Table 1). The top five most CDM-resistant genotypes had LAI scores of <2% and included PI 482312, PI 482301, PI 482273, ‘Carolina Strongback’ (CSB), and PI 482331. All the top five resistant genotypes were from Zimbabwe in Southern Africa except for CSB, which is from the United States (Supplementary Table S1). The most CDM-susceptible genotypes included PI 596692, PI 271767, Grif 17032, and PI 248774. Genotypes PI 596692 and PI 271767 were obtained from South Africa; Grif 17032 is a germplasm material from the United States; and PI 248774 was collected from Namibia in Southern Africa (Supplementary Table S1).

Download as PowerPoint

Table 1. Analysis of variance showing mean squares percentage of total variance explained, and F and P values for the cucurbit downy mildew leaf area infection of the 122 Citrullus amarus genotypes evaluated over two screening tests under growth chamber conditions

Source of variation	dfa	Mean square	% of TSSb	F value	P value
Genotype	121	575.9	40.2	2.20	0.0001
Test	1	8,733.7	5.0	33.34	0.0001
Genotype × test	121	254.2	17.7	0.97	0.5689
Rep (test)	2	413.3	0.5	1.58	0.2085
Residual error	242	261.9	36.6	NAc	NA

^a d_f, degrees of freedom.

^b % of TSS, Percentage of total sum of squares (total variance) explained.

^c NA, not applicable.

View as image HTML

Population structure and linkage disequilibrium.

A neighbor-joining tree and PCA, based on genetic distances using SNP markers, were used to assess clustering of the genetic variation present in the C. amarus collection (Fig. 3). There was subtle genetic stratification within the core collection especially because 90.2% of the genotypes were collected or received from Southern Africa (Supplementary Table S1 and Fig. 3). The first and second PCs accounted for 22.1 and 4.7%, respectively, of the variability within the SNP marker data (Fig. 3). The third, fourth, and fifth PCs explained 3.6, 2.7, and 2.4% of the genetic variability present in the C. amarus core collection. To account for population structure in GWAS, five PCs that together explained 35.5% of genetic variation in the C. amarus germplasm were used for the MLM analysis. The genetic variation explained by the subsequent PCs after the fifth PC only had marginal incremental values, hence the decision to use only the first five PCs. The decay of LD with physical distance varied among the C. amarus chromosomes and ranged from 10.8 Kb on Ca10 to 106.2 Kb on Ca06 (Supplementary Figs. S1, S2, S3, and S4). Overall, the majority of the C. amarus chromosomes (Ca01 to Ca05, Ca07, and Ca09 to Ca11) exhibited high LD decay rates (<50-Kb decay distances). Chromosomes Ca08 and Ca06 had moderate LD decay rates of 61.6 and 106.2 Kb, respectively (Supplementary Figs. S1, S2, S3, and S4). Overall, presence of moderate/high LD decay rates is expected in populations of cross-pollinated plant species such as watermelon (Ferreira et al. 2000; Guo et al. 2019). Moreover, high LD decay rates of 1 to 100 Kb in maize (cross pollinated) versus longer LD decay distances of 0.1 to 1 Mb in rice (self-pollinated) have been reported (Gore et al. 2009; Remington et al. 2001; Zhao et al. 2011).

Download as PowerPoint

Marker-trait associations and allelic effects on leaf area infection.

Significant SNP markers associated with resistance to CDM were identified using GLM and MLM (Table 2; Supplementary Fig. S5). The significant SNPs explained 24 to 30.2% of the phenotypic variability in LAI (Table 2). Overall, GLM identified more significant markers compared with MLM (Supplementary Fig. S5). This is reasonable because GLM is less stringent and accounts for only population structure covariates in the analysis compared with MLM, which incorporates both population structure and kinship (Price et al. 2006; Yu et al. 2006). Both GLM and MLM uncovered five identical significant SNP markers on chromosomes Ca03, Ca05, and Ca07 (Table 2). MLM identified three unique SNP markers on chromosomes Ca05 and Ca11 (S5_1047147, S5_26089384, and S11_5245781). The top two most significant SNP markers were uncovered using GLM and included S1_37497469 and S2_5605535 on chromosomes Ca01 and Ca02, respectively (Table 2).

Table 2. Details of the peak SNPs and candidate genes associated with resistance to cucurbit downy mildew across 122 Citrullus amarus genotypes evaluated over two screening tests under growth chamber conditions

Chr.	SNP	Position (bp)a	Major allele	Minor allele	MAFb	P value	R2 (%)	Methodc	Candidate gene: annotationd
Ca01	S1_37497469	37,497,469	A	T	0.05	5.23E-12	26.2	GLM	CaU01G30010: LRR–RLK
Ca02	S2_5605535	5,605,535	T	C	0.07	7.44E-12	25.3	GLM	–
Ca03	S3_16504881	16,504,881	A	G	0.05	2.19E-07	28.9	GLM, MLM	–
Ca03	S3_17354658	17,354,658	T	C	0.05	2.19E-07	28.9	GLM, MLM	–
Ca05	S5_732202	732,202	G	A	0.25	9.92E-08	30.2	GLM, MLM	CaU05G00920: LRR–RLK
Ca05	S5_1047147	1,047,147	T	G	0.18	2.92E-07	27.9	MLM	CaU05G01280: WRKY transcription factor
Ca05	S5_26089384	26,089,384	G	A	0.05	3.58E-07	27.5	MLM	−
Ca06	S6_12840690	12,840,690	A	G	0.06	1.23E-11	24.0	GLM	−
Ca07	S7_22984409	22,984,409	C	T	0.05	2.15E-07	28.6	GLM, MLM	−
Ca07	S7_22985428	22,985,428	T	A	0.05	2.15E-07	28.6	GLM, MLM	−
Ca11	S11_5245781	5,245,781	A	T	0.05	1.72E-07	29.1	MLM	CaU11G06010: LRR–RLK

^a bp, basepair.

^b MAF, minor allele frequency.

^c GLM, general linear model; MLM, mixed linear model.

^d LRR−RLK, leucine-rich repeat receptor-like protein kinase.

View as image HTML

The “AA” allele of SNP S3_16504881 on chromosome Ca03 was responsible for an increase in resistance to CDM (mean = 11.9% LAI) compared with a mean of 42.3% LAI among genotypes homozygous for the alternative allele “GG” (P value = 0.007; Fig. 4). Genotypes carrying the “CC” allele of S3_17354658 were more susceptible to CDM (mean = 42.3% LAI) compared with the “TT” allele (mean = 11.9% LAI; P value = 0.007). On chromosome Ca05, the “AA” allele for marker S5_732202 resulted in slightly less resistant genotypic mean of 16.5% LAI compared with allele “GG” (mean LAI of 11.3%) at a P value of 0.039 (Fig. 4). Genotypes with allele “AA” for marker S5_26089384 were more susceptible to CDM (mean = 41.1% LAI) compared with allele “GG” (mean = 11.8% LAI; P value = 0.006). On chromosome Ca07, genotypes with allele “CC” for SNP S7_22984409 were more resistant to CDM (11.7% LAI) than allele “TT” (mean = 39.4% LAI; P value = 0.004). The allele “AA” of marker S7_22985428 was associated with more susceptibility to CDM (mean = 39.4% LAI) compared with genotypes with allele TT, which had a mean LAI of 11.7% at P value = 0.004 (Fig. 4). On chromosome Ca11, the allele “AA” for SNP marker S11_5245781 resulted in high resistance to CDM (mean = 11.7% LAI) compared with the alternative allele of “TT” with mean LAI of 44.9% at a P value of 0.002 (Fig. 4).

Download as PowerPoint

Candidate genes.

Four candidate genes associated with peak SNPs and relevant to disease resistance were identified on chromosomes Ca01, Ca05, and Ca11 (Table 2). Candidate gene CaU01G30010 was located at 78.1 Kb downstream of the peak SNP S1_37497469 on chromosome Ca01 and codes for a leucine-rich repeat receptor-like protein kinase (LRR−RLK; Table 2). Candidate gene CaU05G00920 was found at 23 Kb downstream of the peak SNP S5_732202 on chromosome Ca05 and codes for LRR−RLK gene (Table 2). Candidate gene CaU05G01280 was located at 43.5 Kb upstream of peak SNP S5_1047147 on chromosome Ca05 and codes for a WRKY transcription factor (Table 2). Candidate gene CaU11G06010 was found at 33.4 Kb upstream of the peak SNP S11_5245781 on chromosome Ca11 and codes for LRR-RLK gene (Table 2).

Discussion

CDM, caused by P. cubensis, is an important foliar disease within the Cucurbitaceae family causing yield and quality losses in both temperate and tropical regions around the world (Lebeda 1992; Ojiambo et al. 2015). Developing watermelon cultivars with resistance to CDM is one of the most effective strategies of controlling this important disease. Identification of sources of resistance as well as the underlying genetic basis of CDM resistance is critical to this effort. In this research, we investigated phenotypic variability for CDM resistance present within 122 genotypes of the C. amarus U.S. Department of Agriculture’s germplasm collection. Marker-trait associations relevant to CDM resistance were explored using GWAS.

The ANOVA revealed a strong contribution of genotype to the phenotypic variability in LAI (40.2%, P value = 0.0001). This strong genetic variability enabled identification of CDM-resistant lines (PI 482312, PI 482301, PI 482273, CSB, and PI 482331) in the C. amarus collection. Over two CDM screening tests, genotype reactions to inoculations with P. cubensis varied from 1.14 to 57.76% LAI. These resistant accessions can be targeted as parental materials to improve watermelon for resistance to CDM. Broad sense heritability on entry mean basis was moderately high (0.55), indicating presence of adequate genetic variability for CDM resistance within the C. amarus U.S. Department of Agriculture’s germplasm collection. Genotype PI 482331 performed above the population mean for total root length in a large watermelon root traits screening study while CSB offers additional value with its high resistance to soilborne pathogens and suitability for grafting (Katuuramu et al. 2020; Keinath et al. 2019).

Historically, breeding and genetic studies for CDM have been conducted in the commercial vegetables within the Cucumis genus, particularly melon and cucumbers (Barnes 1948; Call et al. 2012; Kenigsbuch and Cohen 1992; Thomas et al. 1988). CDM is an emerging foliar disease in watermelon production. Research efforts on resistance breeding for CDM in watermelon is scarce. In this research, we identified eight significant SNP markers associated with resistance to CDM on chromosomes Ca03, Ca05, Ca07, and Ca11. The SNP markers explained variability in CDM ranging from 27.9 to 30.2%. Allele “GG” for SNP S5_732202 had the lowest mean score for LAI of 11.3% indicating that genotypes with this SNP allele were more likely to be resistant. Identification of these CDM resistance loci with large phenotypic variability explained values will be vital for marker-assisted breeding in watermelon.

Four of the peak SNPs from this research associated with CDM resistance in the C. amarus collection were located close to candidate resistance genes. The functional annotation of the candidate genes was LRR−RLK and the WRKY transcription factor. Upon infection, molecular plant disease resistance mechanisms include several processes such as pathogen detection, signal transduction, and defense response (Glazebrook 2005; Kou and Wang 2010). LRRs are a hallmark feature of the majority of all cloned resistance genes and have been widely described for their roles in plant host defense (Kourelis and van der Hoorn 2018). The LRR−RLK gene family have been reported to provide increased resistance to rice blast and the bacterial pathogen Xanthomonas oryzae pv. oryzae in rice (Caddell et al. 2017; Peng et al. 2009). Additionally, silencing the LRR1/PR4b gene proteins showed that LRR1 and PR4b are necessary for defense response to Pseudomonas syringae in pepper (Hwang et al. 2014) The second family of candidate genes uncovered in this research is a WRKY transcription factor. WRKY transcription factors have been reported to play key roles in signal transduction in the presence of pathogen or insect attack, ultimately boosting plant host immunity. The WRKY transcription factor has been reported to improve resistance to oomycetes pathogens in pepper and soybean (Cheng et al. 2020; Cui et al. 2019) as well as to fungal pathogens in rice (Peng et al. 2012). The functional gene data for the different crop systems above provides additional support to the potential value of the genomic regions uncovered in this research for CDM resistance breeding in watermelon.