Genome-wide association mapping of Sclerotinia sclerotiorum resistance in soybean using whole-genome resequencing data

DOI: https://doi.org/10.21203/rs.2.14709/v1

Abstract

Sclerotinia stem rot (SSR), caused by Sclerotinia sclerotiorum (Lib.) de Bary, is an important cause of yield loss in soybean. Although many papers have reported different loci contributing to partial resistance, few of these were proved to reproduce the same phenotypic impact in different populations. In this study, we identified a major quantitative trait loci (QTL) associated with resistance to SSR progression on the main stem by using a genome-wide association mapping (GWAM). A population of 127 soybean accessions was genotyped with 1.5M SNPs derived from genotyping-by-sequencing (GBS) and whole-genome sequencing (WGS) ensuring an extensive genome coverage and phenotyped for SSR resistance. SNP-trait association led to discovery of a new QTL on chromosome 1 (Gm01) where resistant lines had shorter lesions on the stem by 29 mm. The impact of this QTL was even more significant in the descendants of a cross between two lines carrying contrasted alleles for Gm01. Individuals carrying the resistance allele developed lesions almost 50% shorter than those bearing the sensitivity allele. These results suggest that this region harbors a promising resistance QTL to SSR that can be used in soybean breeding program.

Introduction

Sclerotinia stem rot (SSR) is a significant disease that causes yield and quality loss in soybean in the northern United States and Canada. This disease is caused by Sclerotinia sclerotiorum, a necrotrophic Ascomycota, capable of infecting more than 408 different species (Boland and Hall, 1994). The fungus infects the plant via the flower then spreads through the stem causing bleaching, severe wilting and shredding of tissue (Bolton et al., 2006). SSR was reported as the second most important disease-causing yield losses in Canada in 1994 and in the USA in 1994, 2004 and 2009 (Wrather and Koenning, 2006; Koenning and Wrather, 2010). However, the impact of this disease is very unpredictable from year to another because fungal development is highly influenced by temperature and humidity (Mila and Yang, 2008). Its impact could be reduced by using chemical or biological control, but results can be variable as these methods can fail when disease incidence is higher than 50% (Zeng et al., 2012). The best results can be achieved when several preventive treatments are applied each year even when SSR doesn’t pose a threat due to unfavorable climate conditions. Considering these facts, enhancing the genetic resistance of soybean cultivars seems to be the most effective solution to reduce the detrimental impacts of SSR.

The evaluation of SSR resistance is quite challenging in variable environmental conditions. However, a reliable inoculation method was developed by Bastien et al. (2012) wherein a mycelium suspension is applied on flower buds in controlled greenhouse conditions. It has been shown to produce consistent results and was used to investigate the genetic determinants of SSR resistance is soybean (Huynh et al., 2010; Bastien et al., 2014; Iquira et al., 2015).

To date, complete resistance has yet to be reported in soybean. Partial resistance is controlled by multiple genes or quantitative trait loci (QTL). Numerous mapping studies have been conducted and have identified more than 114 QTL via conventional biparental mapping (Kim and Diers, 2000; Arahana et al., 2001; Huynh et al., 2010; Vuong et al., 2008; Guo et al., 2008; Zhao et al., 2015 and Kandel et al., 2018). Although this method has been widely used for QTL mapping, it is still limited to the genetic diversity present in the two parents. More recently, with the advancement of genotyping technologies, it was possible to screen quantitative partial resistance in multiple soybean lines with thousands of markers using GWAM. Using this method, more than 130 QTLs have also been reported in different populations (Bastien et al., 2014; Iquira et al., 2015; Moellers et al., 2017; Zhao et al., 2015; Wei et al., 2017 and Wen et al., 2018). Such number of loci raise some questions about their credibility especially when fewer of these were proved to reproduce the same allelic effect in different genetic backgrounds. One explanation is that some of these QTLs identified based on different methods of evaluation, could be confused with an escape or avoidance mechanisms and not genuinely related to the real physiological resistance to SSR (Bastien et al., 2014; Kim and Diers, 2000). As a proof, the only QTL proved to reproduce the same phenotyping effect in a biparental cross was identified on chromosome Gm15 based on resistance evaluation under a controlled environment (Bastien et al., 2014). These results suggest that a reliable phenotyping method is a key factor in this study.

