The phenotypic differences between ECW and ECW50R (bs5) were clear and easily distinguishable following inoculation at a relatively low bacterial concentration (10⁵ CFU/ml) ( Figure 1 ). The ECW leaf tissue developed necrotic lesions surrounded by yellow halos while the ECW50R tissue remained mostly green. In the GBS F₂ population, 91 out of 100 F₂s (19 resistant and 72 susceptible) were phenotyped with high confidence and thus were used for GBS step. The ratio of resistant to susceptible F₂s (1:3.8) was slightly lower than the expected ratio of 1:3 for recessive monogenic inheritance, however the difference was not statistically significant (X^{2 =} 0.824 at 1 degree of freedom; p=0.364).

The phenotype of ECW60R (bs6) resistance was not as distinct as bs5 ( Figure 1 ). As expected, bs6 resistance was characterized by extensive chlorosis. Out of 120 F₂s, 92 most clearly phenotyped individuals (29 resistant and 63 susceptible) were selected for GBS analysis. The ratio of resistant to susceptible F₂s (1:2.2) was not statistically different (X^{2 =} 2.087 at 1 degree of freedom; p=0.1486) from the expected 1:3 ratio.

A total of 169,398,995 reads were generated from the bs5 GBS library ( Supplementary Table 1 ). The GBS pipeline discovered 101 high quality SNPs that were polymorphic between the two parents, and those SNPs were selected for further analysis. The linkage analysis of 88 F₂s that could be genotyped identified thirteen linkage groups, and the bs5 resistance mapped to linkage group 1 in chromosome 3 with highest significance ( Figure 2 ; Supplementary Tables 2 , 3 ). SNPs between positions 134,620 and 1,098,542 of chromosome 3 were the most significantly associated with bs5 (p<0.0001). Genotyping of the F₂ population with CAPS markers spanning the linkage region confirmed 100% marker-trait co-segregation in the mapping population ( Supplementary Tables 4 , 5 ). The results indicate that bs5 is located towards the distal end of the short arm of chromosome 3, within a ~1 Mbp interval between 0.1 and 1.1 Mbp position.

A larger ECW × ECW50R F₂ population was developed to fine-map the position of bs5. Out of 1270 F₂s genotyped with flanking markers 3g_C0.134 and 3g_C1.11 ( Supplementary Table 4 ), 16 individuals were identified as recombinants and were phenotyped. Ten informative recombinants and F₃ RILs developed from six non-informative recombinants placed bs5 into an ~546 Kb interval between markers 3g_C0.134 (~0.4 cM) and 3g_C0.68 (~0.95 cM) with tight linkage with marker 3g_C0.26. ( Figure 3 ; Supplementary Tables 6 , 7 ).

An ECW bs5 super-scaffold was developed by concatenating C. annuum ECW scaffolds that align with in C. annuum UCD10X bs5 interval. This super-scaffold consisted of 535 Kbp sequence including gaps and flanking region and provided complete coverage of UCD10X bs5 interval ( Supplementary Table 8 ). Comparison of whole genome polymorphisms between bs5-fixed line (PI 163192 × ECW50R) and ECW identified a total of 1,718 variants in this region under stringent filtration (data not shown). However, only 28 variants were found to alter the protein sequences, which resulted in 14 putative candidate genes for bs5 resistance ( Table 1 ; Alignment S1-S14 ).

As the reference-based GBS pipeline only identified a small number of polymorphic markers, the reference-free UNEAK pipeline was used for mapping bs6. This pipeline discovered 133 SNPs from a total of 173,074,228 reads generated from sequencing ( Supplementary Table 9 ). Nine linkage groups were generated from the linkage analysis using genotyping information from 92 F₂ plants ( Supplementary Table 10 ), out of which the bs6 resistance phenotype was significantly (p < 0.0001) linked to SNPs on linkage group 3 ( Figure 4 ; Supplementary Tables 10 , 11 ). The linkage group was determined to be physically located in chromosome 6. CAPS markers were developed in the bs6-mapped region, and genotyping of the F₂ population validated the linkage between those markers and the resistance phenotype ( Supplementary Tables 12 , 13 ). The results indicated that bs6 was located within an ~21 Mbp interval between positions 168–189 Mbp in C. annuum UCD10X genome.

