Abstract DNA methylation is an epigenetic modification involved in regulatory processes such as cell differentiation during development, X-chromosome inactivation, genomic imprinting and susceptibility to complex disease. However, the dynamics of DNA methylation changes between humans and their closest relatives are still poorly understood. We performed a comparative analysis of CpG methylation patterns between 9 humans and 23 primate samples including all species of great apes (chimpanzee, bonobo, gorilla and orangutan) using Illumina Methylation450 bead arrays. Our analysis identified ∼800 genes with significantly altered methylation patterns among the great apes, including ∼170 genes with a methylation pattern unique to human. Some of these are known to be involved in developmental and neurological features, suggesting that epigenetic changes have been frequent during recent human and primate evolution. We identified a significant positive relationship between the rate of coding variation and alterations of methylation at the promoter level, indicative of co-occurrence between evolution of protein sequence and gene regulation. In contrast, and supporting the idea that many phenotypic differences between humans and great apes are not due to amino acid differences, our analysis also identified 184 genes that are perfectly conserved at protein level between human and chimpanzee, yet show significant epigenetic differences between these two species. We conclude that epigenetic alterations are an important force during primate evolution and have been under-explored in evolutionary comparative genomics.

Author Summary Differences in protein coding sequences between humans and their closest relatives are too small to account for their phenotypic differences. It has been hypothesized that these differences may be explained by alterations of gene regulation rather than primary genome sequence. DNA methylation is an important epigenetic modification that is involved in many biological processes, but from an evolutionary point of view this modification is still poorly understood. To this end, we performed a comparative analysis of CpG methylation patterns between humans and great apes. Using this approach, we were able to study the dynamics of DNA methylation in recent primate evolution and to identify regions showing species-specific methylation pattern among humans and great apes. We find that genes with alterations of promoter methylation tend to show increased rates of divergence in their protein sequence, and in contrast we also identify many genes with regulatory changes between human and chimpanzee that have perfectly conserved protein sequence. Our study provides the first global view of evolutionary epigenetic changes that have occurred in the genomes of all species of great apes.

Citation: Hernando-Herraez I, Prado-Martinez J, Garg P, Fernandez-Callejo M, Heyn H, Hvilsom C, et al. (2013) Dynamics of DNA Methylation in Recent Human and Great Ape Evolution. PLoS Genet 9(9): e1003763. https://doi.org/10.1371/journal.pgen.1003763 Editor: Yoav Gilad, University of Chicago, United States of America Received: March 20, 2013; Accepted: July 18, 2013; Published: September 5, 2013 Copyright: © 2013 Hernando-Herraez et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Funding: TMB is supported by the European Research Council (ERC Starting Grant, StG_20091118) and the Spanish Government (BFU2011-28549). AJS is supported by NIH grants 1R01DA033660, 1R01HG006696, and a grant from the Alzheimer's Association (2012ALZNIRG69983). IHH is supported by the European Social Fund, AGAUR (Generalitat de Catalunya, Spain) and the Spanish National Research Council (CSIC). We also thank the Spanish Government for the grant BFU2009-13409-C02-02 to AN and the Barcelona Zoo (Ajuntament de Barcelona) for an award to JPM. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. Competing interests: The authors have declared that no competing interests exist.

