Insights into the genetic architecture of Phytophthora capsici root rot resistance in chile pepper (Capsicum spp.) from multi-locus genome-wide association study

Abstract

Background

Phytophthora root rot, a major constraint in chile pepper production worldwide, is caused by the soil-borne oomycete, Phytophthora capsici. This study aimed to detect significant regions in the Capsicum genome linked to Phytophthora root rot resistance using a panel consisting of 157 Capsicum spp. genotypes. Multi-locus genome wide association study (GWAS) was conducted using single nucleotide polymorphism (SNP) markers derived from genotyping-by-sequencing (GBS). Individual plants were separately inoculated with P. capsici isolates, ‘PWB-185’, ‘PWB-186’, and ‘6347’, at the 4–8 leaf stage and were scored for disease symptoms up to 14-days post-inoculation. Disease scores were used to calculate disease parameters including disease severity index percentage, percent of resistant plants, area under disease progress curve, and estimated marginal means for each genotype.

Results

Most of the genotypes displayed root rot symptoms, whereas five accessions were completely resistant to all the isolates and displayed no symptoms of infection. A total of 55,117 SNP markers derived from GBS were used to perform multi-locus GWAS which identified 330 significant SNP markers associated with disease resistance. Of these, 56 SNP markers distributed across all the 12 chromosomes were common across the isolates, indicating association with more durable resistance. Candidate genes including nucleotide-binding site leucine-rich repeat (NBS-LRR), systemic acquired resistance (SAR8.2), and receptor-like kinase (RLKs), were identified within 0.5 Mb of the associated markers.

Conclusions

Results will be used to improve resistance to Phytophthora root rot in chile pepper by the development of Kompetitive allele-specific markers (KASP®) for marker validation, genomewide selection, and marker-assisted breeding.

Background

Phytophthora capsici, a soil-borne oomycete first isolated in 1918 from chile pepper plants at New Mexico Agricultural Experiment Station field plots, and formally described in 1922 [1], remains a major problem in chile production worldwide [2, 3]. The pathogen infects multiple parts of the chile pepper plant, including stems, roots, fruits, and leaves at all stages of growth, and can cause significant production losses up to 100% [2, 4]. The most destructive disease symptom of P. capsici is root rot. Damping off due to root rot can cause death within two to five days in seedlings after infection; root rot can result in stunted growth, wilting, and ultimately plant death in older plants within approximately two weeks [2, 5, 6].

Various management practices including crop rotation, water management, fumigation, soil solarization, and fungicide application have been used but none of these strategies have been completely effective in controlling the disease [2, 7]. The most efficient, eco-friendly, and sustainable approach to managing this pathogen is through the utilization of resistant cultivars [8, 9]. The complexity of the Phytophthora-Capsicum pathosystem is a major challenge in disease resistance breeding. Different resistant genes are known to be required for resistance to different physiological races and syndromes [8, 10,11,12]. Moreover, several inheritance models including single dominant [6] and two-gene or multiple genes with dominant and recessive epistasis [2, 8] have been reported for P. capsici resistance, which presents an additional challenge. A dominant inhibitor for P. capsici resistance (Ipcr) gene which inhibits polygenic host resistance has also been identified in chile peppers [13]. Many resistant accessions have been found with different levels of resistance to various isolates. Among the different sources identified, the Mexican landrace, ‘Criollo de Morelos-334’ (‘CM-334’), has been known to have the greatest degree of resistance to all the physiological races of P. capsici and has been widely utilized in various breeding programs [9, 14]. To this date, however, a resistant commercial cultivar with broad-spectrum resistance comparable to the original resistant sources has not been developed [9].

Genetic studies have mainly focused on traditional biparental linkage analyses for various complex traits in chile pepper such as capsaicinoid content [15,16,17], yield-related traits [18, 19], fruit-related traits [20,21,22], plant architecture [21, 23], cucumber mosaic virus (CMV; Cucumovirus) resistance [24, 25], and P. capsici resistance [8, 9, 26,27,28]. This mapping approach depends on the genetic variation found in the parents and can only identify broad genomic regions, making it challenging to detect the specific candidate genes responsible [16]. Genomewide association study (GWAS) is an efficient tool that can be utilized to identify genetic loci linked with various complex traits using natural populations [29] and can be used to complement biparental mapping. Previously, a total of 117 and 30 significant SNP markers were identified across the entire Capsicum genome using single-locus GWAS for P. capsici resistance using 352 [9] and 342 accessions [30].

The current study utilized multi-locus GWAS models which allowed simultaneous calculation and testing of marker effects, improving its power over single-locus GWAS models [31, 32]. Various crop species have already employed multi-locus GWAS models including wheat (Triticum aestivum) [33,34,35], corn (Zea mays) [36,37,38], soybean (Glycine max) [39,40,41], rice (Oryza sativa) [42, 43], and barley (Hordeum vulgare) [44, 45] to understand the genetics of different complex traits. Although multi-locus GWAS models have been applied recently in chile pepper to identify genomic regions associated with yield, yield components, and other agronomic traits [32], there is currently no existing report on their implementation for the genetic dissection of resistance to P. capsici root rot in Capsicum.

The primary goal of this study was to understand the genetic basis of Phytophthora root rot resistance using a diverse population of Capsicum accessions. The specific objectives were to (1) screen a diverse panel of chile pepper for resistance to three different isolates (‘PWB-185’, ‘PWB-186’, and isolate ‘6347’) of P. capsici with varying degrees of virulence, (2) identify significant genomic regions linked to Phytophthora root rot resistance using multi-locus GWAS models, and (3) determine candidate genes associated with these significant loci. Results from this research will be relevant for molecular breeding for improving Phytophthora root rot resistance in chile pepper.

Results

Disease resistance screening