Five of the CAPS markers within the ~21 Mbp bs6 interval were initially used to more precisely determine the position of bs6. In a fine mapping F₂ population of 940 plants, 277 plants were identified as recombinants, 123 of which were homozygous for 60R alleles throughout part of the recombined region and were phenotyped as F₂ plants; genotyping of these F₂s delimited the resistance locus to an ~9.8 Mbp region between markers 6g_C171.79 and 6g_C181.60 ( Figure 5A ; Supplementary Table 14 ). F₃ RILs developed from 61 F₂s that recombined within the region were genotyped with eight new CAPS markers within the interval ( Supplementary Table 12 ); this delimited bs6 within an ~5.1 Mbp interval between markers 6g_C175.02 and 6g_C180.10 ( Figure 5B ; Supplementary Table 15 ). A second ECW60R × ECW F₂ population of 940 plants was developed and genotyped with new HRM markers ( Supplementary Table 12 ), and 41 recombinants between flanking markers 6g_H171.54 and 6g_H183.16 were identified and developed into F3 RILs. All 41 RILs were phenotyped and were genotyped with markers in the 5.1 Kbp interval, thereby delimiting bs6 to an ~656 Kbp region between markers 6g_H178.44 (~0.11 cM) and 6g_H179.10 (~0.11 cM) ( Figure 5C ; Supplementary Table 16 ).

The ECW bs6 super-scaffold spanned three C. annuum ECW scaffolds with a total size of 681 Kb, providing complete coverage of UCD10X bs6 interval ( Supplementary Table 8 ). A total of 1,718 variants were identified between ECW and ECW60R genome in this region after filtration. Annotation of those variants identified protein coding changes in eight genes, which are candidates for bs6 ( Table 2 ; Alignment S15-S22 ). Interestingly, four of those candidates are functionally annotated as ZED1-related serine/threonine kinases, and three have protein polymorphisms within the putative kinase domain ( Table 2 ).

For developing populations segregating for resistance, ECW50R and ECW60R were used as resistant parents for bs5 and bs6, respectively. ECW was used as susceptible parent for both populations. For both resistances, ECW was crossed with respective resistant parent to produce an F₁ population, which was self-pollinated to generate F₂ seeds. F₃ populations were generated by selfing of F₂s when necessary. F₂ recombinant individuals were self-pollinated, and progeny were genotyped to identify plants fixed for the recombined chromosomal segments (recombinant inbred lines (RILs)). A complete outline of all populations is presented in Supp. Image 1 . For all plants, seeds were sown in a seedling flat, and fourteen-day-old seedlings were transplanted to 10-cm pots containing Fafard Mix 4 (Fafard, Inc., Agawam, MA). For fine-mapping F₂ populations, the plants were grown in 242-well trays (Speedling Inc., Sun City, FL) containing Speedling peat-lite soilless media (Speedling Inc., Sun City, FL). The transplants were grown in a greenhouse at temperatures ranging between 20-30 °C.

As the resistant responses due to bs5 and bs6 do not result in HR induction, they are differentiated from the susceptible response by infiltration of bacterial suspension into pepper leaves at a low concentration (Stall, 1981). In contrast to the development of necrotic lesion in susceptible pepper, the bs5 resistance only causes a slight yellowing of the infiltrated area and the bs6 resistant response is characterized by a more intense chlorosis (Vallejos et al., 2010). Xe race P6 strain Xv157 was grown in nutrient broth (BBL, Cockeysville, MD) overnight at 28 °C with constant shaking. Bacterial cells were pelleted by centrifugation, the supernatant was discarded, and the cells were re-suspended in sterile tap water. The bacterial suspension was adjusted using Spectronic 20 Genesys spectrophotometer (Spectronic Instruments, Rochester, NY) to OD₆₀₀ = 0.3, which is approximately 10⁸ CFU/ml, then diluted to 10⁵ CFU/ml in sterile tap water. The resulting bacterial suspension was infiltrated with a syringe and hypodermic needle into the mesophyll of the first and second true leaf of five- to six-week-old pepper plants. Inoculated plants were maintained in a greenhouse for disease development, and the plants were evaluated three weeks after inoculation. Plants showing confluent necrosis were rated as susceptible, else they were rated as resistant for the respective resistance. For bs6 resistance, the disease screen of each RIL was repeated multiple times to obtain accurate phenotypic result.

