Transcriptional profiles of human islet and exocrine endothelial cells in subjects with or without impaired glucose metabolism

In experimental studies, pancreatic islet microvasculature is essential for islet endocrine function and mass, and islet vascular morphology is altered in diabetic subjects. Even so, almost no information is available concerning human islet microvascular endothelial cell (MVEC) physiology and gene expression. In this study, islets and exocrine pancreatic tissue were acquired from organ donors with normoglycemia or impaired glucose metabolism (IGM) immediately after islet isolation. Following single-cell dissociation, primary islet- and exocrine MVECs were obtained through fluorescence-activated cell sorting (FACS) and transcriptional profiles were generated using AmpliSeq. Multiple gene sets involved in general vascular development and extracellular matrix remodeling were enriched in islet MVEC. In exocrine MVEC samples, multiple enriched gene sets that relate to biosynthesis and biomolecule catabolism were found. No statistically significant enrichment was found in gene sets related to autophagy or endoplasmic reticulum (ER) stress. Although ample differences were found between islet- and exocrine tissue endothelial cells, no differences could be observed between normoglycemic donors and donors with IGM at gene or gene set level. Our data is consistent with active angiogenesis and vascular remodeling in human islets and support the notion of ongoing endocrine pancreas tissue repair and regeneration even in the adult human.

In experimental studies, pancreatic islet microvasculature is essential for islet endocrine function and mass, and islet vascular morphology is altered in diabetic subjects. Even so, almost no information is available concerning human islet microvascular endothelial cell (MVEC) physiology and gene expression. In this study, islets and exocrine pancreatic tissue were acquired from organ donors with normoglycemia or impaired glucose metabolism (IGM) immediately after islet isolation. Following single-cell dissociation, primary islet-and exocrine MVECs were obtained through fluorescenceactivated cell sorting (FACS) and transcriptional profiles were generated using AmpliSeq. Multiple gene sets involved in general vascular development and extracellular matrix remodeling were enriched in islet MVEC. In exocrine MVEC samples, multiple enriched gene sets that relate to biosynthesis and biomolecule catabolism were found. No statistically significant enrichment was found in gene sets related to autophagy or endoplasmic reticulum (ER) stress. Although ample differences were found between islet-and exocrine tissue endothelial cells, no differences could be observed between normoglycemic donors and donors with IGM at gene or gene set level. Our data is consistent with active angiogenesis and vascular remodeling in human islets and support the notion of ongoing endocrine pancreas tissue repair and regeneration even in the adult human.
Endothelial cells are remarkably heterogeneous in different vascular beds 1 . Using > 32,000 endothelial cell single cell RNA-seq sets from 11 different murine tissues, Kalucka et al. constructed a murine endothelial cell transcriptome atlas containing 78 endothelial cell subclusters. They found that endothelial cells in mice differ between and within tissues and kinds of vascular beds, and that endothelial cell subsets in multiple tissues express markers associated with angiogenesis 2 .
Microenvironmental cues originating in the neighboring tissue are thought to be a co-regulator of endothelial cell heterogeneity, contributing to microvascular adaptations to local needs 3 . In the pancreas, β-cells and islet endothelial cells signal bidirectionally through soluble factors, extracellular matrix proteins, and cell-bound molecules 4 . β-cells secrete VEGF-A, which is seemingly an important factor for vessel size, density, and fenestration. Pancreatic endothelial cells have shown importance in β-cell differentiation during organogenesis 5 and in stimulating the proliferation of β-cells 6,7 . Endothelial-to-β cell signaling facilitates glucose stimulated insulin secretion 8,9 and increases islet insulin content 9 . Disrupted morphology of the islet microvasculature is correlated with reduced β-cell mass and diabetes in various rodent models 10 .
The endocrine pancreas is highly vascularized, with rat islets receiving roughly 10% of the pancreatic perfusion under basal conditions 11 while only making up 1% of the pancreatic volume. Rat islet endothelium contains ten times as many fenestrae as the exocrine endothelium, with a sharp margin at the islet-exocrine interface 12 . In human organ donors without diabetes, islet microvascular vessels have a larger diameter and lower density than the exocrine microvasculature 13 . An increase in islet capillary density has been found in patients with type 2 diabetes compared with non-diabetic subjects 14 . Insulin has been shown to have general vasodilating effects, acting through local production of nitrous oxide in vascular endothelium 15 among other mechanisms. A recent study 16 found a general proliferation in all pancreatic cell types in organ donors with an increased duration of stay in the intensive care unit, suggesting that a general proliferation can be triggered in the human pancreas through hitherto unknown processes.
Thus, data from human and rodent studies suggest that islet endothelial cells have important roles in islet physiology. Based on the earlier findings of differences between microvascular beds, we expected to find differences between microvascular endothelial cells in the neighboring endocrine and exocrine compartments of the pancreas also at the gene expression level. Identifying such differences could grant insights into the endothelial characteristics necessary for maintaining normal function in the human pancreas. In this study we thus compared gene expression profiles of primary microvascular endothelial cells (MVECs) from human islets of Langerhans and exocrine tissue.

