Simultaneous epigenomic profiling and regulatory activity measurement using e2MPRA

Domain: biology
Published: 2026-01-14
Editorial outcome: minor_revision

Abstract
Using various biochemical assays that identify transcription factor (TF) binding and histone modifications, cis-regulatory elements (CREs) can be annotated in a genome-wide manner. However, these assays are descriptive and require functional validation. To the best of our knowledge, no technology can simultaneously analyze the regulatory function and epigenomic modifications of a specific sequence. Here, we develop an enrichment followed by epigenomic profiling massively parallel reporter assay (e2MPRA). This technique uses lentivirus to enrich for the integration of specific CREs into the genome and applies MPRA, Cut&Tag or ATAC-seq on them enabling simultaneous, high-throughput analysis of regulatory activity, protein binding, and epigenetic modification. We demonstrate that e2MPRA can dissect the epigenetic functions of TF motifs arranged within synthetic enhancers and evaluate the effects of sequence perturbation on epigenetic states. In summary, e2MPRA advances our understanding of the regulatory code, its effect on the epigenome and how its alteration leads to phenotypic effects.

Introduction
Gene expression is regulated in a spatiotemporal manner by active cis-regulatory elements (CREs), such as promoters, enhancers, insulators and silencers. Nucleotide variants in CREs can alter gene expression, which can lead to subsequent changes in phenotypes, including cell differentiation, disease, evolution, and other biological phenomena. Active regulatory elements such as promoters and enhancers are bound by various transcription factors (TFs) that play different roles in transcriptional regulation, such as pioneering activity that induces chromatin accessibility, transcriptional activation and repression, histone modification, and chromatin remodeling activities that mediate higher-order chromatin structure.

Several technologies have been developed to identify CREs in a genome-wide manner. For example, DNase-seq and ATAC-seq can identify open chromatin regions that are accessible to TFs1,2. ChIP-seq identifies the binding of specific TFs, co-factors or histone marks by using chromatin immunoprecipitation followed by sequencing3. Cleavage Under Targets and Release Using Nuclease (CUT&RUN) uses an antibody to target specific TFs, co-factors or histone marks followed by the binding of a Protein A/G fused to micrococcal nuclease (pAG-MNase) to cleave the primary antibody-bound sites allowing to profile candidate CREs (cCREs) using a lower cell number than ChIP-seq4,5 (i.e. at least 100 cells for a histone modification and 1,000 cells for a transcription factor are required). CUT&Tag (Cleavage Under Targets and Tagmentation) has been developed from CUT&RUN by using Tn5-mediated tagmentation, generating fragments ready for PCR enrichment and DNA sequencing6. These technologies have identified millions of cCREs across many different cell types and tissues. However, these technologies are descriptive, as the binding of a protein or specific histone mark to a certain DNA sequence does not ultimately mean it is a functional CRE. In addition, it is hard to assess the effect of nucleotide variants in these cCREs on regulatory activity using these technologies.

Massively parallel reporter assays (MPRA) overcome these hurdles, allowing to test the activity of thousands of sequences and variants within them for their regulatory activity by measuring a transcribed barcode7. We previously developed a lentivirus-based MPRA (lentiMPRA), where cCREs are integrated into the host genome and can be tested in a wide-array of cell types8. By testing a similar episomal MPRA library side-by-side, we showed that this ‘in genome’ readout is more strongly correlated with ENCODE annotations and sequence-based models and later work on numerous MPRA technologies showed that it also provides higher cell-type specificity predictions than episomal MPRA9. However, while MPRA can reveal the regulatory activity of a sequence, it cannot identify the proteins binding to that sequence nor its epigenetic modifications.

Here, we develop a technology, e2MPRA, by combining lentiMPRA with CUT&Tag and ATAC-seq (Fig. 1). We used e2MPRA to analyze synthetic enhancers containing liver-specific transcription factors and dissect their epigenetic roles. Additionally, we examined a perturbation library of nine enhancers containing the POU::SOX motif, which is essential for pluripotency, demonstrating the applicability of this technology for identifying key motifs involved in epigenetic function. This technology can systematically characterize thousands of cCREs and their variants for their functional effect on regulatory activity and epigenetic modification, allowing the dissection of the DNA and epigenetic regulatory code side by side.

Results
e2MPRA development
To develop e2MPRA, we first designed a pilot library consisting of 400 elements, each 100 base pairs (bp) in length (Methods), which have been previously tested for regulatory activity using lentiMPRA in HepG2 cells10,11. The pilot library included six distinct sequence categories: 1) 50 random genomic sequences as negative controls; 2) 100 scrambled sequences as negative controls; 3) 100 active CREs from the human genome; 4) 50 inactive genomic sequences; 5) the top 50 and 6) the bottom 50 synthetic sequences by MPRA activity from our previous paper (Fig. 2a). These sequences were synthesized, amplified, and coupled with a random 15-bp barcode before insertion into a lentiMPRA plasmid containing an EGFP reporter vector to generate a sequence-barcode library (Fig. 1a, Supplementary Fig. 1). The resulting pilot library was packaged into lentivirus and transduced into HepG2 cells at a multiplicity of infection (MOI) of 50 (Fig. 1b). To minimize background signal arising from unintegrated lentiviral DNA, infected cells were passaged twice and cultured for 10 days following infection.

Using these cultured cells, we simultaneously performed lentiMPRA, ATAC-seq and CUT&Tag specifically on the library to quantify transcriptional activity and epigenetic modifications in parallel (Fig. 1c, Methods). We conducted a single round of viral infection followed by three independent replicate experiments per assay, which served as technical replicates for downstream analyses.

We first analyzed the regulatory function of each sequence in the library using MPRAflow12. The regulatory activity of each element was quantified as the ratio of RNA barcode counts to the corresponding DNA barcode counts. Since barcode sequences are randomly associated with each element, we included only barcodes that were observed in both DNA and RNA barcode datasets within each replicate. Additionally, elements with fewer than 5 unique barcodes were excluded from the analysis. After filtering, the final dataset contained on average 79.1 unique barcodes per element, with only one element excluded (Supplementary Fig. 2a). This ensured a robust evaluation of activity and helped correct for site of integration biases. We observed a high correlation for the number of UMIs per barcodes (Spearman’s ρ=0.78–0.89; Supplementary Fig. 2b) and RNA/DNA ratio (Spearman’s ρ>0.98; Supplementary Fig. 2c) across the three replicates. As expected, active genomic CREs and top 50 synthetic sequences exhibited significantly higher regulatory activity compared to random genomic sequences and scrambled sequences (Fig. 2b). Inactive genomic and bottom 50 synthetic sequences showed low activity, recapitulating previous results.

Next, we investigated epigenetic activities of CREs in the pilot library by applying ATAC-seq and H3K27ac CUT&Tag, both well-characterized enhancer marks. In these assays, Tn5-tagmented elements were enriched using primers that specifically amplify sequences in the pilot library and simultaneously labeled with unique molecular identifiers (UMIs) (enriched CRE count; Fig. 1d, e). To estimate the genomic integration frequency for each element, we quantified “non-tagmented” elements from the input genomic DNA (gDNA) using the same primers (inserted CRE count; Fig. 1d, e). We observed an average of 95.9 and 68.8 UMIs per CRE in ATAC-seq and H3K27ac CUT&Tag, respectively, while the inserted CRE count had an average of 306.7 UMIs per CRE (Supplementary Fig. 2d). The enriched CRE counts or inserted CRE counts were highly correlated between replicates (Spearman’s ρ=0.83−0.97; Supplementary Fig. 2e). The epigenetic activity of each CRE was quantified as the ratio of enriched CRE count to inserted CRE count (Fig. 1e). Though we obtained slightly better correlation in UMI counts across replicates for ATAC-seq and H3K27ac CUT&Tag (Supplementary Fig. 2e) compared to the lentiMPRA result (Supplementary Fig. 2b), we observed a lower moderate log2(Activity) correlation between replicates (ATAC-seq, Spearman’s ρ=0.64−0.66; H3K27ac CUT&Tag, Spearman’s ρ=0.48−0.53, Supplementary Fig. 2f) than that observed for lentiMPRA (Supplementary Fig. 2c). These moderate correlations are consistent with the way that the present method directly counts CRE fragments using UMIs, leading to greater error propagation when calculating activity, whereas lentiMPRA uses the average of multiple barcodes. To mitigate noise, we summed the counts across replicates and then computed activity by dividing the total enriched count by the total inserted count for each element.

We found that active genomic and top 50 synthetic sequences exhibited significantly higher activity than random genomic sequences and scrambled sequences (Mann–Whitney–Wilcoxon tests with Benjamini–Hochberg correction, adj-P values<0.0001), consistent with the lentiMPRA result (Fig. 2b) with ATAC-seq and H3K27ac CUT&Tag activities across the library correlating with lentiMPRA activity (Spearman’s ρ=0.67 and 0.64, respectively; Fig. 2c). As expected, the bottom 50 synthetic sequences showed lower lentiMPRA activity than active genomic sequences; however, they displayed relatively high ATAC-seq and H3K27ac signals. This observation is consistent with the fact that all synthetic sequences were composed of active TF motifs, and that the bottom 50 sequences are enriched for binding sites of pioneer factors (e.g., ONECUT1, PPARA), rather than transcriptional activators (e.g., HNF1A, XBP1) (Supplementary Fig. 2g).

