Type 2 diabetes candidate genes, including PAX5, cause impaired insulin secretion in human pancreatic islets

Transcriptomic analysis of human pancreatic islets identified T2D candidates.

To identify regulators of insulin secretion that may contribute to T2D, we generated RNA-Seq data from the large Lund University Diabetes Centre (LUDC) pancreatic islet cohort, comprising islets from a total of 309 donors (Supplemental Table 1; supplemental material available online with this article; https://doi.org/10.1172/JCI163612DS1). From these, RNA from 283 islet samples was successfully sequenced (Figure 1). Before further analysis and to ensure the highest possible quality of data, we filtered our islet samples (Figure 1). ND controls were defined as individuals lacking a diabetes diagnosis and exhibiting HbA1c levels below 42 mmol/mol. Moreover, age of 40 years or older was used as an inclusion criterion for the ND controls, since none of the individuals with T2D were younger than 40 years of age. The characteristics of the individuals in the LUDC islet case-control cohort are summarized in Table 1 and Supplemental Figure 1, A and B. We subsequently compared RNA-Seq data from islets from 138 ND controls with those from 33 individuals with T2D (Figure 1) using a generalized linear model with correction for age, sex, islet purity, and days in culture (DIC). The analysis identified 395 DEGs (FDR below 5%, q < 0.05, Supplemental Table 2, sheet A), which included 228 upregulated and 167 downregulated genes. We further applied a model in which we also added BMI as a covariate; expression of 394 of the 395 DEGs was then associated with T2D (P = 8.6 × 10^–17 to 4.5 × 10^–2, Supplemental Table 2, sheet A). Additionally, we performed an analysis including only donors with an islet purity of 80% or higher (n = 69 ND and n = 18 T2D); expression of 334 of 395 DEGs was then associated with T2D (P = 6.7 × 10^–11 to 4.9 × 10^–2, Supplemental Table 2, sheet A).

For replication, we next compared our 395 identified DEGs with published expression data from studies of human pancreatic islets from T2D case-control cohorts (1–3, 5–8). Previous bulk expression studies identified 75 of the 395 DEGs identified here (Figure 2A and Supplemental Tables 2 and 3) (1, 2, 6, 7). Additionally, published scRNA-Seq studies found differential expression of 28 of our 395 DEGs in islet cells from ND controls versus individuals with T2D (Figure 2B and Supplemental Tables 2 and 3) (3, 5, 8). Eight and 9 of these DEGs were differentially expressed in α and β cells, respectively (Supplemental Table 2, sheet A), and their expression in islets is presented in Supplemental Figure 1, C and D. Taken together, 94 of our 395 DEGs have been identified in one or more of these previous studies, whereas 301 were not (Supplemental Table 2, sheet A). ARG2, GLRA1, IAPP, IGFBP2, PDHX, PPP1R1A, PTEN, RASD2, SMAD9, SYT13, TBC1D4, and UNC5D are DEGs previously identified by the studies included in Figure 2, A and B (Figure 2C and Supplemental Table 2, sheet A), while, for example, CDKN1C, GABRA1, GABRA2, GAD1, IGFBP4, IGFBP6, IL6, PDE7B, SIRT1, SLC2A5, SOCS1, SOCS6, SYT1, SYT12, SYT14, and TET1 are DEGs not identified by those studies (Figure 2D and Supplemental Table 2, sheet A).

To better understand which islet cell types expressed our 395 DEGs, we studied their expression in FACS human α and β cells from the LUDC sorted α/β cell cohort. This cohort included 13 ND individuals, 2 individuals with T2D, and 3 individuals with prediabetes (Table 1). We found that α and β cells expressed 366 and 368 of our DEGs, respectively (Figure 2E and Supplemental Table 2, sheet A). Of these, 90 genes showed higher expression in α versus β cells, while 174 had higher expression in β versus α cells (q < 0.05, Supplemental Table 2, sheet A). This included 9 of the 11 genes that we selected for functional validation (see below and Figure 2F). Using WebGestaltR, we found 17 enriched gene sets (q < 0.05), including divalent inorganic cation homeostasis and positive regulation of cell adhesion, among the DEGs with higher expression in α versus β cells (Supplemental Figure 1E), while there was no enrichment among the DEGs with higher expression in β versus α cells. We also examined the expression of our DEGs in different islet cell types using published human islet scRNA-Seq data (5). We found that among our 395 DEGs, 366 (93%) were expressed in α cells, 347 (88%) in β cells, 331 (84%) in pancreatic polypeptide (PP) cells, 312 (79%) in delta (δ) cells, and 190 (48%) in epsilon (ε) cells (Supplemental Table 4 and Supplemental Figure 2). We found that 97%–98% of DEGs expressed in α and β cells in the scRNA-Seq expression data overlapped with those expressed in α and β cells from the LUDC sorted α/β cell cohort. We further tested whether any of the 395 islet DEGs also exhibited differential expression in the same direction in sorted α or β cells when comparing cells from 5 donors with T2D or prediabetes versus those from 13 ND controls (Table 1). Here, 4 and 28 of our 395 DEGs showed differential expression (P < 0.05) in sorted α or β cells, respectively (Supplemental Table 2, sheets B and C, Supplemental Figure 3, A and B). The DEGs in sorted β cells included BARX1, NEFL, PAX5, and PCOLCE2, genes we later selected for functional follow-up (see below). It should be noted that, given the modest sample size, there was limited power when comparing T2D with ND samples from both the scRNA-Seq and sorted RNA-Seq data, and hence one would not expect to find all 395 DEGs in these analyses. These data suggest that several DEGs in T2D islets may affect insulin secretion or other aspects of β cell function.

Our data clearly show that individuals with T2D exhibit abundant transcriptomic alterations in pancreatic islets. However, it is not clear whether these changes predispose individuals to T2D or whether they are a consequence of disease progression. To determine whether changes in expression levels of the 395 identified DEGs potentially predispose individuals to T2D, we studied their expression in human islets from the LUDC islet HbA1c cohort (Table 1), including individuals with a wide range of HbA1c levels who were not previously diagnosed with diabetes. Interestingly, expression of 142 of the 395 DEGs (36%) was linearly associated with HbA1c levels (Supplemental Table 2, sheet A). Among those genes, all had a β coefficient with the same directionality as that of the expression differences seen in islets from the LUDC islet case-control cohort. These data suggest that the changes in expression may occur before the diagnosis of T2D and potentially contribute to the development of the disease. They also indicate that the expression changes were not due to T2D treatment.

Gene sets, including hormone secretion, are enriched among T2D-associated DEGs.

We used WebGestaltR and gene ontology to discover enriched gene sets among our 395 DEGs. We found 39, 5, and 6 enriched gene sets for biological processes, cellular components, and molecular functions, respectively (q < 0.05, Supplemental Table 5). Interestingly, gene sets for growth factor binding, regulation of lipid transport, negative regulation of Wnt signaling, aging, regulation of cell-cell adhesion, hormone secretion, regulation of inflammatory response, and divalent inorganic cation homeostasis were among the significantly enriched gene sets (Figure 2G).

T2D-associated DEGs display alterations in open chromatin and DNA methylation in human islets.

Open chromatin regions are associated with active gene transcription (15). Of our 395 DEGs in islets from ND controls versus those from individuals with T2D, 346 (88%) were marked by at least 1 open chromatin region, identified by previous ATAC-Seq of human pancreatic islets (15). From this subset, 194 and 152 DEGs were upregulated and downregulated, respectively, in islets from individuals with T2D. The top DEGs BARX1, OPRD1, PAX5, and PCOLCE2 overlapped with several open chromatin regions (4, 4, 5, and 6, respectively) located both upstream and downstream of the respective transcription start sites (Supplemental Table 6, sheet A). Moreover, 24 DEGs, including CAMK4, DIO2, DKK3, FOXP1, GABRA2, PTPRC, SOCS1, and SYNPO, displayed alterations of open chromatin in islets from ND controls versus those from individuals with T2D. Of these, 23 exhibited more open chromatin and higher expression, and 1 (GABRA2) showed less open chromatin and reduced expression, in islets from donors with T2D (Figure 2H, Supplemental Table 6, sheet A). Together, these data suggest an altered chromatin state in islets from individuals with T2D that was associated with altered gene expression for a subset of identified DEGs.

We further tested whether the 395 DEGs showed altered DNA methylation in islets from individuals with T2D versus those from ND controls based on published data (16) and found 732 differentially methylated regions (DMRs) annotated to 262 DEGs (Supplemental Table 6, sheet B). These included DMRs annotated to CHL1, ELFN1, FAIM2, HHATL, OPRD1, PAX5, SFRP1, and SLC2A2, DEGs selected for functional follow-up (see below).

SNPs associate with the identified DEGs, T2D, and metabolic traits.

We next assessed whether SNPs were associated with expression of the 395 DEGs. We first explored expression quantitative trait loci (eQTL) data from human pancreatic islets from the INSPIRE (Integrated Network for Systematic analysis of Pancreatic Islet RNA Expression) study (17) and the T2D Genes portal (https://t2d.hugeamp.org/). We found 148 SNPs associated with islet expression of 120 of our 395 DEGs, including CHL1, FAIM2, GABRA2, HHATL, PCOLCE2, SFRP1, SIRT1, and SMAD9 (Supplemental Table 7, sheet A, and Figure 2I). We then used the GWAS catalog to determine whether any of these 148 SNPs have previously been linked to T2D or other metabolic diseases or traits (Supplemental Table 7, sheet A). Indeed, 2 eQTL SNPs associated with FXYD2 (rs529623) and RPL39L (rs3887925) islet expression have been linked to T2D risk (12) (Supplemental Table 7, sheet A, and Figure 2I). Moreover, 2 eQTL SNPs associated with islet expression of FOXE1 (rs7038480) and ENTR1 (rs11145930) have been linked to blood glucose levels (12, 18–20) (Supplemental Table 7, sheet A, and Figure 2I). Four eQTL SNPs associated with islet expression of ARPC1B (rs3843540), COMP (rs7260000), DIXDC1 (rs10891295), and HSD3B7 (rs4889599) have been linked to BMI and/or waist-to-hip ratio (18, 21) (Supplemental Table 7, sheet A, and Figure 2I). Finally, 6 eQTL SNPs associated with islet expression of ACP2 (rs11039194), CBLC (rs8104467), CD5 (rs7124430), HSD3B7 (rs4889599), PCOLCE2 (rs6794287), and TMED6 (rs113671952) have been linked to triglyceride or LDL cholesterol levels (18, 22–24) (Supplemental Table 7, sheet A, and Figure 2I). These data demonstrate that SNPs associated with islet expression of numerous identified DEGs impact the risk for T2D and metabolic traits (Figure 2I).

As our islet cohort data included variables, e.g., islet purity, that were not available or used in the INSPIRE study (17), and since we used a cutoff for islet purity, we also performed an eQTL analysis adjusting for age, sex, BMI, DIC, purity, and T2D in the LUDC islet case-control cohort. This showed that 26 of our 395 DEGs, including BEST3, HSD3B7, and TMED6, had at least 1 eQTL, based on a q value of less than 0.05. Moreover, 374 DEGs had eQTLs with nominal P values (P < 0.05, Supplemental Table 7, sheet B), and these included 47 of the 148 SNPs and 117 of the 123 genes identified in the INSPIRE analysis (Supplemental Table 7, sheets A and B).

We then assessed whether SNPs that map to any of the 395 DEGs have been associated with T2D and/or glycemic traits (HbA1c, fasting glucose, fasting insulin, homeostatic model assessment of β cell function [HOMA-B], disposition index [DI], or corrected insulin response [CIR]) in GWAS using the Common Metabolic Diseases Knowledge Portal (hugeamp.org, accessed August 2022). A total of 106 DEGs were linked to 149 unique SNPs associated with these traits in GWAS (Supplemental Table 7, sheet C, Figure 2I); 131 SNPs were associated with T2D, 38 with HbA1c, 10 with fasting glucose, 1 with fasting insulin, 4 with DI, and 1 with CIR. This included 1 SNP mapping to each of BARX1, CHL1, ELFN1, and FAIM2, three SNPs mapping to SLC2A2, and 4 SNPs mapping to SFRP1, DEGs we subsequently selected for functional follow-up (see below).

Since genetic and epigenetic mechanisms interact and together affect biological processes, we tested whether SNPs in cis are associated with DNA methylation of CpG sites annotated to the 395 DEGs in human pancreatic islets, so-called mQTLs (25). We found 490 SNPs associated with DNA methylation of 176 unique sites annotated to 90 DEGs (Supplemental Table 7, sheet D). These included mQTLs, 3 each, annotated to OPRD1 and PAX5, two additional DEGs we subsequently selected for functional follow-up (see below).

Panels A–C in Supplemental Figure 4 show the overlap between DEGs that had eQTLs, DMRs, GWAS SNPs, ATAC-Seq peaks, or mQTLs annotated, and the overlap between the SNPs in the different genetic analyses above.

The identified T2D-associated DEGs affect metabolism in vivo.

To explore potential in vivo evidence for protective or susceptible functions in the development of diabetes for the 395 identified T2D DEGs, we systematically searched the IMPC database (https://www.mousephenotype.org/, accessed April 15, 2020). The IMPC database contains KO mouse strain entries for 168 of the 395 DEGs, with phenotype data available for 125 strains (Figure 3A). Most were homozygous viable (n = 85), with fewer that were homozygous subviable (meaning that the number of alive homozygous pups was lower than expected, n = 6), homozygous lethal but heterozygous viable (n = 24), or of unknown viability (n = 10). Supplemental Table 8 lists DEGs found in T2D islets with IMPC metabolic readouts for KO strains (n = 125). Of the phenotyped KO strains, insulin blood (IB) measurement data were available for 36 KO strains, with 16 (44%) displaying altered levels (P < 0.05, Figure 3 and Supplemental Table 8). Intraperitoneal glucose tolerance test (IPGTT) data showed that the fasted blood glucose concentration (FG) and AUC for glucose (AUCG) were altered in 24 of 81 (30%) and 24 of 83 (29%) KO strains, respectively (P < 0.05), while the initial response to glucose (IGR) challenge was altered in 26 of 80 (33%) (P < 0.05) (Figure 3 and Supplemental Table 8). With respect to body composition, dual-energy x-ray absorptiometry (DEXA) phenotypic data were available for 93 KO models, and among these, 33 (35%) and 35 (38%) had altered total fat mass and total lean mass, respectively (P < 0.05, Figure 3 and Supplemental Table 8).

The effects of gene KO on the studied metabolic phenotypes agreed with the direction of altered expression for several of the T2D-associated DEGs. For example, our data show that HHATL expression was significantly decreased in T2D islets (Supplemental Table 2). In line with this, IMPC data showed that Hhatl-KO mice had impaired glucose tolerance, increased total body fat mass, and decreased total body lean mass (Figure 3B). SLC2A2 expression was also decreased in islets from individuals with T2D (Supplemental Table 2), and Slc2a2-KO mice had increased FG levels and IGR. Overall, these rodent in vivo findings indicate an important role for the identified DEGs in the control of glucose homeostasis and body composition.

Functional follow-up shows that changes in expression of T2D-associated DEGs impair β cell function.

We then asked whether any DEGs previously unrecognized in comparable studies have a functional role in β cells. First, we used a strategy to select genes for functional investigation (Figure 4A). We selected DEGs with a greater than 2-fold change in either direction, based on mean expression levels in islets in ND controls versus individuals with T2D, and the same directionality of change in expression as the HbA1c correlation in islets from the LUDC islet HbA1c cohort. This generated a list of 31 genes. We further excluded DEGs that were not expressed in the sorted α or β cells from ND controls, unless their expression was induced in islets from individuals with T2D, and that were not expressed in endocrine cell types in the scRNA-Seq data sets. Finally, on the basis of the lowest q values, we selected the top 9 genes that, to our knowledge, had not previously been studied in β cells (Figure 4A). In addition, we included the genes CHL1 and SLC2A2. Knockdown of Chl1 in clonal β cells was previously found to impair glucose-stimulated insulin secretion (GSIS) (26), and we included this gene as a positive control. SLC2A2, encoding the glucose transporter GLUT2, has also been studied, but we included it as it is debated whether its function is significant in human β cells (27, 28). Of note, 10 of these 11 genes have either T2D-associated DMRs (Supplemental Table 6, sheet B), INSPIRE eQTLs (Supplemental Table 7, sheet A), SNPs in cis associated with T2D and/or glycemic traits (Supplemental Table 7, sheet C), or mQTLs (Supplemental Table 7, sheet D), as summarized in Figure 4B. Moreover, all 11 DEGs have nominal LUDC eQTLs (P = 2.7 × 10^–4 to 2.3 × 10^–2, Supplemental Table 7, sheet B).

To mimic the situation in individuals with T2D, the selected genes that had lower expression levels in islets in T2D, i.e., CHL1, HHATL, OPRD1, and SLC2A2 (Figure 4C), were silenced by siRNA transfection in human islets from ND individuals. This resulted in a 50%–60% reduction in gene expression (Figure 4D). Determinations of insulin secretion showed that islets deficient for CHL1 displayed a reduction in insulin release when exposed to basal glucose (2.8 mM), whereas islets deficient for SLC2A2 showed nominally reduced insulin secretion at the same glucose concentration (P = 0.062, Figure 4E). At stimulatory glucose levels (16.7 mM), OPRD1 or SLC2A2 silencing resulted in approximately 20% and 40% reduced insulin secretion, respectively, while islets deficient for CHL1 or HHATL showed nominally reduced and unaffected secretion, respectively (Figure 4E). The changes did not translate into significant differences in the fold change of insulin secretion (secretion at stimulatory glucose divided by secretion at basal glucose levels) and occurred without differences in insulin content (Supplemental Figure 5, A and B). We also measured glucagon secretion from islets in which CHL1, HHATL, OPRD1, or SLC2A2 had been silenced and found no differences (Supplemental Figure 5C). These data indicate that reduced islet expression of OPRD1 and SLC2A2 may contribute to the insulin secretion defect seen in individuals with T2D.

Selected DEGs that exhibited higher islet expression in T2D, i.e., BARX1, ELFN1, FAIM2, NEFL, PAX5, PCOLCE2, and SFRP1 (Figure 4C), were overexpressed by plasmid transfection in 832/13 INS1 β cells (hereafter called INS1 β cells). We used INS1 β cells for overexpression experiments because of the limited supply of human islets. Furthermore, in our hands, INS1 β cells behaved more like primary mature human β cells than did the fetal human β cell line EndoC-βH1 in gene manipulation experiments (KB and JKO, unpublished observations). Plasmid transfection resulted in impaired insulin secretion in absolute numbers at basal or stimulatory glucose levels, reduced insulin secretion expressed as fold change, or altered insulin content in cells overexpressing Barx1, Nefl, Pax5, Pcolce2, or Sfrp1 (Supplemental Figure 5, D–G). However, the plasmid transfection efficiency was low (20%–30% of cells, Supplemental Figure 5, H–J). This, together with the fact that expression was driven by the very strong CMV promoter, made us wary that the phenotypes may be due to very high overexpression overwhelming the transfected cells. We therefore proceeded to use lentiviral vectors, allowing for transgene delivery to a larger proportion of cells, with the weaker human phosphoglycerate kinase 1 promoter driving expression of HA-tagged Barx1, Nefl, Pax5, Pcolce2, or Sfrp1. Western blot analysis showed that transduction of INS1 β cells with these vectors resulted in expression of proteins of expected sizes (Figure 4F). Overexpression of Pax5 significantly increased basal secretion and strongly reduced insulin secretion at stimulatory glucose levels, without influencing insulin content (Figure 4, G and H), whereas Nefl or Pcolce2 overexpression caused reduced insulin secretion at basal (both) and stimulatory (Pcolce2) glucose levels (Figure 4G) and increased insulin content (Figure 4H). Overexpression of Barx1 or Sfrp1 did not alter insulin secretion or content. The changes described caused a strong reduction in fold change of insulin secretion in Pax5-overexpressing INS1 β cells (Figure 4I). Finally, wells where cells were transduced with the Pax5 or Pcolce2 viruses also contained significantly less total protein than wells with cells transduced with the GFP virus (Figure 4J). This indicates a lower cell number and, hence, suggests that Pax5 and Pcolce2 may affect cell viability and/or proliferation and therefore potentially β cell mass in T2D. Taken together, these functional experiments identified 4, to our knowledge, previously unrecognized regulators of insulin secretion showing either lower (OPRD1) or higher (NEFL, PAX5, and PCOLCE2) expression in islets from individuals with T2 (Figure 4, C–J). Because Pax5 overexpression had the strongest effect on both insulin secretion and protein content (Figure 4, G and J), we decided to further explore its impact on β cells. First, to exclude the possibility that the HA tag renders Pax5 pathogenic, we repeated the secretion experiments with viral vectors conferring expression of untagged Pax5. The results showed that also untagged Pax5 reduced GSIS in INS1 β cells (Supplemental Figure 5K).

Overexpression of Pax5 perturbs mitochondrial function and causes β cell loss.

PAX5 encodes a transcription factor associated with leukemia (29) and, to our knowledge, has not been studied in β cells. First, we used IHC to assess the distribution of PAX5 in human islets. In accordance with the RNA-Seq data, we found that islets from individuals with T2D displayed robust PAX5 expression, with most of the staining found in β cells, whereas PAX5 was barely detectable in islets from ND controls (Figure 5A).

To further characterize the Pax5-induced defects in INS1 β cells, we stimulated insulin secretion with a depolarizing concentration of K⁺. These experiments again showed that GSIS was severely blunted by Pax5 overexpression, while K⁺-stimulated secretion was increased (Figure 5B). This indicates that secretory defects induced by Pax5 overexpression occurred before depolarization of K_ATP channels in the insulin secretion pathway, potentially in mitochondrial metabolism. To determine whether Pax5 overexpression alters mitochondrial metabolism in β cells, we used the Seahorse XF analyzer to measure the oxygen consumption rate (OCR). This showed that Pax5 overexpression caused a reduction in mitochondrial respiration with significantly lower glucose-stimulated respiration, both in terms of increase compared with basal respiration and of absolute values, as well as lower maximal respiration (Figure 5, C–F). These respiratory defects were reflected in a reduced PercevalHR signal, demonstrating a lower ATP/ADP ratio in response to glucose stimulation in Pax5-overexpressing β cells (Figure 5, G–I). To find a possible cause for these mitochondrial defects, we performed Western blotting to measure protein levels of subunits of complex I–V in the electron transport chain and of the citric acid cycle enzyme citrate synthase. We found that the levels of citrate synthase and Sdhb (complex II) were reduced by Pax5 overexpression, while Uqcrc2 (complex III) was slightly increased (Figure 5, J–O, and Supplemental Figure 5, L and M). Together, these data support the idea that increased Pax5 levels perturb insulin secretion by negative effects on mitochondrial function.

In view of the data shown in Figure 4J, we hypothesized that elevated Pax5 levels may affect cell numbers. This was supported by results from an MTT [3-(4,5-dimethylthiazol-2-yl)-2,5-diphenyltetrazolium bromide] assay, which measures metabolic activity and is used to investigate cell viability and numbers (Figure 6A). To investigate whether the reduced cell number is due to increased apoptosis and/or decreased proliferation, we quantified caspase-3/-7 activity and cleaved (i.e., active) caspase-3 as a measure of apoptotic activity, and used an 5-ethynyl-2′-deoxyuridine (EdU) incorporation assay to quantify proliferation. These analyses showed that Pax5 overexpression induced a clear increase in apoptotic activity and a drop in proliferation rates in INS1 β cells (Figure 6, B–D). Hence, we concluded that Pax5 is a regulator of β cell viability and numbers.

In a final effort to characterize the mechanisms underlying the profound effects of Pax5, we performed a global transcriptomic analysis of Pax5-overexpressing INS1 β cells. The analysis revealed that 3,069 genes were differentially expressed when comparing β cells overexpressing Pax5 and GFP, respectively (q < 0.05, Supplemental Table 9). These included 75 of the 395 T2D DEGs (Figure 6E and Supplemental Table 2, sheet A, and Supplemental Table 9). In addition to Pax5 itself, 3 of the other DEGs selected for functional follow-up, Faim2, Pcolce2, and Slc2a2, were altered in Pax5-overexpressing INS1 β cells, and in the same direction as in islets from individuals with T2D (Figure 6E). The 3,069 DEGs were enriched for many gene sets, including those for insulin secretion, glucose homeostasis, response to glucose, positive regulation of cell death, and aging (Supplemental Table 10 and Figure 6F). Overall, these data clearly demonstrated that Pax5 overexpression leads to transcriptomic changes that can have profound effects on β cell function and survival.

PAX5 is a potential dysregulator of gene expression in T2D.

As PAX5 is a transcription factor, we next used the bioinformatics prediction tool Pscan (30) to test whether the 395 DEGs are enriched for genes with a PAX5-binding motif in the promoter. This would suggest that elevated PAX5 may cause dysregulation of other DEGs in T2D. Indeed, our analysis revealed that 196 DEGs had a putative PAX5-binding motif within the promoter, which is a significant enrichment (P = 0.017, Supplemental Table 11). Of note, the 196 genes included 6 of the genes we selected for functional follow-up (BARX1, CHL1, FAIM2, HHATL, OPRD1, and PAX5 itself), and a literature search showed that approximately 25% of the 196 genes have been found to regulate β cell function and/or numbers, or have genetic variants associated with T2D or other metabolic traits in humans (Supplemental Table 11 and Figure 7A). These include IL6 (31), PPP1R1A (32), PTEN (33), SYT13 (34), and TBC1D4 (35). Additionally, in human islets, the expression of PAX5 correlated with 126 (32%) of the other 394 DEGs (Spearman’s correlation, P < 0.05, Supplemental Table 12), including BARX1, CDKN1C, CHL1, ELFN1, GABRA2, GAD1, NEFL, PCOLCE2, PDE7B, SFRP1, SLC2A2, SOCS1, SYT1, and SYT12. As a final piece of evidence supporting PAX5 as a key DEG, we used weighted correlation network analysis (WGCNA) (36) for a coexpression analysis to establish whether PAX5 is part of a dysregulated gene network in T2D. This analysis revealed a network with PAX5 and 86 other T2D DEGs (Supplemental Table 13 and Figure 7B). Together, data from human islets and Pax5-overexpressing INS1 β cells suggest that PAX5 may cause dysregulation of many identified T2D-associated DEGs.

Human islets.

Human islets of Langerhans were obtained from the Human Tissue Laboratory, which is funded by the Excellence of Diabetes Research in Sweden (EXODIAB) network (www.exodiab.se/home) in collaboration with the Nordic Network for Clinical Islet Transplantation Program (www.nordicislets.org). Islets from cadaver donors were prepared by enzymatic digestion and density gradient separation. Islet preparation purity and counts were determined as described previously (69). Details on islet purity, measured by dithizone staining of individual islet preparations, are presented in Supplemental Table 14.

RNA-Seq.

RNA was extracted from human pancreatic islets using the AllPrep DNA/RNA kit or the miRNeasy Mini Kit (QIAGEN). Sample preparation of 1 μg high-quality RNA was done using the TruSeq RNA Library Preparation Kit or the TruSeq Stranded Total RNA Library Prep kit, followed by RNA-Seq on a HiSeq 2000 or NextSeq 500 (Illumina), respectively.

Data and material availability.

All data needed to evaluate the conclusions in this study are presented here and/or in the Supplemental Methods. The human islet RNA-Seq data sets were deposited in the LUDC repository (https://www.ludc.lu.se/resources/repository) under the following accession numbers: LUDC islet case-control cohort (accession no. LUDC2022.07.111), LUDC sorted α/β cell cohort (accession no. LUDC2022.07.112), and LUDC islet HbA1c cohort (accession no. LUDC2022.07.113). (These LUDC data are available upon request.) Data have also been deposited in the European Genome-Phenome Archive (https://www.ebi.ac.uk/ega/) under accession number EGAD00001005512. The ATAC-Seq data set and mRNA expression data from GFP- and Pax5-overexpressing INS1 β cells are available in the NCBI’s Gene Expression Omnibus (GEO) repository under accession numbers GSE129383 and GSE211310, respectively.

Statistics.

The functional data were analyzed with a 2-tailed, paired t test unless stated otherwise. Other analyses were performed as described in the respective sections in the Supplemental Methods. Data are presented as the mean ± SEM, unless stated otherwise. A P or q value (when correction for multiple testing was required) below 0.05 was considered statistically significant.

Study approval.

Written informed consent was obtained from pancreatic donors or their relatives, and all procedures were approved by the Swedish Ethical Review Authority (Permit number 2011263).

Version 1. 01/19/2023

In-Press Preview

Version 2. 02/15/2023

Electronic publication

Version 3. 12/21/2023

Figure 4D-E corrected to one star for the statistical comparisons