The diversity panel was screened for resistance to P. capsici isolates with varying degrees of virulence, viz. ‘PWB-185’ (moderate virulence), ‘PWB-186’ (low virulence), and ‘6347’ (high virulence). The resistant check (‘CM-334’) showed no symptoms of infection when screened against each of the isolates. The susceptible checks (‘Camelot’ and ‘NMCA 10399’) inoculated with ‘PWB-185’ started showing symptoms at 5-days post-inoculation (DPI) and were completely wilted and dead at 14-DPI. For isolate ‘6347’, the susceptible checks ‘Ninja’ and ‘NMCA 10399’, started displaying symptoms at 5-DPI and were wilted and dead by 10-DPI. The susceptible checks ‘Ninja’ and ‘NMCA 10399’ inoculated with ‘PWB-186’ started showing symptoms at 7–10-DPI, but the complete death of plants was not detected at 14-DPI, when the bioassay was completed.

The plants were scored every day on a 0–6 disease rating scale starting at 3-DPI [8, 46] up to 14-DPI (Fig. 1A). For both ‘PWB-185’ and ‘6347’ isolates, 90.4% and 84.1% accessions started showing symptoms at 5-DPI respectively; however, for ‘PWB-186’, only 24.2% of the accessions started showing symptoms at 5-DPI. Out of 157 accessions in the diversity panel, eight were completely resistant (disease score 0 at 14-DPI) to isolate ‘PWB-185’ and 35 accessions were not significantly different from the resistant check (‘CM-334’) (P ≥ 0.05). For ‘PWB-186’, 33 accessions in the diversity panel were completely resistant (disease score 0 at 14-DPI), and 130 accessions were not significantly different from the resistant check (P ≥ 0.05). Seven accessions were completely resistant to isolate ‘6347’ and 27 were not significantly different from the resistant check for this isolate (P ≥ 0.05). Overall, five C. annuum lines in the diversity panel, namely, ‘Chilhuacle Orange’, ‘Tipo Ancho’, ‘NMCA 10237’, ‘13C905-6’, and ‘Tipo Pasilla’ were completely resistant to the three isolates.

Fig. 1

A Disease progression in Capsicum accessions inoculated with three different isolates of Phytophthora capsici. B Distribution of the average disease scores at 14-DPI for isolate ‘PWB-185’, (C) ‘PWB-186’, and (D) ‘6347’. The test for normality overlaid with a red line

Significant differences were revealed using analysis of variance (ANOVA) between entries, replications, DPI, and blocks nested in replications for each isolate (P < 0.0001) (Table 1). Entry and DPI interactions were also tested in ANOVA and significant interactions were detected for isolates ‘PWB-185’ and ‘6347’ (P < 0.0001), but not for isolate ‘PWB-186’. Average disease score distribution was left-skewed (towards susceptibility) for ‘PWB-185’ and ‘6347’ and it was right-skewed (towards resistance) for ‘PWB-186’ (Shapiro–Wilk test, P < 0.0001; Fig. 1 B, C, D).

Table 1 Analysis of variance (ANOVA) for the three Phytophthora capsici isolates used in the study. Significant differences were detected in days post-inoculation (DPI), replications, blocks nested in replications, and entries (P < 0.0001). A significant interaction between entries and days post-inoculation was detected for isolates ‘PWB-185’ and ‘6347’

The highest average AUDPC values for ‘PWB-185’ were 46.2 (‘NMCA 11322’) and 45.9 (‘Mesilla’). Susceptible checks ‘Camelot’ and ‘NMCA 10399’ had average AUDPC values of 40.2 and 37.8, respectively. A total of 46 accessions had disease severity index percentage (DSI%) greater than 90%, and 11 accessions had resistant plants percentages greater than 90% against ‘PWB-185’. For the less virulent ‘PWB-186’ isolate, the highest average AUDPC values were 20.9 (‘20C239’) and 20.7 (‘Mallorca Paprika’), whereas the susceptible checks ‘Ninja’ and ‘NMCA 10399’ had average AUDPC values of 3.4 and 3.2, respectively. There were no accessions with DSI% greater than 90%, when screened using ‘PWB-186’, whereas 49 accessions had resistant plants percentages greater than 90%. The highest average AUDPC values for isolate ‘6347’ were 54.1 (‘NMCA 10406’) and 54.0 (‘Large Red Thick Cayenne’). The susceptible checks ‘Ninja’ and ‘NMCA 10399’ had average AUDPC values of 37.14 and 45.9, respectively. The number of accessions with DSI% and percentage of resistant plants greater than 90%, against ‘6347’ were 77 and 11, respectively. The broad-sense heritability (H²) for Phytophthora root rot resistance varied among the isolates and was related to the level of virulence. H² values for the traits were higher for the more virulent isolates, with an average of 0.78 and 0.90 for isolates ‘PWB-185’ and ‘6347’, respectively. For the ‘PWB-186’ isolate, H² across the traits had an average value of 0.25 (Table 2).

Table 2 Genotypic and phenotypic variance and broad sense heritability values for the disease parameters for the three Phytophthora capsici isolates

GBS-derived SNP markers

The original genotype file obtained after the sequence analysis had a total of 404,188 SNP markers. Excluding minor SNP states and SNPs with < 0.05 minor allele frequency resulted in 61,714 SNP markers. After filtering out the unmapped markers, 55,117 SNP loci with known chromosomal locations were identified and were used for multi-locus GWAS. The most predominant nucleotide was adenine (‘A’) (23.74%), followed by guanine (‘G’) (23.4%), thymine (‘T’) (23.2%), and cytosine (‘C’) (22.8%). The most common substitution types were ‘C/T’ transition (8,700; 14.09%) and ‘G/A’ transition (8,680; 14.06%). The most common transversion was ‘A/T’ (4,296; 6.96%) (Table 3). Chromosome 3 had the most markers with 7,388, followed by chromosomes 1 with 6,113 and 2 with 5,437. Chromosome 11 had the fewest markers with 3,217, followed by chromosomes 5 with 3,517 and 10 with 3,624. Chromosome 2 had the highest SNP density (33.2 SNPs/Mb), followed by chromosomes 8 (30.6 SNPs/Mb) and 3 (28.3 SNPs/Mb); whereas chromosome 11 had the least marker density of 14.6 SNPs/Mb, followed by chromosomes 9 (15.54 SNPs/Mb) and 5 (16.2 SNPs/Mb) (Fig. 2). The overall marker density for the whole genome was (20.8 SNPs/Mb).

