Influence of cis-regulatory elements on expression divergence in human segmental duplications

Human-duplicated genes are differentially expressed

Our previous work characterizing a subset of HSD genes identified divergent expression patterns over a short evolutionary time span (Shew et al. 2021). Leveraging recent work identifying additional duplicated gene families with at least one paralog unique to humans (Soto et al. 2025), we sought to expand our understanding of gene expression and regulation in HSDs. We also generated long-read RNA-seq data (Pacific Biosciences [PacBio] Kinnex) in GM12878 and SH-SY5Y, two cell lines amenable to functional genomic assays. LCLs like GM12878 are accessible and directly comparable between primate species, and SH-SY5Y retains a neuronal progenitor-like state without in vitro reprogramming (Pezzini et al. 2017). Considering 605 duplicate gene families containing 1394 genes, long- and short-read transcript abundance (Supplemental Table S1; The ENCODE Project Consortium 2012; Pezzini et al. 2017; Luo et al. 2020; Hitz et al. 2023) is strongly correlated in both cell types (r = 0.85 for GM12878 and 0.88 for SH-SY5Y) (Supplemental Fig. S1). Although the long-read quantification is expected to be more accurate for distinguishing paralogs from each other, it is less sensitive overall, with a larger proportion of genes not detected in either cell type (17.5% with long reads vs. 9.0% with short reads) (Supplemental Table S2). Across duplicated gene families, we saw an enrichment of highest-expressed paralogs in regions with synteny to the chimpanzee ortholog (P = 2.7 × 10⁻⁴ for GM12878 Kinnex; P = 1.7 × 10⁻³ for GM12878 short-read; P = 2.8 × 10⁻⁶ for SH-SY5Y Kinnex; P = 2.6 × 10⁻³ for SH-SY5Y short-read) (Fig. 1). This is in agreement with our previous work showing ancestral paralogs exhibit highest expression, as opposed to human-specific, derived genes (Shew et al. 2021). To the same end, we also found an enrichment of coding genes as opposed to pseudogenes among highest-expressed paralogs (P = 2.7 × 10⁻⁸ for GM12878 Kinnex; P = 1.5 × 10⁻³ for GM12878 short-read; P = 2.0 × 10⁻¹⁰ for SH-SY5Y Kinnex; P = 7.7 × 10⁻³ for SH-SY5Y short-read) (Fig. 1). This same enrichment was observed in RNA-seq data from primary adult liver, kidney, heart, and lung tissue; additional LCLs; and fetal neocortex (Supplemental Fig. S2; Pickrell et al. 2010; Fietz et al. 2012; Florio et al. 2015; Blake et al. 2020). Taken together, these analyses confirm that human-specific genes exhibit divergent expression patterns, and the ancestral paralogs are most likely to retain a conserved function.

Figure 1.

Human-duplicated genes are differentially expressed. The highest-expressed paralogs in each gene family (top point per column and inset violin plots) are more likely to be syntenic with chimpanzee (blue) and more likely to be protein coding (circle). Comparisons were performed with transcriptomic data from GM12878 lymphoblastoid (top) and SH-SY5Y neuroblastoma cells (bottom) sequenced with short-read Illumina (left) and long-read PacBio Kinnex (right) reads. Gene families are represented on the x-axis, and transcript abundance on the y-axis (transcripts per million [TPM], TS [total transcript]). Gene family expression is ordered from highest (left) to lowest (right) with high-confidence human-specific duplications, which are the focus of the remainder of the study, numbered and demarcated with a vertical line. The corresponding number followed by the name of the highest-expressed paralog of these priority gene families is shown below the plot and colored to indicate syntenic status: likely ancestral (blue) or nonsyntenic, likely derived (red).

Open in new tabDownload PowerPointLink to figure

MPRA identifies active regulatory elements in human duplications

