
==== Front
101573691
39703
Cell Rep
Cell Rep
Cell reports
2211-1247

39163202
10.1016/j.celrep.2024.114640
nihpa2019932
Article
CRISPR screening uncovers a long-range enhancer for ONECUT1 in pancreatic differentiation and links a diabetes risk variant
Kaplan Samuel Joseph 12
Wong Wilfred 13
Yan Jielin 24
Pulecio Julian 2
Cho Hyein S. 2
Li Qianzi 13
Zhao Jiahui 1
Leslie-Iyer Jayanti 2
Kazakov Jonathan 2
Murphy Dylan 1
Luo Renhe 24
Dey Kushal K. 3
Apostolou Effie 5
Leslie Christina S. 3
Huangfu Danwei 26*
1 Weill Cornell Graduate School of Medical Sciences, Weill Cornell Medical College, New York, NY 10065, USA
2 Developmental Biology Program, Sloan Kettering Institute, Memorial Sloan Kettering Cancer Center, New York, NY 10065, USA
3 Computational and Systems Biology Program, Sloan Kettering Institute, Memorial Sloan Kettering Cancer Center, New York, NY 10065, USA
4 Louis V. Gerstner Jr. Graduate School of Biomedical Sciences, Memorial Sloan Kettering Cancer Center, New York, NY 10065, USA
5 Meyer Cancer Center, Division of Neuro-Oncology, Department of Neurology, Sandra and Edward Meyer Cancer Center, New York-Presbyterian Hospital/Weill Cornell Medicine, New York, NY 10065, USA
6 Lead contact
* Correspondence: huangfud@mskcc.org
AUTHOR CONTRIBUTIONS

Conceptualization, S.J.K. and D.H.; methodology, S.J.K., J.Y., W.W., J.P., H.S.C., R.L., D.M., K.K.D., E.A., C.S.L., and DH; investigation, S.J.K., J.Y., W.W., J.P., H.S.C., J.L.I., J.K., J.Z., Q.L., D.M., and R.L.; visualization, S.J.K., W.W., and Q.L.; funding acquisition, D.H.; project administration, D.H.; supervision, S.J.K. and D.H.; writing – original draft, S.J.K. and D.H.; writing – review & editing, all authors.

11 9 2024
27 8 2024
21 8 2024
17 9 2024
43 8 114640114640
https://creativecommons.org/licenses/by-nc/4.0/ This is an open access article under the CC BY-NC license (http://creativecommons.org/licenses/by-nc/4.0/).
SUMMARY

Functional enhancer annotation is critical for understanding tissue-specific transcriptional regulation and prioritizing disease-associated non-coding variants. However, unbiased enhancer discovery in disease-relevant contexts remains challenging. To identify enhancers pertinent to diabetes, we conducted a CRISPR interference (CRISPRi) screen in the human pluripotent stem cell (hPSC) pancreatic differentiation system. Among the enhancers identified, we focused on an enhancer we named ONECUT1e-664kb, ~664 kb from the ONECUT1 promoter. Previous studies have linked ONECUT1 coding mutations to pancreatic hypoplasia and neonatal diabetes. We found that homozygous deletion of ONECUT1e-664kb in hPSCs leads to a near-complete loss of ONECUT1 expression and impaired pancreatic differentiation. ONECUT1e-664kb contains a type 2 diabetes-associated variant (rs528350911) disrupting a GATA motif. Introducing the risk variant into hPSCs reduced binding of key pancreatic transcription factors (GATA4, GATA6, and FOXA2), supporting its causal role in diabetes. This work highlights the utility of unbiased enhancer discovery in disease-relevant settings for understanding monogenic and complex disease.

Graphical Abstract

In brief

Kaplan et al. performed a CRISPRi repression screen in an hPSC pancreatic differentiation system, identifying enhancers of pancreatic lineage transcription factors. They validated a long-range enhancer ~664 kb from the ONECUT1 promoter and linked a T2D-associated variant in this enhancer to ONECUT1 through enhancer deletion and variant editing.
==== Body
pmcINTRODUCTION

A spectrum of genetic causality often influences the onset, penetrance, and severity (expressivity) of a disease.1 This concept is exemplified in metabolic conditions like diabetes, where a single nucleotide alteration in a protein-coding region can cause severe neonatal diabetes mellitus (NDM) or maturity-onset diabetes of the young (MODY), while the aggregation of polygenic effects is thought to confer susceptibility to type 2 diabetes (T2D).2,3 However, most genetic variants have moderate effects and are in unannotated non-coding genomic regions, with the affected genes typically unknown.3,4 Therefore, there is a pressing need to understand the functional impact of disease-associated variants uncovered by genome-wide association studies (GWASs) and manage the burgeoning set of variants of unknown significance identified through the wide adoption of high-throughput DNA sequencing.

It is hypothesized that many genetic variants affect enhancers,5 and alterations in their sequences can affect gene expression, contributing to both Mendelian and complex disease traits. A clear illustration of this impact is shown by the limb malformation caused by a mutation in an SHH enhancer ~1 Mb from the gene promoter.6 This represents one of a handful of clinical genetics examples where enhancer mutations have a pronounced pathological impact.7 On the opposite end of the disease severity spectrum, functional consequences of individual sequence alterations in common complex diseases are generally unknown, but a significant proportion of disease-associated genomic variations identified through GWASs are located within putative enhancers.5 Together, these findings have contributed to the inception of the enhanceropathy concept.7 However, diagnosis of suspected enhanceropathies is hindered by the lack of functional annotation for most enhancers. While putative enhancers can be identified based on genomic features like chromatin accessibility, scalable functional characterization through the perturbation of endogenous genomic loci in disease-relevant cell contexts remains a considerable challenge. One obstacle is the variable distance enhancers can have from their target genes. In our recent interrogation of developmental enhancers, we found that a significant number of identified enhancers were located >100 kb away from the target gene promoter,8 yet practical considerations often limit functional enhancer screens to regions within 100 kb of the target gene promoter.9 Consequently, the overall prevalence of long-distance enhancer-promoter regulations remains unclear despite a growing set of long-range enhancers characterized through genetic studies.10

Human pluripotent stem cell (hPSC) differentiation is a powerful system for disease modeling. The integration of genetic perturbation and genomic approaches has proven instrumental for investigating the cascade of gene and enhancer activation during development and holds promise for understanding long-term health implications. We and others have shown that NDM gene knockout in hPSCs recapitulates key aspects of human genetic conditions.11 In addition, the hPSC platform has been used to study the impact of clinically observed recessive mutations in a PTF1A enhancer, shedding light on its role in pancreatic agenesis and NDM.12,13 However, many NDM cases lack a known genetic cause even with whole-exome sequencing, underscoring the importance of identifying non-coding regulatory regions.14 Significant progress in biochemical characterization has narrowed the search space, but many putative enhancers have not been functionally interrogated. Furthermore, the extensive epigenome rewiring during development further swells the number of potential enhancers, complicating the characterization task.15,16 This challenge is becoming more tractable with the use of CRISPR-Cas9 and dCas9-based CRISPR interference (CRISPRi) tools, enabling the simultaneous perturbation of genes or enhancers en masse.9 We and others have applied large-scale CRISPR-Cas9 screens in hPSCs to uncover genes important for gastrulation and pancreatic development.17-22 Recent CRISPRi screens have further identified enhancers and long non-coding RNAs (lncRNAs) involved in hPSC endoderm and cardiac differentiation.8,23,24 These advances motivated us to leverage CRISPRi tools to interrogate putative enhancers identified through our epigenomic characterization of the hPSC pancreatic differentiation system25 to uncover functional enhancers with potential roles in diabetes.

The hPSC pancreatic differentiation system directs cells progressively through the definitive endoderm (DE), primitive gut tube (GT), and posterior foregut (PFG; also known as PP1) stages that model embryonic development in humans.26 Building on our recent success in using a cell identity gene as a readout for enhancer discovery in cell-state transitions,8 we utilized the developmental expression initiation of the diabetes-associated gene PDX1 at the PFG stage as a readout for the onset of pancreatic differentiation.11,17 We then selected candidate transcription factors (TFs) likely to influence pancreatic specification and hypothesized that repressing enhancers of these TFs would impact PDX1 activation. Our subsequent screen indeed uncovered enhancers of multiple TFs that affected PDX1 expression at the PFG stage. We focused on a distal enhancer of ONECUT1, located ~664 kb from the gene promoter. This choice was motivated by the presence of a high-confidence fine-mapped T2D-associated single-nucleotide polymorphism (SNP) within this enhancer27 and clinical genetics findings linking ONECUT1 coding mutations to NDM.28,29 We show that ONECUT1 enhancer deletion in hPSCs leads to an almost complete absence of ONECUT1 transcript and protein, resulting in impaired pancreatic differentiation. Enhancer deletion had locus-wide effects within the topologically associating domain (TAD), with reduced antisense transcription and decreased H3K27ac levels. Finally, we generated hPSC lines harboring the T2D variant (rs528350911) and found that it significantly impaired the binding of the key pancreatic TFs GATA4, GATA6, and FOXA2. Our study not only identified an essential enhancer of ONECUT1 in pancreatic development but also developed and reduced to practice a framework for prioritizing disease-associated variants for investigation.

RESULTS

A CRISPRi screen for pancreatic enhancers

To discover enhancers involved in pancreatic development and diabetes, we conducted a CRISPRi screen in an hPSC pancreatic differentiation system, using PDX1 expression as a readout. PDX1 is a key marker for the initiation of pancreatic differentiation and is indispensable for mammalian pancreatic development.30-33 Homozygous and heterozygous mutations in PDX1 cause NDM and MODY, respectively.34,35 Additionally, the PDX1 locus is associated with T2D risk.36 To identify putative enhancer regions for study, we first selected 23 candidate TFs (including PDX1) based on their likely impact on PDX1+ induction and diabetes by integrating findings from murine experiments,30,31,37-55 clinical genetics, and our previous CRISPR-Cas9 screens in hPSC pancreatic differentiation (Figure 1B). 11 of these TFs have established roles in human metabolic disorders or NDM.29,34,56-61 15 TFs, including GATA6, RFX6, and ONECUT1, have also been shown to affect pancreatic differentiation and PDX1 expression in hPSC disease models or CRISPR-Cas9 coding screens for regulators of endoderm and pancreatic differentiation.17,18,28,62,63

TADs have been proposed to constrain regulatory enhancer-promoter interactions,64 so we conducted Hi-C assays on PFG-stage cells and used TAD information to delineate an ~2-6-Mb linear window around each TF for interrogation (Table S1). Within each window, we collated all accessible regions identified based on ATAC-seq peaks across the DE, GT, and PFG stages for interrogation (Figure 1C; Table S1).25 These accessible regions also included 758 annotated promoters, facilitating benchmarking of the CRISPRi screen results against our previous genome-scale CRISPR-Cas9 coding screen conducted in the same differentiation context.17 Thus, we designed a library with 37,184 gRNAs targeting 6,443 regions (average length of ~300 bp), 1,100 safe-targeting,65 and 463 non-targeting gRNAs66 and cloned the sequences into a lentiviral backbone (Table S1).

The CRISPRi screen identifies enhancers and genes necessary for human pancreatic differentiation

We conducted a pooled CRISPRi screen in our hPSC pancreatic differentiation system17 with two biological replicates (independent infections and differentiations), referred to as screens 1 and 2 (Figure 1A). For each screen, we infected inducible dCas9-KRAB H1 hPSCs67 with the gRNA library and induced dCas9-KRAB expression. Following differentiation, we enriched for PDX1+ and PDX1− cells at the PFG stage through fluorescence-activated cell sorting (FACS) (Figures 1D and S1A). Subsequently, next-generation sequencing determined gRNA abundance within the sorted populations, and Model-based Analysis of Genome-wide CRISPR/Cas9 Knockout (MAGeCK) analysis68 identified 38 top hits that caused a PDX1 decrease when perturbed (Figures 1E and S1B; Table S1). Both promoter (20) and non-promoter (18) hits have greater sequence conservation compared to all regions investigated, as shown by the average conservation scores calculated from PhyloP10069 (Figure S1C). Furthermore, we examined all interrogated promoter regions and found the gene hits to be largely concordant with our previous Cas9 coding screen conducted in a similar pancreatic differentiation context17 (Figure 1F). These findings support the utility of CRISPRi repression for the discovery of both coding and noncoding regulators of hPSC differentiation.

We next focused on non-promoter hits and examined H3K27ac, a histone mark associated with active enhancers.70 Non-promoter hit regions did not have H3K27ac peaks in undifferentiated hPSCs but gained H3K27ac at various stages during pancreatic differentiation (Figures 1G, 1H, and S1E), indicating stage-specific enhancer activity. We further reasoned that each enhancer hit is likely to correspond to a promoter hit for the regulated gene. Indeed, the majority of discovered enhancers were linked to a single promoter hit within the linear neighborhood, resulting in a total of 15 enhancer-gene pairs (Figures 1I and S1D). This approach demonstrates the specificity of the enhancer-gene regulatory mechanism. In addition, we observed that the enhancer-promoter distances ranged from 2 to 664 kb, with three enhancers found beyond 100 kb from the annotated transcription start site (TSS) (Figure 1J). Thus, unbiased screening in hPSC pancreatic differentiation uncovered many previously undiscovered enhancers and expanded our understanding of gene regulation during human pancreas development.

ONECUT1e-664kb deletion causes allele-specific loss of ONECUT1 expression

