Evidence for APOBEC3G activity in human evolution

We set out to investigate whether A3 proteins have led to enzyme-induced evolutionary changes via introduction of mutation clusters in the germline. For this purpose, we used an approach similar to that established for the detection of A3 activity in human cancer (Roberts et al. 2013). Briefly, the methodology exploits the fact that different A3 paralogs have a unique mutagenesis pattern that consists of preferred substrate nucleotide motifs and a tendency to mutate several cytosines on the same ssDNA strand (Supplemental Fig. S1). We first conducted a multispecies comparative genomic analysis that identified mutations that have occurred along the human lineage since its divergence from chimpanzee, i.e., positions in the reference human genome that differ from their ancestral state (Table 1; Fig. 1A; Methods). Second, we considered clusters of such mutations as groups of at least two mutations within ≤50 base pairs (bp) from each other (Fig. 1B,C). We excluded clusters with potential complex mutations that are ≤10 bp from each other and clusters composed of mutations arising from different DNA strands or different nucleotides in the ancestral state (Fig. 1C). To further minimize the extent to which nearby yet independently occurring mutations are considered as a cluster, we further excluded clusters according to a cluster specific P-value that accounts for the probability of clustering independent events (Roberts et al. 2012). We refer to these clusters by their ancestral nucleotide, e.g., a cluster with all mutations originating from cytosine is termed a C-coordinated cluster (Fig. 1B).

Figure 1.

Clustered mutations in the human lineage. (A) Phylogenetic tree of the species used for computing the human-chimpanzee ancestral state (dashed arrow). (B) Multiple sequence alignment demonstrating an incidence of a C-coordinated cluster in which the ancestral state harbors closely located cytosine residues that were mutated in the human lineage. (C) An example of an N-coordinated cluster where the mismatches originated from different ancestral nucleotides. Alignment matches are represented by dots; omitted parts of the sequence are indicated by tildes (∼). Red letters and arrows indicate mismatches in the cluster. Positions in which an outgroup does not agree with the other outgroups are marked in blue, and positions in which the ancestral state could not be inferred are marked in gray. The boundaries of clusters are highlighted in gray.

Open in new tabDownload PowerPointLink to figure

Table 1.

Mutations along hominid lineages

Open in new tabDownload PowerPointLink to table

Next, we compared the prevalence of C- (or G-) coordinated clusters to that of A- (or T-) coordinated clusters, with C- (or G-) but not A- (or T-) coordinated clusters being potentially a result of A3 activity (on either strand). We observed a significantly higher abundance of C- (or G-) coordinated clusters than A- (or T-) coordinated clusters (243,923 against 210,172—16% higher) (Supplemental Fig. S2) and a higher abundance of C- (or G-) coordinated clusters than expected by chance (Methods; Supplemental Table S1). In order to assess whether a subset of the former clusters were formed by members of the AID/APOBEC family of enzymes, we searched for the A3- or AID-associated mutagenic patterns within the detected C- or G-coordinated clusters while considering the different sequence motif they act on. Namely, we tested the clustered mutations for enrichment in the sequence motifs of the following enzymes: A3B (TC [Hultquist et al. 2011] and TCW [Roberts et al. 2012] motifs; the mutated nucleotide is underlined; W – A or T), A3F (TC [Liddament et al. 2004] and TTC [Liddament et al. 2004; Armitage et al. 2014]), A3G (CC [Hultquist et al. 2011] and CCC [Chelico et al. 2006]) and AID (WRC [Shen and Storb 2004]; R – A or G) by calculating the ratio between the number of mutations falling within the motif to all clustered mutations (e.g., the number of all mutations in the underlined base in CCC divided by the overall number of mutations in C-clusters). Motif enrichment was then determined based on this ratio, following normalization by several similarly calculated backgrounds (Methods). When using the genomic background, we found mutations in the A3G-associated CCC motif to be highly enriched within C- or G-coordinated clusters (Fig. 2A). Additionally, the less stringent of the two A3G motifs, CC, was also found to be enriched, albeit to a lesser extent (Fig. 2A). These findings were robust to the method of normalization or clustering and when considering local variations in nucleotide context (Methods; Supplemental Figs. S3–S6). In contrast, no other APOBEC-related motifs were found to be enriched across clusters by any of the different approaches (Fig. 2A). We further tested all 48 possible trinucleotides that harbor mutated cytosines (i.e., nucleotide triplets with a cytosine in the ancestral state showing a derived mutation in any of the three positions of the triplet) and found, with the exception of CpG-containing triplets, CCC to be the only significantly enriched motif (Fig. 2B; Supplemental Fig. S7).