Compared to biparental mapping, diversity panels offer a lower level of linkage disequilibrium (LD) between markers and QTLs. Hence, for GWAM, a higher marker density is needed depending on population size and diversity. For higher QTL detection power, the LD between the QTL and any flanking markers should be higher than 0.8. To achieve such a coverage, Bastien et al. (2014) estimated that at least 12,900 SNPs in the pericentromeric regions and 55,700 SNPs in the telomeric region would be needed for a total of over 68K well-distributed SNPs to cover the entire genome. For mapping SSR resistance loci in soybean, many attempts were made to achieve such converge using different genotyping approaches like genotyping by sequencing (GBS) (Bastien et al., 2014; Iquira et al., 2015; Wei et al., 2017) or specific locus amplified fragment sequencing (SLAF-seq) (Zhao et al., 2015). To date, the largest number of informative SNPs was achieved using the SoySNP50K array in two studies. One obtained 35,683 SNPs on 466 accessions (Moellers et al., 2017) and the other achieved 31,600 and 35,708 SNPs, respectively, in populations of 915 improved lines and 405 soybean landraces (Wen et al., 2018). It is likely that the marker coverage obtained in these most recent papers still falls short of the number needed to ensure exhaustive genome coverage.

One alternative to the previously used genotyping approaches is whole-genome sequencing (WGS). However, this approach is still expensive, especially when using large populations. In previous work, Torkamaneh et al. (2017a) proposed a two-step approach termed “scanning and filling”. In a first step, a large population can be genotyped at tens of thousands of SNP loci (using GBS or an array). In a second step, WGS can be performed on a subset of these lines (e.g. 20%) and these can serve as a reference panel to impute millions of SNP markers onto the entire set of accessions.

In this work, we used such a combined GBS and WGS genotyping approach to genotype an association panel (comprising elite Canadian soybean lines) at millions of SNPs. We then used this exhaustive marker dataset to perform GWAM in the association panel to identify QTLs responsible for partial resistance to SSR in Canadian soybean.

Materials and Methods

Association mapping panel

The association mapping panel used for this study was composed of 127 lines exhibiting a wide variation in their response to SSR. These were chosen from a larger group of 530 accessions (cultivars/advanced breeding lines) representative of the genetic diversity in Canadian soybean based on previous work (Torkamaneh et al., 2017a). These 127 lines belonged to maturity groups (MGs) ranging from 000 to II except for one line, Williams 82, from MG III. A series of six checks were also included: three cultivars known to offer a good level of SSR partial resistance (Karlo RR, Maple Donovan and S19–90), two moderately resistant cultivars (OAC Bayfield and Williams 82) and one highly susceptible cultivar (Nattosan) (Bastien et al., 2012).

Validation panel

A total of 47 F6:8 lines segregating for the candidate QTL region on chromosome Gm01 were selected to serve as a validation panel. These lines were generated from a cross between the partially resistant Maple Donovan and the susceptible OAC Bayfield.

Phenotyping

Lines were evaluated for SSR partial resistance using the cotton pad method described in Bastien et al. (2012). For the association panel, the phenotypic data are those previously reported by Bastien et al. (2014). Briefly, plants were sown in a greenhouse in a randomized complete block design with four blocks separated in time (25 Sept, 6 Nov,18 Dec 2009 and 29 Jan 2010). Experimental units consisted of a total of six plants grown in three 6-L pots (two per pot). The same experimental design was used to characterize the validation panel but with only two blocks separated in time (4 May and 7 Sept 2017) and four plants per experimental unit.

The potting mix was prepared using a mixture of black earth (50%), perlite (30%) and Promix (20%) (Premier Tech Horticulture, Rivière-du-Loup, QC, Canada). At sowing, seeds were inoculated with RhizoStick® inoculant (Becker Underwood, Ames, IA). Plants were grown under a 16-h photoperiod and the day/night temperature was maintained at 26/22°C.

The inoculum was prepared from strain NB–5 (provided by Dr. S. Rioux of CEROM, Quebec City, QC, Canada) as described in Bastien et al. (2012). Briefly, S. sclerotiorum was cultured in potato dextrose broth (PDA) for three days until almost reaching saturation. Inoculation was performed once the plants started to flower. First, the suspension was homogenized for 30 s in a blender. Then, pieces (2.7 x 5.5 cm) of cotton pad were soaked in the suspension. The inoculum was applied on the petiole of the lowest node bearing flowers. After inoculation, plants were transferred to a different greenhouse where day/night temperatures were 22°C/18°C and high humidity was maintained at 2.5 g/m3 with a fogging system. For the validation panel, all plants were inoculated on the same day, while for the association panel, several days were needed because of differences in flowering date. Lesion length was measured 7 d after inoculation.

SNP genotyping and imputation