Introduction The genomic era is characterized by different comparative approaches to understand the effect of genomic changes upon phenotypes. In the context of human evolution, the genomes of all species of great apes have now been sequenced [1]–[4] allowing nucleotide resolution comparisons to understand the evolution of our genome. However, in contrast to these advances in comparative genomic analyses, there has been relatively little progress in the understanding of the evolution of genome regulation [5]–[9]. DNA methylation is an important epigenetic modification found in many taxa. In mammals, it is involved in numerous biological processes such as cell differentiation, X-chromosome inactivation, genomic imprinting and susceptibility to complex diseases [10]–[13]. Promoter hypermethylation is generally thought to act as a durable silencing mechanism [14]. However, the exact relationship between DNA methylation and gene expression is not clear since recent studies have also linked gene body methylation with transcriptional activity and alternative splicing [15]–[17]. At some loci DNA methylation patterns are influenced by the underlying genotype [18]–[20]. However, due to the fact that patterns of DNA methylation can change during development [16], [21], [22] or as a result of environmental factors [23], [24], the exact mechanisms governing DNA methylation states remain unclear. Most efforts to understand DNA methylation changes in primates have focused on the comparison of human with chimpanzee or macaque [6], [7], [9], [25]. This is largely attributable to the difficulty of obtaining samples from endangered species and the lack of genome sequence for the great apes. The publication last year of draft sequences of the gorilla [2] and bonobo [3] genomes facilitates a more accurate characterization of the species-specific events in all the great ape phylogeny, and interrogation of this epigenetic modification from an evolutionary point of view. Studies to date have found that DNA methylation profiles are, in general, more similar between homologous tissues than between different tissues of the same species [9]. However, differentially expressed genes between human and chimpanzee are often associated with promoter methylation differences, regardless of tissue type, establishing that some differences in the expression rates of genes between the species are associated with differences in DNA methylation. It is estimated that around 12–18% (depending on the tissue) of interspecies differences in gene expression levels could be explained by changes in promoter methylation [9]. Here we present the first comparative analysis of DNA methylation patterns between humans and all great ape species, allowing us to recapitulate the evolution of CpG methylation over the last 15 million years in these species. We used Illumina Methylation450 BeadChips to profile DNA methylation genome-wide in blood-derived DNA from a total of 9 humans and 23 wild-born individuals of different species and sub-species of chimpanzee, bonobo, gorilla and orangutan. We observed that the methylation values recapitulate the known phylogenetic relationships of the species, and we were able to characterize methylation differences that have occurred exclusively in the human lineage and among different great apes species. We also identified a significant positive relationship between the rate of coding variation and alterations of methylation at the promoter level, indicative of co-occurrence between evolution of protein sequence and gene regulation