We observed a moderate correlation between epigenetic activity as measured by e2MPRA and endogenous ATAC-seq and H3K27ac ChIP-seq signals obtained from the ENCODE dataset13 (ATAC, Spearman’s ρ=0.44; H3K27ac, Spearman’s ρ=0.53; Fig. 2d). This could be due to a technical aspect with e2MPRA where specific primers are used to count how many times each element has been tagmented, while genome-wide Cut&Tag and ATAC-seq measure the enrichment based on the number of unique reads in each region. Taken together, our pilot library suggests that although e2MPRA does not quantitatively measure epigenetic activity in a comparable manner to the endogenous genomic context, it can qualitatively capture epigenetic modifications to assess intrinsic regulatory potential.

Dissecting TFs function on epigenetic modifications via e2MPRA
Next, we set out to use e2MPRA to systematically evaluate the effect of TFs on epigenetic modifications. We designed a library of synthetic cCREs that contain nine binding motifs of TFs that are known to be expressed in the liver: CEBPA, CTCF, FOXA1, HNF1A, NR2F2, ONECUT1, PPARA, REST, and XBP111,14 (Fig. 3a). These TF motifs were systematically arranged into three classes on two distinct neutral templates (chr9:83712634-83712733 and chr2:211153273-211153372; hg19) previously shown not to have enhancer activity in HepG2 cells11,14. Class 1 (N=27 per template, total N=54) consisted of homotypic arrangements of 1, 2, or 4 evenly spaced copies of each TF motif, to assess the impact of the number of motifs on their epigenetic activity. Class 2 consisted of heterotypic arrangements of two different motifs to examine their cooperation (N=288). Class 3 comprises heterotypic arrangements of four distinct motifs in all possible combinations to explore further complex interactions (N=6048). We then characterized this library using e2MPRA in HepG2 cells. cCRE activities on the two templates were treated as replicates, yielding two replicates per cCRE, in addition to the three biological replicates of independent library infections.

For the lentiMPRA assay, we observed an average of 157.4 unique barcodes per cCRE (Supplementary Fig. 3a), and the correlations between replicates were high (Spearman’s ρ≈0.9; Supplementary Fig. 3b). For the e2MPRA, we detected on average 406.0 UMIs per cCRE for inserted counts, and 66.5 and 38.8 UMIs per cCRE for enriched counts in ATAC-seq and H3K27ac CUT&Tag, respectively (Supplementary Fig. 3c). The correlations of the activity between replicates for ATAC-seq and H3K27ac were moderate, ranging from 0.4 to 0.5 (Supplementary Fig. 3d), as the e2MPRA data exhibited greater variability compared to lentiMPRA, consistent with our pilot library results. Additionally, sequences containing the NR2F2 motif were rarely amplified or not detected in genomic DNA, likely due to its polyC-like sequence (Supplementary Data 11; see Methods). Therefore, we excluded these data from subsequent analyses, and the remaining eight TF motifs were used for further analyses.

For Class 1, we found that having more CEBPA, FOXA1, HNF1A and XBP1 motifs significantly increased regulatory activity (MPRA), suggesting their primary role as direct transcriptional activators, as shown previously11 (t-test with Benjamini–Hochberg correction, FDR<0.05; Fig. 3b, Supplementary Fig. 4). We found that HNF1A and ONECUT1 motifs were associated with an increase in ATAC signal, fitting with their characterized role as pioneer TFs15,16. The PPARA motif was shown to affect both chromatin accessibility and H3K27ac modification but not transcriptional activation, i.e. having more copies of it did not lead to regulatory activity, consistent with its interaction with chromatin remodeling factors17. Conversely, REST and CTCF, which are known to function as repressors in a context-dependent manner18–20, showed no significant correlation with epigenetic activity (Supplementary Fig. 4).

Next, we compared Class 1 (homotypic arrangements of individual TF motifs) and Class 2 (heterotypic arrangements of two different motifs) to examine synergistic effects of TF motif pairs. For each pair, we performed linear regression to quantify the individual contributions to epigenetic activity, introducing a synergistic interaction term. We observed multiple cooperative effects of two TF motifs on the transcriptional activity measured by lentiMPRA (Fig. 3c, d, Supplementary Fig. 5a), but not for ATAC-seq and H2K27ac activities (Supplementary Fig. 5b, c). For example, PPARA synergistically increased transcriptional activity in combination with CEBPA, CTCF, FOXA1, ONECUT1, and XBP1 (Fig. 3c upper panel), despite its homotypic arrangements increasing only ATAC-seq and H3K27ac activity but not transcriptional activity (Fig. 3b). This suggests that PPARA facilitates transcriptional activation by mediating its epigenetic activity. In addition to the synergistic effects, we also observed repressive effects mediated by REST. REST significantly reduced lentiMPRA activity when combined with any other TF except ONECUT1 (Fig. 3c, d). This suggests that REST represses transcription in the presence of adjacent activators.

To further investigate regulatory grammar, we analyzed Class 3 sequences (heterotypic arrangements of four distinct motifs in all possible combinations) and assessed the effects of motif order on regulatory activity. For each combination of four motifs (N=70), we performed one-way ANOVA to determine the significance of order on regulatory activity. Among the significant combinations affecting transcriptional activity, the sequence containing four TF motifs [HNF1A, PPARA, REST, XBP1] showed the largest differential regulatory activity depending on TF motif order (Fig. 4a, b). Many of these TF motif combinations showed significant effects on transcriptional activity, but not ATAC-seq or H3K27ac activities, which is consistent with epigenetic modifications occurring in broader regions around cCREs (FDR<0.01; Fig. 4b). In addition, we further investigated positional enrichment of each TF motif in the top 200 and bottom 200 sequences, with position 1-to-4 being distal-to-proximal from the minimal promoter (mP) (Fig. 4c). Specifically, transcriptional activators such as HNF1A and XBP1 showed greater enhancement at positions closer to the mP in the top 200 sequences and depletion in the bottom 200. Conversely, the REST repressor exhibited stronger repressive effect when located closer to the mP. Taken together, these results demonstrate that e2MPRA allows the analysis of epigenetic activity of TFs and their cooperative function independent of their transcriptional activity.

e2MPRA analysis of variant effects
We next set out to test the ability of e2MPRA to simultaneously characterize the effect of variants on regulatory activity and epigenetic marks. We selected five sequences positive for enhancer activity from a previous lentiMPRA10 carried out in induced pluripotent stem cells (iPSCs; WTC11 line). These sequences had ATAC-seq and ChIP-seq peaks for POU5F1/SOX2 in human ESCs and iPSCs, as well as being in close proximity (<1 Mb) to genes associated with pluripotency and/or early development (Fig. 5a, Methods). We generated a single-nucleotide substitution library (Fig. 5b), where each nucleotide within the 100-bp sequence was individually mutated to each of the three alternative nucleotides (N=300 variants per CRE). We then performed e2MPRA using WTC11 iPSCs to quantify the impact of these variants on transcriptional and epigenetic activities.

For the lentiMPRA, we observed an average of 211.8 unique barcodes per CRE (Supplementary Fig. 6a), and the activity of variants was highly correlated between replicates (Spearman’s ρ≈0.95; Supplementary Fig. 6b). For the ATAC-seq and CUT&Tag H3K27ac, we detected on average 754.4 UMIs per element for inserted counts, and 43.3 and 56.5 UMIs per element for enriched counts, respectively (Supplementary Fig. 6c). The correlations between replicates were lower for both ATAC-seq and H3K27ac assays (Spearman’s ρ≈0.4–0.5; Supplementary Fig. 6d), consistent with previous observations from the pilot and HepG2 libraries. One of the five active CREs (seq6846_R; Fig. 5a) yielded insufficient read coverage and was therefore excluded from downstream analyses. The remaining four CREs showed high lentiMPRA activities (log2(Activity)=0.66–4.79) for the wild type sequence fitting with their selection as active enhancers (Fig. 5a). Thus, we focused our subsequent analyses on these four CREs. For the single-nucleotide substitution libraries, we quantified the effects of each variant using a linear regression model and found positional clusters of variant effects (Fig. 5c, d, Supplementary Fig. 7). For example, CRE seq68781_R, located near the FOXK1 gene, contains ETS and POU5F1::SOX2 motifs, and their mutations significantly decreased its activity (Fig. 5c).

