A Computational Approach for Defining a Signature of β-Cell Golgi Stress in Diabetes

Diabetes impacts more than 425 million individuals worldwide, and this number is expected to approach 629 million by the year 2045 (1). Type 1 diabetes (T1D) accounts for 5–10% of all diabetes cases and results from autoimmune destruction of pancreatic β-cells, whereas type 2 diabetes (T2D) accounts for 90–95% of cases and results from a combination of peripheral insulin resistance and β-cell dysfunction (2). While the origins of T1D and T2D have been viewed classically as distinct, inadequate insulin secretion plays a central role in the pathophysiology of both forms of diabetes. Accumulating data suggest that common β-cell stress pathways underlie the diminished insulin secretory capacity and reduced β-cell survival observed during T1D and T2D progression.

The primary function of the β-cell is to rapidly sense elevations in blood glucose and respond with carefully titrated levels of insulin secretion (3). To meet the heavy biosynthetic burden of insulin production, β-cells possess a highly developed secretory pathway in which insulin is sequentially folded, processed, and packaged into vesicles for exocytosis. Insulin production begins within the lumen of the endoplasmic reticulum (ER) with conversion of preproinsulin to proinsulin following cleavage of the preproinsulin signal peptide by signal peptidase (3). Proinsulin exits the ER and enters the Golgi apparatus (GA), where it undergoes packaging into early secretory vesicles, and these immature vesicles undergo a series of maturation steps (3). Cleavage of proinsulin by prohormone convertases to produce mature insulin begins within the trans GA subcompartment, and this process is completed within the mature secretory vesicles (3). The mild acidity of the GA aids the maturation of the prohormone convertases needed for insulin maturation (3), and this process continues and accelerates in the more acidic secretory vesicles to form mature insulin granules. The vesicles are then stored near the plasma membrane as a readily releasable pool that is rapidly mobilized during glucose-stimulated insulin secretion (3). As insulin accounts for nearly half of the total protein produced in β-cells (4), disruptions in the fidelity of proinsulin processing and alterations in secretory organelle function can significantly affect β-cell function and health (5,6).

In this regard, ER stress has been well documented to play a role in the pathogenesis of both T1D and T2D (7–10), and the tripartite ER stress response has been studied extensively, resulting in multiple well-validated protein and transcriptional markers of ER stress. The concept of GA stress has been proposed primarily in the context of neurologic disorders (11,12), but has also been suggested to occur in other secretory cells, such as monocytes (13) and Brunner gland cells (14). Several candidate GA stress markers were recently identified in other cell types and linked with GA morphologic changes observed by electron microscopy (15). In contrast with ER stress, the pathways underlying GA stress are incompletely described and may be much more diverse in nature (14). Despite the importance of the GA in insulin processing and maturation, whether GA stress plays a role in β-cells or whether it exists in either T1D or T2D remains unclear.

To this end, we identified publicly available human islet data sets from T1D and T2D models, including a combination of data sets from donors with diabetes and data sets where diabetes pathophysiology was modeled ex vivo, and we analyzed these data sets for changes in the expression of GA-associated genes. In parallel, we treated human islets with the known GA stress inducer brefeldin A (BFA) and used a computational approach to identify overlap between the T1D, T2D, and BFA data sets. From this analysis, we identified a set of genes and pathways associated with β-cell GA stress in diabetes.

Studies of interest were identified by searching specific terms in the National Center for Biotechnology Information Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/gds). Search terms included: “islet and cytokine,” “islet and palmitate,” “islet and type 2 diabetes,” and “islet and type 1 diabetes.” Some data sets were found in more than one search query, but these were weighted equally with data sets appearing in only one query. Data sets were excluded if the primary tissue was not human islets or if there was no nondiabetic control or control condition available for comparison. The experimental procedures of all the publicly available microarray and RNA-sequencing (RNA-Seq) human pancreatic islets samples used in this study have been fully described (16–24).

