Nucleotide-resolution mapping of regulatory elements via allelic readout of tiled base editing

Domain: chemistry
Published: 2026-03-03
Editorial outcome: minor_revision

Abstract
CRISPR tiling screens have enabled the characterization of regulatory sequences but are limited by low resolution arising from the indirect readout of editing via guide RNA sequencing and enrichment analysis. This study introduces an end-to-end experimental assay and computational pipeline, which leverages targeted sequencing of CRISPR-introduced alleles at the endogenous target locus following dense base-editing mutagenesis. As a proof of concept, we studied a putative CD19 enhancer, an immunotherapy target in leukemia, and identified alleles and single nucleotides crucial for CD19 regulation. Our visualization tools revealed transcription factor motifs corresponding to the top-ranked nucleotides. Validation experiments confirmed that mutations in MYB, PAX5, and EBF1 binding sites reduce CD19 expression. Critically, editing MYB and PAX5 motifs conferred resistance to aCD19 CAR-T cell therapy, revealing how noncoding variants can drive immunotherapy escape. Taken together, this approach achieves nucleotide-resolution genotype-phenotype mapping at regulatory elements beyond conventional gRNA-based screens.

Introduction
While chromatin profiling, GWAS, and eQTL analyses propelled significant breakthroughs in identifying potential functional non-coding regions and variants in disease1, these methods do not provide mechanistic understanding of regulatory elements at the nucleotide level. In this context, CRISPR-based methods facilitate in situ perturbations of regulatory sequences, presenting opportunities to explore regulatory functions across various developmental stages and in numerous disease states in their native chromatin context2–6. Current CRISPR screens commonly perform indirect readout of perturbations via guide RNA sequencing7. However, this readout is largely confounded by the editing efficiency of each gRNA, and the enrichment of specific editing outcomes cannot be resolved. Initial attempts were made to uncover allelic contribution to phenotype by comparing disruptive and non-disruptive alleles at regulatory sequences2,8. These experiments however had limited resolution (~10bp), inherent to the repair profiles of nucleases, and low throughput. Further methods to link nucleotide variants to phenotype have been developed but focus mainly on coding variants and rely on scRNAseq for phenotypic read-out9–11. While prime- and base-editing screens provide high-precision interrogation of individual regulatory variants12,13, they remain constrained by editing efficiency and a lack of computational methods to accurately identify the causal variants based on mutational outcomes. While Cooper et al.14, and Nam et al.15 combine scRNAseq with targeted single cell genotyping, these studies did not investigate the effect of non-coding variants on gene expression. As these methods rely on single-cell technologies, high experimental costs impose significant limitations on their scalability.

To address these limitations, we introduce an experimental and computational framework that enables precise determination of genotype-phenotype relationships at the nucleotide level from bulk amplicon sequencing, in the context of non-coding sequences, utilizing information from direct sequencing of alleles produced by base editing at regulatory regions. To interrogate cellular phenotypes, our approach is coupled with fluorescence-activated cell sorting, allowing to link variants to a variety of cellular phenotypes.

As a proof-of-concept, we tile a CD19 enhancer using evoCDSA and ABE8e-SpRY base editors and directly compare gRNA-based and allele-based workflows. While we observe shared hit regions between gRNA-based and allele-based approaches, sequencing of individual alleles uncovers additional hits. Leveraging the information from the identified single nucleotide substitution hits, we identify disruptions of MYB, PAX5 and EBF1 binding motifs and further demonstrate their trans-regulatory functions. Finally, we show that our prioritized enhancer mutations suffice to and render cells resistant to aCD19 CAR T therapy, highlighting the importance in the identification of disease-linked causal variants.

Results
High-resolution dissection of regulatory elements
To demonstrate the efficacy of our framework, we investigated a putative regulatory element upstream of the CD19 gene, a B-cell marker and immunotherapy target in pediatric and adult B-cell acute lymphoblastic leukemia (B-ALL) and lymphoma. The candidate regulatory element was identified based on the presence of a highly evolutionary conserved and open chromatin region in proximity to the CD19 gene (Fig. 1a), as indicated by ATAC-seq and H3K27ac ChIP-seq16. In order to validate this enhancer, we isolated clones with mono- or biallelic deletions of the 346bp element. These heterozygous and homozygous clones showed reduced CD19 mRNA (Fig. 1b) and protein levels (Fig. 1c), confirming the regulatory potential of this enhancer element. To further dissect the enhancer, we performed base editor screens, for which we designed a library of 200 gRNAs tiling the CD19 enhancer and used both cytosine (evoCDA17) and adenine (ABE8e-SpRY18) base editors in NALM6 cells, a pediatric B-ALL cell line (Fig. 1d). The ABE8e-SpRY base editor operates with relaxed PAM requirement, therefore installing edits at high density at the CD19 enhancer (Supplementary Fig. 1a), in contrast to the evoCDA base editor which is restricted by the -NGG proximal sequences (Supplementary Fig. 1b). While ABE8e performs A-to-G substitutions (Supplementary Fig. 1c), evoCDA is a spurious editor primarily introducing C-to-T substitutions along with different types of C-to-N substitutions and short indels, thereby increasing allele diversity (Supplementary Fig. 1d).

After transducing the gRNAs at a 0.3 MOI, we cultured the cells for 7 days to allow for efficient editing and phenotype development. We then stained cells with a fluorescently conjugated antibody and performed fluorescence-activated cell sorting (FACS) in triplicates to sort cells into “CD19 positive” and “CD19 negative” populations. For gating and sorting, we targeted the AAVS1 locus as negative control and the CD19 exon 2 splice acceptor site as positive control (Fig. 1e, Supplementary Fig. 1e). For both the CD19 positive and CD19 negative sorted samples, we split the extracted genomic DNA for both targeted sequencing of the gRNA cassette and the endogenously targeted region. This approach enabled direct comparison between standard gRNA count-based CRISPR screens analyzed with CRISPR-SURF19 (an established computational tool for identifying functional regulatory regions from gRNA enrichment data) and direct allele sequencing. We developed CRISPR-Millipede to achieve single-nucleotide resolution and adapted DESeq2 to identify functionally significant alleles20.

To address the "self-editing" at the gRNA cassette, a feature of base editors with loose PAM requirements and confounding factor in pooled screens, we also developed CRISPR-Correct, a relaxed gRNA mapping method that corrects for substitutions at the gRNA cassette, increases gRNA count recovery and power, and further captures gRNA mismapping due to sequencing or amplification errors. CRISPR-Correct increases the read count of gRNAs by an average of 26.5%-41.5% across all ABE8e-SpRY CD19 positive and CD19 negative samples (Supplementary Fig. 2a) and 13.6%-16.6% across all evoCDA CD19 positive and CD19 negative samples (Supplementary Fig. 2b). We observed an average of 20.3%-22.9% self-editing frequency in ABE8e-SpRY presort samples (Supplementary Fig. 2c) primarily consisting of A>G self-edits (Supplementary Fig. 2d). Additionally, we detected a direct correlation between gRNA self-editing frequency and post-correction percent increase (Supplementary Fig. 2e). We then applied CRISPR-SURF19 to the corrected gRNA count data for both ABE8e-SpRY and evoCDA datasets to establish a baseline for comparison and evaluate the efficacy of standard gRNA count-based CRISPR screens in identifying functional sub regions within the CD19 enhancer supported by multiple gRNAs.

