Single-Cell Meta-Analysis Uncovers the Pancreatic Endothelial Cell Transcriptomic Signature and Reveals a Key Role for NKX2-3 in PLVAP Expression

Abstract

Background: The pancreatic vasculature displays tissue-specific physiological and functional adaptations that support rapid insulin response by β-cells. However, the digestive enzymes have made it difficult to characterize pancreatic endothelial cells (ECs), resulting in the poor understanding of pancreatic EC specialization.

Methods: Available single-nuclei/single-cell RNA-sequencing data sets were mined to identify pancreatic EC-enriched signature genes and to develop an integrated atlas of human pancreatic ECs. We validated the findings using independent single-nuclei/single-cell RNA-sequencing data, bulk RNA-sequencing data of isolated ECs, spatial transcriptomics data, immunofluorescence, and RNAScope of selected markers. The NK2 homeobox 3 (NKX2-3) TF (transcription factor) was expressed in HUVECs via gene transfection, and the expression of pancreatic EC-enriched signature genes was assessed via RT-qPCR.

Results: We defined a pancreatic EC-enriched gene signature conserved across species and developmental stages that included genes involved in ECM (extracellular matrix) composition (COL15A1 and COL4A1), permeability and barrier function (PLVAP, EHD4, CAVIN3, HSPG2, ROBO4, HEG1, and CLEC14A), and key signaling pathways (S1P [sphingosine-1-phosphate], TGF-β [transforming growth factor-β], RHO/RAC GTPase [guanosine triphosphatase], PI3K/AKT [phosphoinositide 3-kinase/protein kinase B], and PDGF [platelet-derived growth factor]). The integrated atlas revealed the vascular hierarchy within the pancreas. We identified and validated a specialized islet capillary subpopulation characterized by genes involved in permeability (PLVAP and EHD4), immune-modulation (FABP5, HLA-C, and B2M), ECM composition (SPARC and SPARCL1), IGF (insulin-like growth factor) signaling (IGFBP7), and membrane transport (SLCO2A1, SLC2A3, and CD320). Importantly, we identified NKX2-3 as a key TF enriched in pancreatic ECs. DNA-binding motif analysis found NKX2-3 motifs in ≈40% of the signature genes. Induction of NKX2-3 in HUVECs promoted the expression of the islet capillary EC-enriched genes PLVAP and SPARCL1.

Conclusions: We defined a validated transcriptomic signature of pancreatic ECs and uncovered their intratissue transcriptomic heterogeneity. We showed that NKX2-3 acts upstream of PLVAP and provided a single-cell online resource that can be further explored by the community: https://vasconcelos.shinyapps.io/pancreatic_endothelial/.

Keywords: capillary permeability; endothelial cells; extracellular matrix; insulin secretion; islets of langerhans; pancreas; single-cell analysis.

Figure 1.

A conserved pancreas endothelial cell (EC)–enriched gene signature exists across species and age. A through C, Uncovering biologically relevant pancreas EC-enriched signature genes. A, Pancreatic EC-enriched signature genes lie at the overlap between pancreatic enriched differentially expressed genes (DEGs) derived from EC atlases and EC-enriched DEGs derived from pancreatic data sets (curated set of EC genes). B, The overlap of significant pancreatic enriched DEGs in the Descartes fetal human EC atlas and Tabula Muris adult mouse EC atlas and the curated set of EC genes. C, The expression of pancreas EC-enriched signature genes across the 2 multiorgan EC atlases. The x and y axes demonstrate the average log2 fold change (log2FC) of the genes in pancreas ECs compared with other tissue ECs for each EC atlas. The normalized percentage difference was calculated by subtracting the percent expression in population 1 (pancreatic ECs) with that of population 2 (the expression in all other tissue ECs) and was weighted by dividing the difference by the percent expression in population 1. The normalized percentage differences are visualized in color scales and dot sizes for each EC atlas. D through F, Pancreas EC-enriched signatures genes were validated independently in the adult human Tabula Sapiens multiorgan atlas. D, Module scores derived by assessing the collective expression of the pancreas EC-enriched gene signature in the different tissue ECs. Each violin is colored based on the average scaled module score for that tissue. E, The individual log-normalized average scaled expression of the top 30 pancreas EC-enriched signature genes in the Tabula Sapiens EC atlas. Genes are ordered based on normalized percentage difference, and hierarchical clustering of tissues is based on Euclidian distances between the average scaled expression of all the signature genes weighted by the percentage expression. F, RNAScope of Nkx2-3, Plcb1, Clec14a, and Heg1 on embryonic day (E) 18.5 mouse pancreas tissue sections. Islet and exocrine regions are identified by immunofluorescence costaining of insulin and glucagon (pink), which stains cells of the islets; vasculature is identified by costaining of Pecam1 and endomucin (neon blue); epithelial cells are stained by E-cadherin (yellow); nuclei are stained by DAPI (4′,6-diamidino-2-phenylindole; blue); and genes of interest are identified by RNAScope probes (red).