The CRISPRi screen identified the ONECUT1 promoter as well as a site ~664 kb from the ONECUT1 TSS, hereafter denoted as ONECUT1e-664kb. We observed that, at the PFG stage, the ONECUT1 promoter showed significantly increased Hi-C contact with a region encompassing ONECUT1e-664kb compared to the DE stage (Figures 2A and 2B). This observation suggests a physical enhancer-promoter interaction and substantial enhancer rewiring during differentiation. We chose ONECUT1e-664kb for further investigation because ONECUT1 coding mutations cause pancreatic hypoplasia in humans,28,29 and ONECUT1 is necessary for efficient hPSC differentiation to pancreatic progenitors.28,71 Consistent with the clinical findings, ONECUT1-null mice had a reduced number of endodermal cells expressing PDX1 at pancreatic budding,55 accompanied by compromised pancreatic organogenesis.72-74 To validate the regulatory impact of ONECUT1e-664kb on ONECUT1, we used CRISPR-Cas9 deletion as an orthogonal method to complement our dCas9-KRAB-based repression approach. Using paired gRNAs flanking an ~750-bp region, we generated two heterozygous and one homozygous clonal enhancer deletion lines (Figure 2C). ONECUT1 transcription is activated during the GT-to-PFG transition (Figure S2A), so we differentiated the homozygous and heterozygous enhancer knockout (eKO) lines to the PFG stage and assessed PDX1 and ONECUT1 protein levels via flow cytometry (Figure 2D). Homozygous eKO led to a near-complete absence of ONECUT1+ cells (Figure 2E). In addition, the percentage of cells activating PDX1 expression was significantly lower in homozygous eKO cells compared to unedited wild-type (WT) cells, showing that the loss of ONECUT1 impaired pancreatic differentiation (Figure 2F). These findings were further confirmed through immunofluorescence staining for PDX1 and ONECUT1 (Figure S2B). Further differentiation to the pancreatic progenitor stage (PP2) (Figure S2C) revealed that homozygous eKO cells remained ONECUT1 negative (Figures S2F and S2G). In addition, they formed a significantly reduced percentage of PDX1+NKX6-1+ cells compared to WT cells (Figures S2D and S2E), recapitulating the pancreatic differentiation phenotype reported for ONECUT1 gene KO.28 To explore the mechanisms by which ONECUT1 may be affecting PDX1 expression, we examined the PDX1 enhancers uncovered in our screen (Table S1) for ONECUT1 binding.17 While all PDX1 enhancers gained ATAC-seq and H3K27ac chromatin immunoprecipitation sequencing (ChIP-seq) signals during differentiation, we observed ONECUT1 ChIP-seq signal primarily at PDX1e-2kb, suggesting that ONECUT1 may activate PDX1 through this enhancer (Figure S2I). Additionally, this enhancer region is bound by key pancreatic TFs, including HHEX, FOXA2, GATA4, and GATA6. This is a potential mechanism by which the homozygous eKO cells could overcome the lack of ONECUT1 and activate PDX1 at the PP2 stage.

Heterozygous eKO lines did not exhibit a significantly lower percentage of PDX1+ or ONECUT1+ cells, but there was a significant decrease in the mean fluorescence intensity (MFI) signal of ONECUT1 among the PDX1+ cells in heterozygous cells compared to WT cells at both the PFG and PP2 stages (Figures 2G, 2H, and S2H). This suggested that the heterozygous cells had lower ONECUT1 expression compared to WT cells. To directly assess the cis effect of heterozygous eKO on ONECUT1 transcription, we pursued a more quantitative approach and designed an allele-specific droplet digital PCR (ddPCR) assay for the ONECUT1 transcript (Table S2). Given the substantial linear genomic distance between ONECUT1e-664kb and the ONECUT1 gene, we first called and phased all SNPs in the ONECUT1 locus after conducting targeted long-read Nanopore sequencing with adaptive sampling (Figure S2J). This approach allowed us to identify a SNP in the 3′ UTR of ONECUT1 and design an allele-specific ddPCR assay (Figures 2I and 2K). We also identified a genotyping SNP close to ONECUT1e-664kb for phasing the deleted enhancer in the heterozygous lines via PCR. Both heterozygous lines had enhancer deletions on the same allele (designated “allele 2”) (Figure 2J). ddPCR assays performed on RNA extracted from cells differentiated to the PFG stage showed that homozygous eKO lines lost >95% of transcription from both alleles (Figure 2L). Heterozygous lines exhibited a specific loss of transcription from allele 2, resulting in a significantly skewed ratio of transcripts from allele 2 over allele 1. This skewed ratio was not observed in control ddPCR assays conducted on genomic DNA (gDNA) (Figure S2K). These findings demonstrate that the transcriptional impact of ONECUT1e-664kb deletion on ONECUT1 expression is due to direct in cis regulation rather than an indirect effect through the downregulation of the ONECUT1 TF or other mechanisms. Notably, one of the heterozygous eKO lines (Het. 1) had the enhancer sequence inverted on allele 1, and ONECUT1 transcription from allele 1 was not significantly affected by this inversion. Together, these results establish ONECUT1e-664kb as a critical cis-regulatory element for ONECUT1 activation in pancreatic development.

ONECUT1 enhancer deletion causes locus-wide transcriptional and epigenomic alterations

The large impact of ONECUT1e-664kb deletion on ONECUT1 transcription led us to investigate its effects on transcription in the broader genomic locus, since some enhancers may serve as “hubs” that simultaneously regulate the expression of multiple genes in a locus.75 We performed TruSeq stranded total (ribodepletion) RNA sequencing (RNA-seq) on WT, heterozygous, and homozygous eKO cells differentiated to the GT and PFG stages (Table S3). Consistent with our ddPCR results, homozygous eKO cells had an approximately 30-fold decrease in ONECUT1 reads at the PFG stage compared to reads in WT cells. We observed a statistically significant decrease in WDR72 expression in homozygous ONECUT1e-664kb eKO cells (Figure 3A), the only other protein-coding gene in the same TAD as ONECUT1, but did not detect any significant transcription changes of genes in the adjacent TADs. We attempted to verify the decrease of WDR72 expression in homozygous eKO cells by ddPCR assays. We observed a noticeable trend toward a decrease, but the difference did not reach statistical significance (Figure S3A), which may be attributed to the low expression of WDR72 at the PFG stage. These results show that ONECUT1e-664kb primarily regulates ONECUT1 transcription during pancreatic differentiation.

We further investigated the impact of ONECUT1e-664kb deletion on ONECUT1 promoter antisense (PAS) RNA transcription. Transcription-associated PAS RNAs have been extensively documented,76,77 and antisense RNAs can be quantified in strand-specific RNA-seq datasets.78 Indeed, we observed an adjacent non-coding region antisense to ONECUT1 that had lower antisense transcription in homozygous eKO cells compared to WT cells (Figure 3B). Quantifying antisense transcripts in this region across three differentiation replicates showed that, during the GT-PFG transition, ONECUT1 transcriptional activation in WT cells was associated with a significant increase in ONECUT1 antisense transcripts. These transcripts were significantly reduced in both heterozygous and homozygous eKO lines at the PFG stage (Figure 3C). In contrast, sense transcripts across the same non-coding region showed no change during the GT-PFG transition and were unperturbed upon enhancer deletion (Figure S3B). Thus, in addition to affecting the ONECUT1 transcript, enhancer deletion also affects PAS RNAs, illustrating an additional transcriptional phenotype of enhancer loss.

The transcriptional changes we observed upon enhancer deletion motivated us to investigate the effect of enhancer deletion on H3K27ac, an epigenetic mark associated with active enhancers70 and transcription.79 Additionally, a previous study documented that PTF1A enhancer deletion causes H3K27ac loss at the PTF1A promoter as well as adjacent regions in the genomic locus.12 However, given that PTF1A encodes a TF, it is challenging to distinguish trans and cis effects of enhancer deletion. We reasoned that analyzing allele-specific H3K27ac ChIP-seq reads could aid in examining cis consequences of enhancer deletion, with a cis effect manifesting as an imbalance in allele-specific reads within heterozygous eKO cells. We performed H3K27ac ChIP-seq on PFG cells differentiated from WT, heterozygous, and homozygous eKO lines (Table S4). First, we compared WT and homozygous eKO cells. The decrease in H3K27ac in the ONECUT1 TAD was the most pronounced of all regions on chromosome 15 (Figure S3C), strongly suggesting cis effects. Among the 16 regions with the most significant H3K27ac reduction (adjusted p value < 1E-15) (Figure 3D), 10 contained heterozygous SNPs. Since these SNPs were phased with long-read sequencing (Figure S2J), we could use the SNPs to quantify changes occurring in cis with enhancer deletion by calculating allele-specific H3K27ac ChIP-seq read ratios in both WT and heterozygous eKO cells (Figure S3D). As controls, we used regions with heterozygous SNPs that did not meet the significance threshold (adjusted P ≥ 1E-15) of the differential analysis between homozygous eKO and WT. For all top differential peaks, H3K27ac was decreased in cis to the allele with enhancer deletion in heterozygous eKO cells, with some positions losing nearly all H3K27ac (Figure 3E), whereas control regions showed a near 1:1 allele read ratio average in both WT and heterozygous eKO cells (Figure 3F). This demonstrates that ONECUT1e-664kb affects H3K27ac at multiple genomic positions within the locus, supporting a direct crosstalk among regions marked by H3K27ac. Thus, the reduction of H3K27ac signal and PAS transcripts in cis with enhancer deletion illustrates the impact of enhancer dysregulation on genome function beyond gene transcription.

A T2D-associated SNP decreases pancreatic TF recruitment to ONECUT1e-664kb

To explore the significance of ONECUT1e-664kb in a clinical genetics context, we examined SNP metabolic disease association statistics in the ONECUT1 locus using data from a recent T2D GWAS meta-analysis27 and other studies.80,81 We observed one distinct T2D disease signal around the ONECUT1 gene with a largely unresolved 99% credible set spanning an ~100-kb genomic region (hg38, chromosome 15 [chr15]: 52777944–52873484). We identified a second independent T2D disease-causal signal ~664 kb away from ONECUT1 and closer to WDR72 that is driven by a single variant, rs528350911 (G/C, reverse strand, minor allele frequency = 0.0068 in individuals of European ancestry27 or 0.0039 in all populations based on gno-mAD v.4.1.0) (Figure 4A). This variant is directly within ONECUT1e-664kb and is highly confidently fine-mapped with the posterior probability of association (PPA) close to 1 (0.994) after linkage disequilibrium correction that is part of Sum of Single Effects (SuSiE) fine-mapping.27,82 To predict whether rs528350911 would affect enhancer activation in pancreatic development, we trained a ChromBPNet83,84 model on PFG-stage ATAC-seq data.25 Then, we computed per nucleotide DeepLift85 attribution scores across the region for both the reference and alternate sequences of rs528350911. The decreased attribution for an entire GATA TF motif suggests that the G in the GATA motif is crucial for the model’s genomic accessibility prediction, and the alternate sequence associated with T2D risk leads to an attribution loss at the site (Figure 4B).

GATA4 and GATA6 are members of the GATA-family TFs known for their important role in pancreatic differentiation62 and are linked to neonatal diabetes.86,87 At the PFG stage, GATA4 and GATA6, along with FOXA2, bind to the ONECUT1e-664kb enhancer region (Figure S4A).25,62 This motivated us to investigate the molecular impact of rs528350911. We engineered three heterozygous hPSC lines with the risk variant: two carried the variant on allele 2 and one on allele 1 (Figure 4C). While all lines differentiated to the PFG stage had a similar percentage of ONECUT1- and PDX1-expressing cells (Figures 4D and S4B) and did not affect WDR72 expression (Figure S4C), we observed a significant reduction in ONECUT1 MFI for two of the three lines (Figure 4E). Allele-specific ddPCR quantification of ONECUT1 transcription (Figure 2K) confirmed that the same two lines had a small but significant allele-specific reduction in ONECUT1 transcription (Figure 4F).

We further investigated the impact of rs528350911 on TF binding. To accomplish this, we developed an allele-specific ddPCR assay targeting the reference and alternate (T2D-associated) nucleotides (Figure 4G) to quantify TF binding following ChIP in an allele-specific manner. By utilizing the heterozygous lines, this allele-specific assay supports an internally controlled experimental system that mitigates technical variabilities arising from in vitro differentiation, ChIP, and PCR. We selected two hPSC lines that contained the risk variant on different alleles; one line had a significant impact on ONECUT1 expression, while the other did not show a significant effect, as mentioned above. After differentiating the lines to the PFG stage, we performed ChIP for GATA4, GATA6, and FOXA2. ddPCR on the precipitated DNA revealed an ~4-fold reduction in GATA4 binding on the allele containing the risk variant in both lines (Figure 4H). We also observed significantly decreased GATA6 and FOXA2 binding. Thus, despite a marginal effect on ONECUT1 transcription, the risk variant substantially reduced TF binding to ONECUT1e-664kb.

DISCUSSION

Pancreatic enhancers have been studied in mice by integrating pre-existing knowledge of important pancreatic genes with annotations such as evolutionary conservation and chromatin accessibility. Using this approach, several murine PDX1 enhancers were identified, and the mouse models showed a spectrum of pancreatic phenotypes depending on the specific region that was removed.88,89 A similar approach also identified a putative GATA4 enhancer, which was subsequently verified by a transgenic reporter assay in mouse endoderm and endoderm-derived tissues, including the pancreas and duodenum.90 These studies set the stage for integration of biochemical marks into the enhancer discovery framework, leading to the discovery of a PTF1a enhanceropathy responsible for pancreatic agenesis and NDM.12,13 The studies also highlight the value of developmental enhancer characterization for elucidating mechanisms of organogenesis and improving disease diagnosis. Building upon these findings, our study expanded the scope of enhancer discovery by systematically targeting developmentally accessible chromatin regions for large-scale functional interrogation using dCas9-KRAB. Our study identified the human PDX1 and GATA4 enhancers corresponding to the aforementioned mouse enhancers (Table S1) and uncovered additional enhancers for GATA6, HHEX, HNF1B, ISL1, MEIS1, ONECUT1, RFX6, and PBX1 that had not been reported previously. Our findings provide a valuable diagnostic resource not only for the large number of monogenic diabetes patients with unknown genetic cause14 but also for comprehending variants identified through population-level association studies.