T1D data sets included cadaveric human islets treated with or without proinflammatory cytokines analogous to autoimmune-mediated β-cell insult in T1D (data sets 1 [18] and 3 [data under review; islets treated with or without 50 units/mL IL-1β and 1,000 units/mL IFNγ for 24 h]) and from organ donors with T1D [data set 2]) (25). T2D data sets included cadaveric human islets from donors diagnosed with T2D (data sets 4 and 5) (17,24) and cadaveric human islets treated with palmitate as a model of T2D (data set 6) (16). The BFA dataset (data set 7) included cadaveric human islets treated with BFA or DMSO as a control. Further information on the data set characteristics, including GEO accession numbers, number of samples, profiling technology, and number of probe sets or genes, is shown in Table 1. Human islet donor characteristics are shown in Table 2.

For the generation of the BFA data set, human cadaveric donor islets were obtained from the Integrated Islet Distribution Program (IIDP) (26). Upon receipt, islets were allowed to recover overnight in phenol red–free DMEM containing 10% FBS, 10 mmol/L HEPES, 2 mmol/L l-glutamate, and 100 units/0.1 mg/L penicillin/streptomycin. Islets were handpicked (250 per condition) and then treated for 24 h with DMSO or 0.1 µg/mL BFA. Islets were collected, washed in PBS, and lysed in TRIzol (Thermo Fisher Scientific, Waltham, MA). Total RNA was isolated using an miRNeasy kit (QIAGEN, Hilden, Germany) for RNA sequencing and RNeasy Micro kit for quantitative RT-PCR (qRT-PCR). Raw data were deposited in the GEO database (GSE152615).

The concentration and quality of total RNA samples were first assessed using an Agilent 2100 Bioanalyzer. An RNA integrity number ≥5 was required to pass quality control (RNA integrity number 8.063 ± 0.4031). Then, 100 ng RNA per sample was used to prepare single-indexed strand-specific cDNA libraries using a KAPA mRNA HyperPrep kit (Roche Sequencing & Life Science, Indianapolis, IN). The resulting libraries were assessed for quantity and size distribution using Qubit and an Agilent 2100 Bioanalyzer. Pooled libraries (300 pmol/L) were sequenced with a 2 × 100 bp paired-end configuration on a NovaSeq 6000 System (Illumina, San Diego, CA). A Phred quality score (Q score) was used to measure the quality of sequencing. More than 90% of the sequencing reads reached Q30 (99.9% base call accuracy).

For microarray data sets, transcriptomic data were background adjusted, summarized, and normalized to probe sets using Robust Multi-Array Average methods from the Bioconductor package affy in R (27) for data sets 4 and 5. Batch correction of data set 5 was performed using the ComBat function from the R/Bioconductor package sva (28). Probe sets were then annotated with their associated gene symbols using annotation GPL files GPL96 (data set 4) and GPL570 (data set 5). Some of the probe sets for each data set had the same gene symbols, and for each of these redundancies, the probe set with the greatest average expression across all samples was chosen to represent each unique gene.

RNA-Seq data sets, including the three GEO data sets (GSE35296, GSE102371, and GSE53949) and the two in-house data sets, were normalized using the reads per kilobase per million mapped reads (RPKM) method. The transcript expressions of all RNA-Seq data were mapped to the hg19 human reference genome using STAR RNA-seq aligner (29). To account for negative values during log transformation in gene RPKM expression, a pseudocount of 1 was added.

Statistical analysis was performed to detect differentially expressed genes (DEGs). For data sets 4 and 5, the empirical Bayes moderated t statistic test was applied using the R/Bioconductor package limma (28), with the P values corrected for age and sex covariates. For data set 7, the analysis was performed using the Bioconductor packages edgeR and limma. DEGs were identified at P values <0.05 and classified as upregulated or downregulated based on the log fold change (FC). For the remaining four data sets, t tests were performed to identify significant DEGs. The P values were corrected for false discovery rate (FDR) using the Benjamini-Hochberg method. The log FC was calculated as the log₂ transformation by dividing the RPKM expression of all genes under treatment/disease conditions by their corresponding RPKM expression under control conditions.