Analyzing ATAC-seq and H3K27ac activities, we observed noisy signals and no clusters of variant effects around putative TF motifs in any CREs (Fig. 5c, d, Supplementary Fig. 7). This could be due to e2MPRA directly using read counts of amplified variants without using barcode counts, thus having weaker statistical power. To overcome this, we introduced 6 bp randomization with 1 bp-sliding windows for each CRE (two randomizations per window, total N=190 per CRE; Fig. 5b). To quantify the impact of these perturbations, we calculated the perturbation effect at each position (k) as the median epigenetic activity across all sequences affected by perturbation at position k, then computed the median absolute deviation (MAD) score using wild-type (WT) enhancer activity as the baseline (set to 0). Finally, we applied a modified Canny edge detection algorithm to the position vs. MAD score data, smoothing and identifying regions with sharp epigenetic activity changes upon perturbation as significant functional peaks (Methods). This approach found that disruptions of POU5F1::SOX2 and ETS motifs in the seq68781_R (Fig. 5e) significantly decreased ATAC-seq and H3K27ac activities, consistent with the single-nucleotide mutation results (Fig. 5c). In the previously characterized POU5F1 enhancer21 (POU5F1_DE_core), we found that disruption of POU5F1::SOX2 motif decreased both regulatory activity and ATAC signals (Fig. 5f). In addition, while the disruption of the YY1 motif increased regulatory activity, it reduced ATAC activity (no significance for H3K27ac), suggesting that the transcriptional repression by YY1 is associated with chromatin opening mediated by its own binding. It was previously shown that both POU5F1::SOX2 and YY1 motifs are conserved among species and required for enhancer activity21. Overall, our results demonstrate that e2MPRA effectively enables high-resolution identification of functional motifs, providing a powerful framework for systematic dissection of regulatory logic.

Methods
MPRAlibrarydesign
HepG2pilotlibrary.To design the pilot library, we selected sequence features listed in Fig. 2a. 100 active CREs were selected from the top 100 sequences in our previous lentiMPRA in HepG210, as well as the overlap with H3K27ac peaks (ENCODE data: ENCFF001SWK). 50 inactive genomic sequences were selected as the bottom 50 sequences.
For synthetic enhancers, we selected the top 50 and bottom 50 sequences based on MPRA signals from the library used in Smith et al.11, which consists of combinations of 12 active TF motifs. 50 random genomic sequences were selected from genomic sequences that do not overlap with out H3K27ac, H3K27me3 marks, and RepeatMasker annotation. Scramble sequences were generated by randomizing nucleotides using 100 active genomic sequences. Each CRE was 200 bp in length, and flanked by 15 bp adapter sequences (5’-AGGACCGGATCAACT and CATTGCGTGAACCGA-3’).
To generate 100 bp and 150 bp libraries, the initial 200 bp CRE sequences were cropped from the center of each sequence and added another adapter sequences (5’-AATGCTAGCGCATGG and CTGCAACCTACGGAA-3’). We tested all the 100 bp, 150 bp and 200 bp libraries with e2MPRA. However, we were not able to effectively obtain PCR amplification followed by ATAC or CUT&Tag for the 150 bp and 200 bp libraries. Therefore, we decided to use the 100 bp library in the following experiments.

HepG2syntheticenhancerlibrary.Seven TF binding motifs (CEBPA, FOXA1, HNF1A, NR2F2, ONECUT1, PPARA, and XBP1) used in Smith et al.11 were selected based on their correlation significance with transcriptional activity. In addition, we included the CTCF and REST sites used in Georgakopoulos-Soares et al.14, resulting in a total of nine TF motifs in this library. These motifs were arranged in various combinations (Class 1, 2, 3; Fig. 3a) on two neutral DNA templates (hg19 chr9:83712634-83712733 and chr2:211153273-211153372), which were the same as those used in the previous study but cropped 100 bp from the center. In addition to these synthetic enhancer sequences, we also included 200 control sequences selected from four sequence features of the pilot library based on lentiMPRA activity: random genomic sequences, scrambled sequences, and active/inactive genomic sequences (50 each). These control sequences were used only for TMM normalization.

WTC11enhancerperturbationlibrary.To select enhancer candidates for perturbation in the WTC11 iPSC line, we first collected POU5F1 and SOX2 peaks from TF ChIP-seq data and H3K27ac ChIP-seq peaks from ChIP-Atlas29 (last accessed: 2023-05-11). These peaks were overlapped with functional WTC11 CREs identified previously10, and the intersected CREs were scanned for POU5F1:SOX2 binding motifs using FIMO30 ver. 5.5.1, resulting in 1000 candidates. Among these candidates, we selected five active (seq6846_R, seq68781_R, seq5934_F, POU5F1_DE_core, NANOG_p) located near genes of interest within 1 Mb and four inactive sequences (seq16767_R, seq34039_F, seq34899_R, seq2846_R). The inactive sequences were included to identify potential variants that relieve repressive elements and lead to their activation. However, their MPRA activity scores were too weak and therefore excluded from subsequent analyses (Supplementary Data 8). We trimmed these sequences to 100 bp, ensuring that the POU5F1:SOX2 motif is included. To generate the single nucleotide substituted subset, we created sequences in which each nucleotide at every position of the target sequence was replaced with each of the other three nucleotides.
To generate the 6-bp window perturbation subset, we created sequences in which random mutations were introduced within 6-bp windows, shifting by 1 bp at a time. To prevent the emergence of TF motifs due to random mutations, we extracted a 6-bp region upstream and downstream of the mutated site (totaling 18 bp) and matched with known binding motifs using FIMO. If known motifs were detected, the introduced mutation was rejected and a new mutation process was performed recursively. The mutation process was repeated twice, and both replicates were included in the library. We also included the same 200 control sequences from the HepG2 library which were also used only for TMM normalization.