The association panel was a part of a larger set of 530 Canadian elite lines on which we had previously performed GBS (ApeKI, as per Elshire et al., 2011) over time (Bastien et al., 2014; Torkamaneh et al., 2017a). To maximize data quality and uniformity, all reads (940M 108-bp single-end Illumina reads) were run on an improved SNP-calling pipeline (Fast-GBS; Torkamaneh et al., 2017b) and on a more recent version of the Williams 82 reference genome (Wm82.a2.v1) (Schmutz et al., 2010). This resulted in a catalogue of 150K SNPs on which all missing data were imputed using BEAGLE v5 (Browning and Browning, 2007) as per Torkamaneh and Belzile (2015). Subsequently, a subset of 102 lines was subjected to whole-genome sequencing (WGS) (Torkamaneh et al., 2017a) and the resulting SNP catalogue (>4M SNPs) was used as a reference panel to perform large-scale imputation of missing loci. The SNP data (4.1M loci) for the 127 lines of the association panel were extracted and filtered using vcftools v0.1.16 (Danecek et al., 2011). We retained SNPs with a minor allele count (MAC) ≥ 1 and a minor allele frequency (MAF) ≥ 0.05. Linkage disequilibrium (LD) was estimated (using r2) for all marker pairs in a sliding window of 50 Kb using PLINK 1.9 (Purcell et al., 2007).

Analysis of population structure

Given the large size of the SNP catalogue (almost 1.5M SNPs), pruning was performed using PLINK 1.9 (Purcell et al., 2007) to remove markers in high LD (r2 ≥ 0.9). The resulting set of 85K SNPs was used to assess population structure using fastSTRUCTURE (Raj et al., 2014) with K set between 1 and 12. The most likely number of subpopulations was estimated using the chooseK tool from fastSTRUCTURE (Raj et al., 2014).

Genome-wide association analysis

In view of GWAM, only SNPs having a minor allele frequency (MAF) ≥ 5% in the association panel were used, and this resulted in a catalog of close to 1.5M filtered SNPs. An association mapping analysis for SSR partial resistance was performed using the phenotypic (mean lesion length) and genotypic data described above with the Genomic Association and Prediction Integrated Tool (GAPIT version 2) (Tang et al., 2016). To correct for false-positive associations, a mixed linear model (Q + K model) taking into account both population structure (Q matrix) and relative kinship (K matrix) was used. The Q matrix (for K = 6) was derived from fastSTRUCTURE while the K matrix was generated in TASSEL. Marker-trait associations were deemed significant when the measured p-values were below a critical p-value corresponding to a false discovery rate (FDR) of 0.1 (Benjamini and Yekutieli, 2005).

QTL validation

A codominant cleaved amplified polymorphic sequence (CAPS) marker was designed to genotype one of the candidate SNPs (Gm01: 5594765) residing in the haplotype block containing the peak SNP on Gm01. Two specific primers (5’-GTTGTATGGAAGTGCAACTAAAGTTCT–3’ and 5’- GGTACTTTTTCTTACCTTAC GATGA–3’) were used to amplify an 800-bp region encompassing the targeted SNP. The two alleles can be distinguished by digesting the resulting amplicon with NmuCI. The PCR product derived from the allele associated with partial resistance to SSR (present in Maple Donovan) will be cut once while the product obtained after amplification of the allele from the susceptible parent (OAC Bayfield) is not cut. All 47 F6:8 lines of the validation panel (described above), along with the two parental lines were genotyped using this CAPS marker.

Genomic landscape around the peak SNP

LD values from Plink (Purcell et al., 2007) were extracted for a 2Mb window around the most significant associated SNP and LD blocks were visualized using Haploview (V4.2) (Barrett et al., 2005). Information about the genes found in the LD block containing the peak association were obtained from SoyBase (www.soybase.org). Functional annotation of nucleotide variation in the region was explored using SnpEFF (Cingolani et al., 2012).

Results

SSR resistance in lines of the association panel

Lesion length was measured seven days after inoculation on young flower buds, and the mean value for each genotype is shown in Supplementary Table S1. As illustrated in Figure 1, lesion lengths were found to range broadly, from as low as 29 mm to a maximum of 192 mm, with lesion length in the population averaging 114 mm. The distribution of lesion lengths was bell-shaped suggesting that several genes control this trait. The resistant checks (Karlo RR, S19–90 and Maple Donovan) ranked among the lines with the shortest lesions (1st, 3rd and 21st out of 127) while the highly susceptible check Nattosan had the second longest lesions (177 mm) and the two moderately susceptible checks (Williams 82 and OAC Bayfield) showed lesions slightly above the population average.

Marker distribution

