The DNA samples include 68 wild-caught unrelated chimpanzees (Becquet et al. 2007; Gonder et al. 2006, 2011), including Central African subspecies, P. t. troglodytes (15 individuals), Western African subspecies, P. t. verus (22 individuals), and the Gulf of Guinea subspecies, P. t. ellioti (31 individuals) (supplementary table S1, Supplementary Material online). We studied the genetic variation in the WFDC locus by Sanger sequencing the coding regions of 18 WFDC and SEMG genes (comprising a total of 66 exons) and a number of intervening noncoding regions (spaced every ∼10 kb). Additionally, we sequenced 47 pseudogenes located in unrelated, neutrally evolving regions across the chimpanzee genome, used as control regions, as previously described (Andrés et al. 2010; Ferreira et al. 2013) (supplementary table S2, Supplementary Material online).

Primers for amplification and sequencing were designed based on the Human Genome Reference Sequence (March 2006 assembly - v36.1), available at the UCSC Genome Browser (genome.ucsc.edu, last accessed December 24, 2013). All samples were polymerase chain reaction (PCR)-amplified and analyzed by bidirectional Sanger sequencing. Further details about PCR and DNA sequencing are available from the authors upon request. The sequences were aligned to the Human Genome v26.1 and polymorphic sites and fixed differences were detected with phred-phrap-consed package (Nickerson et al. 1997). To ensure sequence quality, we discarded variant sites in the first and last 75 bp of each amplicon segment. We manually curated sites found to have discordant genotypes in different amplicons. The ancestral state of each SNP was inferred by comparison with the human, orangutan, and macaque genome sequences (Gibbs et al. 2007; Andrés et al. 2010; genome.ucsc.edu).

We assessed the subspecies differentiation levels by calculating the population differentiation (F_ST statistic) and by performing principal component analysis (PCA) on the SNP data (Excoffier 2002; Patterson et al. 2006). We used a locus-by-locus analysis of molecular variance, using 20,000 simulations, which was performed by Arlequin using its default values (constant model; Excoffier et al. 2005). The EIGENSOFT software package was used for PCA (Patterson et al. 2006). We performed cluster analysis using STRUCTURE version 2.3 (Pritchard et al. 2000), assuming admixture and correlated allele frequencies. Fifty iterations of the data at each K = 1–5 with 500,000 Markov chain Monte Carlo (MCMC) burn-in steps and 500,000 MCMC iterations. STRUCTURE output was processed with CLUMPP and plotted with DISTRUCT (Conrad et al. 2006). We used STRUCTURE harvester to determine the best K estimate. Population structure analyses were performed blinded to a priori population labels.

A model of chimpanzee demography was inferred using BP&P program which implements Bayesian inference with a MCMC and accommodates a species phylogeny as well as lineage sorting due to ancestral polymorphism (Rannala and Yang 2003; Yang and Rannala 2010). We used a proposed phylogenetic tree of the three chimpanzee subspecies (Gonder et al. 2011) and included human as the outgroup. A gamma prior G(2, 1,000), with mean 2/2,000 = 0.001, was used on the population size parameters (θ_s), and the age of the root in the species tree (τ₀) was assigned to the gamma prior G(25, 5,000). This was based on the assumption of divergence time between chimpanzee and humans of 5 Myr and a mutation rate of 10⁻⁹ per site/per year (2 × 10⁻⁸ per site per 20-year generation). The other divergence time parameters were assigned to the Dirichlet prior (Yang and Rannala, 2010). The analyses were run twice to confirm consistency, and parameters of historical demographic events are expressed as the mean and 95% confidence intervals of the posterior distributions.

