Article Identification of constrained sequence elements across 239 primate genomes https://doi.org/10.1038/s41586-023-06798-8 Noncoding DNA is central to our understanding of human gene regulation and Received: 9 February 2023 complex diseases1,2, and measuring the evolutionary sequence constraint can establish the functional relevance of putative regulatory elements in the human genome3–9. Accepted: 30 October 2023 Identifying the genomic elements that have become constrained specifically in Published online: xx xx xxxx primates has been hampered by the faster evolution of noncoding DNA compared Open access to protein-coding DNA10, the relatively short timescales separating primate species11, Check for updates and the previously limited availability of whole-genome sequences12. Here we construct a whole-genome alignment of 239 species, representing nearly half of all extant species in the primate order. Using this resource, we identified human regulatory elements that are under selective constraint across primates and other mammals at a 5% false discovery rate. We detected 111,318 DNase I hypersensitivity sites and 267,410 transcription factor binding sites that are constrained specifically in primates but not across other placental mammals and validate their cis-regulatory effects on gene expression. These regulatory elements are enriched for human genetic variants that affect gene expression and complex traits and diseases. Our results highlight the important role of recent evolution in regulatory sequence elements differentiating primates, including humans, from other placental mammals. Functional genomic elements that have acquired selective constraint regulatory sequences in the human genome. Using comparisons with specific to the primate order are prime candidates for understanding other mammals, we identify an important class of noncoding regulatory the evolutionary changes that have contributed to the uniqueness elements with constraint specific to primates and delineate a role for of our own species13–16. Whereas comparisons between the human these elements in human health by integrating functional genomics genome and those of other mammal and vertebrate species have and population genetics datasets. revealed an extensive catalogue of constrained genes and regulatory elements4–6,17,18, identifying constrained sequence elements that are specific to primates has been particularly challenging owing to the A 239-way primate whole-genome alignment short evolutionary distances separating these species5,18. Compared To identify genomic elements with primate-specific constraint, we with the mammalian lineage, which includes more than 6,000 species constructed a multiple sequence alignment that densely samples the separated by more than 100 million years of evolution19, the primate primate lineage. We identified 187 primate species without an available order only consists of approximately 500 species that are separated reference assembly that had recently reported Illumina whole-genome by a fraction of this time11—around 65 million years. Thus, despite 43 sequencing data11,23, and assembled their genomes using Megahit25 primate species having been aligned in the recent Zoonomia study20 of based on an average coverage of 35× per individual. We combined the 240 placental mammals, the total phylogenetic branch length within resulting contigs together with 52 previously published high-quality these primates is only around 10% that of the placental mammal align- primate reference assemblies to create a reference-free whole-genome ment21. At such short timescales, it is unclear whether the absence of MSA of 239 primate species with Cactus21 (Supplementary Data 1). This genetic changes between species is owing to functional constraints, alignment represents all major primate lineages, including 86% of gen- or simply because insufficient time has passed for random mutations era and all 16 families (Fig. 1a,b). As our goal was to quantify sequence to arise. Consequently, the selective constraints specific to the phy- constraint across the human genome, we confirmed that each base was logenetic branch from which the human species ultimately emerged covered by an average of 174 other primate species, and 85% percent of remain largely unidentified. the euchromatic regions of the human genome were covered by at least We recently reported a catalogue of genetic diversity in primates 100 other primate species (Fig. 1c). To ensure that the per-base error rate based on hundreds of species and individuals, which enabled us to in our de novo assemblies was sufficiently low for subsequent constraint gain insight into evolutionary and population dynamics in the primate analysis, we compared a set of 25 species within our data for which both order11,22. Leveraging the vast new catalogue of benign missense muta- newly generated short-read contigs and previously published reference tions in these species, we further developed and applied models to iden- genomes were available. We found that the rates of mismatches between tify pathogenic variants in protein-coding sequences, which account for these assembly pairs ranged between 0.02 and 0.5% and were largely only 1% of the human genome23,24. Here, we expand on these prior works explained by differences in the species’ heterozygosity (Fig. 1d and by constructing a genome-wide multiple sequence alignment (MSA) Supplementary Table 1). After accounting for intraspecific variation, of 239 primate species to better characterize constraint at noncoding the average remaining mismatch rate attributable to assembly and A list of authors and their affiliations appears at the end of the paper. Nature | www.nature.com | 1 Article sequencing errors was reduced to 0.04% (Methods). Finally, we gener- with shared gene-regulatory functions33,34, the majority of human DNase ated a 441-species mammalian MSA by combining our primate MSA with I hypersensitivity site (DHS) elements and transcription factor bind- the remaining mammalian orders sampled in Zoonomia20. This consti- ing or occupancy sites (TFBSs) currently lack detectable sequence tutes the deepest species sampling for mammals in a whole-genome constraint35,36. This lack of observed constraint in non-primate ances- MSA to date, including 204 primate species unique to this study, and tors might reflect a true divergence in function at these elements, but enables detection of sequence constraint both broadly across mammals could also be owing to recently acquired sequence constraint in the and in the more recent evolution of our own lineage. primate order37. We estimated the average sequence constraint for primates and mammals in high-resolution maps of 1.2 million DHS elements from Primate-constrained protein-coding sequences 438 cell types8 (Methods). At an FDR of 5%, we observed that 35% and Expanding the number of available primate species in the MSA to 239 33% of elements exhibited evidence of constraint across mammals or increased the phylogenetic branch length 2.8-fold over the previously within primates, respectively, and largely overlapped (Supplementary available 43 primate species alignment from the Zoonomia study20. Data 3, OR = 14.1, P < 1.0 × 10−300). After removing DHS elements with We used phyloP26 to estimate genome-wide per-base constraint for ambiguous or contradictory evidence of constraint (Methods), we regions of the MSA without ambiguous alignments and found that 3.1% observed that 42% had evidence of sequence constraint in species of the bases in the human genome were nominally constrained across that had diverged over 100 million years ago (Ma) (42%), and 111,318 all primates (phyloP score > 1.3 or P < 0.05), compared with 7.1% of bases (11%) were significantly constrained in primates but lacked evidence that were constrained in the broader set of 240 mammals at the same of constraint in mammals or vertebrates (Fig. 2a,b, Extended Data thresholds. We additionally detected 157 Mb of constrained sequence Fig. 4a,b and Methods). The identification of these elements was elements in the primate order using phastCons26, comprising 5.1% of the largely consistent regardless of constraint metric (phyloP or phast- human genome. To determine whether primate constraint metrics could Cons, OR of overlap = 12.7, P < 1.0 × 10−300), and sensitivity analyses distinguish functional from neutral sequence, we investigated con- suggested that the identification of primate-specific DHS elements straint scores in annotated sequence elements. First, we observed that was robust to mammalian FDR thresholds, regional differences in protein-coding DNA—including exons, start codons and stop codons— mutation rates and effects of incomplete lineage sorting (Extended was strongly enriched in phastCons elements (Fig. 1e). Noncoding DNA Data Fig. 4c–f). encompassing transcribed regions and cis-regulatory elements (CREs) Within these DHS elements, transcription factor occupancy prevents in accessible chromatin or occupied by a transcription factor was also DNase I cleavage to create footprints of transcription factor binding significantly enriched. We observed periodic patterns of codon con- events at nucleotide resolution8,38. Across 3.6 million TFBS footprints, straint that differentiate exonic from surrounding intronic sequences at we find that 1,034,832 (30%) have evidence of broad constraint in mam- the nucleotide level (Fig. 1e). Primate phyloP also distinguished between mals, whereas 267,410 (8%) show primate-specific constraint (Extended non-synonymous and fourfold degenerate sites, although less well than Data Fig. 5 and Supplementary Data 4). Consistent with previous work, mammal phyloP, which is better powered, given the higher total branch a substantial fraction of footprintable regulatory elements exhibited length in the mammal MSA (Extended Data Figs. 1 and 2). complex architecture (37%) and contain multiple TFBSs with differing We next explored whether we could identify protein-coding genes evolutionary constraints on their binding sequences39 (Methods). Of and exons that are constrained specifically in primates but not in other note, 66% of DHS elements with primate-specific constraint have a TFBS placental mammals27. We estimated primate and non-primate mammal with evidence of constraint in mammals, suggesting that regulatory sequence constraint in canonical protein-coding exons annotated in the function initially evolved in a common ancestor (Fig. 2c). However, human genome, identifying 179,329 exons with evidence of constraint in 19% of mammal-constrained DHS elements contain individual TFBS primates at a false discovery rate (FDR) of 5%. As expected, 99% of these footprints with evidence of primate-specific constraint, suggesting exons were broadly constrained across non-primate mammals and ver- that the function of deeply constrained elements can further evolve. tebrates, but 2,178 were constrained specifically in primates (Extended Furthermore, we find evidence that the number of DHS elements with Data Fig. 3a,b). The majority of primate-constrained exons (72%) are primate-specific constraint is likely to be underestimated by phyloP annotated as protein-coding at orthologous regions in the mouse owing to short branch lengths, including 208,717 DHS elements genome, indicating that they are not newly evolved coding sequences with primate-specific constraint detectable only by phastCons and but instead have been subject to shifts in selective constraint in the an additional 86,987 unconstrained DHS elements with at least one primate order. Genes that had at least one exon constrained among primate-specific TFBS. Overall, we find that a significant fraction of primates but none across other mammals (Supplementary Data 2) putative human CREs have evidence of constraint in primates but not were most highly enriched for involvement in the antibacterial humoral in mammals or vertebrates. response (fold enrichment = 26.4, P = 1.8 × 10−9; Supplementary Table 2). We undertook several studies to validate the biological function The overall structure and splicing of these genes were broadly con- of these putative regulatory elements with evidence of constraint strained across mammals, suggesting that the amino acid sequences specific to the primate order using orthogonal computational and that they encode may have become constrained early on in primate experimental approaches. First, we investigated whether they were evolution as a maintained response to pathogens. Primate-specific more likely to have a regulatory function in humans than elements with- constrained exons were also significantly more likely to undergo alterna- out detectable constraint. Broadly constrained and primate-specific tive splicing (P = 1.3 × 10−7) and had lower levels of transcript inclusion constrained elements had higher chromatin accessibility and were (P = 8.6 × 10−6; Extended Data Fig. 3c,d), hinting at an initially limited accessible in significantly more cell types than unconstrained ele- utilization of recently evolved exons28–31. Our results underscore that the ments (P < 1.0 × 10−300 for both; Fig. 2d). Across massively parallel evolution of new protein-coding genes or exons from existing sequences reporter assays40 (MPRAs) of 148 cis-regulatory sequence elements, is rare, whereas the increased functional importance of pre-existing both mammal and primate constraint at the nucleotide level were exons is a relatively more common, but still infrequent, event32. significantly correlated with transcriptional changes in saturation mutagenesis experiments (49% and 35%, respectively), of which 14% correlated with primate constraint only (Fig. 2e and Supplementary Primate-constrained CREs Data 7). Since elements with primate-specific constraint appeared Although comparative genomic and epigenomic studies of mammals to have more cell-type-specific biochemical activity than broadly and other vertebrates have identified many CREs in the human genome constrained elements, we also tested whether the extent of primate 2 | Nature | www.nature.com a b Position (Gb) c 3.0 Aligned bases (×109) 0 0.5 1.0 1.5 2.0 2.5 2.5 1 2.0 2 1.5 Cercopithecidae (n = 85) 3 1.0 4 5 0.5 6 0 0 50 100 150 200 7 Number of species Hominidae (n = 7) 8 9 d Hylobatidae (n = 10) Pearson r = 0.91 Mismatch rate (×10–3) 10 5 Chromosome Callitrichidae (n = 23) 11 4 12 Aotidae (n = 5) 3 13 Cebidae (n = 11) 2 14 15 1 Atelidae (n = 15) 16 1 2 3 4 5 17 Pairwise heterozygosity (×10–3) Pitheciidae (n = 28) 18 19 Number of aligned species Tarsiidae (n = 4) 20 0 121 239 e 50 Fraction of Lemuridae (n = 19) human genome 21 40 0.1 Indriidae (n = 10) Enrichment (OR) Cheirogaleidae (n = 4) 22 0.3 X 30 0.5 Lepilemuridae (n = 4) Daubentoniidae (n = 1) Y Lorisidae (n = 8) 20 Galagidae (n = 5) f 10 Primate phyloP 1.5 0 1.0 op od ) Ex co on s n U ) 3′ TR fo D TR Pr otp HS TF Tra mo nts oc nsc ters pa ts y St rt c DS 5′ (all nc on do cu rip 0.5 U ri a C St s ( 0 o on 0 25 50 75 100 125 150 175 Ex S DH Base pair Fig. 1 | MSA of 239 primate species. a, Cladogram of primate species included with a corresponding 95% confidence interval ribbon is shown. e, Enrichment in the MSA. The number of sampled species per family is given in parenthesis. of primate phastCons elements for coding and noncoding genomic b, Ideogram of the human genome depicting the average number of species elements. The size of the circle represents the fraction of the human genome. covered by the MSA at 500-kb resolution. Telomeric, centromeric and The dashed grey line indicates an odds ratio (OR) of 1. CDS, coding sequence; heterochromatic regions (light blue) are indicated. c, Cumulative primate TF, transcription factor; UTR, untranslated region. (f) Codon periodicity in species coverage of the human genome in the 239-way primate MSA. d, Per-base the mean primate phyloP scores across 482 protein-coding exons exactly mismatch rate between newly generated short-read contigs and species with 130 nucleotides in length. Coding sequences are shown in dark blue and previously published high-quality reference assemblies. A linear regression fit flanking intronic sequences in beige. constraint at an element was consistent with cell-type-specific regula- deeply constrained and primate-specific constrained elements were tory activity using Enformer41, a deep-learning method that predicts significantly better correlated between human and macaques than gene expression from sequence without using sequence constraint. at elements without evidence of constraint (P = 0.0004 and 0.0001, Across 438 cell types, we observed that primate constraint correlated respectively; Fig. 2i), indicating that constraint on the sequence level better with estimates of gene-regulatory activity when the element was corresponds to constraint of molecular function between species. Nev- accessible in similar cell-type categories to the Enformer predictions ertheless, primate-specific constrained elements also shared functional (Fig. 2f). Together, these results indicate that regulatory elements with similarity between primates and mouse, consistent with the results of evidence of sequence constraint specific to primates have important our TFBS analyses. cis-regulatory functions in humans. Evolutionary constraint estimated in mammals and vertebrates is cor- In addition to the extensive body of human experimental data pro- related with selective constraint estimated in human populations17,44, viding support for the function of primate-constrained regulatory so we explored contemporary human cohorts for evidence of ongoing elements, a limited number of experiments have been conducted in selection against genetic variants that disrupt primate-constrained non-human primates, enabling us to investigate the regulatory activity regulatory elements. Using the gnomAD cohort of 141,456 human indi- of primate-constrained DHS elements in non-human contexts. First, viduals45, we found that predicted target genes of primate-specific we set out to experimentally validate the regulatory capacity of a small elements had significantly fewer loss-of-function mutations than subset of DHS elements with primate-specific constraint. We cloned expected (P < 10−300; Fig. 3a). Moreover, we observed increased muta- orthologous sequences from human, chimpanzee and mouse into tional constraint46 in the noncoding primate-specific constrained ele- luciferase reporter assays, transfected these constructs into human ments themselves (P < 10−300; Fig. 3b). Indeed, polymorphic variants in induced pluripotent stem cells (iPS cells), and measured transcription regulatory elements were more likely to have allele-specific regulatory of the reporter gene for three elements. Of note, two out of three ele- effects by MPRA when there was evidence of constraint in primates at ments drove transcription more strongly from the primate sequences the mutated nucleotide (P = 0.0007) or across the entire regulatory than from the mouse sequence (Fig. 2g and Supplementary Data 6), and element (P = 2.9 × 10−13; Fig. 3c), even after controlling for mammalian we set out to validate this observation more broadly. We investigated constraint (P = 1.1 × 10−5). Together, these results extend previous stud- chromatin accessibility across DHS elements in fibroblasts from four ies44,46 and suggest that regulatory elements constrained specifically non-human primate species, observing that primate-specific con- in the primate order are under purifying selection in human popula- strained DHS elements displayed higher and more consistent chromatin tions and that mutations in these elements are likely to have important accessibility in all four primate species compared to unconstrained regulatory functions. DHS elements42 (Fig. 2h and Extended Data Fig. 6a). We also investi- To explore whether genes expressed in specific tissues were more gated the levels of H3K27ac, a marker of active CREs, in stage-matched likely to be regulated by noncoding elements with primate-specific cell types during corticogenesis at orthologous regions in humans, constraint, we investigated the depth of conservation across 16 rhesus macaques and mice43. We observed that H3K27ac levels at broadly defined cellular contexts47. We confirmed that regulatory Nature | www.nature.com | 3 Article No. of primate- a DHS GRIA4 b c specific TFBSs 15 0.5 3+ Human Mammal- 2 Proportion of DHS elements 0 constrained 1 0.4 ATAC-seq 15 Chimp 0.3 0 42% 30 Mammal Primate Mouse 0.2 0 47% 11% 0.1 2 Primate- specific No evidence 0 PhyloP 0 of constraint constraint co M ain c of o e ined l- tra ce 6 tr ifi a ns am t t ns en tra m ns ec in co vid co -sp e 0 at N im Pr TEAD4 motif d No. of accessible cell types Human GCGTGGAGCAGAGTGAGCATTCCAGAGTCCCAGACAAGAGGAG 1.35 70 Chimpanzee Accessibility (CPM) GCGTGGAGCAGAGTGAGCATTCCAGAGTCCCAGACAAGAGGAG Gibbon GCGTGGGGCAGAGTGAGCATTTCAGAGTCCCAGACAAGAGGAG 1.30 60 Rhesus monkey GCGTGGAGCAGAGTGAGCATTCCAGAGTCCCAGACAAGAGGAG Squirrel monkey GCGCGGAGCGGAGTGAGCATTCCAGAGTCCCGGACGAGAGGAG 1.25 Pacific walrus ACATGGAG--------------------CCGGGCTAAGAAAAG Saiga antelope ACATGGAGTAGAGTGAGTATTCTGGAGCTCCAGTTAAGCAAAG 50 Black flying fox ACATGGAGTAGAGTGAGTATTCCAGCGTC-CAGCTTAAAAAAA 1.20 Mouse ACTGTAAGTAG---GAGTTTTCCAAA--CCCCAAGC--AAAGG 40 1.15 Primate phyloP e r = 0.28, FDR = 0.0004 1 g 0 1.5 Mouse Primates Average substitution effect 3 6 0 *** *** NS Relative luciferase reporter activity –1 2 4 1 Mammal phyloP r = 0.08, FDR = 0.18 1 1 2 0 3 0 0 0 0 –1 f h 3.3 Accessibility (log2 CPM) 3.2 7.5 3.1 Normalized Spearman’s U Kidney 5.0 3.0 2.9 2.5 Brain 2.8 0 Stomach –2.5 Human Chimp Gorilla Orangutan Macaque Primate phyloP Muscle –5.0 i H3K27ac (Spearman's U) Heart Lung 0.7 0.6 Multiple 0.5 Human vs Human vs Macaque vs macaque mouse mouse Enformer prediction Fig. 2 | Identification of noncoding regulatory elements with primate- (Spearman’s ρ) between primate phyloP and sequence-based Enformer specific constraint. a, Example of a primate-specific constrained DHS element predictions of regulatory activity across 438 ENCODE cell types. Categories of in the GRIA4 locus (hg38; chromosome (chr.) 11:105608279–105612792). Assay similar cell types corresponding to specific tissues are indicated. g, Normalized for transposase-accessible chromatin with sequencing (ATAC-seq) insertions luciferase reporter activity in human iPS cells for three selected sets of primate- from human, chimpanzee and mouse iPS cells and phyloP constraint in primates specific constrained DHS elements at orthologous primate and mouse and mammals. A putative TEAD4 binding motif that better matches primate sequences. Colours indicate constraint categories from b. Bars represent mean sequences than non-primate mammal sequences is indicated. b, Proportion and error bars represent 95% confidence intervals; n = 36 across 3 elements. of constrained DHS elements across clades. c, Number of primate-specific P values: 1.4 × 10 −5 (left), 2.8 × 10 −4 (middle) and 0.54 (right). Raw data are constrained footprints (TFBSs) in DHS elements, stratified by constraint across provided in Supplementary Data 6. h, Average chromatin accessibility in the entire DHS. Error bars represent 95% confidence intervals. d, Average fibroblasts for five primate species at orthologous sequence elements chromatin accessibility and the number of accessible cell types is higher at stratified by sequence constraint. Colours indicate constraint categories from more constrained DHS elements. Colours indicate constraint categories from b. Error bars represent 95% confidence intervals; n = 90,827 DHS elements. b. Error bars represent 95% confidence intervals. CPM, counts per million. i, Average Spearman ρ of H3K27ac levels at orthologous CREs for three pairs e, A saturation mutagenesis experiment (MPRA) of a DHS element at of species. Colours indicate constraint categories from b. Error bars represent chr. 2:191049304–191045304 (hg38). Average effects of substitutions at each 95% confidence intervals. n = 12 for human versus mouse, n = 10 for all other nucleotide on transcriptional activity are correlated with phyloP scores from comparisons. ***P < 0.001; NS, not significant. primates but not from mammals. f, Heat map of normalized correlation values 4 | Nature | www.nature.com elements active in multiple cell types—and particularly in neural a 0.64 Mammal-constrained b and musculoskeletal cell types—were most deeply constrained48, Loss of function intolerance in human exomes (LOEUF) 0.8 Primate-specific constraint whereas blood, epithelial, and placental cell types were least con- Non-coding mutation intolerance in human genomes (z-score) 0.68 Primate-specific strained (Fig. 3d). Regulatory elements present in neural, cardiac constraint 0.7 Mammal-constrained and embryonic cell types exhibited higher phyloP scores in primates 0.6 than in mammals (Fig. 3e). We explore the connection between ultra- 0.72 0.5 No evidence conserved elements (UCEs) and neural cell types below. Finally, we No evidence of constraint of constraint investigated whether specific TFBSs were more or less constrained 0.76 0.4 in primates than in mammals, finding that most TFBS motifs in DHS c e Neural footprints had significant, but small, differences (241 out of 282 (85%); 1.5 Fig. 3f). A small number of footprints are over 20% less constrained in 1.4 Common variant allelic effect in MPRA (OR) Cardiac primates than mammals, including the KRAB zinc-finger domain tran- Primate phyloP 1.0 Organ scription factors (KZNFs), ZNF384 and ZNF28. The reduced constraint 1.2 Embryonic development at KZNF binding sites in primates probably reflects the divergence 0.5 of KZNFs themselves, which are among the fastest evolving gene 1.0 0 Cancer/epithelial families in primates49,50. PhyloP PhyloP Phast- Placental (DHS) (nucleotide) Cons 0.5 1.0 1.5 d Mammal phyloP Ultraconserved elements in primates Proportion of DHS elements 0.6 In addition to the elements that we detected as constrained by phy- loP and phastCons, we identified 74.6 million positions in the human 0.4 Mammal-constrained genome that are perfectly conserved without a single substitution Primate-specific constraint across all 239 primate species. These positions were often contigu- 0.2 No evidence of constraint ous, and we catalogued 33,368 primate UCEs that were at least 20 bps in length (Supplementary Data 5), amounting to more than 1 Mb of 0 ue y l Va Org ari . sc an ant /e . St eu . m l lo tro iac ry l B Em ges d na ry e Ly anc c C mp er ce oid ac i. l C lA ss ar ta ro ra ta in dev ar ev N do total DNA sequence including 7,261 coding exons and 22,582 DHS Pl r/ep Re b tiv i i l/c n Di hro Ti on ele en a /e a ye S ard ul d o an h n id m t lm sk v Pu ulo elements. More than half (57%) of the 4,552 recently described mam- c us M malian UCEs18 overlapped our primate UCEs, and 82% overlapped after M allowing for up to 1% of missing species per aligned column within f Not significant the primate alignment. Genes whose protein-coding sequences over- FDR < 0.05 0 lapped primate UCEs were more likely to be involved in nervous sys- footprint constraint (log2 OR) tem development (Supplementary Table 3, fold enrichment = 2.24, Priamte vs mammal P = 8.8 × 10−9). We additionally found that 2.7% of primate UCEs also –0.5 overlapped brain regulatory elements (fold enrichment = 3.1, P < 10−300), FOX Homeodomain consistent with the deep constraint of neuronal protein-coding SOX HOMEZ ZNF28 sequences. –1.0 LIN54 MYB IRF CUX Complex trait variation in constrained CREs –1.5 ZNF384 Genome-wide association studies (GWAS) have identified hundreds Fig. 3 | Characterization of constrained regulatory elements. a, Predicted of thousands of genetic variants associated with complex human dis- target genes have fewer loss-of-function mutations in humans than expected eases and changes in gene expression, the majority of which map to at constrained DHS elements. Error bars represent 95% confidence intervals. noncoding CREs27,33,34,37. We identified DHS elements and footprints b, Constrained DHS elements have fewer mutations in human populations containing fine-mapped GWAS variants (posterior inclusion probability than unconstrained elements. Error bars represent 95% confidence intervals. c, Enrichment of allele-specific regulatory activity (MPRA) for 27,023 common (PIP) > 0.5) for 96 human clinical phenotypes and complex traits from variants, stratified by type of constraint. A colour legend for constraint the UK Biobank8,47, and characterized whether the underlying sequence categories is shown in d. Error bars represent 95% confidence intervals, the was constrained only in primates (65 Ma), placental mammals (100 Ma), central dot represents point estimates; n = 27,023 variants. d, Proportion vertebrates (160–400 Ma), or without evidence of constraint (less than of constrained DHS elements across 16 broad cellular contexts. Error bars 65 Ma; Fig. 4a and Extended Data Fig. 6c). Fine-mapped variants under- represent 95% confidence intervals, centre represents proportion. n = 1,029,688 lying clinical phenotypes and complex traits were enriched across all DHS elements. Dev., development; endo., endothelial; epi., epithelial. e, Scatter classes of distal accessible chromatin element and footprints, includ- plot of mean primate and mammal phyloP scores at DHS elements, stratified by ing those with primate-specific constraint (OR = 2.4, P = 2.5 × 10−13 and cell types. A linear fit is shown with a corresponding 95% confidence interval OR = 4.0, P = 1.8 × 10−7, respectively), with more deeply constrained ele- ribbon. Putative outlier cell types with higher primate phyloP than mammal ments showing greater enrichment51. A heritability enrichment analysis phyloP scores are indicated. f, Differences in the proportion of primate and corroborated the relevance of constrained regulatory elements and mammalian constrained footprints in human DHS elements, for each of 283 primate-specific constraint more generally in complex traits (Extended transcription factor family motifs. Positive values indicate a higher proportion Data Fig. 6d). In comparison, fine-mapped variants underlying changes of constrained TFBSs in primates, negative values indicate a lower proportion in gene expression (expression quantitative trait loci (eQTLs)) from of constrained TFBSs in primates. Transcription factors that are the least constrained in primates compared to mammals are labelled, and significantly the GTEx study showed similar enrichment for elements with recent different transcription factors are coloured in magenta (FDR < 5%). Error bars constraint but were markedly less enriched at elements that are broadly represent 95% confidence intervals. constrained across mammals or vertebrates. After stratifying human genes by selective constraint quantified by loss-of-function observed/ expected upper bound fraction (LOEUF) scores38, we found that vari- (OR = 4.6, P = 1.0 × 10−53 and OR = 8.0, P = 4.3 × 10−24, respectively), ants affecting the expression of highly constrained genes tended to whereas variants affecting the expression of less constrained genes be enriched at more deeply constrained DHS elements and footprints tended to reside at elements with more recent constraint (Fig. 4b). Nature | www.nature.com | 5 Article a eQTL (GTEx v8) d r 2 to rs686030 Complex traits (UKBB) HDL-C (0.8,1] (0.2,0.4] –log10(P value) 6 (0.6,0.8] [0,0.2] 40 (0.4,0.6] Sentinel Noncoding fine-mapped 20 variant enrichment 4 0 Posterior probability 500 kb 1.0 rs686030 2 (PIP) 0.5 0 0 PSIP1 Split time (Ma) 160–400 100 65 <65 Non- FREM1 TTC39B CCDC171 accessible Accessible chromatin chromatin SNAPC3 Human selective constraint b 0–20% (Most constrained genes) 500 kb 20–40% B cell (CD20+) rs686030 40–60% 2 DHS or ChIP (CPM) 6 0 60–80% 2 Monocyte (CD14+) 80–100% (Least constrained genes) Noncoding fine-mapped 0 variant enrichment 2 Small intestine 4 0 2 HepG2 (liver carcinoma) 0 1 CEBPα 0 2 300 bp Primate 2 0 0 Mammal 6 PhyloP Split time (Ma) 160–400 100 65 <65 Non- 0 accessible Vertebrate Accessible chromatin 6 chromatin 0 c 1.00 613 variants 795 variants CEBPα Constraint motif <65 Ma rs686030 65 Ma C fine-mapped variants 0.75 100 Ma Human GGTTGACAGACTTGATCAATCAAGGTCAA 160–400 Ma GGTTGACAGACTTGATCGATCAAGGTCAA Proportion of Chimpanzee Western hoolock gibbon GGTTGACAGACTTGATCAATCAAGGTCAA 0.50 L’Hoests monkey AGTTGACAGACTTGATCAATCAAGGTCAA Ecuadorian squirrel monkey GGTTAAGAGACTTGATCAATCAAGGTCAA Brown rat GGCTA-G--ACTTATC-AGTCAAGGTCAT 0.25 Pacific walrus GAATGAGAGACTTGC----TCAAAATTAC Black flying fox GGTTGAGAGATCAGCTCAATCAAGGTCAT 0 Saiga antelope GGTTGAAAGACTGGGTTAGTCAGGGTAAC Coding Accessible exons chromatin Fig. 4 | Enrichment of complex trait variants at constrained noncoding CREs. genes separated into five bins based on their selective population constraint a, Enrichment of fine-mapped GWAS variants from 96 UK Biobank (UKBB) (LOEUF). Ribbons represent 95% confidence intervals and the centre represents complex traits and clinical phenotypes (red) or eQTLs for 49 GTEx tissues the point estimate. c, Total count of fine-mapped variants for 96 UK Biobank (blue) in DHS elements, stratified by sequence constraint of the element. phenotypes in protein-coding exons or accessible chromatin sites, stratified Approximate split times for vertebrates (160–400 Ma), placental mammals by extent of constraint as in a. d, Example of a fine-mapped variant (rs686030) (100 Ma) and primates (65 Ma) are shown. Enrichments are computed as the for HDL-C and cholelithiasis at a primate-specific constrained DHS element. ratio of the proportion of variants with PIP > 0.5 compared to the proportion of GWAS signal at the locus, fine-mapping probability, DNase signal, CEBPα variants with PIP < 0.01. Ribbons represent 95% confidence intervals and the chromatin immunoprecipitation with sequencing (ChIP–seq) signal, constraint centre represents the point estimate. The grey dashed line indicates an OR of 1. scores and MSAs of primate (blue) and mammal (green) species are shown. b, Enrichment of fine-mapped eQTL variants within DHS elements as in a, with To explore the functional role of primate-specific constrained CREs bile cholesterol secretion. Although 36% of fine-mapped variants at in human complex traits and clinical phenotypes, we partitioned the DHS elements lack significant constraint across primates and other fine-mapped variants from the UK Biobank by protein-coding con- mammals, these elements were also not significantly enriched for sequence and constraint depth. In contrast to 88% of fine-mapped heritability in humans (Extended Data Fig. 6d), suggesting that further protein-coding variants residing within deeply constrained exons data are needed to resolve these loci, some of which might be false that predate the emergence of placental mammals (Fig. 4c and Sup- positives53. Of note, we find residual enrichment for fine-mapped plementary Data 8), only 37% of noncoding variants in accessible variants in DHS elements that lack evidence of constraint by phyloP chromatin were constrained to this extent. 12% of fine-mapped vari- (FDR < 5%) but overlap with phastCons elements in primates (Extended ants in CREs were constrained only in primates and not in placental Data Fig. 6f). Additional sequencing to increase sampling density on mammals, corresponding to 93 probably causal regulatory variants this branch may help to define the selective constraints at the origin of underlying human complex traits and clinical phenotypes (Supple- our own species and their contribution to human clinical phenotypes mentary Data 9 and 10). One example is rs686030, a fine-mapped and complex diseases. noncoding variant in a primate-constrained DHS element near the TCC39B gene, which is associated with high-density lipoprotein (HDL) cholesterol levels (PIP = 0.99) and cholelithiasis (PIP = 0.38) (Fig. 4d). Discussion The derived allele strengthens a motif for the bound CEBPα transcrip- Heritable modifications in genomic sequence are necessary for trait tion factor and is associated with TCC39B gene expression (PIP = 0.43 adaptations and the emergence of new species, but the nature of for liver), and mouse knockout studies of TTC39B showed an increase these sequence changes remains incompletely understood. Although in HDL-C levels52, potentially modulating the risk of cholelithiasis via constrained noncoding elements in mammals have been extensively 6 | Nature | www.nature.com catalogued, less attention has been paid to those in the primate line- 13. Boffelli, D. et al. Phylogenetic shadowing of primate sequences to find functional regions of the human genome. Science 299, 1391–1394 (2003). ages, in part owing to the challenges in detecting constraint at short 14. Gilad, Y., Oshlack, A., Smyth, G. K., Speed, T. P. & White, K. P. Expression profiling in phylogenetic distances with previously available species sampling. primates reveals a rapid evolution of human transcription factors. Nature 440, 242–245 By placing the genomes of 239 primate species, including 187 newly (2006). 15. Orkin, J. D., Kuderna, L. F. K. & Marques-Bonet, T. The diversity of primates: from assembled here, in the context of other mammalian and vertebrate biomedicine to conservation genomics. Annu. Rev. Anim. Biosci. 9, 103–124 (2021). genomes20, we identified hundreds of thousands of constrained 16. Sousa, A. M. M., Meyer, K. A., Santpere, G., Gulden, F. O. & Sestan, N. Evolution of the noncoding sequence elements and catalogued the origins of their human nervous system function, structure, and development. Cell 170, 226–247 (2017). 17. Lindblad-Toh, K. et al. A high-resolution map of human evolutionary constraint using 29 sequence constraint in primates, placental mammals and more distant mammals. Nature 478, 476–482 (2011). vertebrates. Collectively, these CREs are unique evolutionary records 18. Christmas, M. J. et al. Evolutionary constraint and innovation across hundreds of placental that provide a lens through which to view the mechanisms of recent mammals. Science 380, eabn3943 (2023). 19. Wilson, D. E. & Reeder, D. M. Mammal Species of the World: A Taxonomic and Geographic exaptations leading to our species10. Reference (JHU Press, 2005). In keeping with prior work showing that noncoding DNA evolves 20. Zoonomia Consortium. A comparative genomics multitool for scientific discovery and more rapidly than protein-coding sequences17,18,54,55, we find that many conservation. Nature 587, 240–245 (2020). 21. Armstrong, J. et al. Progressive Cactus is a multiple-genome aligner for the thousand- human CREs that previously showed no evidence of sequence con- genome era. Nature 587, 246–251 (2020). straint are in fact constrained exclusively in primates, considerably 22. Sørensen, E. F. et al. Genome-wide coancestry reveals details of ancient and recent expanding the number of known constrained noncoding elements male-driven reticulation in baboons. Science 380, eabn8153 (2023). 23. Gao, H. et al. The landscape of tolerated genetic variation in humans and primates. in the human genome. Indeed, sequence constraint in primates Science 380, eabn8153 (2023). uniquely predicted the function of a subset of regulatory elements, 24. Fiziev, P. P. et al. Rare penetrant mutations confer severe risk of common diseases. and specifically constrained elements had higher and more similar Science 380, eabo1131 (2023). 25. Li, D., Liu, C.-M., Luo, R., Sadakane, K. & Lam, T.-W. MEGAHIT: an ultra-fast single-node regulatory functions in diverse human cell types and across distinct solution for large and complex metagenomics assembly via succinct de Bruijn graph. primate species. These elements are predicted to regulate genes Bioinformatics 31, 1674–1676 (2015). that are more intolerant to deleterious mutations in human popula- 26. Pollard, K. S., Hubisz, M. J., Rosenbloom, K. R. & Siepel, A. Detection of nonneutral substitution rates on mammalian phylogenies. Genome Res. 20, 110–121 (2010). tions and are significantly enriched for common genetic variants 27. Frankish, A., Diekhans, M., Jungreis, I. & Lagarde, J. GENCODE 2021. Nucleic Acids Res. associated with variation in gene expression and complex human 49, D916–D923 (2021). 28. Pan, Q. et al. Alternative splicing of conserved exons is frequently species-specific in traits and diseases. Nevertheless, some functional genomic elements human and mouse. Trends Genet. 21, 73–77 (2005). underlying complex human phenotypes do not show evidence of 29. Merkin, J., Russell, C., Chen, P. & Burge, C. B. Evolutionary dynamics of gene and isoform constraint in either primates or mammals in our analysis, suggesting regulation in Mammalian tissues. Science 338, 1593–1599 (2012). 30. Xiong, J. et al. Predominant patterns of splicing evolution on human, chimpanzee and that they potentially emerged after the initial radiation of primates macaque evolutionary lineages. Hum. Mol. Genet. 27, 1474–1485 (2018). and thus became selectively constrained only in a sub-lineage such 31. Suntsova, M. V. & Buzdin, A. A. Differences between human and chimpanzee genomes as anthropoids or apes, or that functional sequence elements were and their implications in gene expression, protein functions and biochemical properties of the two species. BMC Genomics 21, 535 (2020). selectively lost in one or more lineages. Additional sequencing of the 32. Kondrashov, F. A. & Koonin, E. V. Origin of alternative splicing by tandem exon duplication. remaining species in the primate order, including population-level Hum. Mol. Genet. 10, 2661–2669 (2001). oversampling of key lineages, would help to provide the resolution 33. Mikkelsen, T. S. et al. Comparative epigenomic analysis of murine and human adipogenesis. Cell 143, 156–169 (2010). needed to detect sequence elements under selective constraint in 34. Odom, D. T. et al. Tissue-specific transcriptional regulation has diverged significantly finer detail, especially those specific to clades from which the human between human and mouse. Nat. Genet. 39, 730–732 (2007). species ultimately emerged. 35. Ward, L. D. & Kellis, M. Evidence of abundant purifying selection in humans for recently acquired regulatory functions. Science 337, 1675–1678 (2012). 36. Necsulea, A. & Kaessmann, H. Evolutionary dynamics of coding and non-coding transcriptomes. Nat. Rev. Genet. 15, 734–748 (2014). Online content 37. Villar, D. et al. Enhancer evolution across 20 mammalian species. Cell 160, 554–566 (2015). Any methods, additional references, Nature Portfolio reporting summa- 38. Vierstra, J. et al. Global reference mapping of human transcription factor footprints. ries, source data, extended data, supplementary information, acknowl- Nature 583, 729–736 (2020). edgements, peer review information; details of author contributions 39. Fong, S. L. & Capra, J. A. Modeling the evolutionary architectures of transcribed human enhancer sequences reveals distinct origins, functions, and associations with human trait and competing interests; and statements of data and code availability variation. Mol. Biol. Evol. 38, 3681–3696 (2021). are available at https://doi.org/10.1038/s41586-023-06798-8. 40. Kircher, M. et al. Saturation mutagenesis of twenty disease-associated regulatory elements at single base-pair resolution. Nat. Commun. 10, 3583 (2019). 41. Avsec, Ž. et al. Effective gene expression prediction from sequence by integrating 1. Albert, F. W. & Kruglyak, L. The role of regulatory variation in complex traits and disease. long-range interactions. Nat. Methods 18, 1196–1203 (2021). Nat. Rev. Genet. 16, 197–212 (2015). 42. Edsall, L. E. et al. Evaluating chromatin accessibility differences across multiple primate 2. Lappalainen, T. & MacArthur, D. G. From variant to function in human disease genetics. species using a joint modeling approach. Genome Biol. Evol. 11, 3035–3053 (2019). Science 373, 1464–1468 (2021). 43. Reilly, S. K. et al. Evolutionary genomics. Evolutionary changes in promoter and enhancer 3. Dermitzakis, E. T. & Clark, A. G. Evolution of transcription factor binding sites in activity during human corticogenesis. Science 347, 1155–1159 (2015). mammalian gene regulatory regions: conservation and turnover. Mol. Biol. Evol. 19, 44. Drake, J. A. et al. Conserved noncoding sequences are selectively constrained and not 1114–1121 (2002). mutation cold spots. Nat. Genet. 38, 223–227 (2006). 4. Thomas, J. W. et al. Comparative analyses of multi-species sequences from targeted 45. Karczewski, K. J. et al. The mutational constraint spectrum quantified from variation in genomic regions. Nature 424, 788–793 (2003). 141,456 humans. Nature 581, 434–443 (2020). 5. Boffelli, D., Nobrega, M. A. & Rubin, E. M. Comparative genomics at the vertebrate 46. Chen, S. et al. A genome-wide mutational constraint map quantified from variation in extremes. Nat. Rev. Genet. 5, 456–465 (2004). 76,156 human genomes. Preprint at bioRxiv https://doi.org/10.1101/2022.03.20.485034 6. Margulies, E. H. et al. Analyses of deep mammalian sequence alignments and constraint (2022). predictions for 1% of the human genome. Genome Res. 17, 760–774 (2007). 47. Meuleman, W. et al. Index and biological spectrum of human DNase I hypersensitive sites. 7. Sullivan, P. F. et al. Leveraging base-pair mammalian constraint to understand genetic Nature 584, 244–251 (2020). variation and human disease. Science 380, eabn2937 (2023). 48. Cardoso-Moreira, M. et al. Gene expression across mammalian organ development. 8. Roadmap Epigenomics Consortium. Integrative analysis of 111 reference human Nature 571, 505–509 (2019). epigenomes. Nature 518, 317–330 (2015). 49. Pontis, J. et al. Primate-specific transposable elements shape transcriptional networks 9. ENCODE Project Consortium. An integrated encyclopedia of DNA elements in the human during human development. Nat. Commun. 13, 7178 (2022). genome. Nature 489, 57–74 (2012). 50. Nowick, K. et al. Gain, loss and divergence in primate zinc-finger genes: a rich resource 10. King, M. C. & Wilson, A. C. Evolution at two levels in humans and chimpanzees. Science for evolution of gene regulatory differences between species. PLoS ONE 6, e21553 (2011). 188, 107–116 (1975). 51. Vierstra, J. et al. Mouse regulatory DNA landscapes reveal global principles of 11. Kuderna, L. F. K. et al. A global catalog of whole-genome diversity from 233 primate cis-regulatory evolution. Science 346, 1007–1012 (2014). species. Science 380, 906–913 (2023). 52. Teslovich, T. M. et al. Biological, clinical and population relevance of 95 loci for blood 12. Juan, D., Santpere, G., Kelley, J. L., Cornejo, O. E. & Marques-Bonet, T. Current advances in lipids. Nature 466, 707–713 (2010). primate genomics: novel approaches for understanding evolution and disease. Nat. Rev. 53. Cui, R. et al. Improving fine-mapping by modeling infinitesimal effects. Preprint at bioRxiv Genet. 24, 314–331 (2023). https://doi.org/10.1101/2022.10.21.513123 (2022). Nature | www.nature.com | 7 Article 54. Hardison, R. C. et al. Covariation in frequencies of substitution, deletion, transposition, National Parks, Arusha, Tanzania. 19North Carolina Museum of Natural Sciences, Raleigh, NC, and recombination during eutherian evolution. Genome Res. 13, 13–26 (2003). USA. 20Department of Biological and Biomedical Sciences, North Carolina Central University, 55. Kryukov, G. V., Pennacchio, L. A. & Sunyaev, S. R. Most rare missense alleles are Durham, NC, USA. 21Department of Biological Sciences, North Carolina State University, deleterious in humans: implications for complex disease and association studies. Am. J. Raleigh, NC, USA. 22Department of Evolutionary Anthropology, Duke University, Durham, NC, Hum. Genet. 80, 727–739 (2007). USA. 23Renaissance Computing Institute, University of North Carolina at Chapel Hill, Chapel Hill, NC, USA. 24Copenhagen Zoo, Frederiksberg, Denmark. 25Universidade Federal de Viçosa, Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in Viçosa, Brazil. 26Universidade Federal do Amazonas, Departamento de Genética, Laboratório published maps and institutional affiliations. de Evolução e Genética Animal (LEGAL), Manaus, Brazil. 27Department of Anthropology, University of Utah, Salt Lake City, UT, USA. 28Universidade Federal do Para, Bragança, Brazil. 29 Research Group on Terrestrial Vertebrate Ecology, Mamirauá Institute for Sustainable Open Access This article is licensed under a Creative Commons Attribution Development, Tefé, Brazil. 30Rede de Pesquisa em Diversidade, Conservação e Uso da Fauna 4.0 International License, which permits use, sharing, adaptation, distribution da Amazônia – RedeFauna, Manaus, Brazil. 31Comunidad de Manejo de Fauna Silvestre en la and reproduction in any medium or format, as long as you give appropriate Amazonía y en Latinoamérica—ComFauna, Iquitos, Peru. 32Universidade Federal de Rondônia, credit to the original author(s) and the source, provide a link to the Creative Commons licence, Porto Velho, Brazil. 33Instituto Nacional de Pesquisas da Amazônia, Manaus, Brazil. 34Instituto and indicate if changes were made. The images or other third party material in this article are de Biociências, Universidade Federal do Mato Grosso, Cuiabá, Brazil. 35Department of Biology, included in the article’s Creative Commons licence, unless indicated otherwise in a credit line Trinity University, San Antonio, TX, USA. 36Life Sciences and Environment, Technology and to the material. If material is not included in the article’s Creative Commons licence and your Environment of Mahajanga, University of Mahajanga, Mahajanga, Madagascar. 37Department intended use is not permitted by statutory regulation or exceeds the permitted use, you will of Anthropology, New York University, New York, NY, USA. 38Department of Neuroscience, need to obtain permission directly from the copyright holder. To view a copy of this licence, Washington University School of Medicine in St Louis, St Louis, MO, USA. 39Keeling visit http://creativecommons.org/licenses/by/4.0/. Center for Comparative Medicine and Research, MD Anderson Cancer Center, Bastrop, TX, USA. 40Department of Anthropology, Yale University, New Haven, CT, USA. 41School of © The Author(s) 2023 Interdisciplinary Forensics, Arizona State University, Phoenix, AZ, USA. 42California National Primate Research Center, University of California, Davis, CA, USA. 43Guinea Worm Eradication Program, The Carter Center Ethiopia, Addis Ababa, Ethiopia. 44State Key Laboratory of Genetic Resources and Evolution, Kunming Institute of Zoology, Chinese Academy of Lukas F. K. Kuderna1,84, Jacob C. Ulirsch1,84, Sabrina Rashid1,84, Mohamed Ameen1, Sciences, Kunming, China. 45Center for Evolutionary and Organismal Biology, Zhejiang Laksshman Sundaram1, Glenn Hickey2, Anthony J. Cox1, Hong Gao1, Arvind Kumar1, University School of Medicine, Hangzhou, China. 46Villum Centre for Biodiversity Genomics, Francois Aguet1, Matthew J. Christmas3, Hiram Clawson2, Maximilian Haeussler2, Section for Ecology and Evolution, Department of Biology, University of Copenhagen, Mareike C. Janiak4, Martin Kuhlwilm5,6, Joseph D. Orkin7, Thomas Bataillon8, Copenhagen, Denmark. 47Liangzhu Laboratory, Zhejiang University Medical Center, Shivakumara Manu9,10, Alejandro Valenzuela11, Juraj Bergman8,12, Marjolaine Rouselle8, Hangzhou, China. 48Women’s Hospital, School of Medicine, Zhejiang University, Hangzhou, Felipe Ennes Silva13,14, Lidia Agueda15, Julie Blanc15, Marta Gut15, Dorien de Vries4, China. 49Tanzania Wildlife Research Institute (TAWIRI), Arusha, Tanzania. 50Institute of Ian Goodhead4, R. Alan Harris16, Muthuswamy Raveendran16, Axel Jensen17, International Animal Health/One Health, Friedrich-Loeffler-Institut, Federal Research Institute Idriss S. Chuma18, Julie E. Horvath19,20,21,22,23, Christina Hvilsom24, David Juan11, for Animal Health, Greifswald–Insel Riems, Germany. 51Professorship for International Animal Peter Frandsen24, Joshua G. Schraiber1, Fabiano R. de Melo25, Fabrício Bertuol26, Health/One Health, Faculty of Veterinary Medicine, Justus Liebig University, Giessen, Hazel Byrne27, Iracilda Sampaio28, Izeni Farias26, João Valsecchi29,30,31, Malu Messias32, Germany. 52Department of Environmental Ecology, Faculty of Environmental Sciences, Maria N. F. da Silva33, Mihir Trivedi10, Rogerio Rossi34, Tomas Hrbek26,35, University of Science and Central Institute for Natural Resources and Environmental Studies, Nicole Andriaholinirina36, Clément J. Rabarivola36, Alphonse Zaramody36, Clifford J. Jolly37, Vietnam National University, Hanoi, Vietnam. 53Institut Català de Paleontologia Miquel Jane Phillips-Conroy38, Gregory Wilkerson39, Christian Abee39, Joe H. Simmons39, Crusafont, Universitat Autònoma de Barcelona, Barcelona, Spain. 54Department of Zoology, Eduardo Fernandez-Duque40, Sree Kanthaswamy41,42, Fekadu Shiferaw43, Dongdong Wu44, State Museum of Natural History Stuttgart, Stuttgart, Germany. 55Institució Catalana de Long Zhou45, Yong Shao44, Guojie Zhang44,45,46,47,48, Julius D. Keyyu49, Sascha Knauf50,51, Recerca i Estudis Avançats (ICREA), Barcelona, Spain. 56Centre for Genomic Regulation (CRG), Minh D. Le52, Esther Lizano11,53, Stefan Merker54, Arcadi Navarro11,55,56,57,58, Tilo Nadler59, The Barcelona Institute of Science and Technology, Barcelona, Spain. 57Barcelonaβeta Brain Chiea Chuen Khor60, Jessica Lee61, Patrick Tan60,62,63, Weng Khong Lim62,63,64, Research Center, Pasqual Maragall Foundation, Barcelona, Spain. 58Universitat Pompeu Andrew C. Kitchener65,66, Dietmar Zinner67,68,69, Ivo Gut15, Amanda D. Melin70,71,72, Fabra, Barcelona, Spain. 59Cuc Phuong Commune, Nho Quan District, Vietnam. 60Genome Katerina Guschanski17,73, Mikkel Heide Schierup8, Robin M. D. Beck4, Ioannis Karakikes74,75, Institute of Singapore, Agency for Science, Technology and Research, Singapore, Singapore. Kevin C. Wang76,77,78, Govindhaswamy Umapathy9,10, Christian Roos79, Jean P. Boubli4, 61 Mandai Nature, Singapore, Singapore. 62SingHealth Duke–NUS Institute of Precision Adam Siepel80, Anshul Kundaje81,82, Benedict Paten2, Kerstin Lindblad-Toh3,83, Medicine (PRISM), Singapore, Singapore. 63Cancer and Stem Cell Biology Program, Duke-NUS Jeffrey Rogers16 ✉, Tomas Marques Bonet11,15,53,55,58 ✉ & Kyle Kai-How Farh1 ✉ Medical School, Singapore, Singapore. 64SingHealth Duke–NUS Genomic Medicine Centre, Singapore, Singapore. 65Department of Natural Sciences, National Museums Scotland, Edinburgh, UK. 66School of Geosciences, Edinburgh, UK. 67Cognitive Ethology Laboratory, 1 Illumina Artificial Intelligence Laboratory, Illumina, San Diego, CA, USA. 2UC Santa Cruz Germany Primate Center, Leibniz Institute for Primate Research, Göttingen, Germany. Genomics Institute, University of California, Santa Cruz, CA, USA. 3Science for Life Laboratory, 68 Department of Primate Cognition, Georg-August-Universität Göttingen, Göttingen, Department of Medical Biochemistry and Microbiology, Uppsala University, Uppsala, Sweden. Germany. 69Leibniz ScienceCampus Primate Cognition, Göttingen, Germany. 70Department 4 School of Science, Engineering and Environment, University of Salford, Salford, UK. of Anthropology and Archaeology, University of Calgary, Calgary, Alberta, Canada. 5 Department of Evolutionary Anthropology, University of Vienna, Vienna, Austria. 6Human 71 Department of Medical Genetics, University of Calgary, Calgary, Alberta, Canada. Evolution and Archaeological Sciences (HEAS), University of Vienna, Vienna, Austria. 72 Alberta Children’s Hospital Research Institute, University of Calgary, Calgary, Alberta, 7 Département d’Anthropologie, Université de Montréal, Montréal, Quebec, Canada. Canada. 73Institute of Ecology and Evolution, School of Biological Sciences, University of 8 Bioinformatics Research Centre, Aarhus University, Aarhus, Denmark. 9Academy of Scientific Edinburgh, Edinburgh, UK. 74Cardiovascular Institute, Stanford University, Stanford, CA, USA. and Innovative Research (AcSIR), Ghaziabad, India. 10Laboratory for the Conservation of 75 Department of Cardiothoracic Surgery, Stanford University, Stanford, CA, USA. 76Department Endangered Species, CSIR-Centre for Cellular and Molecular Biology, Hyderabad, India. of Cancer Biology, Stanford University, Stanford, CA, USA. 77Department of Dermatology, 11 IBE, Institute of Evolutionary Biology (UPF-CSIC), Department of Medicine and Life Sciences, Stanford University School of Medicine, Stanford, CA, USA. 78Veterans Affairs Palo Alto Universitat Pompeu Fabra, Barcelona, Spain. 12Section for Ecoinformatics and Biodiversity, Healthcare System, Palo Alto, CA, USA. 79Gene Bank of Primates and Primate Genetics Department of Biology, Aarhus University, Aarhus, Denmark. 13Research Group on Primate Laboratory, German Primate Center, Leibniz Institute for Primate Research, Göttingen, Biology and Conservation, Mamirauá Institute for Sustainable Development, Tefé, Brazil. Germany. 80Simons Center for Quantitative Biology, Cold Spring Harbor Laboratory, Cold 14 Evolutionary Biology and Ecology (EBE), Département de Biologie des Organismes, Spring Harbor, NY, USA. 81Department of Computer Science, Stanford University, Stanford, Université libre de Bruxelles (ULB), Brussels, Belgium. 15Centro Nacional de Analisis Genomico CA, USA. 82Department of Genetics, Stanford University, Stanford, CA, USA. 83Broad Institute (CNAG), Barcelona, Spain. 16Human Genome Sequencing Center and Department of of MIT and Harvard, Cambridge, MA, USA. 84These authors contributed equally: Lukas F. K. Molecular and Human Genetics, Baylor College of Medicine, Houston, TX, USA. 17Department Kuderna, Jacob C. Ulirsch, Sabrina Rashid. ✉e-mail: jr13@bcm.edu; tomas.marques@upf.edu; of Ecology and Genetics, Animal Ecology, Uppsala University, Uppsala, Sweden. 18Tanzania kfarh@illumina.com 8 | Nature | www.nature.com Methods accounting for the fraction of the genome within alignment blocks, resulting in mismatch rates that range from 0.00026–0.00515 mis- De novo assembly and repeat-masking matches per bp. As the genome assemblies produced herein are haploid To maximize the species diversity of primates in our analyses, we newly compressions of diploid organisms, a random allele will be sampled sequenced and assembled the genomes of 187 different primate spe- and incorporated at heterozygous positions, and thus the resulting cies, initially presented in refs. 11,23, for which no other reference differences between two assemblies of the same species should be genome assembly was available. In brief, each individual was sequenced strongly correlated with the species’ intraspecific diversity. We com- with 150 bp paired end reads on the Illumina NovaSeq 6000 platform pared our mismatch rates to the estimates of heterozygosity for the to an average whole-genome coverage of ~35×, and we assembled the same genomes presented in ref. 11, and confirmed that heterozygosity resulting reads into contigs using Megahit25 (version 1.2.9) using default accounts for 83% of the observed variation in mismatch rates across parameters. The resulting assemblies had an average contig N50 of assemblies. We quantified the residual mismatch rate after regressing 34 kb, and the assembly sizes ranged from 2.1–3.0 Gb, thus falling within out it’s the effects of heterozygosity, and found the resulting average the typical range of previously reported genome sizes for primates56 mismatch rate to be 0.0004 mismatches per bp, which we consider to (see Extended Data Fig. 1a). We then combined these assemblies with be sufficiently low for our analyses. We note that the number of base the reference genomes of 52 additional species that had been previously differences due to assembly error is likely lower than this, as residual generated as part of other studies57 and or available through public mismatches also include fixed differences between individuals, which repositories (Supplementary Data 1). The final species sampling densely are not accounted for by heterozygosity. covers the whole primate radiation and includes members of all 16 primate families and 72 primate genera. We identified and soft-masked Detecting selective constraint common genomic repeats within the assemblies using RepeatMasker We measured selective constraint genome wide using the widely used (version 4.1.2-p1; http://www.repeatmasker.org), utilizing the primates phyloP and phastCons algorithms from the PHAST package26,59. To do repeat catalogue as query. so, we extracted the ancestral genomes of primates and of eutherian mammals from our alignment using haltools hal2fasta, and annotated Multiple sequence alignment common genomic repeats in both using ReapeatMasker as described We aligned the assemblies with Cactus21 (version 2.1.1), using the phy- above, but using the mammalian repeat catalogue for the eutherian logeny presented in11 as a guide tree for progressive decomposition, ancestor. We lifted the resulting annotations into human reference and used the previously available high-quality assemblies as alignment space, and randomly sampled 1 Mb of autosomal short interspersed outgroups. All computation was done by running cactus-prepare with nuclear element (SINE), long interspersed nuclear element (LINE), options–wdl–noLocalInputs–preprocessBatchSize 5–defaultDisk long terminal repeat (LTR) and DNA repeats from the alignments as 3000 G–halAppendDisk 9000 G–defaultCores 64–gpu–gpuCount putatively neutrally evolving regions. We used these regions as input for 8–defaultMemory 385 G–alignMemory 450 to produce a script in work- phyloFit together with the general reversible model (“–subst-mod REV”) flow description language (WDL), then uploading it to Terra (https:// as the nucleotide substitution model and expectation maximization app.terra.bio/) where it was executed on Google Cloud Platform. algorithm (“-EM”) to fit it to the data. As our goal is to detect elements GPU-related issues prevented that version of Cactus from executing to with sequence constraint specific to primates, we generated the neutral completion, so the job was resumed using a WDL made without the–gpu background models once for all primates, and once for all mammals and–gpuCount options. An outgroup to primates (Mus musculus refer- after excluding the primate branch. We additionally generated a neutral ence mm10) was manually added to the root alignment job by editing model for the 100-way vertebrate MSA from UCSC Genome Browser the WDL, and the ‘LOCAL’ disk parameter of the hal_append_subtree in our analysis to minimize false negatives on the mammalian track, task was manually increased to 9,000. Cactus has since been fixed for which we also excluded the primate branch containing 11 species (v2.2.3) to resolve all issues encountered during this alignment. and defined neutral background models via alignments at 4D sites as We then combined our resulting primate MSA with the recently putatively neutral regions, due to their easier detection across the much generated mammalian MSA by the Zoonomia consortium20. In brief, larger phylogenetic distances present in this alignment. we used hal2fasta from the haltools21 package to output the ancestral We used the models to estimate constraint in different ways across genome at the root of the primate MSA, and used it to generate a bridge the three clades (primates, mammals, vertebrates): For phyloP, we alignment with the Sunda colugo (Galeopterus variegatus), the closest calculated scores for both constraint and acceleration with the “–mode outgroup to primates in the Zoonomia MSA. We used this bridge align- CONACC” flag, and used the likelihood ratio test “–method LRT” yield- ment to insert the primate MSA into the Zoonomia MSA, and replace ing phyloP scores—that is, the −log10(P value) from the hypothesis the original primate branch with it. test, and the associated scale factor. We scored individual bases by To generate the final, filtered alignment used as input for subsequent outputting them via the “–wig-scores” flags. We additionally scored analyses described below, we output maf files centred on the human element-wide annotations for coding sequences, DHS and TFBS by genome reference using haltools including the “–onlyOrthologs– passing them to phyloP via the “–features” flag, to increase power as noAncestors –noDupes” flags, thus removing any regions with poten- the test is performed across more than a single basepair. Finally, we tially ambiguous mappings at multiple locations. generated discrete constrained elements in primates using phastCons, using primate neutral background model, the “–expected-length 45– Pairwise alignments error rate estimate target-coverage 0.3–rho 0.31” consistent with previous studies18, and To quantify residual error rates within the genome assemblies gen- output constrained elements with the “–most-conserved” flag. erated in this project, we identified 25 species for which a reference To explore the potential impact of regional variation in substitution genome was previously assembled with an orthogonal, state of the rates on our estimates of constraint, we additionally generated regional art combination of technologies (Supplementary Table 1). After intro- neutral background models for primates and other mammals from 1-Mb ducing a minimum contig length cutoff of 1 kb, we generated pairwise sliding windows across the human genome. In each window, we subset alignments between the two assemblies using minimap258 (v. 2.17-r941) the previously identified ancestral repeats and randomly selected using the following flags:–cs -x asm5. We called variants on the resulting 100 kb of sequence after trimming sites with >20% missing data. As alignments by retaining alignment blocks of at least 1 kb within the PAF described above, these sites were used to estimate substitution rates file using paftools.js, by applying the following flags: paftools.js call -l input with phyloFit, and the resulting models were used to run phyloP 1000 -L 1000. We quantified mismatch rates from the resulting output for individual bases and DHSs elements. Article To additionally ensure our estimates of constraint are robust to topo- in our primary analyses38. Cell types and tissues where each DHS ele- logical variation in the underlying phylogeny due to potential sources ment was most strongly associated were previously estimated using of uncertainty such as incomplete lineage sorting, we additionally non-negative matrix factorization with 16 components47. inferred regional phylogenies for primates using a maximum likeli- We defined a core 40-bp window surrounding the summit of the hood approach implemented in IQtree. In brief, we randomly subset peak of each DHS annotation as the input to calculate element-wise. 150 kb of trimmed sequence from each 1 Mb window, which was used Analogous to protein-coding exons, we then calculated constraint in to estimate an appropriate substitution model and infer the phylogeny DHS and TFBS element-wise using phyloP across primates, mammals, including 1,000 bootstraps. We used the topology of the resulting and vertebrates, and define constrained elements in each clade as consensus tree and the ancestral repeat alignments to infer neutral those remaining significant at a 5% FDR60. To define primate-specific models as described, using the same subset of sites as for the regional constraint in DHS and TFBS, we required the elements to be significantly models to minimize additional sources of variation, and assessed the constrained in primates, but not in mammals or vertebrates. Finally, concordance of constraint for DHS elements between regional models DHS elements and TFBSs that did not have primate-specific constraint using the canonical and regional phylogenies. by phyloP but overlapped with a primate PhastCons elements were excluded from the primary analyses for consistency in interpretation, Protein-coding exons since these sequences represent a mixture of primate-specific and To identify protein-coding exons with constrained specifically in the deeply constrained sequences. The depth of constraint for each DHS primate lineage, we used phyloP with protein-coding exons from and TFBS are provided in Supplementary Data 9 and 10. Approximate GENCODE (v 42)9,27 as element-wise input as described above across target genes of each DHS element were based on the closest gene using the primate, mammalian, and vertebrate tracks. We restricted these the ‘nearest’ function the R GenomicRanges package. analyses to exons that are part of ‘Ensembl canonical’ transcript, and additionally excluded any exon that overlaps known human segmental TFBS enrichment analysis duplications, as defined by the segmental duplication track on UCSC We obtained archetypal motifs overlapping each TFBS or DHS footprint Genome Browser. We ran element-wise phyloP tests on these remain- from the annotations presented in ref. 38. Footprints typically had mul- ing coding exons, and defined constrained exons for each clade (pri- tiple motif matches and were considered independently. For each motif, mates, mammals, vertebrates) directly based on the resulting P values. we computed the proportion of footprints in either constraint category We accounted for multiple testing by retaining those that remained (primate or mammal constrained below an FDR of 5%, as described significant at a 5% FDR60. To define exons with primate-specific con- above), where the denominator was the total number of constrained straint, we required them to be significantly constrained in primates, footprints (primate or mammal) regardless of motif match. We then but not in mammals or vertebrates. To detect whether these exons calculated the odds ratio for each motif to test whether the proportion also have coding potential in other mammals, we lifted the underly- of primate-constrained and mammal-constrained footprints were dif- ing coordinates to the mouse genomes (mm10) and checked weather ferent. After observing a small bias where short footprints were more they overlap protein-coding annotations there. To define genes with likely to be detected as constrained in mammals, we split footprints primate-specific constraint, we looked for genes containing one or into 10 equal-sized bins, computed the odds ratio for each motif in more exons with primate-specific constraint, but no mammal differ- each bin, then performed a fixed effects meta-analysis for each motif. entially constrained ones. To calculate differences in the proportion of alternatively spliced exons between broadly constrained and primate Primate UCEs specifically constrained exons, we calculated the mean exon inclu- We defined UCEs across primates analogous to ref. 18: We filtered sion rate across tissues from the GTEx project61, and defined exons regions with ambiguous or multiple alignments using haltools includ- with an inclusion rate different from 1 as alternatively spliced. A list ing the “–onlyOrthologs–noAncestors –noDupes” flags, and parsed the of exons and genes with primate-specific constraint is presented in resulting alignment to exclude any alignment column that is different Supplementary Data 2. from all other species in at least one species. We then kept consecu- tive stretches of 20 bp or more for the final set of UCEs. For a more lax GO-term enrichment definition, we allowed for missing data (“-” or “N”) in the alignment in We used Panther62 to calculate Gene Ontology (GO)-term enrichments at most 2 species (1%). We strictly defined overlap to previous annota- of genes with primate-specific constraint, and those overlapping pri- tions as 1 bp or more. mate UCEs. We used Fishers’ exact to test for statistical overrepresenta- tion on the ‘GO biological process’ annotation, by using the Ensembl Estimates of constraint in human populations identifiers of the underlying genes from either analysis as foreground Gene constraint in the human population was estimated using the set and the human gene annotation as background. To account for LOEUF metric. In brief, this metric conservatively estimates the selec- multiple testing, we report only results that remain significance at a tion against loss-of-function mutations by taking the upper bound of FDR (Benjamini–Hochberg) of 5%. a 95% Poisson confidence interval around the observed to expected ratio of loss-of-function mutations. LOEUF values were obtained DHSs and TFBSs from 141,456 individuals in gnomAD v245. Constraint across noncod- We obtained high-resolution maps of DHSs from 733 human biosamples ing regions of the genome was estimated as a z-score for depletion of encompassing 438 cell and tissue types and states47. The study reported mutations compared to expectation46. Z-scores for non-overlapping 3.6 million DHS elements, and we applied several additional quality 1,000-bp bins were obtained from 71,156 individuals in gnomAD v3. control steps to remove low-quality peaks. First, we excluded all peaks When a DHS element overlapped multiple bins the average z-score without 1-to-1 matches between GRCh38 and hg19. We normalized peaks was used. to 300 bps in size for all analyses, except for the element-wise constraint scoring described below. Finally, we required all peaks to be within Trait-associated variant analyses the top 100,000 in at least one annotated cell type in the datasets, by Fine-mapping results for 96 complex traits and diseases across 366,194 the normalized score provided from the study. After excluding sex unrelated ‘white British’ individuals in the UKBB63 were obtained from chromosomes, this resulted in a set of 1,238,405 peaks that were used https://www.finucanelab.org/data and have previously been described in downstream analyses. We similarly obtained 3,622,316 consensus in detail64. In brief, fine mapping was performed using FINEMAP65,66 DNase I hypersensitivity footprints for the set of DHS elements used and SuSiE67 with GWAS summary statistics from SAIGE/BOLT-LMM and in-sample dosage linkage disequilibrium (LD) computed by LDstore 268. primate-specific constraint for all bio-sample pairs, and computed the Regions were defined by expanding ±1.5 Mb for each lead variant and correlation between the Enformer and phyloP scores for the retained were merged if they overlapped. Up to 10 causal variants were allowed regions, and row and column normalized the final correlation matrix. per region. Posterior inclusion probabilities (PIPs) were averaged across The final matrix was hierarchically clustered on the rows, and the same the two methods and variants where PIPs from the two methods dis­ order was retained for the columns in the heat map. Major cell types agreed by >0.05 were excluded. for each correlation block identified are highlighted as annotations. Fine-mapping results for expression quantitative traits in 49 tissues across 838 individuals were obtained from https://www.finucanelab. Luciferase reporter vector construction org/data and have been described in detail61,64. In brief, fine mapping Mouse, chimp and human CRE with 150 bp in length were synthesized by was performed using SuSiE on cis-eQTL summary statistics from the IDT. The CRE was cloned into the linearized pGL3- Promoter vector (cut GTEx portal (https://gtexportal.org/). Covariates (sex, PCR amplifica- by Nhel and BglII). The fusion product (pGL3-cRE) was subsequently tion, sequencing platform, genotype principal components, and proba- transformed into Mix & Go Competent Cells Strain Zymo 5-a (Zymo bilistic estimation of expression residuals factors69) were projected Research, T3007). Clones were selected by ampicillin and plasmids were out from the genotypes prior to fine mapping. After fine mapping, all prepared using the NucleoSpin Plasmid Transfection-grade (Takara, variants were lifted over from GRCh38 to hg19. 740490). Definition of constraint at DHS and TFBSs was slightly modified such that evidence of constraint out to mammals or vertebrates was Transfection and luciferase assays separated and elements with discrepant estimates of constraint were Human iPS cells were transfected in a 24-well plate using the Lipo- excluded. Specifically, constraint at approximately 100 Ma required fectamine Stem Transfection Reagent (Invitrogen, STEM00001) that mammal and primate phyloP scores were below the FDR thresh- and Opti-MEM Reduced Serum medium (Invitrogen, 31-985- 070). old but vertebrate phyloP was above the FDR threshold. Similarly, On the day of transfection, cell density was 50% confluent. For each constraint at approximately 160–400 Ma required that vertebrate, well, 500 ng of pGL3-enhancer, pGL3-control, or pGL3-promoter was mammal, and primate phyloP scores were below the FDR threshold. co-transfected with 10 ng of pRL-CMV (Promega, E2261) as an inter- Bigwig files for accessible chromatin and transcription factor occu- nal control for the normalization of luciferase activity. Cells were pancy were obtained from the ENCODE project47,70 (ENCFF220IWU, incubated with DNA–lipid complex overnight and media was changed ENCFF659BVQ, ENCFF619LIB and ENCFF842XRQ) or the Sequence for another two days. The firefly and Renilla luciferase activity were Read Archive (SRX097095). Coding variants were annotated as measured respectively using a Dual-Glo Luciferase Assay System loss-of-function, missense, or synonymous using the Ensembl Vari- (Promega, E2920). Human iPS cells were obtained from the Stanford ant Effect Predictor (VEP) v8571. When a variant had multiple coding CVI iPS cell Biobank. annotations, the most severe consequence on the canonical transcript (GENCODE v19) was used. Massively parallel reporter assays We computed the enrichment of fine-mapped variants for different Measured effects of single nucleotide substitution effects from satu- annotations by comparing the proportion of variants with PIP > 0.5 to ration mutagenesis experiments across 29 regulatory elements were the proportion of variants with PIP < 0.01. Distal elements were defined obtained from40 and across 131 elements from9. For each nucleotide, the as DHS elements that did not overlap promoters72. When variants were mean substitution effect across all reported nucleotides was correlated fine-mapped across multiple traits, tissues, or genes, only the highest (Pearson) with phyloP scores that were truncated such that negative PIP variant was used. Confidence intervals and P values were estimated values, which are indicative of possible acceleration, were set to zero. A using Fisher’s exact test. Enrichments were performed in hg19 and Storey FDR60 was used to control for multiple comparisons. Regulatory annotations were lifted over from GRCh38. effects from 27,017 common variants in the DHS elements investigated A similar enrichment analysis was performed using stratified LD score in this study were obtained from9. Variants with a reported FDR below regression (S-LDSC)72 to estimate the heritability in each annotation. 5% were defined as allele-specific. A generalized linear model with a Similar to previous studies7, S-LDSC models were fit using approxi- binomial probability distribution was used to estimate the effects of mately 10 million common variants including the Baseline v2.2 anno- constraint on allele-specific activity. tations. Annotations derived in GrCH38 were lifted over to hg19, and their LD scores were estimated using the EUR sub-population of the Chromatin accessibility and histone modifications in non-humans 1000 Genomes project. Enrichment and average per-SNP heritability Chromatin accessibility from ATAC-seq in fibroblasts obtained from estimates were meta-analysed across 69 mostly independent traits human and 4 non-human primates (chimpanzee, gorilla, orangutan using a random effects model. and macaque) at 89,744 merged peaks with orthologous sequences in The predicted effects of fine-mapped variants on transcription fac- all 5 species were obtained from42,75. Counts were transformed to log2 tor binding was estimated using motifbreakR73 for 426 position weight counts per million (cpm), and FDR values from differential accessibility matrices from HOCOMOCOv1174. A motif match was determined using testing across any primate species were obtained42. the information content (ic) if either allele obtained a P value < 0.0001. Histone modifications (H3K27ac) were also obtained from three A variant disrupted a motif match if there was a difference of >0.4 for matching cell types during corticogenesis for human, macaque, and the scaled motif matrix between alleles. mouse43. First, H3K27ac peaks at orthologous sequences from all species were obtained from the authors and filtered such that at least Enformer analysis 200 bp of these peaks overlapped with a DHS element in this study. We obtained the 733 bio-sample aggregated DNase peak dataset as Next, DHS elements coordinates in GRCh38 were lifted over to each curated by47 and deduplicated the technical replicates by retaining species and the maximum H3K27ac signal (cpm) at each element was the top bio-sample for samples with technical replicates. We retained calculated using the provided bigwig files. Spearman correlations all DHS peaks found in more than two biosamples for downstream between matching cell types were then computed for each pair of spe- analysis, calculated the midpoint for each DHS and scored the regions cies stratified by the type of constraint on the DHS element. using the Enformer model41. To assess the local functional relevance of the Enformer scores, we averaged them across ±128 bp around the mid- Reporting summary point of each DHS. To compute the correlation between the Enformer Further information on research design is available in the Nature Port- score and phyloP in each bio-sample, we pairwise intersected DHS with folio Reporting Summary linked to this article. Article F.E.S. was supported by Brazilian National Council for Scientific and Technological Development Data availability (CNPq) (Processes 303286/2014-8, 303579/2014-5, 200502/2015-8, 302140/2020-4, 300365/2021-7, 301407/2021-5 and 301925/2021-6), by the Fonds de la Recherche Scientifique Primate assemblies have been deposited at the European Nucleotide - FNRS (#40017464), and by the European Union’s Horizon 2020 research and innovation Archive (ENA) under the accession PRJEB67744. The MSA and constraint programme under the Marie Skłodowska-Curie grant agreement no. 801505. Fieldwork for samples collected in the Brazilian Amazon was funded by grants from Conselho Nacional de tracks are available through the UCSC Genome Browser. Desenvolvimento Científico e Tecnológico (CNPq/SISBIOTA Program #563348/2010-0), Fundação de Amparo à Pesquisa do Estado do Amazonas (FAPEAM/SISBIOTA #2317/2011), and 56. Kuderna, L. F., Esteller-Cucala, P. & Marques-Bonet, T. Branching out: what omics can tell Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES AUX # 3261/2013) to us about primate evolution. Curr. Opin. Genet. Dev. 62, 65–71 (2020). I.F. Samples from Amazônia, Brazil were accessed under SisGen no. A8F3D55. Sampling of 57. Shao, Y. et al. Phylogenomic analyses provide insights into primate evolution. Science non-human primates in Tanzania was funded by the German Research Foundation (KN1097/3-1 380, 913–924 (2023). to S. Knauf and RO3055/2-1 to CR). No animals in Tanzania were sampled purposely for this 58. Li, H. New strategies to improve minimap2 alignment accuracy. Bioinformatics 37, study. Details of the original study on Treponema pallidum infection can be requested from 4572–4574 (2021). S. Knauf. This research was funded by the Vietnamese Ministry of Science and Technology’s 59. Hubisz, M. J., Pollard, K. S. & Siepel, A. PHAST and RPHAST: phylogenetic analysis with Program 562 (grant no. ĐTĐL.CN-64/19). A.N. is supported by AEI-PGC2018-101927-BI00 704 space/time models. Brief. Bioinform. 12, 41–51 (2011). (FEDER/UE), FEDER (Fondo Europeo de Desarrollo Regional)/FSE (Fondo Social Europeo), 60. Storey, J. D. A direct approach to false discovery rates. J. R. Stat. Soc. B 64, 479–498 Unidad de Excelencia María de Maeztu, funded by the AEI (CEX2018-000792-M) and Secretaria (2002). d’Universitats i Recerca and CERCA Programme del Departament d’Economia i Coneixement 61. The GTEx Consortium. et al. The GTEx Consortium atlas of genetic regulatory effects de la Generalitat de Catalunya (GRC 2017 SGR 880). A.D.M. was supported by the National across human tissues. Science 369, 1318–1330 (2020). Sciences and Engineering Research Council of Canada and Canada Research Chairs 62. Thomas, P. D. et al. PANTHER: Making genome-scale phylogenetics accessible to all. program. Aotus azarae samples from Argentina were obtained with grant support to E.F.-D. Protein Sci. 31, 8–22 (2022). from the Zoological Society of San Diego, Wenner-Gren Foundation, the L. S. B. Leakey 63. Bycroft, C. et al. The UK Biobank resource with deep phenotyping and genomic data. Foundation, the National Geographic Society, the US National Science Foundation Nature 562, 203–209 (2018). (NSF-BCS-0621020, 1232349, 1503753, 1848954; NSF-RAPID-1219368, NSF-FAIN-1952072; 64. Kanai, M. et al. Insights from complex trait fine-mapping across diverse populations. NSF-DDIG-1540255; NSF-REU 0837921, 0924352, 1026991), and the US National Institute on Preprint at bioRxiv https://doi.org/10.1101/2021.09.03.21262975 (2021). Aging (NIA- P30 AG012836-19, NICHD R24 HD-044964-11). J.H.S. was supported in part by the 65. Benner, C. et al. FINEMAP: efficient variable selection using summary data from genome- NIH under award number P40OD024628 - SPF Baboon Research Resource. K.G. was supported wide association studies. Bioinformatics 32, 1493–1501 (2016). by the Swedish Research Council VR (2020-03398). We thank C. Escudé and B. Bed’Homme, 66. Benner, C., Havulinna, A. S., Salomaa, V., Ripatti, S. & Pirinen, M. Refining fine-mapping: and L. Cacheux and J.-P. Gautier for providing guenon cell culture and tissue samples. This effect sizes and regional heritability. Preprint at bioRxiv https://doi.org/10.1101/318618 research is supported by the National Research Foundation Singapore under its National (2018). Precision Medicine Programme (NPM) Phase II Funding (MOH-000588 to P.T. and W.K.L.) and 67. Wang, G., Sarkar, A., Carbonetto, P. & Stephens, M. A simple new approach to variable administered by the Singapore Ministry of Health’s National Medical Research Council. The selection in regression, with application to genetic fine mapping. J. R. Stat. Soc. B 82, authors thank the Veterinary and Zoology staff at Wildlife Reserves Singapore for their help 1273–1300 (2020). in obtaining the tissue samples, and the Lee Kong Chian Natural History Museum for storage 68. Benner, C. et al. Prospects of fine-mapping trait-associated genomic regions by using and provision of the tissue samples. T.M.B. is supported by funding from the European summary statistics from genome-wide association studies. Am. J. Hum. Genet. 101, Research Council (ERC) under the European Union’s Horizon 2020 research and innovation 539–551 (2017). programme (grant agreement no. 864203), PID2021-126004NB-100 (MICIIN/FEDER, UE) and 69. Stegle, O., Parts, L., Durbin, R. & Winn, J. A Bayesian framework to account for complex Secretaria d’Universitats i Recerca and CERCA Programme del Departament d’Economia i non-genetic factors in gene expression levels greatly increases power in eQTL studies. Coneixement de la Generalitat de Catalunya (GRC 2021 SGR 00177). K.L.-T. is a recipient of PLoS Comput. Biol. 6, e1000770 (2010). Distinguished Professor award from the Swedish Research Council and the Knut and 70. ENCODE Project Consortium. Expanded encyclopaedias of DNA elements in the human Wallenberg Foundation. and mouse genomes. Nature 583, 699–710 (2020). 71. McLaren, W. et al. The ensembl variant effect predictor. Genome Biol. 17, 122 (2016). Author contributions L.F.K.K., J.C.U., S.R., M.A., L.S., G.H., A.J.C., H.G., A. Kumar, F.A., M.J.C., 72. Finucane, H. K. et al. Partitioning heritability by functional annotation using genome-wide H.C., M.H. and K.K.-H.F. performed the analysis and wrote the manuscript. M.C.J., M.K., J.D.O., association summary statistics. Nat. Genet. 47, 1228–1235 (2015). T.B., S. Manu, A.V., J. Bergman, M. Rouselle, F.E.S., L.A., J. Blanc, M.G., D.d.V., I. Goodhead, R.A.H., 73. Coetzee, S. G., Coetzee, G. A. & Hazelett, D. J. motifbreakR: an R/Bioconductor package M. Raveendran, A.J., I.S.C., J.H., C.H., D.J., P.F., J.G.S., F.R.d.M., F.B., H.B., I.S., I.F., J.V., M.M., for predicting variant effects at transcription factor binding sites. Bioinformatics 31, M.N.F.d.S., M.T., R.R., T.H., N.A., C.J.R., A.Z., C.J.J., J.P.-C., G.W., C.A., J.H.S., E.F.-D., S.K., F.S., D.W., 3847–3849 (2015). L.Z., Y.S., G.Z., J.D.K., S.K., M.D.L., E.L., S. Merker, A.N., T.N., C.C.K., J.L., P.T., W.K.L., A.C.K., D.Z., 74. Kulakovskiy, I. V. et al. HOCOMOCO: towards a complete collection of transcription I. Gut, A.D.M., K.G., M.H.S., R.M.D.B., I.K., K.C.W., G.U., C.R. and J.P.B. contributed the primate factor binding models for human and mouse via large-scale ChIP-seq analysis. Nucleic samples, sequencing data and laboratory resources. A.S., A. Kundaje, B.P., K.L.-T., J.R., T.M.B. Acids Res. 46, D252–D259 (2018). and K.K.-H.F. supervised the work. 75. García-Pérez, R. et al. Epigenomic profiling of primate lymphoblastoid cell lines reveals the evolutionary patterns of epigenetic activities in gene regulatory architectures. Nat. Competing interests L.F.K.K., J.C.U., S.R., M.A., L.S., A.J.C., H.G., A.K., F.A., J.G.S. and K.K.-H.F. Commun. 12, 3116 (2021). were employees of Illumina Inc. as of the initial submission of this manuscript. Acknowledgements The authors thank H. Chang, O. Ryder, S. Reilly and E. Karlsson for helpful Additional information discussions. Funding: M.C.J., D.d.V., I. Goodhead, R.M.D.B. and J.P.B. were supported by a UKRI Supplementary information The online version contains supplementary material available at NERC standard grant (NE/T000341/1). H.C. and M.H. were supported by NHGRI U24HG002371. https://doi.org/10.1038/s41586-023-06798-8. M.K. was supported by la Caixa Foundation (ID 100010434), fellowship code LCF/BQ/PR19/ Correspondence and requests for materials should be addressed to Jeffrey Rogers, 11700002, and by the Vienna Science and Technology Fund (WWTF) and the City of Vienna Tomas Marques Bonet or Kyle Kai-How Farh. through project VRG20-001. J.D.O. was supported by la Caixa Foundation (ID 100010434) Peer review information Nature thanks Dimitrios Vitsios and the other, anonymous, reviewer(s) and the European Union’s Horizon 2020 research and innovation programme under the Marie for their contribution to the peer review of this work. Skłodowska-Curie grant agreement no. 847648. The fellowship code is LCF/BQ/PI20/11760004. Reprints and permissions information is available at http://www.nature.com/reprints. Extended Data Fig. 1 | Genome assemblies and constraint metrics. missing data in the filtered MSA, after excluding ambiguous alignments and (a) Distribution of genome assembly span and contiguity for newly assembled duplications for a given species, versus the pairwise phylogenetic distance to primate species in this project. The cluster with assembly spans <2.3 Gb human. The proportion of resolved bases has a strong phylogenetic clustering, corresponds to Strepsirrhines, which have smaller genomes sizes then remaining points are colored by the corresponding primate family following the color primate species. (b) ROC-curves for coding benchmark across mammal and scheme presented in Fig. 1a. (d) Effect of alignment composition on phyloP primate phyloP, comparing codon positions 2 (CD2) as putatively constrained scores for 3 different scenarios: Site 1 contains positions with perfectly positive cases, and human four-fold degenerate sites (4D) as negative cases. matching alignments in 151-171 species and missing alignments in the remaining Both primate and mammal phyloP distinguish well between non-synonymous ones, Site 2 contains positions with perfectly matching alignments in 151-171 CD2 and 4-fold degenerate sites, while mammal phyloP achieves expectedly species but mismatches in over 50 species, Site 3 contains perfect alignments higher performance due to the larger total branch-length covered by the MSA. across all species. Distributions for Site1 and Site 2 are significantly different (c) Scatterplot showing the proportion of bases in the human genome with (P = 1.4 × 10 −66, two-sided Rank Sum Test). Article Extended Data Fig. 2 | Regional and global substitution models. sequences used to calibrate the regional model, resulting in less reliable (a) Comparison of neutral background models with genome-wide random estimates of local substitution rates (<50 kb, annotated as purple crosses). sampling of ancestral repeats from all autosomes (green) versus regional (b) Comparison of performance of global versus regional model at separating modeling of substitution rates at a 1 Mb scale (purple). The upper panel shows codon position 2 (amino acid-altering positions) versus 4-fold degenerate sites median phyloP scores in 1 Mb windows along chromosome 1, the lower panel (synonymous positions), and promoters versus matched distal non-coding the corresponding standard deviations. Median scores and dispersion are very sequence. Global and regional models achieve similar performance on both similar between global and regional neutral models, values of larger discrepancy coding and non-coding benchmarks. tend to fall within windows that containe a limited number of ancestral repeat Extended Data Fig. 3 | Constraint in human protein-coding exons. spliced significantly more often than broadly constrained ones (OR = 1.35, (a) Average per-base mammal and primate phyloP scores for human canonical P = 1.3 × 10 −7, two-sided Fisher’s Exact Test). (d) Mean exon inclusion rates (PSI) protein-coding exons classified by primate-specific constraint. (b) Distribution of alternatively spliced exons across GTEx tissues. Exons constrained specifically of constraint across clades for 185,275 protein-coding exons. Most human in primates have significantly lower inclusion rates than broadly constrained protein coding exons are deeply constrained. (c) Fraction of alternatively ones (P = 8.6 × 10 −6, two-sided Rank Sum Test, n = 28,127 exons). Boxes show spliced exons for exons constrained either specifically in primates, or broadly mean and interquartile range (IQR), whiskers delimit +/− 1.5 x IQR. across mammals. Exons with primate-specific constraint are alternatively Article Extended Data Fig. 4 | Sensitivity analysis of constraint in DHS elements. rates modeled in a regional 1 Mb context and a genome-wide averaged model. (a) Distribution of non-primate mammalian scaling factors for DHS elements (e) Normalized Robinson–Foulds distance between 1 Mb scale phylogeny stratified by clade-specificity of constraint. The dashed gray line denotes and canonical phylogeny along human chromosome 1. (f) Venn diagram where the mammal-constrained and primate-specific constrained distributions intersecting DHS elements on chr1 classified as constrained in primates using intersect. (b) Distribution of primate scaling factors for DHS elements stratified regional substitution rate models and a fixed, canonical topology, or regional by clade-specificity of constraint. (c) Proportion of DHS with primate-specific substitution rate models and a variable, regional topology. Models that constraint for variable FDR cutoffs in mammals excluding primates. Primate accounting for regional differences in topology due to e.g. incomplete lineage FDR is fixed at 5%. (d) Proportion of constrained DHS elements across clades sorting are highly concordant to those that use a single genome-wide topology when modeling substitution rates at a 1 Mb scale, compare to Fig. 2b. The (OR = 806.5, P ≈ 0, two-sided Fisher’s Exact Test). estimated proportions are robust to differences between neutral substitution Extended Data Fig. 5 | UCEs and constrained TF footprints. (a) Overlap between ultraconserved elements as recently defined by Zoonomia (zooUCEs) and primate UCEs allowing up to 1% missing data. (b) Distribution of constraint across clades for TF footprints assessed in this study. Article Extended Data Fig. 6 | Extended characterization of constrained noncoding constraint. n = 3,221 on x-axis and 3,447 on y-axis fine-mapped variants. regulatory elements. (a) Differential chromatin accessibility at orthologous (d) Heritability enrichment as measured by LD Score regression for 6 regulatory sequence elements across 5 primate species. The y-axis indicates the proportion constraint annotations and primate Phastcons. n = 69 traits. Error bars represent of elements where differential accessibility was not detected in (37), stratified 95% CIs. (e) Comparison of noncoding fine-mapped variant enrichment with by sequence constraint. (b) For elements tested by Luciferase reporter in and without adjustment for MAF distributions between the set of variants Fig. 2g, multiple sequence alignments for select primate and mammal species with PIP > 0.5 and the set with PIP < 0.01. Error bars represent 95% CIs, centers are shown for a subsequence of tested elements. Subsequences with high represent point estimates. n = 3,221 fine-mapped variants. (f) Enrichment of DeepLift contribution scores that had matching TF motifs were selected and fine-mapped variants (PIP > 0.5) in DHS elements, related to Fig. 4a,b. Error these data are shown. (c) Comparison between the enrichment of fine-mapped bars represent 95% CIs, centers represent point estimates. Colors indicate variants (PIP > 0.5) in DHS elements or further restricted to TFBSs is shown, sequence constraint, including primate specific constraint as defined by related to Fig. 4a,b. Error bars represent 95% CIs, centers represent point phyloP and by phastCons but not phyloP. n = 3,221 for UKBB and 48,183 for estimates. A grey dashed line indicates y = x. The shape of the point indicates GTEx fine-mapped variants. whether the enrichment is for eQTLs or complex traits. Colors indicate sequence nature portfolio | reporting summary Kyle Farh Jeffrey Rogers Corresponding author(s): Tomas Marques Bonet Last updated by author(s): Oct 9, 2023 Reporting Summary Nature Portfolio wishes to improve the reproducibility of the work that we publish. This form provides structure for consistency and transparency in reporting. For further information on Nature Portfolio policies, see our Editorial Policies and the Editorial Policy Checklist. Statistics For all statistical analyses, confirm that the following items are present in the figure legend, table legend, main text, or Methods section. n/a Confirmed The exact sample size (n) for each experimental group/condition, given as a discrete number and unit of measurement A statement on whether measurements were taken from distinct samples or whether the same sample was measured repeatedly The statistical test(s) used AND whether they are one- or two-sided Only common tests should be described solely by name; describe more complex techniques in the Methods section. A description of all covariates tested A description of any assumptions or corrections, such as tests of normality and adjustment for multiple comparisons A full description of the statistical parameters including central tendency (e.g. means) or other basic estimates (e.g. regression coefficient) AND variation (e.g. standard deviation) or associated estimates of uncertainty (e.g. confidence intervals) For null hypothesis testing, the test statistic (e.g. F, t, r) with confidence intervals, effect sizes, degrees of freedom and P value noted Give P values as exact values whenever suitable. For Bayesian analysis, information on the choice of priors and Markov chain Monte Carlo settings For hierarchical and complex designs, identification of the appropriate level for tests and full reporting of outcomes Estimates of effect sizes (e.g. Cohen's d, Pearson's r), indicating how they were calculated Our web collection on statistics for biologists contains articles on many of the points above. Software and code Policy information about availability of computer code Data collection No software was used for data collection. Data analysis Megahit (v1.2.9) RepeatMasker (v4.1.2-p1) Cactus (v2.1.1) minimap2 (v2.17-r941) PHAST (v1.5) FINEMAP (v1.4) Su5iE (v0.12.10) LDstore (v2) motifbreakR (v3.16) Panther (v16) IQTree (v1.6.12) VEP (v85) GenomicRanges (v3.1.7) April 2023 Haltoos (v2.2) For manuscripts utilizing custom algorithms or software that are central to the research but not yet described in published literature, software must be made available to editors and reviewers. We strongly encourage code deposition in a community repository (e.g. GitHub). See the Nature Portfolio guidelines for submitting code & software for further information. 1 nature portfolio | reporting summary Data Policy information about availability of data All manuscripts must include a data availability statement. This statement should provide the following information, where applicable: - Accession codes, unique identifiers, or web links for publicly available datasets - A description of any restrictions on data availability - For clinical datasets or third party data, please ensure that the statement adheres to our policy Primate assemblies have been deposited at ENA under the accession PRJEB67744. The MSA and constraint tracks are available through the UCSC genome browser. Research involving human participants, their data, or biological material Policy information about studies with human participants or human data. See also policy information about sex, gender (identity/presentation), and sexual orientation and race, ethnicity and racism. Reporting on sex and gender not applicable Reporting on race, ethnicity, or not applicable other socially relevant groupings Population characteristics not applicable Recruitment not applicable Ethics oversight not applicable Note that full information on the approval of the study protocol must also be provided in the manuscript. Field-specific reporting Please select the one below that is the best fit for your research. If you are not sure, read the appropriate sections before making your selection. Life sciences Behavioural & social sciences Ecological, evolutionary & environmental sciences For a reference copy of the document with all sections, see nature.com/documents/nr-reporting-summary-flat.pdf Life sciences study design All studies must disclose on these points even when the disclosure is negative. Sample size Sample size was based on available primate genomes. Data exclusions No data was excluded from analyses. Replication All wet-lab experiments were replicated 3 or more times. All attempts at replication were successful. Randomization Randomization was not necessary for this study, as genomic sequences were derived from a single individual. Blinding Investigators were not blinded to the species names or genome sequences. Blinding was was not necessary for this study, as genomic sequences were derived from a single individual. Reporting for specific materials, systems and methods We require information from authors about some types of materials, experimental systems and methods used in many studies. Here, indicate whether each material, system or method listed is relevant to your study. If you are not sure if a list item applies to your research, read the appropriate section before selecting a response. April 2023 2 Materials & experimental systems Methods nature portfolio | reporting summary n/a Involved in the study n/a Involved in the study Antibodies ChIP-seq Eukaryotic cell lines Flow cytometry Palaeontology and archaeology MRI-based neuroimaging Animals and other organisms Clinical data Dual use research of concern Plants Eukaryotic cell lines Policy information about cell lines and Sex and Gender in Research Cell line source(s) Human iPSCs were obtained from the Stanford CVI iPSC Biobank. Authentication Immunofluresence assay was performed to check the expression of stem cell markers such as NANOG, POUSF1 and SOX2. SNP karyotyping was tested through HuCytoSNP-12 chip (Illumina), and CNV and SNP visualization was performed using KaryoStudio v1.4 (Illumina). Mycoplasma contamination We confirmed that the cell line was negative for mycoplasma contamination using MycoAlertTM PLUS Mycoplasma Detection Kit (Lanza) Commonly misidentified lines No commonly misidentified cell lines were used. (See ICLAC register) Plants Seed stocks No plant materials were used in this study. Novel plant genotypes No plant materials were used in this study. Authentication No plant materials were used in this study. April 2023 3