Corrected P values (<0.05) were taken as the threshold for significant differences in gene expression. In addition, a log FC of greater than ±1.5 was necessary to be considered a DEG. A gene was selected as upregulated or downregulated if at least 80% of islet samples exhibited a significant change in expression in one direction and none in the opposite direction. The numbers of normalized transcripts are listed in Table 1.

For all data sets used, the DEGs were filtered to GA-associated genes using the 1,030 genes identified as localized to the GA as defined by The Human Protein Atlas (http://www.proteinatlas.org/humancell/golgi+apparatus). All DEGs are listed in Supplementary File 1. For identification and visualization of the functional profiles for genes and gene clusters, further analyses were performed to identify the significant gene ontology (GO) terms, including biological processes, cellular components, and molecular functions, from the enriched genes in the T1D, T2D, and BFA data sets in relation to GA stress in β-cells, using the Bioconductor clusterProfiler R package (30).

qRT-PCR was performed as previously described (10). For ARF4, ATF3, CREB3, and ACTB, qRT-PCR was performed using SsoAdvanced Universal SYBR Green Supermix (Bio-Rad, Hercules, CA). For COG6, GOSR2, and ACTB, qRT-PCR was performed using TaqMan Fast Advanced Master Mix (Thermo Fisher Scientific). Primers used are listed in Supplementary Table 1.

The data set generated during the current study is available in the GEO repository (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE152615). All data analyzed during this study are included here and in the Supplementary Material.

Following normalization by RPKM, data sets 1–3 were analyzed for DEGs to gain insight into whether β-cell GA stress exists in human models of T1D. Data sets 1 and 3 were from human islets treated ex vivo with proinflammatory cytokines, which is an established model of T1D (29). In contrast, data set 2 included islets from organ donors with T1D. We identified a total of (upregulated/downregulated) 946 (367/579), 1,168 (162/1,006), and 2,434 (934/1,500) DEGs from these three T1D data sets, respectively.

Because overlap between the three data sets was minimal, we used DEG summation for filtering with Golgi genes. Summation resulted in 5,195 unique DEGs (FC >1.5, FDR <0.05), irrespective of profiling technology. After filtering, 611 unique DEGs were found to be GA associated (Fig. 1A). A complete list of DEGs can be found in the Supplementary Material. To simplify functional enrichment and pathway analyses, we applied a more stringent cutoff (FDR <0.01), yielding 173 GA-associated DEGs that were used to elucidate GA-associated pathways involved in T1D.

Functional enrichment yielded 182 pathways, and the top 5 pathways and their associated DEGs were visualized in a cnetplot (Fig. 1B). These top five pathways were “Golgi-associated vesicle,” “Golgi-associated vesicle membrane,” “MHC protein complex,” “cellular response to interferon-γ,” and “response to interferon-γ.” Among these pathways, common GA-associated genes included HLA-DQA1, HLA-DRA, HLA-DPA1, HLA-F, and HLA-B, which have been identified as components of the class I and II MHC and are integral to the autoimmune response during T1D (31). Of note, HLA-F is known to be localized to the GA and not the cell surface like other HLAs (32). We also identified several genes involved in ER-to-GA and GA-to-ER transport, including COPZ2 and KDELR1 (33,34).

We followed a similar workflow to analyze data sets from cadaveric human donors with T2D (data sets 4 and 5) or human islets from donors without diabetes that had been treated ex vivo with palmitate (data set 6). We identified (upregulated/downregulated) 2,466 (884/1,582), 1,722 (931/791), and 300 (80/220) DEGs from the three T2D data sets, respectively. Similar to the T1D data sets, DEG summation was performed and further filtered for GA-associated genes. Summation resulted in 8,969 DEGs (FC >1.5, FDR <0.05) and 1,013 GA-associated DEGs after the Golgi filter was applied (Fig. 2A). A more stringent cutoff (FDR <0.01) was used to simplify the analysis and resulted in 386 DEGs and 390 pathways. Of the 390 pathways, the top 5 were “endomembrane system organization,” “Golgi organization,” “ER-to-Golgi vesicle-mediated transport,” “Golgi vesicle transport,” and “trans Golgi network,” and these pathways were visualized with their associated DEGs in a cnetplot (Fig. 2B). Three GA-associated genes, VTI1A, RAB30, and COG7, were common among the top five pathways and have been identified as important to GA transport and localization (35).

BFA is a known inducer of GA stress (11); however, no GEO-deposited data set is available for BFA-treated human islets. To this end, cadaveric human islets from four donors who did not have diabetes were treated with 0.1 μg/mL BFA for 24 h, followed by RNA sequencing as outlined in 2research design and methods (data set 7). Principle component analysis showed a clear separation between the DMSO control–treated islets and BFA-treated islets (Fig. 3A). Differential expression analysis of this data set identified 8,667 DEGs, with 4,289 genes showing a pattern of upregulation and 4,378 showing a pattern of downregulation (FC >1.5, FDR <0.05). BFA-treated islet DEGs were visualized in a volcano plot (Fig. 3B), with the top 50 DEGs labeled. Functional enrichment of the DEGs revealed 1,034 significant pathways (FDR <0.05), and the top 5 pathways were visualized as a cnetplot (Fig. 3C). Consistent with the notion that BFA is a GA stress inducer, these pathways included “endomembrane system organization,” “ER-to-Golgi vesicle-mediated transport,” “Golgi organization,” “Golgi vesicle transport,” and “post–Golgi vesicle-mediated transport.”

DEGs were then filtered for Golgi-related genes, yielding 503 DEGs, of which 251 were upregulated and 252 downregulated (Fig. 4A). Functional enrichment of 364 DEGs (FDR <0.01) revealed 110 significant pathways, and the top 5 pathways were visualized in an alluvial plot (Fig. 4B). Pathway analysis identified the top five significant pathways as “intrinsic component of Golgi membrane,” “Golgi organization,” “Golgi vesicle transport,” “retrograde vesicle-mediated transport, Golgi to ER,” and “Golgi stack.” Commonly found genes included COG3, COG2, GOLPH3, GBF1, and GOLGA2—all associated with GA localization, trafficking, and integrity (35).

To identify genes that may serve as markers of GA stress in diabetes, we generated sets of genes that overlapped between the T1D and BFA data sets and the T2D and BFA data sets, with and without GA filtering. The overlap of T1D and BFA genes included 2,289 DEGs (Fig. 5A), and after filtering for GA-associated genes, 121 were found to be common (Fig. 5C). The overlap of T2D and BFA genes included 3,540 DEGs (Fig. 5B), and following GA filtering, 204 DEGs were common (Fig. 5D). Overlap among all three groups yielded 55 GA-associated DEGs (Fig. 5E). GO functional enrichment analysis was performed using clusterProfiler R, which analyzed and visualized the functional profiles of genes and gene clusters (33). GO enrichment analysis returned 475 GO terms for T1D, 261 for T2D, 108 for BFA, and 29 common to all three groups (FC >1.5, FDR <0.05). To visualize the overlapping GO terms, we generated a heat map with GA-associated genes and GO terms from our three groups (Fig. 5F). Common GO terms included “ER-to-Golgi vesicle-mediated transport,” “endosome membrane,” “trans Golgi network,” and “transport vesicle.” Examples of unique GO terms included “clathrin-coated endocytic vesicle” and “regulation of TNF production” in T1D, “protein localization to ER,” and “angiogenesis” in T2D, and “synaptic vesicle” and “ARF protein signal transduction” in BFA.

From our BFA-treated islet GA-associated DEGs, we selected five genes for validation as induced by GA stress in islets: CREB3, ARF4, ATF3, COG6, and GOSR2. CREB3, ARF4, and ATF3 were previously reported as GA stress markers (14). COG6 (FC 2.541, FDR 2.16E−05) and GOSR2 (FC 2.251, FDR 8.89E−06) were among the most highly expressed GA-associated DEGs in our BFA data set. COG6 is a subunit of the GA required for maintaining normal GA structure and activity (36), and GOSR2 is a trafficking protein in the medial and trans GA compartments (37).

Next, islets from human cadaveric donors (Table 2), which were distinct from donor islets used for RNA sequencing, were treated with 0.1 μg/mL BFA or DMSO for 24 h. qRT-PCR revealed significantly higher expression of CREB3, ARF4, ATF3, and GOSR2 in BFA-treated islets (Fig. 6A–D). COG6 expression was also increased in BFA-treated islets but not significantly (P = 0.087) (Fig. 6E).

To maintain blood glucose homeostasis and robust levels of insulin production, pancreatic β-cells are subject to a heavy biosynthetic burden that renders them susceptible to influences that disrupt secretory pathway organelle function. Stress within the ER arising from calcium dyshomeostasis, oxidative stress, lipotoxicity, and inflammatory stress are regarded as disruptive to insulin transcription, maturation, and secretion (5,6,8–10). The impact of GA stress on β-cell function, however, is largely unknown and further complicated by the lack of validated GA stress markers. Herein, we sought to address this knowledge gap by using an informatics-based approach using publicly available data sets from T1D and T2D human islet models and our own sequencing data set from human islets treated with the known GA stress inducer BFA.

Evidence from other cell types suggests that the GA is capable of acting as a cell sensor. The GA’s response to stressors that modify or compromise secretory function was initially identified from morphologic assessments, in which the GA was described as having a swollen or fragmented appearance (11,12). GA dysfunction can occur in response to pH modification and extrinsic stressors that interrupt GA-mediated posttranslational modifications or interfere with vesicle transport. In addition, the secretion of large amounts of proteins has been shown to impair GA glycosylation and transport (38), and GA morphologic changes have also been observed after viral infection, such as GA fragmentation and dispersal resulting from herpes virus infection (39). Much of the evidence for GA stress has been discerned from studies of neurodegenerative diseases (11,12,40), but GA stress has also been observed or implicated in skeletal muscle disorders (41), alcohol-induced liver disorders (42), and disorders affecting monocytes (13) and intestinal Brunner gland cells (43). Notably, monocytes and Brunner gland cells are secretory cell types, secreting cytokines (44) and mucins (45), respectively. With emerging evidence indicating that GA stress also affects other secretory cells, it stands to reason that GA stress may also affect β-cells.

Similar to the process of ER stress, GA stress has been linked with increased transcription of genes or activation of signaling pathways that are thought to increase the functional capacity of the GA as a means of restoring organelle homeostasis (11). However, in contrast to ER stress, markers for these responses are limited. Recently, a few candidate markers have been identified in other cell types, including TFE3, HSP47, and CREB3 (14,38). TFE3 is a transcription factor that promotes transcription of GA structural proteins, glycosylation enzymes, and a transport protein (38). HSP47 is an ER chaperone induced by GA stress that prevents GA-induced apoptosis (14,38). CREB3 is a transcription factor that activates ARF4, which is a GA-localized regulator of GA-to-ER vesicle transport. ARF4 activity is impaired by the GA-stress inducer BFA (14,38).

The goal of this study was to identify putative transcriptional signatures of GA stress in diabetic β-cells. To this end, we examined T1D and T2D data sets from publicly available and in-house sources. Differential expression analysis indicated that DEG overlap within each group was minimal; therefore, DEG summation was used for both the T1D and T2D data sets. Each group was then filtered for GA-associated genes using a list of 1,030 genes reportedly specific to the GA according to The Human Protein Atlas. The T1D group included 611 GA-associated DEGs, whereas the T2D group included 1,013 GA-associated DEGs. Following pathway analysis, both the T1D and T2D data sets revealed multiple pathways associated with GA integrity, organization, and transport. In addition, T1D data sets revealed pathways associated with the immune system, including the notable identification of the MHC complex pathway. MHC molecules are processed in the GA before transport to the cell surface, directly for MHC-II and some MHC-I molecules (46,47). Some MHC-I molecules are subsequently transported back to the ER before transport to the cell surface (48). Islet MHC-I is associated with increased ER stress and elevated cell-surface expression in T1D (49). Interestingly, increased MHC-I production and shuttling through the secretory pathway were shown to also induce GA stress in a myositis model of skeletal muscle MHC-I overexpression (41). T2D data sets centered on pathways associated with either GA organization or transport. In this regard, disruption of either ER-to-GA or GA-to-ER transport has been linked with reduced insulin production and secretion and T2D pathophysiology (50,51).

BFA causes GA stress by collapsing the GA, and BFA has been used to prevent cytokine secretion in immune cells (22). In liver and cancerous cells, prolonged BFA exposure has been shown to induce the unfolded protein response and eventually apoptosis (14,52). Within β-cells, BFA reportedly inhibits glucose-stimulated insulin secretion (53–55). To define the overlap between DEGs from T1D and T2D data sets with a known GA stressor, we generated a data set using BFA-treated human islets. To the best of our knowledge, this is the first RNA-Seq data set from BFA-treated human islets. The BFA data set itself had 8,164 DEGs, with 503 associated with the GA. After overlap of GA-associated DEGs from the BFA data set and the T1D and T2D groups, 121 and 204 DEGs were identified, respectively. These results suggest that GA stress may be more prominent in T2D. To further assess overlap, we examined GO terms from the T1D, T2D, and BFA data sets. Twenty-nine GO terms were common to all three groups. Several GO terms were associated with the Golgi complex and involved in mediating transport within the GA or between the GA and other parts of the cell, whereas several others were associated with GA structure.

For validation of our analysis, we selected five highly expressed DEGs from BFA-treated islets based on candidate markers from the literature (CREB3, ARF4, ATF3) and components associated with GA integrity and trafficking (COG6 and GOSR2). Our validation set also identified increased expression of the previously implicated GA stress candidate markers CREB3, ARF4, and ATF3. These genes were differentially expressed in the analyzed data sets: CREB3 in two T2D data sets, ARF4 in one T2D data set, and ATF3 in one T1D data set and two T2D data sets. The conserved oligomeric Golgi (COG) complex is composed of eight subunits and impacts trafficking, glycosylation, and GA integrity (56). COG6 was upregulated in our BFA-treated islets, and COG2 and COG4 were differentially expressed in at least one data set in our diabetic groups, suggesting that expression of components of the COG complex may be altered by GA or diabetes stress. The COG complex has been implicated in BFA-induced GA collapse (57) but not as a marker of GA stress and thus represents an avenue for future investigation. GOSR2 is a member of the GA SNARE transport complex involved in trafficking and glycosylation (37). GOSR2 was differentially expressed in two T2D data sets and our BFA data set, and we observed a trend toward increased GOSR2 expression in our validation set.

Our results showed increased expression of the ATF3 gene, which encodes activating transcription factor 3 (ATF3). ATF3 is a member of the CREB family of transcription factors and is induced by multiple β-cell stressors, including proinflammatory cytokines, free fatty acids, thapsigargin, and nitric oxide (58–61). Notably, ATF3 is induced by the PERK arm of the ER stress response, and ATF3 deletion in β-cells is protective against apoptosis (59,60). These findings highlight the close connection that is likely to exist between ER and GA stress and suggest that ATF3 may serve as a bridge between these processes. Future studies will be needed to delineate the role of ATF3 in GA stress and to test whether changes in ATF3 following GA stress occur independently of ER stress.

Along these lines, several monogenic syndromes associated with diabetes are known to impact β-cell function and survival through alterations in ER health, including Wolfram syndrome, Wolcott-Rallison syndrome, and mutant INS-gene-induced diabetes of youth (MIDY). Wolfram syndrome is caused by a mutation in WFS1, leading to Ca²⁺ leakage from the ER and subsequent ER stress (62). Wolcott-Rallison syndrome is the result of mutations in the EIF2AK3 gene, which encodes for PERK, and loss of functional PERK induces ER stress, reduces insulin secretion, and increases β-cell death (63). MIDY is a syndrome where mutations in the INS gene lead to misfolded proinsulin, resulting in ER stress and atypical proinsulin distribution through the GA (6,64). Maturity-onset diabetes of the young (MODY) subtypes GCK-MODY and ABCC8-MODY are primarily associated with defects in glucose sensing and insulin release (65). However, these genes are also expressed in the GA, and mutations in these genes could contribute to GA dysfunction. While these are rare, the study of preclinical models as well as tissues from human donors with these monogenic subtypes of diabetes could provide insight into the mechanisms of GA stress as well as the intersection between ER and GA dysfunction.

In summary, our analysis revealed that multiple GA-associated genes are differentially expressed in T1D and T2D. Overlapping the T1D and T2D groups with the BFA data set, we identified 120 and 204 DEGs, respectively. The higher number of DEGs in T2D may be due to the different mechanisms by which T1D and T2D develop and suggests that GA stress may be more prevalent in T2D, but further investigation is needed. In BFA-treated islets, we validated a signature of GA stress, including ATF3, ARF4, CREB3, and COG6. Several members of the COG complex were differentially expressed across data sets and thus may be a new target in investigations of GA stress. Taken together, our data indicate that GA-associated genes are dysregulated in diabetes, potentially opening an avenue of research to identify new markers of GA stress.

Funding. This work was supported by National Institute of Diabetes and Digestive and Kidney Diseases (NIDDK) grants R01 DK093954 and UC4 DK104166 (to C.E.-M.), U.S. Department of Veterans Affairs Merit Award I01BX001733 (to C.E.-M.), JDRF grant 2-SRA-2018-493-A-B (to C.E.-M.), and gifts from the Sigma Beta Sorority, the Ball Brothers Foundation, and the George and Frances Ball Foundation (to C.E.-M.). The authors acknowledge the support of the Indiana Diabetes Research Center Islet & Physiology Core (P30 DK097512). R.N.B. was supported by a National Institute of Allergy and Infectious Diseases, National Institutes of Health (NIH), training grant (T32 AI060519) and by a JDRF postdoctoral fellowship (3-PDF-2017-385-A-N). F.S. was supported by a JDRF postdoctoral fellowship (3-PDF-2016-199-A-N). Human pancreatic islets were provided by the NIDDK-funded Integrated Islet Distribution Program (IIDP) (RRID:SCR_014387) at City of Hope, NIH grant no. 2UC4DK098085.

The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Duality of Interest. No potential conflicts of interest relevant to this article were reported.

Author Contributions. R.N.B. directed the design of the study, the data analysis and interpretation, the collection and assembly of data, and the writing of the manuscript. O.O., S.T., S.S., and P.K. participated in data analysis and interpretation and critical revision of the manuscript. F.S. participated in data collection and critical revision of the manuscript. H.W. participated in study conception and design and critical revision of the manuscript. C.E.-M. directed funding acquisition and study conception and design, participated in the collection and assembly of data, contributed to data analysis, directed the writing of the manuscript, and gave final approval of the manuscript. C.E.-M. is the guarantor of this work and, as such, had full access to all the data in the study and takes responsibility for the integrity of the data and the accuracy of the data analysis.