Next, we focused on allele-based readouts to overcome the limitations of gRNA count-based methods in elucidating genotype-phenotype relation with nucleotide-level precision. First, we developed CRISPR-Millipede to enable fine-mapping and identification of single nucleotide causal variants, and second, adapted the DESeq2 method to identify significant alleles (Fig. 1f). CRISPR-Millipede represents a Bayesian linear regression and feature selection framework that we designed to effectively distinguish functional nucleotides from bystander edits across both ABE8e-SpRY and evoCDA screens. This regression model, which leverages the recent “millipede” statistical tool21 for Bayesian Variable Selection, employs an efficient Markov Chain Monte Carlo (MCMC) scheme for inference, allowing it to scale to scenarios with a large number of potential regulatory nucleotides and assigns both effect sizes (beta coefficient) and posterior inclusion probabilities (PIP) scores to each observed substitution under an additive assumption, based on their presence and frequency in the observed alleles across sorted populations.

CRISPR-Millipede leverages CRISPResso222 to merge, align, and quantify edited alleles. The resulting alleles are subsequently encoded into a feature indicator matrix, focusing exclusively on canonical edits: A>G and T>C transitions for ABE8e-SpRY screens, and G>A and C>T transitions for evoCDA screens. To capture the diversity of editing imposed by evoCDA, we additionally performed CRISPR-Millipede analysis incorporating all edit types, including single base deletions, in the encoding matrix. With these computational tools established, we next applied them to identify functional regulatory elements within the CD19 enhancer.

Pinpointing functional subregions, alleles, and nucleotides
We identified regions, alleles and nucleotides crucial for CD19 enhancer function by comparing CD19 negative and CD19 positive populations and employing CRISPR-Correct followed by CRISPR-SURF, DESeq2 and CRISPR-Millipede. After calculating the enrichment of gRNA read counts between the CD19 positive and CD19 negative populations, CRISPR-SURF identified one significant hotspot region at amplicon position 218-235 that was shared between evoCDA and ABE8e-SpRY screens. CRISPR-SURF also identified a hit sub-region unique to the evoCDA screen at amplicon position 145-166 and a hit unique to the ABE8e-SpRY screen at amplicon position 326-330 (Fig. 2a). When CRISPR-SURF was run without mismatch-correction by CRISPR-Correct, the evoCDA screen yielded only a single hit at position 145-166, and the ABE8e-SpRY screen produced no detectable hits (Supplementary Fig. 3b). Additionally, supporting gRNAs within the CRISPR-SURF sub-regions show low correlation of log fold-change enrichment scores across replicates without CRISPR-Correct (r=-0.08 to 0.45) and improved correlation when using CRISPR-Correct in the ABE8e-SpRY screen (r=0.62 to 0.78) (Supplementary Fig. 3a), highlighting the necessity of correcting for mismatches from self-editing prior to downstream analyses.

In contrast to the gRNA readout model, CRISPR-Millipede demonstrates superior capabilities in elucidating genotype-phenotype relationship with nucleotide-level precision, a significant improvement over the ~5-10 base pair resolution of traditional genetic screens. While gRNA-based analysis is biased by differences in the editing frequency of individual gRNAs and underestimates the signal from low-editing gRNAs, CRISPR-Millipede captures editing rates, and identifies regulatory potential based on relative allele frequencies (Fig. 2b). This awareness of editing efficiency also enables more accurate estimation of effect sizes and allows the method to disambiguate phenotypic effects more precisely than gRNA-based methods. CRISPR-Millipede analysis was performed at three levels: (1) individual replicates, (2) within-editor joint models combining all replicates for each editor separately, and (3) a cross-editor joint model integrating data from both ABE8e-SpRY and evoCDA screens. To visualize CRISPR-Millipede results at nucleotide resolution, we developed 'board plots' that display both posterior inclusion probabilities (PIPs, shown as square size) and effect sizes (shown as color intensity) for each nucleotide position, enabling rapid identification of functional variants.

Strong correlations across replicates for allele read counts in CD19 positive and CD19 negative populations (Supplementary Fig. 4a-d) and differential editing frequencies (CD19 negative minus CD19 positive) per variant (Supplementary Fig. 5a-b) demonstrated the robustness of this approach. CRISPR-Millipede effectively highlighted clustered substitutions with high and reproducible penetrance on CD19 expression at positions 151 and 223, with high editing frequency providing power to detect effects across the tiled regions. Sites identified by both approaches were predominantly located within highly conserved sequences containing binding motifs for transcription factors (TFs) abundantly expressed in NALM6 cells, such as MYB, EBF1, SPIB, IKZF1, and PAX5 (Fig. 2c). CRISPR-Millipede also revealed four regulatory sites with a high PIP (>0.4) which were not detected using the count-based screening method. The ABE8e-SpRY screen identified hits at positions 137, 176, and 311, as well as positions 200-201 identified in the evoCDA screen (Fig. 2b). Among the top 6 CRISPR-Millipede hits in the ABE8e-SpRY screen, position 223 showed consistent PIPs and beta coefficients between individual replicates and the within-editor joint model (Supplementary Fig. 6a). In contrast, the remaining 5 hits displayed discordant individual replicate scores yet were still identified as significant by the within-editor joint CRISPR-Millipede model. This highlights CRISPR-Millipede's ability to leverage cross-replicate information to identify true regulatory variants that individual replicate analyses might miss due to experimental variability, noise, and low read depth. A similar pattern was also observed for the evoCDA PIPs and betas across replicates (Supplementary Fig. 6b).

Furthermore, to take advantage of the diversity of edits installed by evoCDA, we compared CRISPR-Millipede analyses to all edits or restricted the analysis to C>T substitutions (Supplementary Fig. 7). The unrestricted analysis recovered not only significant C>T substitutions but discovered additional variants. Notably, this analysis also detected single nucleotide deletion events that corresponded with certain ABE8e-SpRY variant hits, such as at position 137 and 223, suggesting convergent functional effects across different editing modalities (Supplementary Fig. 7b).

To leverage the endogenous allelic readout, we also adapted the DESeq2 framework20 to study allele-level enrichment and significance. This approach allowed us to quantify the differential abundance of specific edited alleles between the CD19 positive and negative populations, including the combined effect of multiple substitutions in cis. DESeq2 identified numerous single-variant and composite alleles that overlapped with CRISPR-Millipede results across ABE8e-SpRY and evoCDA screens (Fig. 2d). Despite low LFC per-replicate correlation among the DESeq2 hits (r=-0.12 to 0.25), the per-replicate comparisons of LFC values showed consistent directional agreement, with most variants falling in the same quadrant across replicates in both ABE8e-SpRY and evoCDA screens (Supplementary Fig. 6c-d).

Additionally, DESeq2 discovered unique hits not detected by CRISPR-Millipede, while CRISPR-Millipede also detected additional unique hits in the ABE8e-SpRY screen, such as at position 176, due to its ability to “fine-map” composite alleles and pinpoint which specific nucleotides confer functionality (Supplementary Fig. 8a-b). To explore the DEseq2-only hits, we examined two variants, 276T>C and 172G>A (Supplementary Fig. 8c-d). Both variants were significant individually, but none of the composite alleles containing them reached significance, suggesting that some of these signals may represent false positives.

To test CRISPR-Millipede's fine-mapping capabilities, we performed a simulation experiment. We excluded all single-variant alleles from the input encoding matrix and re-ran CRISPR-Millipede to determine whether it could still identify functional nucleotides using only composite alleles containing multiple edits. Using our board plots visualizations, the previously identified hits persisted with comparable PIPs and effect sizes, confirming CRISPR-Millipede's ability to extract functional signals from multi-edit composite alleles (Supplementary Fig. 8e-h). For instance, variant 137T>C in the ABE8e-SpRY screen, a single-variant hit in DESeq2 (p-adj = 1.4e-9), persists as a CRISPR-Millipede hit even after removal of the single-variant “137T>C” allele from the input encoding matrix. This is due to composite alleles such as “131T>C,137T>C” being called as hits (p-adj = 1.9e-5), while composite alleles with the bystander variant 131T>C are not called as functional such as “131T>C,130T>C” (p-adj = 0.65) and “131T>C” (p-adj = 0.32) (Supplementary Fig. 8e,g). Additionally, variant 154G>A in the evoCDA screen, a single-variant hit in DESeq2 (p-adj = 3.9e-3), persists as a CRISPR-Millipede hit after removal of the single-variant “154G>A” allele from the input encoding matrix, since composite alleles such as “154G>A,155G>A”, “154G>A,155G>A,161G>A,162G>A”, and “154G>A,161G>A,162G>A” are also called as significant (p-adj = 4.9e-7 and p-adj = 4.2e-14). Alleles containing only the bystander variants such as “155G>A” or “155G>A,161G>A,162G>A” are not called as functional (p-adj = 0.21 and p-adj = 0.24) despite having differential editing frequencies, highlighting CRISPR-Millipede’s ability to pinpoint the true functional variant (Supplementary Fig. 8f,h). These results demonstrate how CRISPR-Millipede can call variant hits with low base-editing purity due to deconvoluting between composite alleles.

