Identifying candidate genes for Phytophthora capsici resistance in pepper (Capsicum annuum) via genotyping-by-sequencing-based QTL mapping and genome-wide association study

Phytophthora capsici (Leon.) is a globally prevalent, devastating oomycete pathogen that causes root rot in pepper (Capsicum annuum). Several studies have identified quantitative trait loci (QTL) underlying resistance to P. capsici root rot (PcRR). However, breeding for pepper cultivars resistant to PcRR remains challenging due to the complexity of PcRR resistance. Here, we combined traditional QTL mapping with GWAS to broaden our understanding of PcRR resistance in pepper. Three major-effect loci (5.1, 5.2, and 5.3) conferring broad-spectrum resistance to three isolates of P. capsici were mapped to pepper chromosome P5. In addition, QTLs with epistatic interactions and minor effects specific to isolate and environment were detected on other chromosomes. GWAS detected 117 significant SNPs across the genome associated with PcRR resistance, including SNPs on chromosomes P5, P7, and P11 that colocalized with the QTLs identified here and in previous studies. Clusters of candidate nucleotide-binding site-leucine-rich repeat (NBS-LRR) and receptor-like kinase (RLK) genes were predicted within the QTL and GWAS regions; such genes often function in disease resistance. These candidate genes lay the foundation for the molecular dissection of PcRR resistance. SNP markers associated with QTLs for PcRR resistance will be useful for marker-assisted breeding and genomic selection in pepper breeding.