Table 3 Allele summary for the SNP markers obtained for the chile pepper population

Fig. 2

Chromosome wise single nucleotide polymorphism (SNP) density plot with the number of SNP markers within a 1 Mb window. The horizontal axis represents the length of the chromosomes, and the bars represent chromosomes 1 to 12

Population structure and linkage disequilibrium

Principal component analysis (PCA) (Fig. 3A) revealed that the proportion of the variance defined by the principal components dropped after the fifth principal component (PC) and almost plateaued. A few of the accessions did not belong to any of the three major clusters, based on the PCA (Fig. 3B). Thus, to account for the structure of the whole population, five PCs were selected for the GWAS analysis. The PC1, PC2, PC3, PC4, and PC5 explained 17.8%, 11.1%, 4.9%, 2.9%, and 2.3% of the genotypic variance, respectively. Seven phylogenetic clusters were revealed using the neighbor-joining tree (Fig. 3C). Similarly, the Evanno criterion determined the optimum number of clusters, K = 7 (ΔK = 4,005.3), in the population based on STRUCTURE analysis, which suggested that the population could be divided into seven admixed subpopulations (Figs. 3D, E). Cluster 1 (Group A; N = 19) consisted of cultivars with a wide range of fruit types including bell, wax, cayenne, de arbol, poblano, serrano, pasilla, mirasol, jalapeño, and chilhuacle. Cluster 2 (Group B; N = 20) included breeding lines with New Mexican, cayenne, and paprika fruit types. Cluster 3 (Group C; N = 5) consisted of cultivars with jalapeño fruit types and cluster 4 (Group D; N = 13) included breeding lines (8) and cultivars (2) with cayenne fruit type and three cultivars with banana or wax pepper, jalapeño, and serrano fruit types. Cluster 5 (Group E; N = 4) consisted of interspecific hybrid accessions derived from crossing C. annuum with C. frutescens. Cluster 6 (Group F; N = 19) mostly included cultivars with New Mexican (including paprika types), a landrace, and a breeding line with jalapeño fruit type. Cluster 7 (Group G; N = 76) was an admixed of 64 breeding lines and 12 landraces with different fruit types including New Mexican, jalapeño, cayenne, poblano, and the rest included Aleppo, Asian, guajillo, mirasol, paprika, pasilla, pimiento, shipkas, and wax fruit types. The number of clusters, K = 3 also serves as an alternative number of optimum clusters as it displayed a high ΔK value of 752.47 compared to the other clusters (Fig. 3D).

Fig. 3

A Scree plot illustrating the gradual decline in the proportion of variance explained with the increase in the number of principal components. B PCA biplot constructed using the first two principal components demonstrating the genetic stratification in the diversity panel. C Neighbor-joining tree depicts the genetic relatedness of the 157 accessions used in this study. D The Evanno method determined the optimal number of clusters for the population, K = 7 with the highest value for ΔK. E Bar plots of the population structure coefficients of the 157 accessions. Seven colors represent seven subpopulations, respectively. An individual is represented by each vertical column and the percentage of the individual in the population is represented by each colored segment in each column. F Decay plot for linkage disequilibrium (LD). The half-decay value for LD (r.² = 0.236) is represented by the blue line and the non-linear regression curve corresponds to the red solid line. The green line corresponds to the distance in the base pair when the LD drops to half of its maximum value (~ 0.10 Mb)

A total of 1,944,143 intrachromosomal marker pairs were used for LD analysis, and 421,560 pairs (21.7%) were within significant LD (P < 0.05). The largest number of marker pairs within significant LD were found on chromosome 3 (61, 440), whereas chromosome 11 had the fewest (22, 047). A total of 174, 621 marker pairs showed complete LD (coefficient of LD, r² = 1.0). Chromosome 3 had the highest number of marker pairs in complete LD (26,041). The average distance between the pairs of markers in significant LD ranged from 0.63 (chromosome 2) to 1.64 Mb (chromosome 9) (Table S1). The r² values were observed to decrease over genetic distance (in Mb). Genome wide LD decay rate was approximately 0.10 Mb (Fig. 3F).

Significant SNPs for Phytophthora capsici root rot resistance

Overall, 330 SNP markers were detected to be linked to Phytophthora root rot resistance using different multi-locus GWAS models on various disease parameters (Table S2). For the ‘PWB-185’ isolate, 33, 31, 37, 28, and 34 SNPs were found to be significant (LOD score ≥ 3.0) by all six multi-locus GWAS models for AUDPC, disease score at 14-DPI, DSI%, estimated marginal means, and percent resistant plants, respectively. For the ‘PWB-186’ isolate, 30, 19, 24, 23, and 25 SNPs were significant for AUDPC, disease score at 14-DPI, DSI%, estimated marginal means, and percent resistant plants, respectively. In the case of isolate ‘6347’, 33, 39, 40, 37, and 44 significant SNPs were identified for AUDPC, disease score at 14-DPI, DSI%, estimated marginal means, and percent resistant plants, respectively. Lastly, for BLUPs, 33, 50, 55, 52, and 43 significant SNPs were found using all the multi-locus GWAS methods for AUDPC, disease score at 14-DPI, DSI%, estimated marginal means, and percent resistant plants, respectively. There was a total of 95, 76, 94, and 113 significant SNP markers for the four datasets for average disease score at 14-DPI for the four phenotypic datasets (Fig. 4). In the case of disease parameters, the highest number of total significant markers were identified for DSI% (136), followed by percentage of resistant plants (127). The total significant markers identified for disease score at 14DPI, estimated marginal means, and AUDPC were 123, 124, and 114 respectively. A total of 56 SNP markers were common across all four phenotypic datasets (Table S3).

Fig. 4

Manhattan plot for average disease scores at 14-DPI for Phytophthora capsici (A) isolate ‘PWB-185’ BLUEs, (B) isolate ‘PWB-186’ BLUEs, (C) isolate ‘6347’ BLUEs, and (D) BLUPs

