Single-Cell Transcriptomics Links Loss of Human Pancreatic β-Cell Identity to ER Stress

The maintenance of pancreatic islet architecture is crucial for proper β-cell function. We previously reported that disruption of human islet integrity could result in altered β-cell identity. Here we combine β-cell lineage tracing and single-cell transcriptomics to investigate the mechanisms underlying this process in primary human islet cells. Using drug-induced ER stress and cytoskeleton modification models, we demonstrate that altering the islet structure triggers an unfolding protein response that causes the downregulation of β-cell maturity genes. Collectively, our findings illustrate the close relationship between endoplasmic reticulum homeostasis and β-cell phenotype, and strengthen the concept of altered β-cell identity as a mechanism underlying the loss of functional β-cell mass.


Introduction
Type 2 diabetes mellitus (T2D) is a complex metabolic disorder caused by failure of pancreatic β-cells in the presence of peripheral insulin resistance. Studies in mouse models led to the concept of β-cell dedifferentiation as one of the mechanisms of β-cell dysfunction in diabetes [1,2]. While most studies on β-cell dedifferentiation rely on the forced activation or deletion of β-cell transcription factors in murine models, we and others found evidence of alterations in β-cell identity in islets from patients diagnosed with type 2 diabetes, such as an increased frequency of polyhormonal cells, a reduced expression of key β-cell transcription factors (e.g., MAFA, PDX1) as well as an increased proportion of degranulated islet cells [3][4][5] β-cell function (assessed on isolated islets) was found to be inversely correlated with the higher proportion of hormone-negative cells in islets from T2D individuals [5]. Similarly, the higher frequency of islet cells presenting an altered identity was correlated with the presence of islet amyloid deposits, which are associated with a reduced functional β-cell mass in type 2 diabetes [6]. Collectively these findings point towards a causative link between β-cell identity change and reduced functional β-cell mass.
In vitro models for β-cell identity changes in primary human islets are limited. We recently reported the downregulation of β-cell maturity markers, particularly MAFA, in drug-induced diabetes [7]. In addition, we previously described a model, in which disruption of primary human islet integrity triggered severe phenotypic alterations in β-cells including the conversion of part of insulin-producing β-cells into glucagon-positive α-cells [8]. Importantly, this process could be prevented by inactivation of ARX, leading to preserved β-cell identity and function.
In recent years, single-cell RNA sequencing (scRNAseq) has emerged as a powerful tool to study complex cell populations including the human pancreas [9][10][11][12][13]. Here, we combine lentivirus-mediated lineage tracing and single-cell transcriptomics to investigate the molecular mechanisms underlying the response of primary human β-cell to cellular stress.

Cell Culture
Islets were isolated from donor pancreas allocated (after anonymisation) by Eurotransplant for the clinical islet transplantation program of our institute (Leiden University Medical Center). Islets were used for research only if they could not be used for clinical purposes, and if research consent was present, according to Dutch national laws. Human islets preparations from 35 non-diabetic donors were used for this study (Table S7) [14]. Islets with a purity of at least 75% were cultured as previously described [8].

Islet Integrity Disruption Model
Two distinct experimental setups were performed to investigate the β-cell response to stress.
Experimental setup 2: Intact islets were transduced overnight with the two lentiviral vectors followed by tamoxifen exposure, as described above. After 3-4 days, transduced islets were dispersed and prepared for scRNAseq ( Figure S4A).

Other Experimental Models
Human islets or EndoC-βH1 cells were treated with either 0.1 µM thapsigargin (TG, Bio-Connect, Huissen, the Netherlands) for 24 h or 1 µM TG for 5 h. Read-outs were performed 24 h after the start of the treatment with TG, unless stated otherwise in the figure legends.

Cell Sorting and scRNAseq
Islets and aggregates cultured as outlined in experimental setup 1 and 2 were briefly washed in PBS and dispersed into single cells. Cells were sorted using a FACS Jazz or FACS Aria III (BD Biosciences, Franklin Lakes, NJ, USA). Live single cells were selected based on DAPI exclusion and FSC ( Figure S1B), sorted into 384-well hard shell plates (BioRad, Hercules, CA, USA) and further processed with SORTseq (SOrting and Robot-assisted Transcriptome sequencing) as previously described [9].
Islets from donors 1, 2 and 3 were used in experimental setup 1, and donors 4, 5 and 6 for setup 2. Of note, donor 1 is has been studied on days 3 and 7 while the other donors have been processed on days 2, 5 and 7 for practical reasons (limited accessibility to the FACSsorting facility). In addition, scRNAseq data from intact islets from donor 1 and 2 were combined with our previously generated data from the same donors (and using the same SORTseq protocol) (donor D25 and D28, respectively-GEO: GSE85241 [9]) to increase cell number from intact islets as reference.
A detailed bioinformatics analysis is presented in Appendix A. Briefly, the raw sequencing data were mapped to the human reference transcriptome. The resulting count tables were filtered to remove low-quality cells. Dimensionality reduction of the data was performed using tSNE [20] and cells were clustered using graph-based clustering using Seurat [21] to identify all the pancreatic cell types. Differential expression analysis was performed using the Wilcoxon rank-sum test, and p-values were adjusted with Benjamini-Hochberg correction. Significant genes (Padj < 0.05) were subjected to pathway analysis using Ingenuity Pathway Analysis (www.ingenuity.com, accessed on 5 July 2018).
The transcriptional profiles of dispersed β-cells were compared to scRNAseq data of islets from individuals with T2D [10] using the rank-rank hypergeometric overlap (RRHO) algorithm [22]. We used pseudo-temporal ordering of lineage traced β-cells in the context of canonical αand β-cells. Of note, potential changes due to culture time were taken into account in the analysis by comparing data with control cells sorted at the same time point (data not shown).