Of all the discovered enhancers, this work adds ONECUT1e-664kb to a select group of well-documented long-range regulatory sequences that activate transcription over genomic intervals exceeding 500 kb.91-99 These enhancers are valuable models for further investigation into long-range transcriptional regulation and its implications for health and disease. Our findings underscore the value of unbiased discovery, and we speculate that gene regulation by long-range enhancers is more common than appreciated. From a technical standpoint, the substantial genomic distance between ONECUT1e-664kb and its target gene presented a challenge for definitively demonstrating the cis impact of the enhancer. In addition to conducting Hi-C to support enhancer-promoter interaction, we phased the variants in the ONECUT1 locus and employed high-sensitivity variant-specific ddPCR assays in heterozygous lines. These assays revealed a near-complete, allele-specific loss of ONECUT1 transcription resulting from the deletion of ONECUT1e-664kb. In addition, we uncovered two mechanisms beyond transcriptional control by which regulatory regions could influence disease traits: loss of PAS transcription and H3K27ac in a genomic locus. Considering results from ONECUT1 KO studies in mice55 and hPSC pancreatic differentiation,28 we infer ONECUT1 protein loss to be the driver of the pancreatic differentiation phenotype in our eKO study. We predict that loss of ONECUT1e-664kb in an individual would phenocopy the severe pancreatic hypoplasia observed in patients with ONECUT1 protein-coding mutations.28,29 Supporting this prediction, a proband with NDM was recently discovered to have an ~100-kb deletion encompassing ONECUT1e-664kb in trans with a ONECUT1 protein-truncating frameshift mutation.100

The discovery of ONECUT1e-664kb also helped us prioritize a disease risk variant for functional interrogation. In previous GWASs, the T2D-associated variant rs528350911 is assigned to the closest gene—WDR72.27 Considering the strong effect of ONECUT1e-664kb deletion on ONECUT1 expression and the known association of ONECUT1 with diabetes, it is likely that this variant confers disease susceptibility through its effect on ONECUT1 expression. Allelic series are primarily reported in the context of protein-coding variants with variable impact on gene function.101,102 Our findings support an allelic series for ONECUT1 that spans protein-coding mutations,28,29 enhancer deletion, and a point variant within the enhancer with a spectrum of impacts on ONECUT1 transcriptional regulation, pancreas development, and, ultimately, disease risk and diabetes phenotypes.

Our study demonstrates how the integration of CRISPR screening, precise genome editing, deep learning models, and 3D genomic conformation information can improve the functional pairing of variants with genes. The introduction of rs528350911 into hPSCs did not have a strong effect on ONECUT1 transcription, but its robust effect on TF binding to ONECUT1e-664kb suggests that assays other than gene expression could be valuable for characterizing the molecular consequences of variants in endogenous loci. We speculate that rs528350911 may affect ONECUT1 expression when the variant effect is integrated over time. In addition, the variant may confer a latent disease vulnerability that could be exposed by factors not fully modeled in our current differentiation system, such as age and environment.103,104 Given the considerable distance between ONECUT1e-664kb and the ONECUT1 gene as well as the marginal impact of the rs528350911 risk variant on ONECUT1 transcription, it is unlikely that this variant would have been prioritized for investigation in a high-throughput variant modeling screen that relies on gene expression readouts coupled with prime or base editing. In comparison, dCas9-KRAB supports more effective enhancer repression and has a more pronounced impact on gene expression, making it a scalable approach for prioritization. Our work demonstrates that integrating in vitro functional enhancer discovery with population genetics provides a rational approach to efficiently prioritize variants for investigation. This framework holds promise for understanding the exponentially growing set of disease-associated variants discovered through population-wide GWASs and whole-genome sequencing efforts.

Limitations of the study

Our screens may not be sufficiently sensitive to detect weak enhancers or those masked by redundancies. Additionally, we were not able to identify the corresponding promoter hits for three of the enhancer hits, which could be due to variations in gRNA efficiency and screen cutoff stringency.

T2D risk variants may contribute to diabetes over many years and can impact non-pancreatic tissues, such as the liver. While we have provided molecular evidence supporting the role of rs528350911 in embryonic pancreatic differentiation, further studies are needed to fully characterize the roles of this variant in diabetes pathogenesis.

rs528350911 has been linked to LINC02490 (~330 kb away) in some databases. However, the intergenic region between ONECUT1 and WDR72 contains several lncRNAs with inconsistent annotations across different transcriptome assemblies. Together with their low expression levels, this made it challenging to assess the impact of eKO on specific lncRNAs using our RNA-seq data. Therefore, we opted to evaluate them collectively in the PAS analysis, which does not differentiate the effects on individual lncRNAs.

STAR★METHODS

RESOURCE AVAILABILITY

Lead contact

Further information and requests for resources and reagents should be directed to Dr. Danwei Huangfu (huangfud@mskcc.org).

Materials availability

All the hPSC cell lines and plasmids generated in this study are available upon reasonable request.

Data and code availability