Results
evaluation of pcRR resistance in eCRILs and the GWAs core collection. We evaluated 188 F 7:8 ECRILs, parental lines, and additional resistant and susceptible controls in two environments (E1 and E2) for resistance to three P. capsici isolates (Fig. 1). The three isolates showed a significant difference in terms of symptom development and virulence. Susceptible controls (ECW30R, Tean, and Geumsugangsan) infected with the highly virulent isolate KPC-7, showed disease symptoms, including wilting and the appearance of water-soaked lesions at stem collars within 72 h of inoculation, and the plants completely wilted and died at 7 to 10 days post inoculation (DPI). By contrast, in plants infected with the moderately virulent isolate JHAI1-7, we noticed disease symptoms on the susceptible controls at 5 to 7 DPI, and the plants completely wilted and died at two weeks post inoculation. Susceptible controls infected with the isolate MY-1, with low virulence, symptoms were observed at 10 to 14 DPI, but complete plant death was not detected up to three weeks post inoculation. No visual symptoms of collar rot were observed on resistant controls cv. CM334 and Mohanjilju (commercial F 1 hybrid, Syngenta Korea) after inoculation with all three isolates, even at three weeks post inoculation (Fig. 1a). However, PI201234 was partially resistant to isolates KPC-7 and JHAI1-7 and completely resistant to isolate MY-1 (Fig. 1a).
We evaluated PcRR resistance in the ECRIL and GWAS core collection and scored disease severity on a scale of 1-4 at three weeks post inoculation in their respective environments, as described by 36 (Supplementary Fig. S1). Ten plants from each ECRIL were grown and the resistance level of each plant scored individually. We calculated the frequency distribution based on the mean resistance values of each RIL (Fig. 1b,c). The frequency distribution of disease severity showed a normal distribution for KPC-7 and JHAI-7 in both environments, although more ECRILs showed susceptibility to KPC-7 in E1 than E2 (Fig. 1b,c). By contrast, the distribution of disease severity of the RILs skewed towards resistance to the low-virulence isolate MY-1, and no ECRILs reached the disease rating of 4 during the experiment in either environment, except for two ECRILs in environment E2 (Fig. 1b,c and Supplementary Fig. S2). Correlation analysis indicated a significant, positive correlation between the two environments for PcRR resistance; however, the correlations between different isolates varied ( Supplementary  Fig. S2).
We evaluated the GWAS accessions for PcRR resistance against the same three isolates as the ECRILs in two different environments (EA and E1). Among GWAS accessions, CM334 was completely resistant, whereas PI201234, YCM334, and nine other accessions were partially resistant, with slight disease symptoms at the stem collars upon infection with isolates KPC-7 and JHAI-7 (Fig. 1a,d and Supplementary Table S2). Most of these PcRR-resistant accessions include known resistance sources (Supplementary Table S2). The isolate JHAI1-7 showed slightly higher virulence compared to KPC-7 in a few lines of GWAS core collection, but a higher mortality rate was observed for KPC-7 ( Fig. 1d and Supplementary Fig. S3). The disease response to the low-virulence isolate MY-1 was less severe, with 52 accessions showing no PcRR symptoms (Fig. 1d). Correlation analysis revealed significant, positive correlations between the KPC-7 and JHAI1-7 isolates and both environments for PcRR resistance (Supplementary Fig. S3). Highly virulent isolate KPC-7 for ECRILs in E2 and moderately virulent isolate JHAI1-7 for the GWAS accessions showed a minimum coefficient of variation (CV) of 38.55 and 28.09%, respectively (Table 1). Weak isolate MY-1 showed the highest (CV 57.0 to 62.02%) disease score in both screened populations in all environments ( Table 1). The broad-sense heritability (H 2 %) was found to be high for all PcRR isolates. In E1, isolate KPC-7 showed the highest H 2 (82.48%) among ECRILs, whereas isolate JAHI-7 showed an H 2 of 87.40% and 84.16% for ECRILs in E2 and the GWAS accessions, respectively (Table 1). sNp discovery through GBs and bin map construction for eCRILs. We performed genotyping of the ECRILs using GBS after EcoRI and MseI digestion. Sequencing of the prepared GBS libraries of 188 ECRILs and parental lines resulted in 525.62 million raw reads (data not shown). After trimming the raw reads, we obtained an average of approximately 4.4 million reads per sample (Table 2). After aligning the reads to the CM334 reference genome v1.6, 66,405 SNPs were detected (Fig. 2a). The SNPs were distributed more densely at the ends of the chromosomes than in the middle regions (Fig. 2a). After removing more than 90% of missing data and filtering unequally distributed SNPs, we obtained 13,021 high-quality SNPs that were equally distributed across the genome, which were used for bin linkage map construction (Tables 2 and 3 and Supplementary Fig. S4). To impute missing data and genotyping errors, we employed a sliding window approach 23 . Recombination breakpoints were determined by sliding 25 SNPs consecutively as one window. As a result, a high-density bin map consisting of 2,663 bins covering a total genetic distance of 1,428.8 cM was constructed (Table 3 and Supplementary  Fig. S4). The average genetic distance between bins was estimated to be 0.6 cM (Table 3). Among the 12 linkage groups, maximum and minimum genetic distances of 195.6 and 44.7 cM were obtained for chromosomes P1 and P8, respectively (Table 3). www.nature.com/scientificreports www.nature.com/scientificreports/ pcRR resistance QtLs in the eCRILs. We identified QTLs for PcRR resistance to three isolates, KPC-7, JHAI1-7, and MY-1, using 188 ECRILs and the high-density SNP linkage map in two environments, E1 and E2 ( Fig. 3 and Table 4). In total, 14 QTLs, including those with major and minor effects, those commonly detected, those specific to isolates, and those specific to environments were detected across the genome, explaining   Fig. 3 and Table 4). Three minor QTLs were detected in both environments, but only one QTL (on chromosome P1 at 77.61 cM [9.5-9.6 Mb]) was commonly detected in both environments ( Fig. 3 and Table 4 Fig. 3 and Table 4). Four minor QTLs on P8 and P11 were detected in both environments ( Fig. 3 and Table 4).  Table 3. Comparison of the physical length and genetic distance of bins for biparental QTL mapping.  Table 4). One major QTL on chromosome P5, E1My-5.3, was detected in E2 only at 34.61 cM, corresponding to 34.6-37 Mb, which explained 8.90% (LOD 4.38) of phenotypic variation (R 2 ). Two minor QTLs were detected in both environments, but only one minor QTL was commonly detected in both environments (on chromosome P11 at 66.41 cM [189-190 Mb]; Fig. 3 and Table 4). Notably, all QTLs associated with the MY-1 isolate showed lower R 2 and LOD values when compared to the QTLs detected against the two other isolates used in this study ( Fig. 3 and Table 4).