Only a single SNP marker, S3_236038572, present on chromosome 3 was detected in all datasets and accounted for 0.74% to 9.92% of phenotypic variation (Fig. 4A). Out of the 56 significant markers, 11 markers were detected in at least three phenotypic datasets. The maximum number of common significant markers were identified on chromosome 12 (9 SNPs), followed by chromosomes 5 (7) and 6 (7). The maximum phenotypic variation among the common markers was explained by S6_705421 which was detected in three phenotypic datasets (‘PWB-185’, ‘6347’, and BLUP), ranging from 4.6% to 54.4%. Some of the significant SNPs had positive allele effects for AUDPC, DSI%, disease scores at 14-DPI, and estimated marginal means and had negative allele effects for percentage of resistant plants which indicates their potential association with susceptibility genes. For instance, the allele effect for SNP S5_65193705 on chromosome 5 was 21.8 for DSI% and -12.20 for percent resistant plants. A total of 38 out of the 56 common significant markers were found to be potentially linked to susceptibility genes. On the other hand, some SNPs had positive allele effects for percentage of resistant plants and negative allele effects for the rest of the disease parameters indicating their association with resistance genes. For example, SNP S9_167564458 had allele effects ranging from -25.1 for DSI% to 67.6 for percentage of resistant plants. A total of 18 significant SNP markers located on chromosomes 1, 3, 4, 5, 6, 9, 10, 11, and 12 were potentially associated with resistance genes. These potential associations with resistance or susceptibility genes can be confirmed with performing marker and QTL validation studies.

Candidate genes for Phytophthora root rot resistance

The 56 significant markers that were common across all phenotypic datasets were associated with a total of 829 candidate genes (Table S4). These genes are involved in a range of biological functions including cell cycle and division, oxidation–reduction, phosphorylation, defense response, DNA repair, protein serine/threonine kinase activity, metabolic processes, cell wall organization and biosynthesis, among others. A total of 15 genes that code for LRR receptor-like serine/threonine-protein kinase which are well-known for involvement in disease resistance were detected close to six SNPs (S2_137881264, S3_38944749, S4_212420204, S10_122768586, S10_189438330, and S12_5016287). Other genes associated with defense response including late blight resistance protein R1-A, β-1,3-glucanases, mitogen-activated protein kinase, elicitor peptide 6, receptor-like proteins, and ethylene responsive transcription factors, were also identified. These genes play a role in different mechanisms of action in defense response and could be good candidates for Phytophthora root rot resistance in Capsicum.

The SNP, S6_705421, which explained the highest phenotypic variation (R² = 54.4%), had 47 candidate genes in its proximity, which included four putative late blight resistance protein coding genes, one gene coding for leucine rich repeat domain and one gene that codes for NB-ARC domain-containing protein. A systemic acquired resistance (SAR 8.2) protein coding gene was detected 471 kb downstream to the significant SNP, S5_14665044, located on chromosome 5. Another SNP S2_137881264 on chromosome 2 with an LOD value up to 13.74 had a total of 52 genes within 0.5 Mb region upstream and downstream. It included one late blight resistance gene and nine putative LRR receptor-like serine/threonine-protein kinase coding genes.

Discussion

Phytophthora root rot is a highly damaging disease of chile pepper present throughout the globe and is very difficult to manage [2]. With modern molecular plant breeding tools, disease resistance QTL and genes can be detected, which ultimately could contribute to the development of resistant chile pepper varieties. GWAS is a tool used in plant breeding which provides higher mapping resolution than traditional biparental mapping approaches [47]. This study employed multi-locus GWAS approaches to better understand the genetics of resistance to P. capsici root rot in chile pepper. The virulence of P. capsici isolates used for screening exhibited variation, where some genotypes displayed resistance to specific races. A complex genetic relatedness was observed within the Capsicum population used for screening P. capsici root rot. Further, the complex genetic basis of disease resistance is reflected by the diversity of candidate genes identified in this study.

Variation in virulence for the P. capsici isolates and race-specific resistance

The genotypes that are virulent to the same host are referred to as a physiological race of a pathogen [48]. This study involved three isolates of P. capsici, each of which belonged to different physiological races with different levels of virulence. Around 82.8% of the accessions screened with isolate ‘PWB-186’ showed no significant difference in the average disease ratings compared to the resistant control (‘CM-334’). For isolates ‘PWB-185’ and ‘6347’, 22.3% and 17.2% of the accessions, respectively, did not exhibit significant differences from the resistant control. The susceptible checks for ‘PWB-186’ were not dead at 14-DPI, whereas for ‘PWB-185’ and ‘6347’, the susceptible checks were completely dead at 14-DPI and 10-DPI, respectively. Thus, the isolate ‘PWB-186’ was the least virulent among the three; isolates ‘PWB-185’ and ‘6347’ were highly virulent, with ‘6347’ being more virulent than ‘PWB-185’. Isolate ‘6347’ was also described as the most virulent isolate (designated as ‘PWB-175’) among 13 isolates used in a previous study [49].

Breeding pepper for resistance against P. capsici is a challenging task due to the constant emergence of new races that can overcome the existing host resistance [50, 51]. Different resistance genes are required for different physiological races within each disease phase [10]. Race-specific resistance was observed for some of the accessions in this study. For example, ‘Floral Gem’ was resistant to isolates ‘PWB-185’ and ‘PWB-186’ but was susceptible to isolate ‘6347’. In a previous study, this cultivar was screened against 10 isolates of P. capsici and was found to be resistant to eight isolates and susceptible to isolates ‘PWB-24’ and ‘PWB-73’ [52]. ‘Paladin’ was susceptible to all three isolates in this study, but it was resistant to three (‘PWB-53’, ‘PWB-75’, and ‘PWB-54’) out of eight isolates in a previous study [52]. ‘Paladin’ was also found to be resistant to isolates from New York and North Carolina, USA [53, 54]. In this study, ‘NuMex Vaquero’ was susceptible to all three isolates, but a previous study reported that this variety had resistance against Phytophthora root rot races 2 and 3 [55]. Five accessions, namely ‘Chilhuacle Orange’, ‘Tipo Ancho’, ‘NMCA10237’, ‘13C905-6’, and ‘Tipo Pasilla’ exhibited broad-spectrum resistance and were found to be completely resistant to the three isolates used in this study. These accessions can serve as potential resistant sources for future genomic breeding aimed at developing chile pepper with Phytophthora root rot resistance.