Figure 2.

Gene set overrepresentation analysis of pancreas endothelial cell (EC)–enriched signature genes. Gene set overrepresentation analysis was conducted on the signature genes using Enrichr (v3.2) with the following pathway databases: WikiPathway (2023), Reactome (2022), KEGG (Kyoto Encylopedia of Genes and Genomes; 2021), BioPlanet (2019), Panther (2016), and BioCarta (2016). Selected pathways of interest are visualized in a network; genes are shown in blue ellipses and terms in boxes; boxes are colored based on the negative log-adjusted P value from the gene set overrepresentation analysis.

Figure 3.

Integrated atlas of human pancreatic endothelial cells (ECs) identifies specialized capillary subpopulation. Three different human pancreatic EC data were subsetted from their parent data and integrated using Harmony. This includes the fetal human pancreas data set from the Descartes gene expression atlas and the neonatal and adult pancreas data set from Tosti et al. A, The resulting integrated uniform manifold approximation and projection (UMAP). B, The percentage distribution of different subpopulations in each data set. C, The UMAP of the integrated map split by data set. D, The expression of markers used to differentiate the vascular bed subclusters in the integrated EC atlas. E, The expression of selected genes across biological processes of interest in the atlas. F, The UMAP of the integrated atlas colored by pseudotime inferred through trajectory analysis using SCORPIUS. The trajectory is colored in blue. G, The heatmap demonstrating the change in average scaled expression of the top unique 50 markers determined by Wilcoxon rank-sum test for each subpopulation across pseudotime. Genes are hierarchically clustered based on the Euclidian distance between the average scaled expression of the genes. H, The expression of the top 12 genes across pseudotime ranked by importance to trajectory inference using SCORPIUS.

Figure 4.

Pancreas endothelial cell (EC)–enriched gene signature is enriched in islet capillary ECs. A, The module scores of the pancreas EC-enriched gene signature in the different EC subpopulations in the integrated EC atlas. B, Box plot of the collective module score for the pancreas EC-enriched signature in bulk RNA sequencing (RNAseq) islet and exocrine ECs from healthy normoglycemic human donors sequenced by Jonsson et al. Each dot represents 1 biological replicate, defined as 1 donor (n=7). The P value from a paired t test is shown. C, The principal component analysis of the islet and exocrine bulk RNAseq EC samples using only log normalization of counts per million (logCPM) scaled expression values of the signature genes. D, The heatmap of the average scaled expression of the signature genes with samples ordered based on PC1 (principal component 1) scores and genes ordered based on PC1 feature loading scores. E, The individual dotplot of selected genes showing the average scaled expression in the different subpopulations in the integrated atlas. F, Box plots of the same selected genes with the logCPM values of the bulk RNAseq islet and exocrine EC samples. Adjusted P values from the DE analysis of the bulk RNAseq data using DESeq2 are shown. Each N represents 1 donor (n=7).

Figure 5.

