Evolution of the rapidly mutating human salivary agglutinin gene (DMBT1) and population subsistence strategy

Characterization of Copy Number Variation at DMBT1.

Genome-wide analysis using high-resolution tiling arrayCGH (aCGH) identified two CNVs whose location was consistent with the known polymorphic deletions described in the literature, but showing a much more extensive and complex loss and gain of copy number (17) (Fig. 1B). We used these two CNV regions as a starting point for our analysis, by interrogating the two CNVs (CNV1 involving SRCR3–SRCR6 and CNV2 involving SRCR9–SRCR11, Fig. 1B) separately. To this end, we designed several paralog ratio tests (PRTs), a form of quantitative PCR that is particularly robust in accurately calling CNVs (27). These PRTs were used to estimate copy number at each CNV in 270 individuals from HapMap phase 1. To verify our PRTs, we used concordance between PRT assays, clustering of PRT copy number estimates into distinct groups reflecting integer copy number and comparison with aCGH probe intensities (Fig. 1C and Fig. S1). By typing samples in duplicate, we estimated an upper 95% confidence limit for the error rate in determining CNV1 and CNV2 copy number to be 0.37% (366 samples, no discrepancies) and 0.33% (407 samples, no discrepancies), respectively. In addition, we further validated a subset of samples using long-PCR and fiber-FISH (Figs. S2 and S3). We show that CNV1 is a multiallelic CNV with a copy number varying between 0 and 5 per diploid genome, and the copy number variable unit includes four SRCR domains (Fig. S1). For CNV1, zero, one, and two or more copies reflect homozygous deletion, heterozygous deletion, and normal genotype of the deletion described previously (26). Sanger sequencing of homozygous deleted samples suggests that nonallelic homologous recombination (NAHR) between the 98% identical SRCR repeats carrying SRCR2 and SRCR6 is responsible for CNV1 (Fig. 1D and SI Methods). It is also clear that CNV2 is considerably more complex than the small deletion described previously, being a multiallelic CNV ranging between 1 and 11 copies per diploid genome with each repeat unit carrying a single SRCR domain.

Analysis of further samples from the CEPH-Human Genome Diversity Project (HGDP) panel of 971 individuals from 52 populations worldwide (28) showed rare individuals with a CNV2 copy number of zero. Sanger sequencing of PCR products from these individuals showed that all of the zero-copy CNV2 alleles had a breakpoint within 33 bp of sequence identical between SRCR8 and SRCR11 (Fig. S4), just upstream of the exon encoding the SRCR domain, suggesting that this allele was generated by NAHR between these repeats (SI Methods). This finding suggests that other larger CNV2 alleles have also been generated by NAHR between any of the repeats carrying SRCR domains 8–11.

DMBT1 Copy Number Variation Has a High Mutation Rate.

The extensive allelic diversity and repetitive genomic structure of DMBT1, together with the knowledge that NAHR is likely to have mediated generation of new alleles, led us to consider whether CNV1 and CNV2 had a high mutation rate. To study this directly, we used our validated PRT assays to call copy number of DMBT1 CNV1 and CNV2 on 522 samples from 40 large multigenerational families from the Centre d’Etude de Polymorphisme Humain (CEPH) collection. We robustly identified de novo copy number mutations at both loci (Fig. S5 and Datasets S1 and S2). The mutation rate at CNV1 is estimated to be 1.4% per gamete (9 out of 632 meioses, 95% CI 0.7–2.7%) and the mutation rate at CNV2 is 3.3% per gamete (21 out of 632 meioses, 95% CI 2.1–5.0%). These mutation rate estimates place both loci among the most highly mutating loci known, with comparable rates seen only for noncoding minisatellites (29) or for the mitochondrial D-loop (30). Error rates for CNV1 and CNV2 of 0.37% and 0.33% respectively are below the lower 95% CI bound of both mutation rates, showing that these high mutation rates are not due to errors in copy number calling. Furthermore, examination of the copy number calls of the individuals showing de novo mutations indicates a high posterior probability of that copy number call (Fig. S5). All mutations were of a loss or gain of one CNV repeat unit, with no evidence of a bias toward loss or gain. Analysis of our data showed that two CNV1 mutational events and one CNV2 mutational event were associated with a crossover involving flanking marker exchange at the correct position in that individual. This observation suggests that although NAHR events involving homologous chromosomes do contribute to CNV mutation rate, most events are likely to be inter- or intrachromatidal.

Global Diversity of DMBT1 and Agriculture.