Results
Endothelial cells are highly enriched in sorted microvascular endothelial cell samples. Relative expression levels of tissue-specific genes targeted by quantitative real-time PCR (RT-qPCR) are summarized in Fig. 1A. Endocrine and exocrine tissue genes were mainly detected in bulk samples from the corresponding tissue, with some signal also detected in the corresponding sorted samples. Endothelial genes PECAM and VWF were almost exclusively detected in sorted MVEC samples.
Through transcriptomics-based deconvolution of cellular constituents using a previously published data set 17 , a median of 92% (range 70-97%) of the cells in our sorted samples were identified as endothelial cells (Fig. 1C). FACS and transcriptomics-based deconvolution correlated well (data not shown).
The tissue provenance of endothelial cells, but not glycemic status, is reflected in global gene expression. Principal component analysis (PCA) revealed a trend towards islet-and exocrine endothelial cell samples forming two separate clusters ( Fig. 2A). There was no obvious separation of normoglycemic donors and donors with IGM on PC1 or PC2. There was a slight tendency for samples from the same donor to associate (not shown).
For glycemic status, no differentially expressed genes or enriched gene sets were identified. In a heatmap constructed using the 1000 genes with the most variance, islet and exocrine samples tended to create separate, somewhat overlapping clusters (Fig. 2B).

Islet-and exocrine MVECs differ in expression of a large number of genes. A total of 13,451
genes passed the filtering (filterByExpr) step, and 1249 genes were found to be differentially expressed between islet and exocrine samples, using an FDR < 0.05 assuming a log fold change of ± log2 (1.2) as significance threshold. A total of 676 genes were enriched in islet samples, while 573 genes were enriched in exocrine samples (Fig. 3A). The full results of the differential gene expression analysis can be found in GEO (GSE157546).
Pre-ranked gene set enrichment analysis (GSEA) was performed using 2937 "gene ontology: biological processes" (GO:BP) gene sets containing 15-500 genes each. 347 gene sets were enriched in islet MVEC samples, and 138 gene sets were enriched in exocrine MVEC samples. The full GSEA results for all included gene sets are available online at GEO (GSE157546). The most enriched gene sets in islet samples contain several GO terms related to vascular development (Fig. 3B). The most enriched gene sets enriched in exocrine samples contain several GO terms related to ribosomal assembly and protein translation (Fig. 3B). Leading edge analysis was performed to identify gene sets with an overlap in the genes driving gene set enrichment. The overlap between enriched gene sets are visualized in Fig. 3C.
Comparison to murine tissue-specific and vascular-process specific endothelial cell markers. Islet endothelial cell samples were enriched for kidney, small intestine, lung, colon, soleus and spleen (FDR < 0.25) marker sets (Fig. 4). The most enriched tissue set was "kidney" followed by "small intestine" and "colon". No tissue set was significantly enriched in exocrine samples.
Islet endothelial cell samples were also enriched for the "angiogenic" marker set while the exocrine endothelial cell samples were enriched for the "proliferation" marker set (FDR < 0.25).