Importantly, to demonstrate that CRISPR-Millipede identifies functional variants independent of editing frequency, and to test whether increased read depth improves sensitivity at regions with low base-editing purity, we performed an orthogonal validation experiment using a single guide (sg218) targeting the position 218-235 region. We transduced NALM6 cells expressing ABE8e-SpRY with sg218 alone, followed by sorting, amplicon-sequencing, and CRISPR-Millipede analysis. While the pooled ABE8e-SpRY screen identified 223A>G and 230A>G as variant hits, the sg218-only experiment additionally identified 220A>G and 224A>G (Supplementary Fig. 9a,b). These CRISPR-Millipede results are consistent with the per-replicate single-variant allele DESeq2 LFC scores (Supplementary Fig. 9c). Notably, despite the high editing frequency of 228A>G, CRISPR-Millipede did not call this variant as a hit, while 230A>G was detected as a hit despite low editing frequency, demonstrating that increased read depth enables better deconvolution between bystanders. These results confirm that increasing read depth helps CRISPR-Millipede identify additional potentially causal variants and that CRISPR-Millipede scores are driven by functional effects rather than editing frequency alone.

Validation of gRNA-based and allele-based hits
Following our computational analysis, we validated gRNA-based and allele-based hits in vitro. We transduced NALM6 cells expressing different base editors, including ABE8e-SpRY, evoCDA and the cytosine base editor CBE-SpRY, with individual gRNAs in an arrayed fashion and measured CD19 expression 7 days post transduction via flow cytometry. To link alleles to function, we used gRNAs from the initial screening library as well as newly designed orthogonal gRNAs (Fig. 3a, Supplementary Fig. 10a). Importantly, the PAMless CBE-SpRY editor enabled precise installation of specific alleles and targeting of individual cytosines. For normalization, we aimed for a 50% gRNA transduction efficiency and compared the CD19 levels of gRNA-transduced cells (GFP+) versus non-transduced cells (GFP-) within the same well (Supplementary Fig. 10b).

Using this approach, we tested the two hit regions that overlapped between the gRNA- and allele-based analyses, hits uniquely recovered using the allele-based framework, and the region identified by the gRNA-based analysis only. By targeting nucleotides at positions 145-166 and 218-235, we validated the two hits that were identified by both gRNA- and allele-based analysis (Fig. 3b-c, Supplementary Fig. 10c-d). Next, we focused on variants exclusively identified by the allele-based analysis, at positions 137, 176, 200-201 and 311. By editing position 137 in NALM6 ABE8e-SpRY cells, we observed that a single base substitution at this position using different gRNAs reduces CD19 expression by 17-44%, compared to non-transduced cells, (Fig. 3d-e), while three neutral gRNAs, designed to target non-hit regions, showed no effect on CD19 expression. Similarly, editing of position 176 led to a 23% reduction of CD19 levels (Fig. 3f, Supplementary Fig. 10e). Analogously, installing indels at positions 200-201 using Cas9 and targeting position 311 with ABE8e-SpRY further validated the contribution of these regions to CD19 expression (Fig. 3g-h, Supplementary Fig. 10f-g). These results demonstrate that CRISPR-Millipede uncovers truly functional nucleotides within enhancers.

To investigate the hit region at position 326-330, identified by the gRNA-based analysis only (Fig. 2a), we sought to determine whether specific nucleotide substitutions or rather the steric hindrance of the base editor localizing at the hit region led to decreased CD19 expression. For this, we developed a destabilizing domain23 on the N-terminus of our ABE8e-SpRY editor (ddABE8e-SpRY), which ensures the continuous targeted degradation of the editor. Culturing ddABE8e-SpRY expressing cells in the presence of Shield-1 prevents proteasome-induced editor degradation, allowing for editing, while ddABE8e-SpRY is rapidly degraded in the absence of Shield-1 (Supplementary Fig. 10h). We then performed validation experiments by transducing ddABE8e-SpRY cells with gRNAs targeting position 326-330 as well as the shared hit regions 145-166 and 218-235 and culturing them in Shield-1 for 7 days to ensure editing. Following, we cultured the cells in the presence and absence of Shield-1 for 3 additional days before measuring CD19 levels by flow cytometry (Fig. 3i). This approach distinguished the contribution of single nucleotide substitutions from base editor steric hindrance effects by directly comparing CD19 levels in cells subjected to a Shield-1 wash-out and cells continuously grown in Shield-1. Strikingly, the reduction in CD19 levels upon targeting position 326-330 was rescued after Shield-1 withdrawal (Fig. 3j, Supplementary Fig. 10i). In contrast, reduction in CD19 levels persisted upon targeting positions 145-166 (sg 145) and 218-235 (sg 218) regardless of Shield-1 presence (Fig. 3j). These results indicate that while specific mutations at positions 145-166 and 218-235 impair CD19 enhancer function, the effect at position 326-330 results from base editor occupancy at sgRNAs 322 and 325, rather than nucleotide substitutions. Use of such a degradable base editor hence allows to distinguish motif-disrupting edits from sequence-independent TF displacement due to Cas9 steric hindrance.

CRISPR-Millipede boards to identify CD19 trans-regulators
As non-coding variants frequently disrupt TF motifs24,25 and thereby affect regulatory element activity, we investigated whether validated CRISPR-Millipede hits at CD19 colocalize with relevant TF binding sites. To link CRISPR-Millipede variant hits to potential causal TF motifs, we extended our board plot visualizations by plotting recovered effect sizes as logo-style tracks next to position weight matrix (PWM) motif logos (Fig. 4a, Supplementary Fig. 11a). This enables assessing the potential impact of each observed substitution on the binding affinity of each factor. For instance, the logo tracks reveal that observed significant substitutions often disrupt critical nucleotides for TFs like IKZF1, SPIB, PAX5, and EBF1 (Fig. 4a, Supplementary Fig. 11a), providing insight into the potential mechanism of action for these variants. Notably, the board plot shows that 154G>A affects the motifs of MYB, ZBTB18 and ZEB1, while 151A>G only disrupts the MYB motif (Fig. 4a).

To further validate that disruption of the MYB and PAX5 motifs affects CD19 expression, we derived single cell clones of ddABE8e-SpRY cells homozygous for 150A>G and 151A>G mutations or for 228A>G and 230A>G substitutions, respectively. We assessed CD19 expression of these clones upon washing-out Shield-1 and culturing them in normal media for 6 days, thus avoiding steric hindrance effects. In line with our allele-based analysis, motif-mutated cells display lower CD19 levels compared to AAVS1 transduced cells (Supplementary Fig. 11b). ChIP-seq data of PAX5 in NALM6 cells16 further confirmed that this factor localizes at the CD19 enhancer (Supplementary Fig. 12a), suggesting its CD19 regulatory role through enhancer binding. Similarly, IKZF1 and EBF1 were also found to occupy this region (Supplementary Fig. 12a), suggesting that several TFs, whose motifs are enriched in our hit regions, also localize at the CD19 enhancer.