Population structure of the Capsicum spp. accessions

Genetic subpopulations for the genotypes were derived using Bayesian iterative algorithm approach used in STRUCTURE software. Based on the ad hoc Evanno criterion, the optimal number of clusters, K, was determined to be 7 (Fig. 3D), indicating that the population could be divided into seven subpopulations. This finding was also supported by the results of the phylogenetic analysis (Fig. 3C). Values for alpha, α (relative admixture levels between populations) for this K were also consistent among replications, and hence convergence of the implemented algorithm was verified (data not shown). About 97% of the panel used in this study consisted of C. annuum accessions and only four accessions were interspecific hybrids between C. annuum and C. frutescens. These four interspecific accessions formed a separate cluster based on the results from STRUCTURE, PCA, and Neighbour-joining analysis (Group E; Figs. 3B and 3C). The rest of the six clusters consisted of the C. annuum accessions. In previous studies, complex genetic relatedness has been found within the C. annuum group [32, 56], where accessions clustered into different groups. Taranto et al. [57] found clustering in C. annuum based on geographical locations and fruit-related traits. Similar to previous studies [58, 59], clusters of C. annuum accessions were detected on the basis of fruit or pod type in this study. The PCA revealed five groups, whereas the STRUCTURE analysis resulted in K = 7 optimum clusters; this disparity between PCA and STRUCTURE results has been previously noted in other studies [32, 58]. The optimum K derived in performing STRUCTURE analysis may not necessarily represent the ideal number of clusters since it is determined based on a pre-determined sampling method. Therefore, when interpreting the results of STRUCTURE, it is crucial to consider the biological significance of the optimum number of clusters [60].

The knowledge of physical distance of LD decay across the population of interest facilitates the selection of the number of markers for association studies. If LD decays over long distances, fewer markers are necessary for association mapping, resulting in lower resolution. Conversely, if LD decays rapidly within a short distance, a larger number of markers is needed for association mapping and the resolution is higher [61]. Rapid LD decay (~ 0.10 Mb) was observed for the population in this study. The marker density of one single nucleotide polymorphism (SNP) per 48.04 kb, on average, was thus sufficient for detecting the genomic loci associated with P. capsici root rot resistance. Our results were consistent with previous studies which detected rapid LD decay for C. annuum at 0.07 Mb [58] and 0.01 Mb [57]. While rapid LD decay in chile pepper would require an increased number of markers, it would also provide advantages to breeders in performing fine mapping of the genes of interest identified from association studies.

Significant markers associated with P. capsici root rot resistance

Performing multi-locus GWAS resulted in the detection of a total of 330 SNP markers with the highest being detected using the isolate ‘6347’ (114 SNPs). Fifty-six markers which were common in two or more isolates distributed across all 12 chromosomes were detected. Previous studies have also mapped the resistance loci on chromosomes 1–12 of chile pepper [8, 9, 27, 28, 62,63,64,65,66,67,68,69,70]. Chromosome 5 has been previously identified as a prominent region for the occurrence of resistance genes against P. capsici in Capsicum spp. [8, 9, 28, 71]. In this study, seven significant SNPs were identified on chromosome 5 with S5_14665044 explaining the maximum phenotypic variation of up to 29.6%. This marker was located ~ 0.2 Mb upstream of a major QTL peak detected using bulk segregant analysis that covered 18.8 cM and was delimited by two markers, CONTIG6473 and CONTIG1896 [72]. Phyto5SAR marker showing highest LOD at the major QTL detected by Liu et al. [72], was located ~ 1.4 Mb downstream to the S5_25470050 SNP identified in this study. GBS-derived-SNPs on chromosome 5 detected by Siddique et al. [9] at 27.0 to 29.5 Mb, were also located ~ 1.4 Mb downstream to the S5_25470050. Three SNP markers (S5_9345710, S5_14665044, and S5_25470050) detected using GWAS in this study coincide with the extended Pc5.1 (Ext-Pc5.1) region on chromosome 5 that is located between 8.35 Mb and 38.13 Mb [73]. These SNPs also coincide with the MQTL5.1 identified from meta-QTL analysis in chile pepper [71]. SNP S5_207353938 detected in this study was located ~ 0.19 Mb upstream to significant SNP S05_207549766 detected for isolate ‘KPC-7’ in another study with phenotypic variation ranging from 3.5% to 29.5% [9]. The findings of this study along with previous research provide further evidence supporting the significance of chromosome 5 for Phytophthora root rot resistance.

