Genetic risk converges on regulatory networks mediating early type 2 diabetes

Traditional immunohistochemistry.

Multiple sections from pancreatic head, body, and tail regions of 20 T2D and 11 age-matched ND donors were lightly PFA-fixed and prepared for immunohistochemistry and stained as described previously^41,49,51. Primary and secondary antibodies and their dilutions are listed in Supplementary Table 2. Amyloid, a marker of islet pathology in T2D^52,53, was visualized using a 2-min incubation in Thioflavin S (0.5% w/v; T-1892, Sigma) followed by a brief wash in 70% ethanol as described previously^16,47,54. Images were acquired at 20× with 2× digital zoom using a FV3000 confocal laser scanning microscope (Olympus) or a ScanScope FL (Leica/Aperio) and processed using cytonuclear algorithms (HighPlex FL v3.2.1) or tissue classifiers via HALO software (Indica Labs) or morphometric measurement via Metamorph software v7.10 (Molecular Devices). Analyses were run on the entire tissue section or manually annotated islets as indicated in figure legends. Endocrine cell mass was quantified by using pancreas weight and the ratio of hormone positive cells as identified by cytonuclear logarithm within the entire pancreatic section from multiple blocks representing the head, body, and tail regions. To obtain islet capillary measurements, caveolin-1 channel was isolated and colour thresholding was used on a per-image basis to gather object data using the Integrated Morphometry Analysis function (Metamorph). The following analysis metrics represent mean ± standard error: endocrine cells (Fig. 2a and Extended Data Fig. 3e–g) 16,151 ± 1,715 islet cells per donor and 570,508 ± 51,866 total cells per donor; endocrine cell area (Extended Data Fig. 3c,d) 2.34 ± 0.24 mm² per donor; capillary morphology (Fig. 2e) 48 ± 4 islets per donor; macrophage area (Fig. 2i) 0.64 ± 0.07 mm² per donor; amyloid (Extended Data Fig. 3l) 108 ± 19 islets per donor; cilia (Extended Data Fig. 7d) 0.32 ± 0.05 mm² per donor; RFX6 (Fig. 4e) 1,863 ± 362 cells per donor; pseudoislets (Extended Data Fig. 8d) 2,797 ± 508 cells per sample. Islet grafts were stained, imaged, and analysed using the same technique described for tissue above, with graft area annotated by hand to exclude quantification of mouse kidney tissue. The following analysis metrics represent mean ± standard error: endocrine cells (Extended Data Fig. 1r) 2,109 ± 347 cells per donor; endocrine cell area (Extended Data Fig. 1r) 0.1606 ± 0.026 mm² per donor; vascularization (Extended Data Fig. 1s) 0.1728 ± 0.031 mm² per donor.

CODEX multiplexed imaging.

To visualize additional cell types modulating the islet microenvironment⁵⁵, a subset of samples was analysed using CODEX⁵⁶. Antibodies were purchased pre-conjugated from Akoya Biosciences or sourced from other vendors and conjugated in-house using the CODEX Conjugation Kit (Akoya Biosciences) or by Leinco Technologies (Supplementary Table 3). Ten-micrometre lightly fixed⁴¹ pancreas sections were mounted onto 22 × 22 mm glass coverslips (Electron Microscopy Sciences) coated in 0.1% poly-l-lysine (Sigma) and stained with the CODEX Staining Kit (Akoya Biosciences) in uncoated 6-well tissue culture plates (VWR) per manufacturer instructions. Fluorescent oligonucleotide-conjugated reporters were combined with Nuclear Stain and CODEX Assay Reagent (Akoya Biosciences) in light-protected 96-well plates sealed with foil (Akoya Biosciences) and automated image acquisition and fluidics exchange were performed using the Akoya CODEX instrument and CODEX Instrument Manager (CIM) v1.29 driver software (Akoya Biosciences) integrated with a BZ-X810 epifluorescence microscope (Keyence). Tissue was hydrated in 1× CODEX buffer (10× CODEX Buffer diluted in Milli-Q water) and hybridization/stripping of the fluorescent oligonucleotides was performed using dimethyl sulfoxide (Sigma). After loading of coverslip into stage insert, tissue was visualized with Nuclear Stain diluted 1:1,000 in PBS and imaging area was set by centre point and tile number using BZ-X810 viewing software (Keyence). All images were acquired using a CFI plan Apo I 20×/0.75 objective (Nikon) with 30% tile overlap and 5 z-planes (1.5 μm/z).

Processing and annotation of CODEX images.