To achieve extensive genome coverage, we re-analyzed previously obtained sequence data (940M single-end reads from ApeKI GBS libraries prepared from DNA of 530 elite Canadian soybean lines) using an improved SNP-calling pipeline (Fast-GBS) and a more recent version of the soybean reference genome. This yielded nearly 150K SNPs on the panel of 530 lines that included all lines of the association panel. We then used a catalog of 4.1M SNPs obtained from WGS of 102 lines, also included in the set of 530 lines, as a reference panel to impute genotypes at all the missing loci, thus resulting in a full dataset of 4.1M markers. Of these, 3.5M SNPs were polymorphic in our association panel (i.e., carried an alternate allele in at least for one of the 127 lines). After removal of SNPs mapping to scaffolds (49.7K SNPs), 3.4M SNPs mapped onto one of the 20 soybean chromosomes. Finally, we removed markers whose MAF was lower than 0.05, thus resulting in a final catalog of 1,493,960 SNPs with which we performed the GWAM analysis.

Population structure and kinship

To characterize population structure, we pruned SNPs in high LD (r2 ≥ 0.9; windows of 50 SNPs), and the remaining 84,708 SNPs were used in fastSTRUCTURE. The results suggested that the panel was composed of between three and six subpopulations. Based on these two results, we chose to perform the ensuing analysis using K = 6 and the corresponding plot is shown in Figure 2. To further reduce confounding, we estimated the kinship matrix between lines of the association panel.

Genome-wide association mapping for SSR resistance

Marker-trait associations were estimated using the phenotypic data (mean lesion length) and the full set of filtered SNP markers (close to 1.5M markers). These were analyzed using an MLM (Q+K) and associations with p-values corresponding to an FDR < 0.1 were considered significant. In total, only two chromosomal regions were found to have at least one peak SNP exceeding this threshold (on Gm01 Gm15; Figure 3). As detailed in Table 1, the peak SNP on Gm01 was at position 5,594,597, showed a p-value of 5.08 x 10–5 and explained 32% of the phenotypic variation. As shown in Figure 4, accessions carrying the favorable allele (frequency = 0.38) at this locus showed shorter lesions although some accessions still exhibited lesions averaging over 100 mm. A second associated region was found on chromosome Gm15 with a single significantly associated marker at position 13,665,369 (p-value = 9.76 x 10–5; FDR = 0.04) and explained 15% of the variation. Accessions fixed for the minor allele (frequency = 0.32) had lesions that were 15 mm shorter than those fixed for the major allele.

Validation experiment

As the association on Gm15 had already been validated in previous work, we focused here on validating the candidate region for SSR resistance on Gm01. To do this, we used a population of F6:8 lines derived from a cross between OAC Bayfield (S) and Maple Donovan (R). These parents were contrasted for the peak marker on Gm01 as well for SSR resistance; Maple Donovan carries the resistance allele and developed lesions 78.3 mm shorter than those exhibited by OAC Bayfield. The parents were used as checks in the validation trial in addition to 47 recombinant inbred lines (RILs) selected as a validation population. For each line, four plants were genotyped using a CAPS marker developed to assay the QTL on Gm01. Among the 47 RILs, 21 were homozygous for the resistance allele while 26 were fixed for the susceptible allele. These RILs, along with the parents, were then evaluated for SSR resistance in two greenhouse trials. The contrast in lesion length between the parents was still evident (60 mm). Among the RILs, the average lesion length was 63 mm and ranged from 16 to 107 mm. Interestingly, almost all genotypes fixed for the resistance allele developed lesions under the average, ranging between 16 and 78 mm (average of 40 mm), whereas lines homozygous for the susceptible allele averaged 83 mm with lesion length extending from 51 to 107 mm. The phenotypic contrast between the two genotypic classes (43 mm) (Figure 5) was significant (p = 0.007).

Candidate genes near the peak SNP on Gm01

The LD (r2) was estimated between all marker pairs between positions 5.4 and 5.8 Mb to investigate the genomic landscape around the associated region on chromosome Gm01. Then, LD blocks were constructed using an r2 threshold of 0.8. The LD plot (Supplementary figure 1) showed that the associated region was a part of a single LD block spanning nearly 276 kb (from 5,528,405 to 5,805,120). This block contained 21 genes in high LD with the peak SNP on Gm01 (Table 2). Gene ontology analysis showed different predicted molecular functions for these genes such as DNA Binding Transcription Factor activity (Glyma.01g047900, Glyma.01g049100, and Glyma.01g049400) carbohydrate binding (Glyma.01g048500) and oxidoreductase activity (Glyma.01g048800, Glyma.01g048900). Based on the currently available annotation, none of these genes could readily be associated with the resistance phenotype.