All the sequencing data generated in this study are available at GEO under accession code GEO: GSE267330. The Hi-C data from the PFG cells are also available in the 4DN Data Portal (https://data.4dnucleome.org) under accession number 4DNFIQBPRZZD. Previously published data that were reanalyzed here are available at GEO under codes GEO: GSE114102 (ATAC-seq data used for screen design, and ATAC-seq and H3K27ac ChIP-seq visualization in Figure S1E) and GATA4, GATA6, and FOXA2 ChIP-seq visualization in and GEO: GSE181480 (HHEX and ONECUT1 ChIP-seq), and the 4DN Data Portal, accession numbers: 4DNFIKCEWMIA, 4DNFI7RGXYFY, 4DNFIRLG5UWL. All genomic coordinates reported in manuscript are hg38 unless otherwise indicated.

This paper does not report original code.

Any additional information required to reanalyze the data reported in this work is available from the lead contact upon request.

EXPERIMENTAL MODEL AND STUDY PARTICIPANT DETAILS

Primary cells

The parental hESC line H1 (NIHhESC-10-0043) is available from WiCell under a material transfer agreement.

hPSC culture

Experiments were performed with H1 human embryonic stem cells (NIHhESC-10-0043), which were regularly confirmed to be mycoplasma-free by the Memorial Sloan Kettering Cancer Center (MSKCC) Antibody & Bioresource Core Facility. All experiments were conducted per NIH guidelines and approved by the Tri-SCI Embryonic Stem Cell Research Oversight (ESCRO) Committee. hPSCs were maintained in Essential 8 (E8) medium (Thermo Fisher Scientific, A1517001) on vitronectin (Thermo Fisher Scientific, A14700) pre-coated plates at 37°C with 5% CO2. The Rho-associated protein kinase (ROCK) inhibitor Y-27632 (5 μM; Selleck Chemicals, S1049) was added to the E8 medium for 24 h when passaging or thawing hPSCs.

METHOD DETAILS

hPSC directed pancreatic differentiation

hPSCs were maintained in E8 medium for 2 days to reach ~80% confluence. Cells were washed with PBS and differentiated to DE, GT, PFG, and PP2 stages following previously described protocols.17,105 hPSCs were rinsed with PBS and first differentiated into DE using S1/2 medium supplemented with 100 ng/mL Activin A (Bon Opus Biosciences) for 3 days and CHIR99021 (Stemgent, 04-0004-10) for 2 days (first day, 5 μM; second day, 0.5 μM). DE cells were rinsed with PBS and then exposed to S1/2 medium supplemented with 50 ng/mL KGF (FGF7) (PeproTech, 100-19) and 0.25 mM vitamin C (VitC) (Sigma-Aldrich, A4544) for 2 days to reach GT stage. GT cells were then switched to S3/4 medium supplemented with 50 ng/mL FGF7, 0.25 mM VitC and 1 μM retinoic acid (RA) (Sigma-Aldrich, R2625) for 2 days to reach PFG stage. PFG cells were then switched to S3/4 medium supplemented with 2 ng/mL FGF7, 0.25 mM VitC, 0.1 μM RA, 200 nM LDN, 0.25 μM SANT-1, 100 nM TPB and 1:200 ITS-X for 4 days to reach PP2 stage.

gRNA library assembly

23 transcription factors expressed in pancreatic differentiation were selected based on known clinical roles in metabolic disease and demonstrated pancreatic functional requirement from CRISPR screens.29,56-59 Around each transcription factor, loci within hg19 (lifted over in supplemental tables to hg38, Table S1) were defined to include at least the TAD containing the gene, with TADs called on PFG-stage Hi-C data as previously done.67 To select putative enhancers for interrogation in these loci, Genrich (0.6.1)106 was used to call peaks at an FDR-adjusted p value of 0.01 on previously published DE, GT, and PFG stage ATAC-seq data (Table S2).25 For each region in the loci, and some additional regions, CHOPCHOPv3107 was used to design gRNAs. For regions <750 bp, the top 5 ranked gRNAs were taken, and for regions >750 bp, 1 gRNA/150bp of sequence was selected for interrogation. Finally, 1100 safe harbor65 and 463 non-targeting gRNAs66 were included in the library for a total of 38,747 gRNAs. Oligos were synthesized, amplified, and restriction cloned into the lentiGuide-puro108 (Addgene: 52963; RRID: Addgene_52963) backbone by the MSKCC Gene Editing & Screening Core Facility. Cloned plasmid libraries were PCR amplified to incorporate adapters for NGS. Samples were purified and sequenced using Illumina HiSeq 2500 platform. FASTQ files were clipped by position and reads were mapped back to the reference library file to quantify abundance of reads per gRNA. The overall representation of the libraries was charted over a one-log fold change to evaluate representation in the final library.

Lentiviral gRNA library production

gRNA lentiviral library generation was performed as previously described.17 In brief, a total of 9.45 μg of the library plasmid combined with 6.75 μg lentiviral packaging vector psPAX2 and 1.36 μg vesicular stomatitis virus G (VSV-G) envelope expressing plasmid pMD2.G (Addgene plasmids 12260 and 12259) were transfected with the JetPRIME (VMR; 89137972) reagent into 10E6 293T cells in a 10cm plate to produce the lentiviral particles. Fresh medium was changed 24h after transfection and viral supernatant was collected at 48 and 72h after transfection, spun at 1000 rpm for 5 min, 0.45μm filtered, and stored at −80°C.

CRISPR dCas9-KRAB screens

Two independent screens were performed following the same procedure. H1 hPSC (i)dCas9-KRAB cells were used.67 The lentiviral library was transduced into 12E6 hPSCs distributed in three 10 cm plates following the conditions (MOI: 4, protamine sulfate: 6 μg/mL, Ri: 10 μM). After one day of recovery, doxycycline (DOX, 2 μg/mL) was added, and cells were kept under puromycin selection for 2 days. After two days, 30E6 cells were seeded into 6-well VTN-coated plates (5E5 per well) and grown for 48 h in E8 medium. After washing with Phosphate Buffered Saline (PBS) w/o Ca2+ & Mg2+, cells were differentiated, with DOX added throughout the seven-day differentiation process. At the PFG stage, cells were dissociated with TrypLE, stained with LIVE/DEAD reagent, fixed, stained with the PDX1 antibody, and stained with a fluorescent secondary antibody. They were then sorted by the MSKCC flow cytometry facility based on 2ndary antibody fluorescence using FACSAria sorters (BD Biosciences). Positive and negative cells were sorted aiming for a minimum 1000X representation (number of cells sorted*MOI)/(unique gRNAs in library) per condition, with two independent sorts serving as technical replicates for each biological replicate. For biological replicate 1, approximately 10E6 PDX1−and 10E6 PDX1+ cells were sorted. For replicate 2, approximately 20E6 PDX1−and 20E6 PDX1+ cells were sorted. Cells were pelleted and kept at −80°C for downstream gRNA sequencing.

gRNA sequencing and screen analysis

gRNA enrichment sequencing was performed by the MSKCC Gene Editing & Screening Core Facility as previously described, but with modifications to enable DNA extraction from fixed cells.17 Cells were decrosslinked and treated with proteinase K (0.1 mg/mL) overnight followed by genomic DNA extraction with QIAGEN Blood & Cell Culture DNA Maxi Kit (QIAGEN; 13362) and quantified by Qubit (Thermo-Scientific) following the manufacturer’s guidelines. A quantity of gDNA sufficient for 1000X representation of gRNAs was amplified with oligos containing the Illumina adapters and multiplexing barcodes by PCR. Amplicons were quantified by Qubit and Bioanalyzer (Agilent) and sequenced on the Illumina HiSeq 2500 platform. The gRNA library sequences were used to align the sequencing reads and the counts for each gRNA were determined with MAGeCK. Screens were analyzed with MAGeCK RRA (0.5.9.5), and the two biological screen replicates were analyzed independently (Table S1). To calculate a single value for the significance of effect, a ratio was taken between the positive and negative MAGeCK scores and log-transformed. The top 38 hits associated with a PDX1 decrease had a MAGeCK neg.FDR<0.15 in replicate 2 and caused a PDX1 decrease in both replicates (Table S1). Hits overlapped with annotated promoter regions were designated as “Promoter Hits,” with the rest designated as “Non-Promoter Hits.” Promoter overlap was defined as overlap with an interrogated region ±1.5kb of Gencode (v18) TSS sites. For Figure 1F, if there were multiple overlapped regions, the most significant hit was used for plotting. To link enhancers to the putatively regulated genes, the nearest linearly adjacent “Promoter Hit” was inferred to be the regulated gene of the “Non-Promoter Hit.” To determine hit and non-hit overlap with H3K27ac ChIP-seq peaks during differentiation, Genrich (0.6.1)106 was used to call peaks at an FDR-adjusted p value of 0.01 on previously published ES, DE, GT, and PFG stage H3K27ac ChIP-seq data.8,25 Conservation scores of regions were determined by calculating the average of per nucleotide conservation across all nucleotides within each region from PhyloP100 per-nucleotide conservation scores.69

Generation of ONECUT1e-664kb clonal knockout hPSC lines

H1 iCas9 hPSC lines62 were used to generate three enhancer deletion lines carrying ONECUT1e-664kb deletions. Two gRNAs targeting the flanking regions of the enhancer were designed, and knockouts generated as previously described but with some modifications.8 gRNAs and tracer RNA were ordered from IDT (Alt-R CRISPR-Cas9 crRNA, 1072532) and added at a 15 nM final concentration. In brief, gRNA/tracer RNA and Lipofectamine RNAiMAX (Thermo Fisher Scientific, 13778030) were diluted separately in Opti-MEM (Invitrogen, 31985070), mixed together, and incubated for 15 min at room temperature, and then added dropwise to freshly seeded iCas9 hPSCs in a 24-well plate. DOX (2 μg/mL) was added the day before transfection, the day of transfection, and 1 day after transfection to induce Cas9 expression. Three days after transfection, hPSCs were dissociated into single cells and ~1-3E3 cells were plated into a 10 cm tissue culture dish for colony formation. After ~10 days of expansion, single colonies were picked. Genomic DNA from crude cell lysate was used for PCR genotyping (primers: ONECUT1e_KOgenot_F, ONECUT1e_KOgenot_R). Heterozygous lines contained a band corresponding to the WT band (~1 kb) and a smaller band corresponding the knockout band (~300 bp), while homozygous lines only had smaller bands. Further genotyping was done through TOPO cloning of the genotyping amplicons which determined the exact genotypes of both alleles for all lines (Table S1). gRNA target sequences and primers used for PCR and sequencing are listed in Table S2.

Generation of clonal hPSC lines with disease-associated variant

A process similar to enhancer deletion line generation was used to generate heterozygous knock-in SNP lines carrying the rs528350911 variant. During transfection, a single gRNA proximal to the variant site was used, and a ssDNA HDR template containing the disease-associated variant was added to the mixture at the Lipofectamine transfection step (Table S2) at a concentration of 20 nM. Genomic DNA from crude cell lysate was used for PCR (primers: ONECUT1e_KOgenot_F, ONECUT1e_KOgenot_R). All picked clones were sequenced with both primers, and heterozygous lines picked by the presence of a double peak on the spectrogram that corresponded to the insertion of the disease-associated variant. For lines containing the introduced variant, no indels or other changes from the reference genome were detected in the amplicon.

ONECUT1 locus SNP phasing

To phase all SNPs in the ONECUT1 locus, Nanopore sequencing with adaptive sampling of the region chr15:52,707,645-53,586,740 (and controls regions) was conducted on an H1 iCas9 line with the rs528350911 heterozygous variant (line N11). DNA was extracted using Monarch HMW DNA Extraction Kit for Cells & Blood (NEB, T3050), according to the manufacturer’s instructions. DNA quantity was measured using a QuantiT (Thermo Fisher Scientific), purity on a NanoDrop (Thermo Fisher Scientific) and fragment-size distribution on a TapeStation (Agilent). Prior to ONT library preparation, DNA was sheared to ~15-20 kb fragment size using Covaris g-TUBE. Sequencing libraries were prepared from ~1–2μg of DNA using the native library preparation kit SQK-LSK114 according to the manufacturer’s instructions. Each library was loaded onto a PromethION flow cell R10.4.1 and sequenced on an ONT PromethION P24 device. The sample was run for 72 h with 1 nuclease flush and reload performed during the run to maximize sequencing yield.

Reads were aligned with minimap2 (2.26),109 with samtools110 coverage estimating ~138X coverage across the region enriched with adaptive sampling. WhatsHap (2.2)111 was used to phase the variants in the locus as well as assemble a haplotype containing the ONECUT1 gene, enhancer, and surrounding regions (Figure S2J). WhatsHap was also used to phase-tag reads for visualization in Figure S2J. From these results, a heterozygous SNP was identified in the 3′UTR of the ONECUT1 transcript, as well as the enhancer-proximal Genotyping SNP.

ONECUT1e-664kb heterozygous knockout line phasing

To determine which allele contained the knocked-out enhancer in the heterozygous enhancer deletion lines (Het. 1, Het. 2), a forward primer was designed centromeric to the Genotyping SNP and a reverse primer targeting the enhancer sequence that would be present on the non-deletion allele. Since the line designated as Het. 1 contained an inversion of the WT enhancer sequence, two different reverse primers were needed. Thus, for Het. 1, and amplicon was generated with the primer pair ONECUT1e_GenotSNP_F, ONECUT1e_GenotSNP_R_WT_inv_e, and for Het.2, the pair ONECUT1e_GenotSNP_F, ONECUT1e_GenotSNP_R_WT_e. The amplicons were sent for sanger sequencing with the ONECUT1e_GenotSNP_F primer used as the sequencing primer. The Genotyping SNP from this sequencing was used to determine the allele containing the deletion.

ONECUT1e-664kb heterozygous rs528350911 edited line phasing

A PCR amplifying an ~8.1kb fragment containing both the T2D-associated variant and the Genotyping SNP (primers: ONECUT1e_GenotSNP_F, ONECUT1e_GenotSNP_R_WT_e) with the LongAmp Taq polymerase kit was conducted (Table S2). The amplicon was sent to Plasmidsaurus for linear fragment sequencing. Reads were aligned with minimap2 (2.26), and variants were phased with longshot (0.4.5).112

Immunofluorescence staining and imaging

PFG-stage pancreatic differentiated cells were fixed in 4% paraformaldehyde (Thermo Fisher Scientific, 50980495) for 10 min at room temperature. After washing with PBST (PBS with 0.1% Triton X-100) three times, cells were blocked in 5% donkey serum in PBST buffer for 30 min at room temperature. Primary and second antibodies were diluted in the blocking solution (Table S2). Cells were incubated with primary antibodies overnight at 4°C, followed by 1 hr staining for secondary antibodies at room temperature. The cells were then stained with 4′,6-diamidino-2-phenylindole (DAPI) for ~15 min at room temperature. Images were taken using a confocal laser scanning platform (Leica TCS SP5).

Flow cytometry and analysis

Cells were dissociated using TrypLE Select and resuspended in FACS buffer (5% FBS in PBS). LIVE/DEAD Fixable Violet Dead cell stain (Invitrogen, L34955) was used to discriminate dead cells from live cells. LIVE/DEAD staining was performed per manufacturer’s instructions in FACS buffer, followed by fixation and intracellular staining with a FOXP3 staining buffer set (eBioscience, 00-5523-00) following the manufacturer’s instructions. Permeabilization/fixation was performed at room temperature for 30 min. Antibody staining was performed in permeabilization buffer for 30 min. Antibodies for this study are listed in Table S2. Cells were then analyzed using BD LSRFortessa. Flow cytometry analysis and figures were generated using FlowJo v.10.

RNA-seq and analysis

RNA was extracted from differentiated cells with the Zymo Quick-RNA MiniPrep kit (R1055). TruSeq stranded total (Ribodepletion) RNA-seq was performed on three independent differentiations by the MSKCC Integrated Genomics Operation (IGO) Core. Briefly, STAR alignment was performed followed byfeatureCounts113 quantification (Table S3) and normalization with DESeq2 median of ratios prior to visualization in Figures 3A, 3C, and S3B. Differential gene expression was performed and adjusted p values computed with DESeq2 (Table S3). Strand-specific analysis of promoter antisense transcription analysis was performed with featureCounts strand specific quantification of a region telomeric to ONECUT1 (chr15:52805588-53391232). For strand-specific transcript visualization, samtools was used to extract sense and antisense alignments from bam files from which bigwig coverage files were generated.

Droplet digital PCR (ddPCR) and analysis

Assays specific for the detection of SNPs (Table S2) were designed (using Primer3Plus, if needed), ordered through Bio-Rad, and conducted by the MSKCC Integrated Genomics Operation (IGO) Core. Cycling conditions were tested to ensure optimal annealing/extension temperature as well as optimal separation of positive from empty droplets. Optimization was done with a known positive control. A constant amount of RNA or DNA was used within each experiment. Reactions were partitioned into a median of ~17,500 droplets per well using the QX200 droplet generator. Plates were read and files analyzed with the QuantaSoft software to assess the number of positive droplets.

Procedure for RNA: after PicoGreen quantification of RNA, droplet generation was performed on a QX200 ddPCR system (Bio-Rad catalog # 1864001) using cDNA generated from 0.1–2 ng total RNA with the One-Step RT-ddPCR Advanced Kit for Probes (Bio-Rad catalog # 1864021) according to the manufacturer’s protocol with reverse transcription at 42°C and annealing/extension at 60°C.

Procedure for DNA: after PicoGreen quantification of DNA, 0.25–9 ng gDNA of ChIP-ed DNA was combined with locus-specific primers; FAM- and HEX-labeled probes; Hae III (DNA extracted form hPSC SNP lines), Hind III (ChIP-ed DNA), or Mse I (DNA extracted ONECUT1e-664kb hPSC knockout lines); and digital PCR Supermix for probes (no dUTP). Emulsified PCRs were run on a 96-well thermal cycler using cycling conditions identified during the optimization step (95°C 10′; 40 cycles of 94°C 30’ and 53°C (rs528350911 assays only) or 60°C 1’; 98°C 10’; 4°C hold).

Hi-C and analysis

Hi-C was performed as previously described with cells differentiated to the PFG stage processed with the Arima Hi-C kit (Arima, A510008).8 HiC-Pro (3.1.0)114 was used to align the individual Hi-C replicates to GRCh38, GCA_000001405.15, with alternative contigs removed, and then final merged libraries were obtained by combining all *.allValidPairs files where duplicates across biological replicates are maintained. DE libraries consisted of three biological replicates (two HUES8, one H1 with accessions 4DNFIKCEWMIA.hic, 4DNFI7RGXYFY.hic, and 4DNFIRLG5UWL.hic respectively), and the PFG libraries consisted of two biological replicates differentiated to the PFG stage (merged accession 4DNFIQBPRZZD.hic). For visualization, Hi-C data was visualized at 25 kB resolution.

Virtual 4C counts for the regions of interest were extracted using utilities from plotGardener (1.6.4)115 and plotted as a one-dimensional signal. Hi-C matrix counts for the regions of interest were extracted using plotGardener and plotted as a rectangular region. Observed counts are library normalized prior to plotting so that the total number of reads between panels is comparable. TADs that were previously identified were also overlaid on the data and can be found in the Hi-C subseries of GEO dataset published with this work. Counts within the vicinity of the ONECUT1 promoter (all contacts with start or end bins within a 1.2Mb window around the ONECUT1 promoter), were extracted with strawr (0.0.91) and passed to DESeq2 (1.40.1).116 The Wald test was used to assess significance of all differential contacts within that window, and significant hits anchored on the ONECUT1 containing bin are reported at an FDR of 0.1.

For calling topologically associating domains, merged Hi-C files were first ICE normalized through HiTC (1.34.0)117 and then topologically associating domains were identified with HiCDCplus (1.5.2)118 using its interface to TopDom at 50 kb. Window size for scanning was set at 5.

ChIP and ChIP-seq assays

ChIP was performed as previously described.18 For each sample, ~30E6 cells were crosslinked in-plate with 1% formaldehyde for 10 min at 37°C and quenched with 0.125 M glycine. Fixed cells were collected and washed in cold PBS buffer, snap frozen, and saved at −80°C. On the day of immunoprecipitation, the cell pellet was thawed on ice, resuspended in 700uL SDS buffer (1% SDS, 10 mM EDTA, 50 mM Tris–HCl, pH 8), and incubated for 10min on ice. Sonication was performed on a Branson Sonifier 150 set at 30% amplitude for 5.5min (10s on/off pulsing). Supernatant was pre-cleared with Dynabeads Protein G (#10004D) and then incubated overnight with antibodies (10 μg for transcription factors, 5 μg for H3K27ac, Table S2) at 4°C.). On the next day, the ChIP samples were incubated with Dynabeads Protein G for 6 h at 4°C. Then the beads were pelleted and sequentially washed twice with low salt (0.1% SDS, 1% Triton X-100, 2 mM EDTA, 20 mM Tris–HCl, pH 8, 150 mM NaCl), then high salt (0.1% SDS, 1% Triton X-100, 2 mM EDTA, 20 mM Tris–HCl, pH 8, 500 mM NaCl), and finally TE buffer (10 mM Tris–HCl, pH 8, 1 mM EDTA). The DNA was eluted from the beads by incubating in elution buffer (1% SDS, 0.1 M NaHCO3) at 65°C for 15min and decrosslinked with 190 mM NaCl at 65°C overnight. A total of 5 μL 0.5 M EDTA, 10 μL 1 M Tris–HCl, pH 6.5, and 1 μL Proteinase K (20 mg/mL) were added to 260 μL of the de-crosslinked product and incubated for 1 h at 45°C. DNA was isolated by using QIAquick PCR purification kit (Qiagen, 28104; QIAGEN). H3K27ac ChIP samples were submitted to MSKCC Integrated Genomics Operation core for NGS library preparation and sequencing. GATA4, GATA6, and FOXA2 ChIP samples were submitted for ddPCR.

H3K27ac ChIP-seq analysis

Sequenced data were aligned to the hg38 reference genome using bowtie2 (2.5.1).119 Peak calling was performed using MACS2 (2.2.9.1)120 with the respective ChIP input as the control and the default p value cutoff and default extension size. Irreproducible discovery rate (IDR)121 with a cutoff of 0.01 was used to filter reproducible peaks, where at least 2 pairwise comparisons of the three replicates passed the cutoff. Peaks were filtered by ENCODE exclusion list version 2 for hg38 (Table S4).122 The quantification and differential signal intensity of the samples were calculated in DESeq2.116 The integer counts, log(fold-change), and adjusted p values are provided in Table S4. The replicates were combined, and the signal track was generated from MACS2. MACS2 subtract was used to generate the difference of the signal track. All signal tracks were visualized using Integrative Genomics Viewer (IGV).123

For allele-specific analysis, GATK124 ASEReadCounter was used to quantify reads at heterozygous positions phased by long-read sequencing and ratios were calculated for heterozygous positions (score filters: PHASE_QUAL ≥ 140, and QUAL ≥ 200) with more than 100 reads total across all WT and Het. 2 replicates (Table S4). For Figure 3F, the “Phased Region” is the haplotype block containing ONECUT1 called by WhatsHap (chr15:50462878-54550826) and visualized in Figure S2J.

chromBPNet variant modeling

A ChromBPNet83,84 model was trained to predict the local ATAC profile using bulk ATAC-seq data from the PFG stage.25 Default hyperparameters were utilized with an input size of 2114 base pairs and an output size of 1000 base pairs for training both the bias model and the ChromBPNet model. For the bias model training, regions were selected that do not overlap with any ATAC-seq peaks in the four stages of pancreas differentiation,25 applying a bias threshold factor of 0.3 to filter out high-count regions in the sets. This resulted in a training set of 54957 regions. For ChromBPNet training, ATAC-seq peaks only from the PFG stage with IDR values greater than 830 were used, resulting in a training set of 52541 regions. During training for both models, peaks from chromosomes 15 and 13 for validation were held out. For each specific region of interest, the 2114 base pairs around the summit center was used as input and the DeepLift attribution score obtained using the trained ChromBPNet model. Specifically, the original attribution pipeline from the ChromBPNet GitHub repository was employed, which involves computing the DeepLift score over 20 reference inputs shuffled from the original sequence. To find the GATA motifs, Find Individual Motif Occurrences (FIMO)125 was used on the region of interest using the positional weight matrix (PWMs) of all GATA variants from the CIS-BP2.00 dataset,126 and motif hits with p values <0.05 were selected.

QUANTIFICATION AND STATISTICAL ANALYSIS

All datapoints refer to biological repeats. No statistical method was used to predetermine sample sizes. Unless otherwise indicated, p values were calculated with PRISM and indicated with star notation as follows: ns p > 0.05, *p ≤ 0.05, **p ≤ 0.01, ***p ≤ 0.001, ****p ≤ 0.0001. The investigators were not blinded to allocation during experiments and outcome assessment. No data were excluded from the analyses unless the differentiation experiment itself failed. The number of biological and technical replicates are reported in the legend of each figure. Flow cytometry analysis was derived from at least three independent experiments. For ChIP-seq and bulk RNA-seq, quantification and statistics were derived from at least three independent experiments. CRISPR dCas9-KRAB screening was performed twice. Quantification is shown as the mean ± s.d. All the statistical analysis methods are indicated in the figure legends and methods.

Supplementary Material

1

2

3

4

5

ACKNOWLEDGMENTS

We thank Dr. Thomas Vierbuchen for valuable advice. We thank Dr. Ting Zhou, Dr. Hanuman Kale, Aaron Zhong, and Emily DeBitetto for assisting with experiments not included in the manuscript and acknowledge the assistance from the following Memorial Sloan Kettering Cancer Center (MSKCC) Cores: Integrated Genomics Operation, Flow Cytometry, Gene Editing & Screening, Molecular Cytogenetics, and Stem Cell Research. We thank Dr. Ralph Garippa and Sanjoy Mehta for assistance with CRISPR library generation and sequencing, Cassidy Cobbs for providing advice regarding next-generation sequencing experiments, Dr. Stephanie Chrysanthou for conducting Nanopore sequencing with adaptive sampling, and Dr. Andrea Farina for designing and executing the ddPCR assays. This study was funded in part by National Institutes of Health grant U01DK128852 (to C.S.L., D.H., and E.A.), National Institutes of Health grant U01HG012051 (to D.H.), National Institutes of Health grant R01DK096239 (to D.H.), Starr Tri-I Stem Cell Initiative #2019-001 (to D.H. and E.A.), National Institutes of Health MSKCC Cancer Center Support Grant P30CA008748, and National Institutes of Health T32 training grant T32GM008539 (to S.J.K.).

Figure 1. CRISPRi repression screen to discover pancreatic differentiation enhancers

(A) hPSC stepwise pancreatic differentiation protocol schematic. ActA, Activin A; CHIR, CHIR99021; RA, retinoic acid; VitC, vitamin C.

(B) Gene selection rationale for the enhancer discovery screen.

(C) Putative enhancer region selection and gRNA design schematic.

(D) Screening procedure schematic.

(E) Results from two screen replicates; each point represents a genomic region. Highlighted are 38 regions with false discovery rate < 0.15 in replicate 2 that have a decrease in PDX1 in both replicates.

(F) Results from screen replicate 2 compared to the whole-genome PDX1 expression differentiation screen. Genes showing enrichment in both screens are labeled.

(G) Region overlap with H3K27ac peaks at any stage during differentiation (ES, DE, GT, and PFG).

(H) Differentiation stages at which non-promoter screen hits overlap with H3K27ac peaks and the number of hits with each activation pattern.

(I) Enhancer-gene pair assignment.

(J) Enhancer-gene pair linear distances.

Figure 2. ONECUT1e-664kb deletion and pancreatic differentiation characterization

(A) Hi-C data of DE- and PFG-stage differentiated hPSCs around the ONECUT1 locus. Contact between ONECUT1e-664kb and ONECUT1 promoter regions is circled to highlight the increase in contact frequency between the two loci. Red lines demarcate TAD boundaries identified at 50 kb.

(B) 1D slices of the 2D contact map were obtained with the ONECUT1 promoter set as the bin of origin. Hi-C data were library scaled and show the changes observed in (A). Significantly different bins were determined by DESeq2.

(C) hPSC enhancer deletion and genotyping PCR schematic.

(D) Representative flow cytometry plots of WT and eKO cells differentiated to the PFG stage.

(E) Percentage of cells achieving a ONECUT1+ identity.

(F) Percentage of cells achieving a PDX1+ identity.

(G) Representative histograms of ONECUT1 MFI of PDX1+ cells from flow cytometry of unedited and eKO cells differentiated to the PFG stage.

(H) Quantification and statistical comparison of ONECUT1 MFI of PDX1+ cells.

(I) Schematic for identifying ddPCR SNPs and genotyping SNPs.

(J) Spectrogram of genotyping SNP amplicon sequencing and schematic showing the allele that contained the deletion in heterozygous lines.

(K) ddPCR assay design schematic.

(L) Allele-specific ONECUT1 transcript quantification with ddPCR on RNA extracted from PFG-stage cells. Each allele is shown in a separate plot, with the expression ratio shown in the last plot. Each symbol represents one independent differentiation (n = 3 independent experiments) with two averaged technical ddPCR replicates. Data are presented as the mean ± SD. One-way analysis of variance (ANOVA) followed by Dunnett multiple comparisons test versus WT control.

For (E), (F), and (H), each dot represents one independent experiment (n = 4 or 5 independent experiments, 4 for Het. 1), and data are presented as the mean ± SD.

For all panels, ns p > 0.05, * p ≤ 0.05, ** p ≤ 0.01, *** p ≤ 0.001, **** p ≤ 0.0001.

Figure 3. Transcriptional and epigenetic consequences of ONECUT1e-664kb KO

(A) RNA-seq normalized read counts for genes in the ONECUT1 locus at GT and PFG stages from WT, heterozygous (het. 2), and homozygous eKO hPSCs. * Padj < 0.005 computed by DESeq2.

(B) Visualization of antisense RNA relative to the direction of ONECUT1 and read coverage from strand-specific RNA-seq of PFG-stage WT and homozygous eKO cells of representative differentiation replicates.

(C) Quantification of (B). Each dot represents one independent experiment (n = 3 independent experiments). One-way ANOVA followed by Dunnett multiple-comparisons test versus WT control.

(D) H3K27ac ChIP-seq at the PFG stage of WT, heterozygous, and homozygous eKO cells. TADs were called from PFG-stage Hi-C data. Difference was calculated via MACS2 subtract. Top differential peaks between WT and homozygous eKO cells were quantified with DESeq2. The ONECUT1e-664kb deletion region is within peak #12.

(E) H3K27ac is decreased in an allele-specific manner within significantly affected H3K27ac peaks.

(F) Quantification of average H3K27ac read ratios at all heterozygous SNP positions in the phased locus in the top significantly different peaks (WT vs. Homo, padj < 1E–15) and other peaks. Each point represents a heterozygous SNP position; average of three replicates from independent experiments. Ordinary one-way ANOVA with Šidák’s multiple-comparisons test with a single pooled variance. Data are presented as mean ± SD for (A), (C), (E), and (F). For all panels except (A), ns p > 0.05, * p ≤ 0.05, ** p ≤ 0.01, *** p ≤ 0.001, **** p ≤ 0.0001.

Figure 4. T2D-associated SNP hPSC modeling and pancreatic differentiation characterization

(A) LocusZoom plot of variants from the Mahajan et al.27 T2D GWAS meta-analysis unadjusted for BMI. Variants within credible sets are colored by the PPA score.

(B) Decreased DeepLift attribution score at the GATA motif upon introduction of the variant. The red box represents a GATA6 binding motif hit by FIMO (P < 0.001).

(C) Spectrogram of disease-associated SNP amplicon sequencing and schematic showing the allele that contains the introduced SNP in hPSCs.

(D) Percentage of ONECUT1+ cells at the PFG stage.

(E) Statistical comparison of ONECUT1 MFI of PDX1+ cells.

(F) ONECUT1 transcript ratio of edited alleles divided by unedited alleles. For each independent experiment shown, there were two averaged technical ddPCR replicates.

(G) ChIP-ddPCR schematic.

(H) ChIPed DNA ratio analysis. Each symbol represents one independent differentiation and ChIP experiment; two averaged technical ddPCR replicates ± SD. Two-way ANOVA followed by Tukey’s multiple-comparisons test versus WT control.

For (D)–(F), each symbol represents one independent experiment (n = 6 independent experiments), and data are presented as the mean ± SD. One-way ANOVA followed by Dunnett multiple-comparisons test versus WT control. For all panels, ns p > 0.05, * p ≤ 0.05, ** p ≤ 0.01, *** p ≤ 0.001, **** p ≤ 0.0001.

KEY RESOURCES TABLE

REAGENT or RESOURCE	SOURCE	IDENTIFIER	
Antibodies	
SOX17	BD BIOSCIENCES	Cat#561591; RRID: AB_10717121	
CXCR4	R&D Systems	Cat#FAB170A; RRID: AB_357073	
PDX1	R&D Systems	Cat#AF2419; RRID: AB_355257	
ONECUT1	Santa Cruz Biotechnology	Cat#sc-376308; RRID: AB_10988385	
NKX6-1	DHSB	Cat#F55A12; RRID: AB_532379	
Alexa Fluor 488	ThermoFisher Scientific	Cat#A-11055; RRID: AB_2534102	
Alexa Fluor 647	ThermoFisher Scientific	Cat#A-31571; RRID: AB_162542	
LIVE/DEAD™ Fixable Violet Dead Cell Stain Kit	ThermoFisher Scientific	Cat#L34964	
FOXA2	Cell Signaling Technology	Cat#8186; RRID: AB_10891055	
GATA4	Invitrogen	Cat#MA5-15532; RRID: AB_10989032	
GATA6	Cell Signaling Technology	Cat#5851; RRID: AB_10989032	
H3K27ac	EMD Millipore	Cat#MABE647; RRID: AB_2893037	
Bacterial and virus strains	
NEB 5-alpha Competent E. coli	NEB	Cat#C2987	
Chemicals, peptides, and recombinant proteins	
Complete Essential 8 (E8) medium	Thermo Fisher Scientific	A1517001	
Vitronectin (VTN-N)	Thermo Fisher Scientific	A14700	
EDTA	KD Medical	RGE-3130	
TrypLE Select 1X	Thermo Fisher Scientific	12563–029	
ROCK inhibitor Y-27632 2HCl	Selleck Chemicals	S1049	
Doxycycline	Sigma	D9891	
Lipofectamine RNAiMAX Transfection Reagent	Thermo Fisher Scientific	13778150	
Opti-MEM	Thermo Fisher Scientific	31985070	
MCDB 131	Thermo Fisher Scientific	10372019	
Activin A	PeproTech	120-14E	
CHIR 99021	Tocris	442310	
L-Ascorbic acid (vitamin C)	Sigma-Aldrich	A4544	
FGF7 (KGF)	PeproTech	100–19	
SANT-1	Sigma	S4572	
Retinoic acid	Sigma	R2625	
LDN-193189	Reprocell	04-0074	
TPB	EMD Millipore	565740	
ITS-X (100X)	Thermo Fisher Scientific	51500–056	
Dynabeads Protein G	Thermo Fisher Scientific	10004D	
Critical commercial assays	
Herculase II Fusion DNA Polymerase	Agilent Technologies	600679	
LongAmp Taq DNA Polymerase	NEB	M0323	
eBioscience™ Foxp3/Transcription Factor Staining Buffer Set	Thermo Fisher Scientific	00-5523-00	
Arima-Hi-C kit	Arima	A510008	
Monarch® HMW DNA Extraction Kit	NEB	T3050	
TruSeq Stranded Total RNA	Illumina	20020597	
Zero Blunt TOPO PCR Cloning Kit	Thermo Fisher Scientific	450245	
Deposited data	
Raw and processed data	This paper	GEO: GSE267330	
Raw and processed data	(Lee et al., 2019)25	GEO: GSE114102	
Raw and processed data	(Yang et al., 2022)17	GEO: GSE181480	
Raw and processed data	(Luo et al., 2023)8	4DNFIKCEWMIA, 4DNFI7RGXYFY, 4DNFIRLG5UWL, GSM6585573, GSM6585574	
Experimental models: Cell lines	
H1 dCas9-KRAB line	(Pulecio et al., 2023)67	Derived from the parental line (NIHhESC-10-0043); RRID: CVCL_WS14	
H1 iCas9 line	(Shi et al., 2017)62	Derived from the parental line (NIHhESC-10-0043); RRID: CVCL_WS14	
Oligonucleotides	
See Table S1 for screen gRNA spacer sequences, Table S2 for all other sequences	This paper	N/A	
Software and algorithms	
PRISM v10	GraphPad Software	https://www.graphpad.com/; RRID: SCR_002798	
Flowjo-v10	FlowJo LLC	https://www.flowjo.com/; RRID: SCR_008520	
R	4.4.0	https://cran.r-project.org; RRID: SCR_003005	
Genrich	0.6.1	https://github.com/jsh58/Genrich; RRID: SCR_025320	
CHOPCHOPv3		https://chopchop.cbu.uib.no/; RRID: SCR_015723	
MAGeCK	0.5.9.5	https://sourceforge.net/p/mageck/wiki/Home/; RRID: SCR_025016	
Minimap2	2.26	https://github.com/lh3/minimap2; RRID: SCR_018550	
samtools	1.2	http://htslib.org/; RRID: SCR_002105	
WhatsHap	2.2	https://whatshap.readthedocs.io/; RRID: SCR_025319	
longshot	0.4.5	https://github.com/pjedge/longshot; RRID: SCR_025318	
STAR	2.7.11b	https://github.com/alexdobin/STAR; RRID: SCR_004463	
Subread featureCounts	2.0.6	https://subread.sourceforge.net/; RRID: SCR_009803	
QuantaSoft Analysis Pro	1.0.596	https://www.bio-rad.com; RRID: SCR_025321	
HiC-Pro	3.1.0	https://github.com/nservant/HiC-Pro; RRID: SCR_017643	
DESeq2	1.40.1	https://bioconductor.org/packages/release/bioc/html/DESeq2.html; RRID: SCR_015687	
HiTC	1.34.0	http://www.bioconductor.org/packages//2.10/bioc/html/HiTC.html; RRID:SCR_013175	
HiCDCPlus	1.5.2	https://www.bioconductor.org/packages/release/bioc/html/HiCDCPlus.html; RRID: SCR_025317	
Bowtie2	2.5.1	http://bowtie-bio.sourceforge.net/bowtie2/index.shtml; RRID: SCR_016368	
MACS2	2.2.9.1	https://github.com/macs3-project/MACS; RRID: SCR_013291	
IDR	2.0.3	https://github.com/nboley/idr; RRID: SCR_017237	
IGV	2.17.4	http://www.broadinstitute.org/igv/; RRID:SCR_011793	
GATK	4.0.1.2	https://software.broadinstitute.org/gatk/; RRID:SCR_001876	
ChromBPNet	0.1.7	https://github.com/kundajelab/chrombpnet; RRID: SCR_024806	
MEME Suite (FIMO)	5.5.5	http://meme-suite.org/; RRID:SCR_001783	
bedtools	2.31.1	https://github.com/arq5x/bedtools2; RRID:SCR 006646	

Highlights

CRISPRi screen in hPSCs uncovers enhancers functionally important for pancreatic development

ONECUT1e-664kb deletion disrupts ONECUT1 induction and hampers pancreatic differentiation

Local antisense transcription and H3K27ac levels are reduced upon ONECUT1e-664kb deletion

The T2D variant rs528350911 in ONECUT1e-664kb decreases GATA4, GATA6, and FOXA2 binding

SUPPLEMENTAL INFORMATION

Supplemental information can be found online at https://doi.org/10.1016/j.celrep.2024.114640.

DECLARATION OF INTERESTS

The authors declare no competing interests.

AI-ASSISTED TECHNOLOGIES IN THE WRITING PROCESS

During the preparation of this work, the authors used ChatGPT, developed by OpenAI, to improve the clarity ofsome sentences. After using this tool, the authors reviewed and edited the content and take full responsibility for the content of the publication.
==== Refs
REFERENCES

1. Katsanis N . (2016). The continuum of causality in human genetic disorders. Genome Biol. 17 , 233. 10.1186/s13059-016-1107-9.27855690
2. Lemelman MB , Letourneau L , and Greeley SAW (2018). Neonatal Diabetes Mellitus: An Update on Diagnosis and Management. Clin. Perinatol 45 , 41–59.29406006
3. Croucha DJM , and Bodmer WF (2020). Polygenic inheritance, GWAS, polygenic risk scores, and the search for functional variants. Proceedings of the National Academy of Sciences of the United States of America 117 , 18924–18933.32753378
4. Zhang F , and Lupski JR (2015). Non-coding genetic variants in human disease. Hum. Mol. Genet 24 , R102–R110.26152199
5. Edwards SL , Beesley J , French JD , and Dunning M (2013). Beyond GWASs: Illuminating the dark road from association to function. Am. J. Hum. Genet 93 , 779–797.24210251
6. Lettice LA , Heaney SJH , Purdie LA , Li L , de Beer P , Oostra BA , Goode D , Elgar G , Hill RE , and de Graaff E (2003).A long-range Shh enhancer regulates expression in the developing limb and fin and is associated with preaxial polydactyly. Hum. Mol. Genet 12 , 1725–1735. 10.1093/hmg/ddg180.12837695
7. Claringbould A , and Zaugg JB (2021). Enhancers in disease: molecular basis and emerging treatment strategies. Trends Mol. Med 27 , 1060–1073.34420874
8. Luo R , Yan J , Oh JW , Xi W , Shigaki D , Wong W , Cho HS , Murphy D , Cutler R , Rosen BP , (2023). Dynamic network-guided CRISPRi screen identifies CTCF-loop-constrained nonlinear enhancer gene regulatory activity during cell state transitions. Nat. Genet 55 , 1336–1346. 10.1038/s41588-023-01450-7.37488417
9. Shukla A , and Huangfu D (2018). Decoding the noncoding genome via large-scale CRISPR screens. Curr. Opin. Genet. Dev 52 , 70–76.29913329
10. Schoenfelder S , and Fraser P (2019). Long-range enhancer–promoter contacts in gene expression control. Nat. Rev. Genet 20 , 437–455.31086298
11. Burgos JI , Vallier L , and Rodríguez-Seguí SA (2021). Monogenic Diabetes Modeling: In Vitro Pancreatic Differentiation From Human Pluripotent Stem Cells Gains Momentum. Front. Endocrinol 12 , 692596.
12. Miguel-Escalada I , Maestro MÁ , Balboa D , Elek A , Bernal A , Bernardo E , Grau V , García-Hurtado J , Sebé-Pedrós A , and Ferrer J (2022). Pancreas agenesis mutations disrupt a lead enhancer controlling a developmental enhancer cluster. Dev. Cell 57 , 1922–1936.e9. 10.1016/j.devcel.2022.07.014.35998583
13. Weedon MN , Cebola I , Patch A-M , Flanagan SE , De Franco E , Caswell R , Rodríguez-Seguí SA , Shaw-Smith C , Cho CHH , Allen HL , (2013). Recessive mutations in a distal PTF1A enhancer cause isolated pancreatic agenesis. Nat. Genet 46 , 61–64. 10.1038/ng.2826.24212882
14. Franco ED , Flanagan SE , Houghton JAL , Allen HL , MacKay DJG , Temple IK , Ellard S , and Hattersley AT (2015). The effect of early, comprehensive genomic testing on clinical care in neonatal diabetes: An international cohort study. Lancet 386 , 957–963. 10.1016/S0140-6736(15)60098-8.26231457
15. Gasperini M , Tome JM , and Shendure J (2020). Towards a comprehensive catalogue of validated and target-linked human enhancers. Nat. Rev. Genet 21 , 292–310.31988385
16. Yan J , and Huangfu D (2022). Epigenome rewiring in human pluripotent stem cells. Trends Cell Biol. 32 , 259–271.34955367
17. Yang D , Cho H , Tayyebi Z , Shukla A , Luo R , Dixon G , Ursu V , Stransky S , Tremmel DM , Sackett SD , (2022). CRISPR screening uncovers a central requirement for HHEX in pancreatic lineage commitment and plasticity restriction. Nat. Cell Biol 24 , 1064–1076. 10.1038/s41556-022-00946-4.35787684
18. Li QV , Dixon G , Verma N , Rosen BP , Gordillo M , Luo R , Xu C , Wang Q , Soh CL , Yang D , (2019). Genome-scale screens identify JNK–JUN signaling as a barrier for pluripotency exit and endoderm differentiation. Nat. Genet 51 , 999–1010. 10.1038/s41588-019-0408-9.31110351
19. Dixon G , Pan H , Yang D , Rosen BP , Jashari T , Verma N , Pulecio J , Caspi I , Lee K , Stransky S , (2021). QSER1 protects DNA methylation valleys from de novo methylation. Science 372 , eabd0875. 10.1126/science.abd0875.33833093
20. Xu X , Du Y , Ma L , Zhang S , Shi L , Chen Z , Zhou Z , Hui Y , Liu Y , Fang Y , (2021). Mapping germ-layer specification preventing genes in hPSCs via genome-scale CRISPR screening. iScience 24 , 101926. 10.1016/j.isci.2020.101926.33385119
21. Yilmaz A , Braverman-Gross C , Bialer-Tsypin A , Peretz M , and Benvenisty N (2020). Mapping Gene Circuits Essential for Germ Layer Differentiation via Loss-of-Function Screens in Haploid Human Embryonic Stem Cells. Cell Stem Cell 27 , 679–691.e6. 10.1016/j.stem.2020.06.023.32735778
22. Naxerova K , Stefano BD , Makofske JL , Watson EV , de Kort MA , Martin TD , Dezfulian M , Ricken D , Wooten EC , Kuroda MI , (2021). Integrated loss- And gain-of-function screens define a core network governing human embryonic stem cell behavior. Genes and Development 35 , 1527–1547. 10.1101/GAD.349048.121.34711655
23. Armendariz DA , Goetsch SC , Sundarrajan A , Sivakumar S , Wang Y , Xie S , Munshi NV , and Hon GC (2023). CHD-associated enhancers shape human cardiomyocyte lineage commitment. Elife 12 , e86206. 10.7554/eLife.86206.37096669
24. Haswell JR , Mattioli K , Gerhardinger C , Maass PG , Foster DJ , Peinado P , Wang X , Medina PP , Rinn JL , and Slack FJ (2021). Genome-wide CRISPR interference screen identifies long non-coding RNA loci required for differentiation and pluripotency. PLoS One 16 , e0252848. 10.1371/journal.pone.0252848.34731163
25. Lee K , Cho H , Rickert RW , Li QV , Pulecio J , Leslie CS , and Huangfu D (2019). FOXA2 Is Required for Enhancer Priming during Pancreatic Differentiation. Cell Rep. 28 , 382–393.e7. 10.1016/j.celrep.2019.06.034.31291575
26. Balboa D , Iworima DG , and Kieffer TJ (2021). Human Pluripotent Stem Cells to Model Islet Defects in Diabetes. Front. Endocrinol 12 .
27. Mahajan A , Taliun D , Thurner M , Robertson NR , Torres JM , Rayner NW , Payne AJ , Steinthorsdottir V , Scott RA , Grarup N , (2018). Fine-mapping type 2 diabetes loci to single-variant resolution using high-density imputation and islet-specific epigenome maps. Nat. Genet 50 , 1505–1513.30297969
28. Philippi A , Heller S , Costa IG , Senée V , Breunig M , Li Z , Kwon G , Russell R , Illing A , Lin Q , (2021). Mutations and variants of ONECUT1 in diabetes. Nat. Med 27 , 1928–1940. 10.1038/s41591-021-01502-7.34663987
29. Russ-Silsby J , Patel KA , Laver TW , Hawkes G , Johnson MB , Wakeling MN , Patil PP , Hattersley AT , Flanagan SE , Weedon MN , and De Franco E (2023). The Role of ONECUT1 Variants in Monogenic and Type 2 Diabetes Mellitus. Diabetes 72 , 1729–1734. 10.2337/db23-0498.37639628
30. Offield MF , Jetton TL , Labosky PA , Ray M , Stein RW , Magnuson MA , Hogan BL , and Wright CV (1996). PDX-1 is required for pancreatic outgrowth and differentiation of the rostral duodenum. Development 122 , 983–995. 10.1242/dev.122.3.983.8631275
31. Jonsson J , Carlsson L , Edlund T , and Edlund H (1994). Insulin-promoter-factor 1 is required for pancreas development in mice. Nature 371 , 606–609. 10.1038/371606a0.7935793
32. Zhu Z , Li QV , Lee K , Rosen BP , González F , Soh CL , and Huangfu D (2016). Genome Editing of Lineage Determinants in Human Pluripotent Stem Cells Reveals Mechanisms of Pancreatic Development and Diabetes. Cell Stem Cell 18 , 755–768. 10.1016/j.stem.2016.03.015.27133796
33. Wang X , Sterr M , Burtscher I , Böttcher A , Böttcher A , Beckenbauer J , Siehler J , Häring H-U , Häring HU , Staiger H , (2019). Point mutations in the PDX1 transactivation domain impair human β-cell development and function. Mol. Metabol 24 , 80–97.
34. Stoffers DA , Zinkin NT , Stanojevic V , Clarke WL , and Habener JF (1997). Pancreatic agenesis attributable to a single nucleotide deletion in the human IPF1 gene coding sequence. Nat. Genet 15 , 106–110. 10.1038/ng0197-106.8988180
35. Stoffers DA , Ferrer J , Clarke WL , and Habener JF (1997). Early-onset type-ll diabetes mellitus (MODY4) linked to IPF1. Nat. Genet 17 , 138–139.9326926
36. Steinthorsdottir V , Thorleifsson G , Sulem P , Helgason H , Grarup N , Sigurdsson A , Helgadottir HT , Johannsdottir H , Magnusson OT , Gudjonsson SA , (2014). Identification of low-frequency and rare sequence variants associated with elevated or reduced risk of type 2 diabetes. Nat. Genet 46 , 294–298. 10.1038/ng.2882.24464100
37. Smith SB , Qu H-Q , Taleb N , Kishimoto NY , Scheel DW , Lu Y , Patch A-M , Grabs R , Wang J , Lynn FC , (2010). Rfx6 directs islet formation and insulin production in mice and humans. Nature 463 , 775–780.20148032
38. Kropp PA , Zhu X , and Gannon M (2019). Regulation of the Pancreatic Exocrine Differentiation Program and Morphogenesis by Onecut 1/Hnf6. Cell. Mol. Gastroenterol. Hepatol 7 , 841–856, CMGH 7. 10.1016/j.jcmgh.2019.02.004.30831323
39. Kim SK , Selleri L , Lee JS , Zhang AY , Gu X , Jacobs Y , and Cleary ML (2002). Pbx1 inactivation disrupts pancreas development and in Ipf1-deficient mice promotes diabetes mellitus. Nat. Genet 30 , 430–435. 10.1038/ng860.11912494
40. Ahlgren U , Pfaff SL , Jessell TM , Edlund T , and Edlund H (1997). Independent requirement for ISL1 in formation of pancreatic mesenchyme and islet cells. Nature 385 , 257–260. 10.1038/385257a0.9000074
41. Lee CS , Sund NJ , Vatamaniuk MZ , Matschinsky FM , Stoffers DA , and Kaestner KH (2002). Foxa2 controls Pdx1 gene expression in pancreatic β-cells in vivo. Diabetes 51 , 2546–2551. 10.2337/diabetes.51.8.2546.12145169
42. Barbera JPM , Clements M , Thomas P , Rodriguez T , Meloy D , Kioussis D , and Beddington RSP (2000). The homeobox gene Hex is required in definitive endodermal tissues for normal forebrain, liver and thyroid formation. Development 127 , 2433–2445. 10.1242/dev.127.11.2433.10804184
43. Bort R , Martinez-Barbera JP , Beddington RSP , and Zaret KS (2004). Hex homeobox gene-dependent tissue positioning is required for organogenesis of the ventral pancreas. Development 131 , 797–806. 10.1242/dev.00965.14736744
44. Decker K , Goldman DC , Grasch CL , and Sussel L (2006). Gata6 is an important regulator of mouse pancreas development. Dev. Biol 298 , 415–429. 10.1016/j.ydbio.2006.06.046.16887115
45. Xuan S , Borok MJ , Decker KJ , Battle MA , Duncan SA , Hale MA , Macdonald RJ , and Sussel L (2012). Pancreas-specific deletion of mouse Gata4 and Gata6 causes pancreatic agenesis. J. Clin. Invest 122 , 3516–3528. 10.1172/JCI63352.23006325
46. Carrasco M , Delgado I , Soria B , Martín F , and Rojas A (2012). GATA4 and GATA6 control mouse pancreas organogenesis. J. Clin. Invest 122 , 3504–3515. 10.1172/JCI63240.23006330
47. Watt AJ , Zhao R , Li J , and Duncan SA (2007). Development of the mammalian liver and ventral pancreas is dependent on GATA4. BMC Dev. Biol 7 , 37. 10.1186/1471-213X-7-37.17451603
48. Spence JR , Lange AW , Lin SCJ , Kaestner KH , Lowy AM , Kim I , Whitsett JA , and Wells JM (2009). Sox17 Regulates Organ Lineage Segregation of Ventral Foregut Progenitor Cells. Dev. Cell 17 , 62–74. 10.1016/j.devcel.2009.05.012.19619492
49. Lorberbaum DS , Kishore S , Rosselot C , Sarbaugh D , Brooks EP , Aragon E , Xuan S , Simon O , Ghosh D , Mendelsohn C , (2020). Retinoic acid signaling within pancreatic endocrine progenitors regulates mouse and human β cell specification. Development 147 , dev189977. 10.1242/dev.189977.32467243
50. Haumaitre C , Barbacci E , Jenny M , Ott MO , Gradwohl G , and Cereghini S (2005). Lack of TCF2/vHNF1 in mice leads to pancreas agenesis. Proc. Natl. Acad. Sci. USA 102 , 1490–1495. 10.1073/pnas.0405776102.15668393
51. Vanhorenbeeck V , Jenny M , Cornut JF , Gradwohl G , Lemaigre FP , Rousseau GG , and Jacquemin P (2007). Role of the Onecut transcription factors in pancreas morphogenesis and in pancreatic and enteric endocrine differentiation. Dev. Biol 305 , 685–694. 10.1016/j.ydbio.2007.02.027.17400205
52. Nishimura W , Kondo T , Salameh T , El Khattabi I , Dodge R , Bonner-Weir S , and Sharma A (2006). A switch from MafB to MafA expression accompanies differentiation to pancreatic β-cells. Dev. Biol 293 , 526–539. 10.1016/j.ydbio.2006.02.028.16580660
53. Osipovich AB , Long Q , Manduchi E , Gangula R , Hipkens SB , Schneider J , Okubo T , Stoeckert CJ , Takada S , and Magnuson MA (2014). Insm1 promotes endocrine cell differentiation by modulating the expression of a network of genes that includes Neurog3 and Ripply3. Development 141 , 2939–2949. 10.1242/dev.104810.25053427
54. Nissim S , Weeks O , Talbot JC , Hedgepeth JW , Wucherpfennig J , Schatzman-Bone S , Swinburne I , Cortes M , Alexa K , Megason S , (2016). Iterative use of nuclear receptor Nr5a2 regulates multiple stages of liver and pancreas development. Dev. Biol 418 , 108–123. 10.1016/j.ydbio.2016.07.019.27474396
55. Jacquemin P , Lemaigre FP , and Rousseau GG (2003). The Onecut transcription factor HNF-6 (OC-1) is required for timely specification of the pancreas and acts upstream of Pdx-1 in the specification cascade. Dev. Biol 258 , 105–116. 10.1016/S0012-1606(03)00115-5.12781686
56. PanelApp (2024). Genomics England PanelApp. (date accessed), Monogenic diabetes (Version 2.54). https://panelapp.genomicsengland.co.uk.
57. Horikoshi M , Hara K , Ito C , Shojima N , Nagai R , Ueki K , Froguel P , and Kadowaki T (2007). Variations in the HHEX gene are associated with increased risk of type 2 diabetes in the Japanese population. Diabetologia 50 , 2461–2466. 10.1007/s00125-007-0827-5.17928989
58. Shimomura H , Sanke T , Hanabusa T , Tsunoda K , Furuta H , and Nanjo K (2000). Nonsense mutation of Islet-1 (Q310X) found in a type 2 diabetic patient with a strong family history. Diabetes 49 , 1597–1600. 10.2337/diabetes.49.9.1597.10969846
59. Iacovazzo D , Flanagan SE , Walker E , Quezado R , de Sousa Barros FA , Caswell R , Johnson MB , Wakeling M , Brändle M , Guo M , (2018). MAFA missense mutation causes familial insulinomatosis and diabetes mellitus. Proc. Natl. Acad. Sci. USA 115 , 1027–1032. 10.1073/pnas.1712262115.29339498
60. Ng NHJ , Jasmen JB , Lim CS , Lau HH , Krishnan VG , Kadiwala H , Kulkarni RN , Ræder H , Vallier L , Hoon S , and Teo AKK (2019). HNF4A Haploinsufficiency in MODY1 Abrogates Liver and Pancreas Differentiation from Patient-Derived Induced Pluripotent Stem Cells. iScience 16 , 192–205. 10.1016/j.isci.2019.05.032.31195238
61. Haldorsen IS , Vesterhus M , Ræder H , Jensen DK , Søvik O , Molven A , and Njølstad PR (2008). Lack of pancreatic body and tail in HNF1B mutation carriers. Diabet. Med 25 , 782–787. 10.1111/j.1464-5491.2008.02460.x.18644064
62. Shi ZD , Lee K , Yang D , Amin S , Verma N , Li QV , Zhu Z , Soh CL , Kumar R , Evans T , (2017). Genome Editing in hPSCs Reveals GATA6 Haploinsufficiency and a Genetic Interaction with GATA4 in Human Pancreatic Development. Cell Stem Cell 20 , 675–688.e6. 10.1016/j.stem.2017.01.001.28196600
63. Trott J , Alpagu Y , Tan EK , Shboul M , Dawood Y , Elsy M , Wollmann H , Tano V , Bonnard C , Eng S , (2020). Mitchell-Riley syndrome iPSCs exhibit reduced pancreatic endoderm differentiation due to a mutation in RFX6. Development 147 , dev194878. 10.1242/dev.194878.33033118
64. Bolt CC , Lopez-Delisle L , Hintermann A , Mascrez B , Rauseo A , Andrey G , and Duboule D (2022). Context-dependent enhancer function revealed by targeted inter-TAD relocation. Nat. Commun 13 , 3488. 10.1038/s41467-022-31241-3.35715427
65. Morgens DW , Wainberg M , Boyle EA , Ursu O , Araya CL , Tsui CK , Haney MS , Hess GT , Han K , Jeng EE , (2017). Genome-scale measurement of off-target activity using Cas9 toxicity in high-throughput screens. Nat. Commun 8 , 15178. 10.1038/ncomms15178.28474669
66. Sanson KR , Hanna RE , Hegde M , Donovan KF , Strand C , Sullender ME , Vaimberg EW , Goodale A , Root DE , Piccioni F , and Doench JG (2018). Optimized libraries for CRISPR-Cas9 genetic screens with multiple modalities. Nat. Commun 9 , 5416. 10.1038/s41467-018-07901-8.30575746
67. Pulecio J , Tayyebi Z , Liu D , Wong W , Luo R , Damodaran JR , Kaplan S , Cho H , Yan J , Murphy D , (2023). Discovery of Competent Chromatin Regions in Human Embryonic Stem Cells. Preprint at bioRxiv. 10.1101/2023.06.14.544990.
68. Li W , Xu H , Xiao T , Cong L , Love MI , Zhang F , Irizarry RA , Liu JS , Brown M , and Liu XS (2014). MAGeCK enables robust identification of essential genes from genome-scale CRISPR/Cas9 knockout screens. Genome Biol. 15 .
69. Pollard KS , Hubisz MJ , Rosenbloom KR , and Siepel A (2010). Detection of nonneutral substitution rates on mammalian phylogenies. Genome Res. 20 , 110–121. 10.1101/gr.097857.109.19858363
70. Creyghton MP , Cheng AW , Welstead GG , Kooistra T , Carey BW , Steine EJ , Hanna J , Lodato MA , Frampton GM , Sharp PA , (2010). Histone H3K27ac separates active from poised enhancers and predicts developmental state. Proceedings of the National Academy of Sciences of the United States of America 107 , 21931–21936. 10.1073/pnas.1016071107.21106759
71. Heller S , Li Z , Lin Q , Geusz R , Breunig M , Hohwieler M , Zhang X , Nair GG , Seufferlein T , Hebrok M , (2021). Transcriptional changes and the role of ONECUT1 in hPSC pancreatic differentiation. Commun. Biol 4 , 1298.34789845
72. Jacquemin P , Durviaux SM , Jensen J , Godfraind C , Gradwohl G , Guillemot F , Madsen OD , Carmeliet P , Dewerchin M , Collen D , (2000). Transcription Factor Hepatocyte Nuclear Factor 6 Regulates Pancreatic Endocrine Cell Differentiation and Controls Expression of the Proendocrine Gene ngn3. Mol. Cell Biol 20 , 4445–4454. 10.1128/mcb.20.12.4445-4454.2000.10825208
73. Pierreux CE , Poll AV , Kemp CR , Clotman F , Maestro MA , Cordi S , Ferrer J , Leyns L , Rousseau GG , and Lemaigre FP (2006). The Transcription Factor Hepatocyte Nuclear Factor-6 Controls the Development of Pancreatic Ducts in the Mouse. Gastroenterology 130 , 532–541. 10.1053/j.gastro.2005.12.005.16472605
74. Zhang H , Ables ET , Pope CF , Washington MK , Hipkens S , Means AL , Path G , Seufert J , Costa RH , Leiter AB , (2009). Multiple, temporal-specific roles for HNF6 in pancreatic endocrine and ductal differentiation. Mech. Dev 126 , 958–973. 10.1016/j.mod.2009.09.006.19766716
75. Uyehara CM , and Apostolou E (2023). 3D enhancer-promoter interactions and multi-connected hubs: Organizational principles and functional roles. Cell Rep. 42 .
76. Yang F . (2022). Promoter antisense RNAs: beyond transcription by-products of active promoters. RNA Biol. 19 , 533–540.35427206
77. Struhl K . (2007). Transcriptional noise and the fidelity of initiation by RNA polymerase II. Nat. Struct. Mol. Biol 14 , 103–105.17277804
78. Mills JD , Kawahara Y , and Janitz M (2013). Strand-Specific RNA-Seq Provides Greater Resolution of Transcriptome Profiling. Curr. Genom 14 , 173–181. 10.2174/1389202911314030003.
79. Martin BJE , Brind’Amour J , Kuzmin A , Jensen KN , Liu ZC , Lorincz M , and Howe LJ (2021). Transcription shapes genome-wide histone acetylation patterns. Nat. Commun 12 , 210.33431884
80. Slieker RC , Donnelly LA , Fitipaldi H , Bouland GA , Giordano GN , Åkerlund M , Gerl MJ , Ahlqvist E , Ali A , Dragan I , (2021). Distinct Molecular Signatures of Clinical Clusters in People With Type 2 Diabetes: An IMI-RHAPSODY Study. Diabetes 70 , 2683–2693. 10.2337/DB20-1281.34376475
81. Torres JM , Abdalla M , Payne A , Fernandez-Tajes J , Thurner M , Nylander V , Gloyn AL , Mahajan A , and McCarthy MI (2020). A Multi-omic Integrative Scheme Characterizes Tissues of Action at Loci Associated with Type 2 Diabetes. Am. J. Hum. Genet 107 , 1011–1028. 10.1016/j.ajhg.2020.10.009.33186544
82. Zou Y , Carbonetto P , Wang G , and Stephens M (2022). Fine-mapping from summary data with the “Sum of Single Effects” model. PLoS Genet. 18 , e1010299.35853082
83. Brennan KJ , Weilert M , Krueger S , Pampari A , Liu HY , Yang AWH , Morrison JA , Hughes TR , Rushlow CA , Kundaje A , and Zeitlinger J (2023). Chromatin accessibility in the Drosophila embryo is determined by transcription factor pioneering and enhancer activation. Dev. Cell 58 , 1898–1916.e9. 10.1016/j.devcel.2023.07.007.37557175
84. Pampari A , Shcherbina A , Nair S , Schreiber J , Patel A , Wang A , Kundu S , Shrikumar A , and Kundaje A (2023). Bias Factorized, Base-Resolution Deep Learning Models of Chromatin Accessibility Reveal Cis-Regulatory Sequence Syntax, Transcription Factor Footprints and Regulatory Variants.
85. Shrikumar A , Greenside P , and Kundaje A (2017). Learning Important Features through Propagating Activation Differences.
86. Shaw-Smith C , Franco ED , Allen HL , Batlle M , Flanagan SE , Borowiec M , Taplin CE , Velden JVA-VD , Cruz-Rojo J , Nanclares GPD , (2014). GATA4 mutations are a cause of neonatal and childhood-onset diabetes. Diabetes 63 , 2888–2894. 10.2337/db14-0061.24696446
87. Allen HL , Flanagan SE , Shaw-Smith C , De Franco E , Akerman I , Caswell R , International Pancreatic Agenesis Consortium; Ferrer J , Hattersley AT , and Ellard S (2011). GATA6 haploinsufficiency causes pancreatic agenesis in humans. Nat. Genet 44 , 20–22.22158542
88. Yang YP , Magnuson MA , Stein R , and Wright CVE (2017). The mammal-specific Pdx1 area II enhancer has multiple essential functions in early endocrine cell specification and postnatal β-cell maturation. Development 144 , 248–257. 10.1242/dev.143123.27993987
89. Fujitani Y , Fujitani S , Boyer DF , Gannon M , Kawaguchi Y , Ray M , Shiota M , Stein RW , Magnuson MA , and Wright CVE (2006). Targeted deletion of a cis-regulatory region reveals differential gene dosage requirements for Pdx1 in foregut organ differentiation and pancreas formation. Genes Dev. 20 , 253–266. 10.1101/gad.1360106.16418487
90. Rojas A , Schachterle W , Xu SM , and Black BL (2009). An endoderm-specific transcriptional enhancer from the mouse Gata4 gene requires GATA and homeodomain protein-binding sites for function in vivo. Dev. Dynam 238 , 2588–2598. 10.1002/dvdy.22091.
91. Kelly MR , Wisniewska K , Regner MJ , Lewis MW , Perreault AA , Davis ES , Phanstiel DH , Parker JS , and Franco HL (2022). A multi-omic dissection of super-enhancer driven oncogenic gene expression programs in ovarian cancer. Nat. Commun 13 , 4247. 10.1038/s41467-022-31919-8.35869079
92. Long HK , Osterwalder M , Welsh IC , Hansen K , Davies JOJ , Liu YE , Koska M , Adams AT , Aho R , Arora N , (2020). Loss of Extreme Long-Range Enhancers in Human Neural Crest Drives a Craniofacial Disorder. Cell Stem Cell 27 , 765–783.e14. 10.1016/j.stem.2020.09.001.32991838
93. Sagai T , Hosoya M , Mizushina Y , Tamura M , and Shiroishi T (2005). Elimination of a long-range cis-regulatory module causes complete loss of limb-specific Shh expression and truncation of the mouse limb. Development 132 , 797–803. 10.1242/dev.01613.15677727
94. Uslu VV , Petretich M , Ruf S , Langenfeld K , Fonseca NA , Marioni JC , and Spitz F (2014). Long-range enhancers regulating Myc expression are required for normal facial morphogenesis. Nat. Genet 46 , 753–758. 10.1038/ng.2971.24859337
95. Isoda T , Moore AJ , He Z , Chandra V , Aida M , Denholtz M , Piet van Hamburg J , Fisch KM , Chang AN , Fahl SP , (2017). Non-coding Transcription Instructs Chromatin Folding and Compartmentalization to Dictate Enhancer-Promoter Communication and T Cell Fate. Cell 171 , 103–119.e18. 10.1016/j.cell.2017.09.001.28938112
96. Osterwalder M , Barozzi I , Tissiéres V , Fukuda-Yuzawa Y , Mannion BJ , Afzal SY , Lee EA , Zhu Y , Plajzer-Frick I , Pickle CS , (2018). Enhancer redundancy provides phenotypic robustness in mammalian development. Nature 554 , 239–243. 10.1038/nature25461.29420474
97. Choi HI , Chai JC , Lee YS , Jung KH , and Chai YG (2021). Targeting MYC-Inducing Enhancer-Associated Noncoding (MYC-IEANC) RNAs Inhibits the Proliferation of HCC Cells.
98. Luo X , Liu Y , Dang D , Hu T , Hou Y , Meng X , Zhang F , Li T , Wang C , Li M , (2021). 3D Genome of macaque fetal brain reveals evolutionary innovations during primate corticogenesis. Cell 184 , 723–740.e21. 10.1016/j.cell.2021.01.001.33508230
99. Gonen N , Futtner CR , Wood S , Garcia-Moreno SA , Salamone IM , Samson SC , Sekido R , Poulat F , Maatouk DM , and Lovell-Badge R (2018). Sex reversal following deletion of a single distal enhancer of Sox9. Science 360 , 1469–1473. 10.1126/science.aas9408.29903884
100. Merz S , Senee V , Philippi A , Oswald F , Shaigan M , Fuehrer M , Drewes C , Allgoewer C , Oellinger R , Heni M , (2024). A ONECUT1 regulatory, non-coding region in pancreatic development and diabetes. Preprint at medRxiv. 10.1101/2024.07.23.24310605.
101. McCaw ZR , O’Dushlaine C , Somineni H , Bereket M , Klein C , Karaletsos T , Casale FP , Koller D , and Soare TW (2023). An allelic-series rare-variant association test for candidate-gene discovery. Am. J. Hum. Genet 110 , 1330–1342.37494930
102. Plenge RM , Scolnick EM , and Altshuler D (2013). Validating therapeutic targets through human genetics. Nat. Rev. Drug Discov 12 , 581–594.23868113
103. Lascar N , Brown J , Pattison H , Barnett AH , Bailey CJ , and Bellary S (2018). Type 2 diabetes in adolescents and young adults. Lancet Diabetes Endocrinol 6 , 69–80.28847479
104. Franks PW , Pearson E , and Florez JC (2013). Gene-environment and gene-treatment interactions in type 2 diabetes: Progress, pitfalls, and prospects. Diabetes Care 36 , 1413–1421.23613601
105. Rezania A , Bruin JE , Arora P , Rubin A , Batushansky I , Asadi A , O’Dwyer S , Quiskamp N , Mojibian M , Albrecht T , (2014). Reversal of diabetes with insulin-producing cells derived in vitro from human pluripotent stem cells. Nat. Biotechnol 32 , 1121–1133. 10.1038/nbt.3033.25211370
106. Gaspar JM (2018). Genrich: Detecting Sites of Genomic Enrichment.
107. Labun K , Montague TG , Krause M , Torres Cleuren YN , Tjeldnes F , and Valen E (2019). CHOPCHOP v3: Expanding the CRISPR web toolbox beyond genome editing. Nucleic Acids Res. 47 , W171–W174. 10.1093/nar/gkz365.31106371
108. Sanjana NE , Shalem O , and Zhang F (2014). Improved vectors and genome-wide libraries for CRISPR screening. Nat. Methods 11 , 783–784.25075903
109. Li H (2018). Minimap2: pairwise alignment for nucleotide sequences. Bioinformatics 34 , 3094–3100.29750242
110. Danecek P , Bonfield JK , Liddle J , Marshall J , Ohan V , Pollard MO , Whitwham A , Keane T , McCarthy SA , Davies RM , and Li H (2021). Twelve years of SAMtools and BCFtools. GigaScience 10 , giab008. 10.1093/gigascience/giab008.33590861
111. Martin M , Patterson M , Garg S , Fischer SO , Pisanti N , Klau GW , Schöenhuth A , and Marschall T (2016). WhatsHap: fast and accurate read-based phasing. Preprint at bioRxiv. 10.1101/085050.
112. Edge P , and Bansal V (2019). Longshot enables accurate variant calling in diploid genomes from single-molecule long read sequencing. Nat. Commun 10 , 4660. 10.1038/s41467-019-12493-y.31604920
113. Liao Y , Smyth GK , and Shi W (2014). FeatureCounts: An efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics 30 , 923–930. 10.1093/bioinformatics/btt656.24227677
114. Servant N , Varoquaux N , Lajoie BR , Viara E , Chen CJ , Vert JP , Heard E , Dekker J , and Barillot E (2015). HiC-Pro: An optimized and flexible pipeline for Hi-C data processing. Genome Biol. 16 , 259. 10.1186/s13059-015-0831-x.26619908
115. Kramer NE , Davis ES , Wenger CD , Deoudes EM , Parker SM , Love MI , and Phanstiel DH (2022). Plotgardener: Cultivating precise multi-panel figures in R. Bioinformatics 38 , 2042–2045. 10.1093/bioinformatics/btac057.35134826
116. Love MI , Huber W , and Anders S (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15 , 550. 10.1186/s13059-014-0550-8.25516281
117. Servant N , Lajoie BR , Nora EP , Giorgetti L , Chen CJ , Heard E , Dekker J , and Barillot E (2012). HiTC: Exploration of high-throughput ‘C’ experiments. Bioinformatics 28 , 2843–2844. 10.1093/bioinformatics/bts521.22923296
118. Sahin M , Wong W , Zhan Y , Van Deynze K , Koche R , and Leslie CS (2021). HiC-DC+ enables systematic 3D interaction calls and differential analysis for Hi-C and HiChIP. Nat. Commun 12 , 3366. 10.1038/s41467-021-23749-x.34099725
119. Langmead B , Trapnell C , Pop M , and Salzberg SL (2009). Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 10 , R25. 10.1186/gb-2009-10-3-r25.19261174
120. Zhang Y , Liu T , Meyer CA , Eeckhoute J , Johnson DS , Bernstein BE , Nusbaum C , Myers RM , Brown M , Li W , and Liu XS (2008). Model-based analysis of ChIP-Seq (MACS). Genome Biol. 9 , R137–R139.18798982
121. Li Q , Brown JB , Huang H , and Bickel PJ (2011). Measuring reproducibility of high-throughput experiments. Ann. Appl. Stat 5 , 1752–1779. 10.1214/11-AOAS466.
122. Amemiya HM , Kundaje A , and Boyle AP (2019). The ENCODE Blacklist: Identification of Problematic Regions of the Genome. Sci. Rep 9 , 9354. 10.1038/s41598-019-45839-z.31249361
123. Robinson JT , Thorvaldsdóttir H , Winckler W , Guttman M , Lander ES , Getz G , and Mesirov JP (2011). Integrative genomics viewer. Nat. Biotechnol 9 .
124. Van der Auwera GA , and O’Connor BD (2020). Genomics in the Cloud: Using Docker, GATK, and WDL in Terra (O’Reilly Media).
125. Grant CE , Bailey TL , and Noble WS (2011). FIMO: Scanning for occurrences of a given motif. Bioinformatics 27 , 1017–1018. 10.1093/bioinformatics/btr064.21330290
126. Weirauch MT , Yang A , Albu M , Cote AG , Montenegro-Montero A , Drewe P , Najafabadi HS , Lambert SA , Mann I , Cook K , (2014). Determination and inference of eukaryotic transcription factor sequence specificity. Cell 158 , 1431–1443. 10.1016/j.cell.2014.08.009.25215497
