DOI: 10.5061/dryad.zkh18937m
title: QTL mapping and marker development for tolerance to sulfur phytotoxicity in melon (Cucumis Melo)

description: 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.    
description: <b>Experimental
 design</b> 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 <i>C. melo</i> 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<sup>o</sup>C +/-
 5<sup>o</sup>c. 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<sup>3</sup> glass greenhouse for sulfur treatment.
 Temperature of the greenhouse were maintained at
 30<sup>o</sup>C +/- 5<sup>o</sup>c. Elemental sulfur (Soil Sulfur:&gt;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<sub>1</sub> seeds failed to germinate in the original study,
 so an additional test was performed that included the parents and new seed
 of the F<sub>1</sub> 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.
 <b>Statistical
 analysis</b> 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<sup>2</sup>) 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. <b>QTL
 mapping</b> 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 <i>C.
 melo</i> 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 <i>C. melo</i> 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). <b>Parental resequencing</b> 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 <i>C. melo</i> 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.   <b>Marker development</b> 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.
description: <b>Supplemental
 files</b> <b>Supplementary table
 S1</b> 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. <b>Supplementary table
 S2</b> 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).
 <b>Supplementary table S3</b> 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>B</b>) in blue, sensitivity allele (<b>A</b>)
 in yellow, heterozygous (<b>H</b>) in gray and missing
 (<b>NA</b>) in white. The significance and magnitude of
 correlation between sulfur response and genotype are listed for each
 marker   <b>Supplementary
 table S4 </b>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).  
 <b>Supplementary table S5</b> 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).  
 <b>Supplementary figure S1</b> 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. <b>Supplementary figure S2</b> 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: (<b>a</b>) cotyledons in test 1, (<b>b</b>) cotyledons in test 2, (<b>c</b>) leaves in test 1 and (<b>d</b>) 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.   <b>Supplementary figure S3</b> 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. <br> 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).