Number of bins
Taken together, three major QTLs on P5 were commonly detected for all three isolates in two environments, whereas minor QTLs on different chromosomes were isolate or environment specific ( Fig. 3 and Table 4). Among the major QTLs, QTL5.2, which is located at 27.3-29.2 Mb, was commonly detected for all three isolates and was named according to its genomic position on chromosome P5. The major QTL at 18.7-19.5 Mb (common to JHAI1-7 and MY-1) was named QTL5.1. The major QTL at 34.6-37 Mb (common to KPC-7 and MY-1) was named QTL5. 3. Even though minor QTLs were commonly detected on chromosome 11 for all three isolates, they were located at different positions. We detected isolate-specific minor QTLs on other chromosomes, such as the minor QTLs Kpc1 (for KPC-7) on P1, Jha8.1, Jha8.2, Jha11.1, and Jha11.2 (for JHAI1-7) on P8 and P11, and My11 (for MY-1) on P11. epistatic interactions of pcRR resistance QtLs. We detected additive-by-additive epistatic interactions against highly virulent isolate KPC-7 between major QTLs QTL5. 2 Table S3). The detection of epistatic interactions between loci on chromosomes P5 and P11 points to multilocus epistatic control of PcRR resistance in pepper.
Genome-wide association study of pcRR resistance. We aligned the sequences derived from GBS to C. annuum cv. CM334 reference genome v.1.6. GBS genotyping of 352 accessions using two sets of libraries derived from double digestion with restriction enzyme pairs PstI/MseI and EcoRI/MseI generated a total of 6,126,403 of SNPs. The SNPs were evenly distributed using this combination of restriction enzymes. SNPs in the GWAS core collection were relatively uniformly distributed across the chromosomes (Fig. 2b). We filtered out the SNPs with minor allele frequency >0.05, SNP coverage >0.6, and inbreeding coefficient >0.8, resulting in 507,713 high-quality SNPs distributed evenly on chromosomes P1 to P12 (Fig. 2b and Table 2).
Of the SNPs association with resistance to PcRR against isolate KPC-7, 117 were significant SNPs, with -log10 (p) value > 7.0, explaining phenotypic variation (R 2 ) > 0.2 ( Fig. 4a and Supplementary Table S4). An association study revealed significant SNPs linked to PcRR resistance to KPC-7 on chromosomes P2 (1 SNP Table S4). The two significant SNPs were located at 28.2 to 28.6 Mb on chromosome P5. By contrast, GWAS against the low-virulence isolate MY-1 did not detect significant SNPs ( Supplementary Fig. S5b). Notably, only a few or no significant SNPs were associated with isolates JHAI1-7 and MY-1, perhaps due to environment and genotype interactions or to the low virulence levels of the isolates.
Importantly, significant SNPs associated with resistance to isolate KPC-7 on three regions of chromosomes P5, P7, and P11 colocalized with QTLs detected in this study ( Fig. 4a and Table 4). On chromosome P5, 49 SNPs colocalized with commonly detected major QTL region QTL5.2 (27.3-29.2 Mb) ( Fig. 4b and Table 4). One SNP detected on chromosome P7 at 25 Mb corresponded to the location of QTL E2Kpc-7 ( Fig. 4c and Table 4). Another region on chromosome P11 at 75.2 Mb corresponded to QTL E1Kpc-11.2 ( Fig. 4g and Table 4). The GWAS-SNPs at three regions on chromosomes P5, P10, and P11 colocalized with previously detected QTLs and linked markers (Fig. 4b,f). The GWAS-SNPs on chromosome P5 in a 2.5-Mb region at position 27.0-29.5 Mb colocalized with a QTL detected in earlier reports 4,9,10,21 (Fig. 4b). SNP marker S10_14299301 located at 14,299 kb on chromosome P10 colocalized with the previously identified QTL Ph051-10-2 (RAUDPC) for PcRR resistance 19 (Fig. 4f). Regions detected by GWAS with significant SNPs on chromosome P10 at 196 Mb (Fig. 4f) colocalized with the region containing PhR10, a race-specific PcRR resistance locus 17 . We identified 13 significant GWAS-SNPs linked to PcRR resistance at the end of the long arm (206-212 Mb) of chromosome P5 that were not detected in previous studies, emphasizing the involvement of this genomic region in PcRR resistance (Fig. 4a). We used the significant SNP S05_208290460, with a -log10(p) value of 9.72 identified from the GWAS-SNPs (Supplementary Table S4), to compare the resistance levels of the GWAS accessions. As shown in the box plots in Supplementary Fig. S6e, the homozygous genotype (S05_208290460) AA is associated with an increased resistance level against all three isolates compared to the alternative homozygous genotype GG.
We conducted extensive linkage disequilibrium (LD) analysis of the GWAS core collection (N = 352) based on all adjacent marker pairs within a chromosome or within a haplotype block (Supplementary Table S5). We identified 29,269 haplotypes, with an average of 2,439 per chromosome. The average haplotype was 81.2 kb in size, with an average of 16.3 SNPs per haplotype (Table 2 and Supplementary Table S5). Further haplotype block analysis carried out with significant SNPs on chromosome P5 revealed 4 blocks containing 64 significant SNPs (Fig. 4d and Supplementary Table S4). The three blocks on P5 colocalized with QTL 5.2 detected in this study and with previously identified QTLs 4,9,10,21 (Fig. 4b,d). Two haplotype blocks were obtained on each of the chromosomes P7, P10 and P11, including regions with significant SNPs (Fig. 4e,h,i). These results indicate that the identified haplotype regions are strongly associated with the PcRR resistance. predicted candidate genes for pcRR. We looked for candidate PcRR resistance genes in the 1 Mb upstream and downstream flanking regions of all significant GWAS-SNPs and bin markers linked to PcRR QTLs. Detailed descriptions of the predicted genes in the target regions are provided in Supplementary Table S6. A major and colocalized genomic region detected for isolate KPC-7 via GWAS and biparental QTL mapping (QTL5.2) associated with PcRR resistance encompassing a 2.5 Mb region (27.0-29.5 Mb) (Fig. 3 and 4a) contains 20 predicted genes, including genes for three receptor-like kinase (RLK) domain-containing proteins, two nucleotide-binding site-leucine-rich repeat (NBS-LRR) domain-containing proteins, and one SAR8.2 precursor protein known to be associated with disease resistance (Table 5). Two NBS-LRR genes (PHT81215.1 and PHT81216.1) are located 618 and 319 kb upstream of the most significant SNP (S05_27703815), respectively. The most closely located RLK gene (PHT81221.1) was identified 144 kb upstream of a highly significant GWAS SNP. Two RLK genes (PHT81227.1 and PHT81229.1) were identified at 5.5 and 12.4 kb downstream of GWAS (SNP S05_28163627). One candidate gene, encoding the SAR8.2 precursor protein, was detected 148 kb upstream of GWAS SNP (S05_28825990) ( Table 5 and Supplementary Table S6). These genes were detected by both QTL mapping and GWAS analysis, making them strong candidate genes for PcRR resistance.
www.nature.com/scientificreports www.nature.com/scientificreports/ We used the highly significant GWAS SNP (S05_27703815), with -log10(p) value of 11.23 (Supplementary  Table S4), located close to the candidate genes to compare the resistance levels of the GWAS accessions. As shown in the box plots in Fig. 5c, the homozygous genotype (S05_27703815) AA was associated with increased resistance compared to the alternative homozygous genotype GG among the three PcRR isolates examined. We used www.nature.com/scientificreports www.nature.com/scientificreports/ the tightly linked bin marker EC5-bin40 of QTL segment QTL5.2 overlapping with candidate genes to compare the resistance levels of the ECRILs. As shown in the box plots, the homozygous resistant genotype BB is associated with enhanced resistance compared to the homozygous susceptible genotype AA against all three isolates in both E1 and E2 (Fig. 5a,b).
The major QTL region QTL5.1, spanning 0.7 Mb, which was commonly detected against two isolates, JHAI1-7 and MY-1, contains 15 predicted genes (Supplementary Table S6), including genes associated with disease resistance, such as genes for six RLK domain-containing proteins and four NBS-LRR domain-containing proteins (Table 5). Finally, we used the tightly linked bin marker EC5-bin27 of QTL segment QTL5.1 overlapping with candidate genes to compare the resistance levels of the ECRILs. As shown in the box plots, the homozygous resistant genotype BB is associated with enhanced resistance compared to the homozygous susceptible genotype AA against all three isolates in both E1 and E2 ( Supplementary Fig. S6a,b).
Finally, QTL5.3, a QTL region against isolates KPC-7 and MY-1, was detected at 34.6-37 Mb on chromosome P5, with 28 predicted genes. Among these, three predicted genes are linked to disease resistance proteins, including two NBS-LRR domain-containing proteins and one Kunitz trypsin inhibitor 2-like protein (Table 5). We used the tightly linked bin marker EC5-bin51 of QTL segment QTL5.3 overlapping with the candidate genes to compare the resistance levels of the ECRILs. As shown in the box plots, the homozygous resistant genotype BB is associated with enhanced resistance compared to homozygous susceptible genotype AA against all three isolates in both E1 and E2 ( Supplementary Fig. S6c,d).