A total of 16 tissue regions were captured from 6 ND and 10 T2D donors (mean 50 mm² tissue per donor). Image alignment, stitching, background subtraction, and deconvolution were performed using the CODEX Processor v1.7.0.6 (Akoya Biosciences; see https://help.codex.bio/codex/processor/technical-notes for details). Individual channel images (TIFF files) were imported into HALO software v3.1 (Indica Labs) for all analyses as described below. Tissue and islet areas were annotated by hand to exclude out-of-focus regions and poor tissue quality. Islets (estimated diameter ≥50 μm; mean 42 islets per donor) were annotated based on DAPI and CHGA channels. Cell segmentation and cell-type annotations were performed using the HALO HighPlex FL v3.2.1 module with consistent cytonuclear parameters (nuclear contrast threshold 0.456, maximum cytoplasm radius 0.48). Due to marker intensity variability among samples, thresholds were manually set for each marker and donor. Unless otherwise noted, cells were counted positive for a given marker if minimum intensity was reached in 50% of cytoplasm area (see Supplementary Tables 3 and 4 for a complete list of markers and cell types). For cells with more variable morphology, positivity was also counted for nuclear area (30%: ARG1, CD11c, CD14, CD163, CD206, CD31, CD34, CD45, HLA-DR, IBA1, KRT, MCAM). Proliferating cells were counted only if minimum 60% of nuclear area met KI67 intensity threshold. Vascular structures (CD31) were also measured by random forest classification algorithm (HALO Tissue Classifier module).

The following analysis metrics represent mean ± standard error: endocrine cell area (Extended Data Fig. 3h) 0.88 ± 0.10 mm² per donor; islet cell composition (Fig. 2c and Extended Data Fig. 3j) 7,322 ± 852 cells per donor; immune cells (Fig. 2h,j) 309 ± 43 cells per donor; endothelial cell phenotypes (Extended Data Fig. 4e) 460 ± 92 cells per donor; macrophage phenotypes (Extended Data Fig. 4h) 191 ± 29 cells per donor; T cell phenotypes (Extended Data Fig. 4i,j) 40 ± 17 cells per donor.

High-dimensional, spatial and neighbourhood analyses.

The R implementation of the UMAP algorithm (https://CRAN.R-project.org/package=umap) was used for dimensionality reduction. Cell marker percentages obtained through HALO were standardized across islets (n = 255 ND islets and 426 T2D islets; mean 172 cells per islet), and default parameters were used for UMAP reduction (Fig. 2k and Extended Data Fig. 4k) except for nearest neighbours (80) and minimum distance (0.05). For spatial analyses, CD31 area classifications were converted to an annotation layer. A nearest neighbour algorithm (HALO Spatial Analysis module) was applied to obtain average distance of endocrine cells (n = 4,830 ± 692 cells per donor) to islet capillaries (CD31⁺ region) (Fig. 2f and Extended Data Fig. 4d).

For CN analysis, two methods were applied in parallel to CODEX data from annotated islets. In the community detection method, termed Dynamic CF-IDF (Fig. 2l and Extended Data Fig. 5a,d), a weighted undirected heterogeneous graph for each islet was constructed based on the cell types and normalized distance between cells. A greedy-based graph community detection method⁵⁷ was applied to segment the graph into a set of cell communities, then cell communities were stratified into 6 cellular neighbourhoods (n = 5,582 total cellular neighbourhoods with median 11 cells per cellular neighbourhood). Cell-type enrichment was determined by a new proposed scoring function CF-IDF, which is a modification of the widely used text sequence analysis method term frequency–inverse document frequency scoring⁵⁸. Our cell frequency (CF)-inverse dataset frequency (IDF) score emphasizes the cell type that is not only prevailing, but also uniquely representative in a group of target islets. Therefore, it will deemphasize the most dominant cell types (for example, α and β) throughout all the islets while paying more attention to the relative enrichment of less abundant cell types (for example, vascular and immune cells) in the local regions. The downstream analysis not only introduces insightful results on T2D feature analysis but also shows a robust performance across different resolution levels.

The second cellular neighbourhood analysis method, a k-means approach (Extended Data Fig. 5b,e,g), built on a previously published algorithm used to identify cellular neighbourhoods in the tumour microenvironment⁵⁹. For each cell, we first found its 10 nearest neighbours in the islet and assigned the ith nearest neighbour that was an α cell, β cell, macrophage, endothelial cell or γ cell, a score cos(iπ/20). Then we calculated the total score for each cell type, applied L1 normalization to the scores, and standardized them across all cells. The resulting representations of cells were finally used for k-means clustering to form 5 cellular neighbourhoods (n = 5,021 total cellular neighbourhoods with median 5 cells per cellular neighbourhood).

RNA isolation and bulk RNA sequencing.

RNA was extracted from sorted α and β cells (see ‘Purification of α and β cells by FACS’) or from pelleted whole islets using the Invitrogen RNAqueous-Micro Total RNA Isolation kit (Thermo Fisher AM1931). TURBO DNA-free (Ambion) was used to treat any trace DNA contamination. RNA was quantified by Qubit Fluorometer 2.0 and RNA integrity (RIN) was confirmed (RIN > 7) by 2100 Bioanalyzer (Agilent). Amplified cDNA libraries were constructed using SMART-seq v4 Ultra Low Input RNA-kit (Takara) and sequencing was performed on an NovaSeq platform (Illumina) using paired-end reads (100 bp) and 25 million reads per sample.

We processed the raw RNA-seq reads using FastQC (v0.11.8) for broad quality assessment. In brief, we examined the following parameters: (1) base quality score distribution, (2) sequence quality score distribution, (3) average base content per read, (4) GC distribution in thereads, (5) PCR amplification issue, (6) overrepresented sequences, (7) adapter content. Based on the quality report of fastq files, we trimmed sequence reads using fastq-mcf (v1.05) and cutadapt (v2.5) to only retain high-quality sequence for further analysis. The paired-end reads were aligned to the GRCh37/hg19 human reference with GENCODE v19 gene annotation using STAR splice-aware aligner (v2.5.4b;–outSAMUnmapped Within KeepPairs)⁶⁰.

We counted fragments mapping to features type in GENCODE v19 gene annotation using featureCounts from Subread package⁶¹. The gene list was pruned to contain only protein-coding genes mapping to autosome and chrX, resulting in a total of 20,260 genes. We assessed libraries using comprehensive quality metrics generated by QoRTs⁶² as well as computed derived metrics. In brief, on the top of QoRTs reported metrics, we computed: (1) 5′–3′ gene coverage bias (as the ratio of coverage values at the 90th percentile and 10th percentile of the coverage distribution); (2) Kolmogorov–Smirov test statistic between cumulative gene diversity of each library relative to median distribution of all libraries within each cell type and standardized to a mean of 0 and s.d. of 1 to yield a z-score; (3) number of reads mapped mapped to Xist and SRY genes; (4) average number of reads mapped to chrM; and (5) transcript integrity number (TIN)⁶³ for each library. The labelled sex of donors was matched against the gene expression quantified for sex genes to rule out any sample swaps or mislabelling. We also computed principal components for TPM (transcript per million) normalized count matrix for each cell type in order to detect potential outliers.

Differential gene expression analysis.

As collection of these rare tissues spanned more than five years, we used a latent variable analysis to discern biological variation from technical variation and then examined the datasets by both differential gene expression and gene network analyses. We performed differential gene expression analysis between T2D and ND samples for each cell type individually using DESeq2⁶⁴. In order to minimize potential effects of known and unknown confounding factors, we included known covariates in the DESeq2 model as well accounted for unknown covariates using the RUVseq latent variable approach⁶⁵. In brief, we used the following multi-step process: (1) we first removed genes from the raw count matrix that had less than 10 reads in fewer than 25% of the samples for that cell type. (2) We then ran a first-pass differential expression analysis using DESeq2 with age, sex, BMI and batch as known covariates. The output result was filtered for genes that were non-significant (that is, not differentially expressed between T2D and ND samples) and had P values > 0.5. These genes were used as ‘control’ or ‘empirical’ genes for the RUVSeq::RUVg function to estimate latent variables accounting for variation in the data not attributed to disease status. (3) The latent variables estimated from the RUVseq run were then used as additional covariates (on the top of age, sex, BMI and batch where applicable) for the second run of DESeq2. We selected the number of latent variables to provide the most reasonable separation between T2D and ND samples and minimal deviation from mean in the relative log expression plots. The output results from DESeq2 were filtered for 1% FDR to generate the final list of genes that were differentially expressed between T2D and ND for each cell type. We performed functional enrichment analysis using RNA-Enrich⁶⁶ and retained terms with an FDR threshold of 5%. Terms were condensed using the RelSim function in REVIGO⁶⁷ with similarity parameter set to 0.5 and visualized in semantic space using an.xgmml file imported into Cytoscape software⁶⁸ v3.8.2.

Combined analysis of differentially expressed genes (fold change ≥1.5 or ≤−1.5; P < 0.01) was performed using Metascape⁶⁹. For visualization of overlap between samples (Fig. 1i), each gene list (α, β and islet) was independently analysed to calculate all statistically significant (P < 0.01) ontology terms enriched (two-tailed hypergeometric test). Grey curves on the diagram that connect two lists indicate that an ontology term containing ≤100 genes was statistically enriched in both lists (that is, the genes in the two lists are not identical, but they are both members of terms that are statistically enriched in their respective list). For visualization of enrichment across samples (Fig. 1j), all three gene lists (α, β and islet) were provided and Metascape’s heuristic algorithm sampled the 20 top-score clusters, selected up to the 10 best scoring terms (lowest P values) within each cluster, and connected terms pairs with Kappa similarity above 0.3. The resulting network was exported as a.cys file and visualized using Cytoscape, with the most representative term name in each cluster selected manually.

Gene network analysis.

We adopted the WGCNA⁷⁰ approach to create networks from the gene expression data. This approach created modules (eigengenes) of up to 2,000 genes each, labelled by sample type and numbered consecutively (β cells, modules β00–β48; α cells, modules α00–α54; and islets, modules i00-i67). Importantly, collapsing the expression patterns across >14,000 genes into a smaller number of modules reduced gene-level multiple testing burden and enabled association of transcriptomic profiles with sample features including donor traits⁷¹, islet functional parameters from the same donors defined by dynamic islet perifusion, and enrichment of open chromatin peaks to overlap GWAS variants. In brief, we first filtered genes following the same rule established in ‘Differential gene expression’, where we only kept genes that had at least 10 reads in at least 25% of the samples for each cell type. We then processed raw counts using the varianceStabilizedTransformation function in DESeq2 package and used removeBatchEffect from the limma R package⁷² to adjust for effects of age, sex and BMI while protecting for disease status in the design matrix. The normalized and batch-corrected count matrix was then used as input to blockwiseModules to create a ‘signed hybrid’ network with ‘bicor’ as the correlation function. The power (k) parameter was selected such that the scale free topology fit reached at least 80% fit. To examine cell-type modules associated with quantitative traits of interest, we utilized a linear regression-based framework. We: (1) inverse-normalized the raw quantitative trait; (2) adjusted for age, sex and BMI by linear regression; and (3) computed the spearman rank correlation between residuals and eigengene of all modules. Within each network, we also computed the module membership score and network connectivity for each gene. Estimated enrichment of curated gene lists^22,73,74 (Supplementary Table 5) was calculated using Fisher’s exact test. Functional enrichment of genes in each module was performed using gprofiler2⁷⁵, and the results were visualized as a dotplot.

Integration of network analysis with chromatin accessibility.

We integrated chromatin accessibility information with gene network analysis using single-cell combinatorial indexing on ATAC-seq (sci-ATAC-seq) data for α and β cells derived from our previously published study⁴. For each module within each cell type, we selected: (1) accessible sites that were present within a specified distance of the transcription start site (TSS) of the genes within that module; and (2) the distal chromatin peaks that were linked to the peaks within this set based on the Cicero peak interaction results from the same study. This set of TSS proximal and distal peaks for all of the genes within each module and for each cell type were then used for downstream enrichment analyses.

For variant enrichment analysis in the module linked peaks, we collected the latest published summary statistics for selected traits^3,76. Using a threshold of ±10 kb to define our gene TSS boundary for linking peaks with modules, we created a set of accessible sites for each module. The union of peaks across all modules was used as a ‘bulk’ positive enrichment control. We then tested the enrichment of trait-associated variants from multiple GWAS across module peaks using GARFIELD³⁷ and used a P value threshold of 5 × 10⁻⁸ as input parameter for selecting trait-associated variants.

Next, we considered whether specific transcription factor binding motifs (TFBMs) are enriched to occur in certain modules. To test this, we defined module linked peaks for each module as described before but using a threshold of ±1 kb from gene TSS. For each peak within a module, we then identified the peak summit and extended the summit by 50 bp in each direction. Using genomic sequence in this region as our test sequence, we used the Analysis of Motif Enrichment (AME, v5.3.2) tool from MEME suite⁷⁷ (using default parameters) to identify enriched TFBMs represented in cisBP v.2.0⁷⁸. The control set of sequence was generated using –shuffle–parameter in AME that generates a control sequence by shuffling the test sequence but preserving the 2-mer frequency. The enrichment score was computed as scaled log₂-transformed (true positives + 1)/(false positives + 1) for each TFBM.

Isolation of nuclei.

Pseudoislet samples treated with RFX6 shRNA or scramble RNA were pooled together using a randomized study design, so the targeting and scramble conditions were not confounded by batch (Fig. 5a). To accomplish this, samples were allocated into 6 groups (batches) of n = 490–494 pseudoislets for nuclei isolation. A customized protocol was developed based on recommendations by 10x Genomics (https://www.10xgenomics.com/resources/demonstrated-protocols/), which included optimization steps described below. In brief, the samples were suspended in 1× PBS and pelleted at 2,000g for 3 min at 4 °C. The pellet was resuspended in lysis buffer (10 mM Tris-HCl 7.4 pH, 10 mM NaCl, 3 mM MgCl₂, 0.1% Tween-20, 0.1% NP40, 0.01% digitonin, 1% BSA, 1 mM DTT, and 2 U μl⁻¹ RNAse inhibitor) and rocked in an Eppendorf thermomixer C (EP 5382000015) at 300g for 5 min at 4 °C. Keeping the samples on ice as much as possible, tubes were then transferred to a pre-chilled 2-ml glass dounce homogenizer and homogenized with 15 strokes of tight pestle B before being transferred to a 1.5-ml tube and centrifuged at 500g for 5 min at 4 °C. The resulting pellet was then resuspended in 1 ml wash buffer (10 mM Tris-HCL 7.4 pH, 10 mM NaCl, 3 mM MgCl₂, 1% BSA, 0.1% Tween-20, 1 mM DTT and 2 U μl⁻¹ RNAse inhibitor) and centrifuged at 100g for 1 min at 4 °C. The supernatant was collected, filtered through a pre-wetted 30-μm filter, and centrifuged at 500g for 5 min at 4 °C. Nuclei were resuspended in 300 μl of wash buffer, then 300 μl of sucrose cushion (0.88 M sucrose, 1 mM DTT, 1 mM RNAse inhibitor, and 10% wash buffer) was added to the bottom of the tube and the resulting layered solution was centrifuged at 1,000g for 10 min at 4 °C. Both layers of supernatant were removed, and pellet was resuspended in 1 ml wash buffer and centrifuged at 500g for 5 min at 4 °C. Nuclei were then resuspended in 30 μl of nuclei resuspension buffer before counting and quality assessment. The desired concentration of nuclei was achieved by resuspending the appropriate number of nuclei in 1× diluted nuclei buffer for joint (on the same nucleus) snATAC-seq and snRNA-seq multiome profiling. Nuclei were processed by the University of Michigan Advanced Genomics Core using the 10x Genomics Chromium platform at 20,000 nuclei per well.

Multiome sample genotyping and imputation.

Samples were genotyped with the Infinium Multi-Ethnic Global-8 v1.0 kit using 50 ng I⁻¹ DNA samples in two batches. Probes were mapped to Build 37. We merged the.ped files for the two batches along with samples from other projects that were genotyped on the same chip (resulting in a combined 68 samples). We removed variants with multi mapping probes and updated the variant rsIDs using Illumina support files Multi-EthnicGlobal_D1_MappingComment.txt and Multi-EthnicGlobal_D1.annotated.txt (downloaded from https://support.illumina.com/downloads/infinium-multi-ethnic-global-8-v1-support-files.html). We performed pre-imputation QC using the HRC-1000G-check-bim. pl script (version 4.2.9) obtained from https://www.well.ox.ac.uk/~wrayner/tools/ to check for strand, alleles, position, Ref/Alt assignments and update the same based on the 1000 Genomes reference (https://www.well.ox.ac.uk/~wrayner/tools/1000GP_Phase3_combined.legend.gz). We did not conduct allele frequency checks at this step (that is, used the–noexclude flag) since we had 68 samples from mixed ancestries. These filters resulted in 958,427 variants. We performed pre-phasing and imputation using the Michigan Imputation Server⁷⁹. The standard pipeline (https://imputationserver.readthedocs.io/en/latest/pipeline/) included pre-phasing using Eagle2⁸⁰ and genotype dosage imputation using Minimac4 (https://github.com/statgen/Minimac4) and the 1000 Genomes phase 3 v5 (build GRCh37/hg19) reference panel⁸¹. Post-imputation, we selected biallelic variants with estimated imputation accuracy (r²) > 0.3, variants not significantly deviating from Hardy Weinberg equilibrium (P > 10⁻⁶), minor allele frequency (MAF) in 1000 Genomes European individuals > 0.05 and minor allele count (MAC) > 1 in our 12 samples, resulting in 6,665,607 variants.

Data processing for RNA component.

The RNA component of the multiome data was processed using starSOLO (STAR v. 2.7.3a, with GENCODE v19 annotation; options–soloUMIfiltering MultiGeneUMI–soloCBmatchWLtype 1MM_multi_pseudocounts–soloCellFilter None), which outputs the count matrices needed for most of the analyses⁶⁰. Quality control metrics were gathered on a per-nucleus basis using a custom Python script on the corrected gene counts and aligned BAM file.

Following processing with STAR, we constructed a custom count matrix by combining information from the GeneFull and Gene matrices output by STAR. The GeneFull matrix contains per-gene counts based on intronic and exonic reads, while the Gene matrix contains per-gene counts based on exonic reads only. As nuclear RNA may contain introns, the GeneFull matrix should be preferred. However, due to overlapping transcript annotations that render some read gene assignments ambiguous, some genes may receive fewer counts in the GeneFull matrix than in the Gene matrix. The INS gene was an extreme example of this, receiving very low counts in the GeneFull matrix but high counts in the Gene matrix. To salvage counts for such genes, our custom matrix utilized the GeneFull counts for most genes but utilized the Gene counts for the subset of genes that had greater counts in the Gene matrix than in the GeneFull matrix.

Data processing for the ATAC component.

Adapters were trimmed using cta (https://github.com/ParkerLab/cta). We used a custom Python script, available in the Parker laboratory Github repository, for barcode correction. Barcodes were corrected in a similar manner as in the 10x Genomics Cell Ranger ATAC v. 1.0 software. In brief, barcodes were checked against the 10x Genomics whitelist. If a barcode was not on the whitelist, then we found all whitelisted barcodes within a hamming distance of two from the bad barcode. For each of these whitelisted barcodes, we calculated the probability that the bad barcode should be assigned to the whitelisted barcode using the Phred scores of the mismatched base(s) and the prior probability of a read coming from the whitelisted barcode (based on the whitelisted barcode’s abundance in the rest of the data). If there was at least a 97.5% probability that the bad barcode was derived from one specific whitelisted barcode, it was corrected to the whitelisted barcode.

Reads were mapped using BWA-MEM⁸² with flags ‘-I 200,200, 5000 -M’ (v. 0.7.15-r1140). We used Picard MarkDuplicates (v. 2.25.1; https://broadinstitute.github.io/picard/) to mark duplicates, and filtered to high-quality, non-duplicate autosomal read pairs using SAMtools view⁸³ with flags ‘-f 3 -F 4 -F 8 -F 256 -F 1024 -F 2048 -q 30’ (v. 1.10). Quality control metrics were gathered on a per-nucleus basis using ataqv⁸⁴ (v. 1.2.1) on the BAM file with duplicates marked.

Selection of quality nuclei for downstream analysis.

We performed rigorous QC of all RNA nuclei and only included those deemed as high-quality based on the following four definitions: (1) number of unique molecular identifiers (UMIs) >1,000; (2) mitochondrial fraction <0.2; (3) nuclei where the RNA profile was statistically different from the background or ambient RNA signal; and (4) nuclei were identifiable as a singlet and assignable to a sample using genotypes. We considered droplets with UMIs <10 to be empty and therefore representative of the background or ambient RNA profile. Top genes in the ambient RNA included highly expressed genes across prominent islet cell types such as INS, GCG and SST, along with several mitochondrial genes. We used the testEmptyDrops function from DropletUtils (v 1.6.1)⁸⁵, specifying the lower parameter as 10 and selecting droplets with P < 0.05 as droplets significantly different from the ambient RNA profile. To identify singlets and assign to samples, we ran Demuxlet⁸⁶ using using the BAM files and the genotype VCF file considering all post-QC variants in gene bodies with minor allele count (MAC) > 1. We used the command “demuxlet–sam $bam–tag-group CB–tag-UMI UB–vcf ${vcf}–alpha 0–alpha 0.5–field GT”, and selected singlets. To account for ambient RNA contamination while identifying singlets, we also masked the top 1% genes expressed in the ambient RNA and re-ran Demuxlet with the same parameters; nuclei were considered singlets and kept for downstream analysis if they were called as singlets in either Demuxlet run.

We also performed QC of the ATAC component of the multiome data. For ATAC, we required nuclei to have a minimum TSS enrichment (as calculated by ataqv) of 2, minimum filtered read count of 1.000 (ataqv ‘HQAA’ metric), and maximum mitochondrial fraction of 0.5. We also ran Demuxlet on the ATAC component (command: “demuxlet–sam $bam–tag-group CB–vcf ${vcf}–field GT”) and required that a prospective nucleus be called as a singlet. The ATAC component of nuclei in two wells showed low TSS enrichment and all nuclei from these two wells were therefore excluded from analysis.

If the RNA and the ATAC component of a barcode both passed QC and the Demuxlet sample assignment was the same, both modalities were utilized for downstream analysis. If only the RNA component passed QC, only the RNA component was used in downstream analysis. As we performed clustering on the RNA component, we excluded the few (twelve) barcodes that passed ATAC QC and failed RNA QC. The final count was 15,825 (RNA) and 5,706 (ATAC) high-quality nuclei for downstream analysis.

Removal of ambient RNA counts from single-nucleus gene expression UMI matrices.

Prior to clustering and downstream analysis, we used DecontX⁸⁷ (celda v. 1.8.1, in R v. 4.1.1)⁸⁸ to adjust the nucleus × gene expression count matrices for ambient RNA. DecontX was run on a per-batch basis, as the amount of ambient contamination may vary across batches. Decontaminated counts were generated via the decontX() function, passing barcodes with total UMI count ≤10 to the background argument. Rounded decontaminated counts were used for clustering and all downstream analyses. Nuclei with estimated contamination level >0.2 were excluded from downstream analysis.

Clustering of multiome data.

Nuclei were clustered on the RNA component using Seurat^89–91 (v. 3.9.9.9010, in R v. 3.6.3). After normalizing counts with the NormalizeData function, we identified the top 2,000 variable features (FindVariableFeatures function, with selection. method=‘vst’) and scaled with the ScaleData function. We identified neighbours using the top 20 principal components and k.param = 20, and called clusters using resolution = 0.1 with n.start = 100. We used the top 20 principal components for generating the UMAP.

This clustering protocol identified 10 clusters. One of the smaller clusters shows expression of both INS and GCG, suggesting it may consist of doublets that were not caught by demuxlet. To verify this was a doublet cluster, we ran a different, genotype-independent, ATAC-based doublet detection method (AMULET; v. 1.0-beta, run with default parameters separately on data from each multiome well)⁹² on the ATAC nuclei that otherwise passed QC. This method tagged ~40% of the nuclei in the suspected doublet cluster as doublets, while only ~5% of nuclei in any other cluster were tagged as doublets. We therefore removed the small doublet cluster from the clustering and downstream analysis. Data are available via the UCSC Cell Browser⁹³ at https://theparkerlab.med.umich.edu/data/public/cellbrowser/?ds=Pseudoislet10XMultiome for further exploration.

Differential gene expression analysis.

Differential gene expression was performed within each cluster using DESeq2 (v. 1.28.0)⁶⁴ on pseudobulk counts. UMI counts were summed across nuclei within a donor + construct + cluster. Only donors with paired data (RFX6–2896 and scrambled–mCherry constructs) were used, and the analysis was performed in a paired fashion (DESeq2 model: ~donor + construct). We used an FDR threshold of 5% for considering genes differentially expressed. To compare differential gene expression between multiome data and data from sorted β cells (Extended Data Fig. 9g), Metascape was utilized as described above (‘Transcriptional analysis of α and β cells and islets from ND and T2D donors’, ‘Differential gene expression analysis’).

Gene network exploration for the RNA component.

Single-cell regulatory network inference and clustering (SCENIC)⁹⁴ was applied to the RNA modality (~15,300 nuclei) to identify cell-specific regulons (small gene regulatory networks of transcription factors and their target genes) and discern changes induced by RFX6 knockdown. Soft gene filtering was applied to post-QC nuclei to remove genes present in <1% of all the nuclei and used as the input for SCENIC. In brief, the following steps were performed: network inference based on GRNBoost2, regulon prediction (cisTarget), and cellular enrichment (AUCell). Adjacency matrices, listing interactions of transcription factors and a numerical score, were generated, as were lists of regulons, their regulon specificity scores (RSS), and UMAP of nuclei based on RSS. Using previously identified cluster labels for each nucleus, average RSS was calculated for each regulon to generate a ranked list for each cluster. For comparison between shRFX6 and control nuclei in the β cluster, absolute difference of the RSS score for each regulon was computed between the two groups to identify the most differential regulons. Pathway enrichment was performed using Metascape as described above, using the input of genes (n = 220) contained in top 10 differential β regulons.

Per-cluster processing of ATAC component.

All ATAC reads from pass-QC, clustered nuclei were merged within each cluster. To generate per-cluster peaks, these BAM files were converted to single-ended BED format using bedtools bamtobed⁹⁵ before calling ATAC-seq peak summits with MACS2⁹⁶ (flags -g hs–nomodel–shift −37–extsize 73 -B–keep-dup all–call-summits). We removed summits in blacklist regions, filtered to FDR 0.1% summits, and then generated a peak list from the summits by extending the ATAC-seq peak summits for each cluster +/− 150 bps to get 300-bp peaks (within each cluster, if two 300-bp peaks overlapped the one with the greater MACS2 score was kept). We then removed peaks in blacklist regions. To get the ATAC peak counts used in the ATAC principal components analysis (PCA) and differential chromatin accessibility analyses, we determined the number of ATAC fragments overlapping each of these peaks in each of the per-cluster, per-donor, per-construct pseudobulk samples.

For visualization of ATAC signal, we generated a normalized bedGraph file using MACS2 on the single-end BED file (macs2 callpeak command, with options–SPMR–nomodel–shift −100–extsize 200 -B–broad–keep-dup all) and then converted to bigWig format using the UCSC bedGraphToBigWig⁹⁷. For PCA on the pseudobulk ATAC counts, we first removed any peaks on the mCherry or mKate2 contigs. We then converted peak counts to counts per million and removed the bottom 10% of features with the lowest average CPM across samples. For each peak, we filled any zeros with a value equal to half of the minimum non-zero CPM for that peak across samples. We then log transformed prior to performing the PCA.

Differential chromatin accessibility analysis.

Differential chromatin accessibility was performed within each cluster using DESeq2 (v. 1.28.0)⁶⁴ on pseudobulk ATAC peak counts. Only donors with paired ATAC data (RFX6–2896 and scrambled–mCherry constructs) were used, and we additionally excluded donor 17277513 due to very low read counts. The DESeq2 analysis was performed in a paired fashion, with model: ~donor + tss_enrichment + construct. To compute TSS enrichment for each pseudobulk sample, we merged all ATAC nuclei (regardless of cluster) from each donor and computed TSS enrichment with ataqv.

Testing for enrichment of peak subsets near differential genes.

We used a permutation test to determine whether the most significant peaks (‘top peaks’) from the beta cell differential peak analysis were enriched near beta cell differentially expressed genes. First, we assigned each peak to the gene with the nearest TSS (if multiple TSS were equally close, we took the TSS with the smallest chromosomal coordinate). We then calculated the fraction of top peaks whose nearest gene was differentially expressed. To get the null expectation for this value, we permuted the ‘DE/not DE’ gene labels, such that the same number of genes were always labelled as ‘DE’ but the identity of these differentially expressed genes changed in each permutation. While permuting, we split genes into deciles based on the expression of each gene and permuted the labels only within each decile (this controls for the fact that highly expressed genes are more likely to be differentially expressed than lowly expressed genes due to statistical power in the differentially expressed analysis). We performed 10,000 permutations, in each permutation re-calculating the fraction of top peaks whose nearest gene was differentially expressed to build up the null distribution. We then calculated an empirical P value based on our observed value and the null distribution, adding a pseudocount to avoid a P value of 0 (P = (1 + no. of permutations where the test statistic was greater than or equal to our observed value)/10,001).

Motif scanning for multiome motif enrichment analyses.

The motif scans were performed using FIMO (v. 5.0.4) with a background model calculated from the hg19 reference genome⁹⁸ and otherwise default parameters. We used the motifs from Kheradpour and Kellis 2014⁹⁹, excluding “*_disc” motifs; motifs from cisBP v. 2.0⁷⁸; motifs from Jolma et al. (2013)¹⁰⁰; and custom RFX6 motifs generated using mouse Rfx6 ChIP-seq data from Piccand et al. (2014)³³.

The custom RFX6 motifs were generated during a previous project²⁵. Sequencing reads from Piccand et al. (2014)³³ were mapped to the mouse mm9 genome¹⁰¹ using bwa (v. 0.7.12-r1039) and peaks were called using MACS2 (flags: -t MIN6_Rfx6-HA_IP.bam -c MIN6_Control-HA.bam -B–nomodel -g mm–keep-dup 1 -q 1.00e-4). The MEME (v. 4.11.0)¹⁰² and DREME (v. 4.9.1)¹⁰³ tools from the MEME suite¹⁰⁴ were used to discover novel motifs in the resulting peaks. One non-repetitive motif from the MEME tool and two motifs from the DREME tool, bearing similarity to known RFX family motifs, were selected for use in downstream analysis.

Motif enrichment in most significant peaks.

We used logistic regression to measure enrichment of motifs in subsets of ATAC-seq peaks. We ran one model per peak category and motif. For testing for enrichment in the peaks that had the smallest P values and leaned towards higher signal in shRFX6 samples, we modelled:

Where ‘peak_leans_higher_in_shRFX6’ is 1 if the peak was one of the most significant peaks in the ‘up in RFX6 KD condition’ direction and 0 otherwise; peak_gc_content was the GC content of the sequence within the peak; peak_size was the mean DESeq2-normalized count for the peak across the samples in the DESeq2 analysis; and n_motif_hits_in_peak was the number of motif hits in the peak as determined by the FIMO motif scans. The coefficient of the n_motif_hits_in_peak term was taken as the measure of motif enrichment. For testing for enrichment in the peaks that had the smallest P values and leaned towards lower signal in shRFX6 samples, we used the same model except the outcome variable was ‘peak_leans_lower_in_shRFX6’.

Generation of ATAC footprint plots.

To generate the ATAC footprint plots, we first separated the motif occurrences into those within the beta cells ATAC peaks and those outside of peaks. For each of these two groups, we computed an aggregate Tn5 cut matrix for the 500 bps on either side of the motifs, using beta cell ATAC reads from each individual donor plus construct (using the make_cut_matrix script within the atactk package (https://github.com/ParkerLab/atactk); options -a -r 500). The cut matrices were generated separately for each donor+construct, utilizing only donors with paired ATAC data (RFX6–2896 and scrambled–mCherry constructs) and additionally excluding donor 17277513 due to very low ATAC read counts. To reduce the impact of Tn5 insertion sequence bias, we normalized the Tn5 cut frequency at each position for the motifs in peaks by the corresponding frequencies for the motifs outside of peaks. To adjust for technical differences (for example, TSS enrichment) between the donors plus constructs, we then divided these normalized cut frequencies by the average normalized cut frequency between the −500 and −400 bp positions.

GWAS enrichment in most significant peaks.

We considered if β cell ATAC-seq peaks that score highly for differential accessibility, as measured by P value, are specifically enriched to overlap T2D GWAS variants. We compared the enrichment of T2D (adjusted BMI) GWAS variants (n = 3,062,361) to overlap top 5,000 ATAC-seq differential peaks leaning up and down with the remaining peaks for β cell using GARFIELD³⁷. Using a P value threshold of 10⁻⁵, we also performed a conditional analysis where GARFIELD evaluates if both annotations are conditionally independent of each other in the enrichment model. The coefficients corresponding to each annotation from the conditional enrichment model were shown along with the 95% confidence interval. To ensure robustness of our results, we repeated the analysis for top 2,000 (up and down each) and top 10,000 (up and down each) differential peaks.