GenerationofMPRAlibraries
The MPRA libraries were generated as previously described12, with minor modifications. In brief, all library sequences were synthesized as a Twist oligonucleotide pool with adapter sequences at both ends (adapter02 for the pilot library and the WTC11 library; adapter03 for the HepG2 library; Supplementary Fig. 1a). The oligonucleotide pool was amplified by two rounds of PCR using NEBNext High-Fidelity 2X PCR Master Mix (NEB; Supplementary Fig. 1b). The first-round PCR was performed for 10 cycles using the primer set 5BC-AG02/03-f01 and 5BC-AG02/03-r01 (Supplementary Data 11) to attach a minimal promoter. The second-round PCR was performed for 10 cycles using the primer set 5BC-AG-f02 and 5BC-AG-r02 (Supplementary Data 11) to attach 15-bp random barcodes for lentiMPRA. For both PCR reactions, the following cycling program was used: 98°C for 2 min; 10 cycles of 98°C for 15 s, 60°C for 20 s, and 72°C for 30 s; followed by a final extension at 72°C for 5 min. The amplified fragments were then inserted into the SbfI/AgeI site of the pLS-SceI vector (Addgene, #137725) using the NEBuilder HiFi DNA Assembly mix (NEB). The recombinant products were transformed into electrocompetent cells (NEB) following the manufacturer’s protocol and incubated overnight at 37°C on 15-cm LB agar plates with 100 μL of 100 mg/mL carbenicillin (Nacalai). For each library, we collected ~80,000 colonies for the pilot library, 1.3 million colonies for the HepG2 synthetic enhancer library, and 0.9 million colonies for the WTC11 enhancer perturbation library, aiming to obtain on average 200 barcodes per CRE. To determine the association between CREs and random barcodes, CRE-mP-barcode fragments were PCR-amplified from each plasmid library pool, and P5 and P7 flow cell adapters were attached using P5-pLSmP-ass-i# and P7-pLSmP-ass-gfp (98°C for 1 min; 15 cycles of 98°C for 15 s, 60°C for 20 s, and 72°C for 3 min; 72°C for 5 min) (Supplementary Data 11). The fragments were then sequenced with a iSeq (the pilot library) or NextSeq Mid output 300-cycle kit (the HepG2 and WTC11 library) using custom primers (Read 1: pLSmP-ass02/03-seq-R1; Index read 1: pLSmP-ass-seq-ind1; Index read 2: pLSmP-rand-ind2; Read 2: pLSmP-ass-seq02/03-R2, Supplementary Data 11).

Lentiviruspackagingandtitration
Lentivirus packaging, titration, and infection were performed as previously described12, with minor modifications. Lentivirus for each plasmid library was produced by co-transfecting it with the helper plasmids pMD2.G (Addgene, #12259) and psPAX2 (Addgene, #12260) into HEK293T (CRL-3216, ATCC) cells cultured in T175 flasks using the EndoFectin Lenti transfection reagent (GeneCopoeia), according to the manufacturer’s protocol. At 8 h post-transfection, the culture medium was refreshed, and ViralBoost reagent (Alstem) was added. At 48 h after adding the ViralBoost reagent, the culture supernatant was collected, filtered through a 0.45-µm PES filter unit, concentrated 50-fold using the Lenti-X concentrator (Takara), and stored at 4°C for up to three weeks. The purified lentivirus was then used for titration. For the pilot library and the HepG2 library, HepG2 cells (HB-8065, ATCC) were plated at 1 million cells per well in 6-well plates and incubated for 24 hours. Then, serial volumes (0, 4, 8, 16, 32, 64 µL) of the lentivirus were added along with 8 µg/mL Polybrene. For the WTC11 library, WTC11 cells (P59; GM25256, Coriell Institute) were plated in 6-well plates at 20% confluency and incubated for 24 h. Then, serial volumes (16, 32 µL) of the lentivirus were added along with an equal volume of ViroMag reagent (OZ Biosciences). All infected cells were cultured for three days and then washed three times with PBS. Genomic DNA was extracted using the Wizard SV Genomic DNA Purification Kit (Promega). The multiplicity of infection (MOI) was measured as the relative amount of viral DNA (WPRE region, primer set: WPRE.F and WPRE.R; backbone, primer set: BB.F and BB.R; Supplementary Data 11) to genomic DNA (intronic region of the LIPC gene, primer set: LP34.F and LP34.R, Supplementary Data 11) by qPCR using the Thunderbird SYBR qPCR Mix (Toyobo), according to the manufacturer’s protocol.

Lentiviralinfectionsandcellcultures
For the pilot library and the HepG2 library, 4 million or 2 million cells per replicate (the pilot library: 1 replicate; the HepG2 library: 3 replicates) were seeded in 2 wells of a 6-well plate and incubated for 24 h. The cells were then infected with the lentiviral libraries along with 8 µg/mL polybrene, with an estimated MOI of 50. For the WTC11 library, 5 million cells per replicate (3 replicates) were seeded in a 10-cm dish and incubated for 24 h. The cells were then infected with 160 µL of the lentiviral library along with 80 µL of ViroMag reagent, with an estimated MOI of 5. To eliminate residual viral DNA within the cells, infected cells were passaged twice after infection and maintained a total of 8–10 days. EGFP fluorescence was observed to confirm that the infected library remained unsilenced after passages. The cells were collected and cryopreserved at −80°C using 100 µL of BamBanker (Nippon Genetics) per 1 million cells. These cryopreserved cells were used for downstream assays.

gDNA/RNAextractionandDNA/RNAbarcodesequencingforlentiMPRAassay
The lentiMPRA assay was performed as previously described12, with minor modifications. In brief, the frozen cells (pilot library: 1 million cells; HepG2 library: 2.5 million cells; WTC11 library: 10 million cells per replicate) were washed twice with PBS and lysed in RLT Plus buffer (pilot library: 100 µL; HepG2 library: 200 µL; WTC11 library: 1 mL; QIAGEN). Genomic DNA and total RNA were extracted from the lysate using the AllPrep DNA/RNA Mini Kit (QIAGEN). To prevent genomic DNA contamination, the Turbo DNA-free Kit (Thermo Fisher) was used for further purification of the RNA solution. Then, the entire purified RNA was used for reverse transcription to generate cDNA using SuperScript IV (Invitrogen). To prepare the sequencing library, 1.5 µg (pilot and HepG2 libraries) or 6 µg (WTC11 library) of genomic DNA and the entire RT product were PCR-amplified for 3 cycles (primer: P5-pLSmP-5bc-i#/P7-pLSmP-ass16UMI-gfp, Supplementary Data 11; cycling program: 98°C for 1 min; 3 cycles of 98°C for 10 s, 60°C for 30 s, and 72°C for 1 min; 72°C for 5 min) to attach P5 and P7 flow cell adapters and UMIs to the sequencing library. The amplicons were further amplified 18 cycles using the primer set P5/P7. The final product was pooled the DNA and RNA barcode library at a molar ratio of 1:3 and sequenced on a iSeq (the pilot library) or NextSeq High output 75-cycle kit (the HepG2 and WTC11 library) with custom primers (Read 1: pLSmP-ass-seq-ind1, Read 2: pLSmP-bc-seq, i7 index: pLSmP-UMI-seq, i5 index: pLSmP-5bc-seq-R2) (Supplementary Data 11).

ATACandCUT&Tagfore2MPRAassay
ATAC-seq31 (Illumina Tagment DNA Enzyme and Buffer Small Kit; Illumina) and CUT&Tag (Hyperactive In-Situ ChIP Library Prep Kit for Illumina; Vazyme) were performed on frozen cells following the manufacturer’s protocol with minor modifications. For both assays, we used 500,000 cells (pilot and HepG2 libraries) or 2 million cells (WTC11 library). The following procedure was performed for 500,000 cells. For 2 million cells, all reagent volumes were increased fourfold accordingly. For the ATAC assay, the cells were washed with PBS and suspended in 50 µL of lysis buffer for nuclei isolation. The isolated nuclei were then resuspended in 50 µL of transposition mix per 500,000 cells, incubated at 37°C for 30 minutes, and purified using the MinElute Reaction Cleanup Kit (QIAGEN).
For the CUT&Tag assay, the cells were washed with 250 µL of wash buffer, then resuspended in 50 µL of wash buffer and 5 µL of ConA bead solution. After treatment with ConA beads, the cells were washed again and resuspended in 50 µL of antibody buffer, followed by the addition of 0.5 µL of primary antibody and 0.5 µL of secondary antibody. We used H3K27ac (Abcam: ab4729) as the primary antibody and goat anti-rabbit IgG (Abcam: ab6702) as the secondary antibody. The primary antibody incubation was performed overnight at 4°C, while the secondary antibody was incubated for 2 h at room temperature. After incubation, the solution was washed and resuspended in 150 µL of the tagmentation buffer, then incubated at 37°C for 1 hour. The tagmented DNA was extracted by stopping the reaction with 5 µL of 0.5 M EDTA, 1.5 µL of 10% SDS, and 1.25 µL of Proteinase K, followed by phenol/chloroform extraction and ethanol precipitation.
The purified Tn5-cleaved products from both assays underwent size selection using 0.65×/1.2× SPRIselect (Beckman Coulter) to remove genomic DNA contamination. To prepare the sequencing library, half of the total products were amplified for 3 cycles in the first-round PCR (primers: P5-ctMPRA-amp02/03-i# and P7-ctMPRA-amp02/03-UMI (Supplementary Data 11; Supplementary Fig. 1c); cycling program: 98°C for 1 min; 3 cycles of 98°C for 10 s, 68°C for 30 s, and 72°C for 1 min; final extension at 72°C for 5 min) to attach P5 and P7 flow cell adapters and UMIs. The amplicons were then purified using 1.0× AMPure XP (Beckman Coulter). qPCR was performed using the purified product, and the cycle number at which the reaction reached 80% of the plateau was determined as the cycle number for the second-round PCR. The second-round PCR was performed for 17–24 cycles using the primer set P5/P7. To estimate the count of each CRE inserted into genomic DNA, we used 1 µg (the pilot library and the HepG2 library) or 4 µg (the WTC11 library) of gDNA obtained from the lentiMPRA AllPrep extraction and performed the same PCR amplification to construct an additional sequencing library. The sequencing libraries were pooled and sequenced on a iSeq (the pilot library) or NextSeq Mid output 150-cycle kit (the HepG2 and WTC11 library) flow cell with custom primers (Read 1: pLSmP-ass-seq02/03-R1, Read 2: pLSmP-ass-seq02/03-R2, i7 index: ctMPRA-seq02/03-UMI, i5 index: ctMPRA-seq02/03-idx; Supplementary Data 11; Supplementary Fig. 1d) using the following setting: 72+72+15+8 bp.

Sequencingdataprocessingpipeline
lentiMPRAdataprocessing.The analysis of CRE-barcode associations and DNA/RNA barcode sequencing data was performed using MPRAflow v2.3.512 with minor modifications. For CRE-barcode association data, FASTQ files were generated using bcl2fastq with the following parameters: “--no-lane-splitting --minimum-trimmed-read-length 0 --mask-short-adapter-reads 0 --use-bases-mask Y145n,Y15,I10,Y145n”. These FASTQ files were used as input for “association.nf” in MPRAflow, aligning reads to the designed library sequences with the following parameters: ”--min-cov 3 --mapq -1 --cigar 100M”. Since our libraries, particularly the WTC11 library, were highly sensitive to single-nucleotide mismatches, we modified the alignment strategy by replacing “bwa mem” with “bwa aln -n 0” to exclude all reads containing any mismatches. For DNA/RNA barcode quantification, FASTQ files were generated from barcode sequencing data using bcl2fastq with the following settings: “--use-bases-mask Y15n,Y16,I10,Y15n”. These FASTQ files were then processed with count.nf in MPRAflow using the parameters: “--bc-length 15 --umi-length 16 --thresh 5”. The resulting DNA/RNA barcode count files for each replicate were subsequently used in downstream analyses.

ATACandCUT&Tagdataprocessing.For ATAC and CUT&Tag sequencing data, FASTQ files were generated using bcl2fastq with the following parameters: “--no-lane-splitting --minimum-trimmed-read-length 0 --mask-short-adapter-reads 0 --use-bases-mask Y72,Y15,I8,Y72”. This command generated paired-end reads of CREs as R1 and R3, and UMI reads as R2. Using these FASTQ files, we first associated the CRE with UMIs using fastp v0.23.232 with the following parameters: “-i $R1/$R3 -I $R2 --umi --umi_loc=read2 --umi_len=15 -w 1 -Q -A -L -G -u 100 -n $UMI_LENGTH -Y 100”. This command appended the UMI to the end of the first part of the read header for R1/R3. Next, a consensus sequence was generated from the paired-end reads using fastq-join v1.3.133. All consensus sequences were then aligned to the designed sequences using BWA34 with the parameter: “bwa aln -n 0”. The resulting BAM file was processed to remove PCR duplicates using UMI-tools v1.1.435 with the following parameters: “dedup --per-gene --per-contig --umi-separator=“:””. Finally, UMIs per CRE were counted from the deduplicated BAM file using samtools36 idxstats (v1.11). This process was performed separately for each replicate.

Replicates,normalizationandactivityscores.Activity scores were calculated by dividing the enriched counts of CREs (or barcodes) from each assay by the counts detected from genomic DNA. For lentiMPRA data, the activity score for each CRE was calculated as the ratio of “rna_count” to “dna_count” of barcodes. For ATAC and CUT&Tag data, activity scores were determined by dividing the number of UMIs per CRE enriched by the number of UMIs per CRE amplified from genomic DNA. To ensure comparability of activity scores across different assays, replicates, and libraries, we first performed counts per million (CPM) normalization. We then calculated a scaling factor f using the trimmed mean of M-values (TMM) normalization method37, based on the enriched counts and genomic counts of random genomic sequences, which served as a negative control. We assumed that these sequences exhibited no variability between the enriched counts from assays and the genomic counts. All enriched counts for the library sequences were then divided by the calculated scaling factor f. This normalization process ensured that the activity scores were adjusted so that the expected log2(Activity) of the negative control sequences was 0 across all samples. To combine the data from all three replicates, we followed the method described by Wang et al.38, which first summarizes the enriched and genomic counts within each replicate and then calculates the activity score by dividing the enriched count by the genomic count. For regression analysis, we used the activity scores from each replicate.

EndogenousactivityofCREsinthepilotlibrary
The endogenous activities of CREs from genomic sequences in the pilot library were estimated using whole-genome epigenetic assay data from ENCODE13. We followed our previous method (Agarwal et al.)10,11 to obtain and compute the consensus signal from these data. In brief, we extracted three bigWig files for H3K27ac ChIP-seq data (ENCFF084DIM, ENCFF515WSE, ENCFF759SNY) and five bigWig files for ATAC-seq data (ENCFF622FRD, ENCFF024GLW, ENCFF240VVR, ENCFF782GKX, ENCFF029XKY) of HepG2 cells from the hg38 assembly in ENCODE. For each CRE, we first extended the original 100 bp genomic region to 500 bp to mitigate positional biases and then calculated the mean bigWig signal from the corresponding genomic region using bigWigAverageOverBed39. All data were log-transformed, and multiple replicates corresponding to the same CRE were averaged to compute the consensus signal for each assay. These signals were then compared with the log2(Activity) of e2MPRA.

AnalysisofHepG2syntheticenhancerlibrary
All analyses were performed using in-house Python scripts based on statsmodels40 v0.14.2 and scikit-learn41 v1.2.2. The difference between the two template types was considered as an additional replicate, resulting in a total of 2×3 replicates for each sequence. Of the nine TF motifs used, sequences containing NR2F2 motif were either absent or present at very low counts. From our previous experience, we have repeatedly observed that sequences with long stretches of polyC or polyG bases are often not properly assigned to barcodes during the CRE–barcode association sequencing step. Although the exact cause is unknown, it is likely due to synthesis or PCR errors. In our data, CREs containing a single NR2F2 motif were successfully sequenced, whereas those with two or more NR2F2 motifs were not detected in the lentiMPRA dataset (Supplementary Data 6). Moreover, almost none of the NR2F2-containing sequences were detected in the genomic count, ATAC-seq, or H3K27ac CUT&Tag data (Supplementary Data 7); therefore, we excluded these cases from downstream regression analyses.
To analyze the homotypic amplification of each epigenetic activity score, we calculated Spearman’s correlation and its P-value between log2(Activity) and the count of motifs for each motif and assay. All P-values were adjusted using the Benjamini–Hochberg method42 (FDR=0.05).
To evaluate potential synergistic effects between two TFs, we performed pairwise multiple linear regression analysis for each TF pair. The activity score was modeled as a function of the TF motifs present in each sequence. The observed log2(Activity) of Class 1 and Class 2 sequences was used as the response variable, while the individual TF motif count was used as the predictor variable. Additionally, we introduced an interaction term k, defined as: k={0|Class 1, 1 or 2| Class 2}. The regression model was expressed as follows:
log2(Activity) ~ countTF1 + countTF2 + k (1)
Multiple linear regression was performed for all TF pairs, and the statistical significance of the coefficient of k was assessed by applying multiple comparison corrections to the P-values. Synergistic effects were considered significant if the adjusted P-values met the predefined significance threshold (FDR=0.01).
For Class 3 sequences, we first tested whether changes in TF motif order significantly affected epigenetic activity using one-way analysis of variance (ANOVA). The analysis was performed across all 70 unique motif sets generated by selecting 4 TF motifs out of 8 (8C4=70), with each set comprising 24 permutations (4!=24). For each motif set, one-way ANOVA was conducted using the activity scores of the 24 permutations as input. Statistical significance was determined from ANOVA P-values, followed by Benjamini–Hochberg correction.
To further investigate how the position of each motif relative to the minimal promoter (mP) influences activity, we conducted a positional enrichment analysis. Specifically, we selected the top 200 and bottom 200 sequences based on transcriptional activity from all permutations and calculated the frequency of each TFBS at each of the four positions (position 1=farthest from mP, position 4=closest to mP) as odds ratios. Enrichment was assessed using a hypergeometric test, and P-values were adjusted using the Benjamini–Hochberg method.

AnalysisofWTC11enhancerperturbationlibrary
Single nucleotide substitution subset analysis. To infer the effects of single nucleotide variants, we followed previous studies43,44 with minor modifications. This approach enables formal testing of variant effects and provides p-values that quantify their statistical significance. We fitted a linear regression model of the form:
log2(enriched count) ~ log2(inserted count) + N + offset (2)
where enriched count refers to the number of unique molecular identifiers (UMIs) detected in the RNA, ATAC and CUT&Tag fragments, inserted count refers to the number of UMIs detected in the gDNA, representing the initial insertion frequency, and N represents binary indicators for mutated nucleotides and their positions. The estimated coefficients of N and their P-values were reported as the effects for each variant.

6 bp window perturbation subset analysis. To identify significant functional sites from the 6 bp window perturbation result, we adapted the Canny edge detection method45 to the one-dimensional data. This simplified Canny method was used to detect regions where perturbations induced abrupt changes in activity scores, identifying these regions as functional sites.
i) Positional effect calculation.
For each position i, we calculated the positional effect using the median activity score of six sequences perturbing that position (due to the 6 bp window perturbation design). To assess the statistical significance of the positional effect, we computed the median absolute deviation (MAD) score. The MAD score was calculated using the wild-type (WT) activity score as the mean absolute deviation reference. However, in the ATAC and H3K27ac assays, the number of UMIs for the WT was abnormally higher than that of the variants across all target CREs—around 10,000 UMIs per WT sequence, compared to only hundreds for the variants (Supplementary Data 9). Therefore, instead of using the WT activity score directly, we estimated its reference value by taking the median of single nucleotide substitution results.
ii) Edge detection and functional site identification.
The computed MAD scores were smoothed using a Gaussian filter, and their derivatives were obtained for edge detection. Based on these derivatives, we applied the non-maximum suppression (NMS) technique to enhance edge localization. For hysteresis thresholding, we classified the detected edges as follows:
1. Strong edges: Local maxima identified through NMS and exceeding the median of derivative scores
2. Weak edges: Values exceeding the median of derivative scores.
3. Weak edges connected to strong edges were recursively clustered and reclassified as strong edges.
Through this process, regions with substantial relative changes in activity scores were identified as strong edges.
iii) Peak detection and functional site estimation.
To further suppress noise, we selected MAD score peaks based on the extrema of smoothed values, where absolute values exceeded 0.75. Finally, functional site regions were defined as adjacent regions between peak candidates and strong edges.

Statistics&reproducibility
Statistical analyses for e2MPRA results are detailed in each analysis section of the Methods. Statistical tests for the sequencing data are described in the corresponding figure legends. No statistical method was used to predetermine sample size. For lentiMPRA analyses, only barcodes observed in both DNA and RNA within a replicate were retained, and elements with fewer than 5 unique barcodes were excluded. For the HepG2 synthetic enhancer library, sequences containing the NR2F2 motif were excluded due to absent or very low genomic counts. For the WTC11 enhancer perturbation library, one active CRE (seq6846_R) was excluded due to insufficient read coverage, and downstream analyses were restricted to CREs/variants meeting the read-coverage and detection criteria described in the Results and Methods. The experiments were not randomized. Reproducibility was assessed by comparing replicate measurements (e.g., correlations between replicates). Code and primary sequencing datasets are publicly available as described in the Code availability and Data availability sections.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_2": "Zhang et al present a new approach e2MPRA for simultaneous epigenomic profiling and measuring regulatory activity.\nE2MPRA builds on the authors previous work on lentivirus based MPRA by incorporating epigenetic profiling in addition to\nmeasuring RNA from the integrated sequence. This is an interesting method that could be used to identify factors that lead to\nchanges in chromatin state without activating the regulatory element as exemplified in the present study as well. To show\nthe utility of e2PMRA the authors present two use cases: 1) testing synthetic regulatory elements with binding motifs for nine\ntranscription factors (individual or in combination) in HepG2 cells and 2) the impact of genetic variants on enhancer\nactivity and epigenetics marks at the regulatory element.\nWhile the presented approach provides potentially important insight into CRE function, there are several points that need to\nbe addressed:\nIncluding epigenetic information is major advance of the presented technologies. However, it seems not ideal that a) the\ncorrelation with regular epigenetic profiling and b) reproducibility between replicates (0.4-0.5) is low. How do the authors\nenvision that this could be improved? Is the chromatin state of the elements affected by the integration site in the genome?\nInactive synthetic sequences have higher levels of active chromatin marks than active genomic sequences.\nWhat is the\nunderlying reason? Could this impact the assay performance, e.g. the difference between active and inactive at genomic\nregions seems larger than between active and inactive synthetic regions (Fig. 2d)? Could the difference be increased? Is\nthis part of the reason for the noisy signals when testing the impact of genetic variants (Fig 5c,d)?\nThe authors mention that NR2F2 motif needed to be excluded based on sequence based on high GC content. CTCF motif\nalso has high GC content. Was it also impacted or not?\nRelated to Fig 3:\nHow to explain that ATAC and lentiMPRA for HNF1A is increased but unchanged/decreased levels of H3K27ac? Similar for\nXBP1 increase lentiMPRA leads to unaltered/lower H3K27ac?\nHow do the MPRA results compare to CREs in HEPG2 cells? Do endogeneous CRE show similar changes in epigenetic\nmarks depending on presence of different TF motif and motif combination?\nActivity of TFs combinations with REST is comparable to REST only. However, these interactions are stated as synergistic\neffects (red label). How is the function of e.g. REST and together with the activator CEBPA\nsynergistic? Could a possible\nexplanation be that REST establishes a chromatin state incompatible with activation by the other factors?\nPPARA establish chromatin to facilitate expression. What is the function/mechanistic role of the other factors in the tested\nmodel (no 3D genome architecture, no enhancer-promoter-looping) to induce gene expression?", "reviewer_3": "In this study, Zhang et al. developed a new method they called e2MPRA that utilizes lentivirus-based MPRA to enrich for the\nintegration of specific CREs into the genome followed by Cut&Tag or ATAC-seq targeted specifically for these sequences.\nWhile previous methods focused on either measuring regulatory activity or epigenetic modification, this new method\nmeasures simultaneously regulatory activity, protein binding and epigenetic modification of thousands of candidate CREs.\n1. They first tested the method on 400 100bp enhancers (including positive and negative controls). The results show that\ne2MPRA measures accurately the activity of the enhancer, with ATAC-seq and H3K27ac CUT&Tag activities across the\nlibrary correlating with lentiMPRA activity (but correlation is only 0.64-0.67). In fact, there is a moderate correlation between\nepigenetic activity as measured by e2MPRA and endogenous ATAC-seq and H3K27ac. I think this is not only driven by\ndifference between the epigenetic activity measured by e2MPRA and endogenous methods, but rather, by the fact that\ncombinatorial signal of different epigenetic modifications could better explain enhancer activity than single epigenetic\nmodifications alone (e.g., PMID:34749786) or that other epigenetic modifications beyond H3K27ac could predict\nenhancer activity (e.g., PMIDs:37308596,27089178). Authors could discuss this in the manuscript.\n2. They use e2MPRA to measure binding of different TFs, enhancer activity and epigenetic profiles at these enhnacers.\nThey find that HNF1A and ONECUT1 can potentially open chromatin, which is very interesting and supports reports from the\nliterature. However, this is not the case for FOXA1 (sometimes they refer to this as FOX1A, but I assume that is only a typo).\nThis is interesting since previous studies have reported FOXA1 as a pioneer factor (PMID:27259199), but also some studies\nfound that it might have only moderate capacity to bind closed chromatin and only at the strongest sites (e.g.,\nPMID:37486787). These new results seem to support the latter model, and authors could discuss this in the manuscript.\n3. PPARA seems to facilitate transcriptional activation of other TFs by mediating the enhancer epigenetic activity, while\nREST acts as a repressor and the closer it is to the promoter the stronger the repression is. The latter is known and is nice to\nsee that recapitulated. What maybe the authors could explain more is how REST does not seem to achieve this by closing\nthe chromatin or depleting H3K27ac. Maybe they discussed this, but I missed it.\n4. They tested the effects of variants on activity, ATAC-seq and H3K27ac and found that disruptions of\nPOU5F1::SOX2 and ETS motifs in significantly decreased ATAC-seq and H3K27ac activities, which is interesting and supporting of what is\nexpected for POU5F1::SOX2 and ETS. One aspect that was not clear to me it is why the authors used a linear regression\nmodel (I am not saying it is wrong, just that it was not explained).\nOverall, I think this is a important paper that should be published in Nature Communications, given the authors address my\ncomments above."}, "author_response": "We thank the reviewers for the great comments, which have significantly improved our revised manuscript.\nBelow, we provide a point-by-point response to these reviews.\nReviewer #2\nZhang et al present a new approach e2MPRA for simultaneous epigenomic profiling and measuring\nregulatory activity. E2MPRA builds on the authors previous work on lentivirus based MPRA by incorporating\nepigenetic profiling in addition to measuring RNA from the integrated sequence. This is an interesting\nmethod that could be used to identify factors that lead to changes in chromatin state without activating the\nregulatory element as exemplified in the present study as well. To show the utility of e2PMRA the authors\npresent two use cases: 1) testing synthetic regulatory elements with binding motifs for nine transcription\nfactors (individual or in combination) in HepG2 cells and 2) the impact of genetic variants on enhancer\nactivity and epigenetics marks at the regulatory element.\nWe thank the reviewer for the valuable comments and for recognizing the importance of our work.\nWhile the presented approach provides potentially important insight into CRE function, there are several\npoints that need to be addressed:\nIncluding epigenetic information is major advance of the presented technologies. However, it seems not\nideal that a) the correlation with regular epigenetic profiling and b) reproducibility between replicates (0.4-\n0.5) is low. How do the authors envision that this could be improved? Is the chromatin state of the elements\naffected by the integration site in the genome?\nWe thank the reviewer for this comment. We agree that the concordance of e2MPRA with the endogenous\nsignals and replicate correlations leave room for improvement. In the revised version, we have tried to\naddress these comments. For comment a), the low correlation with regular epigenetic profiling (related to Fig.\n2d) mainly stems from the difference of the methods for how tagmented fragments are detected; i.e. e2MPRA\nuses specific primers to quantify epigenetic modifications by counting how many times each element has\nbeen tagmented. In contrast, genome-wide Cut&Tag and ATAC-seq measure the enrichment based on the\nnumber of unique reads in each region. However, more importantly, the use of specific primers in e2MPRA\nallows us to perform a fair comparison between sequence variants in the library, even though it will not\nquantitatively capture the endogenous signals. To further clarify this, we now mention this issue in our revised\nmanuscript:\n“This could be due to a technical aspect with e2MPRA where specific primers are used to count how many\ntimes each element has been tagmented, while genome-wide Cut&Tag and ATAC-seq measure the\nenrichment based on the number of unique reads in each region.”\nIn addition, we also added the following text to the discussion:\n“We further observed a low correlation between e2MPRA and endogenous epigenetic signals. This could be\ndue to differences in measurement of the two methods, i.e. e2MPRA activities were detected by specific\nprimers and measured in the vector context, while endogenous activities were detected by random\ntagmentation. The use of specific primers in e2MPRA allows for unbiased comparisons between sequence\nvariants within the library. Therefore, e2MPRA could be more suitable for the parallel quantitative comparison\nof enhancer activities among sequence variants rather than for recapitulating endogenous epigenetic\nmodifications.”\nIn terms of the reviewer’s comment b), the low reproducibility between replicates; this could be improved by\nusing easier-toinfect cells (see also WTC11 (Fig. S6d) that is harder to infect and showed lower correlation\nthan HepG2 (Figs. S2f and S3d), which is easier to infect, and adding more sequencing reads and more\nexperimental replicates. To address this, we aggregated barcode counts across replicates as explained in\nthe article:\n“These moderate correlations are consistent with the way that the present method directly counts CRE\nfragments using UMIs, leading to greater error propagation when calculating activity, whereas lentiMPRA\nuses the average of multiple barcodes. To mitigate noise, we summed the counts across replicates and\nthen computed activity by dividing the total enriched count by the total inserted count for each element.”\nIn this way, we were able to detect known and unknown epigenetic functions of transcription factors (Figs.\n3b and 5f); therefore, we are confident that our experimental approach and data analysis methods are\nrobust enough.\nInactive synthetic sequences have higher levels of active chromatin marks than active genomic sequences.\nWhat is the underlying reason? Could this impact the assay performance, e.g. the difference between\nactive and inactive at genomic regions seems larger than between active and inactive synthetic regions\n(Fig. 2d)? Could the difference be increased? Is this part of the reason for the noisy signals when testing\nthe impact of genetic variants (Fig 5c,d)?\nWe apologize for the confusion regarding the category labeling. Both the “active” and “inactive” synthetic\nsequences were taken from our previous study (Smith et al. Nature Genetics 2013) and composed of\nmultiple binding sites of transcription activators. The “active” and “inactive” sets correspond to the top 50\nand the bottom 50 sequences based on the MPRA activity, respectively. Hence, although labeled “inactive”,\nthose sequences still contain TF binding sites and potentially carry epigenetic activities. The figure below\nreflects differences in the compositions of active binding sites between these two subsets. In it, we show for\nexample that the “inactive” sequences include a similar number of PPARA binding sites as “active”\nsequences. This TF is known to be involved in histone modification and pioneering activities and could thus\nlead to open chromatin also in the “inactive” sequences. In contrast, the active sequences contain more\nTFBS copies of HNF1A for example, which is a strong activator in HepG2 cells (Smith et al. Nature\nGenetics 2013).\nTo avoid confusion, we re-labeled these groups as “top 50” and “bottom 50” throughout the manuscript.\nWe should also note that this library is independent from other libraries and does not cause noisy signals or\nimpact on other results.\nWe added the following text in the revised article and also the above figure as Extended Data Fig 2g:\n“As expected, the bottom 50 synthetic sequences showed lower lentiMPRA activity than active genomic\nsequences; however, they displayed relatively high ATAC-seq and H3K27ac signals. This observation is\nconsistent with the fact that all synthetic sequences were composed of active TF motifs, and that the\nbottom 50 sequences are enriched for binding sites of pioneer factors (e.g., ONECUT1, PPARA),\nrather\nthan transcriptional activators (e.g., HNF1A, XBP1) (Extended Data Fig. 2g).”\nThe authors mention that NR2F2 motif needed to be excluded based on sequence based on high GC\ncontent. CTCF motif also has high GC content. Was it also impacted or not?\nIn this study, we unexpectedly failed to retrieve sequences containing NR2F2 motif\n(CCCCCTGACCTTTGCCCCCTGCC) in the library, whereas successfully obtained CTCF\n(CGGCCACCAGGGGGCGCCA) sequences. In our previous experiences, we have repeatedly observed\nthat sequences with many polyC or polyG bases were not properly assigned to barcodes during the CRE–\nbarcode association sequencing step. Although we are not certain of the exact cause, it is likely due to\nsynthetic or PCR errors. Barcode counts per CRE are provided in Supplementary Table S6 (column\nn_obs_bc). In our data, CREs containing a single NR2F2 motif were successfully sequenced, but CREs\nwith two or more NR2F2 motifs were not; therefore, we excluded these cases from downstream regression\nanalyses. CTCF did not have this issue and was included in the analysis. To address this, we have added\nthe following text in the Method section:\n“Of the nine TF motifs used, sequences containing the NR2F2 motif were either absent or present at very\nlow counts. From our previous experience, we have repeatedly observed that sequences with long\nstretches of polyC or polyG bases were not properly assigned to barcodes during the CRE–\nbarcode association sequencing step. Although we are not certain of the exact cause, it is likely due to synthesis or\nPCR errors. In our data, CREs containing a single NR2F2 motif were successfully sequenced, whereas\nthose with two or more NR2F2 motifs were not detected in the lentiMPRA dataset (Supplementary Table\n6). Moreover, almost none of NR2F2 containing sequences were detected in the genomic count, ATAC-seq\nand H3K27ac CUT&Tag (Supplementary Table 7); therefore, we excluded these cases from downstream\nregression analyses.”\nRelated to Fig 3:\nHow to explain that ATAC and lentiMPRA for HNF1A is increased but unchanged/ decreased levels of\nH3K27ac? Similar for XBP1 increase lentiMPRA leads to unaltered/lower H3K27ac?\nWe thank the reviewer for this great comment. Indeed, H3K27ac signals remained unchanged regardless of\nthe number of CEBPA, FOXA1, HNF1A and XBP1 motifs (Fig. 3b, Extended Data Fig. 4). Possible reasons\nthat explain this observation could be: 1) Improving e2MPRA method (adding more reads and more\nreplicates) are required to detect H3K27ac signal. 2) these factors activate gene transcription by directly\ninteracting with gene promoters without being involved in H3K27ac modification.\nTo address this, we have added the following text to the discussion:\n“While CEBPA, FOXA1, and XBP1 increase lentiMPRA activity along with the number of their motifs, we did\nnot observe epigenetic changes for them (Fig. 3b, Extended Data Fig. 4). This could be due to these\nfactors being responsible for activating gene transcription by direct interaction with gene promoters without\nepigenetic modification and/or not having a role related to H3K27ac. Another reason could be due to the\nlimitation of e2MPRA technology, requiring additional reads and replicates to obtain more robust epigenetic\nresults. Thus, we cannot exclude the possibility that these factors are involved in epigenetic modifications.”\nHow do the MPRA results compare to CREs in HEPG2 cells? Do endogeneous CRE show similar changes\nin epigenetic marks depending on presence of different TF motif and motif combination?\nWe thank the reviewer for this suggestion. We have investigated endogenous CREs and found that the\nmajority of CREs consist of combinations of multiple TF motifs, rather than containing a single isolated TF\nmotif. This makes it difficult to quantify the function of individual factors and their cooperative effects in an\nendogenous system.\nThe figure below shows an analysis similar to Fig. 3c but performed on endogenous CREs. We reanalyzed\nthe lentiMPRA data from Vikram et al., Nature 2025, which measured transcriptional activities of\nendogenous CREs in HepG2 cells, and selected sequences that contained at least one of the TFBSs\nexamined in our study. For example, in our e2MPRA, PPARA exhibited strong synergistic effects when\ncombined with other transcription factors (a), whereas in the endogenous data, combinations with other\nfactors showed almost no change in activity distribution compared to CREs containing a single TFBS (b).\nSimilarly, although REST functioned as a strong repressor in our e2MPRA, this effect was not observed in\nthe endogenous data. We think that these observations do not indicate inconsistency between the two\nsystems; rather, they reflect the inherent complexity of endogenous CREs, which contain multiple TFBSs\nand whose activities are highly context dependent. This complexity makes it difficult to isolate and quantify\nthe contribution of each single motif.\nFigure Legend: Synergistic transcriptional effects when combining PPARA or REST motifs with other TF motifs for\nsynthetic constructs (a) and endogenous CREs (b). Transcriptional activities of homotypic arrangements (two identical\nmotifs) were compared with heterotypic arrangements (one motif combined in a 1:1 ratio). For synthetic constructs,\nstatistically significant synergistic effects (FDR < 0.01) are indicated by a red background, and non-significant\ncombinations are shown in blue, as in the main Figure 3c. In Figure 3c, we compared homotypic arrangements of four\nidentical motifs with heterotypic arrangements of two motifs combined in a 2:2 ratio. However, because endogenous\nCREs containing four identical motifs could not be identified, we compared homotypic arrangements of two identical\nmotifs with heterotypic arrangements combining one motif in a 1:1 ratio. The corresponding synthetic constructs were\nanalyzed in the same manner for consistency.\nBy contrast, our synthetic framework places TFBSs on an inactive template sequence, thereby removing\ngenomic contextual effects and enables precise measurement of each factor’s contribution and their\ncombinatorial interactions. We therefore view our results as directly informative for understanding the\ncomplex cooperativity among TFBSs that operates within endogenous CREs.\nTogether, these results highlight the importance of analyzing synthetic enhancers composed of TF motifs to\nunderstand molecular epigenetic function of each TFs using e2MPRA, rather than just investigating\ngenomic enhancers.\nActivity of TFs combinations with REST is comparable to REST only. However, these interactions are stated\nas synergistic effects (red label). How is the function of e.g. REST and together with the activator CEBPA\nsynergistic? Could a possible explanation be that REST establishes a chromatin state incompatible with\nactivation by the other factors?\nWe thank the reviewer for catching this mistake. The term “synergistic” should not be used for REST’s\nrepressive effect. We rephrased a sentence in the Fig.3c legend and the Results section as below:\n“c. Representative examples showing synergistic transcriptional effects when combining PPARA motif with\nother TF motifs, as well as the repressive effect of REST over other motifs (Class 2).”\n“In addition to the synergistic effects, we also observed repressive effects mediated by REST. REST\nsignificantly reduced lentiMPRA activity when combined with any other TF except ONECUT1 (Figs. 3c-d).\nThis suggests that REST represses transcription in the presence of adjacent activators.”\nIn this study, we found that REST functions as a strong repressor even when present alongside a\ntranscriptional activator (Fig. 3c). We also added the following text in the discussion (see also the third\ncomment from the reviewer 3):\n“REST is a well-characterized transcriptional repressor that mediates silencing through recruitment of\nchromatin modifiers, including HDACs, LSD1, G9a/GLP, and PRC1/2. In our assay, REST repressed\ntranscription when present together with transcriptional activators (Fig. 3c), suggesting that it exerted the\nstrongest effect among the factors tested alongside it. By contrast, we did not detect a significant\ndepletion of H3K27ac and closed chromatin for REST (Extended Data Fig. 4). This could be due to a variety of\nfactors, including noisy measurements in this specific e2MPRA. Further technical improvement may be\nrequired to better understand its molecular function.”\nPPARA establish chromatin to facilitate expression. What is the function/mechanistic role of the other\nfactors in the tested model (no 3D genome architecture, no enhancer-promoter-looping) to induce gene\nexpression?\nWe thank the reviewer for this comment. In the previous version of the manuscript, we discussed the\nfunction of other factors in the discussion:\n“We showed that CEBPA, FOXA1, HNF1A, and XBP1 motifs increased transcriptional activity (Extended\nData Fig. 4), consistent with previous results [11,14] and their known roles as transcriptional activators. In\naddition, HNF1A and ONECUT1, previously characterized as pioneer factors [16], were shown to increase\nchromatin accessibility”\nWe have also added the below text in the discussion, according to the reviewer’s suggestion.\n“While CEBPA, FOXA1, and XBP1 increase lentiMPRA activity along with the number of their motifs, we did\nnot observe epigenetic changes for them (Fig. 3b, Extended Data Fig. 4). This could be due to these\nfactors being responsible for activating gene transcription by direct interaction with gene promoters without\nepigenetic modification and/or not having a role related to H3K27ac. Another reason could be due to the\nlimitation of e2MPRA technology, requiring additional reads and replicates to obtain more robust epigenetic\nresults. Thus, we cannot exclude the possibility that these factors are involved in epigenetic modifications.”\nReviewer #3\nIn this study, Zhang et al. developed a new method they called e2MPRA that utilizes lentivirus-based MPRA\nto enrich for the integration of specific CREs into the genome followed by Cut&Tag or ATAC-seq targeted\nspecifically for these sequences. While previous methods focused on either measuring regulatory activity or\nepigenetic modification, this new method measures simultaneously regulatory activity, protein binding and\nepigenetic modification of thousands of candidate CREs.\n1. They first tested the method on 400 100bp enhancers (including positive and negative controls). The\nresults show that e2MPRA measures accurately the activity of the enhancer, with ATAC-seq and H3K27ac\nCUT&Tag activities across the library correlating with lentiMPRA activity (but correlation is only 0.64-0.67).\nIn fact, there is a moderate correlation between epigenetic activity as measured by e2MPRA and\nendogenous ATAC-seq and H3K27ac. I think this is not only driven by difference between the epigenetic\nactivity measured by e2MPRA and endogenous methods, but rather, by the fact that combinatorial signal of\ndifferent epigenetic modifications could better explain enhancer activity than single epigenetic modifications\nalone (e.g., PMID:34749786) or that other epigenetic modifications beyond H3K27ac could predict\nenhancer activity (e.g., PMIDs:37308596,27089178). Authors could discuss this in the manuscript.\nWe thank the reviewer for this suggestion. In our previous study (Agarwal et al. Nature 2025), we built a\nregression model that integrated multiple biochemical features at CREs (TF ChIP-seq, histone ChIP-seq,\nDNase-seq, and ATAC-seq) to predict lentiMPRA activity; the Pearson correlation was around 0.7-0.75\n(Figure 3b in Agarwal et al. Nature 2025). As the reviewer suggests, this supports the idea that composite\nepigenomic features explain activity better than any single mark, while also indicating that epigenomic\nfeatures alone have limited explanatory power. This underscores the value of our present approach, which\nquantifies each epigenomic feature separately. To address this, we added the following text to the revised\ndiscussion:\n“It is generally thought that the combinatorial signal of multiple epigenetic modifications contributes to\nenhancer activity more than single epigenetic modifications alone. Indeed, our previous work integrating\nmultiple biochemical features at CREs (TF ChIP-seq, histone ChIP-seq, DNase-seq, and ATAC-seq)\nenabled to predict lentiMPRA activity to a certain extent; the Pearson correlation was around 0.7-0.75\n(Figure 3b in Agarwal et al. Nature 2025), indicating that epigenomic features alone have limited\nexplanatory power. These findings further emphasize the value of our approach to analyze multiple\nepigenetic modifications separately.”\n2. They use e2MPRA to measure binding of different TFs, enhancer activity and epigenetic profiles at these\nenhancers. They find that HNF1A and ONECUT1 can potentially open chromatin, which is very interesting\nand supports reports from the literature. However, this is not the case for FOXA1 (sometimes they refer to\nthis as FOX1A, but I assume that is only a typo). This is interesting since previous studies have reported\nFOXA1 as a pioneer factor (PMID:27259199), but also some studies found that it might have only moderate\ncapacity to bind closed chromatin and only at the strongest sites (e.g., PMID:37486787). These new results\nseem to support the latter model, and authors could discuss this in the manuscript.\nWe thank the reviewer for this great suggestion! Consistent with this suggestion, the role of FOXA1 as a\npioneer factor remains under active debate. In our assay, the number of FOXA1 sites was positively\ncorrelated with lentiMPRA activity, whereas neither ATAC-seq nor H3K27ac showed a positive correlation.\nThese observations suggest that, under our conditions, FOXA1 behaved primarily as a transcription factor\nrather than as a pioneering factor. That said, its low significance (adjusted p = 0.115; Extended Data Fig. 4)\ncould be due to assay noise; therefore, our data still cannot exclude its pioneer activity. We have added the\nfollowing text in the discussion:\n“The pioneering activity of FOXA1 is still being characterized. While early studies proposed that FOXA1\nacts as a pioneer factor capable of binding to closed chromatin and facilitating subsequent transcriptional\nactivation [23], more recent reports have suggested that its pioneering capacity may be limited to a subset\nof strong binding sites or may depend on cellular context [24]. In our e2MPRA assay, FOXA1 increased\ntranscriptional activity without a corresponding increase in chromatin accessibility or H3K27ac modification\n(adjusted p = 0.115 in ATAC-seq; Extended Data Fig. 4), suggesting that under our experimental conditions,\nFOXA1 primarily functioned as a transcriptional activator rather than as a pioneer factor. Nevertheless,\ngiven the moderate statistical significance of the observed effects, we cannot rule out a pioneering activity\nfor FOXA1.”\nWe also thank the reviewer for catching FOXA1 typos. All typos in the text are now fixed.\n3. PPARA seems to facilitate transcriptional activation of other TFs by mediating the enhancer epigenetic\nactivity, while REST acts as a repressor and the closer it is to the promoter the stronger the repression is.\nThe latter is known and is nice to see that recapitulated. What maybe the authors could explain more is\nhow REST does not seem to achieve this by closing the chromatin or depleting H3K27ac. Maybe they\ndiscussed this, but I missed it.\nAs the reviewer mentioned, REST represses transcription by recruiting chromatin modifiers, including\nHDACs, LSD1, G9a/GLP, and PRC1/2, which can lead to depletion of H3K27ac [25, 27]. In our assay,\nREST repressed transcription when present together with other transcription factors on a CRE, suggesting\nthat it exerted the strongest effect among the factors tested. By contrast, we did not detect a significant\ndepletion of H3K27ac, likely because the H3K27ac measurements were particularly noisy in this dataset.\nWe have added the following text to the discussion:\n“REST is a well-characterized transcriptional repressor that mediates silencing through recruitment of\nchromatin modifiers, including HDACs, LSD1, G9a/GLP, and PRC1/2 [25, 27]. In our assay, REST\nrepressed transcription when present together with transcriptional activators (Fig. 3c), suggesting that it\nexerted the strongest effect among the factors tested alongside it. By contrast, we did not detect a\nsignificant depletion of H3K27ac and closed chromatin for REST (Extended Data Fig. 4). This could be\ndue to a variety of factors, including noisy measurements in this specific e2MPRA. Further technical\nimprovement may be required to better understand its molecular function.”\n4. They tested the effects of variants on activity, ATAC-seq and H3K27ac and found that disruptions of\nPOU5F1::SOX2 and ETS motifs in significantly decreased ATAC-seq and H3K27ac activities, which is\ninteresting and supporting of what is expected for POU5F1::SOX2 and ETS. One aspect that was not clear\nto me it is why the authors used a linear regression model (I am not saying it is wrong, just that it was not\nexplained).\nOur primary reason for using linear regression is continuity with our prior work, in which we estimated\nvariant effects from saturation mutagenesis using linear models (Kreimer et al. 2022, Nat. Comm., Ashuach\net al. 2019, Genome Biol.). Linear regression also permits formal testing of the null hypothesis of no variant\neffect and provides p-values that quantify statistical support. This makes it straightforward to report, for\neach variant, the strength of evidence for its effect. To address this, we revised the following text in the\nMethod section:\n“To infer the effects of single nucleotide variants, we followed previous studies [45, 46] with minor\nmodifications. This approach enables formal testing of variant effects and provides p-values that quantify\ntheir statistical significance.”\nOverall, I think this is a important paper that should be published in Nature Communications, given the\nauthors address my comments above.\nWe thank the reviewer for the supportive comments."}, {"round": 2, "reviewer_comments": {"reviewer_2": "The authors have addressed all my comments in the revision. I recommend publication of the manuscript.", "reviewer_3": "The authors addressed my comments."}, "author_response": ""}]

Decision letter