Archaic Hominin Introgression in Africa Contributes to Functional Salivary MUC7 Genetic Variation

Copy Number Variation of MUC7 PTS Repeats Evolved Recurrently in Humans

Discovering and genotyping the copy number of short repeat segments using short read sequences is challenging due to alignment and assembly errors (Mills etal. 2011; Zhao etal. 2016). Indeed, the 1,000 Genomes Project data did not capture the copy number variation affecting MUC7 (Sudmant etal. 2015). To genotype the full range of PTS-repeat copy numbers, we directly genotyped 251 randomly selected samples from European, African, and Asian populations already sequenced by the 1,000 Genomes Project (supplementary table S1, Supplementary Material online). Then, to understand the haplotypic variation of MUC7 subexonic copy number alleles, we constructed a maximum likelihood tree using the nonrepeated sections of MUC7 from 5,008 phased haplotypes out of the 1,000 Genomes Project Phase 3 (fig. 2a, supplementary fig. S1, Supplementary Material online). The haplotypes fell into eight divergent haplogroup clusters (MUC7*A-H; fig. 2b). When we superimposed the number of PTS-repeat copies and the corresponding haplotypes constructed by single nucleotide variants, we observed that two haplogroups (groups G and E) harbor 5 PTS-repeat alleles, whereas the other groups harbor 6 PTS-repeat alleles (supplementary table S1, Supplementary Material online). This suggests that MUC7 PTS-repeat copy number variation evolved independently in at least two different branches of the human phylogenetic tree, indicating recurrent evolution of these variants.

To test whether recurrence indeed shaped MUC7’s genetic variation, we directly sequenced five-copy haplotypes from haplogroups G and E, as well as six-copy haplotypes. To do this, we amplified the PTS-repeat region in heterozygote individuals, separated five- and six-repeat haplotypes from each other by DNA gel electrophoresis and excised bands containing the five- and six-copy haplotypes from the gel for sequencing. We discovered that the repeat sequences of haplogroup G were distinct from haplogroup E, even though they both contain five repeat copies (fig. 2c). Specifically, compared with the human reference genomes, which all have the six PTS-repeat allele haplotype, each repeat has distinctive single nucleotide differences from each other. As such, it was possible to distinguish individual repeats. We found that haplogroup G is missing the fourth PTS-repeat, whereas haplogroup E is missing the fifth PTS-repeat, counted from the 5′ end. Assuming that the five-copy PTS-repeat was ancestral as we suggested in our earlier work (Xu etal. 2016), our results suggest at least three independent mutational events: 1) the duplication event creating the six-copy PTS-repeat must have occurred before haplogroup E diverged from the other human MUC7 haplotypes; 2) haplogroup G has lost the fourth copy; and 3) haplogroup E lost the fifth copy. Independent of the exact scenario, it is clear that two different five-copy PTS-repeat alleles are present in two independent haplotypic backgrounds, which strongly suggests recurrent formation of these alleles and a high mutation rate in this locus. Overall, we conclude that the contemporary variation of MUC7 PTS-repeat copy number has evolved through multiple independent mutational events happening since the divergence of humans and chimpanzees.

Imputation of MUC7 PTS-Repeat Haplotypes Allows Reevaluation of the Prior Association with Asthma Susceptibility