To systematically quantify cis-regulatory activity in HSDs, we designed an MPRA to directly compare paralogous human and orthologous chimpanzee sequences from 67 genes in 26 gene families (Fig. 2A; Supplemental Tables S3, S4). These gene families were chosen based on their correct representation in the GRCh38 assembly (Nurk et al. 2022; Liao et al. 2023), known elevated copy numbers in all humans surveyed compared with nonhuman primates (Soto et al. 2025), and previously established ancestral and derived identities (Dennis et al. 2017). Additionally, our expression analysis shows 72% (18/25) of these gene families detected with short reads (>1 transcripts per million [TPM]) in GM12878 or SH-SY5Y, and 68% (17/25) with long reads (>1 counts per million [CPM]) (Supplemental Table S5). We used chromatin accessibility and histone H3, lysine 27 acetylation (H3K27ac) chromatin immunoprecipitation sequencing (ChIP-seq) data to identify candidate CREs from the fetal and adult prefrontal cortex (Reilly et al. 2015; Vermunt et al. 2016; de la Torre-Ubieta et al. 2018), owing to the known role of HSD genes in brain development, and LCLs (Degner et al. 2012; McVicker et al. 2013; Bryois et al. 2018), whose accessibility has made them foundational for cross-species comparisons in primates (Supplemental Table S1). HSDs are typically excluded by standard bioinformatic approaches, as reads aligning to these regions have poor mapping quality. To circumvent this, we permissively identified candidate CREs by remapping epigenomic data to the human reference (GRCh38) with multiple alignments. A total of 1163 CREs in ancestral loci were divided into 200-mer “tiles” with 100 bp overlap (N = 2675), and homologous tiles in the human and chimpanzee (panTro6) genomes were identified by sequence alignment. In total, we synthesized 8145 test sequences and 500 scrambled negative controls (Supplemental Fig. S3; Supplemental Table S6), which were packaged into lentiviral vectors and assayed in GM12878 and SH-SY5Y cells according to the lentiMPRA protocol (Gordon et al. 2020). Although these are cancer-derived cell lines with abnormal karyotypes, they provide a relevant trans-regulatory environment, and CRE activity is measured at many genomic positions. Transduction efficiencies were comparable to previous lentiMPRA studies, with a multiplicity of infection (MOI) of about 40 for SH-SY5Y and about 13 for GM12878. Technical replicates (N = 3 per cell line) exhibited high reproducibility (mean pairwise r = 0.93 and 0.93 for within-cell-type comparisons of DNA and RNA libraries, respectively) (Supplemental Fig. S4; Nurk et al. 2022; Liao et al. 2023).

Figure 2.

MPRA quantifies paralog-specific regulatory activity of duplicated CREs. (A) MPRA tiles were designed to multimapped chromatin accessibility and histone H3K27ac ChIP-seq data previously published from human primary brain (fetal and adult) and lymphoblastoid cell lines. Peak regions were tiled with 200-mers at 2× density in ancestral loci, and homologous positions were identified by alignment to the human and chimpanzee genomes. All unique sequences were synthesized in an oligonucleotide pool and cloned into a lentiviral vector. MPRA was performed in the neuroblastoma cell line SH-SY5Y and lymphoblastoid cell line GM12878. (B) Pie charts depict the number and proportion of active sequences (left) and sequence families containing at least one active sequence (right) for each cell type. (C) Sequence activity (estimated transcription rate, ɑ) is plotted type for 7750 tiles measured in both assays. Sequences active over negative controls are colored red (SH-SY5Y), blue (GM12878), and purple (both). The Venn diagram indicates the number of sequences in each category.

Open in new tabDownload PowerPointLink to figure

We first quantified enhancer activity for all tested CREs in both cell lines using the estimated transcription rate (ɑ), modeled from the integrated (DNA) and expressed (RNA) barcode counts (MPRAnalyze) (Ashuach et al. 2019). In SH-SY5Y, 1056/7760 (14%) tested sequences were active relative to the negative controls (median absolute deviation P < 0.05), whereas in GM12878, 1850/7750 (24%) sequences were active (Fig. 2C; Supplemental Tables S7, S8). Of these, 472 (19%) of 2434 active sequences were active in both cell types (Fig. 2B; Supplemental Note S1; Supplemental Fig. S6). Notably, higher activity was measured overall in SH-SY5Y (median ɑ = 2.50, vs. 0.73 in GM12878) (Fig. 2C), likely owing to the threefold greater infection rate achieved in this cell line. Activity measurements were correlated for overlapping tiles (r = 0.21 for GM12878 and 0.43 for SH-SY5Y) (Supplemental Fig. S5C), but we opted against combining measurements across tiles because any target nucleotide was covered by at most two constructs, and the 100 bp offset is much larger than transcription factor binding sites.

