ARTICLE https://doi.org/10.1038/s41467-020-15020-6 OPEN Differential DNA methylation of vocal and facial anatomy genes in modern humans David Gokhman et al.# Changes in potential regulatory elements are thought to be key drivers of phenotypic divergence. However, identifying changes to regulatory elements that underlie human- 1234567890():,; specific traits has proven very challenging. Here, we use 63 reconstructed and experimentally measured DNA methylation maps of ancient and present-day humans, as well as of six chimpanzees, to detect differentially methylated regions that likely emerged in modern humans after the split from Neanderthals and Denisovans. We show that genes associated with face and vocal tract anatomy went through particularly extensive methylation changes. Specifically, we identify widespread hypermethylation in a network of face- and voice- associated genes (SOX9, ACAN, COL2A1, NFIX and XYLT1). We propose that these repression patterns appeared after the split from Neanderthals and Denisovans, and that they might have played a key role in shaping the modern human face and vocal tract. # A full list of authors and their affiliations appears at the end of the paper. NATURE COMMUNICATIONS | (2020)11:1189 | https://doi.org/10.1038/s41467-020-15020-6 | www.nature.com/naturecommunications 1 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15020-6 T he advent of high-coverage ancient genomes of archaic collectively referred to as anatomically modern humans (AMHs), humans (Neanderthal and Denisovan) introduced the while the Neanderthal and Denisovan are referred to as archaic possibility to identify the genetic basis of some unique humans. As an outgroup, we produced methylomes of six modern human traits1. A common approach is to carry out chimpanzees: one WGBS, one reduced representation bisulfite sequence comparisons and detect non-neutral sequence changes. sequencing (RRBS) and four 850K methylation arrays (Supple- However, out of ~30,000 substitutions and indels that reached mentary Fig. 1, Supplementary Tables 1, 2). Together, these data fixation in modern humans, less than 100 directly alter amino establish a unique and comprehensive platform to study DNA acid sequence1, and as of today, our ability to estimate the bio- methylation dynamics in recent human evolution (Supplemen- logical effects of the remaining ~30,000 noncoding changes is tary Data 1). very restricted. Whereas many of them are probably nearly neutral, many others may affect gene function, especially those in regulatory regions such as promoters and enhancers. Such Identification of DMRs. We developed a DMR-detection changes in regulatory elements may have a sizeable impact on method for ancient methylomes, which accounts for potential human evolution, as alterations in gene regulation are thought to noise introduced during reconstruction, as well as differences in underlie much of the phenotypic variation between closely related coverage and deamination rates. To minimize the number of false groups2. Because of the limited ability to interpret noncoding positives and to identify DMRs that are most likely to have a variants, direct examination of regulatory layers such as DNA regulatory effect, we applied a strict threshold of >50% difference methylation has the potential to enhance our understanding of in methylation across a minimum of 50 CpGs. This also filters out the evolutionary origin of human-specific traits far beyond what environmentally-induced DMRs which typically show small can be achieved using sequence comparison alone3. methylation differences and limited spatial scope10. Using this In order to gain insight into changes in regulatory elements method, we identified 9679 regions overall that showed methy- that might underlie human evolution, we previously developed a lation differences between any of the high-quality representative method to reconstruct DNA methylation maps of ancient gen- methylomes of the Denisovan, the Altai Neanderthal, and the omes based on analysis of patterns of damage to ancient DNA4. Ust’-Ishim anatomically modern human. These regions do not We used this method to reconstruct the methylomes of a necessarily represent evolutionary differences between the human Neanderthal and a Denisovan, which were then compared to a groups. Rather, many of them could be attributed to variability partial methylation map of a present-day osteoblast cell line. within populations or to factors separating the three individuals However, the ability to identify differentially methylated regions (e.g., Ust’-Ishim is a male whereas the archaic humans are (DMRs) between the human groups was constrained by the females). One common approach to minimize such false positives incomplete reference map (providing methylation information is to match samples in as many parameters as possible: age, sex, for ~10% of CpG sites), differences in outputs of sequencing bone type, technology, and disease state. However, this is rarely platforms, lack of an outgroup, and a restricted set of skeletal fully accomplished, and samples often remain unmatched in one samples (see Supplementary Methods). or more categories. A stricter approach is to leverage the varia- To study the evolutionary dynamics of DNA methylation along bility between samples to estimate the maximum contribution of the hominin tree on a larger scale, here we use 69 skeletal DNA their confounders. To this end, and to minimize the contribution methylation maps from modern humans, archaic humans, and of within-population variability, we used the 59 additional human chimpanzees to identify 588 genes whose methylation state is maps to filter out regions where variability in methylation is unique to modern humans. We then analyze the function of these detected. Importantly, our samples come from both sexes, from genes by investigating their known anatomical effects, and vali- individuals of various ages and ancestries, from patients and date this using over 50 orthogonal tests and controls. We find that healthy individuals, and from a variety of skeletal parts (femur, the most extensive DNA methylation changes are observed in skull, phalanx, tooth, and rib). We adopted a conservative genes that affect vocal and facial anatomy, and that this trend approach, whereby we take only loci where methylation in one appears to be unique to modern humans. hominin group is found completely outside the range of methy- lation in the other groups (Fig. 1a). Hence, as the differences between groups is considerably larger than the contribution of Results these factors, this procedure is expected to account for these Generating DNA methylation maps. We reconstructed ancient potentially confounding factors, and the remaining DMRs are DNA methylation maps of eight individuals: in addition to the expected to represent true evolutionary differences (Fig. 1a–c and previously published Denisovan and Altai Neanderthal methyla- Supplementary Figs. 2, 3, Supplementary Methods). This step tion maps4, we reconstructed the methylomes of the Vindija resulted in a set of 7649 DMRs that discriminate between the Neanderthal (~52 thousand years ago, kya)5, and three anatomi- human groups, which we ranked according to their significance cally modern humans: the Ust'-Ishim individual (~45 kya, Wes- level. For the top candidates (see later), we also ascertain that the tern Siberia)6, the Loschbour individual (~8 kya, Luxemburg)7, differential methylation is observed when matching samples from and the Stuttgart individual (~7 kya, Germany)7. We also the same bone type, age, disease state, and technology. sequenced to high-coverage and reconstructed the methylomes of Next, using the chimpanzee samples, we were able to the La Braña 1 individual from Spain (~8 kya, 22×) (which was determine for 2825 of these DMRs the lineage where the previously sequenced to low-coverage)8 and an individual from methylation change occurred (Fig. 1d). Of these DMRs, 873 are Barçın Höyük, Western Anatolia, Turkey (I1583, ~8.5 kya, 24×), AMH-derived, 939 are archaic-derived, 443 are Denisovan- which was previously sequenced using a capture array9. derived, and 570 are Neanderthal-derived (Fig. 2a and Supple- To this set we added 52 publicly available partial bone mentary Fig. 3, Supplementary Data 2). To study the derived methylation maps from present-day individuals, produced using biology of AMHs, and to focus on DMRs that are based on the 450K methylation arrays (see Supplementary Methods). To most extensive set of maps, we concentrated on the 873 AMH- obtain full present-day bone maps, we produced whole-genome derived DMRs. We found that these DMRs are located 58× closer bisulfite sequencing (WGBS) methylomes from the femur bones to AMH-derived sequence changes than expected by chance of two individuals (Bone1 and Bone2, Supplementary Fig. 1). (0.092 Mb vs. a median of 5.3 Mb, P < 10−5, permutation test, Hereinafter, ancient and present-day modern humans are Fig. 2b). This suggests that some of the methylation changes 2 NATURE COMMUNICATIONS | (2020)11:1189 | https://doi.org/10.1038/s41467-020-15020-6 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15020-6 ARTICLE a Lineage-specific DMR Modern human Denisovan Neanderthal Chimpanzee 100 MH-derived DMR Femur Skull Femora Teeth Methylation (%) Femora Finger Femora Rib Femur Toe 0 450 K WGBS Recon. Recon. Recon. 850 K WGBS b Limb-specific DMR Modern human Denisovan Neanderthal Chimpanzee 100 Femora Femur Femur Methylation (%) Femora Femora Toe Finger Teeth Skull Rib 0 450 K WGBS Recon. Recon. Recon. 850 K WGBS c Sex-specific DMR Modern human Denisovan Neanderthal Chimpanzee 100 Males Males Males Methylation (%) Unknown Female Females Females Female Females 0 450 K WGBS Recon. Recon. Recon. 850 K WGBS d AMH-derived Archaic-derived Neanderthal-derived Undefined might be associated with cis-regulatory sequence variants that up to 5 kb upstream. The 873 AMH-derived DMRs are linked to arose along the AMH lineage. 588 AMH-derived DMGs (Supplementary Data 2). To gain insight into the function of these DMGs, we first analyzed their Face and voice-associated genes are derived in AMHs. We gene ontology (GO). As expected from genes that show differ- defined differentially methylated genes (DMGs) as genes that ential methylation in the skeleton between human groups, AMH- overlap at least one DMR along their body or in their promoter, derived DMGs are enriched with terms associated with the NATURE COMMUNICATIONS | (2020)11:1189 | https://doi.org/10.1038/s41467-020-15020-6 | www.nature.com/naturecommunications 3 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15020-6 Fig. 1 Variability filtering and lineage assignment. a Methylation levels across AMH, Denisovan, Neanderthal, and chimpanzee samples in DMR#278 (chr4:38,014,896–38,016,197). This is an example of a lineage-specific DMR, defined as a locus in which all samples of a group are found outside the range of methylation in the other groups. Chimpanzee samples were used during the following step of lineage assignment. b A putative limb-specific DMR (chr3:14,339,371–14,339,823) which was removed from the analysis, as it does not comply with our definition of lineage-specific DMRs. Femur, toe, and finger samples are hypermethylated compared to other skeletal elements. Toe and finger are found at the bottom range of limb samples, suggesting some variation in this locus within limb samples too. c A putative sex-specific DMR (chr3:72,394,336–72,396,901) which was removed from the analysis. Males are hypermethylated compared to females. d Lineage assignment using chimpanzee samples. Only DMRs that passed the previous variability filtering steps were analyzed. Each bar at the tree leaves represents a locus in a sample. Methylation levels of the locus in each sample are marked with red (methylated) and green (unmethylated). The lineage where the methylation change has likely occurred (by parsimony) is marked by a star. Branch lengths are not scaled. a Chimpanzee Neanderthal Denisovan AMH c d 570 443 390–440 873 kya 939 520–630 kya ~6 mya b 450 P < 10–5 simulations Number of DMRs (0.092 Mb) 0 0 0.5 3.5 7 Median distance to closest AMH-derived single nucleotide change (Mb) ×2.11 ×3.82 ×1 ×1 ×0.47 ×0.35 Fig. 2 Genes affecting voice and face are the most over-represented within AMH-derived DMRs. a The number of DMRs that emerged along each of the human branches. Divergence times are in thousands of years ago (kya). b Distribution of median distances (turquoise) of DMRs to randomized single nucleotide changes that separate AMHs from archaic humans and chimpanzees. Genomic positions of single nucleotide changes were allocated at random. This was repeated 10,000 times. Red arrow marks the observed distance of DMRs, showing that they tend to be significantly closer to AMH-derived single nucleotide changes than expected by chance. This suggests that some of these sequence changes might be associated with the changes in methylation. c A heat map representing the level of enrichment of each anatomical part within the AMH-derived DMRs. Only body parts that are significantly enriched (FDR < 0.05) are colored. Three skeletal parts are significantly over-represented: the face, pelvis, and larynx (voice box, marked with arrows). d Enrichment levels of anatomical parts within the most significant (top quartile, Q statistic) AMH-derived DMRs, showing a more pronounced enrichment of genes affecting vocal and facial anatomy. skeleton (e.g., endochondral bone morphogenesis, trabecula ORGANizer, which links human genes to the organs they morphogenesis, and palate development). Also notable are terms phenotypically affect11. Unlike tools that use GO terms or RNA associated with the muscular, cardiovascular, and nervous sys- expression data, Gene ORGANizer is based entirely on curated tems (Supplementary Data 3). gene-disease and gene-phenotype associations from monogenic To acquire a more precise understanding of the possible diseases. It relies on direct phenotypic observations in human functional consequences of these DMGs, we used Gene patients whose conditions are associated with known gene 4 NATURE COMMUNICATIONS | (2020)11:1189 | https://doi.org/10.1038/s41467-020-15020-6 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15020-6 ARTICLE perturbations. Hereinafter, we refer to genes as associated with an none of the other branches shows enrichment of the larynx or the organ if they have been shown to have a phenotypic effect on that vocal folds; Neanderthal-derived and Denisovan-derived DMGs organ in patients where this gene is dysfunctional. An enrichment show no significant enrichment in any organ, and archaic-derived or depletion in Gene ORGANizer is detected if the group of genes DMGs are over-represented in the jaws, lips, limbs, scapulae, and analyzed shows a significant deviation in the organs they are spinal column, but not in the larynx or vocal folds (Supplemen- known to phenotypically affect, compared to the rest of the tary Fig. 4e, Supplementary Data 4). In addition, DMRs that genome. Using Gene ORGANizer, we found 11 organs that are separate chimpanzees from all humans (archaic and modern, over-represented within the 588 AMH-derived DMGs, eight of Supplementary Data 2) do not show enrichment of genes which are skeletal parts that can be divided into three regions: the associated with the larynx or face, compatible with the notion face, larynx (voice box), and pelvis (Fig. 2c and Supplementary that this trend emerged along the AMH lineage. Data 4). The strongest enrichment was observed in the laryngeal Taken together, we conclude that DMGs that emerged along region (×2.11 and ×1.68, FDR = 0.017 and 0.046, for the the AMH lineage are uniquely enriched in genes associated with vocal folds (vocal cords) and larynx, respectively), followed by the voice and face, and that this is unlikely to be an artifact of (a) facial and pelvic structures, including the teeth, forehead, jaws, inter-individual variability resulting from age, sex, disease, or and pelvis. Interestingly, the face and pelvis are considered the bone type; (b) significance level of DMRs; (c) the reconstruction most morphologically divergent regions between Neanderthals or DMR-detection processes; (d) number of samples used; (e) and AMHs12 and our results reflect this divergence through gene pleiotropic effects; (f) the types of methylation maps used; (g) the regulation changes. To gain orthogonal evidence for the comparison of skeletal methylomes; (h) gene length distribution; enrichment of the larynx and face within these AMH-derived or (i) biological factors such as cellular composition and DMGs, we carried out a number of additional analyses: First, we developmental state. analyzed gene expression patterns and found that the supralar- Our analyses identified 56 DMRs in genes associated with the yngeal vocal tract (the pharyngeal, oral, and nasal cavities, where facial skeleton, and 32 in genes associated with the laryngeal sound is filtered to specific frequencies) is the most enriched body skeleton. The face- associated genes are known to shape mainly part (1.7× and 1.6×, FDR = 5.6 × 10−6 and FDR = 7.3 × 10−7, for the protrusion of the lower and midface, the size of the nose, and the pharynx and larynx, respectively, hypergeometric test, the slope of the forehead. Interestingly, these traits are considered Supplementary Data 3). Second, 44 of the AMH-derived DMRs some of the most derived between Neanderthals and AMHs12. overlap previously reported putative enhancers of human The larynx-associated genes have been shown to underlie various craniofacial developmental genes (5.1× compared to expected, phenotypes in patients, ranging from slight changes to the pitch P < 10−4, permutation test)13,14. Third, Palate development is the and hoarseness of the voice, to a complete loss of speech ability11 third most enriched GO term among AMH-derived DMGs (Supplementary Data 5). These phenotypes were shown to be (Supplementary Data 3). Fourth, DMGs significantly overlap driven primarily by alterations to the laryngeal and vocal tract genes associated with craniofacial features in the GWAS catalog15 skeleton. Methylation patterns in differentiated cells are often (P = 3.4 × 10−4, hypergeometric test). established during earlier stages of development, and the closer To test whether this enrichment remains if we take only the two tissues are developmentally, the higher the similarity between most confident DMRs, we limited the analysis to DMGs where their methylation maps3,16,17. This is also evident in the fact that the most significant DMRs are found (top quartile, Q statistic). DMRs identified between species in one tissue often exist in other Here, the over-representation of voice-affecting genes is even tissues as well16. Importantly, the laryngeal skeleton, and more pronounced (2.82× and 2.26×, for vocal folds and larynx, particularly the arytenoid cartilage to which the vocal folds are respectively, FDR = 0.028 for both, Fig. 2d and Supplementary anchored, share an origin from the somatic layer of the lateral Data 4). plate mesoderm with the cartilaginous tissue of the limb bones Next, we reasoned that skeleton-associated genes might be prior to their ossification. Thus, it is likely that many of the DMRs over-represented in analyses that compare bone DNA methyla- identified here between limb samples also exist in their closest tion maps, hence introducing potential biases. To test whether tissue—the laryngeal skeleton. This is further supported by the this enrichment might explain the over-representation of the observation that these DMGs are consistent across all examined larynx, face, and pelvis, we compared the fraction of genes skeletal samples, including skull, femur, rib, tibia, and tooth. associated with these organs within all skeletal genes to their Furthermore, we directly measured methylation levels in a subset fraction within the skeletal genes in the AMH-derived DMGs. We of the DMRs in primary chondrocytes and show that their found that genes associated with the face, larynx, and pelvis are patterns extend to these cells as well (see below). significantly over-represented even within skeletal AMH-derived DMGs (P = 1.0 × 10−5, P = 1.3 × 10−3, P = 2.1 × 10−3, P = 0.03, for vocal folds, larynx, face, and pelvis, respectively, hypergeo- Extensive changes within face and voice-associated genes. The metric test). Additionally, using a permutation test, we found that results above suggest that methylation levels in many face- the enrichment levels within AMH-derived DMGs are signifi- associated and voice-associated genes have changed in AMHs cantly higher than expected by chance for the laryngeal and facial since the split from archaic humans, but they do not provide regions, but not for the pelvis (P = 8.0 × 10−5, P = 3.6 × 10−3, information on the extent of changes within each gene. To do so, P = 8.2 × 10−4, and P = 0.115, for vocal folds, larynx, face and we scanned the genome in windows of 100 kb and computed the pelvis, respectively, permutation test, Supplementary Fig. 4, fraction of CpGs which are differentially methylated in AMHs Supplementary Methods). Thus, we found that the enrichment (hereinafter, AMH-derived CpGs). We found that the extent of in the facial and laryngeal regions is not a by-product of a general changes within voice-associated DMGs is most profound, more enrichment in skeletal parts, and we hereinafter focus on genes than 2× compared to other DMGs (0.132 vs. 0.055, FDR = 2.3 × associated with these two regions. 10−3, t-test, Supplementary Data 6). Face- associated DMGs also Finally, we ruled out the options that our DMR-detection present high density of AMH-derived CpGs (0.079 vs. 0.055, algorithm, number of samples, filtering process or biological FDR = 2.8 × 10−3, t-test). In archaic-derived DMGs, on the other factors such as gene length, cellular composition, pleiotropy or hand, the extent of changes within voice-associated and face- developmental stage might underlie the enrichment of these associated genes is not different than expected (FDR = 0.99, organs (see Supplementary Methods). Perhaps most importantly, t-test, and Supplementary Data 6). To control for possible biases, NATURE COMMUNICATIONS | (2020)11:1189 | https://doi.org/10.1038/s41467-020-15020-6 | www.nature.com/naturecommunications 5 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15020-6 we repeated the analysis using only the subset of DMRs in genes phenotypes are also known as traits that differ between modern associated with the skeleton. Here too, we found that voice- and archaic humans. We found that four of the top five most associated AMH-derived DMGs present the highest density of differentially methylated genes (XYLT1, NFIX, ACAN, and changes (2.5× for vocal folds, 2.4× for larynx, FDR = 1.4 × 10−3 COL2A1) are in the top 100 genes with the highest fraction of for both, t-test, Supplementary Data 6), and face-associated divergent traits between Neanderthals and AMHs. Remarkably, DMGs also exhibit a significantly elevated density of changes COL2A1, the most divergent gene in its methylation patterns, is (1.4×, FDR = 0.04, t-test). also the most divergent in its phenotypes: no other gene in the We also found that compared to other AMH-derived DMRs, genome is associated with as many divergent traits between DMRs in voice-associated and face-associated genes tend to modern humans and Neanderthals32 (63 traits, Supplementary be 40% closer to candidate positively selected loci in AMHs18 Data 7, see Supplementary Methods). This suggests that these (P < 10−4, permutation test). extensive methylation changes are possibly linked to phenotypic Strikingly, when scanning the genome for hotspots of divergence between archaic and AMHs. methylation changes, all top five skeleton-related loci are found within genes known to associate with lower and midfacial NFIX methylation patterns suggest downregulation in AMHs. protrusion, as well as the voice (ACAN, SOX9, COL2A1, XYLT1, In order to investigate how methylation changes affect expression and NFIX)11,19 (Fig. 3 and Supplementary Fig. 4f). This is levels, we scanned the DMRs to identify those whose methylation particularly surprising considering that genome-wide, less than levels are strongly correlated with expression across 22 human 2% of genes (345) are known to affect the voice, ~3% of genes tissues33. We found 90 such AMH-derived DMRs (FDR < 0.05, (726) are known to associate with lower and midfacial protrusion, Supplementary Data 2). DMRs in voice-associated genes are and less than 1% (182) are known to associate with both11,19. significantly more likely to be correlated with expression com- The three skeletal DMGs with the highest density of AMH- pared to other DMRs (2.05×, P = 6.65 × 10−4, hypergeometric derived CpGs are the extra-cellular matrix genes ACAN and test). Particularly noteworthy is NFIX, one of the most derived COL2A1, and their key regulator SOX9, which together form a genes in AMHs (ranked 5th among DMGs affecting the skeleton, network that regulates skeletal growth, the transition from Fig. 3a, b). NFIX contains two DMRs (#24 and #167, Fig. 5a), and cartilage to bone, and spatio-temporal patterning of skeletal in both, methylation levels are tightly linked with expression development, including the facial and laryngeal skeleton in (correlation of 81.7 and 73.8%, FDR = 3.5 × 10−6 and 8.6 × 10−5, humans19,20 and mouse21. SOX9 was also shown to be one of the respectively, Pearson’s r, Fig. 5b). In fact, NFIX is one of the top top genes associated with variation in craniofacial morphology ten DMGs with the most significant correlation between within-AMHs22. SOX9 is regulated by a series of upstream methylation and expression in human. The association between enhancers identified in mouse and human23. In human skeletal NFIX methylation and expression was also shown previously samples, hypermethylation of the SOX9 promoter was shown to across several mouse tissues34,35. To further examine this, we downregulate its activity, and consequently its targets24. This was investigated a dataset of DNMT3A-induced methylation of also demonstrated repeatedly in non-skeletal tissues of human MCF-7 cells. Forced induction of methylation in this human25,26 and mouse27,28. We found substantial hypermethyla- study was sufficient to repress NFIX expression by over 50%, tion in AMHs in the following regions: (a) the SOX9 promoter; placing NFIX as one of the genes whose expression is most (b) seven of its proximal and distal skeletal and skeletal affected by hypermethylation36 (ranked in the 98th percentile, progenitor enhancers23; (c) the targets of SOX9: ACAN (DMR FDR = 1.28 × 10−6). We further validated the hypermethylation #80) and COL2A1 (DMR #1, the most significant AMH-derived of NFIX across the skeleton by comparing four human cranial DMR, which spans 32 kb and covers almost the entire COL2A1 samples to four chimpanzee cranial samples through bisulfite- gene, from its 1st intron to its 54th exon and 3′UTR region); and PCR (P = 0.01, t-test, Supplementary Fig. 5, Supplementary (d) an upstream lincRNA (LINC02097). Notably, regions (a), (b), Data 1, Supplementary Table 3, Supplementary Methods). and (d) overlap the longest DMR on the AMH-derived DMR list, Together, these findings suggest that the observed hypermethy- spanning 35,910 bp (DMR #11, Fig. 4). Additionally, a more lation of NFIX in AMHs reflects downregulation that emerged distant putative enhancer, located 345 kb upstream of SOX9, was along the AMH lineage. Indeed, we found that NFIX, as well as shown to bear strong active histone modification marks in SOX9, ACAN, COL2A1, and XYLT1 are hypermethylated in chimpanzee craniofacial progenitor cells; whereas, in humans human femora compared to baboon37 (P = 1.4 × 10−5 and P = these marks are almost absent (~10× lower than chimpanzee, 8.1 × 10−9, compared to baboon femora bone and cartilage, suggesting downregulation, Fig. 4b)13. Importantly, human and respectively, t-test). Also, all five genes show significantly reduced chimpanzee non-skeletal tissues (i.e., brain and blood) exhibit expression in humans compared to mice (Fig. 5c). Taken toge- very similar methylation patterns in these genes, suggesting the ther, these observations suggest that DNA methylation is a pri- DMRs are skeleton-specific. Finally, the amino acid sequence mary mechanism in the regulation of NFIX, and serves as a good coded by each of these genes is identical across the hominin proxy for its expression. Interestingly, NFI proteins were shown groups1, suggesting that the observed changes are purely to bind the upstream enhancers of SOX9 (ref. 38), hence sug- regulatory. Together, these observations support the idea that gesting a possible mechanism to the simultaneous changes in the SOX9 became downregulated in AMH skeletal tissues, likely five top genes we report. followed by downregulation of its targets: ACAN and COL2A1. XYLT1, the 4th highest skeleton-related DMG, is an enzyme involved in the synthesis of glycosaminoglycan. Loss-of-function Discussion mutations, hypermethylation of the gene and its consequent We have shown here that genes associated with vocal and facial reduced expression underlie the Desbuquois dysplasia skeletal anatomy have different DNA methylation patterns in recent syndrome, which was shown to affect the cartilaginous structure AMH, compared to Neanderthals and Denisovans. The differ- of the larynx, and drive a retraction of the face29,30. Very little ences in DNA methylation are manifested both in the number of is known about XYLT1 regulation, but interestingly, in zebrafish divergent genes and in the extent of changes within each gene. it was shown to be bound by SOX9 (ref. 31). Notably, the DMRs we report capture substantial methylation To quantitatively investigate the potential phenotypic con- changes (over 50% between at least one pair of human groups), sequences of these DMGs, we tested what fraction of their known span thousands or tens of thousands of bases, including in 6 NATURE COMMUNICATIONS | (2020)11:1189 | https://doi.org/10.1038/s41467-020-15020-6 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15020-6 ARTICLE a AMH-derived AMH-derived (skeletal) Archaic-derived Simulated DMRs P = 1.5 × 10–4 P = 1.4 × 10–4 P = 0.44 P = 0.42 ACAN 0.582 ACAN 0.582 0.212 0.0135 SOX9 SOX9 TNS3 COL2A1 COL2A1 XYLT1 LINC00673 NFIX SYNPO CCRL2 XYLT1 OPTC GALNT18 NFIX Fraction of differentially methylated CpGs 0.002 0.008 0.001 0.0002 b 0.6 chr15 ACAN 0 0 102,531,392 SOX9 0.35 chr17 Fraction of differentially methylated CpGs 0 0 81,195,210 0.35 chr12 COL2A1 0 0 133,851,895 XYLT1 0.25 chr16 0 0 90,354,753 0.25 chr19 NFIX 0 0 59,128,983 promoters and enhancers. Many of these methylation changes are larynx20,29,39–42. The results we report, which are based on associated with changes in expression. We particularly focused on ancient DNA methylation patterns, provide means to analyze the changes that might affect regulation of the five most derived genetic mechanisms that underlie the evolution of the human face skeletal genes on the AMH lineage: SOX9, ACAN, COL2A1, and vocal tract. XYLT1, and NFIX, whose downregulation was shown to underlie One of the limitations of analyzing regulatory maps is their a retracted face, as well as changes to the structure of the tissue-specific, sex-specific, age-specific, disease-specific, and NATURE COMMUNICATIONS | (2020)11:1189 | https://doi.org/10.1038/s41467-020-15020-6 | www.nature.com/naturecommunications 7 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15020-6 Fig. 3 The extent of differential methylation is highest among genes associated with the voice. a Within each lineage, the fraction of differentially methylated CpGs was computed as the number of derived CpGs per 100 kb centered around the middle of each DMR. DMRs were ranked according to the fraction of derived CpG positions in their vicinity. DMRs in genes associated with the voice are marked with red lines. In AMHs, DMRs in voice-affecting genes tend to be ranked significantly higher. Although known voice-associated genes comprise less than 2% of the genome, three of the top five AMH- derived DMRs, and all top five skeleton-related AMH-derived DMRs are in genes known to associate with the voice. In archaic-derived DMRs and in simulated DMRs, voice-associated genes do not show higher ranking compared to the rest of the DMGs. t-test P-values are shown for each group. b The fraction of differentially methylated CpGs along the five chromosomes containing ACAN, SOX9, COL2A1, XYLT1, and NFIX. In each of these chromosomes, the most extensive changes are found within these genes. All five genes control facial projection and the development of the larynx. a 0.07 SOX9 Frequency DMR #11 0 Mean methylation 0 10 20 30 40 50 60 70 80 90 100 Chimp (850 K) Denisovan Neanderthal Chimp (850 K) Chimp (850 K) Chimp (850 K) AMH (bone1) Chimp (WGBS) AMH (Loschbour) MH (bone2) AMH (Stuttgart) AMH (La Brana 1) AMH (I1583) AMH (Ust’-Ishim) 0.5 ACAN Frequency DMR #80 0 Mean methylation 0 10 20 30 40 50 60 70 80 90 100 Neanderthal Chimp (850 K) Chimp (850 K) Denisovan Chimp (850 K) Chimp (850 K) Chimp (WGBS) AMH (bone1) AMH (Stuttgart) AMH (Loschbour) AMH (bone2) AMH (La Brana 1) AMH (I1583) AMH (Ust’-Ishim) 0.09 COL2A1 Frequency DMR #1 0 Mean methylation 0 10 20 30 40 50 60 70 80 90 100 Neanderthal Denisovan Chimp (WGBS) AMH (Stuttgart) AMH (Loschbour) AMH (bone1) AMH (bone2) AMH (I1583) AMH (La Brana 1) AMH (Ust’-Ishim) b Methylation difference AMHs vs archaic: –4 +27 +19 +1 +13 +10 +2 +25 –6 +10 –3 +47 +9 +14 +31 +6 N/A AMHs vs chimp: –9 –16 +18 +19 +23 –12 +48 +42 –16 +44 +20 +38 +20 +4 +29 +5 chr17 69,772,395 LINC01152 LINC02097 SOX9 70,145,395 Fig. 4 Hypermethylation of SOX9, ACAN, and COL2A1 in AMHs. a Methylation levels in the AMH-derived DMRs in SOX9, ACAN, and COL2A1. AMH samples are marked with red lines, archaic human samples are marked with blue lines and chimpanzee samples are marked with gray lines. The distribution of methylation across 52 AMH samples (450K methylation arrays) is presented in red. b SOX9 and its upstream regulatory elements. AMH-derived DMRs are marked with red rectangles. Previously identified putative enhancers are marked with circles: human craniofacial (orange), human craniofacial, chimpanzee-biased (blue), human skeletal (yellow), human non-skeletal (dark gray), and mouse (light gray). Numbers above skeletal enhancers show the difference in mean bone methylation between AMHs and archaic humans (top) and between AMHs and chimpanzee (bottom). Across almost all SOX9 enhancers, AMHs are hypermethylated compared to archaic humans and the chimpanzee. technology-specific patterns. Although the vast majority of loci to measure methylation: deamination of cytosines. However, in show stable patterns regardless of these factors, some loci are WGBS, the use of bisulfite results in ~99% deamination rate, nevertheless affected by them43. Perfectly matching two human whereas in ancient reconstruction the average rate is usually samples across all of these factors is rarely achievable. Therefore, lower than 5%43. To account for this, we used a window of 25 we used a stricter approach, which results in fewer discoveries, CpGs, which increases power at the cost of resolution. Therefore, but is able to leverage the variability introduced by tissue, sex, age, in this study, we are only able to detect DMRs that are larger than disease, and technology to identify loci where these factors are the window used. To detect more subtle changes, genomes with unlikely to underlie the observed methylation differences. For the higher deamination rates, or merging several genomes together top candidates, we also matched samples for these factors and are required. showed that the differential methylation is still observed. Humans are distinguished from other apes in their unique A central confounder of reconstructed methylation maps is capability to communicate through speech. This capacity has their resolution. Although their resolution is substantially better been attributed not only to neural changes, but also to config- than that of methylation arrays, it is still inferior to that of WGBS. urational alterations to the vocal tract44, as humans have a larynx Both WGBS and ancient reconstruction rely on the same process that is located particularly low among primates44 (Fig. 6a), and 8 NATURE COMMUNICATIONS | (2020)11:1189 | https://doi.org/10.1038/s41467-020-15020-6 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15020-6 ARTICLE a AM Hs Ch AH im s pa nz ee NFIX 25 CpGs I1583 I1583 La Brana 1 La Brana 1 Loschbour Loschbour Stuttgart Stuttgart Bone2 Bone2 Bone1 Bone1 Ust-Ishim Ust’-Ishim Denisovan Denisovan Neanderthal Neanderthal Chimpanzee Chimpanzee chr19 DMR #167 chr19 DMR #24 NFIX NFIX 13,145,588 13,168,871 13,175,658 13,198,259 b 70 DMR #167 70 DMR #24 c –6 P = 2.8 × 10 –5 P = 3.0 × 10 –4 P = 7.2 × 10 P = 3.1 × 10 –17 –22 P = 6.1 × 10 –5 –6 rho = –0.74, P = 8.6 × 10 rho = –0.82, P = 3.5 × 10 100 N N Expression percentile Expression (FPKM) Expression (FPKM) Ao D D Ao MC SC BM Hi MC Ov Es Ad SC Hi Ov AN BM Es AN SN Lu SN Ad Lu B1 RA RA B1 B2 UI Ch LV B2 UI Co SI AG Te LV Sp Pa Co Li SI AG Te Pa Li Th Sp Th 0 0 0 0 Methylation (%) 100 0 Methylation (%) 100 NFIX COL2A1 ACAN SOX9 XYLT1 Ao - Aorta Lu - Lung BM - Breast Myoepithelium Ad - Adipose MC - Midfrontal Cortex B1 - Bone1 (AMH) Li - Liver Ov - Ovary SI - Small Intestine Te - Testis Ch - Chondrocytes B2 - Bone2 (AMH) Hi - Hippocampus Pa - Pancreas Th - Thymus Co - Colon N - Neanderthal Es - Esophagus RA - Right Atrium Sp - Spleen AN - Adipocyte Nuclei D - Denisovan LV - Left Ventricle SC - Sigmoid Colon AG - Adrenal Gland SN - Substantia Nigra UI - Ust’-Ishim (AMH) Fig. 5 DNA methylation and expression pattern of NFIX in modern and archaic humans. a Methylation levels along NFIX, color-coded from green (unmethylated) to red (methylated). Methylation levels around the two AMH-derived DMRs (#24 and #167) are shown in the zoomed-in panels. These two DMRs represent the regions where the most significant methylation changes are observed, but hypermethylation of NFIX in AMHs can be seen throughout the entire gene body. Chimpanzee and present-day samples were smoothed using the same sliding window as in ancient samples to allow easier comparison. The inferred schematic regulatory evolution of NFIX is shown using a phylogenetic tree to the left of the top panel. Star marks the shift in methylation from unmethylated (green) to methylated (red). b Methylation levels in DMRs #167 and #24 vs. expression levels of NFIX across 22 AMH tissues (gray). In both DMRs, higher methylation is significantly associated with lower expression of NFIX. Ust’-Ishim, Bone1 and Bone2 methylation levels (red) are plotted against mean NFIX expression across 13 osteoblast lines. Neanderthal and Denisovan methylation levels (green) are plotted against their predicted expression levels, based on the extrapolated regression line (dashed). Standard errors are marked with dotted lines. The expression levels of NFIX in Neanderthal and Denisovan shown in these graphs are extrapolated (green dots). c Box plots of expression levels of NFIX, COL2A1, ACAN, SOX9, and XYLT1 in 89 AMH samples (red) and four mouse samples (green) from appendicular bones (limbs and pelvis). Central line shows mean, box borders show 25th and 75th percentiles, whiskers extend to the most extreme data points not considered outliers, crosses show outliers. Expression levels were converted to percentiles based on the level of gene expression compared to the rest of the genome in each sample. t-test P-values are shown for each gene. because phonetic range is determined by the different config- current study. Thus, it is still to be determined how the positional urations that the vocal tract can produce. The roles of divergent and structural changes to the larynx emerged and what was the anatomy vs. cognition in our speech skills are still debated45,46, role of the reported methylation changes in this process. For an and some propose that even with a human brain, other apes could in-depth review of the anatomy of vocalization and speech, see not reach the human level of articulation and phonetic range44,47. refs. 44,46. Regardless of its potential role in speech, the process of laryngeal A longstanding question is whether Neanderthals and AMHs descent is developmentally and evolutionarily associated with share similar vocal tract anatomy50,51. Attempts to answer this facial retraction48,49. In this regard, the observation that the top question based on morphological differences have proven hard, as five skeletal DMRs are found in genes that are associated with the larynx is mostly composed of soft tissues (e.g., cartilage), both facial protrusion and the anatomy of the larynx suggests that which do not survive long after death. The only remnant from the these two processes might have been genetically linked, though Neanderthal laryngeal region is the hyoid bone, which is detached the interaction between the two is still to be determined, as their from the rest of the skull51. Based on this single bone, or on exact developmental pathways are beyond the scope of the computer simulations and tentative vocal tract reconstructions, it NATURE COMMUNICATIONS | (2020)11:1189 | https://doi.org/10.1038/s41467-020-15020-6 | www.nature.com/naturecommunications 9 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15020-6 a Larynx Hyoid Epiglottis Tongue b NFIX activity level: Low High Healthy AMH Neanderthal Syndromes Marshall-Smith/Malan syndromes Healthy AMH Neanderthal Phenotypes: SKULLCAP HIGH LOW FOREHEAD PROMINENT RECEDING CRANIUM LARGER LARGER ORBITS SHALLOW DEEP NASAL BRIDGE DEPRESSED PROMINENT NOSE SHORT LONG MIDFACE HYPOPLASTIC HYPERPLASTIC ZYGOMATIC BONE HYPOPLASTIC HYPERPLASTIC SUTURES CLOSE EARLY CLOSE LATE PALATE HIGH LOW PREMAXILLA PROMINENT PROMINENT MOUTH SMALL LARGE TEETH SMALL LARGE LOWER FACE RECEDING PROTRUDING MANDIBULAR RAMUS SHORT SHORT MANDIBLE MICROGNATHIC MACROGNATHIC –4 POINTED AND WIDE AND P = 6.0 × 10 CHIN PROMINENT RECEDING Fig. 6 NFIX down-regulation is associated with modern human-derived traits. a Vocal anatomy of chimpanzee and AMH. The vocal tract is the cavity from the lips to the larynx (marked by dashed lines). In AMHs, the flattening of the face together with the descent of the larynx led to approximately 1:1 proportions of the horizontal and vertical portions of the vocal tract. b Craniofacial features of the Neanderthal, healthy AMH, and AMH with Marshall- Smith or Malan syndromes. Each box shows a phenotype that occurs in the Marshall-Smith/Malan syndromes (i.e., when NFIX is partially or completely inactive). The righthand side of each box shows the observed phenotypes of individuals with the syndromes (left), healthy AMHs (middle) and Neanderthals (right). In most phenotypes, the observed phenotypes match the expected phenotypes based on NFIX expression. Binomial test P-value is shown. is difficult to characterize the full anatomy of the Neanderthal lower and upper projection of the face52. In addition, mutations vocal apparatus, and opinions remain split as to whether it was in NFIX were shown to impair speech capabilities42,53. The exact similar to that of AMHs50,51. mechanism is still unknown, but is thought to occur partly Most skeletal disease phenotypes that result from NFIX dys- through skeletal alterations to the larynx42. To investigate if function are craniofacial, as NFIX influences the balance between changes in NFIX expression could explain morphological changes 10 NATURE COMMUNICATIONS | (2020)11:1189 | https://doi.org/10.1038/s41467-020-15020-6 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15020-6 ARTICLE in the AMH face and larynx, we examined its clinical skeletal patterning20,39–41. It was also shown that facial growth patterns, phenotypes. Mutations in NFIX were shown to cause the which shape facial prognathism, differ between archaic and Marshall–Smith and Malan syndromes, whose phenotypes modern humans not only during early development, but also as include various skeletal alterations such as hypoplasia of the late as adolescence60. Moreover, the main differences between midface, retracted lower jaw, and depressed nasal bridge52. In human and chimpanzee vocal tracts are established during post- many patients, the phenotypic alterations are driven by hetero- infant years61. zygous loss-of-function mutations that cause haploinsufficiency. Although the DMRs we report most likely exist throughout the This shows that reduced activity of NFIX, even if partial, results in skeleton, including the larynx, the evidence we present for the skeletal alterations52. Because NFIX is inferred to have been cranium is more direct, as the patterns are observed in modern downregulated in AMHs compared to archaic humans, we human and chimpanzee crania. Importantly, it has been sug- hypothesized that similar phenotypes to the ones that are driven gested that the 1:1 vocal conformation could have been entirely by NFIX loss-of-function may also exist between modern and driven by cranial, rather than laryngeal, alterations49. Once archaic humans. For example, because reduced activity of NFIX archaic human cranial samples are sequenced, these observations results in a more retracted face, we hypothesized that AMHs could be more directly tested. would present a more retracted face compared to archaic humans. The results we presented open a window to study the evolution A similar relationship between facial features and gene dosage has of the human vocal tract and face from genetic and epigenetic been recently shown in the 16p11.2 locus in humans54. We perspectives. Our data suggest shared genetic mechanisms that therefore examined the phenotypes of the Marshall–Smith and shaped these anatomical regions and point to evolutionary events Malan syndromes and found that not only do most of these that separate AMHs from the Neanderthal and Denisovan. The phenotypes exist between Neanderthals and modern humans, but mechanisms leading to such extensive regulatory shifts, as well as their direction matches the direction expected from NFIX if and to what extent these evolutionary changes affected voca- downregulation along the AMH lineage (18 out of the 22 lization and speech capabilities are still to be determined. Marshall–Smith phenotypes, and 8 out of the 9 Malan pheno- types, P = 6.0 × 10−4, binomial test). In other words, from the Methods Neanderthal, where NFIX activity is expected to be highest, Skeletal methylation maps. Previously, our ability to identify differentially through healthy AMHs, to individuals with NFIX haploinsuffi- methylated regions (DMRs) that discriminate between human groups was confined ciency, phenotypic manifestation matches the level of NFIX by three main factors: (i) We had a single DNA methylation map from a present- day human bone, which was produced using a reduced representation bisulfite activity (Fig. 6b and Supplementary Data 8). sequencing (RRBS) protocol, which provides information for only ~10% of CpG Notably, many cases of laryngeal malformations in the positions in the genome. Moreover, the fact that the archaic and present-day Marshall–Smith syndrome have been reported. Some of the methylomes were produced using different technologies—computational recon- patients exhibit positional changes to the larynx, changes in its struction versus RRBS—potentially introduces a bias. (ii) The analyses included only one bone methylation map from each of the human groups, which limited our width, and, more rarely, structural alterations to the arytenoid ability to identify fixed differences between the groups. Although dozens of maps cartilage—the anchor point of the vocal folds, which controls from additional tissues in present-day humans were included in the analyses, this their movement55. In fact, these laryngeal and facial changes are narrowed the DMRs to represent only human-specific changes that are invariable thought to underlie some of the limited speech capabilities between tissues. (iii) The work did not include a great ape outgroup. Thus, when a AMH-specific change was identified, it was impossible to determine whether it observed in various patients42. This raises the possibility that happened on the AMH lineage, or in the ancestor of Neanderthals and NFIX downregulation in AMHs might be associated with changes Denisovans4. in the larynx too. To overcome these obstacles, a major goal of the current study was to SOX9, ACAN, COL2A1, XYLT1, and NFIX are active in early significantly extend the span of our skeletal methylome collection, covering as many individuals, sexes, and bone types as we could. This included the generation stages of osteochondrogenesis, making the observation of differ- of many samples, including the high-coverage sequencing of additional ancient ential methylation in mature bones puzzling at first glance. This genomes, as listed below. could potentially be explained by two factors: (i) The methylome stabilizes as development progresses and remains largely Present-day human bone DNA methylation maps. We generated full DNA unchanged from late development through adulthood. Thus, methylation maps from two femur head bones from present-day humans using adult methylation states often reflect earlier development, and WGBS. Femora were chosen because of their abundance in present-day human DMRs in adult stages often reflect DMRs in earlier activity samples, as well as in ancient DNA samples5,6,62. In addition, we collected 53 publicly available partial skeletal methylation maps. levels3,17,56. Therefore, these DMRs might reflect early methyla- Trabecular bone tissue from femur heads were taken from two patients with tion changes in mesenchymal progenitors that are carried over to osteoarthritis during a total hip replacement surgery, and after filling in a consent later stages of osteogenesis. Indeed, the methylation patterns of form as per Helsinki approval #0178-13-HMO. Importantly, the effects of NFIX, SOX9, ACAN, and COL2A1 were shown to be established osteoarthritis processes on trabecular bone are much less substantial than those on the synovium, cartilage, and subchondral bone. Bone1 was a left head of femur in early stages of human development and remain stable taken on August 11, 2014 from a 66 years old female and Bone 2 was a right head throughout differentiation from mesenchymal stem cells to of femur taken on September 2, 2014 from a 63-years-old female. mature osteocytes57. It is further supported by the observation DNA was extracted from bones using QIAamp® DNA Investigator kit (56504, that osteoblasts and chondrocytes show almost identical methy- Qiagen). Bones were cut to thin slices (0.2–0.5 mm) and then thoroughly washed lation levels in these DMRs, and are all as hypermethylated as the (X5) with PBS, to clean samples from blood. Bones were crushed with mortar and pestle in liquid nitrogen, and 100 mg bone powder was taken to extract DNA adult bone methylation levels we report58. We have reconfirmed according to the protocol Isolation of Total DNA from Bones and Teeth of the this result by measuring methylation in these DMGs in primary DNA Investigator kit. human chondrocytes. Finally, we show that the upstream We followed the protocol described in ref. 63 to carry out whole-genome mesenchymal enhancer of SOX9 (ref. 23) is differentially methy- bisulfite sequencing at the Center Nacional d’Anàlisi Genòmica (CNAG). In short, DNA libraries were built with the Illumina TruSeq Sample Preparation kit. Bisulfite lated in AMHs (Fig. 4b). (ii) Although expression levels of SOX9, treatment was applied in two rounds using the EpiTect Bisulfite kit (QIAGEN) and ACAN, and COL2A1 gradually decrease with skeletal maturation, paired-end sequencing was performed on an Illumina Hi-Seq 2000 instrument. these genes were shown to remain active in later developmental Reads were aligned using the GEM mapper64, with the reads fully converted in stages in the larynx, vertebrae, limbs, and jaws, including in their silico. To this end, we generated two versions of each of the reference genomes— human (GRCh37) and viral: in the first version, C’s were replaced with T’s, and in osteoblasts21,59. Interestingly, these are also the organs that are the second version, G’s were replaced with A’s. Although methylation state should most affected by mutations in these genes, implying that their late not depend on read position, positional biases have been previously reported65. We stages of activity might still play important roles in morphological observed that the first few bases from each read showed a slightly higher NATURE COMMUNICATIONS | (2020)11:1189 | https://doi.org/10.1038/s41467-020-15020-6 | www.nature.com/naturecommunications 11 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15020-6 probability of being called as methylated, and we thus trimmed the first ten bases were purified with Gel/PCR DNA Mini Kit (RBC, YDF100) and quantified with a from each read (M-bias filtering).. Heterozygous positions, positions with a NanoDrop 2000 spectrophotometer. genotype error probability greater than 0.01, and positions with a read depth CloneJET PCR Cloning Kit (Thermo Scientific, K1231) was used to clone the greater than 250 were filtered out. Only cytosines with six or more reads purified PCR products into a pJET1.2/blunt Cloning Vector following the Blunt- informative for methylation status were considered. On average, half of the reads End Cloning Protocol described in the manufacturer’s instructions. Five microliter from either strand will be informative for methylation status at a given position, so of each cloning reaction product were used for transformation of DH5α Competent minimum coverage is typically greater than 12. Methylated and unmethylated Cells (Invitrogen, 18265017). Colonies were grown overnight on LB plates cytosine conversion rates were determined from spiked-in bacteriophage DNA containing 100 μg/ml ampicillin. Positive transformants were picked and grown (fully methylated phage T7 and unmethylated phage lambda). Five samples were overnight in liquid LB medium containing 100 μg/ml ampicillin. Subsequently, excluded based on conversion rates <0.997, supported by visual inspection of CG plasmid minipreps were purified with a RBC Miniprep Kit (YPD100) according to and non-CG methylation plots. The over-conversion rates for all samples based on manufacturer’s instructions. Purified plasmids were quantified with a NanoDrop methylated phage T7 DNA were ~5%. 2000 spectrophotometer and sequenced on an Applied Biosystems 3730xl Genetic Sequence quality was evaluated using FastQC software v0.11.2. Analyzer (Supplementary Fig. 5a, b). TRIMMOMATIC v.0–32 was used to filter low quality bases with the following parameters: -phred33 LEADING:30 TRAILING:30 MAXINFO:70:0.9 MINLEN:70. Paired-end sequencing reads were mapped to bisulfite converted human (hg19) Human primary chondrocyte validation. Primary chondrocyte cultures were reference genome using Bismark v0.14.3 and bowtie2 v2.2.4 not allowing multiple obtained from osteoarthritis (OA) donors in accordance with Hadassah Medical alignments and using the following parameters: –bowtie2 –non_bs_mm –old_flag Center Institutional Review Board approval and in accordance with the Helsinki -p 4. Potential PCR duplicates were removed using Bismark’s Declaration of ethical principles for medical research involving human subjects. deduplicate_bismark_alignment_output.pl Perl program. Bismark’s End-stage OA patients, with a Kellgren and Lawrence OA severity score of 3–4 bismark_methylation_extractor script was used to produce methylation calls with were recruited following receipt of a formal written informed consent (n = 8; 75% the following parameters: -p –no_overlap –comprehensive –merge_non_CpG – female, mean age 73 ± 7.2 years; mean body mass index 30.1 ± 5.4 kg/m2). Hyaline no_header –bedGraph –multicore 2 –cytosine_report. Examination of the M-bias articular cartilage was dissected and human chondrocytes isolated using 3 mg/ml plots led us to ignore the first 5 bp of both reads in human samples (Supplementary Collagenase Type II (Worthington Cat # LS004177) in DMEM medium (Sigma- Fig. 1). Custom scripts were used to summarize methylation levels at CpG sites Aldrich, St Louis, MI) containing 10% FCS and 1% penicillin-streptomycin (Beit- based on the frequencies of methylated and unmethylated mapped reads on both Haemek Kibutz, Israel), 37 °C, 24 h incubation. Isolated cells were filtered through strands. Methylation data were deposited in NCBI’s Gene Expression Omnibus and a nylon cell strainer (40 mm diameter), washed three times with PBS and plated at are accessible through GEO accession number GSE96833. 1.5 million cells per 14 cm2 tissue culture dish (passage 0, passage 2). Cells were cultured in standard incubation conditions (37 °C, 5% CO2) until confluence. Chondrocyte DNA purification was performed using GenElute™ Mammalian Partial skeletal and full non-skeletal DNA methylation maps. Osteoblast RRBS Genomic DNA Miniprep Kit (Sigma, G1N350). map, extracted from the femur, tibia, and rib bones of a 6-year-old female (NHOst- Osteoblasts by Lonza Pharma, product code: CC-2538, lot number: 6F4124), was Chimpanzee bone DNA methylation maps. Overall, we produced six methylation downloaded from GEO accession number GSE27584. Forty-eight 450K methyla- maps from bones of six common chimpanzee (Pan troglodytes) individuals. They tion array maps, extracted from the femora of adult males and females with include one WGBS of a wild chimpanzee, one RRBS of an infant chimpanzee, and osteoarthritis or osteoporosis, were downloaded from GEO accession number four 850K methylation arrays of captive chimpanzees. GSE64490. Four 450K methylation array maps, extracted from unspecified bones of Chimpanzee tissue samples included in this study were opportunistically adult males and females were downloaded from GEO accession number GSE50192. collected at routine necropsy of these animals. No animals were sacrificed for this Chimpanzee and human WGBS blood methylation maps were downloaded from study, and no living animals were used in this study. NCBI SRA accession number SRP059313. Chimpanzee and human WGBS brain maps were downloaded from GEO accession number GSE37202. WGBS of a chimpanzee bone. We used a rib bone of a 47-year-old female Chimpanzee provided from the Biobank of the Biomedical Primate Bisulfite-PCR of human bone. A skull of an adult male from India was obtained Research Center (BPRC), The Netherlands. The postmortem interval was from the teaching anatomy collection of the Department of Anatomy and approximately 10–12 h. The bone was collected during the necropsy procedure and Anthropology at the Sackler Faculty of Medicine, Tel Aviv University, Israel immediately frozen and stored at –80 °C. (Human 1). Additional two skull specimens (Human 2 and 3) were obtained DNA was extracted in a dedicated ancient DNA laboratory at the Institute of directly from the operating room of the Department of Neurosurgery, Shaare Evolutionary Biology in Barcelona, where no previous work on great apes has ever Zedek Medical Center, Jerusalem, Israel and transferred on dry ice for further been conducted. Standard precautions to avoid and monitor exogenous analysis. All study participants provided informed consent according to an insti- contamination such as frequent cleaning of bench surfaces with bleach, use of tutional review board—approved protocol (SZMC 0048-18). sterile coveralls, UV irradiation and blank controls were taken during the process. Human 1: Standard precautions to avoid contamination were taken, including Two hundred milligram of bone powder were obtained by drilling and the sample wearing disposable coats, masks, hair covers and double gloves. All following steps was extracted following the Dabney et al.66 method. A final 25 µl of extract volume were performed in a UV cabinet dedicated for the preparation of ancient bone was used for genome sequencing. samples and located in a physically separated ancient DNA laboratory at the Analysis was performed similarly to Bone1 and Bone2, with the exception that Faculty of Dental Medicine. The skull was cleaned with an excess of 10% bleach the BSreads were mapped to bisulfite converted chimpanzee (panTro4) reference (equal to 0.6% Sodium hypochlorite) and then subjected to UV radiation for genome, and we ignored the first 5bp of read1 and the first 44 bp of read2 in the 30 min The cortical layer on the temporal surface (facies temporalis) of the chimpanzee sample (Supplementary Fig. 1). Methylation data were deposited in zygomatic bone (ZB) was removed by low-speed drilling using a Wolf Multitool NCBI’s Gene Expression Omnibus and are accessible through GEO accession Combitool Rotary Multi Purpose Tool equipped with a sterile dental burr. Another number GSE96833. sterile burr was used to obtain powder of the subcortical trabecular bone within the body of the zygoma. The powder was collected onto a 10 × 10 cm aluminum foil sheet pretreated with a 10% bleach solution and then transferred into a sterile RRBS of a chimpanzee bone. We used two unidentified long bone fragments that 1.5 ml Eppendorf tube for subsequent DNA extraction. Altogether, three samples belonged to a newborn wild chimpanzee infant who died during a documented were obtained: ZB-3 from the right zygoma weighing 20.3 mg, and ZB-3/1 and ZB- infanticide event at Gombe National Park on 9 March 2012. The infant was known 3/2 from the left zygoma weighing 29.5 mg and 30.3 mg, respectively. Bone DNA to be the offspring of a chimpanzee called Eliza and was partially eaten by an adult was purified from the three bone powder samples using QIAamp DNA Investigator female and her family. The sample was collected from the ground about 48 h after Kit (QIAgen, 56504) according to manufacturer’s instructions. the infant's death and stored in RNAlater solution until arrival at Arizona State Human 2 and 3: DNA was extracted from bones using QIAamp® DNA University (ASU). At ASU the sample was stored at 4 °C until extraction. Investigator kit (56504, Qiagen). Bones were thoroughly washed (X5) with PBS, to Sampling and DNA extractions were conducted at the ASU Ancient DNA clean samples from blood. Bones were crushed with mortar and pestle in liquid Laboratory, a Class 10,000 clean-room facility in a separate building from the nitrogen, and 100 mg bone powder was taken to extract DNA according to the Molecular Anthropology Laboratory. Precautions taken to avoid contamination protocol Isolation of Total DNA from Bones and Teeth of the DNA included bleach decontamination and UV irradiation of tools and work area before Investigator kit. and between uses, and use of full body coverings for all researchers. The bone Genomic DNA was bisulfite converted with the EZ DNA Methylation Lightning samples were pulverized together in December 2012 using a SPEX CertiPrep Kit (Zymo Research, D5030) according to the manufacturer’s instructions. Freezer Mill. Three DNA extractions were conducted using 50–100 mg of bone Specifically, each bone sample was bisulfite converted using 500 ng as genomic powder (Supplementary Table 2) and following the extraction protocol by Rohland DNA input for the conversion. and Hofreiter67. Two extraction blank controls were included to monitor Bisulfite treated DNA were amplified with the FastStart High Fidelity PCR contamination of the extraction process. One microliter each of the sample extract System (Sigma, 03553400001) using the primers listed in Supplementary Table 1. and the blank control were used for fluorometric quantification with the Qubit 2.0 PCR conditions were performed according to manufacturer’s instructions and PCR Broad Range assay. All extracts were combined for a total volume of 345 µl and products were visualized on a 1.5% agarose gel. Prior to cloning, PCR products approximately 0.652 µg of total DNA. 12 NATURE COMMUNICATIONS | (2020)11:1189 | https://doi.org/10.1038/s41467-020-15020-6 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15020-6 ARTICLE RRBS libraries were generated according to Boyle et al.68. 100–200 ng genomic sites within the sex chromosomes were removed using the minfi package in R. This DNA was digested with MspI. Subsequently, the digested DNA fragments were filtering retained a final set of 576,505 probes. end-repaired and adenylated in the same reaction. After ligation with methylated adapters, samples with different adapters were pooled together and were subjected to bisulfite conversion using the EpiTect Bisulfite kit (QIAGen) per the Bisulfite-PCR of chimpanzee cranial bones. Postmortem frontal skull bones from manufacturer’s recommendations with the following modification: after first two different chimpanzees (chimpanzee 1 and chimpanzee 2) were provided by the bisulfite conversion, the converted DNA was treated with sodium bisulfite again to Biomedical Primate Research Center (BPRC, The Netherlands). Bones were guarantee that conversion rates were no less than 99%. Two third of bisulfite opportunistically collected during routine necropsy of these animals and stored at converted DNA was PCR amplified and final RRBS libraries were sequenced in an −80 °C. Chimpanzee 3 and chimpanzee 4 samples were obtained from the chim- Illumina HiSeq 2000 sequencer (Supplementary Data 1). Methylation data were panzee cranial collection in the Department of Paleoanthropology in the Senck- deposited in NCBI’s Gene Expression Omnibus and are accessible through GEO enberg Research Institute Frankfurt (DPSF) and Natural History Museum accession number GSE96833. Frankfurt. These two chimpanzee specimens are owned by the Justus Liebig University Gießen. Chimpanzee 1 and chimpanzee 2: For each sample, bone powder was obtained 850K DNA methylation arrays. Four chimpanzee cadavers from captive colonies by crushing the bones with mortar and pestle. Approximately 100 mg bone powder at the Southwest National Primate Research Center in Texas were used. Femora were used to extract DNA using the QIAamp DNA Investigator Kit (Qiagen) were opportunistically collected at routine necropsy of these animals and stored in following manufacturer’s instructions. −20 °C freezers at the Texas Biomedical Research Institute after dissection. These Chimpanzee 3 and chimpanzee 4: Cochlear bone powder was obtained by preparation and storage conditions ensured the preservation of skeletal DNA accessing the petrous bone from the cranial base72. DNA was extracted from about methylation patterns. 50 mg of powder according to the protocol described by Dabney et al.66, but Samples were then transported to ASU and DNA was extracted from the adapted for the use of High Pure Nucleic Acid Large Volume columns (Roche) femoral trabecular bone using a phenol-chloroform protocol optimized for skeletal instead of the Zymo-Spin V column (Zymo Research) MinElute silica spin column tissues69. From the distal femoral condyles, trabecular bone was collected using (Qiagen) combination. coring devices and pulverized into bone dust using a SPEX SamplePrep Freezer/ Genomic DNA was bisulfite converted with the EZ DNA Methylation— Mill. Specifically, bone cores were obtained from a transverse plane through the Lightning Kit (Zymo Research, D5030) according to the manufacturer’s center of the medial condyle on the right distal femur, such that the articular instructions. Specifically, each bone sample was bisulfite converted two times in surface remained preserved. Cortical bone was removed from these cores using a parallel using 500ng as genomic DNA input for the conversion. Dremel (Supplementary Table 3). Tissue collections were performed at the Texas Three microliter of bisulfite treated DNA were amplified with the FastStart Biomedical Research Institute, and DNA extractions were conducted at the ASU High Fidelity PCR System (Sigma, 03553400001) using the primers listed in Molecular Anthropology Laboratory. Supplementary Table 1. PCR conditions were performed according to Genome-wide DNA methylation was assessed using Illumina Infinium manufacturer’s instructions and PCR products were visualized on a 1.5% agarose MethylationEPIC microarrays. These arrays analyze the methylation status of over gel. Prior to cloning, PCR products were purified with homemade SPRI beads 850,000 sites throughout the genome, covering over 90% of the sites on the (chimpanzee 1 and 2) and Gel/PCR DNA Mini Kit (RBC, YDF100, chimpanzee 3 Infinium HumanMethylation450 BeadChip as well as an additional 350,000 sites and 4), and quantified with a NanoDrop 2000 spectrophotometer. within enhancer regions. For each sample, 400 ng of genomic DNA was bisulfite CloneJET PCR Cloning Kit (Thermo Scientific, K1231) was used to clone the converted using the EZ DNA MethylationTM Gold Kit according to the purified PCR products into a pJET1.2/blunt Cloning Vector following the Blunt- manufacturer’s instructions (Zymo Research), with modifications described in the End Cloning Protocol described in the manufacturer’s instructions. Three Infinium Methylation Assay Protocol. These protocols were conducted at the ASU microliter (chimpanzee 1 and 2) and 3 µl (chimpanzee 3 and 4) of each cloning Molecular Anthropology Laboratory. Following manufacturer guidelines reaction product were used for transformation of DH5α Competent Cells (Illumina), this processed DNA was then whole-genome amplified, enzymatically (Invitrogen, 18265017). Colonies were grown overnight on LB plates containing fragmented, hybridized to the arrays, and imaged using the Illumina iScan system. 100 μg/ml ampicillin. Positive transformants were picked and grown overnight in These protocols were conducted at the Texas Biomedical Research Institute. These liquid LB medium containing 100 μg/ml ampicillin. Subsequently, plasmid array data have been deposited in NCBI’s Gene Expression Omnibus and are minipreps were purified with a QIAprep Miniprep Kit (Qiagen, chimpanzee 1 and accessible through GEO Series accession number GSE94677. 2), and RBC Miniprep Kit (YPD100, chimpanzee 3 and 4) according to Raw fluorescent data were normalized to account for the noise inherent within manufacturer’s instructions. Purified plasmids were quantified with a NanoDrop and between the arrays themselves. Specifically, we performed a normal- 2000 spectrophotometer and sequenced on an Applied Biosystems 3730xl Genetic exponential out-of-band (Noob) background correction method with dye-bias Analyzer (Supplementary Fig. 5a, b). normalization to adjust for background fluorescence and dye-based biases and followed this with a between-array normalization method (functional normalization) which removes unwanted variation by regressing out variability Reconstructing ancient DNA methylation maps. In a dedicated clean room at explained by the control probes present on the array as implemented in the minfi Harvard Medical School, powder was extracted from the root of a lower third package in R which is part of the Bioconductor project. This method has been molar of the Mesolithic La Braña 1 individual (5983–5747 calBCE (6980 ± 50 BP, found to outperform other existing approaches for studies that compare conditions Beta-226472)), from which a non-UDG-treated library was previously sequenced with known large-scale differences70, such as those assessed in this study. to 3.5× coverage8. Two UDG-treated libraries from the same individual were later After normalization, methylation values (β values) for each site were calculated generated and enriched for approximately 1.2 million single targeted polymorph- as the ratio of methylated probe signal intensity to the sum of both methylated and isms and sequenced to an average of 19.5× coverage at these position9. In this unmethylated probe signal intensities. These β values range from 0 to 1 and study, we carried out shotgun sequencing of one of the same UDG-treated libraries represent the average methylation levels at each site across the entire population of from this individual on a NextSeq500 instrument using 2 × 76 bp paired end cells from which DNA was extracted (0 = completely unmethylated sites, 1 = fully sequences. Following the mapping protocol described previously9, we trimmed methylated sites). adapter sequences, only processed read pairs whose ends overlapped by at least 15 bp (allowing for one mismatch) so that we could confidently merge them, and then mapped to the human reference sequence hg19 using the command samse in Methylated signal β ¼ BWA (v0.6.1). We removed duplicated sequences by identifying sequences with the ðMethylated signal þ unmethylated signalÞ same start and stop position and orientation in the alignment, and picking the highest quality one. After restricting to sequences with a map quality of MAPQ ≥ Every β value in the Infinium platform is accompanied by a detection P-value, 10, and sites with a minimum sequencing quality (≥20), we had an average cov- and those with failed detection levels (P-value > 0.05) in greater than 10% of erage measured at the same set of approximately 1.2 million single nucleotide samples were removed from downstream analyses. polymorphism targets of 23.0×. This data is available under GEO accession The probes on the arrays were designed to specifically hybridize with human number: GSE96833, with raw reads deposited under SRA accession number: DNA, so our use of chimpanzee DNA required that probes non-specific to the SRX3194436. chimpanzee genome, which could produce biased methylation measurements, be In a dedicated clean room at the University College Dublin, powder was computationally filtered out and excluded from downstream analyses. This was extracted from the cochlear portion of the petrous bone of individual I1583 accomplished using methods modified from Hernando-Herraez et al.71. Briefly, we (archeological ID L14-200) from the site of Barcın Höyük in the Yenişehir Plain of used blastn to map the 866,837 50 bp probes onto the chimpanzee genome the Marmara Region of Northwest Turkey. The Neolithic individual came from a (Assembly: Pan_tro_3.0, Accession: GCF_000001515.7) using an e-value threshold community that practiced farming, and was anthropologically determined to be of e−10. We only retained probes that successfully mapped to the genome, had only a male aged 6–10 years at the time of death (the sex was confirmed genetically). 1 unique BLAST hit, targeted CpG sites, had 0 mismatches in 5 bp closest to and The direct radiocarbon date was 6426-6236 calBCE (7460 ± 50 BP, Poz-82231). In a including the CpG site, and had 0–2 mismatches in 45 bp not including the CpG dedicated clean room at Harvard Medical School, a UDG-treated library was site. This filtering retained 622,819 probes. prepared from this powder, which was previously enriched for about 1.2 million Additionally, β values associated with cross-reactive probes, probes containing SNP targets, sequenced to 13.5× average coverage, and published in ref. 9. We SNPs at the CpG site (either human or chimpanzee), probes detecting SNP shotgun sequenced the same library on nine lanes of a HiSeqX10 sequencing with information, probes detecting methylation at non-CpG sites, and probes targeting 100bp paired reads. On data processing, we merged overlapping read pairs, NATURE COMMUNICATIONS | (2020)11:1189 | https://doi.org/10.1038/s41467-020-15020-6 | www.nature.com/naturecommunications 13 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15020-6 trimmed Illumina sequencing adapters, and dropped read pairs that did not have finally compared them to the original map and identified DMRs. Any differences in sample barcodes (up to one mismatch) or cannot be unambiguously merged. We methylation levels between the simulated map and the original reference map stem then aligned merged reads with BWA against human reference genome GRCh37 from noise. Thus, running the same DMR-detection algorithm on the simulated (hg19) plus decoy sequences, and combined all nine lanes of data and removed map vs. the reference map, enables an estimation of the false discovery rate. We set duplicate molecules, achieving an average of 24.3× coverage evaluated on the 1.2 the DMR-detection thresholds so that FDR < 0.05. million targets. This data is available under GEO accession number: GSE96833, Step 4: Lineage assignment. The chimpanzee methylation maps were used to with raw reads deposited under SRA accession number: SRX3194436. polarize the DMRs. For each DMR, methylation levels in the chimpanzee were compared to those of the three hominin groups. For example, if methylation levels in the chimpanzee samples clustered with the archaic humans, the DMR was The reconstruction procedure. Reconstruction of DNA methylation maps was assigned to the AMH lineage. performed on the genomes of the following individuals: Ust'-Ishim6, Loschbour7, Step 5: Within-lineage variability filtering. To determine whether a DMR Stuttgart7, La Braña 1, I1583, and the Vindija Neanderthal5, as well as on the represents an individual within a group, or is shared by the entire group, we used a previously published Altai Neanderthal and the Denisovan (Supplementary total of 67 AMH, archaic and chimpanzee methylation maps. We used a Data 1). The Vindija Neanderthal reads were downloaded from the Max Planck conservative approach, where DMRs in which methylation levels in one group Institute for Evolutionary Anthropology website: http://cdna.eva.mpg.de/ overlap (even partially) the methylation levels in another group were discarded. As neandertal/Vindija/bam/. Only the UDG-treated portion of the genome (B8744) 59 out of the 67 maps belong to AMHs, our ability to filter out variation within this was used. Additional UDG-treated ancient human full genomes have been pub- group was better, resulting in fewer DMRs along this lineage. Several various lished to date; however, these were sequenced to a relatively low coverage (<5×), measures were used to ascertain that a DMR along a lineage does not represent a and thus, only crude methylation maps could be reconstructed from them. C→T sex-specific, bone-specific, age-specific, technology, or disease-specific DMR. ratio was computed for every CpG position along the hg19 (GRCh37) human genome assembly, for each of the samples4. In order to exclude from the analyses positions that potentially represent pre- DMR-detection algorithm. We developed an algorithm specifically designed to mortem C→T mutations rather than post-mortem deamination, the following identify DMRs between a deamination map and a full methylome reference. Let i filters were applied: (i) Positions where the sum of A and G reads was greater than enumerate the CpG positions in the genome. In the deamination map, let ti be the the sum of C and T reads were excluded. (ii) For genomes that were produced number of T’s at the C position + the number of A’s in the opposite strand at the G using single-stranded libraries (i.e., Ust'-Ishim, Altai Neanderthal, Denisovan, position, i.e., it counts the total number of T’s that appear in a position that is Vindija Neanderthal and ~1/3 of the Loschbour library) positions, where the G → originally C, in the context of a CpG dinucleotide. We similarly use ci to count the A ratio on the opposite strand was greater than 1/(average single strand coverage) total number of C’s that appear in a position that is originally C, in the context of a were excluded. This fraction represents a threshold of one sequencing error allowed CpG dinucleotide. The C → T ratio is defined as ti/ni, where ni = ci + ti. Let φi and per position. For Loschbour, this was performed only on the fraction of reads that ψi (both between zero and one) be the methylation of this position in the reference came from the single stranded library. (iii) For all genomes, positions with a C → T genome and in the reconstructed one, respectively. If we denote by π the deami- ratio >0.25 were discarded. For the Vindija Neanderthal, this threshold was raised nation rate, assumed to be constant throughout the genome, and if we assume that to 0.5, due to its relatively low coverage (~7×). (iv) Finally, a maximum coverage deamination of C into T is a binomial process with probability of success πψi, threshold of 100 reads was used to filter out regions that are suspected to be PCR we get duplicates.   ti  B ni ; πψ i  ð1Þ In all genomes, excluding Vindija, a fixed sliding window of 25 CpGs was used for smoothing of the C → T ratio. This allowed for an unbiased scanning of Our null hypothesis is that the ith CpG is not part of a DMR, namely that ψi = differentially methylated regions (DMRs) that is not affected by the size of the φi. The alternative hypothesis states that this CpG is part of a DMR. The definition window. Due to its relatively low coverage, we extended the sliding window used of this statement is that |ψi − φi| ≥ Δ, where Δ is some pre-specified threshold. In on the Vindija genome to 50 CpGs. This extended window is not expected to other words, under the alternative hypothesis we get that ψi ≥ φi + Δ if the site has introduce a bias, as this genome was not used for DMR detection, but only for low methylation in the reference genome, and ψi ≤ φi − Δ if it has high methylation subsequent filtering that was applied equally to all genomes (see later). in the reference genome. C → T ratio was translated to methylation percentage using linear transformation determined from two points: zero C → T ratio was set to the value Per-site statistic. Let us start with the first option, testing whether ψi ≥ φi + Δ 0% methylation, and mean C → T ratio in completely methylated (100% when φi is low. A log-likelihood-ratio statistic would be methylation) CpG positions in modern human bone reference (hereinafter μ100)         was set to the value 100% methylation. Positions where C → T ratio > μ100 were set Pr ti jni ; π φi þ Δ Δ 1  π φi þ Δ ‘þ i ¼ ln   ¼ ti ln 1 þ  ln to 100% methylation. For genomes that were extracted from bones, the modern Pr ti jni ; πϕi φi 1  πφi Bone 2 WGBS map, which is the one with the higher coverage between the two   1  π φi þ Δ WGBS modern bone maps, was used to determine μ100. For genomes that were þ ni ln  extracted from teeth, there was no available modern reference methylation map, 1  πφi and therefore, we transformed the C → T ratio into methylation percentage based Similarly, we can test whether ψi ≤ φi − Δ when φi is high using the log- on the assumption that the genome-wide mean methylation is similar to bone likelihood-ratio statistic tissue. Thus, the genome-wide mean C → T ratio represents 75% methylation,         which is the genome-wide mean of measured methylation in the Bone 2 reference Pr ti jni ; π φi  Δ Δ 1  π φi  Δ ‘ ¼ ln   ¼ t ln 1   ln map. This was accomplished by setting μ100 to 1.33× mean genome-wide C → i Pr ti jni ; πϕi i φi 1  πφi T ratio.   1  π φi  Δ þ ni ln : 1  πφi DMR detection. The DMR detection algorithm is comprised of five main steps. We used the value Δ = 0.5 for all samples. The value of π, the deamination rate, We hereby provide an overview of the algorithm followed by a detailed description was estimated using the overall C → T ratio in CpG positions whose methylation of each step. The overall goal of this pipeline is to detect differential methylation, level is 1 in the modern human Bone 2 WGBS methylation map, after exclusion of assign it to the lineage on which it arose and filter out within-lineage variation. putative pre-mortem substitutions, as described in the reconstruction procedure Step 1: Two-way comparisons. To avoid artifacts that could potentially be section (Supplementary Data 1). introduced by comparing DNA methylation maps that were produced using different technologies, our core analysis relied on the comparison of the three reconstructed maps of the Altai Neanderthal, Denisovan, and Ust’-Ishim. Each of Detecting DMRs. The statistics ‘þ  i and ‘i quantify how strongly the estimated the samples was compared to the other two in a pair-wise manner, as a raw C → T methylation in position i deviates from φi. Next, we use these values to identify ratio map against a reconstructed methylation map, and vice versa. This reciprocal DMRs using the cumulative-sum procedure explained below. The process is comparison insured that the reconstruction process does not introduce biases to repeated twice: on the statistic ‘þi to identify DMRs where the sample has elevated one of the groups. The minimum methylation difference threshold was set to 50%, methylation with respect to the reference, and on the statistic ‘ i to identify DMRs spanning >50 CpGs. where the sample has reduced methylation with respect to the reference. Step 2: Three-way comparisons. This step classifies to which of the three For convenience, we explain the cumulative-sum procedure in the context of ‘þ  + hominins the DMR should be attributed. This step is done by overlapping the three i , but an essentially identical procedure is used for ‘i . We define a new vector Q lists of DMRs found in Step 1. For example, a DMR that is detected between the by the recursion Neanderthal and Ust’-Ishim and also between the Denisovan and Ust’-Ishim is Qþ 0 ¼ 0; Qþ þ þ i ¼ maxðQi1 þ ‘i ; 0Þ: considered specific to Ust’-Ishim. Step 3: FDR filtering. Various factors could introduce noise to the Under the null hypothesis, ‘þ ihas a negative expectation which produces a reconstruction process, including the stochasticity of the deamination process, the negative drift that keeps Q+ at zero, or close to zero, levels. Under the alternative use of a sliding window, and variations in read depth within a sample. We ran hypothesis the expectation is positive, hence the drift over a DMR is positive, simulations that mimic the post-mortem degradation processes of ancient DNA, leading to an elevation in the values of Q+. Therefore, our next step is to find all then reconstructed methylation maps from the simulated deamination maps and intervals [a, b] such that Qþ þ þ þ a1 ¼ 0, Qbþ1 ¼ 0, and Qi > 0 for a ≤ i ≤ b. Let Qm be 14 NATURE COMMUNICATIONS | (2020)11:1189 | https://doi.org/10.1038/s41467-020-15020-6 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15020-6 ARTICLE the maximum value of Q+ in this interval, where m is the position of the reference map stem from noise. Thus, running the same DMR-detection algorithm maximum. Then, the interval [a, m] would be called a putative DMR. described above on the simulated map vs. the reference map, enables an estimation The statistics ‘þ  i and ‘i are affected by the number of observed cytosine reads, of the false discovery rate. We ran these simulations 100 times for each of the three and thus have higher power to detect hypermethylation (i.e., larger number of genomes (Altai Neanderthal, Denisovan, Ust’-Ishim) and determined the values of cytosine reads) vs. hypomethylation (Supplementary Fig. 3). the Qþ  T and QT thresholds (see the Filtering DMRs section above) such that the mean number of DMRs that are detected in the simulations is < 0.05 the number of real DMRs detected (i.e., FDR < 0.05). Filtering DMRs. Of course, Q+ may increase locally due to randomness, and thus a putative DMR may not reflect a true DMR. To filter out such intervals, we used two strategies. First, we applied a set of filters to assure that the putative DMRs have DMRs separating chimpanzees and humans. To identify DMRs that separate reasonable biological properties. Second, we cleaned the remaining putative DMRs chimpanzees from all human groups (both modern and archaic), we first compared by applying a false discovery rate (FDR) procedure. In the first strategy, we applied the chimpanzee WGBS bone methylome to each of two present-day WGBS maps two filters: (i) Putative DMRs that harbor less than 50 CpG positions, thus are (those of Bone1 and Bone2). This was done by scanning the chimpanzee map using shorter than twice the smoothing window size, were removed. (ii) To avoid a sliding window of 25 CpGs, in intervals of one CpG position. In each window, we situations where two consecutive CpG sites whose genomic locations are remote counted the number of methylated and unmethylated reads in each sample, and appear on the same DMR, we modify the vector Qþ i as follows. Let di,j be the computed a P-value using Fisher's Exact test. We then computed FDR-adjusted P- distance along the genome (in nucleotides) between CpG sites i and j. Then, for values for each window, and discarded windows with FDR > 0.05 or where the every site i such that di,i–1 > δ we set Qþ i ¼ 0. We used δ = 1000 nt for all samples. mean methylation difference (Δ) was below 0.5. We then merged overlapping To further remove putative DMRs that are unlikely to reflect true DMRs, we windows. This left 8040 DMRs between the chimpanzee and Bone1, and 12,666 eliminated all DMRs where Qþ þ þ + DMRs between the chimpanzee and Bone2. Next, we intersected the two lists to m 2 standard deviations from the mean observed ability to estimate variability within archaic human lineages, we added to the methylation in modern humans. This suggests that although the Neanderthal is analysis the reconstructed methylation map of the Vindija Neanderthal. The hypermethylated compared to almost all modern humans, she is not found USER-treated portion of this genome (the portion amenable for methylation completely outside modern human variation. The Denisovan, on the other reconstruction) was sequenced5 to a depth of 7×. Therefore, the methylation map hand, is found even further away, and significantly outside the other popula- that could be reconstructed from this individual has a considerably lower tions. Given this, the HOXD DMR was classified as Denisovan-derived (Sup- resolution compared to the other reconstructed maps used in this study (coverage plementary Data 2). The Ust'-Ishim remains include a single femur, and to our 19× to 52×). Nevertheless, due to the reduced ability to detect variability along the knowledge, it was not compared morphologically to other humans. Thus, archaic human linages, we employed this map for additional variability filtering further analysis is needed in order to determine whether the hypermethylation along these lineages. DMRs where the Vindija Neanderthal clustered with the other of the Ust'-Ishim individual compared to other AMHs is manifested in mor- hominins, and not with the Altai Neanderthal (or not with either of the archaic phological changes as well. Moreover, as this DMR is classified as Denisovan- humans in the archaic-derived DMRs) were discarded. The number of DMRs derived, we cannot rule out the possibility that it is driven to some extent by mentioned throughout this chapter already includes this filtering. age or bone type differences. A general concern in working with DNA methylation data is that DMRs that Compared to the previously reported DMRs4, in this study we found four are specific to one group do not necessarily represent an evolutionary change, but times as many AMH-derived and archaic-derived DMRs (2805 full bone DMRs rather reflect a characteristic such as technology used to measure methylation, compared to 891) and roughly twice as many Neanderthal- and Denisovan- tissue, sex, disease or age that is shared by individuals in this group and not by derived DMRs (440 and 598 compared to 295 and 307 in the Denisovan and others. We take two complementary approaches to ascertain that the DMRs we Neanderthal, respectively). The list of DMRs reported here cannot be directly report are not driven by these factors: (a) for the top DMGs, we match the samples compared to our previous list of DMRs because of several key differences in the for the above factors and test whether the hypermethylation of AMHs is still analysis: (i) The previous study focused on DMRs that are invariable across observed. To this end, we compared Ust’-Ishim (adult femur with no known tissues, whereas here we focused on DMRs in skeletal tissues. In the previous diseases, methylation map produced using our reconstruction method) to the study, we were therefore able to extrapolate and find trends that extend beyond Vindija Neanderthal (adult femur with no known diseases, methylation map the skeletal system, such as neurological diseases. In this paper we focus on the produced using our reconstruction method), and we also compared 52 modern skeletal system, hence the different appearance of the body map (Fig. 2b, c). (ii) human samples (adult femora, methylation array maps) to four chimpanzee the current study used stricter thresholds for DMR detection, including a samples (adult femora, methylation array maps). In all cases, AMHs show minimum of 50 CpGs in each DMR (compared to ten CpGs previously), and a significant hypermethylation compared to the matched samples (Supplementary requirement for physical overlap in the three-way DMR detection procedure. Fig. 3c, d, see the Methylation in AMH, chimpanzee and Neanderthal femora (iii) In this study, the AMH reference is a reconstructed ancient map, whereas in chapter for additional information). (b) throughout the pipeline, we take only 16 NATURE COMMUNICATIONS | (2020)11:1189 | https://doi.org/10.1038/s41467-020-15020-6 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15020-6 ARTICLE the previous study the AMH reference, as well as the other tissues used for analyses are based on skeletal tissues, we tested whether the proportion of genes filtering out noise, were mainly cultured cell lines with RRBS methylation maps. related to the face, larynx, vocal folds, and pelvis within AMH-derived skeleton- When filtering DMRs along the lines of the previous study by taking only related DMGs is higher than expected by chance. Out of 100 skeleton-related DMRs with low inter-tissue variability in humans (STD < 10%), we indeed observe DMGs, 31 genes are known to affect the voice, 34 affect the larynx, 87 affect the similar trends. For example, when taking AMH-derived DMGs and analyzing their face, and 65 affect the pelvis, whereas genome-wide these proportions are sig- expression patterns using DAVID's tissue expression tools, we found that the brain nificantly lower (14.2%, 20.2%, 70.0%, 52.4%, P = 1.0 × 10−5, P = 1.3 × 10−3, P = is the most represented organ, with 51.5% of DMGs expressed in this organ (×1.28, 2.1 × 10−3, P = 0.03, for vocal folds, larynx, face, and pelvis, respectively, hyper- FDR = 2.6 × 10−4), and glial cells are the most over-represented cell type (×20.6, geometric test). For additional validation tests, see main text. but FDR > 0.05, UP_TISSUE DB, Supplementary Data 3). In fact, the brain is the Genes associated with craniofacial features were taken from the GWAS-catalog only significantly enriched organ in this analysis. Similarly, when analyzing the (version 2019-04-21), using a threshold of P < 10−8. The following features were GNF DB, we found that the subthalamic nucleus is the most enriched body part used: dental caries, cleft palate, facial morphology, intracranial volume, cleft palate (×1.60, FDR = 9.2 × 10−4), followed by additional brain regions, such as the (environmental tobacco smoke interaction), cranial base width, craniofacial olfactory bulb (×1.54, FDR = 0.01), globus pallidus (×1.41, FDR = 0.04), and more macrosomia, facial morphology (factor 1, breadth of lateral portion of upper face), (Supplementary Data 3). Similar enrichment patterns of the brain can be observed facial morphology (factor 10, width of nasal floor), facial morphology (factor 11, when analyzing expression patterns of all AMH-derived DMGs (Supplementary projection of the nose), facial morphology (factor 12, vertical position of sublabial Data 3). Finally, we also find that similarly to the previous report, these DMRs are sulcus relative to central midface), facial morphology (factor 14, intercanthal linked to diseases more often (23.1% compared to the genome average of 10.8%, width), lower facial height, nose morphology, nose size, tooth agenesis (maxillary DAVID OMIM_DISEASE DB). third molar), tooth agenesis (third molar), facial morphology traits (multivariate analysis), lower facial morphology traits (ordinal measurement), lower facial morphology traits (quantitative measurement), middle facial morphology traits Computing correlation between methylation and expression. In order to (quantitative measurement), and upper facial morphology traits (ordinal identify regions where DNA methylation is tightly linked with expression levels, we measurement). We then tested their overlap with DMGs. Genes associated with scanned each DMR in overlapping windows of 25 CpGs (the window used for craniofacial features in the GWAS catalog significantly overlapped DMGs smoothing the deamination signal). In each window we computed Pearson’s cor- compared to the fraction expected by chance (5.17×, P = 3.4 × 10−4, relation between DNA methylation levels and expression levels of overlapping hypergeometric test). As a control, we then tested how this 5-fold enrichment genes as well as the closest genes upstream and downstream genes, across 21 compares to non-craniofacial features. We used blood-related GWAS as a tissues33. For each DMR, we picked the window with the best correlation (in representative of general non-craniofacial GWAS. We extracted from the GWAS absolute value) and computed regression FDR-adjusted P-value. DMRs that catalog 22 blood-related traits (the same number as extracted for craniofacial overlap windows with FDR < 0.05 were considered to be regions where methylation features), by taking the first 22 traits that appear in a search for the term blood and levels are significantly correlated with expression levels. Ninety such DMRs were applying a threshold of P < 10−8. We then used these genes as a background found among the skeletal AMH-derived DMRs, 93 among the archaic human- control for the craniofacial enrichment. We observed a 3.86× enrichment of DMGs derived DMRs, 40 among Neanderthal-derived DMRs, and 19 among Denisovan- with regard to craniofacial-associated vs. non-craniofacial-associated genes (P = derived DMRs. 0.01, chi-square test). As no expression data were available for Ust'-Ishim, Bone1 and Bone2, we Additionally, we conducted a permutation test on the list of 129 AMH-derived approximated their NFIX expression level by taking the average of NFIX expression DMGs that are linked to organs on Gene ORGANizer, replacing those that are from three osteoblast RNA-seq datasets that were downloaded from GEO accession linked to the skeleton with randomly selected skeleton-related genes. We then ran numbers GSE55282, GSE85761, and GSE78608. RNA-seq data for chondrocytes the list in Gene ORGANizer and computed the enrichment. We repeated the was downloaded from the ENCODE project, GEO accession number GSE78607 process 100,000 times and found that the enrichment levels we observed within and plotted against measured methylation levels in primary chondrocytes (see the AMH-derived DMGs are significantly higher than expected by chance for the Human primary chondrocyte validation chapter). Notably, even though the laryngeal and facial regions, but not for the pelvis (P = 8.0 × 10−5, P = 3.6 × 10−3, expression and methylation data come from different individuals, plotting them P = 8.2 × 10−4, and P = 0.115, for vocal folds, larynx, face and pelvis, respectively, against one another positions them only ~one standard deviation from the permutation test, Supplementary Fig. 4b–e). expression value predicted by the regression line (Fig. 5b). Future studies providing Potentially, longer genes have higher probability to overlap DMRs. Indeed, RNA expression levels for the laryngeal skeleton and vocal folds might provide DMGs tend to be longer (148 vs. 39 kb, P = 9.9 × 10−145, t-test). We thus checked further information on the methylation-expression links of these genes. the possibility that genes affecting the larynx and face tend to be longer than other genes, and are thus more likely to contain DMRs. We found that length of genes Studying the function of DMGs. Gene ontology and expression analyses were could not be a factor explaining the enrichment within genes affecting the larynx, conducted using Biological Process and UNIGENE expression tools in DAVID, as these genes tend to be shorter than other genes in the genome (mean: 62.5 vs. using an FDR threshold of 0.05. 73.2 kb, P = 0.001, t-test). Genes affecting the face, on the other hand, tend to be longer than other genes (mean: 77.1 vs. 65.6 kb, P = 4.6 × 10−5, t-test). To examine if this factor may underlie the enrichment we observe, we repeated the analysis Gene ORGANizer analysis. Similarly to sequence mutations, changes in regula- using only DMRs that are found within promoter regions (5 kb upstream to 1 kb tion are likely to be unequally distributed across different body systems, owing to downstream of TSS), thus eliminating the gene length factor. We found that the negative and positive selection, as well as inherent traits of the genes affecting each genes where such DMRs occur are still significantly associated with the face (P = organ. Thus, we turned to investigate which body parts are affected by the DMGs. 0.036, Fisher's exact test). We next repeated the promoter DMR analysis for all To this end, we ran the lists of DMGs in Gene ORGANizer11, which is a tool that genes and compared the Gene ORGANizer enrichment levels in this analysis to the links genes to the organs they affect, through known disease and normal pheno- genome-wide analysis. We observed very similar levels of enrichment (2.02×, types. Thus, it allows us to investigate directly the phenotypic function of genes, to 1.67×, and 1.24×, for vocal folds, larynx, and face, respectively, albeit FDR identify their shared targets and to statistically test the significance of such values >0.05 due to low statistical power). Importantly, AMH-derived DMGs also enrichments. We ran the lists of DMGs in the ORGANize option using the default do not tend to be longer than DMGs on the other branches (148 vs. 147 kb, P = parameters (i.e., based on confident and typical + non-typical gene-phenotype 0.93, t-test). Together, these analyses suggest that gene length does not affect the associations). observed enrichment in genes affecting the face and larynx. When we ran the list of skeletal AMH-derived DMRs, we found 11 significantly Additionally, to test whether cellular composition or differentiation state could enriched body parts, with the vocal folds and the larynx being the most enriched bias the results, we ran Gene ORGANizer on the list of DMGs, following the parts (×2.11 and ×1.68, FDR = 0.017 and FDR = 0.046, respectively). Most other removal of 20 DMRs that are found <10 kb from loci where methylation was parts belonged to the face (teeth, forehead, lips, eyelid, maxilla, face, and jaws), as shown to change during osteogenic differentiation57. We found that genes affecting well as the pelvis and nails (Fig. 2c, d and Supplementary Data 4). For archaic- the voice and face are still the most over-represented (2.13×, 1.71×, and 1.27×, derived DMGs, the lips, limbs, jaws, scapula, and spinal column were enriched FDR = 0.032, FDR = 0.049, and FDR = 0.040, for vocal folds, larynx, and face, (Supplementary Fig. 4e, Supplementary Data 4). The Neanderthal-derived and respectively, Supplementary Data 4). Denisovan-derived DMG lists did not produce any significantly enriched organs, We also investigated the possibility that (for an unknown reason) the DMR- but the immune system was significantly depleted within Neanderthal-derived detection algorithm introduces positional biases that preferentially identify DMRs DMRs (×0.67, FDR = 0.040). within genes affecting the voice or face. To this end, we simulated stochastic In order to examine whether such trends could arise randomly from the deamination processes along the Ust'-Ishim, Altai Neanderthal, and Denisovan reconstruction method, we repeated the analysis on the previously described genomes, reconstructed methylation maps, and ran the DMR-detection algorithm 100 simulations. We ran all simulated DMGs (4153) in Gene ORGANizer and on these maps. We repeated this process 100 times for each hominin and found no found that no enrichment was detected, neither for voice-related organs (vocal enrichment of any body part, including the face, vocal folds, or larynx (1.07×, folds: ×0.99, FDR = 0.731, larynx: ×1.02, FDR = 0.966, FDR = 0.966), nor for any 1.07×, and 1.04×, respectively, FDR = 0.88 for vocal folds, larynx, and face, other organ. permutation test). Perhaps most importantly, none of the other archaic branches shows enrichment of the larynx or the vocal folds. However, archaic-derived Validation of face and larynx enrichment in Gene ORGANizer. To test whether DMGs show over-representation of the jaws, as well as the lips, limbs, scapulae, the enrichment of the face and larynx could be attributed to the fact that the and spinal column (Supplementary Fig. 4e, Supplementary Data 4). In addition, NATURE COMMUNICATIONS | (2020)11:1189 | https://doi.org/10.1038/s41467-020-15020-6 | www.nature.com/naturecommunications 17 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15020-6 DMRs that separate chimpanzees from all humans (archaic and modern, Supplementary Fig. 4e, parts of the face are enriched within Archaic-derived Supplementary Data 2) do not show over-representation of genes that affect the DMGs. However, we currently see no substantial evidence supporting this. voice, larynx, or face, compatible with the notion that this trend emerged along the Importantly, the link between genetic alterations and phenotypes related to the AMH lineage. We also sought to test whether the larynx and vocal folds, which we voice is complex. Some brain-related disorders (i.e., clinical disorders that affect the found to be significantly enriched only along the AMH lineage, are also enriched brain) result in alterations to the voice, the mechanism in which is very difficult to when compared to the other lineages. We ran a chi-squared test on the fraction of pin down. Although the mechanism leading to voice alterations (either in its pitch, vocal folds-affecting and larynx-affecting AMH-derived DMGs (25 and 29, timbre, volume, or range) in some of the genes we report is unknown, many of the respectively, out of a total of 120 organ-associated DMGs), compared to the disorders are skeletal, suggesting the mechanism is related to anatomical changes to corresponding fraction in the DMGs along all the other lineages (42 for vocal folds, the vocal tract. Such changes could also affect more primary functions of the 49 for larynx, out of a total of 275 organ-associated DMGs). We found that both larynx, such as swallowing and breathing. However, the enrichment we observe in the larynx and vocal folds are significantly enriched in AMHs by over 50% Gene ORGANizer shows these genes were also shown to drive vocal alterations in compared to the other lineages (1.57× for both, P = 0.0248 and P = 0.0169 for the disorders they underlie11,19. Voice and speech alterations were also shown to be vocal folds and larynx, respectively, chi-squared test). driven by cultural, dietary and behavioral changes affecting bite configuration74. Furthermore, we added a human bone reduced representation bisulfite Here too, these factors are unlikely to underlie the vocal alterations in the genes we sequencing (RRBS) map, and produced a RRBS map from a chimpanzee infant report, as individuals from the same family as the individual with the disorder, who unspecified long bone (Supplementary Data 1). RRBS methylation maps include do not carry the dysfunctional allele, were not reported to present any vocal information on only ~10% of CpG sites, and are biased towards unmethylated sites. phenotypes. Therefore, they were not included in the previous analyses. However, we added The larynx is an organ which is primarily involved in breathing and swallowing them in this part as they originate from a chimpanzee infant and a present-day in mammals. In humans, the larynx is also used to produce complex speech, but human that is of similar age to the Denisovan (Supplementary Data 1), allowing not every change to the larynx necessarily affects speech. Despite these additional sampling from individuals that are younger than the rest. Repeating the Gene functions, the genes reported by Gene ORGANizer and HPO were specifically ORGANizer analysis after including these samples in the filtering process, we associated with voice alterations, directly or indirectly, suggesting that although found that the face and larynx are the only significantly enriched skeletal regions, they could have additional effects, their effect on the voice is their most shared and the enrichment within voice-affecting genes becomes even more pronounced function. (2.33×, FDR = 7.9 × 10−3, Supplementary Data 4). We also examined if pleiotropy could underlie the observed enrichments. To a Overlap with enhancer regions. To further test whether the AMH-derived DMRs large extent, the statistical tests behind Gene ORGANizer inherently account for overlap skeletal regulatory regions, we examined the previously reported 403,968 pleiotropy11, hence the conclusion that the most significant shared effect of the human loci, where an enrichment of the active enhancer mark H3K27ac was AMH-derived DMGs is in shaping vocal and facial anatomy is valid regardless of detected in developing human limbs (E33, E41, E44, and E47)75. Each DMR was pleiotropy. Nevertheless, we tested this possibility more directly, estimating the allocated a random genomic position in its original chromosome, while keeping its pleiotropy of each gene by counting the number of different Human Phenotype original length and matching the distribution of GC-content and CpG density Ontology (HPO) terms that are associated with it across the entire body19. We between the original and permutated lists. GC-content and CpG density matching found that DMGs do not tend to be more pleiotropic than the rest of the genome was done by matching a 10-bin histogram of the original and permutated lists. This (P = 0.17, t-test), nor do differentially methylated voice-affecting and face-affecting was repeated for 10,000 iterations. We found that AMH-derived DMRs overlap genes tend to be more pleiotropic than other DMGs (P = 0.19 and P = 0.27, limb H3K27ac-enriched regions ~2× more often than expected by chance (610 respectively, t-test). overlapping DMRs, compared to 312.4 ± 21.7, P < 10−4, permutation test). Next, we tested whether the process of within-lineage removal of variable DMRs SOX9 upstream putative enhancer coordinates used in Fig. 4b were taken and the differential number of samples along each lineage biases the Gene from13,14,23,76,77. ORGANizer enrichment analysis. To do so, we analyzed the pre-filtering DMRs along each lineage. We detect very similar trends to the post-filtering analysis, with the laryngeal and facial regions being the most significantly enriched within AMH- Computing the density of changes along the genome. We computed the density derived DMRs (1.58×, 1.44×, and 1.21× to 1.31× for the vocal folds, larynx and of derived CpG positions along the genome in two ways. First, we used a 100 kb different facial regions, respectively, FDR < 0.05), and for archaic-derived DMRs, window centered in the middle of each DMR and computed the fraction of CpGs we detect no enrichment of the laryngeal region (FDR = 0.16 and FDR = 0.43 for in that window which are differentially methylated (i.e., are found within a DMR). the vocal folds and larynx, respectively), and the most enriched regions are the face, Second, for the chromosome density plots, we did not center the window around limbs, and urethra. With the exception of the urethra, these results are very similar each DMR, but rather used a non-overlapping sliding 100 kb window starting at to the results reported for the filtered DMRs, suggesting that the process of within- position 1 and running the length of the chromosome. lineage removal of variable DMRs and the differential number of samples along each lineage does not bias the enrichment results. NFIX, COL2A1, SOX9, ACAN, and XYLT1 phenotypes. The vocal tract and Overall, we observe that AMH-derived DMGs across all 60 AMH samples are larynx affecting genes presented in this paper show involvement in laryngeal found outside archaic human variability, regardless of bone type, disease state, age, cartilage and soft tissue phenotypic variation. Clinical phenotypes can be of high or sex, and that chimpanzee methylation levels in these DMGs cluster closer to severity, with substantial impacts on normal breathing functions, to the point archaic humans than to AMHs, suggesting that these factors are unlikely to where the cause of death is due to respiratory distress. SOX9 and NFIX are often underlie the observed trends. Finally, we tested whether the filtering process in associated with laryngomalacia11,19 (Supplementary Data 5), a collapse of the itself might underlie the observed trends. To this end, we re-ran the entire pipeline larynx due to malformation of the laryngeal cartilaginous framework and/or on Neanderthal-derived and Denisovan-derived DMGs, while applying to them all malformed connective tissues, particularly during inhalation. Patients with muta- the filters as if they were Ust’-Ishim DMGs. This resulted in substantially fewer loci tions in COL2A1 often show backwards displacement of the tongue base11,19. Less (89 for the Neanderthal and 50 for the Denisovan), which limits statistical power, severe phenotypes of the reported genes include variation of voice quality in the but can still be used to examine whether there are any trends of enrichment similar form of pitch variation (high in patients suffering from XYLT1 mutations) and to those observed in AMHs. We found no evidence that the filtering process could sometimes hoarseness of the voice (reported for some patients with mutations of drive the enrichment of the vocal or facial areas: within Neanderthal-derived loci, ACAN, Supplementary Data 5)11,19. Whether this is due to variation of the vocal filtered as if they were Ust’-Ishim-derived, we found that the vocal folds were tract and laryngeal anatomy influenced by the ACAN mutation or due to a scaled ranked only 18th, with a non-significant enrichment of 1.27× (FDR = 0.815, down vocal tract size in the case of the XYLT1 mutation which also causes pri- compared to an enrichment of 2.11× within AMH-derived DMGs). The larynx was mordial dwarfism is not yet clear. ranked 76th and showed a non-significant depletion of 0.87× (FDR = 0.783), and the face was ranked 31st, with a non-significant enrichment of 1.09× (FDR = 0.815). Within Denisovan-derived loci, filtered as if they were Ust’-Ishim-derived, NFIX phenotypes. Skeletal phenotypes that are associated with the none of the loci were linked to the vocal folds nor to the larynx (FDR = 0.535 and Marshall–Smith syndrome were extracted from the Human Phenotype Ontology FDR = 0.834, respectively), and the face was ranked 30th (1.29x, FDR = 0.535, (HPO)19. Non-directional phenotypes (e.g., irregular dentition) and phenotypes Supplementary Data 4). This test suggests that the filtering process in itself is very that are expressed in both directions (e.g., tall stature and short stature) were unlikely to underlie the enrichment of the vocal and facial parts within AMH- removed. derived DMGs. Mutations in NFIX have also been linked to the Sotos syndrome. However, Next, we applied the Neanderthal/Denisovan filters to the Ust’-Ishim-derived NFIX is not the only gene that was linked to this syndrome; mutations in NSD1 loci. This resulted in 792 loci. We found that the vocal folds remained the most were also shown to drive similar phenotypes52. Therefore, it is less relevant in enriched body part (1.76×, FDR = 0.032), the larynx was marginally significant assessing the functional consequences of general shifts in the activity levels of (1.53×, FDR = 0.0502), and the facial region was significantly enriched too (e.g., NFIX. Nevertheless, it is noteworthy that in the Sotos syndrome too, most cheek and chin ranked 2nd, 3rd within significantly enriched body parts, 1.66× and symptoms are a mirror image of the Neanderthal phenotype (e.g., prominent chin 1.63×, FDR = 0.031 and FDR = 0.013, respectively, Supplementary Data 4). and high forehead). Importantly, we do not rule out the option that extensive regulatory changes in genes related to vocal and facial anatomy might have occurred along the Comparing gene expression between AMH and mouse. Ninety-three appen- Neanderthal and Denisovan lineages as well. Indeed, as we report in dicular skeleton samples were used to compare expression levels of NFIX, SOX9, 18 NATURE COMMUNICATIONS | (2020)11:1189 | https://doi.org/10.1038/s41467-020-15020-6 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15020-6 ARTICLE ACAN, and COL2A1 in human and mouse: 1. Five Human expression array data of Received: 25 February 2019; Accepted: 13 February 2020; iliac bones78, downloaded from ArrayExpress accession number E-MEXP-2219. 2. Eighty-four human expression array of iliac bones, downloaded from ArrayExpress accession number E-MEXP-1618. 3. Three mouse expression array data of femur and tibia bones, downloaded from ArrayExpress accession number E-GEOD- 61146. 4. One mouse RNA-seq of a tibia bone, downloaded from supplementary data. Expression values were converted to percentiles, according to each gene expression level compared to the rest of the genome across each sample (Fig. 5c). References 1. Prüfer, K. et al. The complete genome sequence of a Neanderthal from the Altai Mountains. Nature 505, 43–49 (2014). Methylation in AMH, chimpanzee, and Neanderthal femora. To check whether 2. King, M. C. & Wilson, A. C. Evolution at two levels in humans and the AMH hypermethylation of SOX9, ACAN, COL2A1, XYLT1, and NFIX could be chimpanzees. Science 188, 107–116 (1975). a result of variability between bone types, we compared the four chimpanzee femur 3. Hernando-Herraez, I., Garcia-Perez, R., Sharp, A. J. & Marques-Bonet, T. 850K methylation arrays to the 52 present-day femur 450K methylation arrays. We DNA methylation: insights into human evolution. PLoS Genet. 11, (2015). took probes within AMH-derived DMRs that appear on both arrays. We found 4. Gokhman, D. et al. Reconstructing the DNA methylation maps of the that these genes are consistently hypermethylated in AMHs (P = 1.6 × 10−7, t-test), Neandertal and the Denisovan. Science 344, 523–527 (2014). with 38 probes showing >5% hypermethylation in AMH, whereas only eight probes 5. Prüfer, K. et al. A high-coverage Neandertal genome from Vindija Cave in show such hypermethylation in chimpanzees (Supplementary Fig. 5d). Therefore, Croatia. Science 358, 655–658 (2017). even when comparing methylation from the same bone, same sex, same devel- 6. Fu, Q. et al. Genome sequence of a 45,000-year-old modern human from opmental stage, measured by the same technology, and across the same positions, western Siberia. Nature 514, 445–449 (2014). AMH show consistent hypermethylation across all of these DMGs. 7. Lazaridis, I. et al. Ancient human genomes suggest three ancestral populations Similarly, when comparing the DMRs in SOX9, ACAN, COL2A1, XYLT1, and for present-day Europeans. Nature 513, 409–413 (2014). NFIX between the Ust'-Ishim and Vindija Neanderthal samples, the Vindija 8. Olalde, I. et al. Derived immune and ancestral pigmentation alleles in a 7,000- Neanderthal sample is consistently hypomethylated compared to the Ust'-Ishim year-old Mesolithic European. Nature 507, 225–228 (2014). individual (P = 1.2 × 10−5, Supplementary Fig. 5c, t-test). Both of these samples were extracted from femora of adult individuals, and methylation was 9. Mathieson, I. et al. Genome-wide patterns of selection in 230 ancient reconstructed using the same technology. This suggests that the hypermethylation Eurasians. Nature 528, 499–503 (2015). of AMHs compared to Neanderthals is unlikely to be driven by age or bone type, 10. Gokhman, D., Malul, A. & Carmel, L. Inferring past environments from and rather reflects evolutionary shifts. ancient epigenomes. Mol. Biol. Evol. 34, 2429–2438 (2017). 11. Gokhman, D. et al. Gene ORGANizer: linking genes to the organs they affect. Nucleic Acids Res. 45, W138–W145 (2017). Scanning SOX9 for mutations altering NFI binding motifs. To examine whether 12. Aiello, L. & Dean, C. An Introduction to Human Evolutionary Anatomy the changes in regulation of SOX9 could possibly be explained by changes in the (Elsevier, 2002). binding sites of NFI proteins, we searched for the NFI motif along the gene body 13. Prescott, S. L. et al. Enhancer divergence and cis-regulatory evolution in the and the 350 kb upstream region of SOX9. We looked for NFI motifs that exist in human and chimp neural crest. Cell 163, 68–84 (2015). the genomes of the Altai and Vindija Neanderthal, as well as in the Denisovan, but 14. Wilderman, A., VanOudenhove, J., Kron, J., Noonan, J. P. & Cotney, J. High- were abolished in AMHs. We did not find any evidence of such substitutions. resolution epigenomic Atlas of human embryonic craniofacial development. Cell Rep. 23, 1581–1597 (2018). Comparison to divergent traits between Neanderthals and AMHs. To further 15. MacArthur, J. et al. The new NHGRI-EBI Catalog of published genome-wide investigate potential phenotypic consequences of the DMGs we report, we association studies (GWAS Catalog). Nucleic Acids Res. 45, D896–D901 probed the HPO database19 and compared these HPO phenotypes to known (2017). morphological differences between Neanderthals and modern humans32. To 16. Hernando-Herraez, I. et al. The interplay between DNA methylation and compile a list of traits in which Neanderthals and AMHs differ, we reviewed key sequence divergence in recent human evolution. Nucleic Acids Res. 43, sources that surveyed Neanderthal morphology summarized in Aiello, L. and 8204–8214 (2015). Dean12. We identified traits in which Neanderthals are found completely outside 17. Schultz, M. D. et al. Human body epigenome maps reveal noncanonical DNA AMH variation, as well as traits where one group is significantly different from methylation variation. Nature 523, 212–216 (2015). the other, but the distribution of observed measurements partially overlap. Non- 18. Peyrégne, S., Boyle, M. J., Dannemann, M. & Prüfer, K. Detecting ancient directional traits (i.e., traits that could not be described on scales such as higher/ positive selection in humans using extended lineage sorting. Genome Res. 27, lower, accelerated/delayed etc.) were not included, as could not be paralleled 1563–1572 (2017). with HPO phenotypes. The compiled list included 107 phenotypes, 75 of which 19. Köhler, S. et al. The Human Phenotype Ontology project: linking molecular have at least one equivalent HPO phenotype (4.8 on average). For example, the biology and disease through phenotype data. Nucleic Acids Res. 42, HPO phenotype Taurodontia (HP:0000679) was linked to the trait Taurodontia, D966–D974 (2014). and the following HPO phenotypes were linked to the trait Rounded and robust 20. Lee, Y. H. & Saint-Jeannet, J. P. Sox9 function in craniofacial development and rib shafts: broad ribs (HP:0000885), hypoplasia of first ribs (HP:0006657), short disease. Genesis 49, 200–208 (2011). ribs (HP:0000773), thickened cortex of long bones (HP:0000935), thickened ribs 21. Ng, L. J. et al. Sox9 binds DNA, activates transcription, and coexpresses with (HP:0000900), thin ribs (HP:0000883), thoracic hypoplasia (HP:0005257). For each skeleton-affecting phenotype, we determined whether it matches a known type II collagen during chondrogenesis in the mouse. Dev. Biol. 183, 108–121 morphological difference between Neanderthals and AMHs32. For example, (1997). Hypoplastic ilia (HPO ID: HP:0000946) was marked as divergent because in the 22. Claes, P. et al. Genome-wide mapping of global-to-local genetic effects on Neanderthal the iliac bones are considerably enlarged compared to AMHs12. We human facial shape. Nat. Genet. 50, 414–423 (2018). then counted for each gene (whether DMG or not) the fraction of its associated 23. Yao, B. et al. The SOX9 upstream region prone to chromosomal aberrations HPO phenotypes that are divergent between Neanderthals and AMHs. We causing campomelic dysplasia contains multiple cartilage enhancers. Nucleic found that four of the top five most differentially methylated skeletal genes Acids Res. 43, 5394–5408 (2015). (XYLT1, NFIX, ACAN, and COL2A1) are in the top 100 genes with the highest 24. Kim, K. Il, Park, Y. S. & Im, G. Il Changes in the epigenetic status of the SOX- fraction of divergent traits between Neanderthals and AMHs (out of a total of 9 promoter in human osteoarthritic cartilage. J. Bone Miner. Res. 28, 1,789 skeleton-related genes). In fact, COL2A1, which is the top ranked DMR 1050–1060 (2013). (Supplementary Data 2), is also the gene that is overall associated with the 25. Claes, P. et al. Modeling 3D facial shape from DNA. PLoS Genet. 10, e1004224 highest number of derived traits (63) (Supplementary Data 7). This suggests that (2014). these extensive methylation changes are possibly linked to phenotypic diver- 26. Cheng, P. F. et al. Methylation-dependent SOX9 expression mediates invasion gence between archaic and AMHs. in human melanoma cells and is a negative prognostic factor in advanced melanoma. Genome Biol. 16, 42 (2015). Reporting summary. Further information on research design is available in 27. Huang, C.-Z. et al. Sox9 transcriptionally regulates Wnt signaling in intestinal the Nature Research Reporting Summary linked to this article. epithelial stem cells in hypomethylated crypts in the diabetic state. Stem Cell Res. Ther. 8, 60 (2017). 28. Pamnani, M., Sinha, P., Singh, A., Nara, S. & Sachan, M. Methylation of the Data availability Sox9 and Oct4 promoters and its correlation with gene expression during All sequencing and methylation data generated in this work have been deposited in testicular development in the laboratory mouse. Genet. Mol. Biol. 39, 452–458 NCBI’s Gene Expression Omnibus under GEO accession number GSE96833. All other (2016). data and materials are contained in the paper and its supplementary information or 29. Hall, B. D. Lethality in Desbuquois dysplasia: three new cases. Pediatr. Radiol. available upon request. 31, 43–47 (2001). NATURE COMMUNICATIONS | (2020)11:1189 | https://doi.org/10.1038/s41467-020-15020-6 | www.nature.com/naturecommunications 19 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15020-6 30. LaCroix, A. J. et al. GGC repeat expansion and exon 1 methylation of XYLT1 63. Kulis, M. et al. Whole-genome fingerprint of the DNA methylome during Is a common pathogenic variant in Baratela-Scott syndrome. Am. J. Hum. human B cell differentiation. Nat. Genet. 47, 746–756 (2015). Genet. 104, 35–44 (2019). 64. Marco-Sola, S., Sammeth, M., Guigó, R. & Ribeca, P. The GEM mapper: fast, 31. Ohba, S., He, X., Hojo, H. & McMahon, A. P. Distinct transcriptional accurate and versatile alignment by filtration. Nat. Methods 9, 1185–1188 programs underlie Sox9 regulation of the mammalian chondrocyte. Cell Rep. (2012). 12, 229–243 (2015). 65. Hansen, K. D., Langmead, B. & Irizarry, R. A. BSmooth: from whole genome 32. Gokhman, D. et al. Reconstructing Denisovan anatomy using DNA bisulfite sequencing reads to differentially methylated regions. Genome Biol. methylation maps. Cell 179, 180–192 (2019). 13, R83 (2012). 33. Roadmap Epigenomics Consortium. et al. Integrative analysis of 111 reference 66. Dabney, J. et al. Complete mitochondrial genome sequence of a Middle human epigenomes. Nature 518, 317–329 (2015). Pleistocene cave bear reconstructed from ultrashort DNA fragments. Proc. 34. Carrió, E. et al. Deconstruction of DNA methylation patterns during Natl Acad. Sci. USA 110, 15758–15763 (2013). myogenesis reveals specific epigenetic events in the establishment of the 67. Rohland, N. & Hofreiter, M. Ancient DNA extraction from bones and teeth. skeletal muscle lineage. Stem Cells 33, 2025–2036 (2015). Nat. Protoc. 2, 1756–1762 (2007). 35. Maunakea, A. K. et al. Conserved role of intragenic DNA methylation in 68. Boyle, P. et al. Gel-free multiplexed reduced representation bisulfite sequencing regulating alternative promoters. Nature 466, 253–257 (2010). for large-scale DNA methylation profiling. Genome Biol. 13, R92 (2012). 36. Ford, E. E. et al. Frequent lack of repressive capacity of promoter DNA 69. Barnett, R. & Larson, G. A phenol-chloroform protocol for extracting DNA methylation identified through genome-wide epigenomic manipulation. from ancient samples. Methods Mol. Biol. 840, 13–19 (2012). Preprint at bioRxiv https://doi.org/10.1101/170506 (2018). 70. Fortin, J. P., Triche, T. J. & Hansen, K. D. Preprocessing, normalization and 37. Housman, G., Havill, L. M., Quillen, E. E., Comuzzie, A. G. & Stone, A. C. integration of the Illumina HumanMethylationEPIC array with minfi. Assessment of DNA methylation patterns in the bone and cartilage of a Bioinformatics 33, 558–560 (2017). nonhuman primate model of osteoarthritis. Cartilage 10, 335–345 (2018). 71. Hernando-Herraez, I. et al. Dynamics of DNA methylation in recent human 38. Pjanic, M. et al. Nuclear Factor I genomic binding associates with chromatin and great ape evolution. PLoS Genet. 9, e1003763 (2013). boundaries. BMC Genomics 14, 99 (2013). 72. Sirak, K. A. et al. A minimally-invasive method for sampling human petrous 39. Tompson, S. W. et al. A recessive skeletal dysplasia, SEMD aggrecan type, bones from the cranial base for ancient DNA analysis. Biotechniques 62, results from a missense mutation affecting the C-type lectin domain of 283–289 (2017). aggrecan. Am. J. Hum. Genet. 84, 72–79 (2009). 73. Hanghøj, K. et al. Fast, accurate and automatic ancient nucleosome and 40. Frenzel, K., Amann, G. & Lubec, B. Deficiency of laryngeal collagen type II in methylation maps with epiPALEOMIX. Mol. Biol. Evol. 33, 3284–3298 (2016). an infant with respiratory problems. Arch. Dis. Child. 78, 557–559 (1998). 74. Blasi, D. E. et al. Human sound systems are shaped by post-Neolithic changes 41. Hoornaert, K. P. et al. Stickler syndrome caused by COL2A1 mutations: in bite configuration. Science 363, eaav3218 (2019). genotype–phenotype correlation in a series of 100 patients. Eur. J. Hum. 75. Guo, M. et al. Epigenetic profiling of growth plate chondrocytes sheds insight Genet. 18, 872–880 (2010). into regulatory genetic variation influencing height. elife 6, e29329 (2017). 42. Shaw, A. C. et al. Phenotype and natural history in Marshall-Smith syndrome. 76. Bagheri-Fam, S. et al. Long-range upstream and downstream enhancers Am. J. Med. Genet. A 152, 2714–2726 (2010). control distinct subsets of the complex spatiotemporal Sox9 expression 43. Gokhman, D., Meshorer, E. & Carmel, L. Epigenetics: it’s getting old. past pattern. Dev. Biol. 291, 382–397 (2006). meets future in paleoepigenetics. Trends Ecol. Evol. 31, 290–300 (2016). 77. Sekido, R. & Lovell-Badge, R. Sex determination involves synergistic action of 44. Lieberman, P. The evolution of human speech: its anatomical and neural SRY and SF1 on a specific Sox9 enhancer. Nature 453, 930–934 (2008). bases. Curr. Anthropol. 48, 39–66 (2007). 78. Varanasi, S. S. et al. Skeletal site-related variation in human trabecular bone 45. Fitch, W. T., de Boer, B., Mathur, N. & Ghazanfar, A. Response to Lieberman transcriptome and signaling. PLoS ONE 5, e10692 (2010). on “Monkey vocal tracts are speech-ready”. Sci. Adv. 3, e1701859 (2017). 46. Boë, L.-J. et al. Which way to the dawn of speech?: Reanalyzing half a century of debates and data in light of speech science. Sci. Adv. 5, eaaw3916 (2019). Acknowledgements 47. Lieberman, P. Comment on “Monkey vocal tracts are speech-ready”. Sci. Adv. We would like to thank Sagiv Shifman, Yoel Rak, Rodrigo Lacruz, Erella Hovers, Anna 3, e1700442 (2017). Belfer-Cohen, Achinoam Blau, Iain Mathieson, Philip Lieberman, Daniel Lieberman, and 48. Lieberman, D. E. The Evolution of the Human Head (Harvard Univ. Press, Terry Capellini for their useful advice, Svante Pääbo, Janet Kelso, Kay Prüfer, Johannes 2011). Krause, and Anne Pusey for providing data, Sjaak Kaandorp and Christine Kaandorp- 49. Nishimura, T., Mikami, A., Suzuki, J. & Matsuzawa, T. Descent of the hyoid in Huber from Safari park Beekse Bergen in Netherlands for their cooperation in animal chimpanzees: evolution of face flattening and speech. J. Hum. Evol. 51, conservation and use for research, and Maayan Harel for illustrations. D.G. is supported 244–254 (2006). by the Clore Israel Foundation. TMB is supported by BFU2017-86471-P (MINECO/ 50. Fitch, W. T. The evolution of speech: a comparative review. Trends Cogn. Sci. FEDER, UE), U01 MH106874 grant, Howard Hughes International Early Career, Obra 4, 258–267 (2000). Social “La Caixa” and Secretaria d’Universitats i Recerca and CERCA Program del 51. Steele, J., Clegg, M. & Martelli, S. Comparative morphology of the hominin Departament d’Economia i Coneixement de la Generalitat de Catalunya. D.R. is an and african ape hyoid bone, a possible marker of the evolution of speech. Investigator of the Howard Hughes Medical Institute and is also supported by an Allen Hum. Biol. 85, 639–672 (2013). Discovery Center for the Study of Human Brain Evolution funded the Paul G. Allen 52. Malan, V. et al. Distinct effects of allelic NFIX mutations on nonsense- Family Foundation. C.L.-F. is supported by FEDER and BFU2015-64699-P grant from mediated mRNA decay engender either a sotos-like or a Marshall-Smith the Spanish government. R.P. was supported by ERC starting grant ADNABIOARC Syndrome. Am. J. Hum. Genet. 87, 189–198 (2010). (263441). R.M.G. and J.M.O. are supported by NYSTEM contract C030133. Funding for 53. Van Balkom, I. D. C. et al. Development and behaviour in Marshall-Smith the collection and processing of the 850K chimpanzee data was provided by the Leakey syndrome: an exploratory study of cognition, phenotype and autism. J. Foundation Research Grant for Doctoral Students, Wenner-Gren Foundation Disserta- Intellect. Disabil. Res. 55, 973–987 (2011). tion Fieldwork Grant (Gr. 9310), James F. Nacey Fellowship from the Nacey Maggion- 54. Qiu, Y. et al. Oligogenic effects of 16p11.2 copy number variation on calda Foundation, International Primatological Society Research Grant, Sigma Xi Grant- craniofacial development. Cell Rep. 28, 3320–3328.e4 (2019). in-Aid of Research, Center for Evolution and Medicine Venture Fund (ASU), Graduate 55. Cullen á, A., Clarke á, T. A. & O’Dwyer, T. P. The Marshall-Smith syndrome: Research and Support Program Grant (GPSA, ASU), and Graduate Student Research a review of the laryngeal complications. Eur. J. Pediatr. 156, 463–464 (1997). Grant (SHESC, ASU) to G.H. Collection of the chimpanzee bone from Tanzania was 56. Ziller, M. J. et al. Charting a dynamic DNA methylation landscape of the funded by the Jane Goodall Institute, and grants from the US National Institutes of human genome. Nature 500, 477–481 (2013). Health (AI 058715) and National Science Foundation (IOS-1052693), and facilitated by 57. Håkelien, A. M. et al. The regulatory landscape of osteogenic differentiation. Elizabeth Lonsdorf and Beatrice Hahn. Stem Cells 32, 2780–2793 (2014). 58. Bernstein, B. E. et al. An integrated encyclopedia of DNA elements in the human genome. Nature 489, 57–74 (2012). 59. Moriarity, B. S. et al. A sleeping beauty forward genetic screen identifies new Author contributions genes and pathways driving osteosarcoma development and metastasis. Nat. D.G. planned and conducted analyses. L.C. supervised the computational and experi- Genet. 47, 615–624 (2015). mental work. E.M. supervised experiments. L.A.T., B.Y., D.G., and L.C. conceived sta- 60. Lacruz, R. S. et al. Ontogeny of the maxilla in Neanderthals and their tistical analyses. D.G., L.C., and E.M. wrote the manuscript. All other authors contributed ancestors. Nat. Commun. 6, 8996 (2015). to the production of data and wrote their respective parts of the manuscript. 61. Nishimura, T. Developmental changes in the shape of the supralaryngeal vocal tract in Chimpanzees. Am. J. Phys. Anthropol. 126, 193–204 (2005). 62. Horvath, S. et al. The cerebellum ages slowly according to the epigenetic clock. Competing interests Aging 7, 294–306 (2015). The authors declare no competing interests. 20 NATURE COMMUNICATIONS | (2020)11:1189 | https://doi.org/10.1038/s41467-020-15020-6 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15020-6 ARTICLE Additional information Open Access This article is licensed under a Creative Commons Supplementary information is available for this paper at https://doi.org/10.1038/s41467- Attribution 4.0 International License, which permits use, sharing, 020-15020-6. adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Correspondence and requests for materials should be addressed to D.G., E.M. or L.C. Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless Peer review information Nature Communications thanks the anonymous reviewers for indicated otherwise in a credit line to the material. If material is not included in the their contribution to the peer review of this work. Peer reviewer reports are available. article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from Reprints and permission information is available at http://www.nature.com/reprints the copyright holder. To view a copy of this license, visit http://creativecommons.org/ licenses/by/4.0/. Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. © The Author(s) 2020 David Gokhman1 ✉, Malka Nissim-Rafinia1, Lily Agranat-Tamir1,2, Genevieve Housman 3,4, Raquel García-Pérez5, Esther Lizano 5, Olivia Cheronet 6, Swapan Mallick7,8,9, Maria A. Nieves-Colón3,4, Heng Li7, Songül Alpaslan-Roodenberg8, Mario Novak10,11, Hongcang Gu7, Jason M. Osinski12, Manuel Ferrando-Bernal5, Pere Gelabert 5, Iddi Lipende13, Deus Mjungu13, Ivanela Kondova 14, Ronald Bontrop 14, Ottmar Kullmer 15, Gerhard Weber6, Tal Shahar16, Mona Dvir-Ginzberg 17, Marina Faerman18, Ellen E. Quillen19, Alexander Meissner 7,20,21, Yonatan Lahav22,23, Leonid Kandel24, Meir Liebergall24, María E. Prada25, Julio M. Vidal 26, Richard M. Gronostajski 12,27, Anne C. Stone 3,4,28, Benjamin Yakir2, Carles Lalueza-Fox 5, Ron Pinhasi 6, David Reich 7,8,9, Tomas Marques-Bonet 5,29,30,31, Eran Meshorer 1,32 ✉ & Liran Carmel 1 ✉ 1 Department of Genetics, The Alexander Silberman Institute of Life Sciences, Faculty of Science, The Hebrew University of Jerusalem, 91904 Jerusalem, Israel. 2Department of Statistics, The Hebrew University of Jerusalem, 91905 Jerusalem, Israel. 3School of Human Evolution and Social Change, Arizona State University, Tempe, AZ 85281, USA. 4Center for Evolution and Medicine, Arizona State University, Tempe, AZ 85287, USA. 5 Institute of Evolutionary Biology (UPF-CSIC), PRBB, 08003 Barcelona, Spain. 6Department of Evolutionary Anthropology, University of Vienna, 1090 Vienna, Austria. 7Broad Institute, Cambridge, MA 02138, USA. 8Department of Genetics, Harvard Medical School, Boston, MA 02115, USA. 9 Howard Hughes Medical Institute, Harvard Medical School, Boston, MA 02115, USA. 10Institute for Anthropological Research, 10000 Zagreb, Croatia. 11Earth Institute and School of Archaeology, University College Dublin, Dublin 4, Ireland. 12Department of Biochemistry, Jacobs School of Medicine and Biomedical Sciences, University at Buffalo, Buffalo, NY 14203, USA. 13Gombe Stream Research Center, Jane Goodall Institute, Kigoma, Tanzania. 14Biomedical Primate Research Centre (BPRC), Rijswijk, Netherlands. 15Department of Palaeoanthropology and Messel Research, Senckenberg Center of Human Evolution and Paleoecology, Frankfurt am Main, Germany. 16Department of Neurosurgery, Shaare Zedek Medical Center, Jerusalem, Israel. 17Laboratory of Cartilage Biology, Institute of Dental Sciences, Faculty of Dental Medicine, Hebrew University of Jerusalem, 91120 Jerusalem, Israel. 18Laboratory of Bioanthropology and Ancient DNA, Institute of Dental Sciences, Faculty of Dental Medicine, The Hebrew University of Jerusalem, 91120 Jerusalem, Israel. 19Department of Genetics, Texas Biomedical Research Institute, San Antonio, TX 85287, USA. 20Harvard Stem Cell Institute, Cambridge, MA 02138, USA. 21Department of Stem Cell and Regenerative Biology, Harvard University, Cambridge, MA 02138, USA. 22Otolaryngology - Head & Neck Surgery Department, Laryngeal Surgery Unit, Kaplan Medical Center, Rehovot, Israel. 23 The Hebrew University Medical School, Jerusalem, Israel. 24Orthopaedic Department, Hadassah – Hebrew University Medical Center, Jerusalem, Israel. 25I.E.S.O. ‘Los Salados’. Junta de Castilla y León, León, Spain. 26Junta de Castilla y León, Servicio de Cultura de León, León, Spain. 27 Genetics, Genomics and Bioinformatics Program, New York State Center of Excellence in Bioinformatics and Life Sciences, Jacobs School of Medicine and Biomedical Sciences, University at Buffalo, Buffalo, NY 14203, USA. 28Institute of Human Origins, Arizona State University, Tempe, AZ 85287, USA. 29Catalan Institution of Research and Advanced Studies (ICREA), 08010 Barcelona, Spain. 30CNAG-CRG, Centre for Genomic Regulation (CRG), Barcelona Institute of Science and Technology (BIST), 08028 Barcelona, Spain. 31Institut Català de Paleontologia Miquel Crusafont, Universitat Autònoma de Barcelona, Edifici ICTA-ICP, c/ Columnes s/n, Barcelona, Spain. 32The Edmond and Lily Safra Center for Brain Sciences (ELSC), The Hebrew University of Jerusalem, 91904 Jerusalem, Israel. ✉email: david.gokhman@mail.huji.ac.il; eran.meshorer@mail.huji.ac.il; liran.carmel@huji.ac.il NATURE COMMUNICATIONS | (2020)11:1189 | https://doi.org/10.1038/s41467-020-15020-6 | www.nature.com/naturecommunications 21