The gene signatures of human alpha cells in types 1 and 2 diabetes indicate disease-specific pathways of alpha cell dysfunction

Glucagon secretion is perturbed in both type 1 and type 2 diabetes (T1D, T2D) the pathophysiological changes at the level of individual pancreatic alpha cells are still largely obscure. Using recently-curated single-cell RNA data from human donors with either T1D or T2D and appropriate controls, we leveraged alpha cell transcriptomic alterations consistent with both common and discrete pathways. Firstly, altered expression of genes associated with alpha cell identity (ARX, MAFB) was common to both diseases. In contrast, increased expression of cytokine-regulated genes and genes involved in glucagon biosynthesis and processing were apparent in T1D, whereas mitochondrial genes associated with reactive oxygen species generation (COX7B, NQO2) were dysregulated in alpha cells from T2D patients. Conversely, T1D alpha cells displayed alterations in genes associated with autoimmune-induced ER stress (ERLEC1, HSP90) whilst those from T2D patients showed changes in glycolytic and citrate cycle genes (LDH, PDHB, PDK4) which were unaffected in T1D. These findings suggest that despite some similarities related to loss-of-function, the alterations of alpha cells present important disease-specific signatures, suggesting that they are secondary to the main pathogenic events characteristic to each disease, namely immune-mediated-or metabolic-mediated-stress in respectively T1D and T2D.


Introduction
Both types 1 (T1D) and 2 diabetes (T2D) are characterised by varying degrees of pancreatic beta cell failure (Eizirik, Pasquali, and Cnop 2020;Marchetti et al. 2020). This is paralleled by dysfunction of alpha cells, which in T1D may contribute to insulin-induced hypoglycaemia and in T2D, at least at the initial phases of the disease, to hyperglycemia (Brissova et al. 2018;Gromada, Chabosseau, and Rutter 2018).
Alpha and beta cells are intermingled in human pancreatic islets (Bosco et al. 2010) and there is a crosstalk between these cells that regulates at least in part their function (Campbell and Newgard 2021). It is thus conceivable that the reduced functional beta cell mass in T1D and T2D impacts alpha cells and contributes to their dysfunction in each disease. Alternatively, it may be that mechanisms inherent to each disease, i.e. a predominance of autoimmunity and consequent islet inflammation in T1D as compared to severe metabolic stress in T2D (Eizirik, Pasquali, and Cnop 2020), directly impair the alpha cells.
Differentiated cells trigger diverse adaptive responses that are determined by the stress to which they are exposed. For instance, beta cells exposed to pro-inflammatory cytokines trigger branches of the unfolded protein response that are different from the ones triggered in response to the metabolic stressor palmitate, and the global gene signatures of islets obtained from patients affected by T1D or T2D are markedly different (Eizirik, Pasquali, and Cnop 2020). The cellular responses to diverse stresses may leave gene expression footprints -particularly in long-lived cells such as human alpha and beta cells -that can be detected by RNA sequencing. Examination of these footprints may enable to differentiate the principal cause(s) of the alpha cell stress present in T1D and T2D. If the leading cause of alpha cell stress is the relative or absolute loss of neighbouring beta cells, and the deficiency of insulin leading to hyperglycemia, we may expect to find similar gene signatures on alpha cells from T1D or T2D patients; on the other hand, if the stress is disease-specific then alpha cells should show different signatures in each case, for instance immune-induced stress in T1D and more metabolic changes T2D.
To test these hypotheses we presently used recently-curated human islet single-cell transcriptomic data from control donors or individuals affected by either T1D or T2D that are publicly available (Kaestner et al. 2019). The results indicate similar patterns, but also major divergences between the gene expression signatures present in alpha cells from T1D or T2D patients, arguing in favour of disease-specific mechanisms leading to alpha cell dysfunction in each case.

