Report

Reconstructing the DNA Methylation Maps of the Neandertal and the Denisovan

See allHide authors and affiliations

Science  02 May 2014:
Vol. 344, Issue 6183, pp. 523-527
DOI: 10.1126/science.1250368

Methylating the Family Tree

DNA sequences show a high level of similarities between humans and ancient hominids but the degree to which there are differences between methylated regions in their genomes that may explain phenotypic differences is unclear. Gokhman et al. (p. 523, published online 17 April) demonstrate that naturally degraded methylated cytosines in ancient DNA are converted to thymines and can be used to reconstruct ancient methylomes. The results suggest differences in methylation in bone tissues between modern humans and ancient hominids in a set of genes important for limb development.

Abstract

Ancient DNA sequencing has recently provided high-coverage archaic human genomes. However, the evolution of epigenetic regulation along the human lineage remains largely unexplored. We reconstructed the full DNA methylation maps of the Neandertal and the Denisovan by harnessing the natural degradation processes of methylated and unmethylated cytosines. Comparing these ancient methylation maps to those of present-day humans, we identified ~2000 differentially methylated regions (DMRs). Particularly, we found substantial methylation changes in the HOXD cluster that may explain anatomical differences between archaic and present-day humans. Additionally, we found that DMRs are significantly more likely to be associated with diseases. This study provides insight into the epigenetic landscape of our closest evolutionary relatives and opens a window to explore the epigenomes of extinct species.

The high-coverage genomes of the Denisovan (1) and the Neandertal (2), from whom present-day humans split 550,000 to 765,000 years ago (2), allowed the comparison of our DNA to that of our closest extinct relatives. However, differences could be sought exclusively at the sequence level and consequently did not provide direct insights regarding epigenetic regulation, which likely underlies many of the traits that distinguish these human groups.

DNA methylation is a key hallmark of gene activity in mammals, where it occurs almost exclusively in cytosines in the context of CpG dinucleotides. In promoters, it is inversely correlated with expression level, and in intergenic regions it is linked to the activity of enhancers and insulators and to silencing transposable elements (3, 4).

Classic methods to measure full methylomes are impracticable to apply to archaic samples because of damage to the ancient DNA. However, the natural process of cytosine deamination, in which unmethylated cytosines decay with time to uracils and methylated cytosines decay to thymines (5), has been shown to carry information regarding general trends of genome methylation (5, 6). Sequencing protocols of ancient DNA trim uracils but leave thymines unprocessed (5). Therefore, positions in which cytosines were methylated premortem are expected to have a higher fraction of reads with thymines compared with unmethylated positions. Consequently, the fraction of CpG→TpG substitutions (hereinafter, C→T ratio) may serve as a proxy for the levels of methylation in ancient DNA (Fig. 1A and fig. S1).

Fig. 1

C→T ratio is a proxy for ancient DNA methylation. (A) Scheme of the processes that enable the reconstruction of ancient methylation by using C→T ratios. (B) C→T ratios in archaic humans are highly correlated with measured methylation in present-day humans. CpGs were binned according to their measured methylation level in human osteoblast cells, and for each bin a mean C→T ratio was computed. Error bars represent one standard deviation.

We computed, for every CpG in the human genome, a C→T ratio in each of the archaic human genomes. Positions that were expected to represent true premortem thymines rather than postmortem deamination were excluded (7). We then compared the Neandertal and Denisovan C→T ratios to experimentally measured methylation in bone (osteoblast) cells of a present-day human (8). The samples were obtained from the femur, costae, and tibia bones of a young female and therefore highly resemble the Neandertal and Denisovan samples [obtained from limb bones of young females (1, 2)]. For both ancient genomes, the correlation between binned present-day methylation and archaic C→T ratios was ~0.98 (7) (Fig. 1B).

In order for deamination to serve as a reliable proxy for DNA methylation, it should display a stochastic behavior that is homogenous along the genome. By using the genome of a ~4000-year-old Eskimo, the sequencing protocol of which provides a direct measure for deamination rates (9), we show that every region in the genome has the same probability to undergo deamination (7). Such homogeneity is expected given the fast postmortem degradation of cellular structure and proteins that are therefore unlikely to bias deamination across tens of thousands of years.

In order to account for the stochastic nature of deamination, we computed for each CpG a mean C→T ratio by using a sliding window over neighboring CpGs (7) (fig. S2 and table S1). The correlation between the C→T ratio and present-day human methylation is apparent also when examining specific genomic regions (Fig. 2A and fig. S3A). Furthermore, our method is highly robust not only to coverage variation, which could potentially bias C→T ratio through reduced coverage in uracil-trimmed regions, but also to potential differences in deamination rate between methylated and unmethylated cytosines (7) (fig. S4).

Fig. 2 Reconstructed archaic methylation relative to present-day human methylation.