Discussion
In this study, we combined traditional QTL mapping with GWAS to increase our understanding of PcRR resistance in pepper and detected common and race-specific QTLs for this trait. Although many QTL studies have been performed to investigate PcRR using biparental QTL mapping 4,9,10,20,21 , this is the first comprehensive study of PcRR resistance loci for multiple P. capsici isolates in multiple environments using biparental QTL mapping and GWAS. pcRR resistance QtLs. Due to the rapid genetic evolution and diversity of Phytophthora, there is an urgent need to breed pepper with resistance to multiple PcRR isolates 4,35 . Analysis of the differential responses of ECRILs and GWAS accessions against the three PcRR isolates provided us with valuable information about the range and frequency distribution of symptom development and virulence levels (Fig. 1). The mean disease scores of the resistant GWAS accessions were lower than those of the ECRILs (Supplementary Table S2  www.nature.com/scientificreports www.nature.com/scientificreports/ allelic diversity for PcRR resistance in the GWAS core collection. Furthermore, the high heritability of resistance (Table 1) points to the additive effects of genes, offering an excellent opportunity for improving this trait through selection.
Breeding for resistance to PcRR is quite challenging; even resistance incorporated into commercial pepper lines can be readily overcome by highly virulent P. capsici isolates 2 . CM334 is a well-known source of PcRR resistance. The inheritance of PcRR resistance in CM334 is complex 9,10 , and the genetic dissection of PcRR resistance would greatly facilitate the improvement of pepper cultivars. Therefore, in this study, we aimed to identify QTLs for different P. capsici isolates categorized based on virulence in two different environments. QTL mapping for PcRR resistance revealed significant QTLs, including three major-effect QTLs on chromosome P5 and several isolate-specific minor-effect QTLs on various chromosomes. Among the several QTLs identified for PcRR resistance, a major QTL on chromosome P5 has been consistently identified in several resistance sources irrespective of race, isolate, and level of virulence 4,10 . In agreement with this finding, two QTLs (QTL5.1 and QTL5.3) specific to single isolates and one QTL (QTL5.2) common to all three isolates of P. capsici were consistently detected on chromosome P5 in two environments. The location of QTL5.2 detected in this study coincides with that of the previously identified locus Pc5.1 4,10,21 . The dominant resistance locus CaPhyto, conferring race-specific resistance to PcRR, was previously mapped to a genetic interval near Pc5.1 in a 3.3-cM region between two SSR markers, ZL6726 and ZL6970, in C. annuum accession PI201234.
Several minor-effect QTLs have previously been reported on various chromosomes 4,11,[18][19][20]36 . The number of QTLs conferring PcRR resistance and their positions have varied among studies depending on the mapping populations used, disease screening methods, inoculum density, and isolate aggressiveness. In this study, we detected some minor QTLs with isolate specificity on chromosome P2 (Kpc1 for isolate KPC-7) and chromosome P8 (Jha8.1, Jha8.2, Jha11.1, and Jha11.2 for isolate JHAI1-7) and with no isolate specificity on chromosome P11 (Table 4). However, we could not compare the exact physical positions of these QTLs with previously detected QTLs due to the limited number of common markers and the lack of sequence information; instead, we used the closest markers for analysis.
In addition, epistatic interactions between PcRR resistance QTLs have been reported 11,36 . Consistent with these reports, we detected epistatic QTL interactions between resistance QTLs on chromosomes P5 and P11 for isolate KPC-7, on chromosomes P5 and P8 and chromosomes P5 and P11 for isolate JHAI1-7, and on chromosomes P5 and P11 for isolate MY-1 in both environments.
Validation of pcRR loci via GWAs. We took advantage of biparental QTL mapping together with GWAS to identify new PcRR resistance-related loci and to validate known QTL regions. This approach proved to be effective for the identification of PcRR resistance loci. We identified 117 significant SNPs related to PcRR isolate KPC-7 across the chromosomes. However, only SNPs identified on chromosomes P7, P5, and P11 colocalized with QTLs detected in this study (Fig. 4), perhaps because the ECRILs included only two parental alleles compared to the www.nature.com/scientificreports www.nature.com/scientificreports/ highly diverse GWAS accessions 33 . Thus, numerous biparental populations developed using different resistant genetic resources were likely captured among all resistant loci.
GWAS is a high-resolution, powerful technique, as it captures chronological recombination events that have accumulated inside an association panel. PcRR resistance loci have previously been mapped onto chromosomes P1, P3, P4, P6, P8, P9, P10, P11, and P12 4,11,[18][19][20]36 . To determine whether the GWAS-SNPs identified in this study are novel PcRR resistance SNPs, we compared their physical positions with the positions of previously reported QTLs. The GWAS-SNPs on chromosome P5 (at 27.0-29.5 Mb) coincided with major QTLs detected in previous studies 4, 10,21 . Two GWAS-SNP peaks on chromosome P10 at 14.2 and 196 Mb colocalized with the previously detected QTL Ph051-10-2 (RAUDPC) and the race-specific locus PhR10, respectively 17,19 . However, due to the limited availability of marker sequences, it was not possible to compare and validate all GWAS-SNPs with previously reported resistance loci linked to PcRR. The GWAS-SNPs associated with PcRR identified in this study could serve as novel resistance loci that could be used to incorporate resistance into ongoing pepper-breeding programs. potential candidate genes for pcRR resistance. In this study, we identified several classes of genes associated with disease resistance in various QTL regions. A cluster of NBS-LRR domain-containing R genes was identified on chromosome P5: four NBS-LRR genes from QTL5.1 (18.7-19.5 Mb), two from QTL5.2 (27.3-29.2 Mb), and two from QTL5.3 (34.6-37 Mb). These candidate genes share high sequence identity with RPP13-like NBS-LRR protein genes and represent potential candidates for PcRR resistance in pepper. RPP13 is a CC (coiled-coil)-NBS-LRR domain-containing R gene that confers resistance to the oomycete pathogen Peronospora parasitica in Arabidopsis thaliana 37 .
Plant genomes generally contain numerous NBS-LRR class R genes, which are often clustered on specific chromosomes due to tandem and segmental duplications [38][39][40][41] . Several dominant R genes, such as R1 from potato 42 , Ph-3 from tomato 43 , and RpsUN1, RpsUN2, and Rpg1-b from soybean 44,45 , belong to NBS-LRR domain-containing complex R gene clusters. R gene clusters could function as hotspots for the generation of novel R genes and enhance the possibility of structural and copy number variation through various DNA recombination and rearrangement mechanisms, including duplications, unequal crossing over, gene conversion, and diversifying selection, which could account for the gain or loss of resistance [45][46][47] . The candidate R gene clusters identified in this study could account for enhanced durability and resistance to PcRR in pepper. In addition to the NBS-LRR R gene clusters, we identified six RLKs from QTL5.1 (18.7-19.5 Mb) and three from QTL5.2 (27.3-29.2 Mb) on chromosome P5. RLKs are extracellular surface (also called pattern recognition) receptors that play crucial roles in plant defense-related processes, including both host and non-host defense responses 48 . Therefore, the RLK genes identified in this study represent strong candidates for PcRR resistance genes.
We also identified SAR8.2 in QTL5.2 (27.3-29.2 Mb), a systemic acquired resistance (SAR)-related gene known as CASAR82A. CASAR82A localizes to the phloem and epidermal cells of pepper leaf and stem tissues infected by Colletotrichum coccodes and P. capsici or treated with salicylic acid 48 . The involvement of the pepper SAR8.2 gene in pathogen infection and environmental stress responses and its use as a marker of abiotic elicitors 49 suggests it might also be an important candidate gene for PcRR resistance.
Kunitz trypsin inhibitor 2-like protein is involved in programmed cell death (PCD) in Arabidopsis. Arabidopsis serine protease (Kunitz trypsin) inhibitor (KTI1) appears to be a functional KTI when induced in bacteria and in planta. Atkti1 expression is induced at a late stage of infection in response to fungal and bacterial elicitors and to salicylic acid. RNAi silencing of AtKTI1 increased the rate of lesion development in leaf tissue after infiltration with the PCD-eliciting fungal toxin fumonisin B1 (FB1) or the avirulent bacterial pathogen Pseudomonas syringae pv. tomato DC3000 carrying avrB (Pst avrB) 50 . In this study, one gene encoding a Kunitz trypsin inhibitor 2-like protein was identified in a major QTL region on chromosome P5, QTL5.3 (34.6-37 Mb), suggesting it is a candidate gene for PcRR. However, to uncover the roles of intra-and extracellular immune receptors in disease resistance, further analyses of the candidate genes should be performed, including targeted knockdown of candidate genes and subsequent molecular analyses.
In this study, we combined the use of biparental QTL mapping and GWAS to identify QTLs and candidate genes for PcRR resistance in pepper. Based on the genetic map and phenotypic data, 14 significant QTLs were identified with R 2 values ranging from 3.73 to 38.99% among the three PcRR isolates in their respective environments. Two of these QTLs on chromosome P5 showed major effects, with a combined R 2 > 60%. GWAS identified 117 highly significant SNPs associated with PcRR across the genome for isolates KPC-7 and JHAI1-7. The additional loci detected through GWAS might enhance the accuracy of selection for PcRR resistance. The use of SNP markers associated with the candidate genes or QTLs could accelerate MAS and genomic selection in pepper breeding programs for PcRR resistance. The predicted candidate genes for PcRR resistance lay the foundation for the molecular dissection of PcRR resistance.

Materials and Methods
plant materials and DNA extraction. The ECRIL mapping population comprised 188 F 7:8 ECRILs developed from a cross between resistant C. annuum line CM334 and susceptible line ECW30R. The single seed descent method was used for ECRIL development, with shuttle breeding. The population was abbreviated as "EC", based upon the first letter of parental accessions. The plants were grown in plastic pots in a greenhouse. GWAS was carried out on a core collection comprising 352 Capsicum accessions 32 . Both parents of the ECRIL population were included in the core collection. Genomic DNA was extracted from young leaf tissues of plants of the ECRILs and the association population at the seedling stage using the CTAB (cetyl trimethyl ammonium bromide) method as described previously 51 Table S1) used in this study were kindly provided by Dr. Choi (KRICT). Isolates MY-1, JHAI1-7, and KPC-7 were reported to have low, moderate, and high levels of virulence, respectively 10,52 . To prepare the inoculum, the P. capsici isolates were cultured on V8-agar medium and incubated at 27 °C for mycelium growth for five days. The plates were flooded with distilled water to harvest the mycelia. The spore concentration was quantified using a hemocytometer, and the spore suspension density was adjusted to approximately 5 × 10 4 sporangia per mL prior to inoculation. The ECRILs (188 F 7:8 plants) and GWAS panel of 352 core collection accessions were scored for disease resistance with 10 replications per line. The plants were grown in 50-cell plastic trays, and inoculation was carried out by drenching plants at the 4-6 true leaf stage with 5 mL of spore suspension at the base of the stem using a dispenser. The inoculated plants were maintained under protected conditions at 24-30 °C and frequently watered to facilitate disease establishment. C. annuum cultivars CM334, PI201234, and Mohanjilju (commercial F 1 hybrid, Syngenta Korea) were used as the resistant controls, whereas C. annuum cultivars ECW30R, Tean, and Geumsugangsan (commercial F 1 hybrid, Takii Korea) were used as the susceptible controls to compare the disease resistance levels. The disease indexes were scored at three weeks post inoculation based on a previously described disease scale of 1-4, where 1 = no visible symptoms, 2 = dark lesion visible at the base of the stem but surviving without wilting, 3 = wilting with a dark lesion at the base of the stem, and 4 = wilting and death of the whole plant 34,53 . Genotyping-by-sequencing. Genotyping-by-sequencing was performed for the ECRILs and GWAS core collection as described previously 23 . Briefly, genomic DNA was diluted, adjusted to a concentration of 50 ng/ µL per ECRIL, and digested with the restriction enzymes EcoRI and MseI using a SBG 100-Kit v2.0 (Keygene N.V., Wageningen, the Netherlands). For the GWAS core collection population, GBS libraries were constructed manually using restriction enzymes PstI/MseI and EcoRI/MseI as described previously 23,54 . Following the ligation of selective adaptors, the libraries were amplified using selective primers containing "GA" for the ECRILs and "TA" for the GWAS core collection population. Amplified libraries consisting of 188 ECRILs and two replicates of susceptible parents were pooled into a single tube. Libraries consisting of the GWAS core collection populations were pooled into five tubes. The pooled libraries were sequenced on the Illumina HiSeq. 2000 platform at Macrogen (Macrogen, Inc. Seoul, Korea). Trimming and quality control of the GBS raw data were performed using a CLC Genomics Workbench v6.5 (Qiagen, Aarhus, Denmark) with a minimum read length of 80 bp and a minimum quality score of Q20. Filtered reads were aligned to the C. annuum cv. CM334 reference genome (chromosome v1.6, http://peppergenome.snu.ac.kr) using the Burrows-Wheeler Aligner (BWA) 55 . Filtering and SNP calling were performed using the Genome Analysis Toolkit (GATK) Unified Genotyper version 3.3-0. SNPs in the ECRILs were filtered with QUAL value > 20 and minimum read depth of 3. SNPs in the GWAS core collection population were filtered with minor allele frequency (MAF) > 0.03, calling rate >0.6, and inbreeding coefficient (F) > 0.8.

Construction of bin map and QtL analysis.
SNPs with unequal segregation and > 10% missing data were excluded from genetic linkage map construction. For genetic linkage map construction, bins were treated as genetic markers. The sliding window approach was used to impute the missing data and to detect recombination breakpoints as described previously 23 . To assign genetic positions to the bins, arranged bins were mapped with a LOD (logarithm of the odds) threshold of 3.0 and a distance threshold of 30 cM using CarthaGene software. The Kosambi mapping function was used to infer genetic distances between markers in centimorgans (cM). QTL analysis was performed using Composite Interval Mapping (CIM) with Windows QTL Cartographer 2.5 56 . The 1000-permutation test (P < 0.05) was performed to determine the LOD threshold for significance of each QTL. Explanations of the percentage of phenotypic variance [R 2 (%)] and additive effects for each QTL were obtained using the same software. Epistatic effects and the interactions of major QTLs were examined using the Multiple Interval Mapping (MIM) option under the command "Scan through QTL mapping results file" and refined further through "Testing for existing QTLs" with a Bayesian Information Criterion (BIC-X) model using default options. The physical locations of the QTLs were compared with the genetic locations of bins linked to the significant QTLs. The C. annuum cv. CM334 reference genome (chromosome v1.6, http://peppergenome.snu.ac.kr) was used to compare the positions of significant SNPs in the QTLs and GWAS core collection identified in this study and with previously reported linked markers. The marker sequences of previous studies were retrieved using https://solgenomics.net/search/markers and BLAST analysis of the reference genomes. If the sequence information for a particular linked marker was not publically accessible, the closest marker was used.
Genome-wide association analyses for pcRR resistance. The 507,713 filtered SNPs discovered from the 352 core collection population were used for GWAS mapping. Population structure (PCA and kinship matrixes) and genome-wide association (GWAS) were estimated with the Compressed Mixed Linear Model (CMLM) using the R package Genomic Association and Prediction Integrated Tool (GAPIT) 57 with default parameters. All probabilities generated in the association runs were transformed by -log 10 P (FDR p-value < 0.05). Scores for an individual chromosome were then inspected in Manhattan plots to determine whether the SNPs reached the significance threshold. The -log 10 p-values of SNPs from GWAS were adjusted by 58 multiple test correction, and the adjusted cutoff value for accepting thresholds was set to -log 10 (P) value ≥ 7.0. The haplotype www.nature.com/scientificreports www.nature.com/scientificreports/ blocks of the GWAS core collection were estimated using PLINK v1.9 with the following settings:-no-parentsallow-no-sex-blocks-max-kb 2000-blocks-inform-frac 0.9-blocks-strong-highci 0.85-blocks-recomb-highci 0.8 59 . Haplotype block analysis for the chromosomes identified with significant SNPs was performed using Haploview 4.2 software 60 . The significant SNP data sets from GAPIT were used to calculate pair-wise LD between SNPs. A 1% cutoff was used. Haplotype blocks in the region were defined with the confidence interval method 61 .
Identification of candidate genes. To identify candidate genes associated with PcRR resistance, the positions of highly significant and linked bin markers within QTLs on the genetic map were compared with their physical positions on the pepper reference genome (v1.6), and 1-Mb upstream and downstream sequences were mined for genes associated with PcRR disease resistance. Similarly, 1-Mb sequence information for the highly significant SNPs associated with the PcRR resistance trait from the GWAS peak was used for BLAST (blastn) analysis against the C. annuum cv. CM334 reference genome database (chromosome v1.6, http://peppergenome. snu.ac.kr) to identify candidate PcRR resistance genes. statistical analysis to estimate genetic variability against pcRR. The mean disease scores of ten plants per ECRIL and GWAS core collection accessions were used to calculate the frequency distribution, coefficient of variation (CV%), range, and genetic variability parameters for three different isolates in their respective environments. Broad-sense heritability (H 2 bs ) was calculated with the formula (suggested previously 62 ) = × , where, σ g 2 is the genetic variance and σ p 2 is the phenotypic variance. The Pearson correlation was calculated using Rstudio.