The SNP markers S6_170978694 and S12_225476747 identified in this study were ~ 0.53 Mb upstream and ~ 1.28 Mb downstream to significant SNP markers, S06_171517874 and S12_224191357 detected on chromosomes 6 and 12, respectively [9]. Another SNP marker, S10_197712360 identified in this study was ~ 1.69 Mb downstream to the SNP, S10_196014162 on chromosome 10 detected in a previous GWAS [9]. The significant marker, S8_133105638 identified in the current GWAS on chromosome 8 was ~ 0.45 Mb downstream to QTL.Pc8.1 detected using a RIL population obtained by crossing ‘CM-334’ (resistant accession) and ‘Early Jalapeno’ (susceptible accession) [8]. The previously identified QTL.Pc9 was also ~ 1.91 Mb upstream to the significant SNP, S9_237909149 detected on chromosome 9 in the present work [8]. In another GWAS, significant SNPs were identified at 9.14 Mb and 208.4 Mb on chromosomes 7 and 12 respectively, which were ~ 0.11 Mb and ~ 0.91 Mb upstream to the SNPs, S7_9249493 and S12_209274913 detected in this study [30]. Among the significant SNPs discussed above, only S10_197712360 and S12_209274913 markers were potentially associated with resistance genes based on their quantitative trait nucleotide (QTN) effect. The rest of the markers have potential association with the susceptibility genes for P. capsici root rot. A novel major SNP marker, S6_705421 was discovered in this study which accounted for up to 54.4% of the phenotypic variation. Most of the significant SNPs identified in this research represent potential novel resistance loci which can be targeted in breeding toward resistance. Though chromosome 5 has been known to be the major region with disease resistance loci against P. capsici, this study as well as previous studies have detected genomic regions on other chromosomes which demonstrates the complexity of Phytophthora root rot resistance in chile pepper [8, 9, 27, 30, 68, 74]. SNP validation will be performed using Kompetitive allele-specific markers (KASP®). KASP represents a homogeneous genotyping technology that relies on allele-specific oligo extension and fluorescence resonance energy transfer (FRET) to generate signals. This method is recognized for its cost-effectiveness and adaptability, making it a widely employed assay on a global scale [75, 76].

Candidate genes for Phytophthora root rot resistance

This study identified candidate genes linked to disease resistance. Nucleotide-binding site leucine-rich repeat (NBS-LRR) class resistance genes are widely recognized for conferring disease resistance. The LRR regions of these genes interact with extracellular ligands while cytoplasmic kinase domains facilitate signal transduction through phosphorylation [77]. These are essential for the proper function of pathogen recognition as they are responsible for identifying effectors that pathogens deliver to host cells during infection which leads to the activation of a strong resistance response known as effector-triggered immunity [78, 79]. A total of nine genes that encode LRR receptor-like serine/threonine-protein kinase protein was detected near SNP S2_137881264 on chromosome 2, two genes near S12_5016287, and one gene near S3_38944749, S4_212420204, S10_122768586, and S10_189438330. Most resistance (R) proteins that play a role in detecting pathogens and triggering innate immune responses possess a central nucleotide-binding domain (NBS found in NBS-LRR proteins) [80]. The NBS domain, also known as the NB-ARC domain [81] is comprised of three subdomains, NB, ARC1, and ARC2. The NB-ARC domain performs as a functional ATPase domain, and the activity of the R protein is regulated by its nucleotide-binding state [80]. A total of six genes that code for the NB-ARC domain were detected near SNP S6_705421, S6_6636269, S10_201485448, and S11_44892.

Receptor-like proteins and receptor-like kinases are extracellular surface receptors that play role as pattern recognition receptors in plants. They detect the molecular patterns derived from both microbes and the host, initiating the first stage of inducible defense for plant immunity, growth, and development [9, 82]. The genes encoding these proteins were detected in this study near SNPs S5_207353938, S8_141949349, S10_201485448, and S12_37348368. These candidate genes were also previously detected and are good candidates for Phytophthora root rot resistance in chile pepper [9]. SAR8.2 gene, also designated as CASAR82A, is a systemic acquired resistance (SAR)-related gene which plays role in pathogen infection, environmental stresses, and abiotic elicitors [83]. This gene was also identified previously in other P. capsici resistance study [9] and was also detected near SNP, S5_14665044 in this study, suggesting its involvement in Phytophthora root rot resistance in chile pepper.

Ethylene plays a vital role in numerous developmental processes and is known to be a critical mediator of abiotic and biotic stress responses in plants [84]. Previous studies have documented the involvement of ethylene responsive transcription factors in pathogen attack [85,86,87,88]. These factors were recently identified to be upregulated in plants infected with P. capsici [89]. In this study, the ethylene responsive transcription factor genes were identified near S4_206485683, S4_212420204, and S12_5016287 and might be important candidate genes involved in P. capsici resistance. Glucan endo-1,3-beta-glucosidases, also known as β-1,3-glucanases, are essential hydrolytic enzymes in pathogenesis-related groups of proteins and abundant in various plant species following pathogen infection [90]. These enzymes contribute significantly to the defense response by degrading the β-1,3-glucans that are present in the cell wall of microbes, especially fungus and generating signaling glucans which trigger the activation of global responses [91]. The gene coding for this enzyme detected in the proximity of SNPs S7_9249493, S12_1728396, S12_1902119, S12_5016287, and S12_209274913 might have a crucial role in degrading the β-1,3-glucans that are present in the cell wall of P. capsici and can be a strong candidate gene for resistance in chile pepper. A few late blight resistance genes were found to be associated with SNPs S2_137881264 (1), S6_705421 (4), S6_6636269 (1), S11_44892 (3), and S12_225476747 (1). These genes also belong to NB-ARC-domains containing resistance genes and encode a specific domain of potato resistance genes [92]. The plant cell wall serves as the initial physical barrier against pathogen intrusion and plays role in detecting external signals and stimulating defense response [93, 94]. The candidate genes associated with cell wall biosynthesis, organization, and modification like hydroxyproline-rich glycoproteins (S2_137881264, S8_140731657), UDP-glucuronic acid decarboxylase (S3_29930687, S5_9345710), and β-D-xylosidase (S1_65866380, S4_212420204), were identified which can play role in defense response against P. capsici. Previous reports have also identified the involvement of epigenetic mechanisms in the defense response to P. capsici infection [8, 71, 73]. In this study, candidate genes for DNA methylation, histone methylation, acetylation, and ubiquitination were detected, supporting the argument for the involvement of the epigenome in defense response. Overall, many candidate genes involved in different mechanisms were detected reflecting the complex genetic architecture of resistance to P. capsici root rot in chile pepper. Further molecular analysis needs to be conducted to better understand and confirm their roles in conferring root rot disease resistance in Capsicum.

Conclusions