A total of 1142 SNP were found in the same LD block and 20 of these were predicted to moderately impact gene function; none were predicted to have a high impact. These SNPs resided in four genes and might alter the proteine structure and change its effectiveness. Eight SNPs were found in Glyma.01g048000, a gene of unknown function, also carrying the peak marker.

Discussion

Number of markers

In this study, we used nearly 1.5 M high-quality SNPs. Based on the study of three chromosomal segments, Hyten et al. (2007) had estimated that the number of SNPs needed to obtain sufficient coverage (at r2 0.8) in elite material was somewhere between 9,600 and 29,000 markers. More recently, a genome-wide estimation of LD in telomeric and paracentromeric regions led Bastian et al. (2014) to conclude that around 60,000 SNPs would be required to achieve extensive coverage. None of the previous GWAM studies investigating SSR resistance achieved such a marker coverage. In two earlier studies using SNP markers, genotyping was conducted with the GoldenGate assay, achieving between 858 (Mamidi et al., 2011) and 1,142 SNPs (Hao et al., 2012). Later GWAM studies used GBS-derived SNP catalogs of 7,864 (Bastien et al., 2014), 8,397 (Iquira et al., 2015) and 11,811 (Wei et al., 2017), while another study achieved higher coverage with 25,179 SNPs obtained using SLAF-seq (Zhao et al., 2015). The highest coverage prior to this work had been achieved recently using the SoySNP50K BeadChip, giving 31,600 and 35,708 SNPs for two AM populations (Wen et al., 2018). Here, combining GBS data, WGS data and imputation for missing genotypes, we significantly increased SNP coverage, ensuring for the first time a marker coverage likely conferring exhaustive genome-wide coverage in our association panel. It was previously shown that such imputed data are of high accuracy, with 96.4% of the imputed missing genotypes being in agreement with those obtained at loci in common with the SoySNP50K array (Torkamaneh et al., 2017a). Given the recent increase in availability of WGS data for numerous collections of soybean germplasm, we feel the two-step genotyping approach used in this work will tremendously enhance our ability to perform genome-wide scans with full marker coverage.

New QTL on Gm01

Within our association mapping panel of 127 Canadian soybean lines, we identified 16 SNPs significantly associated with SSR resistance falling in small segments of only two chromosomes (Gm01 and Gm15). The first associated region on Gm01, extending over 200 kb, was novel as it did not overlap with any of the previously reported QTLs. In fact, this region didn’t carry any GBS-derived SNPs or those supported by the SoySNP50K BeadChip.

The second QTL was discovered on Gm15 and is the exact same marker-trait association reported previously by Bastien et al. (2014). These results were expected given that we exploited essentially the same association panel (except for three lines that were removed) and the same phenotypic data. In this previous work, three other QTLs had been reported (on Gm01, Gm19 and Gm20), but none of them were rediscovered in this work. All three of these associations were characterized by FDR values (ranging between 0.04 and 0.09) that were much higher than the FDR for the QTL on Gm15 (0.01). To investigate possible causes for this apparent lack of reproducibility, we compared the SNP genotypes at all associated SNPs with their corresponding genotypes in our more recent data set resulting from large-scale imputation. Surprisingly, for SNP on chromosome Gm15, 98% of the genotype calls were identical whereas for the other associated SNPs (on Gm01, Gm19 and Gm20), a much lower level of concordance (60 to 82%) was observed for genotype calls among the lines of the association panel. Based on the work of Torkamaneh and Belzile (2015), we believe that the number of SNPs discovered in the work of Bastien et al. (2014) proved insufficient to adequately capture haplotypes in this association panel and had led to inaccurate imputation of missing genotypic data.

QTL validation

Numerous previously reported QTLs for SSR resistance were discovered in different genetic backgrounds, environments, and various sampling populations. Before considering the use of such QTLs for marker-assisted selection, a validation step is highly recommended. The QTL on Gm15 was previously validated by Bastien et al. (2014) where resistance alleles reduced lesion length by 12.3 and 17.6 mm in two populations of RILs segregating for associated marker. Such validation has rarely been reported in other studies. Here, we wanted to similarly validate the novel marker-trait association discovered on Gm01 by evaluating SSR partial (or quantitative) resistance (lesion length) in RILs carrying contrasting alleles at this locus. Results showed that lines homozygous for the resistance allele (inherited from Maple Donovan) developed lesions 43 mm shorter than those homozygous for susceptible allele (derived from OAC Bayfield). This phenotypic contrast was highly significant and explained a substantial portion of the phenotypic variation between the parents. Also, the estimated allelic effect at this locus in the segregating RILs was more important than the one measured in the association panel, suggesting that fewer QTLs conditioned SSR resistance in the biparental population than in the association panel. Taken together, these data suggest that the region of Gm01 is associated with SSR resistance.