Materials and Methods
Analysis of single-cell data and integration with the other datasets Fastq files and the corresponding metadata were downloaded from the database of the Human Pancreas Analysis Program (HPAP) (https://hpap.pmacs.upenn.edu). Reads were aligned using STAR 2.7.3a (Dobin et al. 2013) against the human reference genome GRCh37 (Ensembl 87 annotation) with different parameters according to the technology used for library production.
After the read mapping step, read count tables were analysed with ad-hoc python scripts implementing the toolbox Scanpy (Wolf, Angerer, and Theis 2018). In particular, cell-wise and gene-wise metrics were computed for quality control analyses (QC) to define excluding criteria of low-quality cells, as previously done . The parameters considered for each cell were: i) the number of genes with at least one read mapped (expressed genes); ii) the number of reads mapped to genes (counts); iii) the ratio of reads mapped on mitochondrial genes (mitochondrial fraction). Additionally, genes with detected expression in three or fewer cells were excluded from the analyses. These values were used to exclude cells with high counts and expressed genes (likely representing multiplets), and with high mitochondrial fraction and low expressed genes (indicative of lysed cells). These variables and their covariation were considered separately for each sample, defining separate threshold values to flag cells as "low quality".
After QC, the processed samples were concatenated and the counts normalized to a total of 10,000 for each cell, then log-transformed. The dispersion of each gene with respect to its mean value was computed in the integrated dataset to annotate highly variable genes using the homonymous Scanpy function. The resulting set of genes displaying high variability was used to perform dataset integration with Batch Balanced KNN (BBKNN) (Polański et al. 2020), similarly as done by Park and colleagues . Briefly, BBKNN was first used with the HPAP id as batch key, followed by a Louvain clustering (Lu, Halappanavar, and Kalyanaraman 2015) on the corrected dataset. Then, the obtained clusters were used as biological covariates to perform a ridge regression on the adjusted data, followed by a further BBKNN correction. The transcriptomes of single-cells were visualized using the Uniform Manifold Approximation Projection (UMAP).

Unsupervised cell type annotation
Single-cell clusters were identified with the Louvain modularity algorithm (Lu, Halappanavar, and Kalyanaraman 2015) as implemented in Scanpy (https://github.com/vtraag/louvain-igraph) with resolution=0.5, with further sub-clustering for groups of interest. The genes with cluster-specific expression were found with the "rank_genes_groups" function of Scanpy, curating their association with known cell types using literature information and gene expression markers reported in PanglaoDB (Franzén, Gan, and Björkegren 2019).

Differential expression and gene set enrichment analysis
The differentially expressed genes (DEGs) in T1D and T2D as compared to the corresponding controls were identified using MAST (Finak et al. 2015) with the following mixed effect model: Counts~Diabetes + ngenes + Technology + Race + Sex + BMI + Age, including also a random effect estimated for each individual (Individual). Since the estimation of such an effect may be hampered by a low number of observations, individuals with < 50 cells were excluded from the analysis. Counts is the matrix of raw count data, filtered of genes being expressed in less than 20% of the cells with the filterLowExpressedGenes function, ngenes is a variable encoding the number of expressed genes and Diabetes is a 2-level factor (Disease, ND) indicating the diabetes status of the donor. Genes with corrected p-value (Benjamini-Hochberg, FDR) lower than 0.05 and absolute fold-change (FC) greater or equal than 0.5 were considered as DEGs.
Gene set enrichment analysis (GSEA) was performed for the following datasets: i) Gene Ontology ( i) single-cell technology used to assess transcriptomes, ii) age, iii) BMI and iiii) sex.
Normoglycemic donors showing positivity towards pancreatic auto-antibodies were excluded.
After the selection of control donors on the basis of such criteria, the T1D dataset included 7 diabetic donors and 6 controls, whereas T2D contained 5 affected donors and 5 controls ( Figure   1A). The features of the donors included in this study are reported in Supplementary Table 1.
Analysis of raw sequencing data and exclusion of low-quality cells and genes (see Methods) The number of cells for each cell type in each dataset are reported in Table 1. Alpha cells were then analysed separately in each dataset to find transcriptomic alterations potentially linked to diabetes.