To systematically identify TFs whose motifs are not only disrupted by single nucleotide substitutions in our identified hit regions, but also regulate CD19, we conducted a KO screen. We transduced Cas9-expressing NALM6 cells with a gRNA library targeting 635 TFs of the hematopoietic lineage and chromatin regulators. The library also included 100 GFP-targeting and 100 non-targeting gRNAs, which were used as negative controls. Following transduction, we cultured the cells for 7 days before sorting for the 5% highest and lowest CD19-expressing gRNA-transduced cells (Supplementary Fig. 12b).

We observed that TFs whose motifs are enriched in our hit regions (Fig. 4a, Supplementary Fig. 10a), such as EBF1, PAX5 and MYB, are among the top CD19 regulators (Fig. 4b), with MYB being the only essential TF (Supplementary Fig. 12c). KO of PAX5 and EBF1 in Cas9-expressing cells resulted in a 64% and 37% reduction in CD19 expression, respectively, further validating our findings (Supplementary Fig. 12d). While publicly available ChIP-seq data14 showed IKZF1 binding at the CD19 promoter and enhancer (Supplementary Fig. 12a), depletion of IKZF1 did not affect CD19 expression (Supplementary Fig. 12e). To confirm these results, we transduced IKZF1-GFP reporter NALM6 cells with gRNAs targeting IKZF1 or treated them with Lenalidomide, a potent IKZF1/3 degrader26. In both experiments, loss of GFP indicated effective inactivation of IKZF1, while CD19 levels stayed consistent (Supplementary Fig. 12e-f).

Interestingly, MYB KO led to an increase in CD19 levels (Fig. 4b), while editing of the MYB motif at the CD19 enhancer decreased CD19 expression (Fig. 4a, S12c), suggesting both a direct and indirect roles for MYB in regulating CD19. To investigate this further, we combined MYB motif editing using sg145 with a 24 h Mebendazole treatment (Fig. 4c). Using a MYB PROTAC allowed us to rapidly degrade MYB protein, circumventing the fact that MYB is essential in NALM6 cells (Supplementary Fig. 12g). When treating cells only with Mebendazole, we observed an increase in CD19 expression, while motif-edited cells showed decreased CD19 levels (Fig. 4c). Upon degrading MYB in motif edited cells, we observed a rescue in CD19 expression, which was most prominent in cells expressing very low CD19 levels (Fig. 4c – black box). Together, these results propose a dual role of MYB, where editing of the motif results in CD19 downregulation, pointing to MYB as a positive regulator, and KO or degradation of MYB leads to CD19 upregulation, suggesting that MYB also functions as a negative CD19 regulator, possibly via an indirect mechanism.

Nucleotide variants affect response to aCD19 CAR-T therapy
To demonstrate the utility of our framework in uncovering clinically relevant nucleotide variants, we investigated the impact of enhancer mutations on CAR-T cell therapy efficacy, a critical issue in cancer treatment. CD19 is a target of CAR-T cell therapy in B-cell malignancies. Despite initial high success rates, ~35% of B-ALL cases treated with aCD19 CAR-T relapse27, with a subset due to loss or downregulation of the CD19 antigen28. To test how substitutions in the CD19 enhancer impact physiological expression levels in the context of aCD19 CAR-T, we performed a competition experiment with unedited NALM6 cells (marked with BFP) and genome-edited NALM6 cells (marked with GFP), co-cultured with aCD19 CAR-T or mock T-cells (Fig. 5a). After 3 days of co-culture, we observed that targeting the CD19 exon 2 splice acceptor site (our positive control) resulted in accumulation of cells resistant to CAR-T-mediated killing, in line with known variants found at the CD19 splice site in CAR-T relapses. Targeting a neutral sequence of the CD19 enhancer with sg119 did not confer any resistance advantage. Using CAR-T cells derived from independent donors, we observed substantial outgrowth of cells that were resistant to aCD19 CAR-T after editing with sg 145 and sg 218 (Fig. 5b-c, Supplementary Fig. 13a-d). These finding demonstrate that genetic alterations at the CD19 enhancer can drive resistance to CD19-targeted immunotherapies. Notably, the enhancer variants we identified reduced CD19 expression to levels sufficient to evade CAR-T recognition while potentially maintaining B-cell viability, representing a clinically relevant mechanism of antigen escape distinct from complete CD19 loss. These results highlight the importance of surveying non-coding regions when investigating mechanisms of immunotherapy resistance.