Figure 2.

Identification of APOBEC3G mutagenesis pattern in human evolution. (A) A3G CC and CCC motifs are significantly enriched within the set of C- (or G-) coordinated clusters with cluster P-value ≤0.0001 (q < 0.006, one-tailed Fisher's exact test after Bonferroni correction. P-values were calculated by comparing the frequency of mutations within and outside a given motif while controlling for the frequency of C or G nucleotides within and outside this motif in the genome). Clusters are not enriched with any other APOBEC-related motif. (B) CCC, the A3G stringent motif, is the only CpG-free trinucleotide motif that is significantly enriched in C- (or G-) clusters (Bonferroni-corrected q-values were calculated as in A). (C) As expected from the known mechanism of A3G, mutations found in C- (or G-) coordinated clusters demonstrated higher enrichment relative to nonclustered events (one-tailed Fisher's exact test P = 1.35 × 10⁻⁸). Clustered mutations found in transcribed regions, which have a higher tendency to form ssDNA, exhibit higher enrichment levels compared with all clustered C (or G) mutations (one-tailed Fisher's exact test P = 0.045). Fold-change value between different bars is indicated above them. (D) A strong positive correlation (Pearson's r = 0.93, P < 0.05) is observed between the level of CCC motif enrichment and the significance threshold used to filter cluster sets (−log[P-value]).

Open in new tabDownload PowerPointLink to figure

While a higher mutation rate is expected at CpG dinucleotides due to the spontaneous deamination of 5-methylcytosine (Bird 1980; Sved and Bird 1990), several lines of evidence suggest this effect cannot account for the enrichment of the A3G-associated motif. First, when repeating our analysis with a CpG-filtered mutation set and a CpG-masked genome as a background, we again found an enrichment of the CCC motif in C- and G-coordinated clusters (Supplemental Fig. S8A). Second, the high enrichment value of the A3G motif in the full data set (Fig. 2A) is not restricted to occurrences of a CCCG motif (Supplemental Fig. S8B). Third, we found the CCC motif to contribute additively to mutated cytosines in CCCpG (Supplemental Fig. S8C). In addition, we tested for the enrichment of TTT motifs within T- (or A-) clusters as a negative control for the possibility that CCC enrichment is biased due to any possible mutational effects in homotypic trinucleotide motifs and found no evidence of enrichment (Supplemental Fig. S9).

Since A3G mutagenesis mainly occurs in clusters and is limited to ssDNA, one would expect that mutations within genomic regions with a higher tendency to form a single-stranded intermediate would show increased enrichment of the A3G-associated motif. Transcription is a major physiological cause for the unwinding of genomic double-helical DNA and therefore forms a convenient substrate for A3G (Schumacher et al. 2005). As expected from the biochemical activity of A3G, clustered mutations were found to exhibit elevated motif enrichment in comparison to single point-mutations, while clusters found within transcribed regions were also found to be further enriched with the A3G motif (Fig. 2C).

The strong correlation of CCC motif enrichment and the confidence level of the clusters (r = 0.93) (Fig. 2D) demonstrate tight linkage between an overrepresentation of the A3G motif and the occurrence of clustered mutations. Namely, as the probability that a given cluster is a random group of separated point mutations decreases, the level of enrichment of the CCC motif increases, indicating an A3G-specific mutagenesis pattern. Combined, this body of evidence strongly supports that concurrent A3G-induced mutations have contributed to human evolution.