QTL detection efficiency and potential uses in genomic selection

The two regions discovered in this work, on Gm01 and Gm15, explained 29 and 15%, respectively, of the phenotypic variation for this trait in our association panel. In previous work, QTLs were reported to explain between 3 to 23% of the variance (Bastien et al., 2014; Iquira et al., 2015; Moellers et al., 2017; Zhao et al., 2015; Wei et al., 2017 and Wen et al., 2018). Thus, the QTL reported here on Gm01 constitutes the most impactful QTL reported for this trait to this date. As the favorable allele was present in a minority (38%) of the elite lines in our panel, there is considerable scope to improve SSR resistance by selecting for this partial resistance QTL in future breeding efforts.

Despite this, it seems that a large proportion of variation for this trait remains uncaptured in our association panel despite the extensive marker coverage. Thus, there could be other undiscovered, potentially impactful QTLs contributing to SSR in this population. In such panels, however, it will always be difficult to capture rare alleles (MAF < 5%) potentially contributing to resistance as the associated markers will have been filtered out. Expanding this work to larger and more diverse panels could help in discovering additional QTLs. Also, the remaining uncaptured heritability may be explained by numerous small-effect QTLs that will be difficult to discover (given the difficulties associated with phenotypic measurements) and would be of limited use in breeding programs as the benefits of a small increase in SSR resistance might not justify the cost of marker-assisted selection for such minor QTLs.

The genomic landscape around the QTL on Gm01

When we investigated the associated region on Gm01, we found 21 genes in high LD with the peak marker. Of these, six genes belonged to the zinc finger superfamily of transcription factors. Genes encoding transcription factors belonging to the zinc finger superfamily have shown increased transcript levels in response to Sclerotinia sclerotiorum infection in oilseed rape (Jean et al., 2018) and Brassica napus (Zhao et al., 2007b). Another gene, Glyma.01g049500, encoded a calcium-binding EF-hand family protein. Such gene families have been shown to be involved in soybean response to SSR (Rojas, 2014). A comparative study of gene expression between OAC Salem (partially resistant) and OAC Shire (susceptible) after 0, 3- and 5-days post-inoculation with S. sclerotiorum identified 712 up-regulated genes from a total of 2,316 that were differentially expressed. From these only 14 were always overexpressed including Glyma02g33660 (calcium- binding EF-hand family protein). The most significant associated markers where found in or in the vicinity of one gene Glyma.01g048000 also carrying eight SNPs predicted to impact its function. Although, no information about these gene where found.

Conclusion

We conclude that the knowledge that comes out from this study will promote the addressing of the SSR challenges under sustainable agricultural practices. We believe that this work is another step forward in rendering GWAM data more applicable in plant breeding. We expect that genetic region identified in this study will become a key tool in soybean breeding programs for SSR and enable geneticists and molecular biologists to identify causal resistance genes for SSR in near future.

Declarations

References

Arahana VS, Graef GL, Specht JE, Steadman JR, Eskridge KM. Identification of QTLs for Resistance to Sclerotinia sclerotiorum in Soybean. Crop Sci. 2001; 41:180–188.

Barrett JC, Fry B, Maller J, Daly MJ. Haploview: analysis and visualization of LD and haplotype maps. Bioinformatics. 2005; Doi: 10.1093/bioinformatics/bth457.

Bastien M, Huynh TT, Giroux G, Iquira E, Rioux S, Belzile F. A reproducible assay for measuring partial resistance to Sclerotinia sclerotiorum in soybean. Can J Plant Sci. 2012; Doi:10.4141/cjps2011–101.

Bastien M, Sonah H, Belzile F. Genome Wide Association Mapping of Sclerotinia sclerotiorum Resistance in Soybean with a Genotyping-by-Sequencing Approach. Plant Genome. 2014; Doi:10.3835/plantgenome2013.10.0030.

Benjamini Y, Yekutieli D. Quantitative trait loci analysis using the false discovery rate. Genetics. 2005; 171:783–90.

Boland GJ, Hall R. Index of plant hosts of Sclerotinia sclerotiorum. Can J Plant Pathol. 1994; Doi:10.1080/07060669409500766.