Accurate imputation of copy number variants is essential to incorporate these variants to genome-wide association studies conducted by single nucleotide genotyping platforms. As previously noted, the MUC7 PTS-repeat copy number was associated with asthma susceptibility in three independent locus-specific association studies (Kirkbride etal. 2001; Rousseau etal. 2006; Watson etal. 2009). However, due to the recurrent evolution of MUC7 PTS-repeat copy number variation, the individual single nucleotide variants do not accurately predict the copy number state of PTS-repeats. For example, the haplogroup E, which we showed to carry five MUC7 PTS-repeats, was not tagged by GWAS. As such, GWAS do not have the power to accurately tag MUC7 PTS-repeats. To address this, we identified multiple single nucleotide variants that tag the two independently evolved five PTS-repeat alleles with very high accuracy in all populations tested (supplementary table S1, Supplementary Material online). Combining the allele frequencies of these tag single nucleotide variants would theoretically allow evaluation of PTS-repeat copy number in existing single nucleotide variant based genotyping studies. However, when we attempted to use these single nucleotide variants to impute the PTS-repeat copy numbers in multiple GWAS studies for asthma (Anon 1999; Boushey etal. 2005; Sorkness etal. 2007; Moore etal. 2010), we realized that some of the tag variants were not included in the genotyping platform that was used in these studies. To address this shortcoming, we used only single nucleotide variants that were used in the GWAS study to predict the PTS-repeat copy number states of 1,000 Genomes data. This approach allowed us to find combinations of single nucleotide variants genotyped in the GWAS platform to accurately predict the copy number status of MUC7 PTS-repeats (supplementary fig. S2, Supplementary Material online). A similar approach was recently used for imputing structural variants within the haptoglobin locus (Boettger etal. 2016). This analysis showed that it is possible to tag the PTS-repeat copy number state (and the haplogroups) of MUC7 using only the data available in GWAS studies with high accuracy (99.8%, supplementary table S2, Supplementary Material online). Using this methodology, we were unable to replicate the association between MUC7 PTS-repeat copy number and asthma prevalence (nominal, locus-specific P-value = 0.2402 for African Americans, P-value = 0.1507 for European descendants, Cochran–Armitage test for trend test; supplementary table S2, Supplementary Material online). Also, the directionality of five-copy frequency in cases and controls are opposite in Africans and European cohorts, where African controls having a higher frequency of five-copy alleles as compared with the cases (supplementary fig. S2, Supplementary Material online). It is plausible that there is an indirect, weak correlation that depends on the external factors (e.g., the abundance of allergens in the environment). A more parsimonious explanation however is that the previous, locus-specific studies found spurious associations between MUC7 PTS-repeat copy number and asthma prevalence.

MUC7 Genetic Variation in Humans Is Associated with Oral Microbiome Composition

MUC7 is known to be recognized and bound by bacteria residing in the oral cavity as well as by systemic pathogens that occasionally traverse the mouth environment (Walz etal. 2009; Heo etal. 2013; Thamadilok etal. 2016). Thus, we hypothesized that genetic variation in MUC7 could affect bacterial colonization in the oral cavity. To test this hypothesis, we conducted a locus-specific association study, by analyzing genetic variation within 50 kb upstream and 50 kb downstream of MUC7 from 93 humans and correlated it with the microbial composition in fifteen body sites of the same individuals (Human Microbiome Project Consortium 2012; Blekhman etal. 2015). From all these 15 different body sites ranging from stool to skin, we found significant correlations only in the oral cavity (fig. 3a, supplementary table S3, Supplementary Material online; P-value ≤ 9.66 × 10⁻⁵ in each of the oral cavity sites). More specifically, we found 20 significant associations between genetic variants and microbial taxon abundance in the anterior nares, palatine tonsils, saliva, supragingival plaque, and tongue dorsum. For example, we highlighted the strongest signal for the association between single nucleotide variant rs12498483 and the abundance of the genus Eubacterium in the palatine tonsil (fig. 3b). Other associations include bacterial genera, such as Lautropia and Neisseria, which are known colonizers of mucosal surfaces and the oral cavity (Rossmann etal. 1998; Mager etal. 2003).

To understand haplotype-level variation that defines these genetic-microbial associations around the MUC7 locus, we investigated linkage disequilibrium patterns among the variants significantly associated with microbiome composition. We found one haplotype, defined by rs12498483, that fully includes MUC7 genetic variation (supplementary fig. S3, Supplementary Material online), and is made up of variants that broadly indicate samples that do not belong to haplogroups A, B, and C (fig. 2a). The haplotype is very common (>30% in most populations), reaching close to 50% in African populations. It is linked with Eubacterium, which has been reported to be a prominent and highly diverse genus in the oral cavity (Downes etal. 2001). Even though the host genotype data that we use in our microbiome association analysis comprise a limited number of variants, we were able to further identify dozens of single nucleotide variants that belong to this haplotype using 1,000 Genomes data set (supplementary fig. S3, Supplementary Material online). These include putatively functional variants, including seven that affect coding sequence (supplementary table S3, Supplementary Material online).

