ARTICLE https://doi.org/10.1038/s41467-021-23397-1 OPEN Epigenomic profiling of primate lymphoblastoid cell lines reveals the evolutionary patterns of epigenetic activities in gene regulatory architectures Raquel García-Pérez 1,11 ✉, Paula Esteller-Cucala 1,11, Glòria Mas2,3, Irene Lobón1, Valerio Di Carlo 2,3, Meritxell Riera1, Martin Kuhlwilm 1, Arcadi Navarro 1,4,5, Antoine Blancher 6,7, Luciano Di Croce 2,3,5, José Luis Gómez-Skarmeta8, David Juan 1,11 ✉ & Tomàs Marquès-Bonet 1,5,9,10,11 ✉ 1234567890():,; Changes in the epigenetic regulation of gene expression have a central role in evolution. Here, we extensively profiled a panel of human, chimpanzee, gorilla, orangutan, and macaque lymphoblastoid cell lines (LCLs), using ChIP-seq for five histone marks, ATAC-seq and RNA- seq, further complemented with whole genome sequencing (WGS) and whole genome bisulfite sequencing (WGBS). We annotated regulatory elements (RE) and integrated chromatin contact maps to define gene regulatory architectures, creating the largest catalog of RE in primates to date. We report that epigenetic conservation and its correlation with sequence conservation in primates depends on the activity state of the regulatory element. Our gene regulatory architectures reveal the coordination of different types of components and highlight the role of promoters and intragenic enhancers (gE) in the regulation of gene expression. We observe that most regulatory changes occur in weakly active gE. Remarkably, novel human-specific gE with weak activities are enriched in human-specific nucleotide changes. These elements appear in genes with signals of positive selection and human acceleration, tissue-specific expression, and particular functional enrichments, suggesting that the regulatory evolution of these genes may have contributed to human adaptation. 1 Institute of Evolutionary Biology (UPF-CSIC), PRBB, Barcelona, Spain. 2 Centre for Genomic Regulation (CRG), The Barcelona Institute of Science and Technology, Barcelona, Spain. 3 Universitat Pompeu Fabra (UPF), Barcelona, Spain. 4 National Institute for Bioinformatics (INB), PRBB, Barcelona, Spain. 5 Institució Catalana de Recerca i Estudis Avançats (ICREA), Barcelona, Spain. 6 Laboratoire d’immunologie, CHU de Toulouse, Institut Fédératif de Biologie, hôpital Purpan, Toulouse, France. 7 Centre de Physiopathologie Toulouse-Purpan (CPTP), Université de Toulouse, Centre National de la Recherche Scientifique (CNRS), Institut National de la Santé et de la Recherche Médicale (Inserm), Université Paul Sabatier (UPS), Toulouse, France. 8 Centro Andaluz de Biología del Desarrollo (CABD), Consejo Superior de Investigaciones Científicas-Universidad Pablo de Olavide-Junta de Andalucía, Seville, Spain. 9 CNAG- CRG, Centre for Genomic Regulation (CRG), Barcelona Institute of Science and Technology (BIST), Barcelona, Spain. 10 Institut Català de Paleontologia Miquel Crusafont, Universitat Autònoma de Barcelona, Cerdanyola del Vallès, Barcelona, Spain. 11These authors contributed equally: Raquel García-Pérez, Paula Esteller-Cucala, David Juan, Tomàs Marquès-Bonet. ✉email: raquel.garcia@bsc.es; david.juan@upf.edu; tomas.marques@upf.edu NATURE COMMUNICATIONS | (2021)12:3116 | https://doi.org/10.1038/s41467-021-23397-1 | www.nature.com/naturecommunications 1 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-23397-1 C hanges in chromatin structure and gene regulation play a overlapping regulatory regions in the two replicates of every crucial role in evolution1,2. Gene expression differences species to define species RE. have been extensively studied in a variety of species and We classified the chromatin states of the RE based on a conditions3–9. However, there is still much unknown about how hierarchy of functionally interpretable epigenetic states. This regulatory landscapes evolve, even in closely related species. hierarchy differentiates chromatin states into promoter (P) and Previous work has focused on the dynamics of the addition and enhancer (E) states, with three different levels of activity each: removal of RE with signals of strong activity during mammalian strong (s), poised (p), or weak (w) (Methods and Supplementary evolution—mainly defined from ChIP-seq experiments on a few Fig. 1). We improved these assignments by applying a linear histone marks10–14. These analyses suggested that enhancers discriminative analysis (LDA) with normalized histone and open evolve faster than promoters. The number of active enhancers chromatin enrichments (Methods and Supplementary Figs. 3 and located near a gene—its regulatory complexity—has also been 4). The refined classification results in more similar regulatory reported to influence the conservation of gene expression in landscapes between biological replicates (Wilcoxon signed-rank mammals10. test: P < 0.05 in all species; Supplementary Fig. 5), increasing the Moreover, in a selected group of primates—mostly chimpan- number of RE with the same state in all species (Wilcoxon zees and macaques—changes in histone mark enrichments are signed-rank test: P = 0.03; Supplementary Fig. 6). At large, associated with gene expression differences15–18. Several studies promoters have a positive H3K4me3/H3K4me1 ratio, whereas have also targeted the appearance of human-specific H3K4me1 is more abundant in enhancers than promoters. RE methylation19–23 patterns and active promoters and enhancers with strong activities have an H3K27ac enrichment level similar in different anatomical structures and cell types11–15. All these to that of H3K4me3 in promoters and H3K4me1 in enhancers. studies have proven that comparative epigenomics is a powerful Poised RE have a characteristic enrichment of H3K27me3. tool to investigate the evolution of RE24,25. However, a deeper Finally, weak promoters and enhancers are associated with lower understanding of the evolution of gene regulation requires the intensity enrichments of all epigenetic signals and generally lack integration of multilayered epigenome data. Only such integra- H3K27me3 and H3K36me3 (Supplementary Fig. 3). tion can provide the necessary resolution of regulatory activities On average, we found ~11,000 and ~76,000 RE with promoter for investigating recent evolutionary time frames, as is the case and enhancer states per species, respectively (Fig. 2b), of which within the primate lineage. Here, we provide an in-depth com- 69% and 33% are strong, 8% and 4% are poised, and 14% and 45% parison of the recent evolution of gene regulatory architectures are weak, respectively (Supplementary Fig. 7 and Supplementary using a homologous cellular model system in human and non- Data 1). Strong and poised activities are more associated with human primates. We observe different levels of correlation promoter states, whereas weak activities are more frequently between sequence and epigenetic conservation for RE with dif- associated with enhancer states (Chi-square test: P < 2.2 × 10−16 in ferent activities and highlight the contribution of intragenic all species). We associated RE with genes using 1D gene proximity enhancers (gE) to explain gene expression levels. We also report and existing high-resolution 3D chromatin contact data for one of the role that often understudied epigenetic states may have in the human LCLs (Fig. 2a and Methods). On average, 70% of the recent human evolution. RE are associated with genes, of which an average across species of 93% of the RE are associated with protein-coding genes and 61% of the RE are associated with one-to-one orthologous protein- Results coding genes in all primate species (Fig. 2c). The set of RE Comprehensive profiling of primate lymphoblastoid cell lines associated with a gene defines its regulatory architecture. (LCLs). We have extensively characterized a panel of lympho- Altogether, this catalog of RE provides a comprehensive view blastoid cell lines (LCLs) from human, chimpanzee, gorilla, of the regulatory landscape of LCLs in humans and nonhuman orangutan, and macaque, including two independent biological primates. In contrast to other commonly used definitions of replicates for each species. This characterization includes chro- promoters and enhancers limited to strongly active regions, our matin immunoprecipitation data (ChIP-seq) from five key his- multilayered integration approach allows the additional annota- tone modifications (H3K4me1, H3K4me3, H3K36me3, H3K27ac, tion of weak and poised activities11,13. These activities are of and H3K27me3) and deep-transcriptome sequencing (RNA-seq) particular relevance to improve the definition of elements in (Fig. 1). We integrate these datasets into gene regulatory archi- regulatory gene architectures. In sum, a detailed primate tectures (Fig. 2a and Supplementary Figs. 1, 2) to (1) understand regulatory catalog enables the study of the evolution of these how primate gene expression levels are controlled and how regulatory activities using LCLs as a proxy of their regulatory expression changes between species occur and to (2) study pat- potential in other cell types or conditions. terns of evolutionary conservation of RE in primates. To com- plement this resource, we have also processed high coverage whole-genome and whole-genome bisulfite sequencing data, as The evolutionary dynamics of promoters and enhancers in well as chromatin accessibility data (Supplementary Tables 1–6 primate LCLs recapitulate previous observations in more dis- and Supplementary Data 1–4). Taken together, this is the most tant mammals. Interspecies differences in regulatory regions can extensive collection of great apes and macaque transcriptomic be associated with genomic or epigenetic changes. Inconsistencies in and epigenomic data to date. the quality of genome assemblies make it difficult to distinguish actual interspecies genomic differences, an issue aggravated in multispecies comparisons. To overcome this problem, we restricted Annotation of RE. We used the signal of the ChIP-seq experi- our analyses to unambiguous one-to-one orthologs between all ments from the five histone marks to identify putative regulatory species. We detected 28,703 one-to-one orthologous genomic regions with characteristic marks of promoters or enhancers regions in the five species with a promoter or enhancer state in at (Supplementary Figs. 1 and 2). We defined regulatory regions for least one species (Supplementary Fig. 8). Most of these orthologous each cell line as those containing chromatin states (over- regulatory regions (~76%, Binomial test: P < 2.2 × 10−16) are asso- represented combinations of histone marks detected by ciated with genes (Methods). In downstream analyses, we focused ChromHMM26) enriched in any regulatory-related histone mark on these regions integrating the regulatory architectures of protein- (Methods and Fig. 2a and Supplementary Fig. 1). We merged coding and non-coding genes. 2 NATURE COMMUNICATIONS | (2021)12:3116 | https://doi.org/10.1038/s41467-021-23397-1 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-23397-1 ARTICLE Fig. 1 Overview of the study design and data generated. a One human and eight nonhuman primate lymphoblastoid cell lines (LCLs) were cultured to perform a variety of high-throughput techniques, including whole genome sequencing (WGS), whole genome bisulfite sequencing (WGBS), chromatin- accessibility sequencing (ATAC-seq), chromatin immunoprecipitation sequencing (ChIP-seq) targeting five different histone modifications (H3K27me3, H3K4me1, H3K27ac, H3K4me3, and H3K36me3) and transcriptome sequencing (RNA-seq). We integrated previously published datasets from an extensively profiled human LCL (GM12878) to balance the number of human samples. b Number of sequencing reads generated per sample and experiment. Striped lines indicate data retrieved from previously published experiments79,108. We quantified the conservation of epigenetic states in In agreement with previous studies in more distant species27,28, regulatory regions as the number of primate species with the nearly all (93%) recently evolved promoter states are acquired epigenetic state in the orthologous regions. In the regulatory through repurposing events, whereas the majority (90%) of architectures of protein-coding genes, promoter states are more recently evolved enhancer states are gained (Chi-square test: P < conserved than enhancer states (Supplementary Fig. 9), with 73% 2.2 × 10−16; Methods and Supplementary Fig. 10). The regulatory and 60% of regions with a promoter or enhancer state being fully architectures of protein-coding and non-coding genes—the latter conserved across primates, respectively (two-tailed Fisher’s exact evaluated in human due to underrepresentation of non-coding test: P < 2.2 × 10−16, OR = 1.48). Less than 14% and 8% of gene annotations in nonhuman species—show this same pattern orthologous regulatory regions with a promoter or enhancer state (two-tailed Fisher’s exact test: P < 2.2 × 10−16, OR = Inf, and P = are specific to a primate species, respectively. These results for 6.6 × 10−7, OR = 137 respectively; Supplementary Fig. 10). protein-coding genes fall in line with the higher conservation of Our results confirm those found in more distant species13,27 promoters previously observed in mammals13. In contrast, for and reinforce the generality of these evolutionary dynamics in non-coding genes, promoter states are less conserved than protein-coding genes. The acquisition of regulatory states in the enhancer states (two-tailed Fisher’s exact test: P < 2.2 × 10−16, regulatory architectures of non-coding genes resembles that of OR = 0.39; Supplementary Fig. 9), with 46% and 69% of fully protein-coding genes. However, the lower conservation of conserved and 26% and 3% of species-specific elements, promoter states associated with non-coding genes suggests that respectively. their overall higher conservation is not an intrinsic characteristic Intrigued by the different epigenetic conservation patterns in of promoter states and that it depends on their specific regulatory protein-coding and non-coding genes, we studied the repurpos- relevance in different genes. ing and acquisition of RE. We defined recently repurposed promoters—or enhancers—as regulatory regions with a promoter state in only one species and enhancer states in the remaining The activity of promoters and enhancers influences their epi- species—or vice versa. Similarly, recently gained promoters or genetic and sequence conservation. Taking advantage of our enhancers are those regions with a promoter or enhancer state in classification of promoters and enhancers into three different one species and without regulatory states in any other species. activities (strong, poised, and weak), we further explored the NATURE COMMUNICATIONS | (2021)12:3116 | https://doi.org/10.1038/s41467-021-23397-1 | www.nature.com/naturecommunications 3 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-23397-1 patterns of evolutionary conservation of the different regulatory activities show a similar pattern of evolutionary conservation states. Globally, orthologous regulatory regions conserve their (Fig. 3a). Enhancer states with strong activities are second in regulatory state (Randomization analyses: 1000 simulations, P < conservation after strong promoters. Nearly 40% of the 0.05; Supplementary Table 7), but different promoter and orthologous regulatory regions with strong enhancer states are enhancer activities show characteristic patterns of conservation fully conserved. Poised enhancers follow closely, with 36% of (Kruskal–Wallis test: P < 2.2 × 10−16; Fig. 3a and Supplementary them conserved in the five species. Lastly, around 21% of the Fig. 11). regions with a weak enhancer conserve their activity across Strong promoters are the most conserved activities: 80% of primates. The regulatory regions associated with protein-coding them are fully conserved in primates. On the contrary, poised and and non-coding genes show the same conservation trends weak promoters are poorly conserved (Fig. 3a). All enhancer (Supplementary Fig. 12). However, strong activities in promoter 4 NATURE COMMUNICATIONS | (2021)12:3116 | https://doi.org/10.1038/s41467-021-23397-1 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-23397-1 ARTICLE Fig. 2 Epigenetic and regulatory characterization of RE annotated in primates. a Approach followed to annotate and classify RE. In short, we classify RE in promoter and enhancer states with three activity levels (strong, poised, or weak) based on a combination of chromatin marks and ATAC-seq signals. Bars represent relative enrichment of epigenetic signals in RE, from left to right: H3K4me3 (dark blue), H3K4me (light blue), H3K27ac (green), H3K36me3 (yellow), H3K27me3 (red), and open chromatin (gray). Promoters have a positive H3K4me3/H3K4me1 ratio, whereas H3K4me1 is more abundant in enhancers. RE with strong activities are associated with high H3K27ac levels and poised RE have a robust enrichment in H3K27me3. RE are then linked to genes based on 1D gene proximity and 3D published chromatin maps for LCLs. RE not associated with any gene are referred to as orphan RE. See Methods and extended representation in Supplementary Fig. 1. b Number of RE with promoter and enhancer epigenetic states in each species. aP and aE refer to ambiguous promoters and enhancers, respectively. Ambiguous RE are defined as those with a consistent state but different activities between replicates. Dashed lines indicate the average number of RE with promoter and enhancer states annotated across species. c Number of RE associated with genes and orphan RE in each species. Genes are divided into one-to-one orthologous protein-coding (1–1 orth PC), protein-coding (PC), and non-protein- coding genes. n n Z Z n n Z Z n n Z Z Fig. 3 Different regulatory activities have different patterns of epigenetic and sequence conservation. a Bar plots show the average number of orthologous regulatory regions across species with the corresponding color-coded epigenetic state conserved in 1, 2, 3, 4, or 5 species. Points indicate average values and the error bars represent the s.d. (n = 5 species). b Distribution of the sequence conservation scores (calculated as Z-scores of the distribution of phastCons30way29 values for non-coding regions in the same TAD85; Methods) of human orthologous regulatory regions with different epigenetic states conserved in 1, 2, 3, 4, or 5 of our primate species. Box plots show medians and the first and third quartiles (the 25th and 75th percentiles), respectively. The upper and lower whiskers extend the largest and smallest value no further than 1.5 × IQR. states are less common for non-coding than for protein-coding in this cell type. These regions are expected to show different genes, leading to lower conservation of promoter compared to epigenetic states in other cell types, and so their evolutionary enhancer states. This shows that differences in activity composi- patterns might also be different. To investigate whether changes tion can lead to differences in the conservation of the regulatory in activity are likely to affect the epigenetic conservation of RE, we architectures. assessed the association between epigenetic and sequence The epigenetic states in a given cell type and their evolutionary conservation—which is cell type-independent. First, we observed conservation reflect the specific function of the regulatory regions that epigenetic conservation significantly correlates with the NATURE COMMUNICATIONS | (2021)12:3116 | https://doi.org/10.1038/s41467-021-23397-1 | www.nature.com/naturecommunications 5 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-23397-1 conservation of the underlying sequence—quantified as Z-scores enriched at gE and PiE, whereas weak enhancers are strongly of background normalized PhastCons values29—in all epigenetic associated with prE (Supplementary Fig. 19). states but weak promoter states (Fig. 3b, Methods, and Gene expression levels are positively associated with the presence Supplementary Figs. 13 and 14). These correlations are seen in of strong activities in their regulatory architectures and are the architectures of protein-coding but not in non-coding genes negatively associated with the presence of poised or weak activities (Randomization analyses: 1000 simulations; Fig. 3b and Supple- (Kruskal–Wallis test: P < 0.05 in all species and regulatory mentary Fig. 15). Of note, orthologous regulatory regions with components; Supplementary Fig. 20). These associations are fully conserved epigenetic states show significant differences in particularly strong in gP and gE (Dwass–Steel–Critchlow–Fligner sequence conservation (Kruskal–Wallis test: P < 2.2 × 10−16; test; Supplementary Fig. 21). Despite the consistency between the Supplementary Fig. 16). In particular, strong and weak promoters components’ activities and gene expression, our results suggest that are associated with higher and lower sequence conservation different types of components might contribute differently to gene scores, respectively, whereas all enhancer states range in between expression regulation. these values (Dwass–Steel–Critchlow–Fligner test). The sequence conservation scores associated with strong and poised enhancers are not significantly different. Orthologous regions associated Regulatory components influence gene expression and its with non-coding genes are fewer and less epigenetically conserved evolution differently. To explore the ability of our classification (Supplementary Figs. 12 and 14), which could explain the lack of of components to discriminate different regulatory roles, we correlation between the conservation of the sequence and the disentangled the underlying network of regulatory co- epigenetic state observed in all but strong enhancers (Supple- dependencies between the different regulatory components and mentary Fig. 15). gene expression in our cell-type. For this, we used Sparse Partial These results demonstrate that a detailed classification of Correlation Analysis (SPCA)34 of the normalized RNA-seq and promoters and enhancers with different activities into regulatory histone mark enrichments (aggregated by promoter and enhancer architectures provides a deeper understanding of their evolu- state in every type of regulatory component) (Methods). This tionary constraints and dynamics, expanding previous observa- approach establishes a stringent protocol (Benjamini–Hochberg’s tions in mammals13 that could mostly be made for active correction, P < 1.8 × 10−22 for all selected partial correlations) regulatory activities. The consistent association of epigenetic and that selects informative partial correlations34. sequence conservation also suggests that the epigenetic conserva- To unravel the contribution of each type of component to gene tion observed in LCLs is a good proxy for the conservation of the expression, we defined their consensus signal (or eigencompo- regulatory activity of these elements in our primate species. nents) inspired by the notion of eigengenes35 (Methods). An SPCA based on the eigencomponents shows a consistent global structure of regulatory interactions. gP and gE directly regulate Definition of different types of components in the regulatory gene expression coordinately, PiE are connected with gPs and EiE architectures. To characterize the evolution of RE based on their with PiE (Fig. 4e and Supplementary Data 5). This regulatory specific role in gene expression, we classified RE into five different scaffold is consistently observed for the residuals of the histone components according to their role in the gene regulatory marks for these eigencomponents (Supplementary Fig. 22 and architectures (Fig. 4a, Methods). We first classified RE based on Supplementary Data 6) when SPCA was performed for all the their proximity to a gene into three types of components: genic histone marks together (Supplementary Data 7) and for each of promoters (gP), gE, and proximal enhancers (prE). As gene them separately (Supplementary Data 8). To account for the expression is controlled by a combination of short- and long- possibility of incompleteness in some of our architectures, we distance regulatory interactions30, we used available 3D chro- replicated all the analyses using only genes with full regulatory matin contact maps for human LCLs31–33 to link interacting RE architectures (i.e., genes associated with regulatory components of to their target gene/s and define two additional types of compo- every type), obtaining consistent results (Supplementary Fig. 23 nents: promoter-interacting enhancers (PiE) and enhancer- and Supplementary Data 8). interacting enhancers (EiE) (Fig. 4a). In agreement with the structure of regulatory interactions We were able to link to genes and classify, on average, nearly recovered by our SPCAs, a generalized linear model of gene 3500 otherwise orphan distal RE per species (Supplementary expression based on H3K27ac, H3K27me3, and H3K36me3 signals Fig. 17). We annotated ~12,500 gP, ~35,000 gE, ~6700 prE, ~6200 at gP and gE and their interactions (15 variables) explains ~67% of PiE, and ~1800 EiE per species (Fig. 4b and Supplementary gene expression variability (Supplementary Table 8). Remarkably, Fig. 18), of which 48%, 69%, 40%, 62%, and 61% are associated this is only 6% lower than an exhaustive naive model, including the with one-to-one orthologous protein-coding genes in all primate signal from all histone marks at all types of regulatory components species (Fig. 4c). with all possible interactions (1225 variables) (Supplementary To assess the consistency of our classification of regulatory Data 9). These results suggest the epigenetic activities of putatively components, we focused on one-to-one orthologous protein- strong and poised gP and gE and their interactions likely have a coding genes considering all their associated RE (i.e., 6 epigenetic large influence on gene expression regulation in our regulatory states x 5 components = 30 regulatory subcategories). We found architectures. However, their co-dependency with the other high concordance between the epigenetic state (based on ChIP- components suggests that they are dependent, in turn, on the seq and ATAC-seq data, Fig. 2a) and the component (based on coordination of the whole architecture. Although we cannot infer the type of association with the gene, Fig. 4a) of the RE. On causality from our SPCA analysis, these networks reflect that average, 75% of gP have a promoter state, and 90% of gene- regulatory co-dependencies between components depend on the associated enhancers have an enhancer state (two-tailed Fisher’s distance of the elements in the network of chromatin contacts (with exact test: P < 2.2 × 10−16 in all species, average OR = 64; gP and gE being in the gene locus, PiE interacting directly, and EiE Supplementary Fig. 19). This concordance is also consistent interacting indirectly with it). The robustness of these networks of across species (Chi-square test: P < 2.2 × 10−16 in all species; direct co-dependencies, their ability to explain gene expression, and Fig. 4d). gP are enriched in RE with strong promoter and poised their correspondence with the spatial disposition of the elements promoter and enhancer states. Strong enhancers are mostly show that these components reflect specific regulatory roles. 6 NATURE COMMUNICATIONS | (2021)12:3116 | https://doi.org/10.1038/s41467-021-23397-1 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-23397-1 ARTICLE Fig. 4 Epigenetic signals in gene regulatory architectures explain gene expression levels. a Classification of RE according to their regulatory roles in gene architectures. Of note, EiE may interact with any type of enhancer in a regulatory architecture (prE, gE, PiE, and EiE). b Average number of RE across species associated with one-to-one orthologous protein-coding genes classified as gP, gE, prE, PiE, and EiE. Differently shaped dots show the number corresponding to each species. Points indicate average values and the error bars represent the s.d. (n = 5 species). c Average number of orthologous protein-coding genes associated with each type of regulatory element. Differently shaped dots show the number corresponding to each species. Points indicate average values and the error bars represent the s.d. (n = 5 species). d Proportion of RE with a given epigenetic state associated with one-to-one orthologous protein-coding genes for each type of regulatory component. Dots and error bars show the average proportion and s.d. across species (n = 5 species), respectively. e Sparse Partial Correlation Networks showing the statistical co-dependence of the RNA-seq (Gene expression) and the consensus ChIP-seq signals for the five histone marks in every component represented by the eigencomponents (minimal partial correlation = −0.41; maximal partial correlation = 0.33; all partial correlations Benjamini–Hochberg’s P < 4.1 × 10−303). Edge widths are proportional to absolute partial correlation values within each network. The networks are based on the 5737 one-to-one orthologous protein-coding genes associated with at least one regulatory element in all species. Only nodes for values with significant and relevant partial correlations are represented. Previous studies have found that gene expression evolution is changes that alter the epigenetic state at gP, specifically the associated with changes in the regulatory complexity of a gene presence of either a strong promoter or poised enhancer, as well (the number of nearby RE)10. Since we could classify the RE of a as the number of gE with either strong or poised enhancer states, gene into different components (Supplementary Fig. 24), we were show the most robust associations with gene expression able to investigate the association of gene expression changes differences. The number of prE in any enhancer epigenetic state (Supplementary Fig. 25) with the evolutionary differences in the and strong promoters and strong and poised enhancers in PiE complexity of each type of component. We found that the effect also show significant though modest effects. These results of changes in complexity on gene expression levels depends on highlight that the additive nature of gene regulation depends on the epigenetic state gained or lost and the type of regulatory regulatory architectures. This dependency can be captured either component affected (Supplementary Fig. 25). Evolutionary by the aggregation of histone enrichment signals (as in our NATURE COMMUNICATIONS | (2021)12:3116 | https://doi.org/10.1038/s41467-021-23397-1 | www.nature.com/naturecommunications 7 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-23397-1 SPCAs) or by quantifying the number of regulatory components and Supplementary Data 16). Similar to genes associated with with specific activities. Moreover, they confirm that our their analogous conserved group (conserved weak enhancers regulatory components represent different regulatory roles with states in gE), these genes are typically expressed at low levels a different contribution to gene expression evolution and which with their highest expression in tissues unrelated to LCLs, evolutionary relevance should be investigated separately. particularly brain, tibial nerve, and testis, while having marginal or no expression in numerous other tissues, including LCLs (Wilcoxon–Nemenyi–McDonald–Thompson test: P < 1 × 10−4; Different regulatory components with poised and weak Rank-biserial correlation effect size between brain and LCLs = enhancer states appear in LCL-unrelated genes and functions. 0.633; Fig. 5b, c, Supplementary Figs. 28 and 29, and We next explored the functional implications that both conserved Supplementary Data 17). Remarkably, these genes have higher and species-specific regulatory states in certain components could tissue-specific expression than those with conserved strong have (overrepresented combinations in Supplementary Fig. 26). activities in their components (median τ = 0.84, For this, we examined their functional role and tissue-specific Dwass–Steel–Critchlow–Fligner test: P < 4.5 × 10−14; Supple- expression (GTEx data36, Supplementary Table 9). We find sig- mentary Fig. 31). nificant enrichments for the genes targeted by conserved strong Intrigued by the high tissue-specificity of the genes with promoter states in gP, conserved strong and weak enhancer states hswEgE, we sought to identify the tissues driving this tissue- in gE, and conserved poised enhancer states in gP and prE (one- specificity taking its analogous conserved group as reference. tailed Fisher’s exact test: Benjamini–Hochberg’s correction, false Testis and brain are the tissues with the highest number of tissue- discovery rate (FDR) < 0.05; Methods, Fig. 5a, Supplementary specific genes (τTissue > 0.8), but most interestingly, whereas the Fig. 27 and Supplementary Data 10). Remarkably, among the fraction of testis-specific genes is comparable between gene sets genes associated with any species-specific epigenetic state, only (two-tailed Fisher’s exact test: P > 0.05, OR = 1.20), brain-specific those linked to human-specific weak enhancers in intragenic genes are more than twofold enriched in genes with human- enhancers (hereafter referred to as hswEgE) have significant specific gE (two-tailed Fisher’s exact test: P = 0.02, OR = 2.29; functional enrichments. Supplementary Fig. 31). These enrichments show the expected association of conserved Our results show that conserved strong epigenetic states are strong epigenetic states (strong promoter states in gP and strong involved in the regulation of important genes highly expressed in enhancer states in gE) with genes involved in relevant cellular LCLs and other tissues. However, poised and weak enhancer processes, such as metabolism, chromatin organization, and states, either conserved or human-specific, are involved in the regulation of the cell cycle (Fig. 5a, Methods, Supplementary regulation of genes marginally expressed in LCLs, but with Fig. 27, and Supplementary Data 11 and 12). Functions specific to particular functional roles and tissue-specific expression patterns. LCLs like those involving viral processes are specifically enriched These unexpected associations, along with the evolutionary in strong enhancers30. Moreover and regardless of whether there conservation patterns of poised and weak enhancer states, suggest is a functional enrichment, the components with strong that these epigenetic states might be indicative of putative RE in epigenetic states show similar expression profiles across tissues other cell types different from LCLs. (Fig. 5b and Supplementary Figs. 28–30), with high expression levels in LCLs and most other tissues and wide expression breadth (Fig. 5c and Supplementary Fig. 31). Genes with novel human-specific intragenic weak enhancers In contrast, conserved poised enhancer states in gP and prE are targeted by signals of adaptation and nucleotide changes. target genes enriched in developmental and proliferative func- The unexpected association of the genes targeted by hswEgE with tions, echoing their known implication in these processes37–40 neuronal functions prompted us to study the relationship these (Fig. 5a, Supplementary Fig. 27, and Supplementary Data 13 and genes might have with signals of positive selection or accelerated 14). Surprisingly, genes with conserved poised enhancer states in evolution in the human genome. In fact, among the genes asso- their gP are also enriched in neuronal functions and have higher ciated with hswEgE (153 hswEgE in 134 genes, Supplementary expression levels in brain (Fig. 5b, c and Supplementary Figs. 27– Data 18), we found several genes previously proposed to have 29). Genes associated with both types of conserved poised been subjected to different evolutionary forces. Some of these enhancer states show overall minimal expression levels but high genes are FOXP2, PALMD, and ROBO1, which have known tissue-specificity (median tissue specificity index (τ, Tau) > 0.85 in brain-related functions42–45 or ADAM1846, CFTR47,48, and both; Methods and Supplementary Fig. 31). TBX1549. Protein-coding genes targeted by gE with conserved weak To assess whether genes with hswEgE have been particularly enhancer states are enriched in various functional annotations, targeted by evolutionary forces, we investigated their co- including neuronal ones, such as cell projection and synapse occurrence in genes associated with signals of positive (Fig. 5a, Supplementary Fig. 27 and Supplementary Data 15). selection50–54 and acceleration55–57, hereafter referred to as This gene set shows low expression in LCLs and high expression putatively selected or accelerated regions (Methods). We find in the brain, which is in agreement with its observed that 41% of the genes with hswEgE (56 genes) are also showing functional annotation (Fig. 5b, c). Also, the tissue-specificity of such signatures (one-tailed Fisher’s exact test: P = 3.33 × 10−11, this group is higher than that of conserved strong regulatory OR = 3.5). The results of our nested analysis (Fig. 6a) show a activities both in promoters and enhancers (median τ = 0.72, significant association of genes targeted by gE (but not in any Dwass–Steel–Critchlow–Fligner test: P < 2.2 × 10−16 in the three other component type), genes targeted by gE with weak enhancer tests; Supplementary Fig. 31). This apparent brain-specificity is states (but not strong or poised enhancer states), and genes not found in genes associated with conserved weak enhancer targeted by gE with human-specific weak enhancer states (but not states in other regulatory components, which have overall higher with fully conserved or any other epigenetic states in nonhuman expression levels and not particular tissue-specificity, as would be primate species). In addition, no enrichment is observed for genes expected from weak epigenetic states41 (Supplementary Figs. 30 with gP with conserved poised enhancer states (one-tailed and 31). Fisher’s exact test: P > 0.05, OR = 1.1), even though they also Finally, we focused on genes targeted by hswEgE. These seem to be associated with brain-specific and neuronal functions genes are solely enriched in neuron parts and synapse (Fig. 5a (Fig. 5). 8 NATURE COMMUNICATIONS | (2021)12:3116 | https://doi.org/10.1038/s41467-021-23397-1 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-23397-1 ARTICLE n n n n n Fig. 5 Functional enrichment and gene expression patterns indicate that weak and poised enhancer states might be associated with LCL-unrelated cell types. a Functional enrichment of conserved and human-specific activities in genic promoters and gE. The size of circles indicates the proportion of genes included in each functional category from the total number of genes contained in the corresponding regulatory group. Extended version in Supplementary Fig. 27. b Heatmap of standardized expression across tissues in state/component regulatory groups with functional enrichments. Extended version in Supplementary Fig. 28. c Median expression levels in testes, LCLs, brain, and whole blood of genes groups in b. Box plots show medians and the first and third quartiles (the 25th and 75th percentiles), respectively. The upper and lower whiskers extend the largest and smallest value no further than 1.5 × IQR. Finally, we explored whether hswEgE could also be associated humans are not enriched in W to S mutations compared to the with human-specific sequence changes. For this, we collected a set ancestral state (one-tailed Fisher’s exact test: P > 0.05, OR = 0.9). of over 2.8 million single nucleotide changes fixed in humans Taken together, roughly 10% of the hswEgE-containing genes (hSNCs) that differ from fixed variants in the genomes of the are in putatively selected or accelerated regions and also harbor remaining nonhuman primates (Methods). Around 22% of the hSNCs in these enhancers (14 genes) and the co-appearance of hswEgE-containing genes (30 genes) harbor at least one hSNC in both features is not significantly enriched (one-tailed Fisher’s exact their hswEgE. Also, the density of these human-specific changes test: P > 0.05, OR = 1.3, Fig. 6b). This observation indicates that is higher in hswEgE compared to gE with conserved weak even though putatively selected or accelerated regions and human- enhancer states (Mann–Whitney U-test: P = 0.01; Methods and specific mutations are both associated with hswEgE, they are not Supplementary Fig. 32). This observation is not driven by GC- mutually conditioned. As such, it implies that none of these biased gene conversion (gBGC), as the derived nucleotides in signals alone is required to explain the appearance of hswEgE. NATURE COMMUNICATIONS | (2021)12:3116 | https://doi.org/10.1038/s41467-021-23397-1 | www.nature.com/naturecommunications 9 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-23397-1 Fig. 6 Human-specific intragenic enhancers with weak activities (hswEgE) co-localize with signals of recent human selection and acceleration and human-specific nucleotide changes. a Diagram with the results of the nested enrichment analysis. hswEgE-containing genes are enriched in genes targeted by putative signals of positive selection and acceleration. Epigenetic states for each species and consensus regulatory components are depicted as boxes with the same color code as Figs. 2b and 4a, respectively. b Top: Schematic representation of a hswEgE with an hSNC (nucleotide change in humans shown in red) in a gene with signals of selection and acceleration. Bottom: Venn diagram illustrating the overlap between the 56 genes containing hswEgE and with signals of selection or acceleration and the 30 genes with hswEgE and hSNCs. Among the 14 genes with both human-evolving signals and nucleotide changes are especially concentrated in three of them. hSNCs (Fig. 6b), there are several interesting candidates for adaptive Two of these genes, CLVS1 and ROBO1, exemplify the confluence evolution associated with different traits. Many of these genes are of signals of human adaptation and hSNCs in protein-coding associated with neuronal functions (ROBO1, CLVS1, SEMA5A, genes that are important for normal neuronal structure, KCNH7, SDK1, and ADGRL2), but also with pigmentation migration, and axon guidance in the human brain. (LRMDA) or actin organization in cardiomyocytes (FHOD3). Other interesting genes include instances of hswEgE-containing Discussion genes that also are targeted by putatively selected or accelerated The evolution of human and nonhuman primates is an area of regions (FOXP2, ASTN2, NPAS3, or NTM) or carry hSNCs in these major interest, but ethical, legal, and practical constraints often enhancers (PALMD, VPS13C, IGSF21, or CADM2). Interestingly, limit access to direct biological material. In this study, we have we found only one antisense RNA gene, MEF2C-AS1, showing generated a comprehensive and unified dataset of epigenomic both. This gene has been associated with ADHD58, and its target landscapes in LCLs for human and four nonhuman primate gene MEF2C is a very well-known target of genetic alterations species. Despite the artificial nature of our cellular model65–67, (many of them also affecting MEF2C-AS1) associated with severe previous studies have shown the value of LCLs as an experi- intellectual disability59, cerebral malformation59, or depression59,60. mentally convenient model of somatic cells that accurately Remarkably, three hswEgE accumulate more hSNCs than resembles the phenotype of its cell type of origin68 and which can expected (Randomization test: 10,000 simulations, Bonferroni be robustly used for comparative studies in humans and correction, P < 0.02 in all cases; Methods and Supplementary primates18,69–71. Moreover, its clonality ensures a cell type- Fig. 32), a number of hswEgE which is also significantly higher specific experimental system reducing the confounding factors than expected (Randomization test: 10,000 simulations, P = 8 × associated with cell population diversity in bulk tissue samples. 10−4; Supplementary Fig. 32). Two of these genes are protein- Using this cell model, we reproduced previous observations on coding genes with known functions in brain cell types and with the dynamics of the evolution of RE reported in more distant signals of human adaptation. CLVS1 is a protein-coding gene species10,13,27, which we show can be extrapolated to closely with brain-specific expression (τBrain = 0.964) required for the related species (at least for great apes and macaques). Moreover, normal morphology of endosomes and lysosomes in neurons61. we have expanded these observations to explain how these ROBO1 is a broadly expressed integral membrane protein that dynamics result from the different evolutionary constraints participates in axon guidance and neuronal migration (τ = associated with their epigenetic activities. Therefore, we show that 0.388)62,63 that has also been associated with human speech and considering weak and poised activities is of major relevance to language acquisition since the split from chimpanzees43. The fully understand the evolution of regulatory regions. third enhancer is included in AC005906.2, a long intergenic non- We also observed that different epigenetic activities have protein-coding gene specifically expressed in brain (τBrain = 1). characteristic evolutionary patterns with higher conservation for Interestingly, this gene overlaps with KCNA1, a voltage-gated strong promoter and strong and poised enhancer states. The potassium channel with the same brain-specific expression correlations between epigenetic and sequence conservations are pattern (τBrain = 0.995) and for which mutations have been also different for each epigenetic state, with higher correlations associated with neurological malfunctions64. for strong and poised promoter and enhancer states. These dif- We show that hswEgE target neuron-related brain-specific ferences are likely due to their different influence on gene genes that are associated with putative signals of positive selection expression. Therefore, the previously reported higher conserva- and acceleration and an excess of hSNCs. These human-specific tion of promoters probably reflects its often greater influence on 10 NATURE COMMUNICATIONS | (2021)12:3116 | https://doi.org/10.1038/s41467-021-23397-1 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-23397-1 ARTICLE gene expression. We are also able to confirm this hypothesis with marginally expressed and unrelated to LCLs, we hypothesize that the lower conservation of promoter states in the regulatory functionally relevant regulatory innovations that appear in a architectures of the non-coding genes where strong promoter given cell type could be echoed as weaker activity signals in a states are scarce. different cell type. This could explain our observations in LCLs. We further classified the RE in five different components of As such, our results suggest that human-specific weak enhancer gene regulatory architectures. These components are genic pro- activities could provide an unexpected window for the study of moters (gP) and intragenic, proximal, promoter-interacting, regulatory evolution in the human lineage. Further research will and enhancer-interacting enhancers (gE, prE, PiE and EiE, be needed to clarify the specific role of these elements in different respectively). By studying their regulatory co-dependencies, we tissues and cell types. observe that the epigenetic activity of each type of component Taken together, this work shows that the evolution of gene influences gene expression differently. In brief, coordinated epi- regulation is deeply influenced by the coordination of epigenetic genetic activities in gP and gE form the core of these architectures activities in gene regulatory architectures. Our insights call for the and are strongly correlated with gene expression levels. Reg- incorporation of better integrative datasets and refined definitions ulatory activities in PiE are also coordinated with promoter of regulatory architectures in comparative evolutionary studies to components, and activities in EiE are associated with PiE. With fully understand the interplay between epigenetic regulation and these results, we show that the impact of regulatory components gene expression. on gene expression reflects the structure of the regulatory architecture. Methods We find that gene expression changes are associated with the Cell line acquisition and cell growth. Lymphoblastoid cell line (LCL) GM19150 regulatory complexity of gene architectures and with the epige- (Yoruban male) was purchased from the Coriell Institute. Chimpanzee, gorilla, and netic activity of the different regulatory components. Namely, we orangutan LCLs75 were kindly provided by Dr. Antoine Blancher. Macaque LCLs76,77 were kindly provided by Dr. Gaby Dioxiadis. report that gene expression changes between primate species Cell lines were grown in suspension at a confluency of 200,000 cells/mL to 1 × occur through the addition or removal of strong promoter 106 cells/mL, at 37 °C and 5% CO2, in RPMI 1640 media (Invitrogen 42401-018) activities in promoter components or strong and poised enhancer with 15% fetal bovine serum (FBS Invitrogen 10270-106), 1% penicillin- activities in gE. The remaining components show fewer changes streptomycin (Invitrogen 15140-122), and supplemented with 2 mM L-glutamine (Sigma G7513-100ML). directly linked to expression differences, but they can still be instrumental for gene expression evolution, probably through DNA and RNA extractions. For DNA and RNA extractions, ~5 million cells were their influence on gP and gE. The gene architectures we present pelleted and washed with PBS (Sigma D8537). provide a starting point for future in-depth investigations on the DNA extractions were performed using a phenol-chloroform protocol. Briefly, interdependence of different regulatory regions and mechanisms cell pellets were resuspended in a lysis buffer (10 mM Tris-HCl pH 8, Life in the evolution of gene regulation. In this sense, we stress the Technologies 15568-025; 10 mM NaCl, Sigma 71386; 0.2% IGEPAL CA-630, Sigma I8896) and incubated in ice for 30 min. Next, cells were washed and incubated at importance of considering promoter and enhancer activity states 37 °C for 30 min in NEB2 buffer (New England Biolabs B7002S), supplemented in the different types of gene components to achieve a more with 10% SDS (Life Technologies AM9820) and 5 μL of RNAse A (QIAGEN detailed description of the regulatory processes. 79254), followed by an overnight incubation at 65 °C with 10 μL Proteinase K (Life Despite the large influence of strong activities on gene reg- Technologies AM2546). Two successive phenol-chloroform extractions were performed the day after. DNA was precipitated with NaAc 3 M and washed twice ulation, our results in LCLs suggest that unforeseen information with 70 and 100% EtOH, respectively. DNA integrity was checked in an agarose can be drawn from the analysis of elements with a repressive or gel. Extracted DNA was used for WGS and WGBS libraries. Libraries were negligible regulatory role in our cell model. Functional enrich- prepared using the TruSeq DNA PCR Free Library Preparation Kit and TruSeq ment analyses show that genes with conserved poised activities DNA Methylation Kit, respectively, following Illumina’s standard protocol. One- are enriched in cell proliferation and differentiation, as shown by fifty base pair paired-end (150PE) reads were sequenced in a HiSeqX machine. RNA was extracted using the miRNeasy Mini kit (QIAGEN 217004) following others37–40. In the case of poised enhancer states in gP, these the manufacturer’s instructions. The TruSeq Stranded Total LT Samples Prep Kit genes are specifically expressed in brain. Moreover, we find that with Ribo-Zero was used for library construction, following Illumina’s standard recently evolved weak gE in the human lineage (hswEgE) occur in protocol. Library preparation and sequencing were performed at Macrogen (Seoul, genes showing patterns of brain-specific gene expression and Korea). Three technical replicates (R1, R2, and R3) were sequenced per sample. All nine libraries per technical replicate were pooled and loaded into one lane. About enriched in neuronal functions. This tissue-specific activity is 101 bp paired-end reads (101PE) were sequenced in a HiSeq4000 machine. consistent with what has been recently reported for gE72. These genes are also associated with signals of positive selection and Chromatin immunoprecipitation. For ChIP-seq experiments, LCL primate cell acceleration, suggesting they may have contributed to human lines were harvested and washed with PBS (Sigma D8537). Cells were fixed by adaptation in several traits. However, we acknowledge that these incubating in 1.5% formaldehyde (Sigma F8775) PBS buffer at room temperature evolutionary signals could also be explained by relaxed for 10 min. Formaldehyde was quenched with 125 mM glycine (Invitrogen 15527- 013) at room temperature for 5 min, after which the samples were placed on ice for selection73. These evolutionary patterns are unique to hswEgE 15 min. Cells were then washed twice in cold PBS, aliquoted, and pellets of ~20 and not found in conserved poised enhancers in gP even though million cells were frozen at −80 °C. ChIP experiments were carried out following a these elements also target genes with brain-related functions. standard protocol in the lab78. Briefly, crosslinked cell pellets were resuspended in Moreover, we show that hswEgE accumulate more human- 500 μL ice-cold ChIP buffer (50 mM Tris-HCl pH 8, 100 mM NaCl, 5 mM EDTA, specific sequence changes than expected. These changes could 0.33% SDS, 1.67% Triton X-100) with Protease inhibitors (cOmpleteTM EDTA- free, SIGMA) and 1 mM PMSF. Lysates were sonicated in a Bioruptor (Diagenode) potentially alter the activity of the enhancer, as has been recently to obtain an average fragment size of 250 bp and centrifuged for 20 min at full shown74. We highlight two interesting protein-coding genes that speed. An aliquot of the soluble chromatin was taken, reverse cross-linked, purified show the convergence of both signals and an excess of human- using the Qiagen PCR purification kit (Cat No./ID: 28104), and quantified by specific mutations which have key functional roles in neuronal Nanodrop. For each ChIP, 30 μg of soluble chromatin, 0.75 μg of mouse E14TG2a (as spike-in control), and 5 μg of the corresponding histone modification antibody structure, migration, or axon guidance: CLVS1 and ROBO1. were incubated overnight at 4 °C. ChIP DNA was recovered by incubating with 42 Using LCLs as our model system, we have seen that poised and μL protein A bead slurry (Diagenode) for 2 h, washed three times with low salt weak enhancer activities may carry information about gene reg- buffer (50 mM HEPES pH 7.5, 140 mM NaCl, 1% Triton X-100) and one time with ulation in unrelated cell types. Notably, we find that genes in high salt buffer (50 mM HEPES pH 7.5, 500 mM NaCl, 1% Triton X-100). DNA complexes were decross-linked at 65 °C overnight with proteinase K (Invitrogen) putatively selected or accelerated regions also harbor human- and DNA was purified using the PCR purification kit (Qiagen, 28104) and specific epigenetic and sequence innovations. Given the tissue- quantified by Nanodrop. We used commercially available antibodies against specificity and functional enrichments of these genes, which are H3K4me3 (Diagenode, C15310003), H3K4me1 (Abcam, ab8895), H3K36me3 NATURE COMMUNICATIONS | (2021)12:3116 | https://doi.org/10.1038/s41467-021-23397-1 | www.nature.com/naturecommunications 11 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-23397-1 (Abcam, ab9050), H3K27me3 (Millipore, 07-449), and H3K27ac (Millipore 07- misclassifications are expected due to an insufficient precision of the qualitative 360). classification. Considering the quantitative relationship between co-existing histone A total of 63 ChIP-seq libraries were constructed: 45 libraries for the ChIP-seq modifications can help to accurately annotate epigenetic states in RE. We used experiments (one per sample and histone modification) and 18 libraries for Input linear discriminant analysis (LDA)80 to refine chromatin-state based annotations. samples (Input1 for histones H3K4me3, H3K4me1, H3K36me3, and H3K27me3; This method is commonly applied to pattern recognition and category prediction. Input2 for H3K27ac). Prior to library preparation, ChIP-qPCR validations were LDA is a technique developed to transform the features into a lower-dimensional performed using primer pairs of known active and repressed genes in human LCLs space, which maximizes the ratio of between-class variance to the within-class (active: CEP250, GAPDH, Beta-Actin; repressed: Sox2) and a negative intergenic variance, thereby granting maximum class separation. We performed LDA analysis region as a control for the histone modifications studied. About 2–10 ng of ChIP using the lda function in the R package MASS (version 7.3-47)81. The predictor DNA were used to prepare the sequencing libraries using the NEBNext Ultra DNA variables were the background-noise normalized IP signals from the five different Library Prep Kit for Illumina (NEB, E7370L) as per the manufacturer’s histone modifications profiled and the chromatin accessibility signal at species RE. instructions. ChIP-seq libraries were size selected to remove fragments below 100 The categorical variable to be predicted based on the underlying enrichments was bp and amplified for 10 PCR cycles. Four to five libraries were pooled and loaded the chromatin-state-based annotation. The regulatory state at the species level was into a total of 13 lanes. About 50 bp single-end reads (50SE) were sequenced in an determined based on the regulatory state in each of the biological replicates. Thus, Illumina HiSeq2500 machine, using v4 chemistry. the regulatory state of an element with different epigenetic states in the two replicates (ambiguous) could be aP or aE, when both samples of a given species were annotated as either P or E but differ in their activity; P/E, when a RE was Assay of transposable chromatin. ATAC-seq libraries for nine lymphoblast cell classified as P in one biological replicate and E in the other one; and P/Non-RE or lines were generated following the Buenrostro et al. protocol79. Briefly, 50,000 cells E/Non-RE, when the RE was detected only in one replicate (Supplementary Fig. 6 were harvested, washed in cold PBS and resuspended in 50 mL of cold lysis buffer and Supplementary Data 1). Larger inter-sample variability is expected between (10 mM Tris-HCl pH 7.4, 10 mM NaCl, 3 mM MgCl2, 0.1% (v/v) Igepal CA-630). human replicates as H1 data were obtained from previously published datasets. Samples were spun down for 10 min at 500 × g, 4 °C. Pellets were resuspended in Hence, to control for interindividual variability, only RE with the same activity in the transposition reaction mix and incubated at 37 °C for 30 min. Samples were the two replicates were considered for downstream analyses. purified using the Qiagen MinElute PCR purification kit. Transposed DNA was eluted in 10 μL of elution buffer and subjected to PCR amplification for eight cycles using barcoded primers and NEBNext High Fidelity PCR master mix. ATAC-seq Analysis of evolutionary conservation at orthologous regulatory regions. This libraries were purified using 1.8X volumes of AMPure XP beads to remove frag- analysis was performed to explore the conservation patterns of RE across species. ments below 100 bp. Library quality was assessed using a Bioanalyser High Sen- We studied the evolutionary conservation of epigenetic signals at two levels. First, sitivity DNA analysis kit (Agilent). All nine libraries were pooled and loaded in a we investigated the conservation of RE with promoter or enhancer epigenetic total of three lanes. 50 bp paired-end reads (50PE) were sequenced on a HiSeq 2500 states, regardless of their activity. Second, we investigated the evolutionary con- platform (Illumina), using v4 chemistry. servation of strong, poised, and weak promoters and enhancers. To study the evolutionary conservation of promoters and enhancers, we focused Definition of RE. This analysis was performed to identify putative RE in each on a set of 21,753 one-to-one orthologous regions associated with genes, in which primate species. We used ChromHMM to jointly learn chromatin states across at least one species showed a promoter or enhancer epigenetic state. To measure samples and segment the genome of each sample26. ChromHMM implements a the evolutionary conservation at each orthologous regulatory region, we count the multivariate Hidden Markov Model aiming to summarize the combinatorial number of species (1, 2, 3, 4, or 5) in which the promoter or enhancer state is interactions between multiple chromatin datasets. Bam files from the five histone conserved (Supplementary Fig. 9). modifications profiled were binarized into 200 bp density maps. Each bin was To study how evolutionary changes occur, we define recently repurposed discretized in two levels, 0 or 1, depending on their enrichment computed by promoters as orthologous regulatory regions in which one species shows a comparing immunoprecipitated (IP) versus background noise (input) signal within promoter state while the others show an enhancer state or vice versa. Similarly, we each bin and using a Poisson distribution. Binarization was performed using the refer as novel promoter or enhancer states to those orthologous regulatory regions BinarizeBam function of the ChromHMM software26. A common model across in which a given species showed a promoter or enhancer state while the others species was learned with the LearnModel ChromHMM function for the con- showed no evidence of regulatory activity (classified as nonregulatory). catenated genomes of all samples but O1 (orangutan sample 1) due to its anom- To study the evolutionary conservation of strong, poised, and weak promoters alous methylation profile (Supplementary Fig. 43). Several models were trained and enhancers, we focused on the subset of 10,641 one-to-one orthologous regions with a number of chromatin states ranging from 8 to 20. To evaluate the different in which at least one species showed a strong, poised, or weak regulatory state (we n-state models, for every sample, the overlap and neighborhood enrichments of do not include orthologous regions including elements with ambiguities, i.e., each state in a series of functional annotations were explored. A 16-state model was different activities between biological replicates). To measure the evolutionary selected for further analysis based on the resolution provided by the defined conservation at each orthologous regulatory region, we count the number of chromatin states, which capture the most significant interactions between histone species (1, 2, 3, 4, or 5) in which each regulatory state (strong, poised, and weak marks and the state enrichments in function-annotated datasets (Supplementary promoters or enhancers) is conserved (Supplementary Fig. 12). To statistically Fig. 2). The genomic coordinates of RE were defined for each sample by merging all assess the different evolutionary dynamics for the different regulatory states, we consecutive 200 bp bins excluding elongating (E1 and E2), repressed hetero- first ran randomization analyses. We randomized (1000 randomizations) the chromatin (E16), and low signal (E15) chromatin states. Species RE were defined as regulatory states associated with each species in the orthologous regulatory regions. the union of sample RE. For orangutan, we did not include RE specific to O1. We calculate the P value as the number of randomizations with average conservation equal to or above the observed conservation for each regulatory state. We further explored the different patterns of conservation combining: (1) Assignment of a regulatory state to RE. This analysis was performed to assign Kruskal–Wallis test (kruskal.test R function)82 to test whether the global each RE a regulatory state based on their epigenetic signals. We did so using a two- distributions of the number of species in which each particular state was conserved step approach. First, we classified RE based on the chromatin states found in them. were different for the different regulatory states, (2) Second, we refined this classification based on their underlying histone modifica- Dwass–Steel–Critchlow–Fligner test to assess the significance of every pairwise tions and open chromatin signals. comparison (dscfAllPairsTest function from the R package PMCMRplus version RE were assigned a chromatin-state based annotation. Combining the 1.4.4)83, and (3) Glass rank biserial correlation coefficient for Mann–Whitney U- information gathered through the overlap and neighborhood enrichment analyses test to compute the effect sizes associated with all statistically significant pairwise in functionally defined regions, we established a hierarchy to designate poised (p), comparisons (wilcoxonRG function from the R package rcompanion version strong (s), and weak (w) promoter and enhancer states. Chromatin states E8, E9, 2.3.25)84. and E11 defined promoter states (P); E8 and E9 were strongly enriched at TSSs, To study the patterns of evolutionary conservation of the sequence underlying CGI, UMR (unmethylated regions), and open chromatin regions, while E11 was orthologous regulatory regions, we assigned each orthologous regulatory region a mostly located downstream the TSS; the presence of E14 defined poised promoter conservation score. We computed this score based on the phastCons30way states (pP); absence of E14 and presence of E9 or E11 defined strong promoter sequence conservation track29. To control for background sequence conservation states (sP); remaining P were classified as weak promoter states (wP). Non- levels, we first computed the average and standard deviations phastCons30way in promoter RE were assigned an enhancer state (E). The presence of E14 defined TADs defined in the cell line GM1287885 (Supplementary Fig. 13). Then, we used poised enhancer states (pE); absence of E14 and presence of E3, E4, E5, E6, and E12 these summary statistics to calculate the Z-score for each bp in every orthologous defined strong enhancer states (sE): E5 and E6 were strongly enriched LMRs (low regulatory region, using the average and standard deviations values of the TAD in methylated regions) whereas E3, E4, and E12 were highly abundant at introns; which each orthologous regulatory region was found. We averaged the Z-scores remaining E were classified as weak enhancer states (wE) (Supplementary Figs. 2 within each orthologous regulatory region in bins of 200 bp that overlap 50 bp with and 33). the next bin and assign each orthologous regulatory region the maximum Z-score One of the limitations of chromatin states is that bin assignments are based on values associated with its bins. We computed the Spearman rho correlation the presence or absence of particular epigenetic marks. However, oftentimes, the between the Z-scores and the number of species in which each orthologous lines separating different RE are blurry: e.g., the distinction between promoter and regulatory region was conserved, separately for each regulatory state. To determine enhancer states generally resides in the H3K4me3/H3K4me1 balance. Hence, some the statistical significance of these correlations we used randomization analysis. For 12 NATURE COMMUNICATIONS | (2021)12:3116 | https://doi.org/10.1038/s41467-021-23397-1 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-23397-1 ARTICLE each regulatory state, we created 1000 sets randomizing the Z-score associated with of statistically significant partial correlations with a cross-validation process to filter each orthologous regulatory region and calculated the Spearman correlation in out those relationships leading to overfitted reciprocal linear LASSO models each randomization. We determined the P value as the number of randomizations (significant partial correlations unlikely to be biologically meaningful). In our case, with a Spearman rho correlation value equal to or above the observed correlation in every analysis, we recovered those partial correlations recovered in at least four (Supplementary Fig. 15). of the five species without leading to overfitting when determining the reciprocal explanatory power in the remaining species. This protocol is intended to detect biologically relevant co-dependences out of the set of significant partial correlations Classification of RE in different types of components of gene regulatory and as a result, this approach filters out many significant partial correlations with architectures. This analysis was performed to classify RE based on their putative very low explanatory power. In fact, all the partial correlations recovered in any of role in the regulation of gene expression. the analyses performed showed very low P values (Benjamini–Hochberg’s We pre-classified each RE into a gene regulatory component based on their correction89, P < 1.8 × 10−22). In our case, given the relatively small amount of genomic location with respect to their corresponding species ENSEMBL release data, we focused on recovering those partial correlations that are likely to be 86 91 gene annotations. RE found up to 5 Kb upstream to the nearest TSS were relevant in any species. For these analyses, we used a modified version of the R code classified as gP. Additional RE located up to 10 Kb to the nearest TSS were provided by the authors (http://spcn.molgen.mpg.de/code/sparse_pcor.R/) to classified as prE. RE that overlapped a gene were classified as gE. Other RE that perform fivefold cross-validation analyses separating by species instead of the could not be linked to a gene based on their genomic proximity were initially original tenfold cross-validation protocol suitable for larger datasets. Network classified as distal enhancers (dE). visualizations were performed with Cytospace90. Then, we made use of available interaction data for the cell line GM12878 In a partial correlation model, direct co-dependencies are established between (HiC31, HiChIP-H3K27ac32, and ChIA-PET33) to map interactions between RE. individual variables. However, we know that coordination of the different histone Each interacting pair was mapped independently to hg38 coordinates using the marks within components is important to define the global epigenetic liftOver tool from the UCSCTOOLS/331 suite87, and only interactions for which configuration of a component (also captured in our epigenetic states), which itself both pairs could be mapped were kept. Subsequently, interactions were mapped to could be considered the relevant variable for this analysis. To better address this the nonhuman primate reference genome assemblies. For interspecies mappings, situation in our analysis, we defined a consensus signal for every component coordinates were mapped twice, going forward and backward, and only pairs that following the same approach established by WGCN35 to define eigengenes as could be mapped in both directions were kept. Interacting RE were defined as those representative variables of clusters of co-expressed genes. In brief, we defined that overlapped with each pair of any given interaction. First-order interactions eigencomponents as the variables summarizing the common signals of the different were annotated between promoters and enhancers, allowing the definition of PiE. histone marks in a component (actually calculated as the first PCA component of Second-order interactions were annotated between enhancer components (gE, prE, these five variables). So that eigencomponents keep the meaning of the activities, or PiE), allowing the definition of EiE (Fig. 4a and Supplementary Fig. 1). they were defined as codirectional with H3K27ac signals in each component Considering both their epigenetic state and regulatory component, RE were (eigenvectors negatively correlated with H3K27ac signals were multiplied by -1). separated into 30 (6 × 5) different subcategories. We used a Chi-square test to We performed a Sparse Partial Correlation Analysis of these ten eigencomponents identify the component-epigenetic state combinations enriched in orthologous and RNA-seq that recovers very clearly the structure of direct co-dependencies regulatory regions with fully conserved and species-specific epigenetic states between the epigenetic configuration of the different components and gene (Supplementary Fig. 26). expression (Fig. 4e and Supplementary Data 5). In addition, we defined the remaining unexplained signal of every histone mark by its eigencomponent as the residuals of a linear model of the original variables Gene expression levels and regulatory states in gene components. This ana- and the corresponding eigencomponent. A Sparse Partial Correlation Analysis of lysis was performed to explore whether gene expression levels are associated with these residuals (Supplementary Fig. 22 and Supplementary Data 6) shows that even the activity state of RE in each type of component. these residuals reflect the same inter-component structure and highlights that our To investigate the influence of the activity state of RE in each type of eigencomponents miss some relevant information for the definition of this component on gene expression levels, we classified one-to-one orthologous regulatory coordination (mainly weaker co-dependencies involving promoter states protein-coding genes, separately for each species, into six mutually excluding in intragenic and PiE and enhancer states in promoters). categories, one for each regulatory state within each type of component To assess to what extent eigencomponents reflect the behavior of the whole (component-state combinations). Whereas genes can only be associated with one network of co-dependencies of the histone marks or of each of the specific histone genic promoter and hence, they can only be classified into one category for gP marks, we also performed SPCAs using the actual ChIP-seq enrichment signals. A depending on the corresponding epigenetic state of the RE, genes can be associated global partial correlation analysis considering all 51 variables shows a very clear with more than enhancer component (gE, prE, PiE, and EiE). In those cases, we structure of direct co-dependencies with a strong intra-component contribution for classified genes into a given component-state category accordingly to the presence the two states of every single component and a clear but more modest exclusive of at least one RE with a given epigenetic state in that component using the inter-component contribution (Supplementary Data 7). Analyses to determine the following state hierarchy: pE > pP > sE > sP > wE > wP (Supplementary Fig. 20). To Sparse Partial Correlation Network of each of the histone marks and RNA-seq statistically assess the influence of each state in each component, we used (1) without considering the possible influence of the remaining histone marks Kruskal–Wallis test (kruskal.test function as implemented in R)82 to test whether (Supplementary Data 8) retrieve very similar networks pointing to the common the distributions of the expression levels of genes associated with each component- backbone of inter-component co-dependences reflected in our SPCA of the state combination were different for the different regulatory states, (2) eigencomponents. Dwass–Steel–Critchlow–Fligner test to assess the significance of every pairwise Our dataset of regulatory components shows a quite unbalanced contribution of comparison (dscfAllPairsTest function from the R package PCMRplus version the components to the architectures, with gE being the most abundant type of 1.4.4)83, and (3) Glass rank biserial correlation coefficient effect size for component and promoter-interacting and EiE being the least abundant Mann–Whitney U-test to compute the effect sizes associated with all statistically (Supplementary Fig. 18). These differences could be at least partially related to our significant pairwise comparisons (wilcoxonRG function from the R package inability to recover some of the chromatin interaction-mediated regulatory rcompanion version 2.3.25)84 (Supplementary Fig. 21). associations. More importantly, this imbalance, if not real, could affect the ability of our partial correlation networks to reflect the contribution of those components Partial correlation analysis. This analysis was performed to understand in greater less represented in our datasets. To explore this point, we recovered the subset of detail how the epigenetic signals associated with RE with different epigenetic states genes (an average of 1068 genes per sample) with full architectures (those with at and in different components relate to each other and with gene expression. least one element in every type of component) and repeated all the Sparse Partial To disentangle the network of direct co-dependencies between the different Correlation Analyses explained above with this dataset of genes. In all the cases, we components, regulatory states, histone marks, and gene expression, we performed a obtained very similar results, recovering fewer relevant partial correlations due to series of partial correlation analyses34,88. To tackle the diversity of architectures the smaller number of genes, but with no signal of any relevant difference in the detected for the different genes, we added up the calibrated signal of all the RE with global structure of the coordinated network of components and gene expression a given regulatory state (promoter or enhancer) in a given type of component for (Supplementary Fig. 23 and Supplementary Data 5–8). any gene architecture. This decision was based on the observed relationship All the components of the connected network can be very influential in gene between the number of strong elements in a gene architecture and the expression expression through their direct or indirect connection with gene expression. level of its target gene. Separation of histone signals in each type of component However, our Sparse Partial Correlation Networks point consistently to the direct between those contributing to a promoter or to an enhancer was intended to reflect co-dependency of RNA-seq with the genic promoter and gE components and the the potential differences in their role in gene expression regulation. As a result of co-dependency between them. To quantify the explanatory power of these this design, our system has 51 variables (RNA-seq signal + 5 histone mark signals dependencies for gene expression, we performed a simple generalized linear model x 2 regulatory states x 5 components) and 57,370 cases (5737 genes x 5 species x (glm function as implemented in R82) for RNA-seq using H3K27ac, H3K27me3, 2 samples). and H3K36me3 signals in gP and gE and the interactions between them. This All partial correlation analyses were performed using an adaptation of a recently model was able to explain 67% of the gene expression variance (Supplementary published Sparse Partial Correlation Analysis protocol34 based on the continuous Data 8), a percentage 5% higher than the 62% explained by a naïve model values of the accumulated ChIP-seq signals (instead of their ranks) to take including the signals of all histone marks in all the components but no interaction advantage of their pseudo-quantitative nature. This protocol combines the recovery between them (Supplementary Table 8), supporting that gP and gE contained NATURE COMMUNICATIONS | (2021)12:3116 | https://doi.org/10.1038/s41467-021-23397-1 | www.nature.com/naturecommunications 13 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-23397-1 nearly all the epigenetic information needed to define gene expression levels in TPM and then we excluded from the analyses those genes that had an accumulated our data. expression value in all tissues below 0.1xNumber of tissues (n = 29 tissues). For each component-state combination, differences in median expression across tissues were assessed with the Friedman test using the friedman.test function as imple- Differential gene expression analyses. This analysis was performed to find genes mented in R82. We used the Wilcoxon–Nemenyi–McDonald–Thompson test with expression changes between species and investigate how these expression implemented in the pWNMT function of the R package NSM3 (version 1.14)96 to changes are connected to changes in gene regulatory components with different assess whether expression levels were significantly different for all pairwise tissue epigenetic states. combinations. Then, we made use of the rank-biserial correlation to calculate the We identified genes with differential expression levels across species using the effect sizes for all statistically significant pairwise tests with the wilcoxonPairedRC iDEGES/edgeR pipeline in the R package TCC (version 1.12.1)91,92 at an FDR of function of the R package rcompanion (version 2.3.25)84. 0.1 and testing all species pairwise comparisons. Then, we determined the patterns We then evaluated the tissue-specificity of the genes associated with the of differential expression, species, and direction of the gene expression change, different component-state combinations. For this, we calculated the tissue using a two-step approach. For every gene, the Q-values obtained in species specificity index97 (τ, tau) for each gene, which is defined as: pairwise comparisons were ordered from smallest to largest. Different expression labels were then assigned to each species according to the ordered Q-values. Once N all species had an assigned label, the average normalized expression values between τ ¼ ∑ ð1  xi Þ=N  1 ð1Þ i¼1 groups were compared to determine the directionality of the change. We separate differentially expressed genes into two categories: genes with species-specific expression changes and genes with non-species-specific expression changes. where N is the number of tissues and xi is the expression value normalized by the To investigate the relationship between changes in gene expression and changes maximum expression value. This value ranges from 0, for housekeeping genes, to 1, in the regulatory architecture of a gene, for every type of regulatory component, we for tissue-specific genes (values above 0.8 are used to identify tissue-specific run a Wilcoxon signed-rank test evaluating whether the number of RE with a given genes)98. Tissue-specificity indices were calculated for all genes included in the regulatory state in that particular regulatory component was significantly latest GTEx release36. Gene expression levels (median TMP) of grouped tissue associated with higher expression levels, for strong and weak activities, or lower categories (Supplementary Table 9) were normalized within and across tissues expression levels, for poised activities. P values obtained for each regulatory role before calculating τ as implemented in the R package tispec (version 0.99.0)99. The were corrected for multiple testing using the Benjamini–Hochberg procedure89. calcTau function from this package provides a tau value for each gene and also a tau expression fraction for each tissue (which also ranges from 0 to 1) that indicates the specificity of a given gene for that tissue. Over-representation analyses (ORA) of functional annotations. This analysis After calculating τ values, we compared their distributions between gene was performed to explore the functional role of both conserved as well as species- datasets with the Kruskall–Wallis test and assessed the significance of every specific regulatory states in certain regulatory/gene components (overrepresented pairwise comparison with the Dwass–Steel–Critchlow–Fligner test combinations of regulatory states in specific components, Supplementary Fig. 26). (dscfAllPairsTest function from the R package PMCMRplus version 1.4.4)83. Glass For that, we defined sets of genes associated with fully conserved and species- rank biserial correlation coefficient was used to compute the effect sizes associated specific component-epigenetic state combinations and explored their functional with all statistically significant pairwise comparisons using the wilcoxonRG enrichments. To ensure the representativeness of the functional enrichments, for function from the R package rcompanion version 2.3.2584 (P < 0.05). the gene lists associated with each type of component, we excluded genes associated with components with different epigenetic states activities (i.e., genes associated Association of genes containing gE with signals of positive selection and with both conserved strong and weak gE) or associated with both conserved and accelerated evolution. The main objective of this analysis was to examine whether species-specific components levels (i.e., genes associated with both a conserved and hswEge-containing genes had been targeted by different evolutionary forces. We a species-specific weak gE) and kept those gene lists with a minimum of 15 genes built a database of human genomic regions with previously detected signals of for enrichment analyses. Of note, orangutan-specific component-epigenetic state positive selection in humans50–54 and human acceleration55–57. In brief, for the combinations were excluded from the analysis because they were defined using signals of positive selection, we retrieved Supplemental Files 1 and 2 from only one LCL replicate (see above) and they are likely to be enriched in inter- Peyrégne50 et al. Table S5 from Racimo51 et al. Table S3 (regions >1600 genera- individual variable activities. tions) from Zhou52 et al. Table S19a.3 from Prüfer53 et al. and Table S37 from Over-representation of Gene Ontology (GO) terms was performed using the Green54 et al. For the datasets with different signals of human acceleration we took WebGestaltR function from the R package WebGestaltR (version 0.4.3; http://www. haDHSs (Supplemental Table 2) from Gittelman55 et al. HARs (Supplementary webgestalt.org) with minNum = 25 and remaining default options. This function Info) from Lindblad-Toh56 et al., and CNSs (Table S1) from Prabhakar57 et al. To controls the FDR by applying Benjamini–Hochberg procedure (default threshold account for the potential effect that gBGC could have in these regions, we filtered FDR = 0.05)89,93. Previous analyses have shown that recent enhancers tend to out all those regions overlapping with gBGC tracts in the human genome100. After occur in the same genes that already have highly conserved enhancers10. To control this filtering, BEDtools101 was used to assign these regions to both protein-coding for the particular background of each component, we built different background and non-coding genes following similar criteria to those used for building the gene gene sets, including the set of human genes associated with at least one-to-one regulatory architectures (Methods’ section Classification of RE in different types of orthologous regulatory regions of each type of component. This way, we have components of gene regulatory architectures). We assigned these regions to a specific and different backgrounds for gP, gE, and PiE. To represent and compare protein-coding gene if they were located within the gene or up to 5 Kb upstream of enriched GO terms between component-state combinations, we performed a its TSS. Then, we made use of available interaction data for the cell line GM12878 clustering analysis of all significantly enriched GO terms using REVIGO94. We (HiC31, HiChIP-H3K27ac32, and ChIA-PET33) to assign this composite of posi- associated each GO term with the proportion of genes from each component-state tively selected regions and acceleration regions to their interacting genes. We ended combination that overlapped that GO term. In the case of GO terms enriched in up with a set of 4747 genes associated with at least one positively selected region or more than one gene set, we chose the highest proportion of genes. We used this list one accelerated region in the human lineage. We computed the overlaps between as input for REVIGO. Given that REVIGO only reports the clustering of ~350 GO this gene list and the lists of genes associated with the different component-state terms and our input list was larger than that, we used the R package GofuncR combinations at different degrees and backgrounds. We used one-tailed Fisher’s (version 1.8.0)95 to retrieve the parent GO terms of the remaining unassigned GO exact test to assess the enrichment significance. terms and add them to the corresponding group as defined by REVIGO. REVIGO group names were manually assigned, taking into account the most representative parent term (Supplementary Data 10). Analyses of the density of human-fixed single nucleotide changes (hSNCs) in gE with weak enhancer states. The purpose of this analysis was to assess the co- occurrence of human-nucleotide changes in human-specific RE. In order to study Analysis of tissue-specific gene expression. The aim of this analysis was to the distribution of these human-fixed changes, we first generated a dataset with inspect the expression pattern and tissue-specificity of certain RE in particular human-specific changes. We used sequencing data from a diversity panel of 27 components (both conserved as well as species-specific) across different human cell orangutans, 42 gorillas, 11 bonobos, and 61 chimpanzees102–104, as well as 19 types or tissues. We defined sets of human genes associated with fully conserved modern humans from the 1000 genomes project105, all mapped to the human component-state combinations and human genes associated with human-specific reference assembly hg19. We applied a basic filter for each site in each individual gains/losses of RE. Note that these gene lists are not mutually exclusive since a gene (sequencing coverage >3 and <100) and kept sites where at least half of the indi- can be associated with different types of conserved or species-specific component- viduals in a given species had sufficient data. Furthermore, at least 90% of the kept state combinations (e.g., a gene with both a human-specific intragenic enhancer individuals at a given site in a given species had to share the same allele; otherwise, with weak activity and a fully conserved gE with a strong activity). We obtained the site was labeled as polymorphic in the population. Indels and triallelic sites expression levels (median TPM values) across a collection of different tissues from were removed, and only biallelic sites were kept. We used data from a macaque the latest GTEx release (v8)36. We only included tissues with at least 70 samples diversity panel106, applying the same filters described above. The allele at mono- and grouped tissue subregions into the same tissue category, as stated in Supple- morphic sites was added using bedtools getfasta101 from the macaque reference mentary Table 9. For each component-state combination, we followed a two-step genome rheMac8. Since this panel uses the macaque reference genome, we per- approach to remove consistently low-expressed genes across tissues. For that, we formed a liftover to hg19 using the R package rtracklayer107 and merged the data first assigned a value of 0 to all genes with a median expression level below 0.1 with the great ape diversity panel. 14 NATURE COMMUNICATIONS | (2021)12:3116 | https://doi.org/10.1038/s41467-021-23397-1 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-23397-1 ARTICLE Lineage-specific changes were retrieved as polymorphisms with sufficient 11. Vermunt, M. W. et al. Epigenomic annotation of gene regulatory alterations information. Hence, human-specific changes (hSNCs) were defined as positions during evolution of the primate brain. Nat. Neurosci. 19, 494–503 (2016). where each species carry only or mostly one allele within their respective 12. Reilly, S. K. et al. Evolutionary genomics. Evolutionary changes in promoter population, the majority of individuals in each population have a genotype call at and enhancer activity during human corticogenesis. Science 347, 1155–1159 sufficient coverage, and the human allele differs from the allele in the other (2015). populations. 13. Villar, D. et al. Enhancer evolution across 20 mammalian species. Cell 160, BEDtools101 was used to annotate those hSNCs in conserved or hswEgE, and 554–566 (2015). the density of changes was calculated as the number of hSNCs present in each 14. Prescott, S. L. et al. Enhancer divergence and cis-regulatory evolution in the enhancer divided by the length of the enhancer. human and chimp neural crest. Cell 163, 68–83 (2015). To determine which hswEgE were enriched in human-specific changes, we 15. Castelijns, B. et al. Recently evolved enhancers emerge with high compared their density to what would be expected at random. For that, we first interindividual variability and less frequently associate with disease. Cell Rep. established the number of hSNCs that fall in human gE with weak enhancer states 31, 107799 (2020). associated with one-to-one orthologous regulatory regions (our universe of 16. Danko, C. G. et al. Dynamic evolution of regulatory element ensembles in enhancers). In each simulation, this number of mutations was randomly placed in this primate CD4 T cells. Nat. Ecol. Evol. 2, 537–548 (2018). universe, and we computed the density for each of the hswEgE (10,000 simulations). 17. Trizzino, M. et al. Transposable elements are the primary source of novelty in With this approach, we corrected for the differences in the length of the enhancers. primate gene regulation. Genome Res. 27, 1623–1633 (2017). The P value for each enhancer was computed as the number of simulations with a 18. Zhou, X. et al. Epigenetic modifications are associated with inter-species gene density equal to or above the observed density for that particular enhancer. All P expression variation in primates. Genome Biol. 15, 547 (2014). values were corrected by multiple testing using the Bonferroni method with the 19. Blake, L. E. et al. A comparison of gene expression and DNA methylation number of tests equal to the number of hswEgE. patterns across tissues and species. Genome Res. 30, 250–262 (2020). We then assessed whether the number of enhancers that were statistically enriched in hSNCs (or number of hits) was greater than what would be expected at 20. Gokhman, D. et al. Differential DNA methylation of vocal and facial anatomy random. In order to do that, for each enhancer, we defined its mutation density genes in modern humans. Nat. Commun. 11, 1189 (2020). critical value adjusting by multiple testing and using the simulated values. For 21. Hernando-Herraez, I. et al. The interplay between DNA methylation and example, in a hypothetical case of 100 enhancers and 10,000 simulations, for each sequence divergence in recent human evolution. Nucleic Acids Res. 43, enhancer, we would order its simulated density of hSNCs from smallest to largest 8204–8214 (2015). and take the fifth value as the critical one (given that our chosen alpha equals 5%, 22. Hernando-Herraez, I. et al. Dynamics of DNA methylation in recent human but it has to be corrected by 100 tests; therefore it becomes 0.05%). Once we and great ape evolution. PLoS Genet. 9, e1003763 (2013). established a critical value for each human-specific intragenic weak enhancer, we 23. Zeng, J. et al. Divergent whole-genome methylation maps of human and determined, for each simulation, how many enhancers had a density equal to or chimpanzee brains reveal epigenetic basis of human regulatory evolution. Am. above their corresponding critical value. Finally, we computed the P value J. Hum. Genet. 91, 455–465 (2012). comparing the number of artificial hits in each simulation with the number of 24. Eres, I. E., Luo, K., Hsiao, C. J., Blake, L. E. & Gilad, Y. Reorganization of 3D observed hits. genome structure may contribute to gene regulatory evolution in primates. PLoS Genet. 15, e1008278 (2019). 25. Lowdon, R. F., Jang, H. S. & Wang, T. Evolution of epigenetic regulation in Data availability vertebrate genomes. Trends Genet. 32, 269–283 (2016). The raw fastq files from the genomic, transcriptomic, and epigenomic data generated and 26. Ernst, J. & Kellis, M. Chromatin-state discovery and genome annotation with used for the analyses in this study were uploaded to the Sequence Read Archive (SRA) ChromHMM. Nat. Protoc. 12, 2478–2492 (2017). with the BioProject accession number PRJNA563344. Data from human cell line 27. Carelli, F. N., Liechti, A., Halbert, J., Warnefors, M. & Kaessmann, H. GM12878 used in the study was obtained from SRA: RNA-seq accessions SRR998197 Repurposing of promoters and enhancers during mammalian evolution. Nat. and SRR998198, ChIP-seq for H3K4me3 accession SRR998194, H3K4me1 accession Commun. 9, 4066 (2018). SRR998191, H3K27ac accession SRR998178, H3K36me3 accession SRR998187, and 28. Jacques, P.-É., Jeyakani, J. & Bourque, G. The majority of primate-specific H3K27me3 accession SRR998189; GM12878 ChIP-seq input accession and ATAC-seq regulatory sequences are derived from transposable elements. PLoS Genet. 9, data accession. Interaction data available for GM12878: HiC accession GSE63525, e1003504 (2013). HiChIP-H3K27ac, and ChIA-PET accession GSM1872887. Expression data is available 29. Siepel, A. et al. Evolutionarily conserved elements in vertebrate, insect, worm, from GTEx v8. The remaining data are available within the manuscript, Supplementary and yeast genomes. Genome Res. 15, 1034–1050 (2005). Information, Supplementary Data or from the corresponding authors upon reasonable 30. Ong, C.-T. & Corces, V. G. Enhancer function: new insights into the request. regulation of tissue-specific gene expression. Nat. Rev. Genet. 12, 283–293 (2011). Received: 23 September 2020; Accepted: 21 April 2021; 31. Rao, S. S. P. et al. A 3D map of the human genome at kilobase resolution reveals principles of chromatin looping. Cell 159, 1665–1680 (2014). 32. Mumbach, M. R. et al. Enhancer connectome in primary human cells identifies target genes of disease-associated DNA elements. Nat. Genet. 49, 1602–1612 (2017). 33. Tang, Z. et al. CTCF-mediated human 3D genome architecture reveals chromatin topology for transcription. Cell 163, 1611–1627 (2015). References 34. Lasserre, J., Chung, H.-R. & Vingron, M. Finding associations among histone 1. Britten, R. J. & Davidson, E. H. Gene regulation for higher cells: a theory. modifications using sparse partial correlation networks. PLoS Comput. Biol. 9, Science 165, 349–357 (1969). e1003168 (2013). 2. Britten, R. J. & Davidson, E. H. Repetitive and non-repetitive DNA sequences 35. Langfelder, P. & Horvath, S. WGCNA: an R package for weighted correlation and a speculation on the origins of evolutionary novelty. Q. Rev. Biol. 46, network analysis. BMC Bioinforma. 9, 1–13 (2008). 111–138 (1971). 36. Consortium, T. G., The GTEx Consortium. The GTEx Consortium atlas of 3. Cardoso-Moreira, M. et al. Gene expression across mammalian organ genetic regulatory effects across human tissues. Science 369, 1318–1330 development. Nature 571, 505–509 (2019). (2020). 4. Xu, C. et al. Human-specific features of spatial gene expression and regulation 37. Bogdanovic, O. et al. Dynamics of enhancer chromatin signatures mark the in eight brain regions. Genome Res. 28, 1097–1110 (2018). transition from pluripotency to cell specification during embryogenesis. 5. Zhu, Y. et al. Spatiotemporal transcriptomic divergence across human and Genome Res. 22, 2043–2053 (2012). macaque brain development. Science 362, eaat8077 (2018). 38. Bonn, S. et al. Tissue-specific analysis of chromatin state identifies temporal 6. Blake, L. E. et al. A comparative study of endoderm differentiation in humans signatures of enhancer activity during embryonic development. Nat. Genet. and chimpanzees. Genome Biol. 19, 162 (2018). 44, 148–156 (2012). 7. Khan, Z. et al. Primate transcript and protein expression levels evolve under 39. Rada-Iglesias, A. et al. A unique chromatin signature uncovers early compensatory selection pressures. Science 342, 1100–1104 (2013). developmental enhancers in humans. Nature 470, 279–283 (2011). 8. Merkin, J., Russell, C., Chen, P. & Burge, C. B. Evolutionary dynamics of gene 40. Creyghton, M. P. et al. Histone H3K27ac separates active from poised and isoform regulation in Mammalian tissues. Science 338, 1593–1599 (2012). enhancers and predicts developmental state. Proc. Natl Acad. Sci. USA 107, 9. Brawand, D. et al. The evolution of gene expression levels in mammalian 21931–21936 (2010). organs. Nature 478, 343–348 (2011). 41. Ernst, J. et al. Mapping and analysis of chromatin state dynamics in nine 10. Berthelot, C., Villar, D., Horvath, J. E., Odom, D. T. & Flicek, P. Complexity human cell types. Nature 473, 43–49 (2011). and conservation of regulatory landscapes underlie evolutionary resilience of 42. Enard, W. et al. Molecular evolution of FOXP2, a gene involved in speech and mammalian gene expression. Nat. Ecol. Evol. 2, 152–163 (2018). language. Nature 418, 869–872 (2002). NATURE COMMUNICATIONS | (2021)12:3116 | https://doi.org/10.1038/s41467-021-23397-1 | www.nature.com/naturecommunications 15 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-23397-1 43. Mozzi, A. et al. The evolutionary history of genes involved in spoken and 75. Ducos, J. et al. Polyreactivity of human monoclonal antibodies: human anti- written language: beyond FOXP2. Sci. Rep. 6, 22157 (2016). Rh monoclonal antibodies of IgM isotype are frequently polyreactive. Vox 44. Kalebic, N. et al. Neocortical expansion due to increased proliferation of basal Sanguinis 61, 196–204 (1991). progenitors is linked to changes in their morphology. Cell Stem Cell 24, 76. Doxiadis, G. G. M. et al. Haplotype diversity generated by ancient 535–550.e9 (2019). recombination-like events in the MHC of Indian rhesus macaques. 45. Kuhlwilm, M. & Boeckx, C. A catalog of single nucleotide changes Immunogenetics 65, 569–584 (2013). distinguishing modern humans from archaic hominins. Sci. Rep. 9, 8463 77. Groot, N. Gde et al. A specialist macaque MHC class I molecule with HLA- (2019). B*27–like peptide-binding characteristics. J. Immunol. 199, 3679–3690 (2017). 46. Finn, S. & Civetta, A. Sexual selection and the molecular evolution of ADAM 78. Morey, L. et al. Nonoverlapping functions of the polycomb group Cbx family proteins. J. Mol. Evol. 71, 231–240 (2010). of proteins in embryonic stem cells. Cell Stem Cell 10, 47–62 (2012). 47. Riordan, J. R. Identification of the cystic fibrosis gene: cloning and 79. Buenrostro, J. D., Giresi, P. G., Zaba, L. C., Chang, H. Y. & Greenleaf, W. J. characterization of complementary DNA. Trends Genet. 5, 363 (1989). Transposition of native chromatin for fast and sensitive epigenomic profiling 48. Poolman, E. M. & Galvani, A. P. Evaluating candidate agents of selective of open chromatin, DNA-binding proteins and nucleosome position. Nat. pressure for cystic fibrosis. J. R. Soc. Interface 4, 91–98 (2007). Methods 10, 1213–1218 (2013). 49. Racimo, F. et al. Archaic adaptive introgression in TBX15/WARS2. Mol. Biol. 80. Tharwat, A., Gaber, T., Ibrahim, A. & Hassanien, A. E. Linear discriminant Evol. https://doi.org/10.1093/molbev/msw283 (2016). analysis: a detailed tutorial. AI Commun. 30, 169–190 (2017). 50. Peyrégne, S., Boyle, M. J., Dannemann, M. & Prüfer, K. Detecting ancient 81. Venables, W. N. & Ripley, B. D. Modern Applied Statistics with S (Springer, positive selection in humans using extended lineage sorting. Genome Res. 27, 2002). 1563–1572 (2017). 82. R Core Team. R: a language and environment for statistical computing (The R 51. Racimo, F. Testing for ancient selection using cross-population allele Foundation, 2015). frequency differentiation. Genetics 202, 733–750 (2016). 83. Pohlert, T. PMCMRplus: calculate pairwise multiple comparisons of mean 52. Zhou, H. et al. A chronological atlas of natural selection in the human genome rank sums extended (The R Foundation, 2020). during the past half-million years. Cold Spring Harb. Lab. https://doi.org/ 84. Mangiafico, S. rcompanion: functions to support extension education program 10.1101/018929 (2015). evaluation (The R Foundation, 2020). 53. Prüfer, K. et al. The complete genome sequence of a Neanderthal from the 85. Wang, Y. et al. The 3D Genome Browser: a web-based browser for visualizing Altai Mountains. Nature 505, 43–49 (2013). 3D genome organization and long-range chromatin interactions. Genome Biol. 54. Green, R. E. et al. A draft sequence of the Neandertal genome. Science 328, 19, 151 (2018). 710–722 (2010). 86. Aken, B. L. et al. Ensembl 2017. Nucleic Acids Res. 45, D635–D642 (2017). 55. Gittelman, R. M. et al. Comprehensive identification and analysis of human 87. Kuhn, R. M., Haussler, D. & Kent, W. J. The UCSC genome browser and accelerated regulatory DNA. Genome Res. 25, 1245–1255 (2015). associated tools. Brief. Bioinforma. 14, 144–161 (2013). 56. Lindblad-Toh, K. et al. A high-resolution map of human evolutionary 88. Juan, D. et al. Epigenomic co-localization and co-evolution reveal a key role constraint using 29 mammals. Nature 478, 476–482 (2011). for 5hmC as a communication hub in the chromatin network of ESCs. Cell 57. Prabhakar, S., Noonan, J. P., Pääbo, S. & Rubin, E. M. Accelerated evolution of Rep. 14, 1246–1257 (2016). conserved noncoding sequences in humans. Science 314, 786 (2006). 89. Benjamini, Y., Drai, D., Elmer, G., Kafkafi, N. & Golani, I. Controlling the false 58. Demontis, D. et al. Discovery of the first genome-wide significant risk loci for discovery rate in behavior genetics research. Behav. Brain Res. 125, 279–284 attention deficit/hyperactivity disorder. Nat. Genet. 51, 63–75 (2019). (2001). 59. Meur, N. L. et al. MEF2C haploinsufficiency caused by either microdeletion of 90. Shannon, P. et al. Cytoscape: a software environment for integrated models of the 5q14.3 region or mutation is responsible for severe mental retardation biomolecular interaction networks. Genome Res. 13, 2498–2504 (2003). with stereotypic movements, epilepsy and/or cerebral malformations. J. Med. 91. Bolstad, B. M., Irizarry, R. A., Astrand, M. & Speed, T. P. A comparison of Genet. 47, 22–29 (2010). normalization methods for high density oligonucleotide array data based on 60. Hyde, C. L. et al. Identification of 15 genetic loci associated with risk of major variance and bias. Bioinformatics 19, 185–193 (2003). depression in individuals of European descent. Nat. Genet. 48, 1031–1036 92. Sun, J., Nishiyama, T., Shimizu, K. & Kadota, K. TCC: an R package for (2016). comparing tag count data with robust normalization strategies. BMC 61. Katoh, Y. et al. The clavesin family, neuron-specific lipid- and clathrin- Bioinforma. 14, 219 (2013). binding Sec14 proteins regulating lysosomal morphology. J. Biol. Chem. 284, 93. Liao, Y., Wang, J., Jaehnig, E. J., Shi, Z. & Zhang, B. WebGestalt 2019: gene set 27646–27654 (2009). analysis toolkit with revamped UIs and APIs. Nucleic Acids Res. 47, 62. Long, H. et al. Conserved roles for slit and robo proteins in midline W199–W205 (2019). commissural axon guidance. Neuron 42, 213–223 (2004). 94. Supek, F., Bošnjak, M., Škunca, N. & Šmuc, T. REVIGO summarizes and 63. Andrews, W. et al. Robo1 regulates the development of major axon tracts and visualizes long lists of gene ontology terms. PLoS ONE 6, e21800 (2011). interneuron migration in the forebrain. Development 133, 2243–2252 (2006). 95. Grote, S. GOfuncR: gene ontology enrichment using FUNC (The R 64. Yin, X.-M. et al. Familial paroxysmal kinesigenic dyskinesia is associated with Foundation, 2020). mutations in the KCNA1 gene. Hum. Mol. Genet. 27, 757–758 (2018). 96. Grant Schneider, Chicken, E. & Becvarik, R. NSM3: functions and datasets to 65. Carter, K. L., Cahir-McFarland, E. & Kieff, E. Epstein-Barr virus-induced accompany (The R Foundation, 2020). changes in B-lymphocyte gene expression. J. Virol. 76, 10427–10436 (2002). 97. Yanai, I. et al. Genome-wide midrange transcription profiles reveal expression 66. Hansen, K. D. et al. Increased methylation variation in epigenetic domains level relationships in human tissue specification. Bioinformatics 21, 650–659 across cancer types. Nat. Genet. 43, 768–775 (2011). (2005). 67. Sugawara, H. et al. Comprehensive DNA methylation analysis of human 98. Kryuchkova-Mostacci, N. & Robinson-Rechavi, M. A benchmark of gene peripheral blood leukocytes and lymphoblastoid cell lines. Epigenetics 6, expression tissue-specificity metrics. Brief. Bioinform. 18, 205–214 (2017). 508–515 (2011). 99. Condon, K. tispec: calculates tissue specificity from RNA-seq data. (2020). 68. Hussain, T. & Mulherkar, R. Lymphoblastoid cell lines: a continuous in vitro 100. Capra, J. A., Hubisz, M. J., Kostka, D., Pollard, K. S. & Siepel, A. A model- source of cells to study carcinogen sensitivity and DNA repair. Int. J. Mol. Cell based analysis of GC-biased gene conversion in the human and chimpanzee Med. 1, 75–87 (2012). genomes. PLoS Genet. 9, e1003684 (2013). 69. Khaitovich, P., Enard, W., Lachmann, M. & Pääbo, S. Evolution of primate 101. Quinlan, A. R. & Hall, I. M. BEDTools: a flexible suite of utilities for gene expression. Nat. Rev. Genet. 7, 693–702 (2006). comparing genomic features. Bioinformatics 26, 841–842 (2010). 70. Pai, A. A., Bell, J. T., Marioni, J. C., Pritchard, J. K. & Gilad, Y. A genome-wide 102. Prado-Martinez, J. et al. Great ape genetic diversity and population history. study of DNA methylation patterns and gene expression levels in multiple Nature 499, 471–475 (2013). human and chimpanzee tissues. PLoS Genet. 7, e1001316 (2011). 103. de Manuel, M. et al. Chimpanzee genomic diversity reveals ancient admixture 71. Shibata, Y. et al. Extensive evolutionary changes in regulatory element activity with bonobos. Science 354, 477–481 (2016). during human origins are associated with altered gene expression and positive 104. Nater, A. et al. Morphometric, behavioral, and genomic evidence for a new selection. PLoS Genet. 8, e1002789 (2012). Orangutan species. Curr. Biol. 27, 3576–3577 (2017). 72. Borsari, B. et al. Intronic enhancers regulate the expression of genes involved 105. Consortium, T. 1000 G. P. & the 1000 Genomes Project Consortium. A global in tissue-specific functions and homeostasis. Cold Spring Harb. Lab. https:// reference for human genetic variation. Nature 526, 68–74 (2015). doi.org/10.1101/2020.08.21.260836 (2020). 106. Xue, C. et al. The population genomics of rhesus macaques (Macaca mulatta) 73. Lahti, D. C. et al. Relaxed selection in the wild. Trends Ecol. E 24, 487–496 based on whole-genome sequences. Genome Res. 26, 1651–1662 (2016). (2009). 107. Lawrence, M., Gentleman, R. & Carey, V. rtracklayer: an R package for 74. Uebbing, S. et al. Massively parallel discovery of human-specific substitutions interfacing with genome browsers. Bioinformatics 25, 1841–1842 (2009). that alter enhancer activity. Proc. Natl. Acad. Sci. USA 118, e2007049118 108. Kasowski, M. et al. Extensive variation in chromatin states across humans. (2021). Science 342, 750–752 (2013). 16 NATURE COMMUNICATIONS | (2021)12:3116 | https://doi.org/10.1038/s41467-021-23397-1 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-23397-1 ARTICLE Acknowledgements Competing interests R.G.-P. was supported by a fellowship from MICINN (FPU13/01823). P.E.-C. was The authors declare no competing interests. supported by a Formació de Personal Investigador fellowship from Generalitat de Cat- alunya (FI_B00122). M.K. was supported by a Deutsche Forschungsgemeinschaft (DFG) fellowship (KU 3467/1-1) and the Postdoctoral Junior Leader Fellowship Program from Additional information Supplementary information The online version contains supplementary material “la Caixa” Banking Foundation (LCF/BQ/PR19/11700002). D.J. was supported by a Juan available at https://doi.org/10.1038/s41467-021-23397-1. de la Cierva fellowship (FJCI2016-29558) from MICINN. T.M-B. is supported by funding from the European Research Council (ERC) under the European Union’s Correspondence and requests for materials should be addressed to R.G.-P., D.J. or Horizon 2020 research and innovation program (grant agreement EC-H2020-ERC-CoG- T.M.-B. ApeGenomeDiversity-864203), BFU2017-86471-P (AEI/FEDER, UE), “Unidad de Excelencia María de Maeztu”, funded by the AEI (CEX2018-000792-M), Howard Hughes Peer review information Nature Communications thanks Chris Ponting and the other, International Early Career, NIH 1R01HG010898-01A1, Obra Social “La Caixa” and anonymous, reviewers for their contribution to the peer review of this work. Peer Secretaria d’Universitats i Recerca and CERCA Program del Departament d’Economia i reviewer reports are available. Coneixement de la Generalitat de Catalunya (GRC 2017 SGR 880). G.M., V.D.C., and L. D.C. were supported by grants from the Spanish of Economy, Industry, and Competi- Reprints and permission information is available at http://www.nature.com/reprints tiveness (MEIC) (BFU2016-75008-P) and G.M. was also supported by the “Convocatoria de Ayudas Fundación BBVA a Investigadores, Innovadores y Creadores Culturales”. J.L. Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in G.-S. was supported by the Spanish government (grants BFU2016-74961-P), an insti- published maps and institutional affiliations. tutional grant Unidad de Excelencia María de Maeztu (MDM-2016-0687) and the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 740041). A.N. was supported by Fondo Europeo de Desarrollo Regional (FEDER) with project grants BFU2016-77961-P and Open Access This article is licensed under a Creative Commons PGC2018- 101927-B-I00 and by the Spanish National Institute of Bioinformatics (PT17/ Attribution 4.0 International License, which permits use, sharing, 0009/0020). Figures 1a, 2a, 4a, 6b, and Supplementary Fig. 1 were fully or partially adaptation, distribution and reproduction in any medium or format, as long as you give created using BioRender.com. We dedicate this manuscript to the memory of Jose Luis appropriate credit to the original author(s) and the source, provide a link to the Creative Gomez Skarmeta: a leader, an inspiration, and a friend. 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 indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory Author contributions regulation or exceeds the permitted use, you will need to obtain permission directly from T.M.-B. and J.L.G.-S. conceived the study; D.J. designed and supervised the analyses; L.D. the copyright holder. To view a copy of this license, visit http://creativecommons.org/ C. supervised the work of G.M. and V.D.C.; A.B. procured nonhuman great ape cell lines; licenses/by/4.0/. A.N. provided helpful insights; R.G.-P., G.M., and V.D.C. performed the experimental work; R.G.-P., P.E.-C., D.J., I.L., M.R., and M.K analyzed the data; D.J., R.G.-P., and P.E.- C. wrote the manuscript with input from all co-authors. © The Author(s) 2021 NATURE COMMUNICATIONS | (2021)12:3116 | https://doi.org/10.1038/s41467-021-23397-1 | www.nature.com/naturecommunications 17