Foliar tissue from young leaves was lyophilized and used for DNA extraction. Genomic DNA was extracted using the Qiagen Plant DNeasy Mini Kit (Qiagen, Germantown, MD) according to the manufacturer’s instructions. The DNA was normalized to 5 ng/µL based on quantification with a Synergy 2 multimode microplate reader (Biotek Instruments, Winooski, VT) with the Quant-iT PicoGreen double-stranded DNA quantification assay (Thermo Fisher Scientific, Waltham, MA). A 96-plex (ninety one F₂s, a single F₁, and two each of ECW and respective resistant parent) ApeKI GBS library was constructed using a previously published protocol (Elshire et al., 2011). Barcode-adapter titration indicated that 0.9 ng µL^-1 of each barcode-adapter per 50 ng of genomic DNA produced satisfactory libraries without dimer formation. The barcode-adapter titration mixture and the final GBS library were analyzed on an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA) to ensure acceptable fragment size distribution and quantities. The GBS library was diluted to 3.6 pM and sequenced on one lane (single end, 101 base pair read length) of an Illumina HiSeq 2500 (Illumina Inc, San Diego, CA) at the Genomics Resources Core Facility (Weill Cornell Medicine, NY).

The raw sequencing reads were processed in TASSEL version 3.0 (Bradbury et al., 2007) using either the reference genome-reliant TASSEL-GBS pipeline (Glaubitz et al., 2014) or the reference-free UNEAK pipeline (for bs6) (Lu et al., 2013). For both pipelines, high quality sequencing reads that contained a barcode-adapter, an ApeKI restriction site, and an inserted genomic sequence (hereafter termed GBS tags) were identified and selected based on polymorphism between parents. In TASSEL-GBS pipeline, the reads were aligned with the bwa v0.7.8 (Li and Durbin, 2009) to the C. annuum UCD10X reference genome, release 1.1 (Hulse-Kemp et al., 2018) to identify polymorphisms ( Supplementary Table 1 ). For the UNEAK pipeline, reference genome information was not necessary, and SNPs were identified by pairwise alignment of all unique sequence tags across the entire dataset ( Supplementary Table 2 ). Raw read files from sequencing of GBS libraries are deposited in NCBI SRA under bioproject PRJNA863731.

Polymorphic SNPs identified between the parental lines were employed for linkage analyses using MapDisto v1.7 (implemented within Microsoft Excel 2007), (Lorieux, 2012). The parameters in linkage analyses were a minimum LOD=5, a maximum r=0.3, and the ‘Kosambi’ mapping function. The loci were ordered within each linkage map using the auto-order function. QTL analysis was conducted for each population to determine the association between the SNPs within a linkage group and resistance to race P6. Single marker analysis was performed using the R/qtl package in R v3.3.1 (Broman et al., 2003).

Cleaved Amplified Polymorphic Sequence (CAPS) markers were designed for validating the mapping results from GBS and for fine mapping. Primers for the markers were designed using Primer 3 software (Untergasser et al., 2007) utilizing SNPs identified from GBS. DNA was extracted using a Cetyltrimethylammonium Bromide (CTAB) method (Doyle and Doyle, 1987) and polymerase chain reaction (PCR) was carried out with Phire Hot Start II DNA polymerase (Thermo Fisher Scientific, Waltham, MA) in a 10 μl volume, which consisted of 2 μl of DNA (adjusted to ~20 ng/μl), 4.89 μl of HPLC-H₂O, 2 μl of 5X Phire Reaction Buffer, 1 μl of dNTPs, 0.03 μl each of forward and reverse primers, and 0.05 μl of polymerase. The amplicons were digested with appropriate restriction enzymes according to the manufacturer’s recommendations (New England Biolabs, Ipswich, MA). Results were detected using electrophoresis on 3% agarose gels stained with ethidium bromide.

High Resolution Melting curve (HRM) markers were developed from SNPs identified from GBS. Primers were developed using the IDT PrimerQuest (idtdna.com/Primerquest). DNA was extracted using a NaOH rapid DNA extraction method (Lee et al., 2017). The 5 μl PCR reactions were mixed with 2x AccuStart II PCR SuperMix (Quantabio, Beverly, MA), 0.5 μM of each primer, and 20x EvaGreen Dye (Biotium, Hayward, CA) and run as follows: (95 °C @ 60s) + 40 × ((94 °C @ 5s) + (Tm @ 10s) + (72 °C @ 15s)) + (72 °C for 60s), where Tm is the annealing temperature. For allele determination, melting curve analysis was performed by scanning the PCR product in a LightCycler 480 Instrument II (Roche, Pleasanton, CA).