The transcriptional signatures of alpha cells in T1D
The T1D dataset comprised 2,362 alpha cells, collectively expressing 6,265 genes (after filtering). A more detailed break-down of the number of cells at the level of individual, single-cell technology used and diabetes status is reported in Supplementary Table 3. After removing cells from individuals with less than 50 transcriptomes, the resulting dataset, including 2,225 alpha cells, was analysed with MAST to identify genes differentially expressed in T1D cells.
A total of 346 DEGs were identified, of which 193 were up-regulated and 153 were down-regulated versus cells from the control, normoglycemic group (Supplementary Table 4).
A pattern observed among the DEGs was the over-expression of genes involved in Reactive Oxygen Species (ROS) response (PRNP, NDUFA6, NDUFB4, GLRX and TXN), and protein folding stability via chaperone activity (HSP9, HSP90AA1 and HSP90AB1). We also observed overexpression of IL-8 and HLA-A, genes downstream of the transcription factors NF-κB and STAT1/STAT2 that are regulated by the pro-inflammatory cytokines IL-1β and types I and II interferons (IFNs) and participate in the immune system-islet cell dialogue present in T1D (Eizirik, Pasquali, and Cnop 2020). Importantly, these genes were not overexpressed in alpha cells from patients affected by T2D (see below). The over-expression of DDIT3 (also known as CHOP), a key mediator of ER stress induced-beta cell death (Eizirik, Cardozo, and Cnop 2008), fits with this scenario. Another trend of pathophysiological relevance was the significant down-regulation of genes implicated in endocrine function, namely PCSK2, CHGA, SRP14 and PAK3. These genes contribute to proglucagon peptide maturation and eventual exocytosis, and their decreased expression could contribute to reduced glucagon secretion from alpha cells in human T1D (Gerich et al. 1973). Thus, in PCSK2 gene-null mice lacking the prohormone convertase in alpha cells, proglucagon maturation to the mature hormone is blocked (Furuta et al. 2001). Of interest, there was also down-regulation (-2.5 fold change) of PCSK1N, a specific inhibitor of PCSK1.