Immunofluorescence Microscopy and Flowcytometry
For whole mount immunofluorescence microscopy, formalin-fixed islets were permeabilised using 0.3% triton-X for 1 h and blocked using goat serum for 1 h. Primary and secondary antibodies were sequentially incubated for 24 h and 12 h respectively with occasional shaking. After counterstaining with Hoechst (BD Biosciences, Franklin Lakes, NJ, USA), samples were mounted using DABCO-glycerol on microscopy slides and confocal imaging on the whole islets was conducted using SP8 WLL (Leica, Wetzlar, Germany).
For flowcytometry, islets were dispersed prior to formalin-fixation and permeabilisation with 0.1% triton-X. After blocking with 2% goat serum for 30 min, single cells were sequentially incubated with primary and secondary antibodies for 20-30 min. Cells were analyzed using an LSRFortessa (BD Biosciences, Franklin Lakes, NJ, USA).

RNA/qPCR
Human islets were washed in PBS and total RNA was extracted using the (Micro) RNeasy kit (Qiagen, Hilden, Germany) according to the manufacturer's instructions. Total RNA was reverse transcribed using M-MLV reverse transcriptase (Invitrogen, Waltham, MA, USA) and oligo (dT). Quantitative PCR was performed using a CFX system (Bio-Rad, Hercules, CA, USA). Fold change was calculated using the delta CT method with human β-actin or GAPDH as reference gene. Primers used are listed in Table S8.

Glucose-Stimulated Insulin Secretion
Approximately 50 IEQ per well were placed in a 96-well transwell plate. Glucosestimulated insulin secretion was performed as previously described [8]. Insulin secretion was assessed using a human insulin ELISA kit (Mercodia, Uppsala, Sweden) following the manufacturer's instructions.

Statistical Analysis
All data are expressed as means ± SEM, unless stated otherwise. For analysis of qPCR data, statistical significance of differences between two groups was determined by a paired Student's t test on the delta CT values calculated from the reference gene (β-actin or GAPDH) and the gene in question. A p-value below 0.05 was considered statistically significant.

