Elemental sulfur is an effective, inexpensive fungicide for many foliar pathogens, but severe phytotoxicity prohibits its use on many melon varieties. Sulfur phytotoxicity causes chlorosis and necrosis of leaf tissue, leading to plant death in the most sensitive lines, while other varieties have little to no damage. A high-density, genotyping-by-sequencing (GBS)-based genetic map of a recombinant inbred line (RIL) population segregating for sulfur tolerance was used for a quantitative trait loci (QTL) mapping study of sulfur phytotoxicity in melon. One major (qSulf-1) and two minor (qSulf-8 and qSulf-12) QTL were associated with sulfur tolerance in the population. The development of Kompetitive Allele-Specific PCR (KASP) markers developed across qSulf-1 decreased the QTL interval from 239 kb (cotyledons) and 157 kb (leaves) to 97 kb (both tissues). The markers were validated for linkage to sulfur tolerance in a set of melon cultivars. These KASP markers can be incorporated into melon breeding programs for introgression of sulfur tolerance into elite melon germplasm.

Experimental design

A previously described RIL population (Branham et al. 2018) consisting of 170 lines generated from a cross of MR-1 and Ananas Yok’neum (AY) was evaluated for elemental sulfur tolerance. The Israeli cantaloupe cultivar Ananas Yok’neum was the sulfur tolerant parent and the inbred C. melo line MR-1 (Thomas 1986) was the sensitive parent (Figure 1). Two independent greenhouse tests of the parents and population were initiated in May and June 2017. Each test was planted in a randomized complete block design with two replicates of 5 plants each. Lines were seeded into Metromix 360 (Sun Gro Horticulture, Agawam, MA) in 50-cell propagation trays (Hummert International, Earth City, MO) and allowed to grow to the 2-3 fully expanded leaf stage in a sulfur-free glass greenhouse. Temperature of the greenhouse were maintained at 30^oC +/- 5^oc. Seedlings were fertilized the day prior to sulfur treatment by soaking trays in a liquid fertilizer solution (3g Peters water soluble fertilizer per liter) (Scotts, Maryville, OH, USA). Trays were transferred into a temperature-controlled, 650 m³ glass greenhouse for sulfur treatment. Temperature of the greenhouse were maintained at 30^oC +/- 5^oc. Elemental sulfur (Soil Sulfur:>99% purity, National Garden Wholesale, Vancouver, WA, USA) was vaporized using two sulfur burners (Wilmod Sulfur Evaporator WSE75; Zoetermeer, Netherlands) for 4 hours nightly. The sulfur burners were ~2 m above the work benches, suspended 0.75 m below a circulation fan.  The two burners were on adjacent ends of the greenhouse. Each sulfur burner vaporized approximately 1.2 g of sulfur per night. On the fifth day, lines were evaluated for sulfur tolerance by recording percent necrosis for both the most damaged cotyledon and true leaf on every plant (Figure 2). The percent necrosis for each RIL (cotyledon and true leaf) was averaged from evaluations of twenty plants (2 tests x 2 reps x 5 plants). F₁ seeds failed to germinate in the original study, so an additional test was performed that included the parents and new seed of the F₁ hybrid. Two replicates of ten seeds per line were planted in a greenhouse trial in June 2018. An additional test of thirty melon accessions (cultivars and PIs) were evaluated to test the utility of the sulfur markers in a variety of germplasm. Two replicates of five seeds each were planted in a greenhouse trial in March 2019. These additional studies followed the same protocols described above.

Statistical analysis

Pearson’s correlation (r) of line means between tests and between tissue types was calculated with the stats package of R version 3.4.1 (R core team 2017). Broad-sense heritability (H²) of sulfur tolerance, measured as percent affected leaf area (chlorosis and/or necrosis), was determined separately for each tissue type as the RIL variance divided by the total variance in percentage affected tissue area using variance components estimated with a linear mixed model in ASReml-R v3.0 (Gilmour et al. 2009). The model included RIL, test, interaction of RIL and test, replicate, and tray nested within test as random effects.

QTL mapping