Our microbiome analysis primarily involved individuals with Eurasian ancestry. In addition, the bacterial associations are at the genus and occasionally class level, limiting the routes to direct functional testing. As such, our results are by no means definitive, and an adaptive impact of MUC7-microbiome interactions has yet to be resolved. However, our results suggest that the genetic variation within MUC7 influences the oral microbiome. These results, together with the fact that microbiome composition varies across human populations (Yatsunenko etal. 2012; Morton etal. 2015; Gomez etal. 2016), motivate the need to design more focused studies to tackle specific gene–microbe interactions in diverse human groups and generate novel hypotheses with regard to the adaptive role of MUC7 genetic variation.

Haplogroup E Is Unusually Divergent

During our analyses, we noticed an apparently divergent haplotype group, which is also one of the haplotypes carrying five PTS-repeat allele (haplogroup E, fig. 2a and b). We asked whether this haplotype is unusually divergent as compared with genome-wide expectations. To investigate this, we devised a statistical test inspired by Wall’s B (Wall 1999), which measures the maximum number of SNPs in a locus that are in perfect LD (r² = 1.0) between themselves. Similar to Wall’s B, our statistical method obtains high values if nonoverlapping subgroups of the sample are sufficiently distinct genetically. Wall’s B measures the number of single nucleotide variants that are in perfect LD with each other and are subsequent (i.e., no other variant in between). Our statistical approach is different from Wall’s B in that it measures the maximum number of single nucleotide variants that are in perfect LD with each other within a given locus independent of their location within that locus. We also normalize this measure by the total number of single nucleotide variants in that locus (see Materials and Methods for details). Thus, compared with Wall’s B, our approach is more robust to noise generated by rare variants or sequencing errors. Our results showed a significantly higher than expected number of variants that are in linkage disequilibrium with each other in the MUC7 locus compared with 1,000 randomly selected regions of the genome for African (YRI) populations (P-value = 0.050, Wilcoxon rank-sum test), but not for Eurasian populations (P-values > 0.1, Wilcoxon rank-sum test, supplementary table S4, Supplementary Material online, fig. 4a). We also confirmed this observation using other statistics such as theta (θ) and Wall’s B (supplementary fig. S4, Supplementary Material online). These results support that haplogroup E indeed constitutes a distinct clade separated by an unexpectedly long branch in the gene tree.

The unusually high divergence of haplogroup E from other haplogroups suggests an increased coalescence time to the most recent common ancestor in the MUC7 locus. To investigate the coalescent depth of the haplogroups, we used the BEAST software (Drummond and Rambaut 2007) to construct a phylogeny of two haplotypes from each haplogroup described above, including the chimpanzee and rhesus macaque reference haplotypes as outgroups (fig. 4b). This allowed us to estimate the most likely coalescence time for haplogroup E and the other human MUC7 haplotypes to ∼4.5 million years before present (95% confidence interval ranges from 3.2 to 6.3 million years before present). This indicates that for MUC7 both tree depth among humans and divergence between humans and chimpanzees date back further than most other parts of the genome (Schiffels and Durbin 2014). Our analysis also confirms that the depth of the tree is primarily due to the divergent haplotype group E.

On the basis of the coalescent tree, we noticed that haplogroup E originated prior to the assumed Human–Neanderthal population divergence, previously dated to between 260 and 765 thousand years before present (Prüfer etal. 2014). To quantify this, we first calculated the number of nucleotide differences between each human haplotype in a pairwise fashion. Then we plotted these differences as a histogram (fig. 4c). On the basis of this plot, nucleotide differences between haplotypes generate “clusters” of similar haplotypes, whereas more distant haplotypes are further apart. We found that human and chimpanzee haplotypes clustered separately. We also found that human–Neanderthal and human–Denisovan haplotypes, as well as most human–human haplotypes clustered together. Strikingly, haplogroup E haplotypes when compared with other haplotypes form a distinct cluster, suggesting that haplogroup E is unusually divergent from other human haplotypes. On the basis of linkage disequilibrium analysis, we further showed that the haplogroup E covers the entire MUC7 gene and corresponds to an unbroken haplotype block of ∼20 kb (fig. 4d).

The Divergent MUC7 Haplotype Likely Introgressed from an Archaic African Hominin Population

