A representative sample of 377 pepper accessions (283 C. baccatum and 94 diverse C. annuum accessions) collected from 32 countries across the world were obtained from the USDA-ARS, Germplasm Resource Information Network, Plant Genetic Resources Conservation Unit, Griffin, GA and World Vegetable Center (AVRDC, Shanhua, Taiwan) (Table S1). The C. annuum collection was comprised of 90 domesticated cultivars and 4 wild accessions. The C. baccatum collection had 218 lines of C. baccatum var. pendulum and 17 wild accessions (C. baccatum var. baccatum). Peduncle length (cm) was measured for 5 plants each of 217 accessions belonging to C. baccatum var. pendulum grown in a greenhouse in three replications.

Genomic DNA was isolated from the seedlings using the DNeasy plant mini kit (QIAGEN, Germany), and GBS was as described (Elshire et al., 2011). DNA was treated with the restriction enzyme ApeKI, a type II restriction endonuclease, barcoded by accession, and sequenced on an Illumina HiSeq 2500 as described (Elshire et al., 2011). SNPs were identified using the TASSEL-GBS Discovery/Production pipeline (https://bitbucket.org/tasseladmin/tassel-5-source/wiki/Tassel5GBSv2Pipeline). Chromosomal assignment and position on the physical map of various SNPs were deduced from the C. annuum whole genome sequence (WGS) draft at http://peppergenome.snu.ac.kr. SNPs were designated by chromosome number and position (e.g., S10_172735351, which indicates an SNP located at position 172735351 on chromosome 10).

Genetic diversity values were calculated by a neighbor-joining algorithm using TASSEL 5. In a second approach, we utilized IBS and principle component analysis (PCA) with the SNP & Variation Suite (SVS v8.1.5) (Golden Helix, Inc., Bozeman, MT, USA; www.goldenhelix.com). Observed nucleotide diversity (π) and Tajima's D were estimated by using TASSEL v5.0 with a sliding-window approach as described (Korneliussen et al., 2013). The fixation index (F_ST) was estimated on the basis of the Wright F statistic (Weir and Cockerham, 1984) with use of SVS v8.1.5.

For GBS data, we considered only SNPs successfully mapped to the C. annuum WGS draft, because knowing the chromosome location of SNPs helps prevent spurious LD and thereby unreliable association mapping. Mapped SNPs were further filtered by call rate >90%. Before studying LD decay, haplotype blocks were calculated for all markers by using the default settings in SVS v8.1.5. Adjacent and pairwise measurements of LD for GBS data were calculated separately for SNPs in each chromosome. For computing LD, we used the expectation-maximization (EM) algorithm (Dempster et al., 1977) as an iterative technique for obtaining maximum likelihood estimates of sample haplotype frequencies.

The PC matrix was constructed with the program “EIGENSTRAT” (http://genetics.med.harvard.edu/reich/Reich_Lab/) and the PCA correction technique; the method of stratification was as described (Price et al., 2006). IBS was calculated as described (Purcell et al., 2007). GWAS involved a single-locus mixed linear model (SLMM), a method that uses a forward and backward stepwise approach to select markers as fixed-effects covariates in the model (Segura et al., 2012), and implemented in SVS v8.1.5. We used a PC matrix to correct for population stratification and an IBS matrix to correct for a polygenic background. Manhattan plots for associated SNPs were visualized by using GenomeBrowse v1.0 (Golden Helix, Inc.). The SNP P-values from GWAS underwent false discovery rate (FDR) analysis (Storey, 2002).

A total of 77,407 SNPs were isolated from the nucleotide sequence obtained for the 283 C. baccatum and 94 C. annuum accessions studied; a total of 8661, 8086, 9843, 6197, 5688, 7410, 5588, 5086, 4472, 5336, 5079, and 5961 SNPs were mapped to the WGS draft and located on chromosomes 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, and 12, respectively. We noted the presence of one SNP at every 35.6 kb across the genome, with average gap size of 31.7 kb and one SNP at every 104.4 kb in the coding regions. A total of 36,621 SNPs had minor allele frequency [MAF] ≥0.05, identified collectively for C. annuum and C. baccatum, and were used for various analyses in the current study. For C. baccatum, 13,129 SNPs had MAF ≥0.05; their chromosome distribution is listed in Table 1. In addition, we identified 26,697 SNPs located in various exons. SNP counts in exons of various genes were 2985, 3308, 3630, 2032, 1837, 2474, 1897, 1758, 1406, 1799, 1550, and 2021 on chromosomes 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, and 12, respectively.

We used PCA of the 36,621 SNPs identified from C. baccatum and C. annuum with MAF ≥0.05 to characterize domesticated and wild C. annuum and C. baccatum peppers. PCA with first and second eigen vectors that explained 80% of the total variation produced two clusters of C. baccatum and C. annuum accessions (Figure 1). Tepin and Tepin Guatemala, two wild peppers belonging to C. annuum var. glabriusculum that are native to southern North America and northern South America, were close to CB-77, a wild C. baccatum pepper. Similarly, three other wild C. baccatum peppers, CB-93, CB-92, and CB-40, were intermediate between the major C. annuum and C. baccatum clusters. A third cluster comprised the remaining wild, semi-domesticated and crown shaped fruit type C. baccatum accessions that grouped with the domesticated large-fruited C. baccatum peppers. A separate PCA with 13,129 SNPs that were polymorphic for C. baccatum accessions resolved the population structure comprised by this group of C. baccatum accessions. This PCA identified 283 C. baccatum accessions in 3 distinct clusters (Figure 2). The middle cluster (cluster II) was parallel to the C. annuum cluster, and the wild species Tepin, Tepin Guatemala, CB-77, CB-93, CB-92, and CB-40 were found in the middle, which indicates intercrossing between wild C. annuum and C. baccatum peppers while or before domestication. PCA placement of various accessions are noted in Tables S2, S3.

F_ST was estimated with 95% confidence intervals between wild and domesticated C. annuum and C. baccatum. The F_ST between wild (C. annuum + C. baccatum) and domesticated (C. annuum + C. baccatum) accessions was 0.09 and 0.05, respectively. The F_STbetween domesticated C. annuum and C. baccatum was 0.78, which indicates genome-wide divergence. The F_ST between wild C. baccatum and wild C. annuum was 0.66. Crown-shaped fruited C. baccatum types are unique for this species group, and pairwise F_ST values with wild, semi-domesticated and domesticated were 0.10, 0.06, and 0.03, respectively, which indicates their closeness to domesticated types. F_ST-values for semi-domesticated with wild and domesticated C. baccatum types were 0.03 and 0.01, respectively. We present an overall F_ST distribution in a Manhattan plot for all chromosomes showing important chromosomal regions with the highest F_ST as peaks (Figure 3, Table S4). Based on F_ST values, peaks on chromosomes 1, 2, 3, 4, 5, 6, and 9 in the Manhattan plot might be the regions of positive selection and important for improvement.

Because of the strong population structure, we assessed patterns of variation separately for each group of domesticated accessions from the respective species when making inferences about the evolutionary dynamics of domestication. Crop domestication is often associated with “population bottlenecks” because of the limited number of founding individuals experiencing domestication events. These bottlenecks may be evident in pepper when comparing diversity between cultivated forms of C. annuum and C. baccatum. We estimated nucleotide diversity (π) and Tajima's D across various chromosomes to understand genome-wide bottleneck effects. The frequency of segregating SNPs as reflected by various chromosomal measures of mean π and Tajima's D is presented in Figure 4. For cultivated C. baccatum, chromosome 4 was positive for π and Tajima's D which indicates accumulation of rapid mutations on this chromosome. The remaining chromosomes were negative or nearly negative for Tajima's D, which indicates bottlenecks in domestication. In contrast, C. annuum chromosomes were positive for Tajima's D on all chromosomes except chromosome 8, which indicates differential evolution after the domestication or the influence of diverse breeding.

We conducted an extensive LD analysis on the entire dataset of C. baccatum accessions on all adjacent marker pairs within a chromosome or within a haplotype block. Haplotype distribution is important to understand patterns of genetic variation of C. baccatum gene pools and has a wide range of applications. The 2 major processes that shape haplotype structure are the domestication process and breeding history. We used “minimize historical recombination,” a block-defining algorithm developed by Gabriel et al. (2002). The upper confidence boundary was set to 0.98 and the lower boundary to 0.70. SNPs with MAF <0.05 were omitted. Maximum block length was set to 160 kb. The expectation maximization (EM) algorithm was used for haplotype estimation, with convergence tolerance 0.0001, and frequency threshold 0.01. Maximum EM iterations were set to 50. We identified 1742 haplotypes containing 4420 SNPs, with a range of 9–2 SNPs per haplotype (Table S5). The results provided values for both the EM algorithm (Dempster et al., 1977) and composite haplotype method (CHM; Weir and Cockerham, 1996). Squared-allele frequency correlations (r²) and LD estimate (D′) for the EM and CHM methods are in Table S6. We created LD plots by using marker-pair associations of adjacent SNPs within a chromosome, within a haplotype block, and within genes (Figure 5). The length of individual LD blocks varied among chromosomes, with regions of high and low LD interspersed (Table 2). The genome-wide average LD block was 99.1 kb. The largest LD block, of 13,021 kb, was on chromosome 11. Pairwise LD was estimated by r² and we compared the pattern of decay at different levels. With pair-wise analysis considering adjacent SNPs across chromosomes, most SNP associations were within 50 kb (Figure 5). The second analysis based on adjacent SNPs within haplotypes revealed most associations within 20 kb (Figure S1). The third analysis of SNPs located in genes revealed most associations within 5 kb (Figure S2).

Peduncle length is the prime differentiating trait between wild and domesticated forms of C. baccatum. Mean peduncle lengths for respective accessions are listed in Table S7. The cultivated form of C. baccatum, var. pendulum, is named based on the epithet related to pendant fruits. In our GWAS, 36 SNPs located on chromosomes 1, 2, 3, 4, 6, 7, 8, 9, 10, and 11 were identified as significantly associated with peduncle length and cumulatively explained 21% of the total variation (Figure 6). Four SNPs located in the intergenic space between the oxidoreductase family protein/arogenate dehydrogenase on chromosome 7 explained 10.6% of the total variation. Chromosome number, map position, P-value, regression beta, FDR correction, variance explained, call rate, and minor/major allele frequencies for all significantly associated SNPs are in Table S8.

The predicted gene set from the annotated C. annuum cv. CM334 reference genome (Kim et al., 2014) was used to characterize the genes containing SNPs or nearby SNPs. Eleven candidate genes containing SNPs in exons or promoters were significantly associated with peduncle length, and 12 more SNPs in introns or intergenic regions of candidate genes were proposed. GWAS details and strengths of association of SNPs are in Table S8. Details of annotation for various associated SNPs, their location in various genes and type of mutation (synonymous or non-synonymous) are in Table 3.

Wild C. baccatum has a relatively restricted distribution confined to southern Peru, Bolivia, and southern Brazil (Eshbaugh, 1970). C. baccatum var. pendulum is a widely distributed cultivated plant found throughout western South America and now spreading worldwide (Eshbaugh, 1970). Wild C. baccatum has red, erect, and non-persistent fruits, and C. baccatum var. pendulum has red, orange, yellow, green, or brown fruits that are pendant and persistent. Because the peduncle is the most differentiating trait between domesticated and wild C. baccatum species, we performed GWAS for peduncle length. We associated 36 SNPs with the trait peduncle. Four of these SNPs clustered with candidate genes on chromosome 7. Annotation for some of these associated SNP-containing sequences revealed their location in various genes, so these genes might play a role in peduncle length, peduncle architecture and C. baccatum domestication.

Length of peduncle is determined by the cell number or cell size, although it is indirectly regulated by hormones and multiple pathways. Kinases play important roles in plant growth and development. Peduncle associated SNPs in the current study were located in leucine-rich repeat receptor like kinases (LRR-RLKs), serine/threonine protein kinase, ABC transporter gene and RING finger protein, which may play important roles in growth and development as well as cell wall integrity and elongation as has been shown in other plants (Lally et al., 2001; Arunyawat et al., 2007; Guo et al., 2009; Gish and Clark, 2011; Ghosh et al., 2013). Plant cell walls contain a glycoprotein component rich in the otherwise rare amino acid hydroxyproline and accumulation of this amino acid was positively correlated with cell elongation in pea epicotyls (Flint-Garcia et al., 2003). In the current study, we also associated a marker S11_725918 on GABA (γ-aminobutyric acid), a ubiquitous non-protein amino acid. An Arabidopsis GABA gene mutant pop2 exhibited defects in hypocotyl cell elongation and pollen-tube elongation via influence on cell-wall–related genes (Bulmer, 1971).

Our study describes the utility of SNPs generated by GBS for genome-wide divergence and LD patterns between C. annuum and C. baccatum. Mapping all the SNPs to the C. annuum reference genome helped to identify homologous SNPs between the two incompatible cultivated pepper genomes, which was further useful to reduce ascertainment bias, so this SNP set was useful in estimating genome-wide population differentiation and allele sharing between the two genomes. Furthermore, the SNPs anchored to the C. annuum genome may not be in the same order in the C. baccatum genome because some genomic regions may not be co-linear to the C. annuum genome because of genome rearrangements. In a comparison of C. baccatum and C. annuum linkage maps, Lee et al. (2016) identified two major reciprocal translocations between chromosomes 3 and 5 and between chromosomes 3 and 9, as well as translocations between chromosomes 1 and 8.

Such uncertain positions of SNPs can be corrected only when the whole genomesequence is available for C. baccatum genome. This SNP panel and the results pertaining to population structure, IBS and LD decay analyses will facilitate routine use of GWAS for identification of genes associated with various economically important traits in Peruvian peppers. Our identification of SNPs associated with fruit peduncle length demonstrates opportunities for utilization of GWAS in crop improvement.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.