The dataset for islet atlas was obtained from Camunas-Soler et. al. (islet1) (Camunas-Soler et al., 2020), Segerstolpe et. al. (islet2) (Segerstolpe et al., 2016), Enge et. al. (islet3) (Enge et al., 2017), Xin et. al. (islet4) (Xin et al., 2016), Dorajoo et. al. (islet5) (Dorajoo et al., 2017), Lawlor et. al. (islet6) (Lawlor et al., 2017), Wang et. al. (islet7) (Wang et al., 2016), and Marquez-Curtis et. al. (islet8) (Marquez-Curtis et al., 2022). For islet1, we excluded 3 T1D patients due to the small sample size. We excluded the R269 sample since the raw data was not available. For islet3, we only used adult samples. For islet5, we excluded “subject 3” sample since the raw data was not available. For islet7, we excluded HP15041 due to genotyping problems. We used cell-type information from the original papers. For islet6, we used marker genes threshold for cell typing based on Supplementary Table S3 in the original paper. We only considered major islet cell types: alpha, beta, gamma, delta, acinar, and ductal cells.

We used the Trinity Cancer Transcriptome Analysis Toolkit (CTAT) which is based on the GATK Best Practices pipeline (Liu et al., 2019; Brouard and Bissonnette, 2022) for the variant calling of single cells in scRNA-seq data (Fangal, 2020). We ran the pipeline with either “single-end” or “paired-end” for each corresponding dataset. We set the parameter “boosting_method”, which speeds up during the variant calling, to be “none”. We excluded SNVs from the known RNA-editing sites from RADAR (Ramaswami and Li, 2014) and RediProtal (Picardi et al., 2017) databases, which were provided by the CTAT pipeline. During the process, 13/1,464 islet1 cells, 24/2,038 islet2 cells, 13/1,578 islet3 cells, 13/1,507 islet4 cells, 222/333 islet5 cells, 9/593 islet6 cells, 59/268 islet7 cells, and 1/147 islet8 cells were excluded due to lack of confident SNVs or absence of raw data from a certain sample. After obtaining all the SNVs from single cells, we merged per donor using bcftools (1.13) (Danecek et al., 2021). Since somatic mutations could exist, we ran a germline mutation calling pipeline from Strelka2 (2.9.2) (Kim et al., 2018) with merged bam files (using Samtools; v1.11) (Danecek et al., 2021) for each donor to identify germline mutations. We only keep the “passed” germline mutation from Strelka2 overlapping the merged mutation profiles at the donor level from the CTAT pipeline. We only keep the genotype information if there is a read at that region. We excluded X, Y, and M chromosomes. We excluded variants with minor allele frequency (MAF) < 0.1 (genomicMateSelectR v0.2.0) (Wolfe, 2022) from both healthy individuals and T2D patients. Locational information for SNP was obtained from CTAT result.

We obtained the transcript and gene expression matrix by running RSEM (Li and Dewey, 2011; Westoby et al., 2018) (1.3.0) using rsem-calculate-expression--estimate-rspd--no-bam-output -p 16 with GRCh37 version 19 (Ensembl 74) (transcript file: *.isoforms.results, gene file: *.genes.results). As a result, we could align reads to 38,381 transcripts and 20,392 genes. The mean number of expressed genes per cell was 6,196 and the mean number of expressed transcripts per cell was 6,835.

We extracted genomic region of each gene by GenomicFeatures package (v 3.17) with [type = = “gene”, gene_status ! = “PUTATIVE”, no duplicated gene name] (Lawrence et al., 2013).

We used log2 (TPM + 1) matrix of transcript from RSEM for running Seurat (v 4.1.0) (Butler et al., 2018)pipeline (FindVariableFeatures, ScaleData, RunPCA, RunUMAP [dims = 1:30]). We corrected batch effect of different dataset by running Harmony (v 0.1.0) (Korsunsky et al., 2019). We performed clustering by FindNeighbor (reduction = “harmony”, dims = 1:30) and FindClusters (resolution = 1.2). We excluded misaligned cell type from each cluster to remove batch outliers. We generated pseudobulk cell type-specific transcript matrix by averaging the transcript expression level per sample. After that, we used inverse log transformation (2^x −1) for sQTLseeker2 input since they suggested using a non-log transformed matrix to better calculate Hellinger distance for transcript expression.