We focused on three possible scenarios to explain the unusual divergence of haplogroup E. First, our null hypothesis was that ancient structure within the human lineage may be responsible for this divergent haplotype in the MUC7 locus. Ancient structure refers to the remnants of genetic structure of populations ancestral to modern humans, which have been maintained across the human genome. It has been reported previously that such variation is especially detectable in regions with low recombination and high mutation rate (μ; Plagnol and Wall 2006; Lin etal. 2015). Second, we considered long-term balancing selection as a possibility to explain the maintenance of highly divergent lineages within the MUC7 region. Indeed, our previous work as well as that of others have shown evidence of such loci where ancient lineages have been maintained within human populations (Gokcumen, Zhu, etal. 2013; Leffler etal. 2013; DeGiorgio etal. 2014; Key etal. 2014; Pajic etal. 2016). However, based on the phylogenetic tree (fig. 4b) and distribution of pairwise differences (fig. 4c), both ancestral structure and balancing selection scenarios are unlikely. To be more specific, the depth of the branches within haplogroup E is relatively shallow despite the fact that they are highly divergent from other haplogroups. This indicates that haplogroup E’s origin within the human population is relatively recent as compared with other haplogroups. Corroborating this observation is that the haplotype block that defines haplogroup E remained intact (i.e., did not undergo recombination or gene conversion events, fig. 4d), also indicating that the haplogroup did not have the time to recombine within the human lineage. These observations fit neither to ancestral structure nor to balancing selection scenarios very well (Charlesworth 2006; Green etal. 2010; Yang etal. 2012).

Third, we considered introgression from an archaic hominin as a source of this haplotype, which will explain the divergence of haplogroup E from other haplogroups, as well as the apparently short coalescence of this haplogroup within the human lineage. This assumption was based on recent observations of similarly divergent haplotypes that were due to archaic introgression (Huerta-Sánchez etal. 2014; Dannemann etal. 2016). However, introgression of haplogroup E from Neanderthals or Denisovans is unlikely, because haplogroup E is exclusively found in sub-Saharan African populations, whereas introgression events between Neanderthals or Denisovans with modern humans happened after modern humans migrated out of Africa (Vernot and Akey 2014; fig. 5a, supplementary fig. S5, Supplementary Material online). As such, we hypothesize that an archaic hominin, at that time still roaming in Africa, contributed haplogroup E to the ancestors of extant Africans. The possibility of such an introgression has been discussed in recent studies (Hammer etal. 2011; Hsieh etal. 2016).

To distinguish between these scenarios and test our hypothesis in a more rigorous and quantitative way, we simulated the effects of the hypothesized ancient introgression to present day African genetic variation. It is important to note that such an introgression will have a genome-wide and not locus-specific effect. As such, before testing whether MUC7 haplogroup E was indeed the result of introgression, we first asked whether such introgression has happened and, if it did, we wanted to simulate its signatures in present day African genomes. To do this, we first assumed the demographic scenario described in Gravel etal. (2011), which is characterized by a population expansion that took place 148,000 years ago to its final present-day effective population size of 14,474 (supplementary table S5, Supplementary Material online). On the basis of these parameters, we considered two competing scenarios (fig. 5b). First, we considered a scenario with no introgression and where all the variation in observed summary statistics can be explained by drift, and variation of mutation and recombination rates in the genome. Second, we considered a scenario where some of the observed variation can be explained by invoking introgression from a “ghost” species that shares a common ancestry with the modern human lineage. To distinguish between these two competing scenarios, we used an Approximate Bayesian Computation (ABC) approach (Beaumont etal. 2002; Csilléry etal. 2012) to find the posterior probability distribution of multiple parameters (supplementary table S5, fig. S6A–F, Supplementary Material online), including the rate and duration of the putative introgression, that would explain the genetic variation observed in each of the 50 random 10 kb fragments along chromosome 4. Using this approach, we generated 300,000 simulations for each of the 50 fragments. Our results showed evidence for introgression from an unsampled archaic population to the genomes of present-day Africans (supplementary fig. S7, Supplementary Material online, Bayes Factor > 19), and as such, are concordant with recent studies describing such an introgression in Africa (Hsieh etal. 2016).

