A fundamental question in biology is how a single genome gives rise to the great diversity of cell types that make up organs and tissues. A key goal is to map all cell types of developing and mature organs such as the pancreas, an essential organ at the basis of multiple human disorders including diabetes and cancer (Kahn, Cooper and Del Prato, 2014; Kamisawa et al., 2016; Han et al., 2020) . Single-cell RNA sequencing (scRNA-seq) provides a powerful tool to resolve cellular heterogeneity, identify cell types and capture highresolution snapshots of gene expression in individual cells (Grün et al., 2016) . With the advent of singlecell transcriptomics, great progress has been made toward the creation of a reference cell atlas of the pancreas (Baron et al., 2016; Segerstolpe et al., 2016; Wang et al., 2016; Xin et al., 2016; Enge et al., 2017; Lawlor et al., 2017; Augsornworawat and Millman, 2020; Han et al., 2020) . Work from several groups provided cellular atlases of the pancreas during mouse development (Stanescu et al., 2017; Byrnes et al., 2018; Scavuzzo et al., 2018) , in adult mice (Muraro et al., 2016) and in human fetal and adult pancreas (Liu et al., 2014; Baron et al., 2016; Muraro et al., 2016; Segerstolpe et al., 2016; Wang et al., 2016; Enge et al., 2017; Han et al., 2020) . Efforts have also been made to map cellular identity during pancreas development starting from human pluripotent stem cells Hrvatin et al., 2014; Zhu et al., 2016; Han et al., 2020; Hogrebe et al., 2020; Peterson et al., 2020) . Taken together, these studies provide an opportunity to better understand the maintenance and establishment of cellular identity among different pancreatic cell types.
Work over the past decades indicated that cellular identity is established by combinations of transcription factors (TFs) that recognize and interact with cis-regulatory elements in the genome. These transcription factors, together with chromatin modifiers, give rise to gene expression programs. A small number of core TFs are thought to be sufficient for the establishment and maintenance of gene expression programs that define cellular identity during and after development (Ohno, 1979) . Studies conducted in both mouse and human have successfully identified TFs that are pivotal for the acquisition and maintenance of pancreatic cell fates (Dassaye, Naidoo and Cerf, 2016) . These include PDX1 (Zhou et al., 2008; Shih et al., 2015) , MAFA (Nishimura, Takahashi and Yasuda, 2015) , NGN3 (Gradwohl et al., 2000) , NKX2.2 (Sussel et al., 1998) , PAX4 (Sosa-Pineda et al., 1997) , NKX6.1, NEUROD1 (Mastracci et al., 2013) , ARX (Collombat et al., 2003) , MAFB (Artner et al., 2006) , RFX6 (Smith et al., 2010) , GATA4 (Ketola et al., 2004) , FOXA2 and SOX9 (Shroff et al., 2014; Shih et al., 2015) . Conditional deletion of TFs such as FOXA2 and PDX1 in adult beta cells results in the loss of cellular identity and function (Sund et al., 2001; Gao et al., 2014) . Genetic evidence for the role of these TFs in establishing human pancreatic cell identity is provided by the identification of TF loss-of-function mutations that cause pancreatic agenesis (Stoffers et al., 1997; Sellick et al., 2004; Allen et al., 2012; Shaw-Smith et al., 2014; De Franco et al., 2019) or neonatal or young-onset diabetes (Senée et al., 2006; Solomon et al., 2009; Rubio-Cabezas et al., 2010 Smith et al., 2010; Bonnefond et al., 2013; Flanagan et al., 2014) . In addition, TF overexpression can reprogram somatic cells to adopt alternative identities (Takahashi and Yamanaka, 2006; Zhou et al., 2008; Vierbuchen et al., 2010; Lima et al., 2016) . For example, the induced expression of Ngn3, Pdx1, and Mafa was shown to reprogram mouse alpha cells into beta-like cells in vivo (Zhou et al., 2008) . However, how key TFs underlie the maintenance of cellular identity in the human pancreas remains incompletely understood.
Over the past decade, multiple approaches to reconstruct gene regulatory networks (GRNs) from bulk and single-cell omics data have been developed (Ghazanfar et al., 2016; Lim et al., 2016; Matsumoto et al., 2017; Fiers et al., 2018) . In particular, it is now possible to combine single-cell transcriptomic data with cisregulatory information to infer GRNs (Janky et al., 2014; Aibar et al., 2017; Van de Sande et al., 2020) .
Because TFs recognize DNA motifs in the genome, one can measure if inferred target genes are expressed within single cells, and therefore quantify the activity of TFs. Such approaches have revealed the regulatory programs in distinct systems including the Drosophila brain (Davie et al , 2018) , cancer (Wouters et al., 2020) , during early mouse embryonic development (Peng et al , 2019) , in a mouse cell atlas (Suo et al., 2018 ) and a human cell atlas (Han et al., 2020) .
Analysis of GRNs in the human adult pancreas has identified distinct endocrine and exocrine regulatory states with multiple stable cell states for alpha, beta and ductal cells (Kumar and Vinod, 2019) . Type 2 diabetes and body mass index (BMI) were shown not to impact GRN activity of alpha and beta cells (Kumar and Vinod, 2019) . Previous data shows that type 2 diabetic (Faerch et al., 2015; Dennis et al., 2019; Dybala and Hara, 2019; Zaharia et al., 2019) and non-diabetic human islet preparations vary greatly depending on age (Enge et al., 2017; Westacott et al., 2017) and BMI (Henquin, 2018) warranting the exploration of GRNs in larger cohorts. Hence, it remains unclear whether previous GRN findings can be extrapolated to a broader, highly heterogeneous non-diabetic and type 2 diabetes patient population. The development of integration methods provides an opportunity to analyze multiple scRNA-seq studies from multiple laboratories and patients (Butler et al., 2018; Luecken et al., 2020) . Additional knowledge on how GRNs maintain cellular identity in the human adult pancreas may further the understanding of disease states as well as improve efforts to convert patient cells into functional, mature beta cells for diabetes treatment.
Here, we build an integrated human pancreas gene regulatory atlas. In this resource, we use single-cell transcriptomes of the human adult pancreas, taking advantage of integration strategies and computational tools to reconstruct GRNs. Our analysis identifies the GRN landscape and candidate regulators that are critical for cellular identity in the human adult pancreas.
Integrating multiple human adult pancreas scRNA-seq datasets can further improve the power of scRNAseq analyses to create a human adult pancreas cell atlas. We set out to analyse and integrate five publicly available datasets covering a total of 35 non-diabetic, 15 type 2 diabetic and 1 type 1 diabetic individuals using Seurat v3.0 CCA integration tools ( Figure 1A , detailed donor information can be found in Table S1 ) (Segerstolpe et al., 2016; Wang et al., 2016; Xin et al., 2016; Enge et al., 2017; Lawlor et al., 2017; Hafemeister and Satija, 2019; Stuart et al., 2019) . After filtering out low quality transcriptomes and data integration, uniform manifold approximation and projection for dimension reduction (UMAP) visualisation revealed that 7393 cells localize into distinct clusters ( Figure 1B) . Cells from each original dataset localize together suggesting that the location of cells on the UMAP is not driven by the dataset of origin ( Figure 1C ).
We next sought to identify pancreatic cell types ( Figure 1D ). Clustering analyses based on the expression of well-established cell type specific markers led to the identification of eight cell types in the human adult pancreas: beta, alpha, gamma, delta, acinar, ductal, stellate and endothelial cells ( Figure 1B , Table S2 ).
UMAP visualization allowed for the segregation of endocrine, exocrine and other lineages ( Figure S1A ).
Beta cells grouped together, away from other clusters and were marked by INS expression (Figure 1E ).
Other distinct clusters corresponded to alpha, gamma and delta cells based on global transcriptional similarity and elevated expression of GCG, PPY and SST, respectively, and other markers ( Figure 1E , Table S2 ). Using a similar approach, we reliably detected other, previously described, major pancreatic cell types (acinar, ductal, endothelial and stellate, Figure 1B ). All cell types were detected in both non-diabetic and type 2 diabetic pancreases ( Figure 1F ). An additional four rare cell populations, that cannot be robustly identified through clustering analyses, were identified manually by assessing the expression of GHRL (epsilon cells), TPS1AB (schwann cells), CD86 (mast cells) and SOX10 (Major histocompatibility complex (MHC) class 2 cells) (Wierup et al., 2002; Segerstolpe et al., 2016) (Figure 1G ). These rare cell types often cluster with other common cell types. Importantly, our annotation recapitulated previous annotations to a large extent (Figure S1B-C). In summary, we reconstructed an integrated single-cell atlas of the human adult pancreas, and annotated 12 pancreatic cell types.
Next, we set out to comprehensively reconstruct GRNs for all pancreatic cell types from single-cell transcriptomic data, applying single-cell regulatory network inference and clustering (pySCENIC) (Aibar et al., 2017; Van de Sande et al., 2020) . PySCENIC links cis-regulatory sequence information together with single-cell transcriptomes in three sequential steps by 1) co-expression analysis, 2) target gene motif and ChIP-seq track enrichment analysis, and 3) regulon activity evaluation (Figure 2A) . Each regulon consists of a TF with its predicted target genes (co-expressed genes with an enriched TF motif), altogether forming a regulon. pySCENIC identified 142 regulons that characterize the GRNs of the human adult pancreas 5 ( Figure 2B /C, Table S3 ). Multiple regulons identified here as active in the pancreas correspond to TF binding motifs enriched in accessible chromatin in the pancreas (assessed by ATAC-seq in FACS-purified pancreatic cells, (Arda et al., 2018) ), supporting the validity of the approach ( Figure S2A ).
UMAP visualization based on the activity of 142 regulons in non-diabetic and type 2 diabetic pancreata revealed groups of cells that differ from one another based on their regulatory activity ( Figure 2B -C). In particular, there are distinct regulatory states for exocrine and endocrine pancreatic lineages, stellate and endothelial cells (Figure 2B /C). Endocrine cell types clustered together, indicating shared regulatory states, while exocrine cell types formed two distinct clusters. Stellate and endothelial cells differed most from other cell types in their regulatory states. These results are consistent with previous analyses (Baron, Veres, Samuel L. Wolock, et al., 2016; Lawlor et al., 2017; Kumar and Vinod, 2019) and are also in line with our findings above based on gene expression analysis ( Figure 1D ). As expected, regulons active in endocrine cell types include RFX6, PAX6 and NEUROD1 ( Figure 2D -G). These TFs have reported roles in endocrine cell fate commitment and maintenance of cell identity throughout adult life (Smith et al., 2010; Hart et al., 2013; Mastracci et al., 2013) . Using iRegulon for visualization, many of the NEUROD1 target genes identified here have been previously linked to beta cell survival and function, such as SNAP25, TSPAN2, (Gierl et al., 2006; Haumaitre et al., 2008; Hart et al., 2013; Churchill et al., 2017) . Clustering all cells based on the activity of all regulons identifies regulatory modules ( Figure 2D , red squares). In the exocrine pancreas, one regulatory module, containing NR5A2, was shared between acinar and ductal cells, although with a tendency for increased regulon activity in ductal cells ( Figure 2D /H). Other exocrine regulons included ONECUT1, REST and HNF1B with reported roles in exocrine development (Nissim et al., 2016; Kropp, Zhu and Gannon, 2019) and the adult exocrine pancreas (Quilichini et al., 2019; Bray et al., 2020 ) ( Figure 2D ). In summary, this analysis confirms the expected separation of exocrine and endocrine cells with distinct gene regulatory programs, and identifies known and novel candidate regulators of pancreas cell states.
Several regulatory modules are shared between different cell types within the endocrine and exocrine pancreas. Additionally, each cell type is defined by cell-type specific regulatory modules ( Figure 1D ). In the endocrine pancreas, alpha and beta cells shared endocrine regulons (MAFB, MEIS2), whereas we observed distinct activities for ARX and IRX2 regulons in alpha cells and RXRG and PDX1 in beta cells ( Figure 2D /H), expanding previous findings (Kumar and Vinod, 2019) . Using iRegulon for visualization, PDX1 target genes include SLC6A17, PDIA6 and ABHD3, which have been reported to control insulin release (Thomas, Brown and Brown, 2014; Eletto et al., 2016; Rorsman and Ashcroft, 2018 ) ( Figure S2C ).
Interestingly, gamma and delta cells overlapped with alpha and beta cells, respectively, suggesting a shared regulatory state ( Figure 2C/D) . This includes shared regulon activity for ARX in gamma and alpha and PDX1 in beta and delta cells ( Figure 2D /H), consistent with their reported expression in published scRNA-seq studies (Baron et al., 2016; Segerstolpe et al., 2016; Lawlor et al., 2017) . GATA4 and RBPJL, known acinar-specific TFs, were highly active in acinar cells (Masui et al., 2010; Carrasco et al., 2012) ( Figure 2H) . Similarly, ductal cells were characterised by highly active SOX9 and POU2F3 regulons, in line with previous literature (Shroff et al., 2014; Yamashita et al., 2017) (Figure 2H ). In sum, this analysis confirms that alpha, beta, acinar and ductal cells are defined by the activity of distinct TF combinations that form gene regulatory modules.
In conclusion, the network approach recovers many of the expected regulators of pancreatic cellular identity allowing for the comprehensive characterisation of the gene regulatory state of all major human adult pancreatic cell types.
A comprehensive network analysis provides an opportunity to predict and identify critical regulators of cell identity. To identify regulons with highly cell type-specific activities within the human adult non-diabetic pancreas, we calculated regulon specificity scores (RSS) (a complete list of RSSs can be found in Table S4 ) (Suo et al., 2018) . The RSS utilises Jensen-Shannon divergence to measure the similarity between the probability distribution of the regulon's enrichment score and cell type annotation wherein outliers receive a higher RSS and are therefore considered cell type-specific (Suo et al., 2018) . It can therefore be used to rank the activity of TFs within specific cell types.
Among the top regulons identified in alpha cells, we recover well known regulators of alpha and endocrine cell fate such as ARX, IRX2, PAX6, MAFB, NEUROD1 and RFX6 ( Figure 3A /B) (Collombat et al., 2003; Artner et al., 2006; Delporte et al., 2008; Smith et al., 2010; Dorrell et al., 2011; Mastracci et al., 2013) . In addition, we identified JUND, EGR4, SREBF1 and STAT4, which have not yet been implicated in alpha cell identity. EGR1 (but not EGR4) has been shown to transcriptionally regulate glucagon expression (Leung- Theung-Long et al., 2005) as well as the PDX1 promoter in beta cells (Eto, Kaur and Thomas, 2007) . STAT4 and JUND have been described in pancreatic tissue in general and beta cells, respectively, but not in alpha cells (Yu and Kim, 2012; Good et al., 2019) (Figure 3B /C). These TFs respond to the JNK and EGFR signalling pathways and may have important physiological functions. Both JUND and the JUND/JNK signaling pathway have been implicated in pancreatic cancer (Shin et al., 2009; Recio-Boiles et al., 2016) .
Immunocytochemistry of the human adult pancreas confirmed the presence of nuclear JUND in islets ( Figure 3Di ). We also detected nuclear JUND protein in a subset of human induced pluripotent stem cells (iPSCs) subjected to beta cell differentiation (Figure 3E/F) . Surprisingly, we also detected JUND protein in ductal cells, despite lower JUND regulon activity in this cell type (Figures 3C, 3Dii) . Thus, protein 7 expression does not always predict regulatory activity. Nevertheless, these results show that JUND is present and active in a subset of pancreatic cell types in the human adult pancreas and human pluripotent stem cell derived islet cells. Altogether, this analysis predicts TFs active in human alpha cells, recovering known as well as new candidate TFs.
Among the top regulons identified in beta cells, we retrieved well-known as well as new candidate regulators of beta and endocrine cell identity. Known TFs include RXRG, PDX1, NEUROD1, PAX6 and RFX6 (Zhou et al., 2008; Miyazaki et al., 2010; Smith et al., 2010; Mastracci and Sussel, 2012; Hart et al., 2013) (Figure   3G /H). In addition, we found that ZNF705D, ASCL2 and HOXD13 were highly ranked regulons ( Figure   3H /I). HOXD13 and BHLHE41 have been shown to be present in the exocrine pancreas (Cantile et al., 2009; Sato et al., 2012) . Interestingly, ASCL2 has been reported to interact with β-catenin of the Wnt pathway, the latter has an established role in endocrine fate specification during in vitro differentiation (Schuijers et al., 2015; Sharon et al., 2019; Vethe et al., 2019) . Many putative target genes of ASCL2 including PDX1, INS, ABCC8, FOXA1, KCNK16, FXYD2 are directly related to glucose sensing and beta cell identity, in line with the beta cell-specific regulatory activity of ASCL2 ( Figure S3A ) (Gao et al., 2010; Arystarkhova et al., 2013; Vierra et al., 2015; Park, Lee and Park, 2016) . FXYD2γa, a regulatory subunit of the Na + -K + -ATPase, is a transcript exclusively expressed in human beta cells (Flamez et al., 2010) .
Immunohistochemistry of human adult pancreas sections showed that ASCL2 is expressed in INS + beta and islet cells ( Figure 3J) . Surprisingly, ASCL2 was mainly localized to the cytoplasm ( Figure 3J) , which is unexpected for TFs which tend to localize to the nucleus (Baranek, Sock and Wegner, 2005) . Cytoplasmic localization of ASCL2 has been reported in the context of colon and breast cancer (Zhu et al., 2012; Xu et al., 2017) . These results implicate additional TFs including ASCL2 in the regulation of beta cell identity.
They also illustrate the value of network analyses to increase our understanding of the biology of the human pancreas.
In summary, GRN analysis and regulon ranking allowed us to pinpoint both known and novel candidate regulators of pancreatic endocrine cell identity, providing a resource for further investigation of their roles in cellular identity and function.
Similarly, the comprehensive network analysis provides an opportunity to predict and identify critical regulators of exocrine cell identity.
We also identify known and new TFs in acinar cells. Among the top acinar-specific regulons, we recovered well known regulators of acinar and exocrine cell identity such as PTF1A, RBPJL, GATA4 and NR5A2 (Ketola et al., 2004; Masui et al., 2010; Nissim et al., 2016; Sakikubo et al., 2018) (Figure 4A/B) . These findings are in line with a recent study that used single-nucleus RNA-seq on pancreatic acinar tissue (Tosti 8 et al., 2019) . Furthermore, we identified MECOM, HEYL and TGIF1 as highly ranked regulons ( Figure   4B /C). Interestingly, aberrant Mecom expression has been linked to the induction of gastric genes in acinar cells, which disrupts acinar cell identity and increases susceptibility to malignancy (Hoang et al., 2016) . The loss of TGIF1 has been linked to pancreatic ductal adenocarcinoma progression making further exploration of these regulons interesting in the context of cancer biology (Weng et al., 2019) . Ectopic expression of Tgif2 (but not Tgif1) reprograms mouse liver cells towards a pancreas progenitor state (Cerdá-Esteban et al., 2017) . HEYL is a reported Notch signalling target gene in NGN3 + exocrine cells (Gomez et al, 2015) .
We confirmed nuclear expression of HEYL in human acinar and islet cells (Figure 4Di /ii, detailed donor information can be found in Table S1 ) by immunohistochemistry, in agreement with elevated HEYL regulon activity in acinar cells ( Figure 3C ). further functional studies in the healthy pancreatic context (Matsumoto et al., 2004; Coleman et al., 2013) .
In summary, GRN analysis and regulon ranking allowed us to pinpoint both known and novel candidate regulators of pancreatic exocrine cell identity. Specifically, we identified HEYL as a candidate TF that might be important for acinar cell identity, warranting further investigation.
To enable users to easily navigate the human pancreatic cell network atlas, we provide a loom file that allows for the visualisation and exploration of the data using the web-based portal SCope (Davie et al., 2018) (.loom file and tutorial available at http://scope.aertslab.org/#/PancreasAtlas/*/welcome and https://github.com/pasquelab/scPancreasAtlas). Features such as cell type annotation as defined in this paper, gene expression and regulon activity can be explored on the regulon and gene expression based UMAP. This resource enables users to select and visualize up to three genes or regulons simultaneously and select subsets of cells for downstream analyses. For example, the expression of COVID-19 related genes can be interactively explored (Yang et al., 2020) . Target genes of a specific regulon can be downloaded to facilitate further exploration, for example in iRegulon or gene ontology analysis (Janky et al., 2014) . A list of predicted target genes of all 142 regulons can also be found in Table S5 . Furthermore, a list of target genes can be manually defined to compute the activity of a custom regulon. This resource can be used to further study cell identity and gene regulation in the context of the pancreas, diabetes and cancer.
In this resource, we take advantage of integration strategies and new computational tools to reconstruct an integrated cell and GRN atlas of the human adult pancreas from single-cell transcriptome data. This approach provides a comprehensive analysis of the gene regulatory logic underlying cellular identity in the human adult pancreas in a broad range of individuals, limiting the influence of inter-donor variability. We recovered known regulators of pancreatic cell identity and uncovered novel candidate regulators of cell identity that can be further investigated for their roles in cellular identity and function. By validating regulon analyses and creating an easily accessible interactive online resource which allows for the exploration of the gene regulatory state of 7393 cells from 51 individuals, this approach extends beyond previous gene regulatory studies in the human adult pancreas (Augsornworawat and Millman, 2020) .
The present analysis identified regulators of pancreatic development, function and survival that are known to be critical in humans because loss-of-gene function causes pancreatic agenesis or young onset diabetes.
For example, PTF1A (Sellick et al., 2004; Weedon et al., 2014) and GATA4 (Shaw-Smith et al., 2014) , whose loss of function are linked to pancreatic agenesis and neonatal diabetes, were among the top acinarspecific regulons (Figure 3 ). In addition, monogenic diabetes related genes PDX1 (Nicolino et al., 2010 ), NEUROD1 (Rubio-Cabezas et al., 2010 , PAX6 (Solomon et al., 2009) , RFX6 (Smith et al., 2010; Patel et al., 2017) and GLIS3 (Senée et al., 2006) were among the top beta cell-specific regulons (Figure 3 and Table S4 ).
stress signalling (Table S4) . CREB3 and CREB3L2 are non-canonical ER stress transducers that are induced in human islets and clonal beta cells upon exposure to the saturated fatty acid palmitate (Cnop et al., 2014) . Interestingly, SREBF1 and -2 undergo similar ER exit and proteolytic processing in the Golgi as these ER stress transducers, but they do so in response to changes in ER cholesterol content; both also have high regulon activity in alpha and beta cells. XBP1 is abundantly expressed in the exocrine and endocrine pancreas (Cnop et al., 2017) , but the XBP1 regulon has its highest specificity in beta cells. ATF3 and ATF4 are TFs that are activated upon eIF2α phosphorylation, an ER stress response pathway to which no less than 5 monogenic forms of diabetes belong (Eizirik, Pasquali and Cnop, 2020) . Our data underscores the importance of these TFs for endocrine pancreatic cell identity.
Given that we predict novel regulators of cell identity in the human pancreas, it will be interesting to also expand this analysis to pancreas embryonic development. Our work may also be beneficial in guiding improvements of and better understanding the in vitro derivation of pancreatic cell types. For example, the emergence of SST-positive cells together with beta-like cells at the end of in vitro differentiation could be explained by the overlap in regulatory states between beta and delta cells (Baron et al., 2016) . GRN analyses are particularly interesting for in vitro derived beta cells since a better understanding of the regulatory logic underlying control of beta cell fate may improve or facilitate future applications in regenerative medicine (Pagliuca et al., 2014; Rezania et al., 2014; Nostro et al., 2015; Russ et al., 2015; Baeyens et al., 2018) . Alternatively, many observed GRNs such as ASCL2, MECOM, PPARD, GATA6 and CDX2 are linked to pancreatic cancer making the additional exploration of GRNs interesting in the context of cancer biology (Matsumoto et al., 2004; Zhu et al., 2012; Coleman et al., 2013; Hoang et al., 2016; Xu et al., 2017; Weng et al., 2019; Brunton et al., 2020) . Finally, recent reports have stratified type 2 diabetes patients based on age at diagnosis, BMI, HbA1c and insulin secretion and sensitivity, and identified subtypes with different genetic predisposition, treatment response, disease progression and complication rates (Ahlqvist et al., 2018) . Hence, it would be interesting to assess differences in gene regulatory state and gene expression profiles of alpha and beta cells between different type 2 diabetic subgroups.
It is important to note that pySCENIC is a stochastic algorithm that does not produce precisely the same regulons for repeated applications, limiting reproducibility when comparing different datasets (Huynh-Thu et al., 2010; Van de Sande et al., 2020) . To mitigate this uncertainty, we ran the full pySCENIC pipeline five times and only kept consistent regulons with the highest regulon activity. The performance of pySCENIC, and other GRN inference methods, suffers due to the large amount of drop-out events in scRNA-seq data warranting caution when interpreting results (Chen and Mar, 2018) . This could explain the absence of wellestablished pancreas TFs such as MAFA (Olbrot et al., 2002) , MNX1 (Flanagan et al., 2014) , NEUROG3 (Krentz et al., 2017) , FOXA2 (Lee et al., 2002) and NKX2-2 (Mastracci et al., 2011) in this analysis.
Nevertheless, in support of the validity of our findings, ATAC-seq, literature and immunohistochemistry of human pancreas sections corroborate several pySCENIC predictions. Chen and colleagues underline the importance of using large sample sizes to derive the most accurate network inference possible (Chen and Mar, 2018) , highlighting the importance of dataset integration to increase the number of cells analysed. In the future, it will be interesting to extend our analyses to include many more cells and patients. In spite of current caveats, GRN analysis has enabled the capture of biological relevant information (Butte et al., 2000) .
One additional limitation of this study is the assumption that all TFs bind their binding motifs in the promoters of expressed genes. However, TF binding can be restricted to a subset of TF motifs in the genome due to influence of chromatin processes including the presence of nucleosomes as well as DNA methylation.
Therefore, additional approaches such as single cell multi-omics that capture additional layers of genome regulation will be helpful to increase our understanding of gene regulation in the context of the human pancreas.
Taken together, our GRN atlas, containing 51 individuals, provides a valuable resource for future studies in the human pancreas development, donor variability, homeostasis and disease including type 2 diabetes 11 and pancreatic cancer. Finally, our results provide new insights into the activity of TFs and gene regulation in the human adult pancreas from a gene regulatory perspective. questions about data analysis should be directed to the Lead Contact, Vincent Pasque (vincent.pasque@kuleuven.be). The reviewer tokens for this GEO repository is cfsjciumfferjkd.
Motif discovery of bulk ATAC-seq data Paired-end raw reads for bulk ATAC-seq (see Key Resource Table) were downloaded from SRA using SRA toolkit (v2.9.4). Reads were aligned and further analyzed using the ENCODE ATAC-seq pipeline with default parameters using the ENCODE human reference genome GRCh38.15 (Lee, 2016) . Bed files
containing the global open chromatin landscape of Adult Alpha (Alpha_1; EA4 and Alpha_2; EA28), Beta (Beta_1; EA5 and Beta_2; EA29), Acinar (Acinar_1; EA7 and Acinar_2; EA27) and ductal (EA11) cells or cell type specific differentially accessible regions were used as input for motif discovery by HOMER (v4.10.4) using the 'findMotifsGenome.pl' with options using hg38 with size given (Heinz et al., 2010) . The TFs whose motifs identified by HOMER correspond with TFs identified by pySCENIC are visualized in Figure S2A .
Analysis of publicly available scRNA-seq data Raw reads for five publicly available scRNA-seq datasets (see Key Resource Table) were downloaded from SRA using SRA toolkit (v2.9.4). Afterwards, reads were aligned to the human reference genome
GRCh38.95 using STAR (v2.5.3a) with default parameters followed by the conversion to the coordinate sorted BAM format. Next, the featureCounts command from the "Rsubread" (v1.5.2) package in R (v3.6.1) was used to assign mapped reads to genomic features. Low quality transcriptomes with a mitochondrial contamination greater than 5% and less than 200 expressed genes per cell were excluded from subsequent analyses. The resulting raw count matrix was batch corrected using the FindIntegrationAnchors and IntegrateData functions from the "Seurat" package (v3.1.1) after which subsequent analyses were carried out in the R package "Seurat" (v3.1.4). Gene expression was used to cluster all 7393 cells with UMAP, using Seurat's function RunUMAP.
Clusters for cell type annotation were defined using Seurat's shared nearest neighbour algorithm FindClusters function after which differential expression analysis was performed using Wilcoxon's rank sum test with a minimum cutoff of 0.25 average log fold change and min.pct of 0.25.
pySCENIC GRNs were inferred using pySCENIC (python implementation of SCENIC, v0.9.15) in Python version 3.6.9 (Aibar et al., 2017) . Integrated read counts were used as input to run GENIE3 (Huynh-Thu et al., 2010) which is part of arboreto (v0.1.5). GRNs were subsequently inferred using pySCENIC with the hg38_refseq-r80 motif database and default settings. To control for the stochasticity, which is inherent to pySCENIC, a consensus GRN was generated by merging results from five repeat pySCENIC runs. If regulons were identified in multiple pySCENIC runs, only the regulon with the highest AUC value was retained. Regulon activity represented by AUCell values was used to cluster all 7393 cells with UMAP, using Seurat's RunUMAP function.
All 142 regulons within non-diabetic cell types were visualized using the 'clustermap' function of the Python package "seaborn" (v0.9.0). The z-score for each regulon across all cells was calculated using the z-score parameter of the seaborn 'clustermap' function.
Extended analysis of the target genes of specific regulons was conducted in Cytoscape (v3.7.1) using the iRegulon application (v1.3). The list of target genes of a specific regulon was downloaded from the loom file through the SCope platform (https://github.com/pasquelab/scPancreasAtlas) (Davie et al., 2018) .
To quantify the cell-type specificity of a regulon, we utilized an entropy-based strategy as described previously (Suo et al., 2018) using the AUCell matrix as input in MATLAB R2019b. The top 10 most specific regulons were subsequently visualized using the R package ggplot2 (v3.1.1). The complete regulon ranking list is available in Table S2 .
Control iPSC line HEL115.6 (Cosentino et al., 2018) was differentiated into beta cells using a previously published 7-step protocol (Cosentino et al., 2018) . At the end of the stage 4, cells were seeded into 24-well Aggrewell 400 microwell plates (Stem Cell Technologies) at a density of 0.9 ·10 6 cells per well after which differentiation was carried out as described previously (Cosentino et al., 2018) .
Stage 7 differentiated beta cells were washed twice with PBS containing 0.5 mM EDTA and incubated in Accumax (SIGMA #A7089) for 8 min at 37 °C after which 50% volume of KnockOut Serum Replacement (ThermoFisher #10828028) was added to stop the reaction. After centrifugation at 400 g for 5 min at room temperature, cells were resuspended in 1 mL HAM's F-10 medium, supplemented as indicated before (Demine et al., 2020) . 75,000 cells in 500 μL medium were seeded per square of a Nunc Lab-Tek II ICC chamber (ThermoFisher).
Immunohistochemistry analyses were carried out largely as described previously (Demine et al., 2020) , Immunohistochemistry analyses were carried out largely as described previously (Ceulemans et al., 2017) , using primary antibodies against the following proteins: ASCL2 (Merck, MAB4418, clone 7E2, 1/5000), JUND (Atlas Antibodies, HPA063029, 1/50), HEYL (Atlas Antibodies, HPA076960, 1/200) and INS (Agilent, IR002, 1/100). Pictures were taken using a Leica DMLB (Leica Microsystems).
The integrated single-cell RNA-seq data and pySCENIC results can be explored interactively in SCope (Davie et al., 2018) . Loompy (v2.0.17) (Linnarsson Lab., 2015) was used to create the loom files which were uploaded to SCope. The embedding of the regulon and integrated gene expression based UMAP clustering, as seen in this article, were added to the loom file. This table is related to Figure 1 , 2, 3 and 4. This table is related to Figure 1 . This table is related to Figure 2 , 3 and 4. List of regulon specificity scores for all 142 regulons of non-diabetic alpha, beta, acinar and ductal cells.
This table is related to Figure 3 and 4. List of putative target genes for each regulon. This table is related to Figure 2 , 3 and 4.