Discussion The primary focus to date for understanding human evolution from a comparative genomic perspective has been the study of changes in DNA sequence and gene expression levels [43]–[45]. Our study of DNA methylation profiles among human and great apes adds to this wealth of information, reinforcing the view that epigenetic changes contribute significantly to species divergence, and therefore they should be considered in studies of human evolution. In this study, one of the main challenges was the technical limitation stemming from the use of arrays designed against the human genome to profile methylation patterns in great ape species with divergent genomes. We utilized a set of filters to account for these differences, and were also able to replicate the results even after limiting our analysis to those probes with 100% identity in each of the non-human reference genome assemblies. Supporting a biological role for our findings, we observed that the clustering of differential methylation within each species was highly non-random, and showed significant enrichments within functional genomic elements. From a biological perspective, it is conceivable that differences in the constitutive fractions of whole blood between species might introduce a bias due to the fact that different cell types possess distinct epigenomes [27]. This limitation is shared by nearly all comparative molecular studies of primary tissues from endangered species (i.e. great apes) due to the difficulty of obtaining relevant samples, especially in the case of wild-born individuals as the ones used in this study. However, in order to minimize this problem we removed all CpG sites that vary significantly between whole blood and the most abundant cell populations in blood. We further required a minimum threshold of 10% change in global methylation between sites in these species in order to identify differentially methylated sites, meaning that changes in the prevalence of minor cell fractions would not influence the results. Finally, while all samples were obtained from adult individuals, we could not match the ages perfectly among all samples, so in order to compensate for this effect, and to minimize the effects of intraspecific polymorphism, we focused our study on sites with low intragenus variance. Our results show that ∼9% of the CpGs we assayed showed significant methylation differences between human and chimpanzee, including the promoter regions of 745 genes (10% of those tested). We estimate that over 2,500 genes present at least some methylation changes between human and chimpanzees (≥2 differentially methylated sites separated by ≤1 kb), a higher number than that known to be affected by copy number variation or under positive selection in the same species [46]–[49]. Although the arrays we used do not provide a complete and unbiased coverage of the genome, these data suggest that epigenetic changes have been frequent during recent primate evolution and represent an important substrate for adaptive modification of genome function. Underlining this idea, the changes we observed among primates are highly enriched for sites showing intermediate DNA methylation levels. Previous studies have shown that such methylation values are often a hallmark of distal regulatory elements [34], suggesting that many epigenetic changes occurring among human and great ape species impact transcriptional regulation. Consistent with these findings, we detected global enrichments for epigenetic change within known regulatory regions, including distal regions upstream of gene transcription start sites and regions flanking CpG islands (termed ‘CpG shores’). We observed that the great ape phylogeny can be recapitulated from methylation data alone. Potential explanations for this are that methylation values could be driven by proximal DNA changes that were not controlled in this study, or that epigenetic changes have occurred independently of DNA sequence but are subject to similar rates of change either through selective pressures or neutral drift. Interestingly, we also identified a significant positive relationship between the rate of coding variation within genes and alterations of promoter methylation, suggesting a co-occurrence between changes in protein sequence and gene regulation that may be related to expression changes in fast evolving genes [50]. In contrast, and consistent with previous analysis indicating the importance of regulatory changes in evolution [51], our study also identified scores of genes that are perfectly conserved at the amino acid level between human and chimpanzee, yet showing significant epigenetic change between these two species. Furthermore, gene ontology analysis of this set showed that they are significantly enriched for the functional category of gene expression. These observations highlight the evolutionary importance of epigenetic changes that affect gene regulation, and also demonstrate that sequence-based studies are insufficient to capture the full spectrum of evolutionary change. Overall our analysis identified >800 genes with significantly altered methylation patterns specifically within each species of human and great apes, including 171 with a methylation pattern unique to humans. Analysis of these 171 genes identified interesting enrichments for a number of functional categories that could suggest a relationship to human-specific traits. For example, we observed that genes involved in the regulation of blood pressure and development of the semicircular canal of the inner ear among others, were all highly enriched for DNA methylation changes specifically in the human lineage. While major changes in circulatory physiology are required for upright locomotion, the inner ear provides sensory input crucial for maintaining balance. Furthermore, a previous study of primates and other mammals has shown that the size of the semicircular canals is correlated with locomotion and with relatively larger canals found in species that utilize fast or agile movement [52]. While these trends hint at the potential importance of epigenetic changes in the evolution of several human-specific features, we caution that at this stage they should be considered as preliminary, as our studies were performed using DNA derived from whole blood, and it is well known that epigenetic patterns often vary widely between different tissues of an organism [27]. Therefore further studies in physiologically relevant tissues will be required to confirm the significance of these findings. However, we note that comparison with previously published data [9] suggests that many of the changes in DNA methylation that we detected between blood of human and chimpanzee appear to be conserved across several other tissues, suggesting that inter-specific differences observed in blood can in some cases be informative for other tissues. Although sequencing studies have undoubtedly provided major advances in our understanding of primate evolution, our analysis of primate epigenomes unveils many novel differences among the great apes that are not apparent from purely sequence-based approaches. Of particular note is the fact that we identify enrichments in multiple independent functional gene categories which suggests that regulatory changes may have played a key role in the acquisition of human-specific trait. Therefore, epigenetic alterations likely represent an important facet of evolutionary change in primate genomes. Future studies that integrate epigenetic data with recent detailed maps of functional elements, selective constraint and chromatin interactions in the human genome [53]–[55] will likely provide many novel insights into genomic and phenotypic evolution.