For each simulated data set we recorded the parameter values, as well as the number of introgression events that occurred during the genealogical history of the sample. Using these simulated data sets, we were able to analyze the impact of such an introgression event to the modern day genomic variation (supplementary fig. S7, Supplementary Material online). For example, Wall’s B and Watterson’s θ Estimator are both clearly higher for simulated fragments that were affected by introgression (“introgressed”) than for those that were not (“nonintrogressed”; fig. 5c). This set-up allowed us to specifically test whether MUC7’s genetic variation is most likely explained under a scenario which involves introgression or under one that does not. Indeed, in all statistics that were influenced by introgression events, the respective values for MUC7 are best explained by invoking introgression (fig. 5c, supplementary fig. S7, Supplementary Material online).

To further interrogate whether introgression best explains the data, we used S* statistics (Plagnol and Wall 2006) as applied recently to detect introgression(s) from Neanderthal and Denisovan genomes into the modern gene pool (Vernot and Akey 2014; Vernot etal. 2016). We first calculated S* for the 50 kb region that encompasses MUC7 in five replicates (each with 20 randomly chosen YRI genomes) using Northwestern European genomes as outliers (see supplementary methods, Supplementary Material online). High S* scores indicate haplotypes that harbor higher numbers of derived variants and are shown to indicate introgressed regions in the genome (Vernot and Akey 2014; Hsieh etal. 2016). To empirically contextualize our results, we first compared the S* values we got from our replicates to S* values from other 50 kb regions across chromosome 4 (fig. 5d). We found that the MUC7 region is a clear and significant outlier (P-value < 0.0001). To verify this, we simulated the S* distribution where we use the number of segregating sites and recombination rate matching what is observed for MUC7 (0.17 Cm/Mb, 495 segregating sites). In these simulations, we assume no introgression (see Materials and Methods for details). Note that the number of segregation sites is very high and the recombination rate is very low. This set of parameters is the most conservative allowing derived haplotypes without invoking introgression. Even then, our results show that MUC7 is a clear outlier and show unusually high S* scores (P-value < 0.0001), indicating the effect of introgression in this locus (fig. 5d).

It is plausible that long-term balancing selection may have maintained old lineages, mimicking the signatures of introgression. To test this possibility, we used SLiM (Messer 2013), which allows simulation of complex evolutionary scenarios, to model long-term heterozygote advantage using the posterior distribution of the models that we described above alongside with the Yoruba demographic history described in Gravel etal. (2011) (see supplementary methods, Supplementary Material online). Two new parameters included here that are different from previous simulations are the selection coefficient (s) and dominance parameter (h). The former is the strength of the selection (here we considered s = 0.001, s = 0.01, and s = 0.1) and the latter is the equilibrium frequency at which the allele is balanced (here we considered equilibrium frequencies ranging from 0.5 [h = 100] to 0.9 [h = 1.125]). Overall, we performed 1,000 forward simulations for each of the different balancing selection scenarios with different combinations of s and h. None of these scenarios explains the observed variation in MUC7 fully, where multiple summary statistics do not fit the simulated results invoking balancing selection (supplementary fig. S8, Supplementary Material online). We conclude based on these results that balancing selection in the form of heterozygote advantage with different s and h values cannot explain the observed variation as well as the introgression model.

Last, to test which of the above-described scenarios explain MUC7 genetic variation in a collective and comprehensive manner, we conducted a principal component analysis of all summary statistics for all simulation results. Briefly, we constructed a matrix of summary statistic distributions for simulated sequences that are generated under no introgression, introgression and balancing selection scenarios. On the basis of this matrix, we were able to calculate principal components of the multivariate data where principal components 1 and 2, explain 49.1% and 26.5% of the data (fig. 6). We showed that MUC7 is placed within the introgressed sequences and away from simulated sequences belonging to other scenarios. To quantify this observation, we have conducted an ABC-based categorization of MUC7, where we showed that indeed MUC7 genetic variation is best explained by the introgression scenario with robust statistical support (Bayes Factor > 36 vs. no introgression; Bayes Factor > 12 vs. balancing selection, see supplementary methods, Supplementary Material online). The simulation results show with high confidence that introgression is a better model to explain the polymorphism patterns of the MUC7 region than the alternative models which do not account for such an introgression.

Taken together, our analyses collectively support the hypothesis that ancient introgression has influenced modern African genomes and that genetic variation in MUC7 locus is influenced by this introgression event.