A Combinatorial Approach of Biparental QTL Mapping and Genome-Wide Association Analysis Identifies Candidate Genes for Phytophthora Blight Resistance in Sesame

Phytophthora blight, caused by pathogen Phytophthora nicotianae, is responsible for a huge reduction in sesame (Sesamum indicum L.) crop yields. In this study, we utilized a combinatorial approach involving biparental QTL mapping and genome-wide association (GWAS) analysis to identify genes associated with Phytophthora blight resistance in sesame. Evaluation of resistant of the parental varieties (Goenbaek, Osan and Milsung) and the RILs of both the populations in greenhouse conditions suggested the qualitative nature of the trait.. The genetic map comprised thirteen LGs covering a total map length of 887.49 cM with an average inter-marker distance of 4.69 cM. Significant QTLs explaining phenotypic variation in the range of 2.25% to 69.24% were identified on chromosomes 10 and 13 (Chr10 and Chr13). A resistance locus detected on Chr10 was found to be highly significant. The association of this locus to PBR was also identified through BSA and single marker analysis in Goenbaek × Milsung cross and through genome-wide association mapping of 87 sesame accessions. The GWAS analysis identified 44 SNP loci significantly associated with Phytophthora disease-resistant traits on Chr10. Further, the haplotype block analysis conducted in order to find whether the SNPs associated with resistance in this study showed that the SNPs are in high LD with the resistance QTL. We obtained a total of 68 candidate genes, which included a number of defense-related R genes. One of the genes, SIN_1019016 (At1g58390) showed high expression in the resistant parent. The results from this study would be highly useful in identifying genetic and molecular factors associated with Phytophthora blight resistance in sesame.


Introduction
from the GWAS sesame panel to develop two bi-parental mapping populations, which were 145 later genotyped using GBS to identify the QTL associated with Phytophthora blight resistance.  (Table S2). For the Goenbaek and 156 Milsung F2 plants, the segregation ratio also followed the expected segregation ratio for a bulks (Table S3). Out of 54 primers used, six primer pairs were polymorphic in two susceptible cvrs. Milsung and Osan and resistant cvr. Goenbaek parental lines. The similar polymorphic 179 pattern showed in Population-I five primer pairs 180 and SiSSM-85023 and SiSSM-86906), in Population-II six primer pairs 181 SiSSM-84958, SiSSM-85023, SiSSM85354, SiSSM-85570 and SiSSM-86906) were 182 polymorphic (Table 1; Table S3). In total seven SSR markers, 183 SiSSM-84996, SiSSM85354, discriminate between resistance and susceptible genotypes and their resistant and susceptible 185 bulks in both the populations (Table S3). Among the surveyed primers in BSA analysis, six 186 primers were able to distinguish parents and F2 bulks derived from the Goenbaek and Osan, 187 Goenbaek and Milsung varieties ( Figure S1) correlated with phenotypic disease evaluation 188 scores. As a result, these primers were scored in order to know the correlation with the  (Table 1).

Single marker analysis to identify markers associated with Phytophthora blight in 192
Population-II 193 Single marker analysis was performed to identify markers associated with Phytophthora blight 194 in Population-II. From the BSA, we identified six polymorphic SSR markers in which two from Goenbaek (Table 4). A single QTL qPhn-13_1 detected on Chr13 showed a negative 239 additive effect, meaning allele derived from susceptible parent 'Osan' (Table 4). The QTL 240 qPhn-10_2 showed the highest LOD score and explained 69.24% of total PVE. The physical 241 positions of the QTLs detected in Goenbaek ⅹ Osan 90 F5:7 RILs were compared using the S. 242 indicum 'ver.2' reference genome. Two genomic regions on  Mbp) and 13 243 (6.74-10.27 Mbp) were associated with resistance to Phytophthora blight, and the QTLs located 244 on Chr10 were considered the most significant. Linkage analyses revealed Phn-10 was linked 245 to polymorphic markers from the BSA F2 polymorphism survey ( Figure S1). Goenbaek was 246 the same resistance parent in the two RIL populations, the two loci were the same locus. 248 To validate the QTLs detected from the biparental populations, a GWAS study was 249 performed using sesame accessions using the greenhouse-based layer test assay (Table S1).