We used the sQTLseekeR2 (v1.1.0) (Garrido-Martin et al., 2021) package for cell type level sQTL. We filtered out the lowly expressed transcript and gene by “min.transcript.exp = 0.1” and “min.gene.exp = 0.1” using prepare.trans.exp function. We used sqtl.seeker function to obtain nominal p-value with “svQTL = TRUE” (to remove a potential false positive sQTL by inhomogeneity in the variance from different genotype groups), “genic.window = 50,000 (bp)”, “min.nb.ind.geno = 5 (at least five samples for a certain genotype)”, “ld.filter = 0.8 (Linkage disequilibrium threshold)”. Age and gender were used for covariables to regress the transcript expression level. For multiple hypotheses correction, we used sqtls function with “FDR = 0.1” and “FDR.svQTL = 0.01” with default parameters (md.min: the minimum MD (Maximum Difference) in relative expression: 0.05).

For DEG analysis, we normalized the gene expression from the Seurat package (v 4.1.0) (Butler et al., 2018) and performed the FindMarkers function with logfc.threshold = 0.25 and min.pct = 0.1. We further used p.adjusted < 0.05.

We used FastQTL (v2.184) (Ongen et al., 2016) to obtain cell type level eQTL with only the SNPs tested for sQTL analysis. The same covariables were used with “--window 1e7” parameter setting. The gene expression matrix was obtained by the same process as the transcript matrix. We only kept eQTLs with q-value < 0.1.