(A) An exemplary region (chromosome 7, 6,606,884 to 6,687,956) demonstrating that reconstructed archaic methylation patterns follow present-day human patterns. C→T ratio and measured methylation are color-coded with a gradient from green (unmethylated) to red (methylated). The bottom line shows genes in the region. For visual clarity, osteoblast methylation in present-day humans appears twice, below each of the C→T ratio lanes. (B) Normalized heat maps of present-day human methylation versus archaic reconstructed methylation for all single CpG positions across the genome. The hot region along the diagonal line shows that, in most CpG positions, reconstructed methylation matches measured present-day methylation.

We then scaled the C→T ratios to represent premortem methylation percentages (7). The reconstructed archaic methylation in single CpGs shows strong correlation to measured present-day methylation (Embedded Image in the Denisovan; Embedded Image in the Neandertal) (Fig. 2B). Next, we investigated CpG islands, which tend to be hypomethylated in mammalian genomes (3). The reconstructed methylation in CpG islands falls within the range of values in present-day humans (fig. S3B). We also investigated 3804 housekeeping genes (10), whose promoters are constitutively hypomethylated (4). The reconstructed methylation shows strong hypomethylation and falls within the range of 16 present-day human tissues taken from 20 individuals (11) (Fig. 3A and fig. S3C). Taken together, these results suggest that the C→T ratio can be used as a proxy for premortem methylation and that genome-wide methylation maps can be reconstructed after tens of thousands of years.

Fig. 3

Characteristics of DMRs. (A) Distribution of methylation in housekeeping promoters in various human tissues (blue) and in the Neandertal and the Denisovan (red). The distribution of methylation in the Neandertal and the Denisovan falls within the variability across human tissues. (B) The numbers of present-day human–, Neandertal-, and Denisovan-specific DMRs. (C) Intergenic DMRs colocalize with enhancer and insulator marks. Overlap of intergenic DMRs with insulator mark (CTCF) and with five different enhancer marks (histone modifications H3K4me1 and H3K27ac, histone variant H2A.Z, binding by p300, and DNase I sensitivity). DMRs are divided according to the number of enhancer marks they overlap.

As expected from groups that diverged recently, over 99% of both archaic genomes show no significant methylation differences compared to the present-day human. However, in each archaic human we found ~1100 differentially methylated regions (DMRs) (7). We divided the DMRs according to where the methylation change has likely occurred, that is, into Neandertal-specific, Denisovan-specific, and present-day human–specific DMRs. DMRs that could not be confidently assigned to a specific lineage were tagged as unclassified DMRs (7) (Fig. 3B and table S2). Because housekeeping genes serve fundamental functions in all cells, we expected them to present low rates of methylation changes. Reassuringly, despite constituting ~10% of genes, none of the housekeeping genes harbors DMRs (Embedded Image, hypergeometric test). Last, when comparing differences in methylation to differences in sequence, we observed a significant overlap between the two, as expected if some of the methylation changes were driven by sequence changes (Embedded Image, randomization test).

A difference in methylation between an archaic and a present-day human does not necessarily imply a fixed difference between the human groups. This difference can stem from variability within a population or from the comparison of close, yet not identical, cell types (osteoblasts versus whole bones). Hence, we compared the archaic methylation to 37 bones samples taken from osteoblasts and whole bones (1214). We sought reliable DMRs, in which the archaic methylation significantly differs from that in modern bones. These samples were measured with 27K arrays and provided information for ~5% of DMRs. In most DMRs, archaic methylation was significantly different from the 37 bones and therefore classified as reliable [false discovery rate (Embedded Image) < 0.01, z-test] (tables S2 and S3).

