Introduction {#Sec1}
============

Chili pepper is an economically important horticultural crop in Solanaceae family that also includes potato, tomato, eggplant, petunia and tobacco. The Solanaceae family includes more than 3,000 varied species with the similar numbers of chromosomes (n = 12) but significantly different genomic sizes. Peppers have been used as a vegetable, condiment, spice, medicine, coloring agent and source of vitamins^[@CR1]--[@CR3]^. The most common cultivated pepper species are *Capsicum annuum*, *Capsicum frutescens*, *Capsicum chinense*, *Capsicum pubescens*, and *Capsicum baccatum*^[@CR4],[@CR5]^. Though pepper consists of several potential economic values, fungi, bacteria and viruses cause heavy losses in pepper fruit production. Powdery mildew (PM) is the most common devastating fungal disease in pepper and is caused by *Leveillulataurici*. In an agricultural setting, this disease could be controlled using agrochemicals or genetic resistance lines. The selection of good PM resistance varieties through traditional breeding potentially requires more than 10 years. Hence, molecular marker-assisted breeding is the current plant breeding method of choice, and the most frequently used markers include single-nucleotide polymorphisms (SNPs). DNA-based molecular markers are employed in plant breeding for genetic diversity and genome association analyses^[@CR6]--[@CR9]^. Over the last decade, major innovations in sequencing technologies and bioinformatics have been achieved, prompting a transition from classical conservation genetics to conservation genomics^[@CR10]--[@CR13]^. Rapid innovations in genome sequencing platforms, such as next generation sequencing (NGS), provide numerous opportunities for transcriptome assembly, functional annotation of genes, and identification of molecular markers^[@CR14],[@CR15]^. New software tools in NGS technology enable the cost effective identification, confirmation, and evaluation of genetic markers on a large scale.

SNPs have been accepted as potential selection markers in genome-wide studies given the high density of markers near loci of interest^[@CR6]^. NGS technologies have identified genome-wide SNPs in several crops, such as bean^[@CR16]^, barley^[@CR17]^, cassava^[@CR18]^, cabbage^[@CR19]^, grape^[@CR20]^ and maize^[@CR21]^. In pepper, several thousand genetic markers, especially SNPs have been discovered^[@CR22]--[@CR28]^. Recently, Kim *et al*.^[@CR29]^ sequenced and assembled the pepper genome (*Capsicum annuum cv*. CM334) at a genomic size of 3.48 Gb. This reference genome will provide the opportunity to improve quality, cultivation, and disease resistance in *Capsicum* species. The aim of this research is to discover SNP variants for future marker-assisted breeding studies related to PM resistance using *Capsicum annuum cv*. CM334 as a reference for data mining. Thus, in the present study resequencing of two pepper varieties, *Capsicum baccatum* (PRH1- PM resistant line) and *Capsicum annuum* (Saengryeg - PM susceptible line), using the HiSeq. 4000 Illumina platform and the genome wide identification of SNPs have been implemented

Results {#Sec2}
=======

Genome sequencing, pre-processing and alignment of reads to the reference genome {#Sec3}
--------------------------------------------------------------------------------

A summary of the sequencing, sequence preprocessing, and alignment to the read mapping were presented in Table [1](#Tab1){ref-type="table"}. In total, 130,370,103 and 118,588,231 paired-raw reads were discovered for PRH1 and Saengryeg, respectively, with an average length of 151 bp. A total of 19.69 and 17.91 Gb paired-end raw reads were recorded for both pepper varieties. The total genome coverages were ≒ 11.31× and ≒ 10.29× of the reference genome. The Solexa QA (v.1.13) package was used to generate high-quality clean reads. Raw reads were assessed for quality, and impractical parts were discarded. After the removal of adaptor sequences, ambiguous and low-quality reads (Q value \<20), a total of 97,216,537 and 88,964,871 reads were discovered for PM resistant and susceptible pepper varieties, with ≒ 5.61× and ≒ 5.17× of genome coverage respectively. After the removal of non-specific reads, the remaining reads were mapped to the reference genome. A total of 194,523,074 and 177,929,742 clean, high-quality reads were recorded for PRH1 and Saengryeg, respectively, compared with the reference genome, covering 88,448,386 (45.47%) and 1,080,500,795 (39.24%) of mapped reads, respectively.Table 1Summary of sequencing, sequence pre-processing and alignment of reads to the reference genome.SampleRead parametersPRH1SaengryegRaw read dataNo. of reads130,370,103118,588,231130,370,103118,588,231Avg. length (bp)151151151151Total length (Gb)19.6917.9119.6917.91Genome coverage\#≒11.31X≒10.29XCleaned dataNo. of reads97,261,53788,964,87197,261,53788,964,871Avg. length (bp)1201218181Total length (Gb)11.6910.797.837.21Trimmed/raw\*59.4160.2439.7840.25Genome coverage\#≒5.61×≒5.71XRead mappingNo. of total reads194,523,074177,929,742No. of mapped reads (%)88,448,386 (45.47)1,080,500,795 (39.24)Mapped region\*\* (%)1,080,500,765 (39.24)2,514,912,154 (91.34)^\*^Trimmed/raw: total length of trimmed read / total length of raw read.^\#^Genome coverage: Total length of all reads divided by reference genome size (3.48 Gb).^\*\*^Mapped region: Coverage of read mapping relative to the reference genome.

Identification and distribution of SNP markers {#Sec4}
----------------------------------------------

Genome-wide SNPs were identified using an improved BWA-SAMtools workflow. The high-quality filtered reads of PRH1 and Saengryeg were mapped to the reference genome. A total of 6,213,009 and 6,840,889 SNPs were identified for both pepper varieties. Based on the SNP ratio to the read map, SNPs were classified into homozygous, heterozygous and other types. Among the identified SNPs, 88.59% homozygous, 3.65% heterozygous, and 7.76% other types of SNPs were determined in PRH1. Likewise, in Saengryeg, 95.04% homozygous, 1.91% heterozygous, and 3.05% other type SNPs were identified. The occurrence of low percentage of heterozygous SNPs in both lines was due to the relatively low sequence depth and rigid SNP calling requirement. Capsicum consists of 12 chromosomes, and the SNPs are distributed evenly across all chromosomes. Our further analysis revealed that the number of SNPs differed in chromosome 1 to 12 for the two pepper varieties (Fig. [1](#Fig1){ref-type="fig"}). The greatest number of homozygous SNPs were noted in chromosome 10 (1,096,754) in Saengryeg whereas chromosome 1 consisted of maximum number of homozygous SNPs (601,032) in PRH1. Similarly, chromosome 1 in PRH1 possessed higher number of heterozygous SNPs (23,932) and chromosome 12 consisted of maximum heretozygous SNPs (15,942) in Saengryeg. However, the least number of SNPs was discovered on chromosome 8 in both pepper varieties. The detailed dataset for the chromosomal distribution of SNPs is listed in Table [2](#Tab2){ref-type="table"}.Figure 1SNP distribution observed per 1 Mb chromosome. Thedistribution of SNPs detected with resequencing of pepper varieties along 12 chromosomes. The horizontalx-axis denotes the length (Mb) of chromosome and y-axis represents number of SNPs254x190 mm.Table 2Distribution of SNPs in the chromosomes of PRH1 and Saengryeg.Chromosome No.PRH1SaengryegHomozygousHeterozygousHomozygousHeterozygous1601,03223,932692,3269,9772400,51315,849388,7708,4283557,18518,211466,43911,4414420,71320,654296,1839,3455405,74421,116577,30412,7616448,88620,823516,71010,5637552,12017,025620,30614,0058304,39511,915153,5187,2299384,00916,739890,1359,78710460,91223,5171,096,75410,23611499,17614,391467,19310,27812469,52322,456301,55615,942Total5,504,208226,6286,467,194129,992

Annotation of SNPs based on their position in the pepper genome {#Sec5}
---------------------------------------------------------------

The SNPs were classified into two main categories (intergenic or genic region) according to their position in the pepper genome sequence. Further genic SNPs were sub-classified as intron and coding DNA sequences (CDS). A total of 6,213,009 and 6,804,889 genome-wide SNPs were discovered for PRH1 and Saengryeg, respectively. Of the discovered SNPs, 5,781,951 (93.06%) and 6,695,385 (93.39%) of intergenic SNPs were recognized for PRH1 and Saengryeg, respectively. Further, these SNPs were classified into homozygous, heterozygous and other type depending upon the ratio to read map. In addition, 82.28% and 93.58% of homozygous type SNPs were identified in the intergenic region for PRH1 and Saengryeg, respectively. We discovered that the number of SNPs in intron was greater than that of CDSs in the genic regions. Most of the SNPs were located in the intergenic regions and were classified as homozygous type (Table [3](#Tab3){ref-type="table"}). All the identified SNPs were analyzed for polymorphisms between PRH1 and Saengryeg. A total number of 15,941,182 SNP loci were identified with respect to the reference genome. Of the identified SNP loci, 4,887,031 polymorphic and 469,978 non-polymorphic loci were identified between PRH1 and Saengryeg. The genomic distribution of polymorphic SNP markers is presented in Fig. [2](#Fig2){ref-type="fig"}. High-resolution melting (HRM) marker primers were identified by targeting SNPs to discriminate between two lines. Among the polymorphic SNPs, 4,164,456 HRM candidates were identified, and 597,434 primer sets were selected. A total of 306,871 HRM primer markers were recommended for further breeding purposes (Supplementary file [S1](#MOESM1){ref-type="media"}). These sets of HRM primers possibly discriminate between the two lines.Table 3Summary of SNP classification by genome structure.SampleTotal no. of SNPRegionTotalHomozygousHeterozygousOtherPRH16,213,009Introns280,076258,9047,82413,348CDS150,932133,4367,02510,471Genic region431,058392,38314,85123,824Intergenic region5,781,9515,111,825211,777458,349Saengryeg6,804,889Introns69,54264,9941,9382,610CDS39,95534,2122,4493,294Genic region109,50499,2104,3885,906Intergenic region6,695,3856,367,984125,604201,797Figure 2Genomic distribution of polymorphic SNP markers (PRH1 Vs Saengryeg) 254 × 190mm.

Chromosome-wise characterization of polymorphic SNPs {#Sec6}
----------------------------------------------------

In order to gain deeper insight into the SNPs associated with the genes involved in disease resistance and stress tolerance process, chromosome-wide functional annotation of polymorphic SNPs were performed. The distribution of SNP markers were analyzed in each chromosomes and the functional characterization of genes with higher polymorphic SNPs have been carried out. Overall, the majority of the genes with high polymorphic SNPs widely involved in carbohydrate metabolism, transcription regulation, ion binding, nucleotide binding, protein transport, fatty acid metabolism, receptors, photosynthesis, post-translational modifications, stress response, regulatory elements, proteolysis, secondary metabolism, biosynthesis, diseases resistance, and others. However, in each chromosome the genes with various functions displayed the major proportion (Fig. [3](#Fig3){ref-type="fig"}). For instance, in chromosome 1 the SNPs were highly identified in genes involved in carbohydrate metabolism followed by transport related genes. Transcription regulation related genes consisted of numerous polymorphic SNPs in chromosome 2 and 8. In chromosome 3, the genes associated with post-translational modifications consisted of more polymorphic SNPs. Likewise the diseases resistances genes with high polymorphic SNPs dominated the chromosome 4. Moreover, nucleotide/ion binding and ion transport genes with polymorphic SNPs were identified in chromosomes 5, 6, 7, 9, and 10. Genes involved in biosynthesis consisted of vast number of SNPs in chromosome 11 and 12.Figure 3Chromosome wise annotation of polymorphic genic SNPs associated with important functionsinPRH1 and Saengryeg.

Identification of polymorphic SNP markers associated with pathogen resistance genes {#Sec7}
-----------------------------------------------------------------------------------

In total, 6281 SNPs associated with 46 pathogen resistance genes with nucleotide binding site-leucine rich repeat (NBS-LRR) motif were identified in the introns and coding regions of the genes (Supplementary file [S2](#MOESM1){ref-type="media"}). The occurrence of SNPs related to NBS-LRR genes in each chromosome has been listed in Fig. [4](#Fig4){ref-type="fig"}. The maximum number of SNPs was distributed in chromosome 4, whereas the least number of SNPs was observed in chromosome 8. Moreover, the PM resistant line PRH1 consisted of greater number of NB-LRR linked SNPs in comparison with the susceptible line Saengryeg. Overall, the occurrence of higher number SNPs particularly associated with the NB-LRR resistance genes could play a vital role in the attribution of PM resistance.Figure 4Chromosome wise occurrence of SNPs associated with NBS-LRR genesinPRH1 and Saengryegincomparison with reference genome.

Phenotypic evaluation for PM resistance and validation of SNP markers {#Sec8}
---------------------------------------------------------------------

In order to assess the disease resistance indexes, the parental types and the F~4~ population were co-cultivated with the powdery mildew pathogen. The infection range observed in the plants has been categorized from 1--5 scale from PM resistance to susceptible (Supplementary Table file [3](#MOESM1){ref-type="media"}). The parents of F~4~ population exhibited contrasting degree of resistance to the PM disease. The *C. baccatum* variety (PRH1) displayed high resistance scale of 1, whereas the *C. annuum* variety (Saengryeg) exhibited resistance score of 5. However, among the 45 individuals in F~4~ population, 11 exhibited the resistance score of 1 followed by 22 plants resulted in the moderate disease resistance level of 3 and 12 plants displayed the severity with the index of 5. Further, to validate the identified SNP markers, HRM assay in both the parental types along with the F~4~ population of 45 progenies has been performed. Among the 36 HRM primers employed, 19 primers significantly distinguished the resistant and susceptible progenies in the F~4~ population. The HRM primers employed in this study have been listed in the Table [4](#Tab4){ref-type="table"}.The representative HRM melt curves obtained for the parents with the heterozygous SNP variation of G/A and C/A have been illustrated in Fig. [5](#Fig5){ref-type="fig"}. Moreover the majority of heterozygous SNPs were observed to be prominent among the population studied. Thus, the current HRM platform provided a suitable approach for the validation of SNP markers among the population.Table 4List of HRM primers designed for genotyping polymorphic genic SNPs from each chromosome.**Primer no**.**Locus nameGene name5′-3′primer sequence3′-5′primer sequenceResult of HRM analysis validated in F** ~**4**~ **population**1CA01g00370Serine/threonine protein kinase%2C putativeCGGCCAATGTATCAAGACTCGAACGAATTCAACAACCGCGTPositive2CA01g02310Xpa-binding protein%2C putativeTCCCTTCTGCGGTTTTCCTCTGTTGCAAACTTCTCCTTGTAGGPositive3CA01g04020Kinesin heavy chain%2C putativeCCCACTGGTGAAAGCAGTGTTGGAGAGAAGGCCTCAATGGPositive4CA02g00020DNA-repair protein UVH3%2C putativeTGGTCAGGTAATGGTGGTTCTCTCTCCCTCATCTGGCAAACANegative5CA02g00720Pentatricopeptide repeat-containing protein%2C putativeAGAGCACTAACCTCTTTAGCAGACTGCAAAGACCCCACAGAPositive6CA02g02750MYBR domain class transcription factorACAGTCATACTAGATGAAGGCGGTGATGCAATGTGGTCAGATGANegative7CA03g00110Beta-galactosidaseAGTAACTGATGGAATTTCGGAATGGATGCGTTTTAGCCTGACTPositive8CA03g00740Small subunit processome component-like proteinTCCCAGCATACTCGTCCAACCCTCAACCTAGGCATGCCAANegative9CA03g15330PREDICTED: Golgi to ER traffic protein 4 homologTGGTTAGTCTTTCCTAATCCGGTCTATTTCTTTTTCCATTCCATTGCPositive10CA04g00830Phosphatidylinositol 4-kinase%2C putativeGGGGGCTAGTCTTCTCTTCTGGCAACAAGGTGGAAAGACGNegative11CA04g00250PREDICTED: transmembrane emp24 domain-containing protein p24beta3-likeCGGATCATCCCGGCATTGATTCACCTCCGATTCACAACTCANegative12CA04g00360Protein transport protein sec. 23%2C putativeGCACGCCCATACCTTGTCAAATCAATGCCAAGCCCATCCAPositive13CA05g00010RNA polymerase II transcription mediators isoform 1CAACGAGGCTGACCGAAAGACTCCACTCGCCCATCTTCTCPositive14CA05g00320Folylpolyglutamate synthaseGGTGGGGGCTTTTGTCTTCTACTACATCTTCTGAGGTAACACCNegative15CA05g15050PREDICTED: mediator of RNA polymerase II transcription subunit 33ACCACCGTTTCAATCCCTTGCACGTGTCAGGATTCATAAGCTPositive16CA06g00010Kinesin heavy chain%2C putativeTGAAGCCGCCTCGAATTTCTAATGAGACTTCGAGGGGCACNegative17CA06g01280Myosin XI%2C putativeATAGACCCCGGCTCAGGAATGCAAAGGTAGCTCCACCACTPositive18CA06g01570PREDICTED: TBC1 domain family memberGGCAGGAAGATACAATAAATGTACAGCAGTATCGTGATTTCATTTGGTNegative19CA07g03700PREDICTED: synaptotagmin-5-likeAGTAAGGTCAAATGTGGAGCCAAGAACGTTAATACTGGCCATCGNegative20CA07g04200Transducin family proteinTGCGAACTTAAGGAAAAAGAAGCAGTAATGCTTGTCGGGAGCCTPositive21CA07g12460ForminGGGATAACGCTCTTCCATATGGACATGTCTGACAGAGGGTGCANegative22CA08g00950Transcription cofactor%2C putativeACACTGAGATGCATGCACCATACCTGGTTTTGGCTGTGTTPositive23CA08g08740DNA-directed RNA polymeraseACAACAGGGACATGATTTCATCAACACTAAACCCTTCTGTGCACAPositive24CA08g09730PREDICTED: protein ZINC INDUCED FACILITATOR-LIKE 1-like isoform X3TGTGTGTCGAAGCAATTGATCTGTTGGAAGATTTGTCAATATCAPositive25CA09g00140O-linked n-acetylglucosamine transferaseCTGCACATAGAATTCTTGCCCATGGGATTGTTTCGTGCTTTTNegative26CA09g01180Vacuolar protein sorting-associated proteinTTGTCCTCCTCCTCAGATGAACCACCAGCAAGAACGTCAANegative27CA09g14940Beta-amyrin synthaseTGGCACCATTTTTAAACAACAACAGTCAGAAGCACACTGTGAPositive28CA10g01250PREDICTED: heterogeneous nuclear ribonucleoproteinTGATGAGCTCGGAGGAGTCAAAGTGGCTGGGATTCAAGGGNegative29CA10g01280Protein binding protein%2C putativeGGGTGAGTTTCCTAAGAGGTCCCAAATCACATGGCCAAACGCPositive30CA10g07870Amidase%2C putativeGCTGCAGCAATGTAATTGGACCTCTGACCATCATCGCTGANegative31CA11g11870Xanthine dehydrogenaseACCTTGACTGGTACACTTTTTCAAGTGATGACGGACAATTGTGTPositive32CA11g15420Tubulin family proteinGGCCTCATAACACCGTGGAATTACCAGCAGCATTGATCGANegative33CA11g15430ABA aldehyde oxidaseTTAATGGAGGCTTCAGAGAGAGCTTGGGACTCTTGAAAGAAGCPositive34CA12g01070IsoleucyltRNA synthetase%2C putativeACAACACCCATCGACTTCCCTGCAGAGCCAGATTTCAGGTPositive35CA12g02370N-like proteinTGGTGTTTTTCCATTTGCCTTCTCTAGAACGTAAGGGTATTCANegative36CA12g22510PREDICTED: pleiotropic drug resistance proteinACCGAGTCGAAAGAGGAAGCAAGGGCAGAGTCGAGCTTTCNegativeFigure 5HRM melt curve and temperature peaks obtained from candidate SNPs between *C. baccatum* (AR1) and *C. annum* (TF68) illustrating the G/A and C/A SNP variation.

Discussion {#Sec9}
==========

In general, a primary requisite of genotyping of all the individuals in a population is necessary for trait mapping in traditional approaches of breeding, which is a highly expensive, labor intensive and time consuming process. Moreover, the occurrence of mere levels of variations or polymorphism also acts as a vital challenge during molecular marker discovery. In order to address these difficulties, next generation sequencing (NGS) strategies have been widely applied in genomics based on breeding of important agricultural and horticultural crops. Recent advancements in NGS technology have facilitated the routine use of high-throughput, low-cost markers for plant breeding programs. New software tools enable the discovery, validation, and assessment of genetic markers on a large scale. Among different marker systems, SNPs are the most important and attractive DNA-based molecular markers used for genetic diversity and genome association analyses and comparative genetics in plant breeding^[@CR6]--[@CR9]^. SNP markers are highly polymorphic, co-dominant, precise, reproducible, high-throughput, economical and informative^[@CR28]^. Moreover, the discovery of genome-wide SNPs aids in the improvement of marker assisted selection, particularly for the identification of traits associated with disease resistance. In this study, a complete genome resequencing of two pepper varieties with contrasting powdery mildew (PM) tolerance ability, PRH1 (PM resistance) and Saengryeg (PM susceptible), has been examined for the identification of SNP markers associated with powdery mildew resistance. The available whole genome sequence information of *Capsicum annuum* cv CM334 has been utilized as the reference genome to enable the comparison between the *C. annum* and *C. baccatum* lines used in this study.

In the current endeavor,intersepecific breeding of sexually incompatible pepper species has been performed due to their potential traits. For instance, the *C. baccatum* is well-known for fruit quality, disease resistance, and high contents of valuable secondary metabolites^[@CR29]^. Therefore, the interspecific breeding of peppers results in progenies with high fruit quality and disease resistance. The *C. baccatum* variety used in this study displayed resistance to powdery mildew and anthracnose diseases. Hence, the whole genome re-sequencing (WGRS) based on discovery of SNPs in the variable pepper varieties could enhance the understanding of SNPs associated with disease resistance. The resequencing and SNP discovery resulted in the identification of 6,213,009 SNPs for PRH1 and 6,840,889 SNPs for Saengryeg. The SNPs identified in the present study were higher than the SNPs discovered by Nimmakayala *et al*. in *C.annuum* and *C.baccatum* varieties using genotyping by sequencing approach^[@CR30]^. The report suggested the collective identification of 36,621 potential SNP markers linked to various genomic regions in in *C. annuum* and *C. baccatum* that can be utilized for the genome wide association studies in pepper varieties^[@CR30]^. Moreover, the identified SNPs in the present study have been majorly categorized into homozygous type with 88.59% and 95.04% for PRH1 and Saengryeg, respectively. This suggests that the sequence of reference genome could be generated from homozygous loci. Further, the chromosomal distribution of SNPs in the pepper genome revealed that a total of 10.92% of homozygous SNPs were located on chromosome 1, and 16.96% of the homozygous SNPs were located on chromosome 10 for PRH1 and Saengryeg, respectively.

In addition, the distribution of SNPs within the pepper genome illustrated the occurrence of higher percentage of SNPs in intergenic regions compared with genic regions. Likewise, several SNPs were identified in the intronic region than in CDSs. Similar results were also reported in tomato by Kim *et al*.^[@CR31]^. Furthermore, the location of SNPs plays a vital role, particularly SNPs should be located in intragenic regions to implicate the phenotypic traits. These SNPs are expected to be applied to marker assisted selection because they could be considered as functional markers. A total of 5,941,182 SNP loci have been detected between Saengryeg and PRH1. Of them, 30.63% SNPs were distributed in polymorphic loci. Potential polymorphic homozygous SNPs were filtered to discover breed-specific markers in both of the pepper varieties. HRM analysis has been applied to identify precise, cost-effective and efficient tool to detect sequence variations, such as SNPs^[@CR32]^. This technique has been successfully implemented to identify SNPs that have been used for genotype discovery, genetic mapping and mutation scanning^[@CR33]--[@CR36]^. Among the discovered homozygous type polymorphic SNPs, 597,434 HRM marker primers were identified that potentially discriminate between two lines. Of them, 306,871 HRM primers were recommended for further experimental research related to PM-based melting patterns and amplification efficacy.

The numerous amount of polymorphic SNPs identified in the genic region were functionally annotated in each chromosome to gain deeper insight into the SNPs associated with the genes involved in disease resistance. A comparative genetics study on the resistance genes in Solanaceae family has shed light on to the potential loci in different chromosomes linked with disease resistance^[@CR37]^. The vital R genes associated with disease resistance were conserved among the related species such as pepper, tomato, and potato^[@CR37]^. The current results revealed that each chromosome consisted of several SNPs associated with the genes involved in vital metabolic processes. However, chromosome 4 consisted of larger set of SNPs associated with disease resistance in comparison with other chromosomes. According to Grube *et al*.^[@CR37]^, the diseases resistance gene loci located in the chromosome 4 of pepper could render resistance against fungal pathogens. Correspondingly, chromosome 4 could play a vital role in encompassing the genes required for disease resistance in pepper. Moreover, the chromosome 5--10 consisted of SNPs related to genes involved in ion and metal binding. The roles of ion/metal binding genes are inevitable particularly under stressed conditions in pepper plants. The uptake and transportation of nutrients and water from the environment to the plant is a complex as well an important process for the improvement of physiological functioning of plants in stress. Hence, the SNPs related to these genes could act as a vital marker under stress.

Furthermore, higher number of polymorphic SNPs associated with disease resistance genes such as NBS-LRR were also identified in chromosome 4. Among the two varieties, the resistant PRH1 possessed higher distribution of polymorphic SNPs related to NBS-LRR genes. In plants, NBS-LRR is a large family of proteins encoded by the resistance genes and NBS-LRR proteins involved in the recognition of pathogens^[@CR38]^. Several reports suggested the importance of NBS-LRR proteins in the resistance against numerous diseases including powdery mildew in plants^[@CR39]--[@CR41]^. In the present study, polymorphic SNPs were identified in genes encoding for LRR receptor-like serine/threonine-protein kinase, F-box/LRR, TIR-NBS-LRR resistances protein, CC-NBS-LRR resistance protein, and TIR1 like protein, etc. Hence, the identification of SNPs associated with the disease resistance genes could aid in the enhancement of screening processes in the molecular breeding of pepper with powdery mildew resistance.

The identified SNPs were validated using HRM primers in the parents and F~4~ population derived from the *C. annuum* and *C. baccatum* varieties. The HRM primers were selected from all the chromosomes and evaluated in the parents and the population. Among the tested primers, 19 primers were able to distinguish the population and the results were correlated with the phenotypic disease evaluation scores for each individual. Overall, the polymorphic SNPs discovered in this study can be utilized for the identification of powdery mildew resistance and susceptible cultivars in pepper breeding. However, in future the present investigation will be extended to evaluate large populations with more number of HRM primers corresponding to important SNPs associated with powdery mildew resistance in pepper.

In summary, the present endeavor reports the discovery of numerous SNP markers with potential applications in population genetics, molecular breeding, linkage mapping, and comparative genomics on gene-based association studies. For the first time, polymorphic SNPs were discovered from *C. annuum* and *C. baccatum* varieties of pepper with different powdery mildew resistance property. The SNP information obtained from the current WGRS approach in pepper can be utilized for the genomics assisted breeding of Capsicum with powdery mildew resistance.

Methods {#Sec10}
=======

Isolation of genomic DNA from pepper plants {#Sec11}
-------------------------------------------

Young leaves of PRH1 and Saengryeg were used for genomic DNA isolation. Briefly, 300 mg of leaves were ground into fine powder using liquid nitrogen. High-quality DNA was extracted using the cetyltrimethylammonium bromide (CTAB) extraction method^[@CR42]^. Powdered samples were mixed with CTAB buffer and incubated at 65 °C for 10 minutes. Sample mixtures were cooled to room temperature, and chloroform was then added to the sample mixture. Chloroform sample mixtures were mixed thoroughly and centrifuged at 13,000 rpm for 5 minutes at 4 °C. The supernatant was transferred into a new tube, and an equal volume of absolute ethanol was added. The solution was centrifuged at 13,000 rpm for 5 minutes at 4 °C, and the supernatant was discarded. Then, 70% ethanol was added to the sample, which was then centrifuged at 13,000 rpm for 5 minutes at 4 °C. Once again, the supernatant was discarded, and precipitated DNA pellets were dried at room temperature. The precipitated DNA pellets were then used as a starting material for purification using the Sigma Genelute plant DNA isolation kit (G2N70, Sigma). The DNA quality was assessed by electrophoresing the DNA on 1% agarose gel. The concentration of the extracted DNA was estimated using a GE Healthcare Bio-Science NanoVue via assessment of a single absorbance peak at 260 nm, a 260/280 absorbance ratio of 1.8 to 2.0 and no evidence of substantial band shearing or contamination (either RNA or polysaccharide).

DNA library construction and massively parallel sequencing {#Sec12}
----------------------------------------------------------

Purified whole genomic DNA was randomly sheared using a Covaris S2 (Covaris, Woburn, MA) to yield DNA fragments in the target range of 400 to 500 bp, and average molecular sizes were assessed using an Agilent Bioanalyzer 2100 (Agilent Technologies, Palo Alto, CA). Subsequently, the resulted overhangs were converted to blunt ends using a TruSeq DNA Sample Preparation Kit v2 (Illumina, CA, USA) followed by a clean-up protocol using AMPure XP Beads (Beckman Coulter Genomics, Danvers, MA). To enhance the ligation between the fragmented DNA and index adapters and to avoid self-ligation, the 3′ ends were adenylated. After adenylation, the index adapters were ligated to the fragmented genomic DNA, and the ligated products were purified using the AMPure XP Beads. The ligated products were size-selected on a 2% agarose gel followed by gel elution and column purification. The selected ligated DNA fragments with adapter sequences were enhanced through PCR using adapter-specific primers. Further, the DNA was re-isolated and the average molecular sizes of the libraries were evaluated using the Agilent Bioanalyzer 2100 (Agilent Technologies, Palo Alto, CA) to assess a sharp peak in the expected 500--600 bp range. Each library was loaded on the HiSeq. 4000 platform, and the high-throughput sequencing was performed to ensure that each sample met the 10-fold average sequencing depth.

Preprocessing {#Sec13}
-------------

After sequencing, the raw reads were trimmed using the Solexa QA v.1.13 package (Cox *et al*., 2010). The quality of bases from either end of Illumina reads commonly drop in, therefore either end of the reads were trimmed when the Phred quality score dropped below Q = 20 (or 0.05 probability of error). In addition, all 5′ and 3′ stretches of ambiguous 'N' nucleotides were also clipped. Trimming resulted in reads with a mean length of 101 bp across all samples, and a minimum length of 25 bp was applied during sequence trimming. These data were used for downstream analysis. The reference genome sequence of *Capsicum annuum* cv. CN334 was downloaded from Sol Genomic Network (SGN) at <http://www.sgn.cornell.edu/>.

Alignment, detection, and annotation of SNPs {#Sec14}
--------------------------------------------

To align the reads to the pepper reference genome, the Burrows-Wheeler Aligner (BWA 0.6.1-r104) program^[@CR43]^ was applied. The BWA default values for mapping were used, except for seed length (−l) = 30, maximum differences in the seed (−k) = 1, number of threads (−t) = 16, maximum number of gap extensions (−e) = 50, mismatch penalty (−M) = 6, gap open penalty (−O) = 15, and gap extension penalty (−E) = 8. Mapped reads were extracted from the resulting BAM file using SAMtools 0.1.16^[@CR44]^ for further analyses. The high mapping quality ensures reliable (unique) mapping of the reads, which is important for variant calling. Using the varFilter command, SNPs were called only for variable positions with a minimal mapping quality (−Q) of 30. The minimum and maximum of read depths were set as 3 and 100, respectively. An *in-house* script considering biallelic loci was used to select significant sites in the called SNP positions^[@CR31]^. Depending on the ratio of SNP reads to mapped reads, variant types were classified into three categories: homozygous SNP (more than 90%), heterozygous SNP (morethan 40% and less than 60%), and other SNPs for the remaining types. The polymorphic SNPs between two samples with sufficient sequences on both sides of the SNP site, without structural variation were noted adjacent to the SNP site and selected for primer design. To design primers flanking the SNP, an *in-house* script and Primer3 (v2.3.5) software were used^[@CR45]^. The parameters employed for the primer designing areas follows, primer length 18--24 bp, with 20 bp as the optimum; primer GC% = 20--80%, with the optimum value being 50%; primer Tm 55--65 °C, with 60 °C as the optimum; and product size range of 80--600 bp. After the designed primers were mapped to the genome sequence, only the primers that aligned were selected as candidates for SNP markers.

Functional annotation of genic SNPs {#Sec15}
-----------------------------------

The functional annotations of polymorphic SNPs were determined using the information acquired from gene ontology consortium ([www.geneontology.org](http://www.geneontology.org)) and Gene Ontology (UniProt) ([www.uniprot.org/help/gene_ontology](http://www.uniprot.org/help/gene_ontology)). The number of SNPs associated with each gene was identified manually.

Genotyping of SNPs using high resolution melt assay (HRM) {#Sec16}
---------------------------------------------------------

For the SNP validation, HRM primers were designed from each chromosome and evaluated in 46 F~4~individuals and compared with the parental lines. The HRM analyses were performed in 20 μl of total reaction mixture containing 2 μl of DNA extract (200 ng), 1× of SsoFastEvagreenSupermix (Bio-Rad Laboratories, Hercules, CA, USA), and 200 nM of forward and reverse primers. The reactions were performed in a fluorometric thermal cycler CFX96 real-time system (Bio-Rad Laboratories, Hercules, CA, USA), following program: 98 °C for 2 min, 45 cycles at 98 °C for 5 s and 60 °C for 10 s. The peaks obtained were normalized and analyzed for the difference in the melt curve.

Physiological disease resistance evaluation {#Sec17}
-------------------------------------------

The HRM results were correlated with the physiological evaluation of disease resistance. For the infection of powdery mildew, the parental lines used in this study, *C. annuum* - TF68 and *C. baccatum* -- ARI are the close relatives of PRH1 and Saengryeg. The parents as well as the F~4~ populations were maintained in a polyvinyl house along with disease infected plants under a normal day light condition with night/day set temperatures of 27/15 °C and 60--70% RH. The experiment was performed in triplicates in random block design. The disease severity has been assessed in 1--5 scale (1-resistant, 3-moderate and 5- sensitive) after two weeks.

Electronic supplementary material
=================================

 {#Sec18}

Supplementary Information. Supplementary Dataset 2 Supplementary Dataset 1

**Electronic supplementary material**

**Supplementary information** accompanies this paper at 10.1038/s41598-018-23279-5.

**Publisher\'s note:** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

This research work was supported by the Cooperative Research Program for Agriculture Science and Technology Development \[Project No. PJ012671022018\], Rural Development Administration, Republic of Korea.

Y.K.A., S.K., T.H.J. and E.Y.Y., designed and conceived the experiments; S.K., S.C., J.H.K., A.M., and E.S.L. performed the experiments. Y.K.A., A.M., S.K. and D.S.K. analyzed the data. A.M. and S.K. wrote the paper. Y.K.A. proofread and finalized the manuscript.

Competing Interests {#FPar1}
===================

The authors declare no competing interests.