We obtained RNA binding protein (RBP) site from POSTAR3 (http://111.198.139.65/index.html) (Zhao et al., 2022) with only the “experimental results”. We obtained RBP overlapping SNP if the SNP is located within the RBP binding region. We only considered an RBP if it is expressed in the cell type (log2 (exp+1) > 0.25 and fraction of cells expressing > 0.1).

To analyze the association between GWAS related to pancreas disease and SNPs from sQTL, we collected “type 2 diabetes mellitus” (EFO ID: MONDO_0005148), “type 1 diabetes mellitus” (EFO ID: MONDO_0005147), and “pancreatitis” (EFO ID: EFO_0000278) from GWAS catalog (https://www.ebi.ac.uk/gwas/) (Sollis et al., 2023). We obtained “glycemic traits” from Ji et al. (Chen et al., 2021). Since GWAS SNPs were in hg38 coordinates, we converted hg38 genome coordinate into hg19 by LiftOver (https://liftover.broadinstitute.org/) (Hinrichs et al., 2006). Unknown alleles were subsequently removed. We obtained LD > 0.6 SNPs from each GWAS by SNiPA (setting: GRCh37.1000 genomes, Phase 3 v5, European, Ensembl 87) (https://snipa.org/snipa3/) (Arnold et al., 2015).

We used the “splicing events classification” pipeline from sQTLseekeR (v2.1) (Monlong et al., 2014).

We used EnrichR (Xie et al., 2021) with KEGG (Kanehisa, 2019), NCI-Nature (Schaefer et al., 2009), and Gene ontology (Gene Ontology et al., 2023). We required FDR < 0.1 for a pathway to be considered significant.

We collected eight full-length scRNA-seq pancreatic islet datasets (Segerstolpe et al., 2016; Wang et al., 2016; Xin et al., 2016; Dorajoo et al., 2017; Enge et al., 2017; Lawlor et al., 2017; Camunas-Soler et al., 2020; Marquez-Curtis et al., 2022; Supplementary Table S1), comprising of 4,841 cells from 54 healthy individuals and 2,574 cells from 21 type 2 diabetes (T2D) patients divided amongst six cell types (acinar, alpha, beta, gamma, delta, and ductal cells) (Supplementary Table S1). We used the CTAT pipeline (Fangal, 2020), which is based on the STAR aligner and the GATK-best practice variant calling pipeline. Importantly, it has been shown to perform well in a benchmark for genotype inference from scRNA-seq (Liu et al., 2019). To infer the genotype we combined all the mutation profiles from each cell within each individual. Since these profiles may contain somatic mutations, we only keep the germline variants predicted by Strelka2 (Kim et al., 2018). We further filtered out the SNPs with minor allele frequency (MAF) < 0.1. As a result, we could obtain ∼24,000 SNPs per individual, similar to the ∼35,000 SNPs that would be expected from exome sequencing (Marian, 2014) (Supplementary Tables S1, S2) The lower number of SNPs is expected since we are restricted to the genes that are expressed in our sample. For transcript quantification we used RSEM (Li and Dewey, 2011) since it has been reported as the best performing tool for isoform level expression quantification with full-length scRNA-seq data (Westoby et al., 2018). Prior to clustering we ran Harmony (Korsunsky et al., 2019) to correct the batch effect from healthy donors or T2D patients, respectively. Then, we removed 416 cells from the healthy and 283 cells from the T2D group that had been aligned to different cell types from the original paper, assuming they were batch outliers (Supplementary Figure S2; Supplementary Table S3). For cell type level sQTL analysis, we created pseudobulks for each cell type following the same strategy used in single-cell eQTL analysis (Cuomo et al., 2021; Yazar et al., 2022). However, due to the small number of cells, we only analyzed alpha, beta, delta, and ductal cells (Figure 1A). We also created a pseudobulk of pancreatic islets called “collapsed” to mimic bulk pancreatic islet data (Supplementary Table S4).

After preprocessing, we used sQTLseeker2 (Garrido-Martin et al., 2021) to obtain cell type level sQTLs (Figure 1B). We followed the recommendations by the authors of sQTLseeker2 and filtered out the lowly expressed transcripts. We applied a linear model to adjust the expression levels, using age and gender as covariates. A 50 kbp window centered around each splice junction was used to identify candidate variants since previous studies have shown that most sQTLs occur within this distance, while larger windows reduce power due to multiple-hypothesis correction (Garrido-Martin et al., 2021; Atla et al., 2022; Brotman et al., 2022). The number of variants tested was further reduced by merging those that had a linkage disequilibrium > 0.8. Not all cell types were annotated in all datasets, and for each cell type we had at least 15 individuals per genotype.

Our analysis yielded 136 sQTLs from alpha cells, 92 from beta cells, while only 1 or 2 sQTLs were significant in other cell types at an FDR of 0.1. For the collapsed dataset we identified 268 sQTLs (Supplementary Tables S5, S6). Previous analyses using bulk data in the GTEX database (Shan et al., 2019) have shown that the number of significant eQTLs is strongly correlated with the sample size (Pearson correlation coefficient: 0.86) (Supplementary Figure S1A). We also observed a similar pattern from our data by down-sampling (Supplementary Figure S1B). We assumed a similar relationship exists for sQTLs, and extrapolating from a brain study, the expected number of sQTLs from 54 samples is ∼300 (63 samples: 387, 124 samples: 849) (Zhang et al., 2020), similar to the 268 sQTLs that we observed in the collapsed data. Similarly, a study of sQTLs using bulk RNAseq from pancreatic islets (Atla et al., 2022), reported 4,858 sQTLs from 400 samples using an FDR of 0.01. Assuming that the number of sQTLs scales linearly with the sample size, then 268 sQTLs for 54 samples is an expected result. Taken together, these comparisons show that the number of sQTLs from the “collapsed” data was within the expected range. For T2D patients, due to the sample size, we could not get any suitable number of significant results (Supplementary Table S5).

We also note that the effect sizes are often substantial. For both alpha and beta cells, > 92% of sQTLs result in a change in a maximum difference of relative abundance > 0.2 which is considered a large effect size in the sQTLseeker2 paper (Garrido-Martin et al., 2021). Comparing the distribution of sQTLs across cell types we found a substantial number of sQTLs and genes showing the alteration in the splicing affected by these variants (sgenes) that do not overlap with the “collapsed” result (Figure 2A). This result is consistent with single cell eQTL studies which have shown that bulk transcriptomes can mask the signal from individual cell types. For both alpha and beta cells, the overlap was only around 15%–20%, indicating that AS in different cell types is regulated by different SNPs (Figure 2B). At the same time, the splicing of a given gene was also affected differently between different cell types with only 15.8% of sgenes found in both alpha and beta cells, suggesting that the set of sQTLs affecting a given pathway will differ between cell types.

Interestingly, the expression levels of sgenes are typically not higher in the cell type where they are active. We performed differentially expressed gene (DEG) analysis between alpha and beta cells to obtain highly expressed genes for each cell type. The 32% for alpha cells and 20% for beta cells is the ratio of the DE genes among sgenes. However, we only found 1.7% of highly expressed genes in the alpha cell (27/1,626), and 0.32% of highly expressed genes in the beta cell (9/2,780) were cell type-specific sGenes. Hence, this implies that DE genes are distinct from sgenes.

Next, we characterized the sQTLs’ genomic locations (Figure 3A; Supplementary Table S7; Fangal, 2020), and this revealed that most sQTLs were located outside of exonic regions, with only 26/137 alpha cell variants and 15/92 beta cell variants overlapping exons. Comparing each region to the total number of detected variants suggested that SNPs from intergenic regions were not enriched above chance level since the odd ratios were close to 1. A previous study of sQTLs in the pancreatic islet showed that 5′ and 3′ UTRs were the most enriched regions (Atla et al., 2022). Although we could not find as many sQTLs from the 5′ UTR as in the previous report, we observed substantial numbers of sQTLs in 3′ UTR (odd ratios, alpha: 2.14, beta:1.77, collapsed: 1.57). We also observed a higher number of sQTLs in intronic regions which was not surprising given its key role in splicing (Jo and Choi, 2015), but the odd ratios were surprisingly low (alpha: 0.33, beta: 0.7, collapsed: 0.57). Although the splice site has been shown to be enriched for sQTLs (Garrido-Martin et al., 2021; Atla et al., 2022; Qi et al., 2022), we did not observe this in our study. One explanation could be that the splice junction can only be observed from unspliced transcripts. We also cannot rule out that the discrepancy from previous studies is due to variant calling from the transcripts. Second, we investigated whether sQTLs are associated with eQTLs (Supplementary Table S8). A previous bulk pancreatic islet study only reported a weak association between sQTLs and eQTLs (∼25%) (Atla et al., 2022), suggesting that regulation of gene expression levels and splicing are largely independent. Similarly, when we compared cell type level sQTLs and eQTLs, we found a low overlap (<27%) (Figure 3B), suggesting that gene expression and splicing regulation are largely independent also at the cell type level. Third, we hypothesize that sQTLs would be enriched at RNA binding protein (RBP) (Zhang et al., 2020; Garrido-Martin et al., 2021; Qi et al., 2022) sites. We analyzed whether sQTLs are enriched at RBP sites using the POSTAR3 database (Zhao et al., 2022), and as expected odd ratios were higher than 1 for alpha, beta, and “collapsed” (Figure 3C; Supplementary Table S9). Fourth, since many sQTLs have been associated with disease (Zhang et al., 2020; Garrido-Martin et al., 2021; Atla et al., 2022; Brotman et al., 2022; Qi et al., 2022), we investigated overlap with GWAS variants related to type 1 diabetes (T1D), type 2 diabetes (T2D), and glycemic traits (Atla et al., 2022) within LD > 0.6 SNPs using SNiPA (Arnold et al., 2015). However, we could not observe any GWAS association for our sQTL results, and we conjecture that this is due to the small number of variants identified.

Next, we characterized the splicing events involved in sQTL using the definitions from sQTLseeker (Monlong et al., 2014). Most of the splicing events were complex 5′ events, complex 3′ events, and skipped exon events, except unknown events (Monlong et al., 2014; Figure 4A; Supplementary Table S10). This finding was consistent for alpha, beta, and “collapsed”. Compared to previous reports, we observed a similar number of complex 5′ and 3′ events, but we observed more skipped exons (Monlong et al., 2014). We also carried out functional analysis for sgenes using EnrichR (Xie et al., 2021) to quantify the overlap with KEGG (Kanehisa, 2019), NIC-Nature (2016) (Schaefer et al., 2009), and Gene ontology (2023) (Gene Ontology et al., 2023; Figure 4B; Supplementary Table S11). A previous bulk pancreatic islet sQTL study (Atla et al., 2022) showed enrichment for terms related to virus infection (herpes simplex virus 1 infection and viral myocarditis) from the KEGG pathway. Our “collapsed” result was also enriched in bacterial response (yersinia infection and shigellosis), indicating an association with the inflammation pathway. Previously, it has been reported that the association between inflammation and pancreatic islets such as microbial DNA mediated pancreatic islet inflammation in obesity (Gao et al., 2022) or bacterial impact on T1D (Abdellatif et al., 2019). Hence, this result indicates that one of the mechanisms through which the genotype impacts obesity or T1D is via changed splicing pattern in response to inflammation. However, only a few pathways were enriched in the “collapsed”, indicating a dilution effect by averaging the signals from different cell types.

On the other hand, there were several enriched pathways from both alpha and beta cells that were found in only one cell type (Figure 4B; Supplementary Figure S3). One interesting finding was that the alpha cells are enriched for inflammation and fatty acid pathways, while the beta cells were enriched for the mTORC1 pathway. Both the association between the alpha cell and fatty acid pathway (Armour et al., 2023), and the association between the beta cell and mTOR signaling (Blandino-Rosano et al., 2012; Blandino-Rosano et al., 2017; Asahara et al., 2022) were reported previously from a bulk study (Atla et al., 2022), indicating cell type specific regulation. The alpha cells were also enriched in TGF-beta signaling, which was surprising since a tight association between beta cells and TGF-beta signaling has been reported previously (Lee et al., 2021; Wang et al., 2022). This result implies an additional functional association of TGF-beta signaling in alpha cells.

Another interesting pathway was the cell morphology pathway mediated by CDC42 which was enriched in both alpha and beta cells. It has been reported that CDC42 is related to cytoskeleton and cellular structure (Filipenko et al., 2005; Watson et al., 2017; Pichaud et al., 2019). This suggests that specific variants can alter pancreatic islet morphology via splicing alterations in a cell type specific manner. In addition, it has been reported that the morphology of the pancreatic islet is associated with diabetes (Bonner-Weir et al., 2015; Okajima et al., 2022). Especially, CDC42 was reported as a regulator for not only diabetes, but also insulin secretion and development (Martinelli et al., 2018; Veluthakal et al., 2018; Huang et al., 2019; He et al., 2020). This implies a further association between splicing alteration via different genotypes and pancreatic islet diseases like diabetes. One example of cell type specific splicing regulation of a pathway is given by BCAR1 which was associated with CDC42 only in alpha cells (Figures 4B, 5A; Supplementary Figure S4A). Reassuringly, it was previously reported to affect cell differentiation via CDC42 (Bauer et al., 2021). By contrast, RAC1 was associated with CDC42 only in beta cells (Figures 4B, 5B; Supplementary Figure S4B), consistent with previous reports linking it to pancreatic islet morphogenesis via beta cells (Greiner et al., 2009). However, none of these three genes were observed in previous bulk pancreatic islet sQTLs (Atla et al., 2022). We investigated CDC42 more carefully and found four SNPs associated with this sQTL from alpha and beta cells. Interestingly, even though they influence the same gene, the perturbation in each cell type was different resulting in varying ratios between isoforms depending on the combination of genetic variants (Figure 5C). This result indicates an additional layer of heterogeneity for sQTLs, which will lead to more diversification in the transcriptome in a cell type specific manner.