Methods Ethics statement The non-human research has been approved by the ethical committee of the European Research Union. No living animal has been used and DNA has been obtained during standard veterinary checks. Methylation profiling of human subjects was approved by the Institutional Review Board of the Icahn School of Medicine at Mount Sinai (HS#: 12-00567 HG). Hybridization and normalization We obtained methylation data from peripheral blood DNA extracted from 9 adult humans, 5 chimpanzees, 6 bonobos, 6 gorillas and 6 orangutans. All individuals were unrelated adults and the non-human primates were all wild born. DNA samples were bisulfite converted, whole-genome amplified, enzymatically fragmented, and hybridized to the Infinium HumanMethylation450 BeadChip which provides quantitative estimates of methylation levels at 482,421 CpG sites distributed genome-wide. The assay was performed according to the manufacturer's instructions. The BeadChip array data discussed in this publication have been deposited in NCBI's Gene Expression Omnibus and are accessible through GEO Series accession number GSE41782. Due to the low density of probes targeting non-CpG dinucleotides (<0.7%) on the array, we focused our study on CpG sites. Since the 50 bp probes on the array were designed against the human reference genome but we performed hybridizations utilizing DNA from different great ape species, we first mapped the probe sequences to the chimpanzee (panTro3), bonobo (panPan1), gorilla (gorGor3) and orangutan (ponAbe2) reference genomes using BWA [56], allowing a maximum edit distance of 3. We then assessed probe performance as a function of the number and relative location of sequence differences at the probe binding site in each primate genome (Figure S1 and Figure S2). Based on this analysis, in each species we only retained those probes that had either a perfect match, or had 1 or 2 mismatches in the first 45 bp but no mismatches in the 3′ 5 bp closest to the CpG site being assayed. We also removed all probes that contained human SNPs with minor allele frequency ≥0.05 within the last 5 bp of their binding site closest to the CpG being assayed [57]. Using published SNP data [26] for each species we removed probes containing SNPs with minor allele frequency ≥0.15 within the last 5 bp of their binding site closest to the CpG being assayed. We also removed all probes that contained more than two SNPs with minor allele frequency ≥0.15 in the first 45 bp. Methylation values for CpG sites in each sample were obtained as β-values, calculated as the ratio of the methylated signal intensity to the sum of both methylated and unmethylated signals after background subtraction (β-values range from 0 to 1, corresponding to completely unmethylated and fully methylated sites, respectively). Within each individual, probes with a detection p>0.01 were excluded. We performed a two color channel signal adjustment and quantile normalization on the pooled signals from both channels and recalculation of average β-values as implemented in “lumi” package of R [58]. The Illumina Infinium HumanMethylation450 BeadChip contains two assay types (Infinium type I and type II probes) which utilize different probe designs. As the data produced by these two assay types shows distinct profiles (Figure S6), to correct this problem we performed a BMIQ (beta mixture quantile method) [59] on the quantile normalized data sets. Using a published human data set [27] we identified differentially methylated sites between whole blood and CD4+ T-cells, and between whole blood and CD16+ neutrophils, representing the two most abundant cell fractions of blood (comprising ∼13% and ∼65%, respectively) (Wilcoxon rank-sum test, FDR-adjusted p<0.05 and mean β-value difference in each case ≥0.1). These sites (n = 10,151) were removed to mitigate potential confounders due to differing proportions of blood cell types among primates, leaving for comparison only those sites that do not significantly vary among the most abundant cell types of blood. β-values can be interpreted as the percentage of methylation at a given site. A β-value of 0.1 indicates that there has been a change in methylation in 10% of the molecules tested. Because our analyses required a mean β-value difference >0.1 to achieve significance, this threshold means that changes in blood cell fractions representing <10% of whole blood will be unlikely to affect our results. The final dataset after all filtering steps comprised 114,739 probes shared across all great ape species, and 291,553 probes shared between human and chimpanzee. Phylogenetic relationships To investigate the global correspondence of DNA sequence differences between species and the degree of methylation changes, we examined the Enredo-Pecan-Orthus (EPO) whole-genome multiple alignments of human, chimpanzee, gorilla, and orangutan [Ensemble Compara.6_primates_EPO] [28], [29]. Considering only those blocks with alignments for all great apes, we first excluded regions containing gaps or indels and then calculated pairwise distances between these four species based on the frequency of single nucleotide substitutions. To calculate the global changes in methylation we used a distance matrix, we first averaged the β-values per probe within a species and then calculated the difference between two species using Euclidean distances. We built phylogenetic trees based on the methylation states of 114,739 filtered probes (perfect match probes and probes containing 1 or 2 mismatches in the first 45 bp) (Figure S3). We used the “ape” R package to construct the phylogenetic tree using the Neighbor-Joining algorithm and 1,000 bootstraps of the resulting tree [60]. We repeated the analysis using only the subset of probes with a perfect match to each of the primate reference genome assemblies (n = 31,853) (Figure 1B). Differentially methylated sites To identify only those methylation differences that represent fixed changes between genera, we retained only those CpGs with low methylation variance within each genus (intragenus standard deviation <0.1). This filtering step resulted in the removal of 1,377 CpGs in human, 5,224 in the Pan genus, 5,289 in Gorilla and 5,740 in Pongo, with the resulting final set being 99,919 CpGs shared across all five species. We performed six pairwise comparisons among groups (Human-Pan species/Human-Gorilla species/Human-Pongo species/Pan species-Gorilla species/Pan species-Pongo species/Gorilla species-Pongo species). We defined a site to be genus-specific differentially methylated if all three comparisons with other groups were significant (Wilcoxon rank-sum test, FDR-adjusted p<0.05) and mean β-value difference in each case ≥0.1. We also tried other statistical approaches (linear modeling, limma package, [59]) and obtained very similar results (concordance for 98% of the sites). All coordinates quoted are based on hg19. We intersected human probe coordinates provided by Illumina with RefSeq genes, retaining CpG sites overlapping genes (−1500 bp from TSS to 3′UTR). We defined a gene to be differentially methylated if there were at least two differentially methylated CpG sites separated by ≤1 kb. To assess significance of these observations we performed a permutation test, as follows. Based on the number of differentially methylated sites detected in each species (Human = 2,284; Pan = 1,245; Gorilla = 1,374; Orangutan = 5,501) we randomly sampled from the 99,919 CpGs and then determined the number of clusters (at least two differentially methylated CpG sites separated by ≤1 kb), repeating this process 10,000 times to create the null distribution. The p-value corresponded to the number of times that differentially methylated clusters appeared within the null distribution divided by the number of permutations (n = 10,000). The Genomic Regions Enrichment of Annotations Tool (GREAT version 2.0.1) [31] was utilized to identify significant enrichments (FDR-corrected p<0.05) for Gene Ontology biological processes. While tools for identifying enriched GO terms are usually based on genes, GREAT permits the assignment of biological function to non-coding genomic regions by analyzing the annotations of nearby genes. For this analysis regulatory regions were associated to the single nearest gene situated within 10 kb. The background data set was the 99,919 CpG sites interrogated in all great ape species. In order to evaluate the positional context of the differentially methylated sites, we compared the distribution of these 10,404 sites detected among the primate species with all 99,919 CpGs tested. Permutation p-values were calculated as described above using 10,000 iterations. X-chromosome inactivation We performed two color channel signal adjustment and quantile normalization on males and females separately. Due to the different methylation pattern in females no BMIQ normalization was done in this data set. For studies of DNA methylation on the X-chromosome that might be linked with XCI between species, we searched for CpG sites presenting no significant changes between males and females in a specific lineage (mean β-value difference <0.1) but showing significant changes in all the other species (mean β-value difference between sexes >0.1). Human-chimpanzee analysis The number of probes shared between human and chimpanzee after applying our mapping and SNP filters was 291,554. Based on this set of probes, we performed a separate two color channel signal adjustment and quantile normalization of the raw data using only human and chimpanzee samples. We performed a BMIQ normalization to correct the probe design bias. After excluding probes with a standard deviation within either species >0.1 we retained a total of 289,007 probes. Differentially methylated sites were those with p<0.05 (Wilcoxon rank-sum test, FDR-adjusted p<0.05) and a mean β-value difference ≥0.1. From the total set of 13,454 human:chimpanzee orthologous genes [1], we removed genes with <150 or >1500 amino acids, and then compared the number of amino acid changes and the K A /K I ratio of genes with robust alterations of promoter methylation (mean β-value difference of top 2 probes within promoter ≥0.1, considering CpGs located ≤1,500 bp upstream of Refseq gene TSSs, in the 5′UTR or the 1st exon, n = 745) versus those without methylation changes (n = 6,507). The Gene Ontology enRIchment anaLysis and visuaLizAtion tool (GOrilla) [41], [42] was utilized to obtain the functional enrichments within the 184 genes conserved at amino acid level, yet having significant epigenetic differences at their promoter. The data set containing 7,252 human∶chimpanzee 1∶1 orthologs was used as a background.

Acknowledgments We wish to thank G. Santpere for constructive comments and useful discussions.

Author Contributions Conceived and designed the experiments: TMB AJS IHH. Performed the experiments: IHH JPM. Analyzed the data: IHH ME HH TMB AJS JPM. Contributed reagents/materials/analysis tools: MFC PG CH AN. Wrote the paper: IHH TMB AJS.