We used the previously published (Branham et al. 2018) high-density genetic map developed for this population for all QTL mapping analyses, which included 5,663 imputed, binned SNPs across the 12 chromosomes (=linkage groups) of the C. melo genome (Garcia-Mas et al. 2012). Haley-Knott regression (Haley and Knott 1992) was used for multiple QTL mapping (MQM) with the stepwiseqtl function (Zeng et al. 1999; Broman and Speed 2002; Broman and Sen 2009) of Rqtl (Broman et al. 2003). The optimal QTL model based upon penalized LOD score (Manichaikul et al. 2009) was chosen through an automated forward and backward search algorithm. One thousand permutations of a two-dimensional, two QTL scan were used to calculate penalties and the genome-wide significance threshold for QTL detection. Multiple QTL models were visualized through LOD profile plots generated from forward selection using standard interval mapping with Haley-Knott regression (Haley and Knott 1992). Distributions of necrosis percentage of both cotyledons and true leaves did not meet the assumptions of parametric interval mapping, therefore the non-parametric model of the scanone function (Kruskal and Wallis 1952; Kruglyak and Lander 1995) was used for QTL verification. Genes within the 1.5-LOD interval of the major QTL were identified using the functional annotation of the C. melo reference genome v3.5.1 (Garcia-Mas et al. 2012), which was obtained through batch query at http://cucurbitgenomics.org/ (Zheng et al. 2019). In addition to using the functional annotation provided with the reference genome to search for candidate genes, conserved domains of genes were identified using the National Center for Biotechnology Information’s batch CD search (CDDv3.16 database) (Marchler-Bauer and Bryant 2004; Marchler-Bauer et al. 2011, 2015, 2017).

Parental resequencing