APOBEC3G-induced genome evolution is common to all hominids

Our findings indicate that APOBEC3G has contributed to human genome evolution by the introduction of clustered mutations. We next examined whether this phenomenon is limited to the modern human lineage, or—as could be expected by the timing of the APOBEC family expansion—may be widespread across hominids. We applied our approach to test for APOBECs activity in chimpanzee and archaic humans (Neandertal and Denisovan), using recently published sequencing data (Meyer et al. 2012; Prüfer et al. 2014). This allowed us to distinguish between mutations specific to chimpanzee, those arising early in the branch common to both modern and archaic humans after the human-chimpanzee split, and more recent mutations specific to Neandertals, Denisovans, or modern humans (Table 1; Fig. 3). We found that all lineages exhibit an enrichment in the A3G signature within clustered mutations (compared to the genomic background) (Fig. 3) and a tight positive correlation between CCC motif enrichment and the cluster's level of confidence given by their P-value (mean r = 0.96) (Supplemental Fig. S10). Combined, these results illustrate that A3G-related mutagenesis is a common feature across hominids.

Figure 3.

APOBEC3G clustered mutations are observed across all hominid lineages. An analysis identical to that performed on mutations in the human lineage was used to analyze chimpanzee mutations (from the point of divergence from the human lineage) and lineage-specific mutation sets from modern human and the archaic hominins: Denisovan and Neandertal. All sets contain clusters enriched with mutagenic patterns associated with A3G activity (cluster P < 0.01; one-tailed Fisher's exact test P < 6.8 × 10⁻¹¹). Additionally, mutations common to all three Homo lineages (marked in gray on the phylogenetic tree) also showed evidence of enrichment with a value of 1.40 ± 0.03. The lengths of branches in the phylogenetic tree are not drawn to scale. Error bars represent the standard errors that were calculated using a block bootstrap approach.

Open in new tabDownload PowerPointLink to figure

Given that A3G is a primate-specific enzyme, we sought to confirm our results by using a set of nonprimate genomes as a negative control (Keane et al. 2011). A similar analysis performed on a genomic set of mutations called between different mouse strains yielded no enrichment in the CCC motif (Supplemental Fig. S11). Moreover, the mouse A3 motif (TYC; Y – C or T) was also not found to be enriched within clustered mutations (Supplemental Fig. S11). These findings point to APOBEC-induced evolution being a primate-specific phenomenon. This nonprimate control also supports the argument that the effect observed is not the result of nonenzymatic processes, such as spontaneous cytosine deamination or guanine oxidation, which are common to both primates and the species considered for the nonprimate control.

APOBEC3G activity is overrepresented in functional regions

We defined the set of potentially A3G-induced mutations as all the mutations within significant C- or G- coordinated clusters (P < 0.01) that have at least one CCC mutation (Table 1; Supplemental Tables S2, S3). We found a significant overrepresentation of such potential A3G clusters in all hominid lineages (Supplemental Fig. S12). In order to test for characteristics unique to the set of A3G-induced mutations, we compared the genomic distribution of this set of mutations to a control set composed of all other mutations, independently for each of the lineages. As shown in Figure 2C, clustered mutations within transcribed regions were highly enriched in the A3G motif. This evidence and the fact that transcribed regions have an ssDNA intermediate make it reasonable to expect that actively transcribed regions will be more susceptible to A3G-induced mutagenesis. Indeed, the A3G mutations set includes 1.7-fold more mutations in transcribed regions relative to the control set (Supplemental Fig. S13). Furthermore, a robust correlation was found between transcription levels to the normalized abundance of A3G-induced mutations (Fig. 4A). As genes are the main transcriptional units, we investigated whether our A3G set is enriched with intra-genic mutations. Indeed, exonic (but not intronic or inter-genic) regions were found to be enriched in most of the hominin lineages tested (Fig. 4B). When distinguishing between different exon types, enrichment in coding sequences (CDS) and 5′ untranslated regions (5′ UTR) was detected in several lineages (Supplemental Fig. S14). We next tested whether the enrichment in A3G mutations within CDS has any evolutionary impact in terms of potential altered functionality. In fact, more than a third (51/147) of the A3G-induced mutations in CDS were found to cause amino acid substitutions across 31 genes in the Homo lineage alone (Supplemental Table S4). In some cases, a cluster is only partially localized in a CDS, but in other instances the whole A3G-induced cluster, constituting up to five different amino acid substitutions, can be found in the coding sequences of a single exon within the genes: ERICH3, CEMP1, MADCAM1, PLIN4, and IBSP.

Figure 4.

APOBEC3G mutations are associated with transcription and overrepresented in exons and regulatory regions. (A) The fold enrichment of the CCC motif in the first nine deciles (solid circles) of the A3G set relative to other mutations is strongly and positively correlated with expression levels (Pearson's correlation coefficient r = 0.87, P = 0.0012). Notably, the top decile (open circle) behaves differently in terms of fold-change and mutation counts. This can be expected from the increased evolutionary conservation of the most highly transcribed regions, which tend to show a lower tolerance to mutations, let alone clusters of mutations. (B) We distinguished between inter-genic regions, introns, and exons. An overrepresentation of A3G mutations within exons relative to the proportion of non-A3G mutations was found in most lineages (q < 0.05 for all regions that were tested in Homo and Pan and both, exons and inter-genic regions in modern and archaic humans. Two-tailed Fisher's exact test with Bonferroni correction). Error bars represent the standard errors that were calculated using a block bootstrap approach. (C) Regulatory regions are highly enriched in A3G mutations. Several transcriptional regulatory regions were inspected: open chromatin, transcription factor binding regions and binding sites, enhancers, and DNase I hypersensitive sites. All regions were found to be enriched in most of the lineages (q ≤ 0.05, two-tailed Fisher's exact test after Bonferroni correction; N indicates nonsignificance). Error bars represent the standard errors that were calculated using a block bootstrap approach. (D) The binding regions of various transcription factors were tested, and several of them showed differential A3G-related mutagenic activity before and after the split of archaic and modern humans. Red bars indicate enrichment before the split, while green bars mark enrichment in the derived lineages after the split. Only results with P < 0.05 (two-tailed Fisher's exact test) are shown. (*) FDR-adjusted P-value < 0.05. Error bars represent two-sided 95% confidence intervals for the exact test.

Open in new tabDownload PowerPointLink to figure

The overrepresentation of A3G mutations in 5′ UTRs led us to hypothesize that A3G mutagenesis could also have functional effects in regulatory regions. Indeed, we found all regulatory elements such as enhancers, transcription factor binding sites, and DNase I hypersensitive regions (DHS) to be enriched for A3G-induced mutations (Fig. 4C). The additive effect between overlapping CCC and CG motifs (Supplemental Fig. S8C) suggests that the A3G-mutations set contains mutations formed by both A3G-induced mutagenesis and nonenzymatic CpG deamination. To rule out that enrichment in functional genomic regions is driven by the latter, we excluded CCCG mutations from the set and obtained comparable results (Supplemental Fig. S15). Next, we focused on transcription factor binding regions (ENCODE ChIP-seq peaks). Interestingly, when comparing ancestral and modern human branches (i.e., before and after the split of modern and archaic humans), several transcription factor binding regions were found to exhibit different tendencies to accumulate A3G mutations between ancestral and derived branches (Fig. 4D). Namely, we found that some of the binding regions were significantly enriched in A3G mutagenesis in early human evolution while others showed a significant increase in A3G activity only in modern humans. Notably, this provides evidence of a shifting landscape of A3G activity along different regions of the genome at different evolutionary time periods and suggests that A3G may have played a role in the evolution of functional regulatory differences among hominids (Fig. 4D).

In conclusion, after carefully controlling for possible confounders, C- (and G-) coordinated clusters across hominids show an enrichment of the hallmark signature of A3G mutagenesis, suggesting that its activity may have been responsible for the introduction of tens of thousands of mutations with unique features in hominid evolution and probably many more in the entire primate tree. Specifically, these mutations are localized in transcribed regions, especially exons, and in regulatory elements, where they show signs of varying rates of occurrence at different evolutionary time periods. Genomic changes in both regions individually could have contributed to the evolution of phenotypes of both modern humans and closely related species.

Calling mutations

Lineage-specific substitutions in human were computed by parsimony-based approach using alignments from the UCSC Genome Browser for the human reference (hg19) and three outgroups: chimpanzee (panTro2), orangutan (ponAbe2), and rhesus macaque (rheMac2) (Kent et al. 2002). Sites where the chimpanzee allele was available (see “filtering” below) and its state was confirmed by matching that in either orangutan or rhesus macaque were kept for further analysis. In these cases, the chimpanzee allele was assigned as the ancestral human state, and human divergences were then called for positions where the human reference was different from the inferred ancestral state.

For filtering, we labeled as “missing” all positions for each reference genome individually, falling in an alignment gap, with an ambiguous allele (e.g., “N”), or falling within a region of poor synteny with human (as determined by the UCSC human/chimpanzee, human/orangutan, or human/rhesus macaque syntenic nets). Positions in the alignment where human or chimpanzee, or where both orangutan and rhesus macaque, were missing were filtered out from our data set. Mutations within simple and low-complexity repeats and segmental duplications were filtered out due to their low reliability (Jurka 2000). For all human/chimpanzee syntenic regions, the Altai Neandertal and Denisovan states were determined using the sequencing data of Prüfer et al. (2014) and Meyer et al. (2012), respectively, and filtered as described in Prüfer et al. (2014). Using the inferred human/chimpanzee ancestral state, we distinguished between mutations occurring before and after the split of modern and archaic human populations. Genomic variation data for different mouse strains was obtained from Keane et al. (2011). Mutations specific to the C57BL/6J mouse strain (mm9 reference genome) were called as previously described using the NZO/HlLtJ, NOD/ShiLtJ, and A/J strains as outgroups. Mutation calls performed using these strains resulted in a set of mutations which is comparable in size to that of the tested hominid branches (∼2.6M mutations) (Table 1). Additionally, these strains, while sufficiently distant on an evolutionary scale to provide similar numbers of mutations as hominids, are sufficiently close as to ensure a similarly low likelihood for multiple mutations at a site, making them effectively comparable under the parsimony-based approach utilized in our pipeline.

Clustering

Mutations were clustered based on previously published methodology (Roberts et al. 2012, 2013) with a single modification where clusters were defined as a group of at least two mutations that are separated by 50 bp or less from each other. Clusters were named by the ancestral nucleotide state on which the mutation arose (e.g., clusters composed solely from mutation derived from cytosines were classified as C-coordinated clusters) (see Supplemental Data). Clusters composed of mutations originating from different ancestral nucleotide states were classified as N-clusters. The interval between pairs of mutations was limited to 50 bp in order to allow a sufficient recovery of non-N clusters, when taking into account the large number of point mutations generated along the evolution of primate species and the high probability of a random recent substitution to disturb a previously formed long coordinated cluster. The same clustering methodology was used in Supplemental Figure S6, except the maximal distance between mutations were set to either 100, 300, or 1000 bp. Cluster P-values were computed using the negative binomial distribution:

Motif enrichment

We tested for an enrichment in APOBEC mutagenesis patterns within mutation clusters. Both positive and negative strands were tested for each motif using its reverse complement sequence strands (e.g., to detect CCC motif in the negative strand, we tested for the GGG motif within G-coordinated clusters). Motif-associated signals were estimated as the fraction of clustered mutations falling in a certain motif out of the total number of clustered mutations. For example, for the CCC motif we estimate the signal as

Genomic background: The original signal compared against the ratio between the number of occurrences of a motif in the genome to the occurrences of the ancestral nucleotide in the genome, e.g., ${\displaystyle \frac{{\text{Genome}}_{CCC}}{{\text{Genome}}_C}}$. Therefore the enrichment is defined as:

Random sets of mutations: For each set of mutations, we generated 100 random sets of mutations as an arbitrary control. Each set was of the same size and with the same mutation types as the original set. Namely, for each mutation in the original set we randomized its genomic location using BEDTools shuffle command (Quinlan and Hall 2010). We then compared the original signal against the fraction of mutations falling within motifs out of the total number of clustered mutations for all 100 sets, e.g., ${\displaystyle \frac{{\text{Randomized clustered mutations}}_{CC\underset{\_}{C}}}{{\text{Randomized clustered mutations}}_{\underset{\_}{C}}}}$. Therefore, the enrichment is defined as:

Random sets of clusters: For each set of clusters, we generated 100 sets of random clusters. Each set included the same size of clusters from each type, with the same cluster length but in random genomic loci. In addition, the random clusters preserve the same number of mutations, with the same mutation types and the inner mutations' distances from clusters boundaries. Then, we calculated a ratio similar to the one that was computed for random mutations sets.

Local nucleotide context: We calculated the ratio between the number of occurrences of a motif in the 10-kb region centered on the middle of each cluster to the occurrences of the ancestral nucleotide in the same context, e.g., ${\displaystyle \frac{{\text{Cluster context}}_{CCC}}{{\text{Cluster context}}_C}}$. Therefore the enrichment is defined as:

Calculating the expected number of clusters

The expected number of C- (or G-) coordinated clusters was calculated as the probability that all mutations in a given cluster originated from the same ancestral nucleotide using the prevalence of mutations derived from each ancestral state: $\text{Expected}=\sum _{i=2}^k\left(N_i\cdot p^i\right)$, where k is the maximal number of mutations in a single cluster, N_i is the number of clusters with i mutations, and p is the prevalence of mutation from a specific ancestral state.

Definition of A3G clusters

A3G clusters were defined as C- (or G-) coordinated clusters with at least one CCC (or GGG) mutation. The expected number of A3G clusters was calculated by the probability of having at least one CCC mutation in a cluster:

Error estimation

Standard errors were calculated using a block bootstrap approach (Kunsch 1989; Liu and Singh 1992; Lahiri 2003; Keinan et al. 2007). Resampling was carried out by dividing the genome into 35,165 equal-sized blocks of ∼90 kb. The blocks were then randomly sampled with replacement to create 10,000 bootstrap genome samples.

Functional regions

A3G-induced sets of mutations were defined as all mutations found within A3G clusters. The fold-change for mutations in different genomic regions between A3G-induced mutations to all other mutations (control set) was computed as the ratio between the fraction of mutations falling in a given region relative to those in the control set. The enrichment in a specific functional region was calculated as follows:

Transcribed regions and open chromatin data from H1-hESC, DHS peaks, and chromatin immunoprecipitation (ChIP)-seq peaks from ESC were used as transcription factor binding regions and were obtained from the ENCODE Project (Rosenbloom et al. 2013). The respective ENCODE tables—H1-hESC (under transcription track, for transcription levels data), wgEncodeRegTfbsClusteredV3 (under Txn Factor ChIP track, for TF binding regions data) and H1-hESC Syn Pk (under Open Chrom Synth track for open chromatin data)—were downloaded from the UCSC Table Browser (Karolchik et al. 2004). Transcription factor binding site data were obtained from factorbook version 3 (Wang et al. 2013). CpG islands data were retrieved from the UCSC Table Browser (Karolchik et al. 2004). Data for the location of transcriptional enhancers were downloaded from the ORegAnno database, picking entries annotated as regulatory regions (Griffith et al. 2007). To ensure that the enrichment in genes is not due to A3G activity in the retrotransposition phase of processed pseudogenes, we computed the number of Homo lineage A3G mutations within processed pseudogenes and found negligible contribution (83 out of 13,011 mutations). Processed pseudogenes data were downloaded from the UCSC Table Browser (Karolchik et al. 2004).