ER Stress Is Associated with Loss of β-Cell Identity in a Model of Islet Integrity Disruption
Disruption of primary human islet integrity triggers severe phenotypic alterations in β-cells including the conversion of part of insulin-producing β-cells into glucagon-positive α-cells [8]. In order to decipher the dynamic process underlying these events, we combined β-cell-specific lineage tracing [8,17] and scRNAseq on primary human islet cells. Briefly, human islets were dispersed into single cells, and transduced with a lentivirus-based β-cell lineage tracing system and all cells were then reaggregated in microwells. From these aggregates, single cells were sorted from the entire live cell population ('AGG-all') or enriched for lineage-traced β-cells ('AGG-GFP+') at 2, 5 and 7 days post-reaggregation (setup 1 outlined in Figure S1A). Single cells from untransduced intact islets ('ISL-all') were used as control. In total, 4093 cells from three donors were processed for scRNAseq using setup 1 ( Figure 1A,B and Figure S1 and Table S1). In line with our previous findings [8], the detection of GFP+ cells sorted from reaggregated islet cells in the α-cell cluster confirmed the identity change of part of β-cells upon dispersion and reaggregation ( Figure S2A). The lineage-traced AGG-GFP+ cells detected in the α-cell cluster express all classical α-cell genes (e.g., GCG, ARX, CHGB), and at similar levels as canonical α-cells, while most β-cell genes are found to be downregulated ( Figure S2B and Figure 2C). Yet, converted α-cells differ notably from canonical α-cells, displaying higher levels of HLA genes (HLA-A, -B, -C and -E), IFI6 as well as remnant β-cell genes INS and IAPP, and lower levels of α-cell markers ALDH1A, LOXL4 and GC ( Figure S2D,E and Table S2). Altogether these data indicate that converted α-cells are highly similar, but not identical to canonical α-cells.  Table S1. (C) In order to monitor the conversion process in the scRNAseq data, we devised an algorithm to define a cell identity score. Cell identity scores of both canonical (Islet all 'ISL-all'-blue bars) and AGG-GFP+ (Aggregates GFP+ green bars) βand α-cells are represented in a histogram (donor 3). The identity score of canonical β-cells is centred around zero and of canonical α-cells is centred around one. The mean conversion scores of both canonical βand α-cells (µ) and two standard deviations from the mean ( Next, we aimed at monitoring this dynamic process in the scRNAseq data. We developed an algorithm based on the expression of βand α-signature genes to define a cell identity score, thereby allowing a pseudo-temporal ordering of the GFP-labelled cells in relation to canonical αand β-cells ( Figure S2F,G, Table S3). The scores of intact islet cells identify two distinct populations (native β-cells: 0, native α-cells: 1) ( Figures 1C and S3A,B). Interestingly, the β-cell peak observed in the AGG-GFP+ cells is off-centred, indicating alterations in β-cell identity. In addition, pseudo-temporal ordering indicates that 3-10% of AGG-GFP+ cells fell into the category of cells with an 'intermediate' identity between βand α-cells ( Figure  The pathway analysis performed on the intermediate cell-specific genes revealed that 30% of the affected genes are coding for components of the ribosome complex and are upregulated in this cell subpopulation. Furthermore, EIF2 signalling and translationrelated genes, as well as genes characteristic of the unfolded protein response (UPR) are most affected (Table S4 and Figure 1E). Importantly, the quality of these intermediate cells is similar to other AGG-GFP+ cells ( Figure S3F). However, the low cell number in one of the donors hampered the statistical power of the comparison. Nonetheless, the combined analysis of intermediate cells across the three donors confirmed the affected pathways found in donor 3 (data not shown).
We specifically evaluated the UPR across identity score pseudo-temporal ordering. Genes involved in all three arms (IRE1α, PERK and ATF6) are mostly affected mostly around the transition between β-cells and intermediate cells ( Figure 1F). Genes associated with the IRE1α arm are affected earlier, as shown by increased XBP1 expression just before the intermediate zone. This is followed by genes associated with the PERK arm; ATF4, which is preferentially translated upon EIF2α-induced translation arrest, is affected at gene expression level around the intermediate zone. ATF4 target genes (TRIB3, DDIT3 and PPP1R15A) are upregulated in early intermediate cells, whereas the expression of the target gene ATF3 is increased later. ATF6 expression is lost when an α-cell identity is acquired across identity score, concomitant with the increased expression of the negative regulator of ATF6, WSF1. DDIT3, which can also be a target gene in the ATF6 arm, seems to have a more similar pattern to other PERK-related genes. islet cells from T2D individuals. Pathways affected upon dispersion in β-cells across different donors and time points in comparison to the pathways affected by T2D. Pathway categories were assigned to pathways with similar genes responsible for the overrepresentation of the pathway, and such redundant pathways were grouped and displayed together. The columns represent the β-cells and α-cells from the different donors and timepoints post-dispersion or type 2 diabetes vs ND. For representation purposes, the -log10 (p-values) were capped at 3. The timepoint of sorting is annotated above and the donor is annotated below in the figure. (B) The ER stress response upon islet integrity disruption is evidenced by the increased ratio of protein levels of phosphorylated-eIF2α (p-eIF2 α) and eIF2α, shown along with the housekeeping protein GAPDH (donor 9). (C) The dispersion-induced transcriptional profiles significantly overlap with gene expression profiles in βand α-cells from subjects with type 2 diabetes. Rank-Rank Hypergeometric Overlap map between dispersed cell profiles (x-axis) and type 2 diabetes cell profiles (y-axis). The colour codes the -log10 transformed hypergeometric p-value corrected for multiple testing and shows the strength of the overlap. (D-F) Dispersion of isolated human islets leads to increased mRNA expression of the ER stress-related genes XBP1s/XBP1u, ATF3 and CHOP as measured by qPCR. Data are presented as means ± SEM of fold change over intact islets at day 0 (t = d0). n = 3-16 donors; numbers depicted in the bars indicate the number of donors per time point. * p < 0.05, ** p < 0.01, *** p < 0.0005, **** p < 0.0001 vs. intact islets at day 0 (t = d0).
Altogether, these data indicate that the β-cells displaying high Endoplasmic Reticulum (ER) stress lose the expression of typical β-cell markers post-disruption of islet integrity, and that the activation of the three arms of the UPR is associated with the transition from mature β-cells to cells displaying severe alterations in identity.
Importantly, in order to validate that the lineage tracing system is faithfully targeting native β-cells at the start of the experiment, we performed a control experiment allowing examination of day 0 (i.e., before dispersion) GFP + cells: we transduced intact human islets (instead of dispersed islet cells as in the first experimental setup); when GFP fluorescence was visible at 3 to 4 days post-transduction, the islets were dispersed and processed for scRNAseq analysis (experimental setup 2; Figure S4A-E). As expected, most of the lineage-traced cells were β-cells at this time point, as shown at mRNA level (scRNAseq, Figure S4F), and at the protein level, as virtually all analysed ISL-GFP + were C-peptidepositive ( Figure S4H,I). Moreover, any potential effect of the lentiviral transduction itself could be ruled out as none of the classical ER stress-responsive genes were found to be differentially expressed in transduced GFP-labelled β-cells, when compared to unlabelled β-cells from intact islets ( Figure S4G). This is in line with an earlier study where we showed that lentiviral transduction does not alter primary human β-cell function [23].

Stress Signature in the Islet Integrity Disruption Model Displays Similarities with Hallmarks of β-Cell Stress in Type 2 Diabetes
We next examined the cell populations identified as β-cells, and α-cells for comparison. In order to determine the response of these cells to the loss of cell-cell contact and disruption of islet integrity, we compared the gene expression profiles of β-cells from aggregates (AGGall) to the profiles of canonical β-cells (ISL-all). Islet cell dispersion induces stress-response pathways such as the UPR, protein ubiquitination and NRF2-mediated oxidative stress response, and an attenuation of the translation in a more pronounced manner in β-cells than in α-cells (Figure 2A, Tables S5 and S6). At the protein level, we confirmed the increased phosphorylation of the α-subunit of the translation initiation factor eIF2 early after islet cell dispersion ( Figure 2B). Moreover, the UPR is quickly activated after dispersion with the immediate upregulation of XBP1s, followed by ATF3 and later CHOP in the whole islet cell population ( Figure 2D-F).
In addition, we found iron homeostasis to be affected in β-cells with the downregulation of ferritin subunits FTH1, FTL and FTH1P3 gene expression, and the upregulation of the hypoxia inducible transcription factor EPAS1 and of the transferrin receptor TFRC. Other alterations, although more donor-specific, were observed post-islet integrity disruption. Among these, immune response pathways with overlapping genes related to phagosome maturation and antigen presentation were upregulated in β-cells (donor 3). The significant overrepresentation of this pathway cluster is driven by the upregulation of HLA genes (HLA-A, -B and -C), B2M and FOS upon dispersion. Data from this donor also showed an upregulation of interferon-induced genes (e.g., IFI6, IFIT3, MX1) specifically in β-cells. Finally, the upregulation of glycolytic enzymes in β-cells was observed in donor 1. Mitochondrial function was attenuated in α-cells (donor 1), with the downregulation of subunits of the respiratory complexes.
Prolonged ER stress has been associated with β-cell dysfunction, β-cell death, and the development of both type 1 and type 2 diabetes. We therefore investigated whether the altered gene expression profile in our model of β-cell identity loss displays some similarities with the transcriptional profiles of β-cells from individuals with type 2 diabetes. scRNAseq profiles of islets from patients with type 2 diabetes previously published by Segerstolpe et al. were used for comparison [10]. The dispersion-induced transcriptional profiles from our study significantly overlapped with gene expression profiles in βand α-cells from subjects with type 2 diabetes ( Figure 2C). In β-cells, both up-and down-regulated genes largely correlated with each other, whereas, in α-cells, mostly upregulated genes in intact and non-diabetes cells overlapped.
Furthermore, we specifically checked whether the pathways affected upon dispersion were also altered in cells from type 2 diabetic donors ( Figure 2A). Genes involved in translation (EIF2 signalling, mTOR and 'regulation of eIF4 and P70S6K' signalling) and stress response (UPR, protein ubiquitination-related pathways) are differentially regulated in β-cells from individuals with type 2 diabetes. Noteworthy, the ER stress and the NRF2mediated oxidative stress pathways are not detected in the differentially expressed genes in α-cells. Pathways involved in phagosome maturation and iron homeostasis are affected in both cell types. Interestingly, the cell-adhesion-related pathways are also differentially regulated in β-cells from subjects with type 2 diabetes.

ER Stress Leads to β-Cell Dysfunction through Loss of β-Cell Identity
We next asked the question of whether ER stress, caused by islet integrity disruption, is a major driver of the dispersion-induced altered gene expression profile in human β-cells. We subjected primary human islets to treatment with thapsigargin (TG), a typical ER stress inducer, and assessed the effect on key β-cell gene expression and on β-cell function.
Different TG conditions were evaluated. Treatment with 0.1 µM TG for 24 h [24] induced a 3-, 4-, and 6-fold upregulation of the classical ER stress markers XBP1s/XBP1u, ATF3 and CHOP mRNA, respectively ( Figure S5A). In addition, we evaluated the effect of 1 µM TG for 5 h [25] at 0, 24 and 48 h after the end of the treatment ( Figure S5B-D). As expected, the UPR response was transitory, slightly decreasing over time. The ER stress response at 24 h after the end of the 5 h treatment with 1 µM TG was similar to the response seen after the 24 h treatment with 0.1 µM TG (and to the one found in the model of islet integrity disruption ( Figure 2D), so these two TG setups were chosen as a model for ER stress in the rest of the study.
We first assessed the effect of ER stress on key β-cell genes. Gene expression of the maturity marker MAFA was strongly decreased upon ER stress induction, both at mRNA ( Figure 3A) and protein level ( Figure 3C). Similarly, gene expression of the key transcription factor PDX1 and of the regulator of β-cell fate PAX4 was reduced by 30% and 60%, respectively, further indicating the loss of β-cell identity after ER stress induction. In contrast, NKX6.1 and insulin gene expression remained unaffected, hence ruling out that the effects described above may be resulting from β-cell-specific death ( Figure 3A).
Next, the effect of ER stress on endocrine progenitor cell markers was assessed. Gene expression of SOX9, HES1 and C-MYC was increased by 55%, 80% and 65% upon TG treatment, respectively, whereas NEUROG3 expression was unaffected ( Figure 3B). Of note, the expression of the typical duct marker KRT19 was unaffected, ruling out that the upregulated expression of SOX9 and HES1, which mark not only progenitor cells but also adult duct cells, results from an increased exocrine fraction in the islet preparations. Furthermore, TG treatment (1 µM, 5 h) of the human β-cell line EndoC-βH1 confirmed an upregulation of both SOX9 (90%) and HES1 (60%) ( Figure S6). No upregulation of α-cell genes was observed (data not shown). genes was observed (data not shown).
Glucose-stimulated insulin secretion was decreased after 5 h treatment with 1 μM TG in two of three donors ( Figure 3D), suggesting impaired β-cell function due to the occurrence of ER stress.
Overall, these data establish a causal link between ER stress induction and loss of β-cell identity and function, in line with our findings from the islet integrity disruption model.  Glucose-stimulated insulin secretion was decreased after 5 h treatment with 1 µM TG in two of three donors ( Figure 3D), suggesting impaired β-cell function due to the occurrence of ER stress.
Overall, these data establish a causal link between ER stress induction and loss of β-cell identity and function, in line with our findings from the islet integrity disruption model.

Altering Actin Cytoskeleton Affects Human β-Cell Identity and Function
Our findings so far indicate that islet integrity disruption, and thus the loss of cell-cell contact, affect β-cell identity via ER stress. To further investigate the importance of islet integrity for the maintenance of β-cell identity and function, we focused on the role of the actin cytoskeleton in β-cells.
Cell adhesion molecules are connected to the actin cytoskeleton and the microtubule network. Actin cytoskeleton remodelling is necessary for insulin granule trafficking and release. Here we hypothesised that actin cytoskeleton remodelling is crucial for the maintenance of β-cell identity and function. To test this hypothesis, we treated human islets with jasplakinolide (JP), a compound with actin polymerizing-and stabilizing capacities, for 24 h. Firstly, JP treatment reduced β-cell function in islets of two out of three donors ( Figure 4A). Furthermore, the treatment of primary human islets with JP induces ER stress as seen by a 10-fold upregulation of ATF3 gene expression, a 3-fold upregulation of CHOP gene expression and a mild increase in XBP1s/XBP1u ( Figure 4B). In line with our previous observations, ER stress induction is correlated with a clear reduction (50%) in MAFA gene expression ( Figure 4C). PDX1, NKX6.1 and insulin gene expression were not significantly altered, ruling out the toxicity of JP on β-cells. No major change in gene expression of SOX9, HES1 and CK19 was observed (data not shown). In addition, these findings were confirmed in a mature human β-cell line (EndoC-βH3), where increased expression of ER stress-related genes is correlated with a clear reduction in MAFA gene expression ( Figure S7).
Collectively, these data show that altering the cytoskeleton dynamics induces an ER stress response in β-cells, followed by altered β-cell identity and some degree of impairment in β-cell function.

Discussion
The current study shows a β-cell adaptation mechanism to ER stress through altered identity and therefore reduced function.
Using a model of loss of islet integrity, we identified a stress signature highly similar to hallmarks of β-cell stress found in type 2 diabetes. As an adaptive response to mild ER stress, the UPR has been shown to modulate β-cell proliferation in response to increased insulin demand [26]. However, more severe ER stress has also been linked to β-cell dysfunction, β-cell death and to the presence of type 1 [27][28][29] and type 2 diabetes [30][31][32][33].
Here we propose that identity change constitutes a novel mechanism of adaptation for βcells to survive irremediable ER stress, and as an alternative to a more definitive path that would be programmed cell death ( Figure 5).

Discussion
The current study shows a β-cell adaptation mechanism to ER stress through altered identity and therefore reduced function.
Using a model of loss of islet integrity, we identified a stress signature highly similar to hallmarks of β-cell stress found in type 2 diabetes. As an adaptive response to mild ER stress, the UPR has been shown to modulate β-cell proliferation in response to increased insulin demand [26]. However, more severe ER stress has also been linked to β-cell dysfunction, β-cell death and to the presence of type 1 [27][28][29] and type 2 diabetes [30][31][32][33].
Here we propose that identity change constitutes a novel mechanism of adaptation for β-cells to survive irremediable ER stress, and as an alternative to a more definitive path that would be programmed cell death ( Figure 5). Figure 5. Proposed model. We present a model in which ER stress occurs as a result of the loss of cell-cell contact (in this particular experimental setup) and the subsequent remodelling of the actin cytoskeleton. When ER stress is resolved, the β-cell can fully recover. However, high or persistent ER stress can lead to either β-cell death or altered β-cell identity, thereby leading to a reduced functional β-cell mass. Overall, we propose β-cell identity changes as a cell-intrinsic mechanism to survive irremediable cellular stress. This adaptation mechanism may contribute to the development of diabetes.
Despite the fact that single-cell transcriptomics analyses are limited to a few donors only, this type of technology allows generating new hypotheses that can then be tested in other types of assays and using more donors/cells. Using this approach, we showed that activation of ER stress is associated with the loss of key transcription factors. Mouse model studies suggested that β-cell failure can be attributed to β-cell dedifferentiation to a progenitor stage and reprogramming to other endocrine cell types [1]. These observations were supported by evidence of altered β-cell identity found in human islets from donors with type 2 diabetes and type 1 diabetes [34], albeit without the detection of a progenitor state [4,5,35,36]. Our findings that MAFA and PDX1 are more susceptible to a stress-induced decrease in mature β-cells are in line with previous findings both in human cells [7] and in mice, where the failure of the adaptive UPR in a diabetic environment is associated with decreased expression of β-cell maturity markers [37]. Furthermore, upon induction of ER stress by TG, we found some evidence of dedifferentiation, with a reactivation of the endocrine progenitor markers SOX9, HES1 and C-MYC in β-cells. This is in line with earlier reports in other in vitro models of β-cell stress [38,39]. In contrast, both islet integrity disruption and JP treatment did not lead to increased expression of progenitor markers. This apparent discrepancy is likely to reflect the effect of additional mechanisms triggered in the latest two conditions.
In the islet integrity disruption model, as in type 2 diabetes, the stress response appears to be more prominent and more persistent in β-cells than in α-cells. β-cells display specifically altered expression of ER stress, UPR, and NRF2-related pathways, while these pathways are not differentially affected in α-cells. Interestingly, α-cells are thought to cope better with ER stress and related apoptosis because of the expression of specific antiapoptotic proteins [40]. In addition, α-cells have been shown to display higher levels of Figure 5. Proposed model. We present a model in which ER stress occurs as a result of the loss of cell-cell contact (in this particular experimental setup) and the subsequent remodelling of the actin cytoskeleton. When ER stress is resolved, the β-cell can fully recover. However, high or persistent ER stress can lead to either β-cell death or altered β-cell identity, thereby leading to a reduced functional β-cell mass. Overall, we propose β-cell identity changes as a cell-intrinsic mechanism to survive irremediable cellular stress. This adaptation mechanism may contribute to the development of diabetes.
Despite the fact that single-cell transcriptomics analyses are limited to a few donors only, this type of technology allows generating new hypotheses that can then be tested in other types of assays and using more donors/cells. Using this approach, we showed that activation of ER stress is associated with the loss of key transcription factors. Mouse model studies suggested that β-cell failure can be attributed to β-cell dedifferentiation to a progenitor stage and reprogramming to other endocrine cell types [1]. These observations were supported by evidence of altered β-cell identity found in human islets from donors with type 2 diabetes and type 1 diabetes [34], albeit without the detection of a progenitor state [4,5,35,36]. Our findings that MAFA and PDX1 are more susceptible to a stressinduced decrease in mature β-cells are in line with previous findings both in human cells [7] and in mice, where the failure of the adaptive UPR in a diabetic environment is associated with decreased expression of β-cell maturity markers [37]. Furthermore, upon induction of ER stress by TG, we found some evidence of dedifferentiation, with a reactivation of the endocrine progenitor markers SOX9, HES1 and C-MYC in β-cells. This is in line with earlier reports in other in vitro models of β-cell stress [38,39]. In contrast, both islet integrity disruption and JP treatment did not lead to increased expression of progenitor markers. This apparent discrepancy is likely to reflect the effect of additional mechanisms triggered in the latest two conditions.
In the islet integrity disruption model, as in type 2 diabetes, the stress response appears to be more prominent and more persistent in β-cells than in α-cells. β-cells display specifically altered expression of ER stress, UPR, and NRF2-related pathways, while these pathways are not differentially affected in α-cells. Interestingly, α-cells are thought to cope better with ER stress and related apoptosis because of the expression of specific anti-apoptotic proteins [40]. In addition, α-cells have been shown to display higher levels of antioxidant defence mechanisms than β-cells, which make the latter ones more sensitive to the overproduction of reactive oxygen species [41]. These intrinsic differences between α-cells and β-cells may explain a distinct adaptation mechanism to stress from these two cell types.
Our data from the islet integrity disruption model and JP treatment indicate that loss of cell-cell contacts triggers ER stress in β-cells. Importantly, the islet architecture is disrupted in both type 1 and type 2 diabetes, due to the specific destruction of β-cells and the development of amyloid [4], respectively. Cell-cell and cell-ECM interaction, mediated by cadherins and integrins, are critical for the maintenance of β-cell identity and function [42,43]. Recent studies have suggested that cytoskeleton status is crucial for the correct development of β-cells from human pluripotent stem cells [44] and the correct insulin granules movements [45]. Cadherins and integrins participate in adherens junctions where they are connected via catenins and focal adhesions to the actin cytoskeleton and the microtubule network. As the cytoskeleton guides granule movement and exocytosis, dispersion-mediated disruption of the cell-to-cell contacts between the islet cells potentially interferes with β-cell function [46,47]. This may be followed by a disruption of the autoregulatory feedback loop that helps to reinforce the β-cell phenotype [48]. Moreover, independent of their role in the UPR, both PERK and IRE1α has been shown to control actin remodelling and dynamics by serving as a scaffolding protein for the actin crosslinking factor filamin A [49,50]. The actin cytoskeleton is connected to the nuclear envelope and as such extracellular signals can be rapidly transmitted and induce structural changes in the nucleus [51,52]. Altogether, this series of events may explain the altered identity of β-cells. Furthermore, detachment from the extracellular matrix has been associated with elevated ROS levels in epithelial cells. ROS constitute signalling molecules sensing environmental cues [53,54]. Elevated ROS levels have been associated with reduced activity of key transcription factors PDX1 and MAFA in β-cells [55]. As found in our scRNAseq data, oxidative stress and ER stress responses are intertwined, and both processes are known to alter β-cell function. Oxidative stress can induce protein misfolding by disturbing the redox state in the ER, while PERK activation in response to protein misfolding also activates the master regulator NRF2 and the oxidative stress response [56].
In conclusion, we present a mechanism of adaptation of primary human β-cells to cellular stress through the loss of β-cell maturity. This mechanism may play a role in the loss of functional β-cell mass that is associated with the onset and the development of diabetes, and therefore may allow new therapeutic opportunities. First of all, alleviating ER stress may not only help to promote β-cell survival, but also to preserve β-cell identity and therefore function. Additionally, if we assume that cellular plasticity is a reversible process, during a specific time window at least, novel (endogenous) β-cell regeneration therapies could be envisaged for diabetes.

Supplementary Materials:
The following are available online at https://www.mdpi.com/article/ 10.3390/cells10123585/s1, Figure S1. Single-cell RNA sequencing to study the adaptive response of human pancreatic β-cells to disruption of islet integrity. Figure S2. Pseudo-temporal ordering of β-cells shows that loss of identity is associated with upregulation of stress markers (1/2). Figure S3. Pseudo-temporal ordering of β-cells show that loss of identity is associated with the upregulation of stress markers (2/2). Figure S4. Lineage tracing control. Figure S5. Characterisation of the TG model. Figure S6. Treatment of EndoC-βH1 cells with thapsigargin leads to an upregulation of ER stress markers and endocrine progenitor markers. Figure S7. Treatment of mature EndoC-βH3 cells with jasplakinolide leads to an upregulation of ER stress markers and a downregulation of β-cell-specific markers. Table S1. scRNAseq data composition per donor across cell types and conditions. This table summarizes the number of cells per cell types and conditions for each donor. Table S2. Differential expression analysis of converted AGG-GFP+ α-cells vs. canonical α-cells. This table lists the genes identified as significantly differentially expressed in AGG-GFP+ α-cells vs. ISL-all α-cells for each donor (as separate sheets in the Excel file). Table S3. Signature genes for β-cells. This table lists the signature genes identified by differential expression between canonical βand α-cells sorted from intact islets (ISL-all) for each donor (as separate sheets in the Excel file).   For the mapping of the sequenced reads and quantification of the number of transcripts reads 2 from the fastq files were mapped to the human reference transcriptome with BWA [56] and reads 1 were used to extract the cell-barcode sequences and unique molecular identifiers (UMI). UMI-correction was performed to the read and barcode counts resulting in unique transcript count tables. Spike-in reads were removed from the dataset.
The data from the two experimental setups, i.e., donors 1-3 and 4-6, were normalised and filtered separately in R (version 3.3.3). For the first dataset (experimental setup 1, donors 1-3), cells with less than 1500 genes and 5000 transcripts detected in at least 5% of the cells were discarded. The 4093 remaining cells were normalised by downsampling to 5000 transcripts/cell ( Figure S1C). Taking into account the donor variation (D1-D3, as observed in Figure S1D) in downstream analyses, cells from each donor were filtered and normalised separately using the abovementioned methods prior to assessment of differential expression per donor. The second dataset yielded less complex sequencing results as fewer reads/cells were detected in comparison to the first dataset. Therefore, for the second dataset (experimental setup 2, donors 4-6), cells with less than 1500 genes and 2000 transcripts detected in at least 5% of the cells were excluded from normalisation and the remaining 813 cells were used for further analysis ( Figure S4B).

Appendix A.2. Clustering, Visualisation and Cell Type Annotation
Dimensionality reduction of the data was conducted using tSNE based on all genes (Rtsne package, version 0.13) [20]. Of note, the transgene GFP was annotated and detected but removed from the dataset prior to dimensionality reduction and cell clustering. Cells were clustered using graph-based clustering as implemented in Seurat (version 1.4.0.16) [21]. The identified clusters were assigned to cell types based on the expression levels of hor-
To ensure the sorting of individual cells for SORTseq, the gate settings of the FACS were purposely very stringent ( Figure S1B). We then tested a few in silico doublet removal methods but with limited success. Therefore we chose to estimate the doublet rate using simultaneous high expression of two islet cell type marker genes listed above. Using this method, we detected a low (potential) doublet rate in donor 1 to 4 (1.3%, 2.8%, 1.6%, 0.63%, respectively). However, we found that donors 5 and 6, which were sorted in a different facility, showed a much higher doublet rate (4.6% and 7%). Therefore, we chose to exclude all estimated doublets from these datasets for donors 5 and 6.

Appendix A.3. Differential Expression Analysis
Differential expression was assessed using the Wilcoxon rank sum test and p-values were adjusted with Benjamini & Hochberg correction as implemented by the "Stats" R package (version 3.3.3). Fold changes were calculated by the log 2 difference in mean expression. Genes with an adjusted p-value < 0.01 and absolute fold change >2 were considered significantly differentially expressed, unless stated otherwise.
Signature genes were identified for canonical βand α-cell types for each donor separately by differential expression between canonical βand α-cells sorted from intact islets using a |log2FC| > 2.
The effect of transduction on the transcriptional profiles of β-cells was addressed by evaluating the differentially expressed genes, as described above, between untransduced (/unlabelled) and transduced (/labelled) β-cells from intact islet sorted at day 0.

Appendix A.4. Pseudotime Ordering
Pseudotime ordering for each cell of interest was calculated based on the expression of βand αsignature genes (Table S2). The pseudotime for a cell x is defined as: where µ α and µ β are the average expression vectors of signature genes across canonical αand β-cells, respectively. The obtained score of zero represents the average canonical β-cell and one represents the average canonical α-cell respectively ( Figure S2F). The distribution of the pseudotime scores of canonical βand α-cells was used to determine intermediate populations of βto α-cell conversion. More specifically, two standard deviations from the average pseudotime scores of both canonical βand α-cells, sorted from intact islets, were used as cut offs for the pseudotime scores of AGG-GFP + cells to determine mature and intermediate cell identities. Differential expression in intermediate cells was addressed by comparison to non-converted β-cells on the one hand and to converted α-cells on the other hand. Expression of genes of interest was visualised across conversion score using the smooth expression curves as implemented in Monocle ("plot_pseudotime_heatmap", version 2.8.0) [57]. The results obtained using the conversion score were compared to the pseudotemporal ordering algorithm as implemented by the Monocle package (version 2.12.0). The comparison is shown for donor 3 in Figure S3C.

Appendix A.5. Assessing the Effect of Dispersion and Reaggregation on Gene Expression
The transcriptional differences induced by dispersion and reaggregation were addressed by comparing cells from aggregates (combining both GFP-labelled and non-labelled cells for β-cells for each donor separately) to cells from intact islets. In order to restrict the analysis to native α-cells, we excluded GFP-labelled α-cells from the α-cell population. Therefore, even though we cannot rule out that part of the α-cells result from conversion from (non GFP-labelled) β-cells, we estimate that the majority of the considered α-cells are of native origin.

. Pathway Analysis
Lists of differentially expressed genes were subjected to pathway analysis using Ingenuity Pathway Analysis (www.Ingenuity.com, accessed on 5 July 2018). Pathways with Benjamini-Hochberg adjusted p-values below 0.05 were considered significantly overrepresented. For the pathway analysis of the intermediate cells, we manually assigned directionality (up-or down-regulated) to each enriched pathway; direction was assigned when at least 75% of the overrepresented genes in a pathway had the same direction of regulation in the comparison of interest. For the effect of dispersion, pathway categories were assigned to pathways with similar genes responsible for the overrepresentation of the pathway. As such, redundant pathways were grouped and displayed together.