DMRs for which data from bones was not available were compared with 25 samples of 23 tissues from 20 individuals representing both sexes, different ethnicities, and ages ranging from 6 to 88 (8). Over 50% of DMRs had stable methylation levels across modern samples, but were significantly different in the archaic humans, and are thus likely to represent reliable DMRs (7) (Embedded Image (tables S2 and S3).

A key regulator of limb development is the 5′ HOXD cluster (15), which includes five genes (HOXD9 to HOXD13). We identified three DMRs within the HOXD locus (Fig. 4) and two intergenic DMRs located 20- and 25-kb upstream to the cluster in a putative enhancer region (16) (table S2). In both archaic humans, we detected hypermethylation in the HOXD9 promoter and in the HOXD10 gene body, whereas in present-day humans these regions are hypomethylated in all 37 bone samples (methylation = 5 ± 2%; Embedded Image, z-test). In addition, we found the gene body of HOXD9 to be hypermethylated in the Denisovan (Fig. 4 and table S2). These DMRs imply linked regulation of HOXD9 and HOXD10, consistent with observations in mouse (17). Distinctions between present-day and Neandertal limbs include larger femur articulations (18), robust hands and fingers (19), broader elbow and knee joints (18, 20), more-curved radius and femur (18, 20), and shorter limbs (18, 20). All of these features were shown in mouse and human to result from changes in the HOXD cluster in general and in HOXD9 and HOXD10 in particular (15, 2123). Together, these findings suggest that the HOXD cluster might have played a key role in the recent evolution of human limbs.

Fig. 4 The HOXD cluster is hypermethylated in archaic humans.

C→T ratio in archaic humans, and present-day human methylation (see caption of Fig. 1C). Gray rectangles mark DMRs. A minimum of 10 positions with measured osteoblast methylation was required for a region to be denoted a DMR.

Enhancers are often deoxyribonuclease I (DNase I)–sensitive regions bound by p300 and histone variant H2A.Z and marked by the histone modifications H3K4me1 and H3K27ac (24, 25). For insulators, CCCTC binding factor (CTCF) plays a major role as a buffer between active and inactive regions. By using a randomization test, we found a significant enrichment for enhancer marks and for CTCF binding in the intergenic DMRs, with over 80% having these marks (Fig. 3C, fig. S5, and tables S2 and S4). We show that this cannot be explained by length, CG-content, CpG density, or coverage (7). This suggests that many of the intergenic DMRs may affect gene activity through distal elements. Additionally, both of the intergenic DMRs upstream of the HOXD cluster possess enhancer marks (table S2), suggesting that they might have affected the methylation of the HOXD cluster through distal regulation.

DNA methylation affects, as well as is affected by, the binding of transcription factors (TFs) (3, 4). To test whether multiple DMRs may be associated with a change in a single TF, we looked for TF-encoding genes that are differentially methylated and whose targets are enriched with DMRs (7) (tables S5 and S6). In the Neandertal, we identified sterol regulatory element binding protein 1 (SREBP-1), a regulator of lipid homeostasis (Embedded Image; Fisher’s exact test), and aryl hydrocarbon receptor (AHR), a receptor of environmental toxins (Embedded Image). In present-day humans, we identified early growth response 3 (EGR3), which regulates biorhythm and development of muscles, lymphocytes, and neurons (Embedded Image); and Meis homeobox 1 (MEIS1), a regulator of limb development that forms complexes with both HOXD9 and HOXD10 (Embedded Image). This enrichment is independent of CpG density, GC content, and coverage (7). Additionally, the number of TFs that passed this test is significantly higher than would be expected by random (P = 0.018, hypergeometric test). These TFs could provide a lead to the mechanism behind the methylation changes and possibly offer a model of burstlike evolution, where multiple genes changed their activity after a single evolutionary event.

Last, we found that DMR-containing genes in present-day humans are almost twice as likely as genes that do not contain DMRs to be disease-related (18.1% compared with 10.8%, Embedded Image, Fisher’s exact test) (fig. S6 and table S7). More than a third of the disease-linked genes (30/81) are involved in neurological and psychiatric disorders. Here, too, the enrichment is independent of gene length, CpG density, GC content, and coverage (7). This could potentially point to a link between recent changes in regulation and the emergence of diseases or to a scenario where disease-linked genes (i.e., genes in which a change is not embryonic lethal) are possibly more prone to methylation changes.

Current techniques for measuring genome-wide methylation result in the destruction of the examined DNA (3) and are therefore impracticable in the study of ancient DNA. Our method allows for the reconstruction of methylation maps directly from ancient DNA reads without the need to use additional ancient material. We hope that the differences in methylation elucidated here will help to uncover the epigenetic basis for phenotypic differences between present-day and archaic humans and shed light on the role of epigenetics in the recent evolution of our lineage.

Supplementary Materials

www.sciencemag.org/content/344/6183/523/suppl/DC1

Materials and Methods

Figs. S1 to S6

Tables S1 to S10

References (2643)

References and Notes

  1. Methods are available as supplementary materials on Science Online.
  2. Acknowledgments: We thank S. Shifman, Y. Stelzer, M. Choder, O. Amster-Choder, U. Ben-David, M. Meyer, Y. Aaronson, and Y. N. Harari for critical comments. Supported by the Israel Science Foundation FIRST individual grant (ISF 1430/13 to L.C. and E.M.), the Israel Science Foundation (ISF 1252/12, 657/12 to E.M.), a European Union Marie Curie Reintegration grant (IRG-248639 to L.C.), the Abisch-Frenkel Foundation (to E.M.), the European Research Council (ERC-281781 to E.M.), the Presidential Innovation Fund of the Max Planck Society (to S.P.), and the Spanish Health Ministry (Instituto de Salud Carlos III-FIS PI09/539, PI12/615). Software (RoAM, after Reconstruction of Archaic Methylation) can be found at http://carmelab.huji.ac.il/software/RoAM/roam.html and data at http://carmelab.huji.ac.il/data.html. The authors declare no competing financial interests. Requests for data should be addressed to Liran Carmel (liran.carmel@huji.ac.il) and Eran Meshorer (meshorer@huji.ac.il).
View Abstract

Navigate This Article