Specialized capillary population is located within the islets. Signature genes for vascular subpopulations identified from the integrated atlas were assessed on spatial transcriptomics data from the study by Olaniru et al. Spearman correlations were calculated for the signature scores with deconvoluted islet and exocrine proportions. A, The Fisher Z score of the difference between Spearman correlation of the signature with islet proportions and exocrine proportions. Each N represents 1 biological replicate, defined as 1 tissue sample sequenced (n=4). The line inside each box represents the median, with the box edges indicating the interquartile range (IQR). Whiskers extend to 1.5× the IQR, and outliers are shown as points beyond the whiskers. P values were obtained using ANOVA, followed by Tukey’s HSD (honestly significant difference) post hoc test for pairwise comparisons. B and C, Islet capillary and exocrine capillary scores in mutually exclusive islet and exocrine regions for postconception week 20 sample 2, characterized as >10% of cells in each voxel identified as islet cells and 0% identified as exocrine or vice versa. B shows the corresponding voxels identified as islet (pink) and exocrine (purple). C visualizes the signature scores for islet capillary and exocrine capillary endothelial cells (ECs) in exclusively islet and exocrine voxels. Arrows indicate selected regions that are different between capillary islet and capillary exocrine. D through H, In situ validation of islet EC capillary population in embryonic day (E) 18.5 mouse pancreas. D, The expression of 3 different DEGs that were enriched in the islet capillary population in the integrated pancreas EC atlas. These genes include PLVAP (left), NKX2-3 (middle), and SPARC (right). E, Immunofluorescence (IF) staining of Plvap. F and G, RNAScope images of Sparc (F) and Nkx2-3 (G). For the in situ images, pink shows IF costaining of glucagon and insulin, which highlights the islets from the exocrine tissue for all images; green shows PLVAP IF staining; red shows Sparc and Nkx2-3 RNAscope punctae; white shows Pecam1/VE-cadherin (vascular endothelial cadherin) costaining; and neon shows Pecam1/endomucin costaining to highlight vessels; yellow shows E-cadherin staining. D, Quantification of fluorescence intensity in arbitrary units (AU) for Plvap (left) and RNAscope area of probe binding (μm²) in islet and exocrine capillaries for Sparc (middle) and Nkx2-3 (right). Each N represents 1 biological replicate, defined as 1 vessel (n=59 for islet and n=35 for exocrine regions for Plvap IF, and n=10 for Sparc and Nkx2-3 RNAscope). The line inside each box represents the median, with the box edges indicating the IQR. Whiskers extend to 1.5× the IQR, and outliers are shown as points beyond the whiskers.

Figure 6.

NKX2-3 regulates pancreas endothelial cell (EC) gene expression. A, TF (transcription factor) motif analysis was conducted to find the NKX2-3 DNA-binding motif obtained from JASPER, an open-access database for TF binding profiles (MA0672.1), within +200 to −1000 bp of transcription start site (TSS) of pancreas EC-enriched signature genes. B, Twenty-four signature genes identified with the NKX2-3 motif are shown. C and D, NKX2-3 tagged C terminally to GFP (green fluorescent protein; NKX2-3-GFP) and GFP only was induced transiently in HUVECs (human umbilical vein endothelial cells). C, Immunofluorescence image of HUVECs transfected 24 hours prior with an NKX2-3-GFP at ×40 magnification. NKX2-3 is shown in red, GFP in green, and DAPI (4′,6-diamidino-2-phenylindole) in blue. D, Gene expression was assessed by RT-qPCR (reverse transcription quantitative polymerase chain reaction) 24 hours after transfection. The mean log2 fold change over GFP-only induced HUVEC genes with P<0.05 based on 1-tailed Wilcoxon rank-sum test is shown. Each dot presents 1 biological replicate (n), defined as the result from RT-qPCR conducted on 1 transfection experiment. The error bars represent the SEM (n=5 for PLVAP; n=4 for SMAD7, HSPG2, CLEC14A, and EDN1; and n=3 for SPARCL1).