Summary statistics of population genetic variation were calculated using SLIDER (http://genapps.uchicago.edu/slider/index.html, last assessed December 24, 2013). We assessed statistical significance of summary statistics using an empirical comparison to the control regions, by calculating the upper and lower 2.5 percentiles of each distribution. Specifically, we used the sequenced control regions to perform an empirical comparison of nucleotide diversity (π) and Tajima’s D values for each WFDC gene in each population. For the WFDC genes, we ran 10⁵ coalescent simulations using “ms” (Hudson 2002) and mutation rate parameters estimated from the sequenced data with SLIDER. For the population recombination parameter, we used the PANMAP estimates of chimpanzee recombination. We assumed demographic models that included constant population size, and historic events as previously described (Hey 2010; Wegmann and Excoffier 2010), and as inferred by us for each subspecies. For every model, we calculated a null distribution of summary statistics values and calculated the 2.5th and 97.5th percentiles.

We performed HKA tests considering all subspecies in DNAsp 5.1 and using a maximum-likelihood method that incorporates values for multiple neutrally evolving regions (Hudson et al. 1987). McDonald–Kreitman test (MKT) was calculated in DNAsp v5.1 using humans as the outgroup and assuming two types of sites: putatively neutral sites (Syn) and functional sites (NSyn) (Rozas and Rozas 1995; Rozas 2009).

Haplotype phasing for all samples was inferred separately for the WFDC centromeric and telomeric subloci (WFDC-CEN and WFDC-TEL, respectively; see fig. 1 for WFDC locus substructure) using PHASE2.1 (Stephens et al. 2001; Stephens and Donnelly 2003). Haplotypes were independently inputted in Haploview 4.2 (Barrett 2009) to calculate linkage disequilibrium (LD) statistics, r² and D′, and to identify LD and haplotype blocks (Gabriel et al. 2002). The potential functional effects at the protein level of non-synonymous (NSyn) SNPs and fixed differences were inferred using PolyPhen v2 (Adzhubei et al. 2010) and SIFT (Kumar et al. 2009).

Maximum-likelihood estimates of d_N/d_S (ω; d_S – synonymous substitution rate and d_N – non-synonymous substitution rate) were carried out using the codeml program from the software package Phylogenetic Analysis by Maximum Likelihood – PAML version 4.2 (Yang 2007). To run PAML, we first reconstructed a phylogenetic tree (DNAml from Phylogeny Inference Package [PHYLIP]; http://evolution.genetics.washington.edu/phylip.html, last accessed December 24, 2013). We used the genomic sequences from human (Homo sapiens), chimpanzee (Pan troglodytes), gorilla (Gorilla gorilla), orangutan (Pongo pygmaeus), gibbon (Nomascus leucogenys), rhesus monkey (Macaca mulatta), baboon (Papio anubis), and marmoset (Callithrix jacchus). They were retrieved from public databases using EPPIN isoform 1 (Uniprot: O95925) and WFDC6 isoform 1 (Uniprot: Q9BQY6) as BLAT templates. Pan paniscus was not included in the analysis as the cDNA sequence is equal to P. troglodytes. The phylogenetic tree diverged from the known primate phylogeny in the position of orangutan and gibbon branches. To test for variable selection pressures among branches, we performed the branch model using either the null model (one ω ratio for the entire tree) or nested models (two-ratios, three-ratios, four-ratios for the tree) (Yang 1997; Bielawski and Yang 2003). The values of ω > 1 were considered as evidences of positive selection, and the values ω < 1 were considered as an indication of purifying selection. The test statistic was constructed as twice the difference in the log of the likelihoods (−2Δl), and significance was assessed by comparing this to the χ² statistic.

We first characterized the levels of differentiation between subspecies at the control regions using F_ST (Excoffier 2002) and PCA (Patterson et al. 2006). The pairwise F_ST values show lower differentiation between P. t. troglodytes and P. t. ellioti (0.104) than each of them compared with P. t. verus (0.392 and 0.400, respectively). The first two principal components (PC1 and PC2) separate P. t. verus from the other two subspecies, and the third principal component appears to not separate completely P. t. troglodytes from P. t. ellioti (supplementary fig. S2, Supplementary Material online). Contrary to previous studies, we could not separate the three subspecies using PCA (Gonder et al. 2011; Bowden et al. 2012). We further examined the shared ancestry levels of the individuals by performing a Bayesian model-based clustering approach available in the STRUCTURE software (Pritchard et al. 2000; Falush et al. 2003). The analysis was performed blinded to a population label and two groups (K = 2) were recovered (supplementary fig. S3, Supplementary Material online). The subspecies P. t. troglodytes and P. t. ellioti could not be confidently distinguished even after the inclusion of the 849 SNPs from the WFDC locus (results not shown).

We inferred the demographic history of the three chimpanzee subspecies from control regions data using a Bayesian MCMC approach, based on the phylogenetic tree proposed by previous studies (Gonder et al. 2011; Bowden et al. 2012). The estimated effective population sizes were ∼52,000 and ∼21,000 for P. t. troglodytes and P. t. verus, respectively (table 1). These values are within the range of previous models of chimpanzee demography as indicated by the overlaps between confidence intervals (table 1), and the 0.31 Myr estimate for divergence from P. t. verus is in the same time frame from previous studies (table 1; Becquet et al. 2007; Hey 2010; Wegmann and Excoffier 2010; Gonder et al. 2011). For P. t. ellioti, we estimated a recent divergence from P. t. troglodytes around 0.173 Ma and an effective size of ∼43,000 individuals, demonstrating values that were of similar order of magnitude to P. t. troglodytes. These two subspecies appear to have preserved a similar effective population size to their ancestral population (fig. 2 and table 1). Conversely, the origin of P. t. verus is associated with a bottleneck and effective population size reduction (fig. 2), providing a good fit to previous models of chimpanzee demography (Becquet et al. 2007; Hey 2010; Wegmann and Excoffier 2010; Gonder et al. 2011). Fig. 2.—

Schematic representation of the inferred demographic history of the three subspecies: Pan troglodytes troglodytes (Ptt); P. t. verus (Ptv); and P. t. ellioti (Pte). NA, ancestral effective population size; N, effective population size; TDIV, divergence time.

The patterns of variation in SEMG1 among individuals include a polymorphic highly repetitive region (9–13 modules encoded by SEMG1 exon 2) (Jensen-Seaman and Li 2003). Because the modular nature of SEMG1 precluded a consistent and unambiguous sequence alignment, the SNPs located in this repetitive region were removed from the analysis for quality purposes. A total of 766 fixed differences were identified in comparison with the human genome reference. Twenty-five indels (insertions/deletions) were found, 24 of which were located in introns, UTRs, and intergenic regions. One remaining indel located in WFDC6 was identified in a single chromosome (f = 0.02). Indels were excluded from all analyses, due to their distinct mutation rate and the low overlap with functional regions, which are unlikely to affect protein function or expression. Additionally, from the 456 human–chimpanzee fixed differences located in the WFDC locus, only 19 were within coding regions and chimpanzee-specific. Of these, 17 were non-synonymous (NSyn) and most of them were classified as benign by Polyphen v2 and SIFT (supplementary table S3, Supplementary Material online).

To characterize the within-subspecies variation, we analyzed the folded site frequency spectrum (SFS) for all SNPs and Syn and NSyn sites among WFDC genes (fig. 3). Additionally, we analyzed the deleterious effects of each coding substitution using SIFT and Polyphen (Kumar et al. 2009; Adzhubei et al. 2010; supplementary table S3, Supplementary Material online). Despite the higher number of NSyn sites in the WFDC genes, these are maintained at low frequencies in the overall species, consistent with the predicted mildly deleterious effects (supplementary table S3, Supplementary Material online). For each WFDC gene, we calculated summary statistics such as nucleotide diversity (π), Tajima’s D (Tajima 1989), Fu and Li’s D (Fu and Li 1993), Fay and Wu’s H (Fay and Wu 2005; Zeng et al. 2006), and HKA (Hudson et al. 1987) (table 2). Fig. 3.—

Folded SFS for the species that were resequenced. The x axis depicts the frequency of the allele frequency bin in the generated data set, whereas the y axis represents the number of alleles found within each frequency bin. Syn, synonymous changes; NSyn, nonsynonymous changes. (A) Folded SFS in WFDC locus; (B) folded SFS of WFDC locus highlighting coding mutations.

To determine whether specific WFDC genes have been under selective pressures in one or all chimpanzee subspecies, we started by comparing the summary statistics for each WFDC gene with the empirical distribution of Tajima’s D in the control regions (fig. 4). Only WFDC6 and EPPIN show unusual patterns in P. t. troglodytes (fig. 4). Although in this subspecies the allele frequency spectrum is generally skewed toward rare alleles, the Tajima’s D values of WFDC6 (−2.073; P value = 1⁻⁴) and EPPIN (−1.811; P value = 0.025) present the lowest values of control and WFDC regions. WFDC6 also presented a low Tajima’s D value (−2.1039) and significant HKA (P = 0.013) when combining all individuals sequenced (supplementary table S4, Supplementary Material online). The other WFDC genes did not show strong significant P values pointing to a neutral evolution based on subspecies genetic diversity (table 2). To confirm WFDC departures from neutrality, we performed 10⁵ coalescent simulations for each subspecies under different demographic scenarios: constant model (P. t. troglodytes, P. t. ellioti, P. t. verus), our best-fit model (P. t. troglodytes, P. t. ellioti, P. t. verus), Hey 2010 model (P. t. troglodytes), and Wegmann and Excoffier 2010 model (P. t. verus). WFDC6 and EPPIN present significantly negative Tajima’s D value compared with all models, but while WFDC6 shows always values below 1st percentile, EPPIN values lie between the 1st and 2.5th percentile (table 2). Fig. 4.—

Empirical comparisons generated from the 47 control regions. Tajima’s D (Tajima 1989) was calculated for each region using SLIDER and plotted with the 2.5 and 97.5 percentiles represented as dashed lines.