The transcriptional signatures of alpha cells in T2D
After filtering out low quality cells, in the T2D dataset there were 1,757 alpha cells expressing a total of 6,872 genes. Following the approach used for the T1D dataset, we detailed the number of alpha cells at the levels of individual, single-cell technology and diabetes status (Supplementary  (Yamaguchi and Wang 2004;Ohoka et al. 2005;Eizirik, Cardozo, and Cnop 2008). TXNIP and ATF3 are instead modulated by glucose concentration and play a role in both apoptosis and hormone secretion. TXNIP controls a pathway (miR-204/MafA/insulin) that reduces insulin expression (Xu et al. 2013), whereas ATF3 increases the expression of proglucagon, as well as regulating expression of genes related to apoptosis (such as GADD45, BNIP3 and NOXA).
The over-expression of such genes, and others related to ROS (SBNO2, EGLN2, MBP) and the unfolded protein response (UPR; EEF2, ATF4, HERPUD1) suggests that alpha cells in T2D are subject to oxidative stress, a hallmark of toxicity induced by elevate glucose levels (Kawahito, Kitahata, and Oshita 2009;Gromada, Chabosseau, and Rutter 2018). The down-regulation of genes involved in pyruvate metabolism (LDHA, PDHB, PDK4) and oxidative phosphorylation (COX7B, NQO2, SUCLA2, UQCR10, SLC25A4) is also consistent with oxidative stress pathways, in that these may affect ROS production by mitochondria.
Of functional relevance to hormone biosynthesis, the genes encoding prohormone convertases

Contrasting the alpha cell signatures of T1D and T2D
By analysing the changes in gene expression in alpha cells separately in T1D and T2D the corresponding transcriptomic signatures were identified and described, revealing a number of similarities, such as upregulated DEGs involved in the stress response or the downregulation of genes relevant for the secretory function. To quantify more precisely the extent at which the two series overlap and to underscore their differences, a systematic comparison of the obtained results was undertaken next.
In T1D and T2D there were a total of 770 DEGs (Supplementary Table 8 for GSEA is reported in Supplementary Table 9. Although the signatures retrieved for T1D and T2D indicate a partial overlap of stress response genes, we found some remarkable differences between the two series. In T1D, genes linked with ER stress and UPR are either significantly up-regulated (ERLEC1, HSP90) or with a significant change that is below the fold-change threshold we used (ATF6, ER mannosidase I, TRAM1), an effect that was not observed in T2D (Supplementary Figure 3).
Another major difference between T1D and T2D involved metabolism i.e. glycolysis, citrate cycle and oxidative phosphorylation (Supplementary Figure 4). Notably, in T1D we observed little changes in expression of genes involved glycolysis and citrate cycle, but increased levels of those contributing to oxidative phosphorylation, whilst these pathways were repressed in T2D: most of the genes involved in glycolysis were down-regulated, four with high significance (PKM, ENO1, ALDH2, TPI1) and three with high significance and fold-change (LDHA, ALDH9A1 and PDHB); strikingly, most of the citrate cycle genes were down-regulated in alpha cells from T2D donors, with six genes displaying a significant change; among the complexes involved in oxidative phosphorylation, genes contributing in complex III were the most negatively affected, with almost all genes being under-expressed in diseased alpha cells.

Discussion
Analysis of single-cell data sets provides a unique opportunity to better understand the changes at the level of the alpha cell which drive dysregulated glucagon secretion in T1D and T2D. In the present work we used publicly-available single-cell transcriptomic data from human islets to test whether the dysfunction of alpha cells in T1D and T2D had common bases or if the transcriptional signatures are more disease-specific. This revealed both shared features (i.e. down-regulation of genes related to alpha cell identity and function, and the up-regulation of stress response mechanisms and inflammation signatures) and, more importantly, disease-specific stress pathways. In T1D, several of the alpha cell alterations were traced back to ER stress, a process associated with chronic inflammation and autoimmune diseases. In T2D, up-regulation of ROS defense mechanisms (SBNO2, EGLN2, MBP) was accompanied by modification of the central metabolism, with the repression of genes involved in glycolysis (LDHA, PDHB, PDK4), citrate cycle and mitochondrial respiration (COX7B, NQO2, SUCLA2, UQCR10, SLC25A4). We note that some of these gene expression changes in T2D alpha cells are akin to findings reported in a recent study using a separate data set (Dai et al. 2022) which highlighted alteration of mitochondrial respiratory complex genes in T2D and inhibitory role of H 2 O 2 in glucagon secretion. However, no data were provided in the latter report on alpha cells from patients with T1D, nor was any comparison made of changes in the two disease types.
We speculate that in T2D down-regulation of glucose catabolism is adaptive to ROS stress, not only by leading to a decrease of radicals from oxidative phosphorylation but also by redirecting glucose flux to Pentose Phosphate Pathway, increasing NADPH to improve ROS scavenging (Mullarky and Cantley 2015). However, the relationship here is complex since lowered LDHA (favouring pyruvate flux into mitochondria) and PDK4 (favouring PDH dephosphorylation and activation) would seem likely to oppose the impact of lowered PDHB activate (decreasing PDH-E1 levels and conversion of pyruvate into acetyl-CoA) (Supplementary Figure 3).
Measurements of the corresponding gene products at the protein level will be important to substantiate these mRNA-based findings (Supplementary Figure 3). Of note, glucose oxidation is lowered in islets from T2D subjects (Del Guerra et al. 2005) though the relative contribution of changes in beta and alpha cells is not established. We also report that, among the respiratory chain genes, the most negatively down-regulated were those involved in mitochondrial complex III, a major producer of ROS with implications into cellular transduction (Bleier and Dröse 2013). Additionally, and complementing the above mechanism, a decreased ability of mitochondria to synthesise ATP in response to elevated glucose concentrations (Ravier and Rutter 2005) may contribute to the failure of the alpha cell in T2D to efficiently shut down glucagon secretion as glucose concentrations rise, consistent with other recent findings (Knudsen et al. 2019).
In T1D we observed signatures of ER stress, with most of the genes involved in ER protein processing being up-regulated (Supplementary Figure 4). Considering that in T1D we did not observe repression of the central metabolism (Supplementary Figure 3), the partial overlap of ROS response between T1D and T2D could be due to the fact that ER stress increases the ROS cellular levels (Zeeshan et al. 2016). Related to ROS production, oxidative phosphorylation is enhanced in T1D, a striking difference with respect to T2D alpha cells in which the pathway is repressed. Finally, alpha cells in T1D display signatures of inflammation, including markers of cytokine exposure, with multiple related pathways among the most significantly enriched.
Besides a stress response, we reported for T1D the down-regulation of a number of genes involved with hormone secretion, with a likely relevance in the pathophysiological context. The apparent repression of PCSK2 transcription is consistent with reduced glucagon production, as well as the decreased expression of PCSK1N, a repressor of PCSK1. As this gene is involved with glucagon-like peptide 1 (GLP-1) production from glucagon (Rouillé et al. 1997), the pattern we observed would be consistent with increased endogenous production of GLP-1 and lower glucagon maturation. Also, the down-regulation of CHGA, encoding chromogranin A, implies impaired hormone secretion with a role in secretory granule biogenesis (Kim et al. 2001).
In conclusion, despite their similarities, the alterations of alpha cells present important disease-specific signatures, suggesting that they are, at least in part, secondary to the main pathogenic events characteristic to each disease, namely immune-mediated-or metabolic-mediated-stress, respectively, in T1D and T2D. Nevertheless, we note that our findings are derived only from single-cell transcriptomics, and are therefore subject to caveats over sensitivity and other limitations associated with this approach (Mawla and Huising 2019).
Future studies will be required to dissect the molecular mechanisms involved in the observed changes, their relevance for disease pathogenesis and the potential for targeted therapeutic interventions.
from the Union's Horizon 2020 research and innovation programme, "EFPIA", "JDRF" and

Figure 2
Cell type annotation of T1D and T2D datasets The single-cells in the UMAP plots of T1D (left) and T2D (right) datasets have been clustered and annotated to a pancreatic cell type (top) based on the expression of known marker genes (bottom).

Figure 3
Comparison of differentially expressed genes in T1D and T2D alpha cells

Figure 4
Comparison of most significantly enriched pathways in T1D and T2D Gene-set enrichment analysis was performed for T1D (left) and T2D (right) datasets. The barplots report the normalized enrichment score (bar length) and the significance (bar color) of the most significant 15 positively and negatively enriched pathways from KEGG (top) and Reactome (bottom).
Supplementary Figure 1

UMAP representation of T1D with confounders
The single-cells from the T1D dataset are reported with cells color-labeled according to the single-cell technology used for sample preparation (left) and the HPAP ID of the donor (right).
Supplementary Figure 2

UMAP representation of T2D with confounders
The single-cells from the T2D dataset are reported with cells color-labeled according to the single-cell technology used for sample preparation (left) and the HPAP ID of the donor (right).

Expression of genes involved in Glycolysis in T1D and T2D datasets
The KEGG map of the Glycolysis pathway (hsa00010) is reported with the enzymes colored according to the fold-change of the corresponding genes in T1D (left) and T2D (right) datasets.

Expression of genes involved in Unfolded Protein Response in T1D and T2D datasets
The KEGG map of the Unfolded Protein Response pathway (hsa04141) is reported with the enzymes colored according to the fold-change of the corresponding genes in T1D (left) and T2D (right) datasets.
Supplementary Table 1 Features of donors in T1D and T2D datasets The sheet "Donor Features" reports the information of each donor included in the study. This Supplementary Table 2 Thresholds used to filter low-quality single-cell transcriptomes The table reports the filters used for each single-cell sample to flag low-quality cells.
Supplementary Table 3 Alpha cell count table The Supplementary Table 4 Alpha cell genes differentially expressed in T1D The table reports the results of the differential expression analysis of alpha cells in the T1D dataset. The reported DEGs are sorted according to their differential expression adjusted significance (FDR).
MAST enrichment analysis, with the effect and significance of each component of the statistical model; the terms are sorted according to their enrichment direction (1=positively enriched in T2D, -1=negatively enriched in T2D) and significance ("combined_adj"). The results are reported for each dataset in a separate sheet, whereas the aggregated results for all the datasets are included in the "All" sheet.
Supplementary Table 8 Differentially expressed genes shared between T1D and T2D Partitioning DEGs in T1D and T2D according to their direction (up-or down-regulated) allowed the comparison between up-regulated genes in T1D and in T2D, down-regulated genes in T1D and T2D, up-regulated in T1D and down-regulated in T2D, and down-regulated in T1D and up-regulated in T2D. The results of such comparisons are reported in each sheet, with a gene classified as "Shared" if it is differentially expressed in both conditions, "unique to T1D" or "unique to T2D" if it is differentially expressed only in one condition.
Supplementary Table 9 Comparison of GSEA enriched terms between T1D and T2D Partitioning enriched functional terms in T1D and T2D according to their direction (up-or down-regulated) allowed the comparison between up-regulated genes in T1D and in T2D, down-regulated genes in T1D and T2D, up-regulated in T1D and down-regulated in T2D, and down-regulated in T1D and up-regulated in T2D. The results of such comparisons are reported in each sheet, with a gene classified as "Shared" if it is differentially expressed in both conditions, "unique to T1D" or "unique to T2D" if it is differentially expressed only in one condition. G C G I N S I A P P P P Y S S T P T P R C P R S S 1 K R T 1 9 S P A R C