Multi-locus GWAS was used to detect the genomic regions linked to Phytophthora root rot resistance. A total of 330 significant markers were identified using various multi-locus GWAS models on five disease parameters and three P. capsici isolates (‘PWB-185’, ‘PWB-186’, and ‘6347’). Overall, 56 significant markers distributed across all 12 chromosomes were found to be common across the isolates. Candidate genes including nucleotide-binding site leucine-rich repeat (NBS-LRR), receptor-like kinase (RLKs), system acquired resistance (SAR8.2), β-1,3-glucanases, late blight resistance genes, and genes involved in epigenetic mechanisms like histone methylation were identified. Results from this study will be further used for validation using Kompetitive allele-specific markers (KASP®) and can be used for marker-assisted selection and genomic selection for P. capsici resistance in Capsicum.

Materials and methods

Plant material

In this study, a diverse panel of Capsicum spp., comprising 157 accessions, including 44 cultivars, 99 breeding lines, and 14 landraces, was used (Table S5). The plants were grown in trays with 15 accessions and three checks in each tray at the Fabián García Research Center Greenhouse, New Mexico State University (NMSU), Las Cruces, NM (32°16′46.7″N, 106°46′24.7″W), under standard conditions for growing chile pepper [95]. Three replications were planted to screen for each isolate, where a replication consisted of up to four plants per accession. The checks used for the isolate ‘PWB-185’ were ‘CM-334’ (landrace; resistant check), ‘Camelot’ (bell pepper; susceptible check 1), and ‘NMCA 10399’ (jalapeño; susceptible check 2). For isolates ‘PWB-186’ and ‘6347’, all checks were the same as for ‘PWB-185’ isolate but ‘Ninja’ (bell pepper) was instead used in place of ‘Camelot’ as susceptible check 1 due to the lack of seeds.

Inoculum preparation and infection

Zoospore production was performed based on previous methods [8, 46, 96]. The isolates were grown on V8^® media petri plate. The V8 media (500 mL) consisted of calcium carbonate (1.5 g), technical agar (10 g), V8 juice (81.5 ml), distilled water (418 ml), rifampicin (0.2 g), ampicillin (0.25 g) and primarcin (100 µL). These plates were used to grow the P. capsici isolates for 5–7 days at room temperature. Then, under a fume hood, 5 cm plugs were excised from these petri dishes using a cork borer. A total of 5–7 of these plugs were placed in deep petri dishes (~ 200) with 25 mL of sterile distilled water. The formation of sporangia was induced by transferring the deep dishes to an incubator at 25 °C for 48 h with fluorescent lighting. Zoospore release was induced by transferring these dishes to a 4 °C cold chamber for 60 min and then placing them back in the incubator for 30 min. A cheesecloth was used to filter the contents of the deep petri dishes to remove agar plugs. The resulting solution was subjected to zoospore quantification using a hemocytometer (Hausser Scientific, PA, USA). The final concentration of zoospores was adjusted using sterile distilled water to ~ 2,000 zoospores per ml.

At the 4–8 leaf stage, inoculation was done by dispensing 5 mL of the final solution of 2,000 zoospores per mL to the soil in each plant cell resulting in ~ 10,000 zoospores for each plant. Disease scoring for the inoculated plants was done after every two or three days starting at three days post-inoculation (DPI) up to 14-DPI. A disease scoring system of 0–6 was used, where 0 = no visible disease symptoms, 1 = stem necrosis with no girdling, 2 = stem necrosis with girdling, 3 = stem necrosis with less than 50% defoliation, 4 = stem necrosis with greater than 50% defoliation, 5 = wilted, and 6 = dead (Fig. 5) [8, 46].

Fig. 5

Disease score ratings (0–6) for screening chile pepper plants for resistance to Phytophthora capsici root rot. Ratings: 0 = no visible disease symptoms, 1 = stem necrosis with no girdling, 2 = stem necrosis with girdling, 3 = stem necrosis with less than 50% defoliation, 4 = stem necrosis with greater than 50% defoliation, 5 = wilted, and 6 = dead [8, 46]

Phenotypic data collection and analysis

Initially, a total of 174, 166, and 164 accessions were screened for ‘PWB-185’, ‘PWB-186’, and ‘6347’ respectively, but only 157 accessions that were common among the three isolates were used for further analysis. The accessions were screened using a randomized complete block design (RCBD) with three replications for each isolate. Each block contained 15 accessions and three checks (one resistant and two susceptible checks), and each accession had a total of four plants. The four disease scores obtained for each replication of each accession were averaged to obtain a single disease score value. Screening for the three isolates was done separately and disease scores for each isolate were used to calculate the area under disease progress curve, disease severity index percentage (DSI%), estimated marginal means (least square means), and percent resistant plants. The ‘audpc’ function in ‘agricolae’ [97] package in R was used to calculate the area under disease progress curve (AUDPC) using the formula AUDPC=$\sum_{i=1}^{n-1}\frac{{y}_{i}+{y}_{i+1}}{2} {\text{x}}\left({t}_{i+1}- {t}_{i}\right),$ where, y_i = disease score of the plants at i^th observation (i = 1 being the first observation at time zero), t_i = number of days post infection at i^th observation, and n = number of observations.

The estimated marginal (least square) means were also calculated in R using the ‘emmeans’ function in the ‘emmeans’ package and disease score at 14-DPI [98]. Disease severity index percentage and percentage of resistant plants were calculated using the SAS® OnDemand for Academics (SAS Institute, Inc., Cary, NC). The formula from Chiang et al. [99] was used to calculate the disease severity index percentage for each line: DSI% = $\frac{[sum\left(class frequency\times score of rating class\right)]}{[\left(total number of plants\right)\times \left(maximal disease index\right)]}\times 100.$ According to the rating scale used in this study, the maximal disease index was a score of 6 for all genotypes. Percent resistance was calculated based on using the formula from Muhyi and Bosland [100]: $\frac{Number of plants with disease score 1 or lower}{Total number of plants}\times 100$.

The ‘aov’ function in R was used to perform analysis of variance (ANOVA) considering DPI, block, entry, and replication as fixed effects. Dunnett tests were also performed to compare the accessions against the resistant and susceptible checks for each isolate using R [101]. The disease score distribution of each isolate for normality was tested using Shapiro–Wilk Test. Broad-sense heritability, H², was calculated for resistance to each isolate of P. capsici using the formula: H² = σ²_G/σ²_P, where σ²_G is the variance associated with the genotype and σ²_P is the variance associated with the phenotype.