Bolton MD, Thomma BPHJ, Nelson BD. Sclerotinia sclerotiorum (Lib.) de Bary: Biology and molecular traits of a cosmopolitan pathogen. Mol Plant Pathol. 2006; Doi:10.1111/j.1364–3703.2005.00316. x.

Browning S, Browning B. Rapid and accurate haplotype phasing and missing-data inference for whole genome association studies by use of localized haplotype clustering. Am J Hum Genet. 2007; 81:1084–1097.

Danecek P, Auton A, Abecasis G, Cornelis AA, Banks E, et al. The Variant Call Format and VCFtools. Bioinformatics. 2011; https://doi.org/10.1093/bioinformatics/btr330.

Elshire RJ, Glaubitz JC, Sun Q, Poland JA, Kawamoto, K, et al. A robust, simple Genotyping-by-Sequencing (GBS) approach for high diversity species. PLoS ONE. 2011; Doi: 10.1371/journal.pone.0019379.

Guo X, Wang D, Gordon S, Helliwell E, Smith T, Berry S, St Martin S, Dorrance A. Genetic mapping of QTLs underlying partial resistance to Sclerotinia sclerotiorum in Soybean PI 391589A and PI 391589B. Crop Sci. 2008; 48:1129–1139.

Hao D, Cheng H, Yin Z, Cui S, Zhang D, et al. Identification of single nucleotide polymorphisms and haplotypes associated with yield and yield components in soybean (Glycine max) landraces across multiple environments. Theor Appl Genet. 2012; Doi:10.1007/s00122–011–1719–0.

Huynh TT, Bastien M, Iquira E, Turcotte P, Belzile F. Identification of QTLs Associated with Partial Resistance to White Mold in Soybean Using Field-Based Inoculation. Crop Sci. 2010; Doi:10.2135/cropsci2009.06.0311.

Hyten DL, Choi IY, Song QJ, Shoemaker RC, Nelson RL, et al. Highly variable patterns of linkage disequilibrium in multiple soybean populations. Genetics. 2007; Doi:10.1534/genetics.106.069740.

Iquira E, Humira S, Belzile F. Association mapping of QTLs for Sclerotinia stem rot resistance in a collection of soybean plant introductions using a genotyping by sequencing (GBS) approach. BMC Plant Biol. 2015; https://doi.org/10.1186/s12870–014–0408-y.

Kandel R, Chen C, Grau C, Dorrance A, Liu J, Wang Y, Wang D. Soybean Resistance to White Mold: Evaluation of Soybean Germplasm Under Different Conditions and Validation of QTL. Front Plant Sci. 9:505. 2018; https://doi.org/10.3389/fpls.2018.00505.

Kim HS, Diers BW. Inheritance of partial resistance to Sclerotinia stem rot in soybean. Crop Sci. 2000; 40:55–61.

Koenning SR, Wrather J. A. Suppression of Soybean Yield Potential in the Continental United States by Plant Diseases from 2006 to 2009. Plant Health Prog. 2010; doi:10.1094/PHP–2010–1122–01-RS.

Mamidi S, Chikara S, Goos RJ, Hyten DL, Annam D, et al. Genomewide association analysis identifies candidate genes associated with iron deficiency chlorosis in soybean. Plant Genome. 2011; Doi:10.3835/plantgenome2011.04.0011.

Mila AL, Yang XB. Effects of fluctuating soil temperature and water potential on sclerotia germination and apothecial production of Sclerotinia sclerotiorum. Plant Dis. 2008; 92:78–82.

Moellers TC, Singh A, Zhang J, Brungardt J, Kabbage M, Mueller DS, et al. Main and epistatic loci studies in soybean for Sclerotinia sclerotiorum resistance reveal multiple modes of resistance in multi-environments. Sci Rep. 2017; Doi:10.1038/s41598–017–03695–9.

Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MAR, et al. PLINK: a toolset for whole-genome association and population-based linkage analysis. American Journal of Human Genetics, 81. 2007; http://pngu.mgh.harvard.edu/purcell/plink/

Raj A, Stephens M, Pritchard J. K. fastSTRUCTURE: Variational Inference of Population Structure in Large SNP Data Sets. Genetics. 197:573–589. 2014; Doi:10.1534/genetics.114.164350.

Rojas V. Physiological, anatomical and molecular characterization of partial resistance against Sclerotinia sclerotiorum in soybean. 2014; Doctoral dissertation, The University of Guelph, Guelph, Ontario, Canada.

Schmutz J, Cannon SB, Schlueter J, Ma J, Mitros T, Nelson W, et al. Genome sequence of the palaeopolyploid soybean. Nature. 2010; 463:178–183.