Methods
Cell culture and cell lines
NALM6 cells (CVCL_0092) were cultured in RPMI 1640 (Thermo Fisher Scientific, # A41923-01) with 10% Fetal Bovine Serum (Fisher Scientific, #10500064), 100 U/ml Penicillin-Streptomycin (Fisher Scientific, #11548876) and 2 mM L-Glutamine (Thermo Fisher Scientific, #11510626) at 37 °C, 5% CO2. HEK293T cells were cultured in DMEM (Thermo Fisher Scientific, #11584486) with 10% Fetal Bovine Serum (Fisher Scientific, #10500064), 100 U/ml Penicillin-Streptomycin (Fisher Scientific, #11548876) and 2 mM L- Glutamine (Thermo Fisher Scientific, #11510626). To achieve stable expression of NGG-evoCDA, ABE8e-SpRY, ddABE8e-SpRY, CBE-SpRY, and Cas9 (Supplementary Data 1) 5x105 cells were seeded in 1 ml RPMI in 24 well plates, transduced with lentiviruses, and selected with 10 µg/ml Blasticidin (Gibco, # A1113903) for five days. Expression of base editors was confirmed by Western Blot using a Cas9 antibody (Thermo Fisher Scientific #MA1-202, dilution 1:1000).

Healthy donor buffy coats were purchased from the Austrian Red Cross, primary human T cells were purified (RosetteSep Human T cell Enrichment kit; STEMCELL Technologies, #15021) and cryopreserved until further use. To generate CAR-expressing T cells, primary human T cells were thawed and activated (Dynabeads Human T-Activator αCD3/αCD28 beads; Thermo Scientific, #11131D), followed by transduction with purified lentiviruses (Lenti-X Concentrator; Takara, #631232) and Puromycin (Gibco, # A1113803) selection29. T cells were cultured at a density of 0.3 - 2x106 cells/ml in AIMV medium (Thermo Fisher Scientific, #12055091) supplemented with 2% Octaplas (Blutspendezentrale Wien), 2.5% HEPES (Thermo Fisher Scientific, #15630056), 1% L-glutamine (Thermo Fisher Scientific, #11510626) and 200 U/mL recombinant human IL-2 (Peprotech, #200-02-250UG).

Generation of genome-edited NALM6 clones
1x106 NALM6 Cas9 cells were electroporated with two GFP-labeled gRNAs targeting sequences flanking the CD19 enhancer using a NEPA21 electroporator. Three days post-electroporation, GFP+ cells were sorted and plated at a density of 5 cells/well in 200 µl in a 96 well plate. Emerging clones were genotyped by PCR. Sequences of gRNAs and genotyping oligos are listed in Supplementary Data 2.

RNA isolation and RT-qPCR
RNA was isolated from cells using RNeasy Mini Kit (Qiagen, #74134), and cDNA was prepared using iScript cDNA Synthesis Kit (Bio-Rad, #1708891). RT-qPCR was performed using iQ SYBR Green Supermix (Bio-Rad, #1708885) on Bio-Rad iCycler RT-PCR detection system. Primers are listed in Supplementary Data 3. Each sample was normalized to the housekeeping gene GAPDH and compared using the ddCt method.

Flow cytometry and cell sorting
CD19 expression was measured using a conjugated antibody (APC Clone SJ25C1; BD Biosciences; #340722). Briefly, 1-2x105 cells were incubated with a 1:25 antibody dilution in DPBS for 25 minutes at 4°C. APC Mouse IgG1 κ Isotype Control (BD Biosciences; #554681) was used as control. After staining, cells were washed twice, resuspended in DPBS and analyzed using a BD LSRFortessa™ Cell Analyzer. Biological triplicates were measured, and experiments were repeated at least twice. CD19 levels were measured in 6-8 different CD19 eWT, CD19 e+/-, and CD19 e-/- clones.

For sorting of the base editor screens, 20-50x106 cells per replicate were washed in DPBS and stained using the APC-conjugated antibody (APC Clone SJ25C1; BD Biosciences; #340722) from above in 200-300 µl staining mix. For the TF KO screen, CD19 expression was measured using a FITC-conjugated antibody (APC Clone SJ25C1; BD Biosciences; #345788) and FITC Mouse IgG1 κ Isotype Control (BD Biosciences; #556649) was used as a control. 130x106 cells were stained, washed, resuspended in DPBS and strained through a 35 µm nylon mesh (Fisher Scientific, #10585801). FACS was performed with BD FACSAria™ Cell Sorter system and sorted cells were collected in 0.5-1 ml DPBS.

Cloning of gRNAs and tiling library
gRNAs were cloned by Golden Gate Cloning in pLKO5.sgRNA.EFS.GFP plasmid (Addgene #57822) and verified by Sanger Sequencing. gRNA sequences are listed in Supplementary Data 2. The coding sequence of evoCDA-be4b (gift of David Liu) was cloned by Gibson Assembly into XbaI-BamHI-digested LentiCas9-Blast (Addgene #52962). The sequence of ABE8e-SpRY was cloned in a lentivector under the SFFV promoter. The ddABE8e-SpRY construct was cloned by amplifying the destabilizing domain (DD) sequence from DD-Cas9-venus (Addgene #90085) and performing Gibson Assembly into AgeI-EcoRI-digested ABE8e-SpRY. For CBE-SpRY (TdCBEd) cloning, a mutated TadA30 (E27K, V28A, M61I, H96N) was fused to the N-terminus of D10A SpRY Cas9, while a P2A-2xUGI-P2A-Blasticidin cassette was fused to its C-terminus. This large DNA cassette was driven by the SFFV promoter and cloned into a lentiviral vector.

For the library, all possible 20-mer sequences tiling the top and bottom strand of the CD19 enhancer were ordered as oligonucleotide pool (IDT), PCR-amplified and cloned by Gibson Assembly in Esp3I-digested pLKO5.sgRNA.EFS.GFP (Addgene #57822). Library gRNAs and oligonucleotides used for cloning are listed in Supplementary Data 4-5, respectively. Quality and completeness of the gRNA library was assessed by NGS.

Lentivirus production
HEK293T cells were seeded in 15 cm dishes 24 h prior to transfection. Cells were transfected at 80% confluency in 16 ml of media with 8.75 μg of VSVG, 16.25 μg of psPAX2, and 25 μg of the lentiviral vector, using 150 μg of linear polyethylenimine (PEI) (Sigma-Aldrich, #408727-100ML). Media was changed to fresh media 16–24 h post-transfection. Lentiviral supernatant was collected at 48 h and 72 h post-transfection and subsequently concentrated by ultracentrifugation (24000 rpm, 4 ºC, 2 h) in a 20% sucrose gradient.

Base editor screen in NALM6 cells
For base editor screens, 5–7.5x106 NALM6 cells stably expressing evoCDA or ABE8e-SpRY were seeded in 10-15 ml RPMI and transduced at a 0.3 MOI to ensure single viral integrations per cell. A titration experiment was performed to determine the amount of virus required to achieve a transduction rate of 30%. A 1,000x gRNA representation was maintained throughout the screen. Seven days post transduction, 20-50x106 cells were collected, counted, and processed for FACS sorting. Cells with lentiviral library integrations were selected by gating for GFP. For the screen conducted in evoCDA expressing cells, 1.5-3.5x105 CD19 low and 3x106 CD19 high cells were sorted per biological replicate. For the screen conducted in ABE8e-SpRY expressing cells, 5.5-6.7x105 CD19 low cells and 1.8-2.3x106 CD19 high cells were sorted per biological replicate. Genomic DNA was isolated with DNeasy Blood & Tissue Kit (Qiagen, #69506) and next-generation sequencing (NGS) libraries for the count-based read-out were prepared as previously reported31. Distributions of edits within a 346 bp fragment of the CD19 enhancer were obtained by targeted amplicon sequencing using MiSeq at PE250 and PE300 configuration.

Oligonucleotides used for count- and allele-based library preparation are listed in Supplementary Data 5-6, respectively. Screens were performed in biological triplicates.

For the validation of sg 218, cells were treated equally to the base editor screen but transduced only with sg 218. 6.7x105-1.2x106 CD19 low cells and 7.8x105-1.2x106 CD19 high cells were sorted per biological replicate. Three biological replicates were used for amplicon sequencing and CRISPR-Millipede analysis.

TF KO screen in NALM6 cells
A gRNA library targeting 635 TFs and chromatin regulator genes, as well as 100 non-targeting and 100 GFP-targeting gRNAs was designed (Supplementary Data 7). Validated gRNAs from the Avana library32 were used for the library design, and if not available, gRNAs for few specific genes were designed using CRISPOR33. For the TF KO screen, 20x106 clonal NALM6 Cas9 cells were transduced at a 0.3 MOI to ensure single viral integrations per cell. A titration experiment was performed beforehand to determine the virus amount necessary to achieve a 30% transduction and a 1000x gRNA representation was maintained throughout the screen. Seven days post transduction, 130x106 cells were collected, counted, stained and processed for FACS sorting. Cells with lentiviral library integrations were selected by gating for iRFP. 5% of cells with highest (3.9-4.7x105 cells) and 5% of cells with lowest (4.7-8.4x105 cells) CD19 levels were sorted. Genomic DNA was isolated with DNeasy Blood & Tissue Kit (Qiagen, #69506) and next-generation sequencing (NGS) libraries for the count-based read-out were prepared as previously reported31. The screen was performed in biological triplicates.

Validation experiments using individual sgRNAs
5x104 NALM6 cells stably expressing evoCDA, ABE8e-SpRY, ddABE8e-SpRY, CBE-SpRY, or Cas9 were seeded in 100 µl RPMI per well in 96 well plates and transduced with lentiviruses. Seven days post transduction, 1-2x105 cells were used for flow cytometry (see above) and the rest was collected for Sanger sequencing.

For experiments using the ddABE8e-SpRY base editor system, NALM6 cells were seeded and transduced in media containing 200 nM of Shield-1 (TargetMol, #T13884). Seven days post transduction, cells were washed once in DPBS to remove Shield-1 (TargetMol, #T13884). After 3 days, 1-2x105 cells were used for flow cytometry (see above) and the rest was collected for Sanger sequencing.

For all validation experiments, data was normalized by calculating the ratio of CD19 MFI in transduced (GFP positive) vs. non-transduced (GFP negative) cells. PCR and Sanger oligonucleotides are listed in Supplementary Data 2.

IKZF1 validation experiment
5x104 NALM6 Cas9 IKZF1-GFP reporter cells were seeded in 100 µl in a 96 well plate and transduced with gRNAs targeting AAVS1 or IKZF1 itself. Transduced cells were selected using Puromycin (Gibco, #A1113803) for 7 days and later GFP levels and CD19 levels were assessed via flow cytometry.

Treatment with IKZF1 and MYB degraders
5x104 NALM6 Cas9 IKZF1-GFP cells in 100 µl in 96 well plates were treated for 48 h with 2 µM Lenalidomide (Merck, #901558) and 1x105 NALM6 ABe8e-SpRY cells in 100 µl in 96 well plates were treated for 24 h with 10 µM Mebendazole (MedChemExpress, #HY-17595) prior to assessment of IKZF1 and/or CD19 levels by flow cytometry. Cells treated with DMSO (Fisher Scientific, # BP231-100) were used as controls.

Western Blot
2x106 cells were lysed using NP-40 buffer (50 mM Tris-Hcl (pH 8.0), 150 mM NaCl, 1% NP-40) supplemented with 1:100 Protease Inhibitor Cocktail (Sigma, #P8340) and 1:200 Benzonase (Merck, #70664). Following a 15 min incubation at 4 °C, samples were centrifuged at 15,000 g for 10 min. The supernatant containing the protein lysates was transferred to a new tube. Protein quantification was conducted using Pierce BCA Protein Assay Kit (Thermo Scientific, #23227) according to manufacturer’s instructions and 70 µg of protein was used for Western Blot. Proteins were denatured using XT Sample Buffer (BioRad, #1610791) for 5 min at 95 °C and separated in 4-12% Criterion XT Bis-Tris Protein Gels (BioRad, #3450125). Proteins were then transferred to a low-fluorescence PVDF membrane (BioRad, #1704275) using the High MW program of the Trans-Blot Turbo Transfer System (BioRad, #1704150). The membranes were blocked for 1 h, incubated overnight at 4 °C in primary antibody (1: 1,000 MYB – Proteintech, 17800-1-AP; 1: 20,000 b-Actin – Santa Cruz, sc-47778) and at RT for 2 h in secondary antibody (1: 10,000 anti-rabbit – Cell Signaling, #7074) in EveryBlot blocking buffer (BioRad, #12010020). The membranes were developed using the Clarity Max™ Western ECL Substrate (BioRad, #1705062) and a BioRad ChemiDoc imager.

Quantification of CD19 molecules per cell
For the quantification of CD19 molecules per cell, the BD Quantibrite™ Beads PE Fluorescence Quantitation Kit (BD Biosciences #340495) was used and CD19 expression in cells was assessed using a PE-conjugated antibody (PE anti-human CD19; clone HIB19; Biolegend, #302208). Briefly, 1-2x105 cells were washed in FACS buffer (DPBS containing 1% Albunorm (200 g/l; Octapharma, #14512) and 0.2% NaN3), incubated in blocking solution (10% human serum (Sigma, #H4522-20ML) in FACS buffer) for 10 min at 4 °C, and stained with a 1:100 antibody dilution for 25 min at 4 °C. PE Mouse IgG1 κ Isotype Control (clone MOPC-21; Biolegend, #400111) was used as control. After two more rounds of washing with FACS buffer, cells were analyzed using a BD LSRFortessa™ Cell Analyzer. Molecule per cell quantifications were done via linear regression according to manufacturer’s instructions. Technical triplicates were measured.

NALM6 T cell co-culture experiments
5x105 NALM6 ABE8e-SpRY cells were seeded in 1 ml RPMI in 24 well plates, transduced with concentrated lentiviruses expressing BFP alone or GFP in combination with a sgRNA, and sorted based on BFP or GFP expression 72 h post transduction. For co-culture experiments, 2.5x104 BFP and 2.5x104 GFP-positive cells were mixed with either mock (non-transduced) or aCD19-expressing CART cells in 96 well plates in an effector-target-ratio of 0.2:1 and cultured in a total of 200 µl RPMI media containing 200 U/ml recombinant human IL-2 (PeproTech, #200-02-250UG), in U-shaped 96 well plates for 3 days. For flow cytometry, cells were washed with FACS buffer and measured using a BD LSRFortessa™ Cell Analyzer. Eight to ten biological replicates were measured and the experiment was performed four times using aCD19 CAR-T cells from four different donors.

Sequencing and demultiplexing
Amplicon-sequencing of endogenous allele samples was conducted on Illumina MiSeq using paired-end 250 and 300 bp and 30% PhiX. Initial demultiplexing for these samples was carried out by the Biomedical Sequencing Facility (BSF) at the Research Center for Molecular Medicine (CeMM). The count samples were sequenced on a NovaSeq SP system, utilizing half of a flow cell with paired-end 100 bp reads. The sequencing core facility performed the primary demultiplexing based on the P7 index. Subsequently, our laboratory conducted a secondary demultiplexing step using the P5 index, which is located internally within the read.

TF KO screen analysis
MAGeCK MLE v0.5.9.534 was used for read mapping as well as comparing plasmid vs presort and low vs high CD19 sorted populations. For the analysis, 100 non-targeting gRNAs were used as negative controls. Genes were ranked based on descending beta scores and significant hits (FDR ≤ 0.05) were plotted in green (positive regulators) and red (negative regulators).

gRNA Mapping
For each demultiplexed FASTQ, “raw” gRNA mapping was performed using CRISPR-SURF Count (docker version pinellolab/crisprsurf:crispr_clear_v1) in base-editing mode. To address possible mismatches in the gRNA cassette, which may include sequencing errors, and self-editing in the ABE8e-SpRY samples, we developed CRISPR-Correct (version 0.0.41). The latter maps the observed gRNA sequence to an inferred gRNA sequence that possesses the lowest sequence hamming distance, therefore allowing for possible mismatches between the observed and mapped sequence. CRISPR-Correct mapping was performed using the following function:
crispr_ambiguous_mapping.get_guide_counts_from_fastq(GUIDE_SERIES,INPUT_FASTQ_FN,hamming_threshold_strict=11,parse_left_flank=False,parse_flank_sequence="GTTTTA", cores=50)

CRISPR-SURF analysis of gRNA counts
To identify significant regions based on the mapped gRNA counts, we ran CRISPR-SURF Deconvolution on both raw and mismatch-corrected ABE8e-SpRY and evoCDA datasets. CRISPR-SURF Deconvolution was run in base-editing mode with a lambda value of 1. For the evoCDA data, gRNAs with an -NGG PAM and a cytosine within the expected editing window (positions 1 to 13 of the protospacer) were designated as observation gRNAs, while gRNAs not meeting these criteria were designated as negative control gRNAs. The perturbation range for evoCDA was set to 13 based on available data on evoCDA editing activity across six HEK293T sites15. For the ABE8e-SpRY data, gRNAs that lack an adenine in the expected editing window (positions 3 to 8) were designated as negative controls, and the perturbation range was set to 6, consistent with the known perturbation window for ABE8e35.

CRISPR-Millipede analyses of endogenous alleles
The CRISPR-Millipede pipeline consists of five steps taking the demultiplexed FASTQs from amplicon-sequencing of the endogenous alleles as input.

1. CRISPResso2
After targeted amplicon-sequencing of the endogenous alleles, the paired-end FASTQs were quality controlled with FastQC36 version 0.12.0 and MultiQC37 version 1.15 to ensure sufficient sequencing quality of our samples. Next, the paired-end FASTQs were passed into CRISPResso2 version 2.1.3 for read merging, alignment, and quantification of the edited alleles by running the following command:
CRISPRessoBatch -bs {FASTQ_FILENAME} -a {AMPLICON_SEQUENCE} -an cd19 -q 30 --exclude_bp_from_left 3 --exclude_bp_from_right 3 --no_rerun -n {SCREEN_NAME} --min_frequency_alleles_around_cut_to_plot 0.001 --max_rows_alleles_around_cut_to_plot 500 -p 20 --plot_window_size 4 --base_editor_output -w 0 -bo {OUTPUT_DIRECTORY}

2. Processing and encoding of alleles
To identify functional variants from the CRISPResso2 outputs, we developed a Python package called CRISPR-Millipede (version 0.1.89) to process the CRISPResso2 allele tables and perform statistical modeling. Specifically, the CRISPResso2 allele frequency tables were encoded into a feature indicator matrix containing columns for all possible REF>ALT variants across the amplicon and rows for all unique alleles. The first and last 20 positions of the amplicon were removed from the encoding matrix due to elevated sequencing error background. Only canonical REF>ALT variants were included in the encoding matrix depending on the editor used (A>G and T>C for ABE8e samples, C>T and G>A for evoCDA samples). Since evoCDA is only expected to edit for protospacers with an -NGG PAM, only variants within the canonical editing window (protospacer position 6±7) of the -NGG PAM guide RNAs were included in the encoding matrix.

3. Per-variant modeling of encoded alleles
The encoding matrices were then fed as input into millipede21 version 0.1.2 (PyPi package millipede-regression), a PyTorch-based package for Bayesian variable selection, in order to infer beta coefficients and posterior inclusion probabilities for each feature. Encoded alleles that do not possess at least two replicates with a read count greater than 3 in at least one population (presort, CD19 positive, or CD19 negative) were filtered out to reduce potential noise. In general, replicates within the population had higher correlation of allelic read counts than across populations (Supplementary Fig. 4 and S5c-d). The model was fit using NormalLikelihoodVariableSelector with the (continuous) response variable set to Y = (xn,pos - xn,neg) / (xn,pos + xn,neg) where xn,pos and xn,neg are the normalized read count in the CD19 positive and CD19 negative populations for each encoded allele n, respectively.
Normalization was performed by dividing the raw read count by the total read count for each sample. To account for the higher standard deviation in the response variable expected for alleles with lower read counts (Supplementary Fig. 14a), we introduced a scaling factor, called sigma scale, which adjusts for the increased uncertainty associated with lower read counts. The sigma scale factor is s = (decay(xn,pos)+decay(xn,neg))/2 where decay(x) = 1+(K-1)e^(-ln(0.01)/a * x), therefore each encoded allele n will possess a different relative standard deviation depending on xn,pos, xn,neg, a (representing the normalized read count where the response noise stabilizes), and K (representing the level of response noise for low read counts vs. a read counts). For the ABE8e and evoCDA screen samples, K is set to 5 and a is set to 0.0005. Alleles in both ABE8e and evoCDA screens exhibited varying scale factors (Supplementary Fig. 14b). For the ABE8e sg218-only sample, K is set to 10 and a is set to 0.0001. This adjustment is made to account for the difference at which the standard deviation approximately stabilizes for the sg218 sample. Replicates are modeled jointly by including multiple intercept variables indicating which replicate an encoded allele/response pair belongs to. Millipede infers posterior parameters under the assumption that Y ~Normal(β0 + Xβ, (s * σ)^2) where the β vector of coefficients corresponds to the set of included features after Bayesian variable selection, thereby providing inferred effect sizes and posterior inclusion probabilities for each feature. The inclusion of covariates is governed by a Bernoulli prior. Moreover, the overall noise level σ is inferred from the data where σ^2 is governed by an inverse Gamma prior. The coefficients including the intercept are governed by normal priors. MCMC is used to sample from the posterior distribution, with a total of T = 5000 iterations and a burn-in period of 500 iterations to allow for convergence. For all other parameters default values were used. See the millipede19 documentation or publication for additional details on the Bayesian variable selection model used.

4. DESeq2 analysis of endogenous alleles
To calculate the differential enrichment of each allele in the CD19 positive and CD19 negative populations, the encoding matrix generated by CRISPR-Millipede together with a design matrix describing the condition of each sample (i.e., replicate and sorting population) were passed into DESeq2 (version 1.48.1). Allele counts were normalized for library size using the unedited allele as the sole control gene in the size factor estimation step, so that the log fold change (LFC) is centered such that an LFC of 0 corresponds to the differential enrichment of the unedited allele. Dispersion estimates were obtained across replicates, and differential testing between CD19 positive and CD19 negative populations was performed using the Wald test. Outliers were assessed using Cook’s distance. No log fold-change shrinkage was applied, and results are reported as the unshrunken log2 fold change estimates. P-values were adjusted for multiple testing using the Benjamini–Hochberg procedure. The resulting log2 fold changes and adjusted p-values were used for downstream visualization and interpretation of allele distributions.
For each allele, DESeq2’s baseMean statistic (the mean of normalized counts across all samples) was also recorded and used downstream as a proxy for overall allele abundance.
We also calculated per-replicate log2 fold changes by applying DESeq2’s variance stabilizing transformation to the count matrix and subtracting the transformed values of the CD19-negative samples from their matched CD19-positive samples within each replicate.

5. Plotting
CRISPR-SURF and CRISPR-Millipede plots were generated using jdgenometracks version 0.1.60 (PyPi package jdgenometracks), available at https://github.com/justin-delano/jdgenometracks. Board plots were generated using a custom function in the Millipede package. Logomaker was used to generate PWM plots38.

Motif analysis
All motif PWMs were obtained from the JASPAR database39 (v24) and CIS-BP40. Motifs were scanned using MOODS41.

Statistics and Reproducibility
For the analysis of the base editor screen gRNA counts and the TF KO screen published statistical frameworks were used as referenced in the respective main text and methods sections. Otherwise, the Prism software (version 9.3.1, GraphPad, USA) was used for statistical analyses, and data are shown as mean ± SD. Statistical tests, p-value cut-offs, and sample sizes are indicated in the respective Figure legends. Additional details on statistics can be found in the respective source data files.

Sample size is the number of biological replicates used for each experiment. All experiments were conducted using at least three biological replicates. No statistical method was used to predetermine sample sizes, but the sample size was sufficient to ensure reproducibility across experiments and in previous studies. No data were excluded from the analyses. Experiments were repeated at least twice on independent days using different batches of reagents and aCD19 CAR T cell donors. Randomization was ensured as cells were plated in different plate formats (6, 24 and 96 well plates) and treated in different wells. Investigators were not blinded to allocation during experiments and outcome assessment, since samples had to be identified for analysis.

Schematic figures were generated with Adobe Illustrator (Adobe, USA).
Figure 1 and Supplementary Figure 12 contain elements created with BioRender.com (https://BioRender.com/8u19bsf; https://BioRender.com/q3kyk09).

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The authors have provided a thoughtful reply to comments and addressed all major points. I am satisfied by the revision.", "reviewer_2": "The authors have substantially revised their manuscript. I appreciate that substantial effort has gone into this and think that the transfer of the manuscript to Nature Communications is sensible. While the statistical ingenuity of the manuscript is still burdened by lacklustre data, I appreciate the added transparency and the enhanced explanations of the methods in the text. Therefore, I tentatively recommend publication with a few minor changes if the other referees concur.\n\nHere, I note some of the limitations of the data set from my technical review and suggest some textual changes to provide a more accurate description.\n\nI was concerned about the reproducibility of the original BE screen data. The addition of Figs S3-S6 has made it clear that reproducibly was, indeed, generally low across replicates. This is important information to know, so I’m glad the authors have included this. Some points from these added figures to note:\n\n1. In FigS3 there are generally poor correlations between all replicates, and even a negative correlation between replicates 2 and 3 in (FigS3a). This is despite the fact the authors are only looking at their hits in this figure and not all regions tested.\n2. Again in FigS6, replicate correlations are low (or even negative), and often effects of variants are driven by the variant scoring as a hit in only 1 of 3 replicates (a,b). I thank the authors for their transparency about this. It seems this may be driven by the relatively low number of cells sorted, which is now more transparently reported.\n\nIn light of these major limitations, it is essential to make sure these figures are included in the final form of the manuscript and that transparency is extended to the main text. I do not see how correlations that are sometimes negative can be referred to as “high” (line 170). It would be more accurate to say that the replicate correlations were weak prior to correction and report the range. Please fix this.\n\nI applaud the authors for adding a considerable amount of validation (Fig. 3). In this work, the ddABE8e-SpRY experiments are particularly interesting. I note that two of the initial hits’ effects were likely driven by BE-mediated steric hindrance, and not allele-specific effects. This is an important point to make. It would be worth commenting on this caveat in the Discussion of the manuscript, as this is something that the field needs to better appreciate.", "reviewer_3": "The authors have justified the context of the field and how their work adds. Meanwhile, of course, more work has been published in this space and I suggest the authors carefully discuss the novelty of their approach.\n\nMy impression is that they have done a good job addressing this point and other points raised by the other reviews. I am therefore supportive of publication of this study."}, "author_response": "Response to Reviewer’s comments concerning “CRISPR-CLEAR: Nucleotide-Resolution Mapping of Regulatory Elements via Allelic Readout of Tiled Base Editing” [redacted] revision, transferred to Nature Communications\nWe thank the reviewers for their thoughtful evaluation of our revised manuscript. We are pleased that all three reviewers support publication: Reviewer #1 stated they are \"satisfied by the revision,\" Reviewer #2 \"tentatively recommend[s] publication,\" and Reviewer #3 is \"supportive of publication of this study.\"\nIn this final revision, we have addressed the remaining minor points raised by Reviewer #2 regarding transparency in reporting replicate correlations, and by Reviewer #3 regarding discussion of recent work in the field.\nWe believe the manuscript is now ready for publication and we thank the reviewers and editor for their constructive guidance throughout the review process\nPlease note that answers to the reviewers are marked in blue and changes in the manuscript are marked in green. Addition of new text is indicated as ‘’[NEW TEXT]”.\nReviewer #1\nThe authors have provided a thoughtful reply to comments and addressed all major points. I am satisfied by the revision.\nWe thank the Reviewer for their positive assessment and support for publication\nReviewer #2\nThe authors have substantially revised their manuscript. I appreciate that substantial effort has gone into this and think that the transfer of the manuscript to Nature Communications is sensible. While the statistical ingenuity of the manuscript is still burdened by lacklustre data, I appreciate the added transparency and the enhanced explanations of the methods in the text. Therefore, I tentatively recommend publication with a few minor changes if the other referees concur.\nHere, I note some of the limitations of the data set from my technical review and suggest some textual changes to provide a more accurate description.\nI was concerned about the reproducibility of the original BE screen data. The addition of Figs S3-S6 has made it clear that reproducibly was, indeed, generally low across replicates. This is important information to know, so I’m glad the authors have included this. Some points from these added figures to note:\n1. In FigS3 there are generally poor correlations between all replicates, and even a negative correlation between replicates 2 and 3 in (FigS3a). This is despite the fact the authors are only looking at their hits in this figure and not all regions tested.\nIn light of these major limitations, it is essential to make sure these figures are included in the final form of the manuscript and that transparency is extended to the main text. I do not see how correlations that are sometimes negative can be referred to as “high” (line 170). It would be more accurate to say that the replicate correlations were weak prior to correction and report the range. Please fix this.\nWe thank the reviewer for their careful evaluation of the newly added supplementary figures and for recognizing the increased transparency of the revised manuscript. We agree that reproducibility across replicates in the original base-editing screens was limited, and we appreciate the emphasis on clearly communicating this limitation.\nIn response to the reviewer’s comments, we have ensured that Supplementary Figures 3–6 are included in the final manuscript and revised the corresponding text to more accurately reflect replicate concordance. In particular, we have removed the characterization of replicate correlations as “high” (line 170) and now explicitly state that correlations were weak prior to correction, reporting the full observed range.\nRegarding the point about Supplementary Fig. 3, mismatches arising from self-editing impair gRNA read mapping and contribute to poor replicate agreement. We now explicitly describe this limitation in the main text and emphasize that mismatch-correction with CRISPR-Correct is a necessary preprocessing step prior to downstream analyses.\n[NEW TEXT, Results] When CRISPR-SURF was run without mismatch-correction by CRISPR-Correct, the evoCDA screen yielded only a single hit at position 145-166, and the ABE8e-SpRY screen produced no detectable hits (Supplementary Fig. 3b). Additionally, supporting gRNAs within the CRISPR-SURF sub-regions show low correlation of log fold-change enrichment scores across replicates without CRISPR-Correct (r=-0.08 to 0.45) and improved correlation when using CRISPR-Correct in the ABE8e-SpRY screen (r=0.62 to 0.78) (Supplementary Fig. 3a), highlighting the necessity of correcting for mismatches from self-editing prior to downstream analyses.\n2. Again in FigS6, replicate correlations are low (or even negative), and often effects of variants are driven by the variant scoring as a hit in only 1 of 3 replicates (a,b). I thank the authors for their transparency about this. It seems this may be driven by the relatively low number of cells sorted, which is now more transparently reported.\nRegarding Supplementary Fig. 6, we agree and this is now mentioned explicitly in the text as follow:\n“Among the top 6 CRISPR-Millipede hits in the ABE8e-SpRY screen, position 223 showed consistent PIPs and beta coefficients between individual replicates and the within-editor joint model (SupplementaryFig. 6a). In contrast, the remaining 5 hits displayed discordant individual replicate scores yet were still identified as significant by the within-editor joint CRISPR-Millipede model.”\nWe also want to highlight that to extend transparency for Supplementary Fig. 6c–d, where per-replicate LFC correlations among DESeq2 hits are low and can be negative when subsetting to hits, we now report the correlation range and include a complementary directional-agreement summary in the main text.\n[NEW TEXT, Results]Despite low LFC per-replicate correlation among the DESeq2 hits (r=-0.12 to 0.25), the per-replicate comparisons of LFC values showed consistent directional agreement, with most variants falling in the same quadrant across replicates in both ABE8e-SpRY and evoCDA screens (Supplementary Fig. 6c-d).\nI applaud the authors for adding a considerable amount of validation (Fig. 3). In this work, the ddABE8e-SpRY experiments are particularly interesting. I note that two of the initial hits’ effects were likely driven by BE-mediated steric hindrance, and not allele-specific effects. This is an important point to make. It would be worth commenting on this caveat in the Discussion of the manuscript, as this is something that the field needs to better appreciate.\nWe thank the Reviewer for appreciating our validation experiments. This experiment indeed helped us to disentangle true hits due to motif disruption (e.g. altering PAX5 and MYB binding) versus sgRNAs whose effect is due to steric hindrance (e.g. at sgRNAs 322 and 325); we reinforced this concept in the manuscript by adding remarks in the results and discussion.\n[NEW TEXT, Results] [...] the effect at position 326-330 results from base editor occupancy at sgRNAs 322 and 325, rather than nucleotide substitutions. Use of such a degradable base editor hence allows to distinguish motif-disrupting edits from sequence-independent TF displacement due to Cas9 steric hindrance.\n[NEW TEXT, Discussion] Furthermore, we observed confounding effects due to chromatin binding of base editors, that might displace transcription factors. We addressed this issue by developing a degradable base editor, ddABE8e-SpRY, that allows testing the relevance of individual nucleotides within regulatory sequences without the issue of steric hindrance, a source of false positives.\nReviewer #3\nThe authors have justified the context of the field and how their work adds. Meanwhile, of course, more work has been published in this space and I suggest the authors carefully discuss the novelty of their approach.\nMy impression is that they have done a good job addressing this point and other points raised by the other reviews. I am therefore supportive of publication of this study.\nWe thank the Reviewer for their positive assessment and support for publication. In our previous revision, we expanded the Introduction to position our framework relative to existing approaches (adding five new citations in response to Reviewer 1). We agree that additional studies have emerged in this area during the review process, and have now incorporated two recent works to discuss how our approach complements these methods (Martyn et al. Cell 2025; Martin-Rufino et al. Cell 2023). The following text has also been added to the introduction.\n[NEW TEXT] While prime- and base-editing screens provide high-precision interrogation of individual regulatory variants 12,13, they remain constrained by editing efficiency and a lack of computational methods to accurately identify the causal variants based on mutational outcomes."}]

Decision letter