Discussion
Transcriptional profiles of highly enriched human pancreatic MVECs from donors with normoglycemia and donors with IGM are presented. The most notable findings pertain to the identification of transcriptional differences between islet and exocrine MVECs. As determined by GSEA, multiple gene sets involved in general vascular development and extracellular matrix remodeling were enriched in islet endothelial cells.
Data from RT-qPCR, flow cytometry, and MuSic demonstrate that endothelial cells constituted an absolute majority of cells in the sorted samples. Tissue was processed immediately after islet isolation in order to minimize the impact of isolation-and culture associated changes to RNA expression. In conjunction with average logCPM being a factor in the significance testing, these quality control measures argue for the majority of the differences in gene expression reflecting true phenotypical differences between islet-and exocrine endothelial cells. It is important to note that there could exist subtypes of endothelial cells within the islets, and that the subtypes may have different transcriptional profiles. Notably, in the context of angiogenesis, endothelial tip cells would have a different transcriptional profile compared with other endothelial cells 2,18 . While single-cell RNA sequencing is a powerful tool, it is limited by a low fidelity of detection. Many genes found to be expressed  www.nature.com/scientificreports/ by bulk RNA-sequencing cannot be reliably detected in any given cell using single-cell RNA sequencing 19 . To elucidate general differences between endothelial cells in islets and exocrine tissue we therefore instead opted for bulk RNA-seq. Due to the paired nature of our islet-and exocrine MVEC samples, we performed pre-ranked GSEA for comparisons between tissues. Many of the enriched gene sets in islet EC compared with exocrine EC were related to angiogenesis and vascular development, endothelial cell development and migration, and cardiovascular system regulation. Several gene sets related to nerve development and differentiation were also found to be upregulatedthis likely represents genes with dual functions. Endothelial tip cells and axonal growth cones share multiple features, with a large overlap in their signaling molecules. The gene sets related to axonal guidance could thus   www.nature.com/scientificreports/ be enriched due to an upregulation of vascular guidance molecules. One of the most significantly upregulated genes in the islet samples was UNC5B, which is almost exclusively expressed in vasculature 20,21 , lending credence to the interpretation of genes with dual functions driving upregulation in nerve cell-related genes. Gene sets related to maintaining and modifying the extracellular matrix were upregulated in islet samples, supporting the notion of ongoing remodulation of the islet vasculature. The findings at the level of gene sets are reflected by an upregulation of several genes that have previously been described to be involved in angiogenesis and morphological processes. PLVAP is known to be important for endothelial fenestration and has been suggested to be upregulated by VEGF-A 22 . Although mature vascular beds are thought to be largely quiescent and independent of VEGF-A signaling 23,24 , VEGF-A has been shown to influence vascular survival and fenestration in various endocrine organs in the adult mouse 23 . Both PLVAP and VEGF-A signaling appears to be upregulated in islets. This is consistent with a greater fenestration in human islet capillaries compared with exocrine capillaries. In mice, β-cells produce VEGF-A which is an important modulator of the islet vasculature and contributes to normoglycemia 6,14 . Notably, while inactivation of VEGF-A at any time point leads to reduced islet vascularization, inactivation during pancreatic development leads to reductions in β-cell mass and function 24 . Both VEGF-A and VEGF-C, as well as FLT1 and KLK, the genes encoding the VEGF receptors, were upregulated in islets in our study. This is consistent with a similar mechanism in humans as that previously described in mice. VEGF-A is an important mediator in angiogenesis, and the VEGF-A receptor FLT1 is typically upregulated in sprouting angiogenesis 21 .
It has previously been hypothesized that islets have a natural progression where they grow until the metabolic demand exceeds the vascular supply and a vascular catastrophe occurs 25 . The results presented herein are consistent with this hypothesis. The increased levels of vascular guidance molecules, such as VEGF-A and other markers of angiogenesis, is consistent with a higher level of angiogenesis in islets compared with exocrine tissue. If islets exist in a state of perpetual growth, a high degree of angiogenesis and vascular remodeling would be necessary to avoid vascular catastrophe and islet death. It is important to remember, however, that genes are often pleotropic and that genes known to be important for angiogenesis may have a different function in the pancreatic islet endothelial cells. There is some support in the literature for the feasibility of ongoing tissue repair in the islets of the adult human pancreas. In a recent report, Smeets et al. found an increase in M2 macrophages, microvascular density, and islet cell proliferation in the pancreases of human organ donors with prolonged intensive care unit stay 16 . This expands further on an earlier report of increased β-cell proliferation in a subset of organ donors 26 . Thus, there is previous immunohistochemical evidence of ongoing tissue repair and cellular proliferation in human islets, and our data provide further evidence for this at the transcriptome level.
Although we found ample differences between islet-and exocrine tissue endothelial cells, we did not identify any differences between normoglycemic donors and donors with IGM at gene or gene set level. This is a surprising result, considering earlier findings of morphological differences between normoglycemic and diabetic subjects in both human and rodent studies. Previous work using RNA-Seq has identified differences in the expression of several genes in cubital vein endothelial cells derived from non-diabetic patients and patients with type 2 diabetes 27 . However, the study by Beckman et al. involved macrovascular endothelial cells, which are not directly comparable to MVEC. One interpretation that reconciles our data with previous studies would be that the reported differences seen between normoglycemia and dysglycemia reflect vasodilatory or vasoconstrictory www.nature.com/scientificreports/ processes, which occurs in some rodent models 24,[28][29][30] . Thus, the reported differences in diabetes could reflect alterations at another level of the microvascular system than the endothelial cells. There was, however, a heterogeneity in clinical characteristics among the subjects in our study, such as differences in HbA1c, treatment regime, and presumably type 2 diabetes duration. It is possible that this heterogeneity could mask subtle or conditional differences between microvascular endothelial cells from normoglycemic donors and donors with IGM. Acinar cells secrete digestive enzymes vital for proper gastrointestinal function. They have high levels of protein synthesis and, during fasting, a correspondingly high protein catabolism 31 . In our exocrine MVEC samples, we found multiple enriched gene sets that relate to biosynthesis and biomolecule catabolism. These findings could reflect an increased turnover of proteins also in these cells. No statistically significant enrichment was found in gene sets related to autophagy or endoplasmic reticulum (ER) stress, albeit there was a tendency towards enrichment in sets related to ER stress.
A murine endothelial cell transcriptome atlas, containing 78 endothelial cell subclusters using single cell RNA-seq sets from 11 different murine tissues 2 , identified marker genes for angiogenesis, interferon signaling and proliferation, as well as tissue specific markers for capillary endothelial cells. Corresponding angiogenesis markers were found to be enriched in our islet samples while corresponding proliferation markers were enriched in our exocrine samples. Our findings may indicate some similarity between endothelial cells from human pancreatic islets and murine kidney, small intestine, and colon.
Our study of the human endothelial cell transcriptome has some limitations. Primarily, the number of organ donors is relatively small. This is inherent to the complicated procedure of obtaining human islets within less than one hour after isolation. However, this is of utmost importance since a period of islet culture could potentially influence the endothelial cell transcriptome. Further, human organ donors with type 2 diabetes are rare and not commonly available for islet research, adding to this limitation in the comparisons between normoglycemia and IGM. Organ donors are subjected to intensive care and critical illness before islet isolation, which may have different effects across tissues and individual donors. It is also possible that the putative differences we have found are affected by preferential acquisition of different endothelial cell subsets in different tissues. For example, younger endothelial cells may be more likely to survive the isolation and sorting process in islet samples. Losses are common in single-cell dissociation and sorting, with a recent study of the human pancreas recovering 4% of the estimated total number of β-cells 32 . Thus, our and other bulk RNA-seq studies require an assumption that cells are not lost in a biased fashion. Another important limitation of our study is the lack of confirmatory results on the protein level supporting the reported transcriptomic findings. However, based on the findings reported in the present communication we will be entitled to obtain a larger number of human islets, allowing proper characterization of the endothelial cell proteome.
To summarize, we present transcriptional profiles from an important cellular component of the human pancreatic islet. Our data demonstrates that islet capillaries differ markedly from their exocrine counterparts, with an enrichment in gene sets related to vascular development and extracellular matrix remodeling. The most parsimonious interpretation would be an ongoing remodulation of the islet vasculature even in the adult human pancreas.