Genotypic data

DNA extraction was performed using single leaf tissue from each chile pepper accession aged 30–45-days. The extraction was conducted using Qiagen DNEasy® plant extraction kits at the University of Minnesota Genomics Center (UMGC) DNA extraction facility (https://genomics.umn.edu/service/dna-extraction). DNA quantification was done using Picogreen (ThermoFisher Scientific, MA, USA). The samples were normalized to a concentration of 10 ng/µl for genotyping-by-sequencing (GBS). The UMGC performed GBS utilizing a single-enzyme digestion protocol as previously described [8, 58]. In this protocol, 10 units of restriction enzyme ApeKI (New England Biolabs®, Inc. MA, USA) were used to digest 100 ng of DNA per sample for 2 h at 75 °C, followed by 20 min of heat inactivation at 80 °C. The DNA samples were then heat killed after being ligated for 1 h at 22 °C with 200 units of T4 ligase (New England Biolabs®, Inc. MA, USA) and phased adaptors. Purification of the ligated samples was done followed by the amplification for the addition of the barcodes. Purification, quantification, and pooling of the GBS libraries was subsequently performed. DNA fragments within the size range of 300–744 bp size region were then specifically selected and diluted to a concentration of 1 nM. These selected fragments were subjected to sequencing with the Illumina NovaSeq™ 6000 (Illumina, CA, USA) using single end 1 × 100 reads.

Illumina ‘bcl2fastq’ software was used for the demultiplexing of the FASTQ files and first 12 bases from the beginning of each read were removed using Trimmomatic to get rid of the adapter sequences [102]. The alignment of FASTQ files to the ‘Zunla-1’ reference genome [103] was done using BWA (Burrows-Wheeler Aligner) [104]. To remove variants with genotype rates < 95%, minor allele frequency < 1%, and samples with genotype rates < 50%, raw VCF files were processed using VCFtools [105]. TASSEL [106, 107] was used to convert VCF files to HapMap format. Further, filtering of the HapMap file based on taxa was done and SNP markers with MAF < 0.05 and minor states were removed. Missing data imputation was done using the LD k-nearest neighbor genotype imputation (LD-knni) function [108] in TASSEL. ‘CMplot’ package in R was used to construct the SNP density plot per chromosome [109].

Population structure and linkage disequilibrium analysis

A total of 55,117 markers were used to identify the subpopulations within the diversity panel using an admixture model in STRUCTURE 2.3.4 [60]. The length of burn-in period and number of Monte Carlo Markov Chain (MCMC) replicates after burn-in were selected as 10,000 and the number of replications per K was set to 5. The number of clusters, K, were defined between 1 and 10 and the optimal K for the population was inferred using the STRUCTURE HARVESTER program (https://github.com/dentearl/structureHarvester/) [110, 111]. Convergence of the algorithm was verified by comparing the average values for alpha, α (i.e., relative admixture levels between populations for each replication per K). TASSEL 5.2.86 software was used to perform principal component analysis (PCA) to analyze the population structure where five major components were considered for further analysis. PCA biplot using the first two PCs and scree plot were constructed using R. TASSEL was also used to construct a neighbor-joining (NJ) phylogenetic tree which was further edited in iTOL (Interactive Tree of Life) website (https://itol.embl.de/) [112]. The PCA biplot and NJ tree were based on the optimum clusters obtained from STRUCTURE results. The “Sliding window” type function in TASSEL with a window size of 50 (0.05 Kb) was used to perform linkage disequilibrium (LD) analysis [106]. The “ggplot2” package in R was used to plot LD decay over distance using the TASSEL results [113].

Genomewide association study

A multi-locus genomewide mapping approach was used to detect significant loci linked with area under disease progress curve (AUDPC), disease severity index percentage (DSI%), percent resistant plants, and estimated marginal means for the three P. capsici isolates. The ‘lme4’ package in R was used to calculate the Best Linear Unbiased Estimates (BLUEs) [114]. Best Linear Unbiased Predictors (BLUPs) were computed considering the three isolates as different “environments” and performing combined analyses using the same package in R. Overall, three different BLUE phenotypic values (‘PWB-185’, ‘PWB-186’, ‘6347’), and a BLUP (‘PWB-185_PWB-186_6347’) for the traits were used for GWAS.

The multi-locus random-SNP effect mixed linear model (‘mrMLM.GUI’) package in R was used to perform GWAS [115]. Six multi-locus GWAS models, namely, Multi-locus random-SNP-effect mixed linear model (mrMLM) [116], Iterative modified-sure independence screening Expectation–Maximization-Bayesian least absolute shrinkage and selection operator (ISIS EM-BLASSO) [117], Polygenic-background-control based least angle regression plus empirical Bayes (pLARmEB) [118], Fast multi-locus random-SNP-effect efficient mixed model association (FASTmrEMMA) [119], Fast multi-locus random-SNP-effect mixed linear model (FASTmrMLM) [120], and Polygenic Kruskal–Wallis method with Empirical Bayes (pKWmEB) [121] were used with the default parameters. The GWAS model included five principal components and a kinship matrix (K-PC) to address population structure. A threshold of 3.0 LOD score was employed to declare markers as statistically significant in the GWAS models.

Candidate gene analysis

Candidate gene analysis was performed using EnsemblPlants [122], whereby genes located within 0.5 Mb of the SNP marker were considered as the candidate genes. An annotation file for ‘Criollo de Morellos-334’ (‘CM-334’; Genome assembly (GA)/ ASM512225v2) (C. annuum L.) was downloaded from the EnsemblPlants website (https://ftp.ensemblgenomes.ebi.ac.uk/pub/plants/release-55/gff3/capsicum_annuum/) [122, 123]. The list of candidate genes and their annotations for each significant SNP marker was subsequently derived from the annotation file using R. Candidate gene analyses were performed for significant markers that were common across all phenotypic datasets.