250
Using CMLM in the GAPIT R package PCA+K Bonferroni correction has captured 251 associations for three P. nicotianae isolates with major effects on Chr10 (Table 6; Figure 4).

269
The raw data of flanking sequences (600 bp) obtained from GBS sequencing of associated 270 SNPs were BLASTed for identifying candidate genes. The analysis of physical locations of 271 associated SNP markers revealed that all markers were located in intergenic regions. Five SNPs 272 (S10-14860929, S10-14860930, S10_14860891 and S10_14860890, S10_14841897) were 273 associated with P. nicotianae isolates KACC48120, KACC48121 and No2040 were found to 274 be in intron part within major resistance gene SIN_1018978 and SIN_1018979 275 (LOC105171865, probable disease resistance protein RF45) on sesame Chr10 using BLAST 276 (Table S5). However, we noticed a single gene with 5 SNPs markers (S10_14860930, 277 S10_14860929, S10_14860891, S10_14860890 and S10_14841897) found inside the gene as 278 a candidate R gene and all are in intergenic region (  (Table S5; Figure 4 a,b,c). Haplotype block possessing marker(s) identified were 285 significantly associated with a Phytophthora blight resistance. 286 We compared the linkage mapping and association mapping results. The QTL intervals qPhn-287 10_1, qPhn-10_2, and qPhn-10_3 in Population-I (Table 5), QTLs qPhn-10_1, qPhn-10_2 in 288 Population-II were overlapped with GWAS QTL interval on Chr10 ( Figure 5). These SSR 289 markers not only exhibited the highest PVE value and had significant marker-trait correlation 290 (P < 0.0001), but were also appeared in haplotype groups of No2040 (28 genes), KACC48120 291 (38 genes) and KACC48121 (56 genes) genomic region ( Figure 4). We survey a number of 292 candidate genes between SSR markers, 539,477 bp to 293 14,765,975 bp) ( Table S4). The genomic region between 44 associated SNP markers is 0.81 294 Mbp (14.31 -15.13 Mbp) and is not large (Table 6). This QTL was considered major associated 295 with four SSR markers mapped at 0.0 -9.0 cM interval in a strong agreement with a physical 296 distance of significant SNP markers associated with GWAS confirming tightly linked to disease 297 resistance. Therefore, we focused on the genes in this interval. The corresponding genomic sequences of the region after haplotype block analysis were extracted. An annotation analysis 299 showed that Chr10 from 14.31-15.13 Mb contains 80 genes (Table S5). These genes starting 300 from the letter "SIN" were annotated according to the S.indicum v.2 reference annotation 301 database (http://www.sesame-bioinfo.org/Sinbase2.0/).

302
Identifying candidate genes associated with SSR markers 303 All the SNPs associated with Phytophthora blight were further evaluated by aligning the tag 304 sequence having the associated SNP against all the known sequences to infer position and the 305 potential candidate genes carrying a true causative polymorphism. Each SNP in intergenic 306 position was considered for potential functional annotation based on the actual proximity of 307 nearby located genes (Table S5; Table 6). The combined analysis involved BSA, linkage and 308 association mapping approaches for identifying the possible R genes as strong candidates in 309 this overlapped region (Table S4 and

349
Molecular tagging of the gene(s) or QTL conferring Phytophthora blight resistance has not 350 been previously reported in sesame. It is known that screening for disease resistance genes 351 using molecular markers offers the additional advantage of permitting selection for resistance 352 in the absence of the pathogen. Therefore, the availability of PCR-based markers will facilitate 353 breeding efforts to efficiently develop Phytophthora blight-resistant sesame cultivars. In the 354 present study, we employed an integrated approach involving QTL mapping using biparental 355 mapping RIL populations, and GWAS to identify candidate genes governing Phytophthora 356 blight resistance in sesame. In this study, using QTL mapping and single-marker analysis 357 approaches, we identified a novel gene, designated as phn-10 on Chr10. This QTL explained 358 approximately 2.93-69.24% of the PVE, providing additional evidence of one minor and major 359 gene coding for P. nicotianae resistance. The QTL was found to colocalize with the markers 360 from BSA and with the QTL mapped using the Goenbaek ⅹ Osan (Population-I) and

361
Goenbaek ⅹ Milsung (Population-II) populations. Several significant SSR markers 362 associated with this major QTL found to cosegregates with Phytophthora blight resistance in 363 the RILs used in this study. Thus, our study shows the success of GBS-generated SNP markers 364 in gene-tagging in sesame.

365
From the GWAS results, a total of 44 SNPs were associated with Phytophthora blight resistance.

366
Among 44 associated SNPs, 20 were common in all three isolates, 17 were co-localized in two 367 isolate strains KACC48120 and No2040, three were common between two strains 368 (KACC48121 and KACC48120) and 5, 5, and 3 were only specific to KACC48121, No2040, 369 and KACC48120, respectively. It is possible that these common SNPs are highly correlated 370 alleles conferring R-genes carrying plausible candidate genes could be referenced in QTL fine 371 mapping of RILs with a favored combination for crossing with a hypersensitive parent to 372 generate a fine-mapping population. However, associated SNPs to Phytophthora resistance for 373 mining putative candidate disease resistance genes were located in intergenic regions not in 374 exon regions. It is reported that the number of SNPs called from the GBS approach may not be 375 large enough for GWAS and that SNP density in genic regions is much lower than that in 376 intergenic regions (Ahn et al., 2018;Kim et al., 2014). Although, our germplasm panel 377 subjected to GWAS identified sufficient and common SNPs associated with Phytophthora blight resistance to the three main strains (isolates KACC48121, KACC48120, No2040) which SIN_1019016 gene, other candidate genes from the QTL regions we identified in this study 410 may also have significant role in Phytophthora blight disease resistance. Such genes when 411 introgressed into the elite cultivars have been shown to impart resistance to respective 412 pathogens. Therefore, precise nature, role, and organization of the R-genes have been a subject 413 of research for long time. In-depth characterization of the genetic regions shown to be 414 associated with the disease resistance to three races may reveal further information related to 415 the molecular mechanisms involved in Phytophthora blight disease resistance in sesame.  In this study, two RIL populations were developed by crossing a sesame cultivar Goenbaek as 430 a resistant parent, with Osan (P2), and Milsung (P2), both being the susceptible parents. The to produce F2:3 family seeds, which were advanced to the F5 generation by single-seed descent.

515
Depending on the ratio of SNP/InDel reads to mapped reads variant type are classified into 516 three categories; homozygous SNP/InDel for more than 90%, heterozygous SNP/InDel for 517 more than 40% and less than 60%, and the rest of them as etc. To control the quality of markers, 518 missing proportion (MSP) <0.3 and minor allele frequency (MAF) > 0.05 were selected.

519
Phenotype-genotype association analysis 520 The association panel used in this study (n = 87 sesame accessions) (Table S1) were subjected 521 to GBS after filtering steps, a total of 8883 SNPs covering the whole genome with minor allele 522 frequency >0.03 were used for GWA analysis. Before fitting the model, each marker was coded 523 with the value 0 used for the reference allele and the value 1 for the alternative allele. probabilities generated in the association runs were transformed by ─log10 P (FDR p value < 0.05). Scores for an individual chromosome were then inspected in Manhattan plots to 530 determine whether the SNPs reached the significance threshold. After multiple test Bonferroni 531 correction significance was defined at a uniform threshold of p < 5.54ⅹ10 -7 (─log10 (p) > 6).

532
The haplotype blocks and the identified significant SNPs for the chromosome were performed 533 using the software PLINK v.1.9 and Haploview 4.2 ( Barrett et al., 2005;Chang et al., 2015).   (Table S6)  The qRT-PCR thermocycling program was 50°C for 2 min, 95°C for 10 min, 40 cycles of 10s 608 at 95°C, 10s at 60°C and 30s at 72°C. Then the melting curve temperature profile was obtained 609 by heating to 95°C for 15s, cooling to 60°C for 1min, and slowly heating to 95°C for 15s with     The target interval on Chr10. The blue bar represents the gene which was identified to exhibit 816 a significant correlation between phenotype variation of Phytophthora disease resistance and 817 gene expression levels.