Tang YLiu XWang JLi MWang Q, et al. GAPIT Version 2: An Enhanced Integrated Tool for Genomic Association and Prediction. Plant Genome. 2016; doi:10.3835/plantgenome2015.11.0120.

Torkamaneh D, Belzile F. Scanning and filling: ultra-dense SNP genotyping combining genotyping-By-sequencing, SNP array and whole genome resequencing data. PLoS ONE. 2015; 10, e0131533.

Torkamaneh D, Laroche J, Bastien M, Abed A, Belzile F. Fast-GBS: a new pipeline for the efficient and highly accurate calling of SNPs from genotyping-by-sequencing data. BMC Bioinformatics. 2017b; DOI 10.1186/s12859–016–1431–9.

Torkamaneh D, Laroche J, Tardivel A, Donoughue L O’, Cober E, et al. Comprehensive Description of Genome-Wide Nucleotide and Structural Variation in Short-Season Soya Bean. Plant Biotechnol. J. 2017a; https://doi.org/10.1111/pbi.12825.

Vuong T, Diers B, Harman G. Identification of QTL for Resistance to Sclerotinia Stem Rot in Soybean Plant Introduction 194639. Crop Sci. 2008; 48:2209–2214.

Wei W, Mesquita ACO, Figueiró A de A, Wu X., Manjunatha S, Wickland DP, et al. Genome-wide association mapping of resistance to a Brazilian isolate of Sclerotinia sclerotiorum in soybean genotypes mostly from Brazil. BMC Genomics. 2017; Doi:10.1111/tpj.12810.

Wen Z, Tan R, Zhang S, Collins PJ, Yuan J, Du W, et al. Integrating GWAS and gene expression data for functional characterization of resistance to white mould in soya bean. Plant Biotechnol J. 2018; Doi:10.1111/pbi.12918.

Wrather JA, Koenning SR. Estimates of disease effects on soybean yields in the United States 2003 to 2005. J Nematol. 2006; 38:173–180.

Zeng W, Kirk W, Hao J. Field management of Sclerotinia stem rot of soybean using biological control agents. Biol Control. 2012; https://doi.org/10.1016/j.biocontrol.2011.09.012.

Zhao X, Han Y, Li Y, Liu D, Sun M, et al. Loci and candidate gene identification for resistance to Sclerotinia sclerotiorum in soybean (Glycine max L. Merr.) via association and linkage maps. Plant J. 2015; Doi:10.1111/pbi.12918.

Tables

Table 1. Characteristics of the markers most highly associated (peak SNPs) with lesion length

Chromosome (Gm)

Position

p-value

FDR

MAF

R2

Allelic effect (mm)

01

5,594,597

5.08 x 10-5

0.02

0.38

0.32

29

15

13,665,369

9.76 x 10-5

0.04

0.32

0.15

15

 

  

Table 2. Genes in LD with the peak SNP at position 5,594,597 on Gm01 and their predicted function.

Gene Name

Start

End

Domain/Function

Glyma.01g047500

5528742

5532260

BED zinc finger, hAT family dimerisation domain

Glyma.01g047600

5545562

5545978

F-box family protein

Glyma.01g047700

5549606

5554165

RING/FYVE/PHD zinc finger superfamily protein

Glyma.01g047800

5559528

5560353

unknown

Glyma.01g047900

5585129

5586802

ZF-HD protein dimerisation region

Glyma.01g048000

5589570

5601987

Protein of unknown function

Glyma.01g048100

5605305

5607452

Pentatricopeptide repeat (PPR) superfamily protein

Glyma.01g048200

5615672

5619252

Ribonuclease, T2 family

Glyma.01g048300

5620361

5624143

unknown

Glyma.01g048400

5627602

5631187

Ribonuclease, T2 family

Glyma.01g048500

5659553

5665638

Galactosyltransferases

Glyma.01g048600

5668640

5670346

unknown

Glyma.01g048700

5680982

5681671

RING ZINC FINGER PROTEIN

Glyma.01g048800

5681940

5684496

Glucose-methanol-choline (GMC) oxidoreductase family protein

Glyma.01g048900

5695082

5697378

Glucose-methanol-choline (GMC) oxidoreductase family protein

Glyma.01g049000

5710637

5711191

glutathione S-transferase THETA 3

Glyma.01g049100

5717569

5719993

Homeodomain-like superfamily protein

Glyma.01g049200

5736667

5737599

Protein of unknown function

Glyma.01g049300

5740729

5741566

F-box family protein

Glyma.01g049400

5763611

5767687

zinc finger superfamily protein

Glyma.01g049500

5802428

5803205

Calcium-binding EF-hand family protein)