A modified microprep protocol was used for DNA extraction for whole genome sequencing of ECW60R (Fulton et al., 1995; Sharma et al., 2022b). DNA concentration and purity was verified using NanoDrop (Thermo Fisher Scientific, Waltham, MA). Subsequently, DNA was cleaned using DNeasy PowerClean Pro Cleanup Kit (Qiagen, Germantown, MD) following the manufacturer’s recommendations. Illumina sequencing library was prepared using a Nextera DNA Flex Library Prep Kit (Illumina Inc, San Diego, CA) using the protocol recommended by the manufacturer. The DNA was sequenced to produce 100 base-pairs (bp) paired end reads in one lane of Illumina HiSeq 3000 at University of Florida Interdisciplinary Center for Biotechnology Research.

The C. annuum ECW whole genome sequence (GCA_011745845.1) was only assembled to scaffold level at the time of analysis (Kim et al., 2017). To produce contiguous sequence, the bs5 or bs6 fine mapped intervals were blasted against the reference genome C. annuum UCD10X (GCF_002878395.1). All ECW scaffolds with query coverage greater than 2% and matching to unique regions were identified and concatenated together in correct order and orientation to produce ECW super-scaffolds for bs5 and bs6. The super-scaffolds also consisted of 5 Kbp region up- and down–stream from flanking markers and 3 Kbp gap between stitched scaffolds. The super-scaffolds were aligned with C. annuum UCD10X resistance intervals to verify complete coverage.

The ECW genes were predicted de-novo to overcome differences in gene annotations between reference genomes. ECW gene prediction model was developed using BRAKER v2.1.6 (Brůna et al., 2021). Within BRAKER, three publicly available ECW RNAseq sequences (SRR13488414, SRR13488423, and SRR13488424) were aligned to C. annuum ECW genome sequence (GCA_011745845.1) and supplied to Genemark-ET v4.68 (Lomsadze et al., 2014) to generate hints for training Augustus v3.4.0 (Stanke et al., 2008). The resulting ECW gene prediction model was used to identify potential protein coding regions in the bs5 and bs6 super-scaffolds. The genes were validated based on their posterior probability and annotation of homologous regions in C. annuum UCD10X or C. annuum CM334 annotation.

Polymorphisms for bs5 were identified using whole genome bulk sequences of PI 163192 × ECW50R F₂ population, which is fixed for bs5 gene (Sharma et al., 2022b). For bs6, the whole genome sequence of ECW60R was used. The sequences were analyzed using an in-house pipeline. The quality of the reads was verified with FastQC 0.11.7 (bioinformatics.babraham.ac.uk/​projects/​fastqc) and the adapters were trimmed using trim_galore v0.6.5 (Krueger et al., 2021). The trimmed reads were aligned to C. annuum ECW genome using Bwa-mem2 v2.2.1 (Vasimuddin et al., 2019). The resulting alignment file was used for variant calling with the HaplotypeCaller tool in GATK 4 (DePristo et al., 2011). The variants were filtered under high stringency as follows: depth ≥ 12, quality-normalized depth ≥ 10, mapping quality ≥ 50, and reference allele depth ≤ 0.1 × alternate allele depth. The sequencing data for PI 163192 × ECW50R F₂s has previously been deposited in NCBI/​ENA/​DDBJ database under bioproject PRJNA789991. ECW60R whole genome sequence is deposited under bioproject PRJNA863893.

The coordinates and allelic sequence of high-quality polymorphisms in bs5/​bs6 super-scaffolds were derived from variant calling of C. annuum ECW scaffolds with an in-house script. The polymorphism were annotated with snpEff v5.0 (Cingolani et al., 2012) using a custom super-scaffold variant annotation database built using previously described sequences and protein coding regions. Only the variations that result in protein coding changes were selected to identify potential candidate genes. Potential homologs of candidate genes in other C. annuum genomes were identified by blasting the predicted amino acid sequences of those genes, which also provided the functional annotations of the candidates. Finally, protein domains containing the polymorphisms between ECW and ECW50R/​ECW60R were identified by Pfam (Mistry et al., 2021) and InterPro search (Blum et al., 2021).