Sequences within 1 kb of transcription start sites (TSSs) were enriched for activity in both cell types (1.35-fold higher proportion of active sequences for promoters in GM12878; 3.46-fold in SH-SY5Y; Fisher's exact test P < 0.02). Although we observed a correlation between promoter activity (mean activity score for all active tiles) and mRNA expression driven by the contrast between active and inactive promoters (Supplemental Fig. S7A), MPRA promoter measurements did not predict differences in ancestral/derived gene pairs (Supplemental Fig. S7B). This suggests a role for distal elements in shaping expression patterns.

Species- and paralog-specific CRE activity

We first examined sequence activity across species. Considering all human–chimpanzee sequence pairs with at least one active element, 35% (357/1017) were differentially active in SH-SY5Y, and 11% (220/1935) were differentially active in GM12878 (Supplemental Fig. S8; Supplemental Tables S9, S10). We found an expected relationship between sequence divergence and the magnitude of differential activity between species in both cell types (r = −0.07, P < 1 × 10⁻³) (Supplemental Fig. S9). However, the strongest differentially active sequences tended to have the highest similarity with chimpanzee, demonstrating that a small number of substitutions can alter sequence activity (Supplemental Fig. S9). Based on our findings that ancestral genes generally exhibit higher expression and increased concordance with chimpanzee orthologs (Fig. 1; Shew et al. 2021), we next hypothesized that activity losses preferentially occur in derived CREs. We partitioned human sequences by ancestral or derived status but found no enrichment of activity gains or losses (considering chimpanzee a proxy for ancestral activity) in either group (hypergeometric test). However, the measured activities of ancestral sequences correlated significantly better than derived sequences with their chimpanzee orthologs (Fisher's z-test, P < 0.01 in both cell types) (Fig. 3A). Thus, although cis-regulatory activity in HSDs is largely conserved between humans and chimpanzees, human-specific derived regions are more diverged, potentially experiencing relaxed evolutionary constraint on CRE activity.

Figure 3.

Differential activity between homologous sequences. (A) Scatterplots depict the activity of human and orthologous chimpanzee sequences, partitioned by human evolutionary status (ancestral above, derived below). Each sequence pair is colored by differential activity relative to chimpanzee: green (human lower), purple (human higher), gray (not differential), black (neither sequence active). The dotted line marks the identity y = x, and Pearson's correlation coefficient is printed on each plot. (B) Volcano plots depict differential activity between paralogous human sequences, with each unique derived–ancestral pair plotted separately. Points are colored by differential activity relative to the ancestral human paralog: blue (derived lower), red (derived higher), gray (not differential), and black (neither sequence active). The marginal histograms depict the distribution of log-fold differences for differentially active sequences.

Open in new tabDownload PowerPointLink to figure

Comparing activities between human duplications, we found a similar proportion of derived–ancestral pairs were differentially active by the same definition: 216/640 (34%) of pairs in SH-SY5Y, with a median 1.13-fold difference; in GM12878, 147/1206 (12%) were differentially active, with a median 1.20-fold difference (Fig. 3B,C). Although most paralogs exhibited modest divergence in activity (1.13-fold and 1.20-fold differences in SH-SY5Y and GM12878, respectively), SRGAP2 and FRMPD2 contained strong differential human-specific sequences (more than twofold) (Supplemental Tables S11, S12) that were also concordant in both cell types, suggesting possible single CRE drivers of paralog differences. Targeting the putative CRE from SRGAP2 (Chr 1: 206,333,093–206,333,293; GRCh38), we used a luciferase reporter assay to validate that a larger homologous region containing the human-specific SRGAP2B (0.43-fold) and SRGAPC (0.36-fold) tile sequences was significantly less active than the ancestral SRGAP2 in SH-SY5Y, with no differential activity observed in GM12878, matching MPRA results (Supplemental Fig. S10; Supplemental Table S13).

Taking this same approach to test larger regions with luciferase reporters, we used 1 kb windows in HSDs to score total measured MPRA activity, secondarily ranking by the maximum fold difference of any contained tile and manually curating select differential elements in both cell types (Supplemental Fig. S11). In the CHRNA7 locus, we tested two candidate regions. The first (Chr 15: 32,183,316–32,184,316; GRCh38) was predicted to have higher activity in the derived CHRFAM7A, which was true for GM12878 (more than 1.3-fold difference) (Fig. 4; Supplemental Fig. S12; Supplemental Table S13). The second was expected to have lower activity in SH-SY5Y, which was also verified by luciferase (less than 0.7-fold difference). In addition, an intronic region within ancestral HYDIN was predicted to have higher activity in the derived HYDIN2 in both cell types. Reporter activity was significantly higher for both homologs in GM12878 and for HYDIN2 in SH-SY5Y (more than 3.8-fold difference) (Supplemental Figs. S12, S13; Supplemental Table S13). Altogether, these experiments demonstrate paralog-specific activity for enhancers in a more biologically relevant sequence context. The reporter activities from these constructs comprising expanded putative CREs (1 kb) were generally in agreement with the MPRA, which measured activity of 200-mer tiles in isolation.

Figure 4.

Paralog-specific regulation at the CHRNA7 locus (Chr 15q13.3). A 143 kb region (Chr 15: 32,153,206–32,460,660; GRCh38) is shown. Segmental duplications (SDs) are colored orange for >99% sequence identity and gray for <98% identity. No ENCODE candidate cis-regulatory elements (cCREs) are annotated in this locus. The MPRA data are visualized as activity relative to the human ancestral locus and plotted on these coordinates. Each point represents a 200-mer tile, and tiles are colored by which homolog they correspond to (teal: chimpanzee CHRNA7; orange: human-specific CHRFAM7A; gray: identical in both) and filled in if they scored as differentially active. Two 1 kb regions are highlighted in the insets; the left panel shows a zoomed in plot of the MPRA measurements per tile, and the right shows a relative luciferase activity of the entire sequence, as assayed in the corresponding cell type. Asterisks on the box plots mark significantly different activity with respect to the human ancestral sequence. See also Supplemental Table S13.

Open in new tabDownload PowerPointLink to figure

Collective analysis of CREs on paralog expression

Because of modest effect sizes of individually tested CREs in the MPRA, we next considered whether their collective activity might explain paralog-expression divergence. Building on the activity-by-contact (ABC) score (Fulco et al. 2019), we developed a similar metric using activity levels directly measured from the MPRA (Methods) (Supplemental Fig. S14A). For expressed derived–ancestral gene pairs (either TPM > 1), expression divergence in the short-read data was significantly correlated with the summed ABC score ratio in SH-SY5Y (P < 0.05, r = 0.70; Supplemental Fig. S14B). This comparison was robust to various cutoffs for considering CREs active or differential and showed the best correlation for more stringent criteria. No relationship was observed using the long-read expression data, however, perhaps reflecting reduced sensitivity of Kinnex to detect transcripts (Supplemental Fig. S14C; Supplemental Note S2; Supplemental Fig. S15). The SH-SY5Y MPRA ABC score ratio also correlated with paralog expression divergence quantified in certain primary brain cell types (i.e., the cortical plate, cortical neurons without basal contact, and the outer subventricular zone) (Supplemental Fig. S16), reflecting functional relevance beyond the cell lines assayed.

We searched for trends in CRE activity that could explain this, namely, in strength, number, or location. In both cell types, neither measured CRE activity nor the number of CREs differed between ancestral and derived regions, but active CREs were slightly closer proximally to ancestral TSSs than to derived ones (P = 0.03, Wilcoxon rank-sum test) (Supplemental Fig. S17). This suggests a trend toward loss of active elements closer to derived paralogs. Because promoter strength alone did not predict expression differences (Supplemental Fig. S7), considering the full complement of CREs may better explain differential mRNA levels across HSD genes.

Indeed, the fraction of contribution of proximal CREs (within 1 kb) as determined by the ABC metric was generally concordant across paralogs (Supplemental Fig. S18A). There were a few exceptions, namely, ARHGEF34P and ARHGAP11B in SH-SY5Y and PTPN20CP and HYDIN2 in GM12878 (Supplemental Fig. S18B,C). Across all elements, the median fraction of contribution of the top CRE to each HSD gene was only 5% (Supplemental Fig. S18D). Together, these results indicate that changes to promoter activity are not the primary driver of HSD expression divergence in most gene families and that many distal elements may influence this process.

For comparison, we also calculated the predicted difference using the ABC score as published, with chromatin accessibility and H3K27ac ChIP-seq. We aligned data sets from GM12878 (Zhang et al. 2020) and SH-SY5Y (Zimmerman et al. 2018; Gartlgruber et al. 2021) to the human reference (GRCh38), allocating multiple alignments probabilistically (Chung et al. 2011; Zhang et al. 2020), enabling the calculation of the ABC difference metric across HSDs and adjacent, unique regions (Methods). For SH-SY5Y, we found a correlation of the ABC score ratio with short-read expression divergence (P = 0.05, r = 0.42), but this did not hold when subsetting to only the duplicated or unique space (Fig. 5A–C). The long-read expression data for both cell types also trended in the expected direction for all and nonduplicated peaks. We also scored the exact regions tested by our MPRA, in place of peaks derived from the epigenomic data, and found a weaker correlation than the MPRA analysis on its own (P = 0.06, r = 0.39) (Fig. 5D). Thus, the measured regulatory activity from the MPRA appears more informative than proxies based on chromatin state. The whole-genome chromatin activity data, however, suggest that more distal, nonduplicated CREs have a greater impact on expression divergence.

Figure 5.

Comparison of SH-SY5Y expression divergence to predictions from epigenomic data. For each derived–ancestral (der./anc.) gene pair, expression divergence was correlated to the ratio of summed ABC scores calculated from multimapped ATAC-seq and H3K27ac ChIP-seq. The comparison was performed on the following sets of regions: all CREs within 5 Mb, as defined by our ChIP-seq analysis (A); CREs within HSDs (B); CREs outside of HSDs (C); or the exact positions tested by MPRA (D). Pearson's correlation and significance of the linear relationship is printed on each plot, and the 95% confidence interval for the regression is shaded in gray.

Open in new tabDownload PowerPointLink to figure

Finally, recent advances in long-read epigenomic assays have better enabled differentiation of regulatory information between human duplications (Altemose et al. 2022; Gershman et al. 2022). Using published directed methylation with long-read sequencing (DiMeLo-seq) data sets from GM12878 (Maslan et al. 2024), we identified expected relationships between DNA methylation (5mC), histone modifications (H3K4me3, H3K27ac), and HSD gene expression, including derived–ancestral differences (Supplemental Fig. S19). These data corroborated our findings using short-read ChIP, identifying a similar correlation between regulatory activity (ABC score) and expression divergence that was primarily driven by nonduplicated elements (Supplemental Note S3; Supplemental Fig. S20). Moving forward, we expect long-read epigenomic data sets will be critical for disentangling the mechanisms underlying paralogous expression divergence in other loci and cell types.

RNA library preparation

Total RNA was extracted from GM12878 and SH-SY5Y cell lines using the Qiagen RNeasy plus kit following the manufacturer's protocol. RNA integrity was confirmed using the Agilent 2100 Bioanalyzer, ensuring RIN values ≥7.0. Full-length (FL) RNA libraries were prepared using the PacBio Kinnex FL RNA kit (103-238-700), including cDNA synthesis, Kinnex PCR amplification, and library cleanup. Libraries were quantified with the Qubit DNA HS assay kit and size-checked on the Agilent Femto Pulse system, which showed average insert sizes of ∼15 kb for both GM12878 and SH-SY5Y, before sequencing on the PacBio Revio platform.

RNA-seq analysis

RNA-seq data were obtained for GM12878 (ENCODE ENCSR000AEC, ENCSR000AEE, and ENCSR000CVT) and SH-SY5Y (Pezzini et al. 2017). Transcripts were quantified with Salmon v1.9.0 (Patro et al. 2017) with the flags “‐‐validateMappings ‐‐gcBias,” the telomere-to-telomere CHM13 v2.0 CAT/Liftoff transcriptome, and the CHM13 v2.0 assembly as decoy sequence. All identical transcripts were removed from the transcriptome prior to index construction. TPM values were summed to the gene level using tximport (Soneson et al. 2015). For relative expression analyses, derived–ancestral gene pairs were considered if at least one paralog was expressed >1 TPM. Expression divergence was calculated as the log₂ ratio of derived/ancestral expression, with a pseudocount one order of magnitude smaller than the smallest nonzero value added to each. For direct comparison of count-normalized data sets, raw reads were downsampled with seqtk v 1.4.

Long-read PacBio Kinnex analysis

Initial processing of PacBio Kinnex data was done with the Iso-Seq workflow (https://isoseq.how/). HiFi reads were first segmented using PacBio concatenated read splitter Skera at adapter positions to generate segmented reads. Lima is then used to identify and remove barcoded Iso-Seq primers. The output of lima are FL reads. These FL reads are then processed with Iso-Seq refine to identify and trim poly(A) tails, as well as identify and remove concatemers. The output of Iso-Seq refinement is full-length, nonconcatemer (FLNC) reads. These FLNC reads are then quantified using Oarfish (Zare Jousheghani et al. 2025) using the same identical transcript-filtered CHM13 v2.0 CAT/Liftoff transcriptome described above. The estimated transcript counts are then summed to get gene-level counts, which were converted to CPM for expression analyses. No length normalization was performed because each transcript detected is expected to yield one FLNC read.

MPRA oligo design

Candidate CREs were identified from H3K27ac ChIP-seq, ATAC-seq, or DNase-seq data in LCLs, fetal cortex, and adult prefrontal cortex with previously published data (Degner et al. 2012; McVicker et al. 2013; Reilly et al. 2015; Vermunt et al. 2016; Bryois et al. 2018; de la Torre-Ubieta et al. 2018). In the case of Bryois et al. (2018), only control samples were used. Reads were trimmed with Trimmomatic (Bolger et al. 2014) aligned to GRCh38, allowing for multiple mapping with Bowtie v1.1.2 (Langmead et al. 2009) with the flags “-a -v2 -m 99.” Peaks were called using MACS2 v2.1.2 (Liu 2014) with a 5% FDR and the following additional settings for chromatin accessibility data: “‐‐nomodel ‐‐shift -100 ‐‐extsize 200 ‐‐broad.” Peaks were called on all samples in each study and combined by reporting genomic regions with nonzero coverage in a minimum number of samples using BEDTools genomecov (Quinlan and Hall 2010): 3/10 LCL H3K27ac, 52/204 LCL DNase, 2/3 fetal cortex H3K27ac, 3/3 fetal ATAC-seq (cortical plate or germinal zone), 1/1 adult prefrontal cortex H3K27ac (peaks called jointly from samples HS1 and HS2), and 102/137 adult prefrontal cortex ATAC-seq. Cutoffs were chosen based on manual inspection of reproducibility and to represent sequences from LCLs and brain approximately equally.

To design the test oligos, all peaks within 400 bp were merged, the maximum distance at which peaks could become adjacent after coordinate expansion. The resulting intervals were padded to the nearest 100 bp, plus an additional 100 bp on either side to allow full 2× coverage of the initial target regions using 200-mer tiles. Tiles were generated with BEDTools make windows. To match orientation of promoters, tiles were assigned to the strand of the nearest annotated feature (GENCODE v32) or randomly selected if multiple features were nearest. This resulted in random strand assignments to 140/661 promoter and 501/7099 distal tiles. Tiles were intersected with ancestral HSD regions and aligned with BLAT (Kent 2002) to a custom reference consisting of GRCh38 with non-HSD sequence masked, in addition to the missing portions of the DUSP22B and GPRIN2B contigs (ancestral loci and contigs from Dennis et al. 2017). Alignments were then filtered in a two-step process to keep only paralogous positions found in human-specific SDs: (1) keep if there were four or fewer hits >190 bp with >95% identity; (2) else keep if there were four or fewer hits >195 bp with >98% identity. This logic maximized target regions by allowing greater divergence for individual tiles on the first pass but increasing stringency if the number of hits was still greater than the expected number of paralogs. To focus on paralogous differences, sequences for which all alignments were identical were discarded. To identify chimpanzee orthologs, human tiles were mapped to the panTro6 assembly using liftOver (UCSC Genome Browser utilities) and filtered for being with 5% (10 bp) of the original 200 bp size (93% success rate). All tiles with this limit were trimmed or expanded to 200 bp; 91% of alignments were within 1 bp difference. FASTA sequences were extracted in a strand-aware manner from GRCh38 or panTro6 and deduplicated. Finally, a universal priming sequence (5′-AGGACCGGATCAACT-[200-mer]-CATTGCGTGAACCGA-3′) to the pLS-SceI vector (Addgene 137725) was added to all 200-mers. Sequences were filtered to remove AgeI and SceI restriction sites, as well as those made singletons after filtering, leaving 8145 test oligos. In addition, 500 ancestral tiles were sampled and scrambled to generate negative controls. All 230mer sequences were synthesized by Agilent Technologies.

MPRA library preparation and sequencing

MPRA was carried out according to the lentiMPRA protocol as previously described (Gordon et al. 2020) with three technical replicates. We generated the lentivirus library by transfecting HEK293T cells with the plasmid library. To titrate the lentivirus, GM12878 and SH-SY5Y were plated into a 24-well plate and infected with varying volumes of the virus in each well. In these experiments, transduction efficiency was ∼10% for GM12878 and ∼30%–40% for SH-SY5Y cells, as determined by GFP-positive cell counts on the Countess II FL. Following a 3 day incubation, genomic DNA was extracted from each well. Then, the MOI was calculated for each viral concentration via qPCR using plasmid backbone and genomic DNA control primers, as described in the protocol. From these titration results, the final infection conditions were set as follows: For SH-SY5Y experiments, 2.5 million cells were infected at a MOI of 40 in DMEM/F12 GlutaMAX medium (Thermo Fisher Scientific 10565018) with 10% FBS (Thermo Fisher Scientific 10-438-026) and 2.5 µg/mL protamine sulfate (Spectrum Chemical TCI-P0675-1G). For GM12878, 55 million cells were transduced at an MOI of 13 in RPMI (Fisher Scientific 11875135) with 10% FBS and 8 µg/mL polybrene (Millipore Sigma MILL-TR-1003-G). Cells were incubated at 37°C in 5% CO₂, and nucleic acid extractions were conducted 48 h after infection. Barcode association was performed on a PE150 NextSeq mid-output run. Barcode counting from DNA and RNA libraries was performed on three PE15 NextSeq high-output runs.

MPRA data analysis

Barcode-insert association was performed with MPRAflow (Gordon et al. 2020) using “‐‐mapq 1” and all other parameters set to default. We manually confirmed that promiscuous barcodes had a majority assignment to one insert and were not mixed between paralogs. Barcodes were counted with the count utility and passed to MPRAnalyze (Ashuach et al. 2019). Counts were depth-normalized with the total sum, and active sequences were defined per cell type against the negative controls with the following model: “dnaDesign = ∼ barcode + replicate, rnaDesign = ∼1.” To identify differential activity of homologous sequences, counts matrices were reformatted to combine all homologs in the same row for model fitting to use homolog identity used as a covariate. Models were fit in “scale” mode with the following design: “rnaDesign = ∼ homolog, reducedDesign = ∼ 1.” The alpha-value used as a measure of sequence activity refers to the estimated transcription rate from the model. Differentially active homologs were defined in two ways: (1) against the chimpanzee ortholog to define relative gains/losses in humans and (2) against the human ancestral sequence to quantify divergence of ancestral-derived pairs. Differential activity was assessed in a pairwise manner with the testCoefficient() function (Wald test), and differential pairs were defined at a 5% FDR using the Benjamini–Hochberg procedure.

ABC score

Regulatory activity for each TSS was calculated using the ABC score framework defined in work by Fulco et al. (2019), substituting the estimated transcription rate ɑ from MPRAnalyze for activity and using (genomic distance)^−0.7 as a proxy for contact frequency. CREs within 5 Mb of each TSS were considered for the calculation. Although ABC values per CRE are typically normalized and interpreted as a relative contribution to a given TSS, we compared summed ABC scores between paralogous genes. Regulatory divergence was scored in a similar manner to expression divergence: A pseudocount of one order of magnitude below the smallest nonzero value was added to all ABC values, and the log₂ ratio of scores was calculated for each derived–ancestral gene pair.

ChIP-seq and ATAC-seq analyses

For calculation of ABC scores as published (using H3K27ac ChIP and ATAC data) and comparison with the MPRA, we aligned short reads (Supplemental Table S1) to GRCh38 with Bowtie v1.1.2 (Langmead et al. 2009), allowing for multiple mapping with the flags “-a -v2 -m 99.” CSEM was used to assign a probability to each alignment position on the basis of the local unique mapping rate (Chung et al. 2011; Zhang et al. 2020). Reads were then allocated to their most likely position by selecting the alignment with the highest posterior probability or at random if all mappings were equally likely. Peaks were called on the ChIP data using input controls at a 5% FDR using MACS2 (Liu 2014). For ATAC data, the additional parameters “‐‐nomodel ‐‐shift -100 ‐‐extsize 200 –broad” were added. Read counts per peak were determined with BEDTools coverage (Quinlan and Hall 2010).

Luciferase reporter assays

One kilobase fragments were synthesized by Azenta Life Sciences and cloned into the pE1B vector using the NEBuilder HiFI assembly master mix (NEB E2621). GM12878 and SH-SY5Y cells were cultured as described above. SH-SY5Y cells at 70%–90% confluence were cotransfected with 50 ng pRL-TK and equimolar test construct using Lipofectamine 3000 (Thermo Fisher Scientific L3000001). Cells were lysed 48 h posttransfection and assayed with the dual-luciferase reporter assay system (Promega E1910) and the Tecan Spark plate reader according to the manufacturer's instructions. GM12878 cells were split 48 and 24 h prior to transfection and electroporated with 12.5 µg pRL-TK and equimolar test construct, brought to 100 µL with RPMI. For each transfection, 12.5 million cells were electroporated with the Neon system (Thermo Fisher Scientific) using buffer E2 and the following program: 1200 V, 20 msec, three pulses. Cells were recovered in 4.3 mL of prewarmed media (RPMI + 15% FBS, no antibiotics). GFP controls were used to monitor transfection efficiency (∼15%). 24 h posttransfection, cells were washed in PBS, plated at 500,000 cells per well, and lysed in an optical plate. Luminescence measurements were performed as above. The activity of each construct was quantified relative to the empty vector, and differential activity was relative to the human ancestral sequence. Significant differences were determined by ANOVA, using batch as a covariate, followed by a Tukey post hoc test.

DiMeLo-seq analyses

DiMeLo data for H3K4me3 and H3K27ac in GM12878 were obtained from GSE208125 (Maslan et al. 2024) as BAM files aligned to T2T CHM13 v2.0. Base modifications (5mC and 6mA) were thresholded at 155 (corresponding to probability of ∼60%), and the frequency of base modification per site was calculated using modkit pileup (Nanoporetech). For comparison with short-read ChIP data and MPRA tiles, coordinates were mapped to T2T CHM13 v2.0 using liftOver (UCSC Genome Browser Utilities). Then, the mean fraction of modified bases was calculated for each interval in all peak sets of interest using BEDTools map. The summed ABC score was calculated per HSD gene as above, using only the H3K27ac signal. Finally, differences in signal (mean fraction modified bases) for all assays (5mC, H3K27ac 6mA, and H3K4me3 6mA) were compared between chimpanzee-syntenic and -nonsyntenic regions, also converted to CHM13 v2.0 coordinates.