To examine global diversity of DMBT1, we determined diploid copy number of CNV1 and CNV2 on the CEPH-HGDP panel (Dataset S2). We observed a similar range for CNV1 (0-5 copies per diploid genome) and for CNV2 (0–11 copies per diploid genome) as in the HapMap samples (Fig. S5). Although there was no linkage disequilibrium between copy number alleles at CNV1 and CNV2 (r² = 0 for CEU parents, r² = 0.01 for YRI parents, r² = 0.01 for all HGDP, Fig. 2A), there was a clear negative relationship between average copy number at CNV1 and CNV2 at the population level (r² = 0.11, Fig. 2B) and at the continental level (r² = 0.43). Because, across populations, the increase in CNV2 is mirrored in part by a decrease in CNV1, the total predicted number of the SRCR domains for the two copies of DMBT1 on homologous chromosomes does not mirror this trend. Nevertheless the total number of SRCR domains per diploid DMBT1 is highly variable, and the number of SRCR domains in a given DMBT1^SAG molecule is predicted to range between 7 and 20, at least (Fig. S5G).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 2.

Distribution of DMBT1 copy number values in the CEPH-HGDP panel. (A) Across individuals. Each point represents the mean unrounded PRT copy number of an individual, with the histogram on each axis representing the distribution of CNV1 copy numbers (x axis) and CNV2 copy numbers (y axis). (B) Across populations. The means of CNV1 and CNV2 in each population are plotted, colored according to continent of origin. The red dashed lines represent the value above which 99.5% of mean copy numbers of simulated populations fall. (C) Average CNV1 and CNV2 copy number in agricultural and nonagricultural populations. Populations are colored according to region (legend in B), thick line indicates median value and thin lines are 25th and 75th centiles, and P values from logistic regression, with distance from Africa as a covariate (Table 1). Agricultural population definition is according to ref. 33.

Given the fact that the SRCR domain is known to bind S. mutans, we considered that the development of agriculture and consequent increase carbohydrate-rich foods, oral S. mutans and dental caries might be a selective pressure influencing the frequency distributions of CNV1 and CNV2. To test this, we correlated CNV1 and CNV2 mean copy number for each HGDP population with an index of extent of agricultural practice for each population, as published (31). We used two statistical approaches. First, a regression analysis corrected for population effects by using distance from East Africa as a covariate. Secondly, a partial Mantel analysis corrected using a population pairwise geographical distance matrix, a population pairwise genetic distance matrix, or distance from East Africa as covariates (32). We found a negative relationship between agricultural populations and the mean CNV1 copy number and a positive relationship between agricultural populations and CNV2 (Table 1). By resampling from an empirical distribution of partial Mantel r correlations with agriculture, we estimated the genome-wide significance of this observation to be P = 0.0467, using the geographic distance matrix as a covariate, and P = 0.0410 using the genetic distance matrix as a covariate. We also tested the association between mean copy number of a population and a subsistence strategy based on carbohydrate-rich foods, such as cereals, roots and tubers, as defined previously (Fig. 2C and ref. 33). We found an association for CNV1 and a weak association for CNV2, both in the expected directions (Table 1). This finding suggests that the subsistence history of a population has affected the frequency distribution of both CNVs within DMBT1.

To investigate this association further, we called DMBT1 CNV1 and CNV2 copy number genotype using sequence read depth analysis on three published ancient DNA samples (refs. 34⇓–36, Table S1, and Fig. S5 E and F). Both Denisova and Neanderthal hominins show a high CNV1 copy number of 3 and a low CNV2 copy number of 3, within the range of hunter–gatherer human populations. Analysis of an 8,000-y-old hunter–gatherer from Loschbour in Luxembourg provides a more recent directly ancestral calibration point. He had a CNV1 genotype of 1, which is common in modern Europeans, and a CNV2 genotype of 4, which is less common but still present in modern Europeans. This observation tentatively suggests that a reduction in CNV1 copy number had occurred or was occurring, but an increase in CNV2 was yet to occur; however, further samples are required before any firm conclusions can be drawn, as the observed genotype is consistent with a copy number distribution that is unchanged from modern Europeans.

To provide further support for natural selection in shaping frequencies of CNV1 and CNV2 copy number, we used forward simulations to model CNV1 and CNV2 in situations of population expansion but no selection. Using a mutation rate derived from our pedigree analysis as well as initial allele frequency distributions based on our largest hunter-gatherer population sample (Biaka) we simulated 1,000 populations at CNV1 and CNV2 using a stepwise mutation model and a realistic demographic model (SI Methods). The resulting distribution of mean copy number for both CNVs provides an empirical test of the departure of our observed populations from this neutral stepwise mutation model. We interpreted a departure from the model as evidence of natural selection occurring on the locus on ancestors of individuals in that population. Using an estimate of mutation rate of the lower 95% confidence limit for both CNV1 and CNV2, 99.5% of simulated populations had a CNV1 mean copy number above 2.4 and 99.5% of populations had a CNV2 mean copy number above 7.4 (Fig. 2B). This observation shows that CNV2 copy numbers are lower than expected given a neutral stepwise mutation model for all populations, probably reflecting selective constraints on DMBT1^SAG protein length. For CNV1, however, six populations show mean copy number consistent with a neutral stepwise mutation model. Four of those populations are from Africa, one from East Asia and one from South America, yet all six have been classified as nonagricultural (33).