Genomic DNA was extracted from young leaf tissue of both parental lines (MR-1 and AY) using a DNeasy Plant Mini kit (Qiagen, Venlo, Netherlands) and sent to the Roy J. Carver Biotechnology Center at the University of Illinois at Urbana-Champaign for whole-genome resequencing. A Hyper Library construction kit (Kapa Biosystems, Roche, Basel, Switzerland) was used to prepare shotgun libraries for each parental DNA. Libraries were quantified by qPCR, pooled, and sequenced on one lane of a NovaSeq 6000 (Illumina, San Diego, CA) with a NovaSeq S2 reagent kit. Paired-end reads (150bp) were demultiplexed with bcl2fastq v2.20 Conversion software (Illumina). Adaptors were trimmed from the 3’ end of the reads. Duplicated read pairs were removed with perl scripts (https://github.com/Sunhh/NGS_data_processing/blob/master/drop_dup_both_end.pl). Low quality reads were removed with trimmomatic v0.38 (Bolger et al. 2014). The remaining high-quality reads were aligned to C. melo reference genome v3.5.1 (Garcia-Mas et al. 2012) with BWA v0.7.17 (Li and Durbin 2009). Picard v2.18.7 (http://broadinstitute.github.io/picard) was used to assign reads to a read group, tag reads originating from a single DNA fragment, and to create a reference sequence dictionary. The reference genome was indexed with Samtools v0.1.8 (Li et al. 2009). The Genome Analysis Toolkit (GATK v3.6) was used for SNP calling following the best practices for variant discovery (McKenna et al. 2010, DePristo et al. 2011, Van der Auwera et al. 2014). SNPs were filtered with Vcftools v0.1.15 (Danecek et al. 2011) to remove those with any missing data, heterozygous genotypes for either inbred parent, and/or genotype quality score of less than 30. SNPs within the major QTL region were functionally annotated with ANNOVAR version 2017 Jul 16 (Wang et al. 2010). Genes with missense or nonsense mutations and mutations to the promotor (less than 1 kb upstream of the start codon) were considered candidate genes.

Marker development

The parental whole-genome resequencing data was used to design markers to saturate the region of the major sulfur tolerance QTL. Eighteen SNPs from across the major QTL region were developed into KASP markers (Supplementary table S1) using “KASP™ by design” services from LGC Genomics (Teddington, Middlesex, UK). PCR reactions (5 µL volume) consisted of 0.07 µL of primer mix (LGC Genomics; fluorophore-labelled allele-specific forward primers and a reverse primer), 2.5 µL of 2x master mix (LGC Genomics) and 20 ng of sample DNA. A standard thermal cycler was used for a touchdown PCR reaction with a 94°C hot-start activation step for 15 min, then ten cycles of 94°C (20 s) and a starting annealing temperature of 61°C that dropped by 0.6°C each cycle. Twenty-six additional cycles of 94°C for 20 s and 55°C for 60 s followed the touchdown steps. Fluorescence was quantified with a Stratagene Mx3005P (Agilent Technologies, Santa Clara, CA) quantitative PCR system at 25°C. Fluorescence values were used to cluster samples into genotypes with MxPro v4.10 software associated with the qPCR machine. Marker linkage to sulfur tolerance in the RIL population was assessed through QTL mapping both alone (KASP markers only) and combined with the binned GBS SNPs following the same procedures as described above. Thirty accessions (cultivars and plant introductions) were evaluated for sulfur tolerance and genotyped with the KASP markers. Correlation between the markers and sulfur phenotype of the accessions was assessed

through analyses of variance (ANOVA) with the aov function (Chambers et al 1992) in R.

Supplemental files

Supplementary table S1 Sequence information for the KASP markers, including: SNP ID, physical position of the SNP, primer sequences, SNP flanking sequence, and nucleotides of the sulfur tolerant (T) and sensitive (S) alleles.

Supplementary table S2 Phenotypic data used for QTL mapping: RIL means of percentage of damaged (chlorotic or necrotic) area after vaporized elemental sulfur treatment of cotyledons (cot) or leaves (leaf) in test 1 (_t1), in test 2 (_t2), and across tests (cot or leaf).

Supplementary table S3 Melon accession genotypes at KASP markers developed across qSulf-1. Markers are named according to physical position (bp) on chromosome 1. Genotypes are color coded, with individuals homozygous for the sulfur tolerance allele (B) in blue, sensitivity allele (A) in yellow, heterozygous (H) in gray and missing (NA) in white. The significance and magnitude of correlation between sulfur response and genotype are listed for each marker

Supplementary table S4 Chromosomal location and functional information for genes that collocated with the major QTL for sulfur tolerance, including the position of the start and stop codons within the chromosome (cs), gene ontology (GO) code and term, and conserved domains and features found within the gene (NCBI).

Supplementary table S5 Functional annotation of candidate gene polymorphisms between the parents in the QTL interval of qSulf-1, including the chromosome, physical position (in bp), parental alleles, polymorphism location relative to the gene (ie. upstream, exonic, downstream, etc.), distance from the gene (in bp), type of exonic mutation (synonymous or nonsynonymous), and detail (which exon and the nucleotide and amino acid changes).

Supplementary figure S1 Interaction plot showing evidence for epistasis between qSulf-1 and qSulf-12. Alleles from the sulfur sensitive parent (MR-1) are ‘AA’ and from the sulfur tolerant parent (AY) are ‘BB’. The circle represents the mean percent damage of individuals in the population with the labelled genotypes. The plus signs indicate the standard error.

Supplementary figure S2 Logarithm of odds (LOD) scores for forward model selection of up to seven QTL associated with mean percentage of damaged (chlorotic or necrotic) area after vaporized elemental sulfur treatment of: (a) cotyledons in test 1, (b) cotyledons in test 2, (c) leaves in test 1 and (d) leaves in test 2. The initial scan shows the likelihood of the first QTL being located at each SNP in the genome (linkage group=chromosome) with subsequent scans showing the LOD of an additional QTL with the effects of the previous QTL(s) controlled for in the model. The dashed line marks the genome-wide significance threshold.

Supplementary figure S3 Interaction plot showing evidence for epistasis between qSulf-1 and qSulf-8. Alleles from the sulfur sensitive parent (MR-1) are ‘AA’ and from the sulfur tolerant parent (AY) are ‘BB’. The circle represents the mean percent damage of individuals in the population with the labelled genotypes. The plus signs indicate the standard error.

sulfQ1_annovar.csv Functional annotation of all gene polymorphisms between the parents in the QTL interval of qSulf-1, including the chromosome, physical position (in bp), parental alleles, polymorphism location relative to the gene (ie. upstream, exonic, downstream, etc.), distance from the gene (in bp), type of exonic mutation (synonymous or nonsynonymous), and detail (which exon and the nucleotide and amino acid changes).