Methods
Organ donors and Islet isolation. Transplant-grade human pancreata were collected from brain dead multi organ donors through the Nordic Network for Islet Transplantation. Tissue from 8 donors with normoglycemia (defined as HbA1c < 42 mmol/mol and no diagnosis of diabetes) and 5 donors with impaired glucose tolerance (IGM-defined as diagnosis of type 2 diabetes and treatment for type 2 diabetes, or HbA1c > 42 mmol/ mol) was included in the study. Available islet and organ donor data is presented in Table 1. Islets and exocrine tissue were isolated at the islet isolation facility (Rudbecklaboratoriet, Uppsala) using a previously described method 33 . Single cell dissociation was initiated immediately after completion of islet isolation.
Ethics statement. Consent for organ donation (for clinical transplantation and for use in research) was obtained via online database (https ://www.socia lstyr elsen .se/en/apply -and-regis ter/join-the-swedi sh-natio naldonor -regis ter/) or verbally from the deceased's next of kin by the attending physician and documented in the medical records of the deceased in accordance with Swedish law and as approved by the Regional Ethics Committee in Uppsala (Dnr 2009/371/2). The present study protocol was approved by the Regional Ethics Committee in Uppsala (Dnr 2015/444). All experiments were carried out in accordance with applicable regulations and ethical guidelines. No organs were procured from prisoners. Table 2. From each donated pancreas, 10,000 IEQ were collected from an endocrine tissue isolate (endocrine content of 78-100% based on automated microscopic evaluation after Dithizone staining) and a similar amount of tissue was collected from an exocrine tissue isolate from the same donor (endocrine content of 0-4%). Single cell suspensions were generated through incubation in Accutase for 20 min at room temperature.

Single-Cell dissociation and fluorescence-activated cell sorting of endothelial cells. Reagent details are presented in
The disassociated cells were then stained for CD31-BB515, CD34-BV421, and CD45-APC-H7 for 30 min at 4 °C, washed and resuspended in a sorting buffer (Tyto running buffer). Small aliquots of the antibody-stained cells were collected for later use in flow cytometric and RT-qPCR analysis of bulk samples. Cell suspensions were passed through a cell strainer before sorting. Cells double positive for CD31 and CD34 and negative for CD45 (CD31 + CD34 + CD45 − ) were considered endothelial and sorted at 8 °C using a microchip-based fluorescence activated cell sorter (MacsQuant Tyto, Miltenyi) which allows for sterile, low pressure sorting of cells. Purple laser (CD34-BV421) was set as trigger channel and blue laser (CD31-BB515) as speed channel. Aliquots were RNA extraction and RT-qPCR. RNA was extracted from the cells lysed in buffer RLT+ using the RNeasy Plus Mini kit according to the manufacturer's instructions.
RT-qPCR was performed as a quality control measure of successful cell-sorting, targeting endocrine (INS, GCG ), exocrine (AMY2A) and endothelial (PECAM, VWF) genes. The cDNA used in RT-qPCR was synthesized using the SuperScript VILO cDNA Synthesis Kit. Reactions were run for 40 cycles on a StepOne Plus instrument with settings following the manufacturer's recommendations. Reactions were run in duplicate and specificity was confirmed with melt curve analysis of the amplified products. Relative expression levels were calculated using the 2 ΔCT method, with GAPDH and 18S as reference genes. In one exocrine MVEC sample, relative expression was calculated using 18S only due to multiple temperature peaks in the GAPDH duplicates. Melt curves were not available for one donor due to hardware failure. In another donor, cDNA was not available from bulk exocrine Table 1. Organ donor data. CIT, Cold Ischemia Time. IA2A, insulinoma-associated protein 2 autoantibody. GAD glutamic acid decarboxylase autoantibody, ICH intracerebral hemorrhage, CVL cerebrovascular lesion, SDH subdural hematoma, CA cerebral anoxia following cardiac arrest, TBI traumatic brain injury, CI cerebral infarction, SAH subarachnoid hemorrhage. a As entered on the organ donor form. b Islet sample excluded from all analyses due to low quality reads in islet EC library. Exocrine sample used in PCA, heatmap, and comparisons between donors with normoglycemia and donors with IGM.  www.nature.com/scientificreports/ material, which was thus substituted by exocrine flow-through material. These aberrations did not noticeably affect the summarized data.
Library preparation and transcriptome sequencing. Libraries were prepared for each sample with the Ion AmpliSeq Transcriptome Human Gene Expression Kit. Sequencing of sample pools was then done using the Ion 550 Kit and the S5 XL (IonS5XL, Thermo Fisher) instrument with a 550 chip in 3 separate runs. The primer pool hg19_AmpliSeq_Transcriptome_21K_v1, containing primers for 20,813 amplicons, was used. Average mean read length for each run was 112, 112 and 113 bp respectively. The produced reads were aligned to hg19 AmpliSeq Transcriptome ERCC v1. The mapping rates were 97%, 96% and 96%.
One sample (the islet EC sample of donor Normo8) produced low quality reads. This sample was thus excluded from data analysis. 13 other samples from pancreatic arteries were sequenced alongside the MVEC samples (possibly affecting the average read lengths and mapping rates of the three runs), but were not included in any data analyses.
Data normalization, exploration and differential gene expression testing. The raw counts were analyzed using the edgeR 34,35 package (v3.28.1) in R (v.3.5.3). The filterByExpr function was used to remove any genes that did not have the equivalent of at least 10 CPM (counts per million) in at least 5 samples (the size of the smallest experimental group). Raw counts were normalized using TMM (trimmed mean of M) values 36,37 . TMM adjusted and log-normalized CPM were used for principal component analysis using the R-package PCAtools v.1.2.04 38 . A heatmap with hierarchical clustering was constructed for the 1000 genes with the highest variance using the pheatmap package in R and the calculated z-scores of the logCPM-TMM normalized counts.
Differential gene expressions were assessed for tissue origin and glycemic status using the edgeR functions glmQLFit and glmTreat (quasi-likelihood, QL). Genes with a fold change of > ± 1.2 and a FDR-value < 0.05 were considered differentially expressed at a significant level. 24 samples from 12 donors were used for differential gene expression analysis using a paired generalized linear model. The exocrine sample from Normo8 was excluded as it lacked a paired islet sample.
Gene set enrichment analysis. Gene set enrichment analysis (GSEA) 39,40 was performed with the GSEA software v4.0.3 (Broad Institute, USA), using a gene list ranked with the negative log of the p-value multiplied by the sign of the fold change for the TREAT-analysis between islet and exocrine samples as the ranking metric. 1000 gene permutations and a weighted gene enrichment statistic was used for the analysis. Gene sets containing between 15 and 500 genes from GO 41,42 biological processes (ftp.broadinstitute.org://pub/gsea/gene_sets/c5.bp. v6.2.symbols.gmt). To test for differences between donors with normoglycemia and donors with IGM, TMM and CPM normalized counts were used and 1000 phenotype permutations were performed.
Gene sets with an FDR q-value < 0.25 and a nominal p-value < 0.01 were considered significantly enriched. Leading edge analysis was performed on gene sets with FDR < 0.05.
Visualization, clustering and annotation of the GSEA result was performed using the Enrichment map suite 43 within Cytoscape 3.7.2 software 43,44 . A nominal p-value cutoff of 0.005 and a FDR-cutoff of 0.01 was used with a Jaccard + Overlap combined model. Gene set annotations were manually curated. Comparison with murine endothelial cell atlas. Single-cell RNA-Seq data sets from a murine endothelial cell atlas 2 were converted from MGI to HGNC using the R-package biomaRt 46 . 14 Gene sets corresponding to 11 tissues of origin (heart, brain, liver, extensor digitorum longus, kidney, small intestine, colon, lung, soleus, spleen, testis) for capillary endothelial cells and 3 different biological processes (angiogenic, proliferation, interferon) were constructed from the genes in each set that overlaps with the genes remaining in y after the filtering step. Pre-ranked GSEA was then used to allow comparison with our data sets. Statistics. Differential gene expression analysis: Differential gene expression was assessed using a paired generalized linear model and a quasi-likelihood function test that tests for a fold change of > ± 1.2. FDR-values were calculated using the Benjamini-Hochberg method as applied in edgeR. Genes with a FDR-value < 0.05 were considered to differentially expressed at a significant level.
GSEA: A ranked gene list was created using the negative log of the p-value multiplied by the sign of the fold change for the TREAT-analysis between islet and exocrine samples as the ranking metric. 1000 gene permutations and a weighted gene enrichment statistic was used for the analysis. Gene sets with an FDR q-value < 0.25 and a nominal p-value < 0.01 were considered significantly enriched. Leading edge analysis was performed on gene sets with FDR < 0.05. For clustering and visualization, a nominal p-value cutoff of 0.005 and a FDR-cutoff of 0.01 was used with a Jaccard + Overlap combined model.

Data availability
Raw counts and TMM-CPM-log2-normalized counts are available on GEO (GSE157546).