The increase in CNV2 copy number and decrease in CNV1 copy number might be due to selection for a particular phenotype more favored following the transition to agriculture. Our favored hypothesis is that a particular SRCR domain that binds S. mutans or hydroxyapatite of the tooth more weakly, thereby reducing the likelihood of caries. Analysis of the S. mutans binding domain sequence shows that there is no difference between SRCR domains in CNV1 and CNV2 (Fig. 1D). However, when manually inspecting the human GRCh37/UCSC hg19 reference sequence but also 10 HGDP samples sequenced to high-depth (ref. 34 and Table S1), the SRCR domains in CNV2 share a serine to tyrosine change which disrupts the hydroxyapatite-binding domain (37) and abolishes a strong potential mucin-type O-linked glycosylation site (Fig. 1D). Replacement of CNV1-type SRCR domains with CNV2-type SRCR domains has therefore allowed this tyrosine substitution to propagate rapidly through the DMBT^SAG molecule. This observation suggests that the transition to agriculture has been accompanied by a partial replacement of canonical SRCR domains with SRCR domains that either bind the tooth less strongly or, because glycosylation is important for binding, bind S. mutans less strongly.

Evolutionary Relationship Between DMBT1 and its Ligand in S. mutans.

If the interaction between DMBT1^SAG and S. mutans is coevolving, we might expect to see a relationship between variation of S. mutans and CNV1 and CNV2 of DMBT1 across different individuals. Adaptation of S. mutans to the DMBT1^SAG phenotype of different mouths reflects a more recent evolutionary time scale than adaptation of DMBT1 in humans. However, given that S. mutans colonizes the mouth in early childhood (38) and the doubling time of biofilm-attached S. mutans is in the order of hours (39), it is likely that it can adapt genetically to the oral environment. We genotyped 125 adult individuals resident in Leicester U.K., 92 of European origin, for CNV1 and CNV2, and sequenced part of the S. mutans gene spaP from DNA isolated from matched saliva. spaP encodes AgI/II which is the ligand for human DMBT1^SAG (40). We focused on a 1-kb region of the spaP gene encoding 336 amino acids from the C-terminal region known to contain two binding domains for human DMBT1^SAG, namely Ad1 and Ad2 (41). For 98 of our cohort (78%), only one S. mutans strain (as defined by homozygosity of the sequenced region) was found. This observation suggests very low within-mouth diversity of S. mutans, and that most people are colonized by only one strain. However, alignment of sequences showed 136 single-nucleotide polymorphisms, 82 of which altered amino acid sequence, reflecting very high levels of diversity between individuals (Fig. S6).

A difference between the allele frequency spectra (AFS) of nonsynonymous and synonymous polymorphisms can indicate selection, if synonymous polymorphisms are assumed to be neutral. Comparison of the AFS of the S. mutans spaP gene shows a difference in both the total and European-only cohort (P = 0.015 and P = 0.018 respectively), with an enrichment of polymorphisms with rare alleles (minor allele frequency < 1%, P = 0.003 for total cohort, P = 0.022 for European-only cohort, Fig. S6 E and F). This result indicates that weak negative selection is acting on spaP, consistent with previous genome-wide approaches (4).

If a particular spaP allele was adapted to a particular DMBT1^SAG phenotype, we might expect the spaP allele and the DMBT1 genotype to be associated across a number of individuals. Across our cohort, we identified 13 polymorphisms that changed amino acid sequence in spaP Ad1 or Ad2. Given the low derived allele frequencies of these polymorphisms, we had limited power to detect an association with CNV1 or CNV2 copy number. However, two of these polymorphisms affect an amino acid highly conserved across oral Streptococci (Fig. S6 G and H), and in one the derived allele (A1090D, PDB P11657) was associated with lower copy number at CNV1 (nominal P = 0.0416) and CNV2 (nominal P = 0.0169) in the European-origin cohort, which leads to an overall association with low DMBT1 copy number (P = 0.046, corrected for multiple comparisons). Only CNV2 remains associated with low copy number (nominal P = 0.0196) in the full cohort, suggesting an effect of ethnicity on this interaction.