Human Pancreatic Cancer Single-Cell Atlas Reveals Association of CXCL10+ Fibroblasts and Basal Subtype Tumor Cells

Abstract

Purpose: Pancreatic ductal adenocarcinoma (PDAC) patients with tumors enriched for the basal-like molecular subtype exhibit enhanced resistance to standard-of-care treatments and have significantly worse overall survival compared with patients with classic subtype-enriched tumors. It is important to develop genomic resources, enabling identification of novel putative targets in a statistically rigorous manner.

Experimental design: We compiled a single-cell RNA sequencing (scRNA-seq) atlas of the human pancreas with 229 patient samples aggregated from publicly available raw data. We mapped cell type-specific scRNA-seq gene signatures in bulk RNA-seq (n = 744) and spatial transcriptomics (ST; n = 22) and performed validation using multiplex immunostaining.

Results: Analysis of tumor cells from our scRNA-seq atlas revealed nine distinct populations, two of which aligned with the basal subtype, correlating with worse overall survival in bulk RNA-seq. Deconvolution identified one of the basal populations to be the predominant tumor subtype in nondissociated ST tissues and in vitro tumor cell and patient-derived organoid lines. We discovered a novel enrichment and spatial association of CXCL10+ cancer-associated fibroblasts with basal tumor cells. We identified that besides immune cells, ductal cells also express CXCR3, the receptor for CXCL10, suggesting a relationship between these cell types in the PDAC tumor microenvironment.

Conclusions: We show that our scRNA-seq atlas (700,000 cells), integrated with ST data, has increased statistical power and is a powerful resource, allowing for expansion of current subtyping paradigms in PDAC. We uncovered a novel signaling niche marked by CXCL10+ cancer-associated fibroblasts and basal tumor cells that could be explored for future targeted therapies.

Figure 1.

scRNA-seq. A, UMAP visualization of (n = 6) normal donor pancreata, (n = 26) adjacent/normal pancreas tissues, (n = 172), primary pancreatic tumors, and (n = 25) metastatic biopsies from the (n = 21) liver, (n = 1) lung, (n = 1) vaginal apex, (n = 1) omentum, and (n = 1) peritoneum. TNK denotes T and natural killer cells. B, Bar plot depicting the relative distribution of coarsely labeled cell types across the different disease states. C, UMAP visualization of the KRT18⁺/KRT19⁺ ductal cluster split by disease state. D, Bar plot showing the distribution of ductal cell types by disease state. E, UMAP visualization showing the expression of KRT18 and AMY2A in all disease states. F, Violin plots of selected (among top 30) significantly upregulated genes comparing each cluster. G, UMAP visualization showing the distribution of CNSs across the ductal cell clusters in all disease states, with a putative tumor threshold of 0.15. H, UMAP visualization showing reclustering of acinar cells by cluster (top) and by initial label (bottom). I, UMAP visualizations showing expression for acinar, ductal, ADM, and PanIN markers on the acinar cell reclustering.

Figure 2.

Tumor-enriched ductal clusters show high levels of within and across sample heterogeneity which correlates with survival differences in bulk RNA-seq. A, UMAP visualization showing the reclustered ductal cells from the PT- and Met-enriched clusters, including all disease states. B, Dot plot showing selected (among the top 30) significantly upregulated genes comparing each cluster. C, UMAP visualization showing the inferred CNSs across the different tumor enriched ductal clusters in all disease states, with a putative tumor threshold of 0.15. D, Bar plot showing the relative frequency of the tumor-enriched ductal clusters per patient. E, Violin plots signature scores from the basal and classic gene lists from Raghaven and colleagues (7) highlighting the enrichment of the basal clusters 4 and 7. F, UMAP visualization showing select classic genes and basal genes in combined disease states. G, Dot plot showing top enriched pathways for each tumor enriched ductal cluster using ReactomeDB pathways. H, Violin plot showing the expression of TP63 in the nine ductal clusters. I, Forest plot showing HRs and 95% confidence intervals for gene signatures for each ductal cluster in four bulk RNA-seq datasets. SLC, solute carrier family; ER, endoplasmic reticulum; NGF, nerve growth factor; MET, mesenchymal-epithelial transition; TNPS, TnpB nuclease dead repressors.

Figure 3.

Fibroblast and myeloid heterogeneity identified in the PDAC microenvironment. A, UMAP visualization showing fibroblast subtypes in combined disease states. Due to the inclusion of donor and AdjN pancreata samples, we refer to fibroblasts as “Fb” instead of cancer associated fibroblasts (CAFs). B, UMAP visualization showing two pan-fibroblast genes (DCN, LUM, COL1A1, and PDPN). C, UMAP visualization showing inflammatory markers (C7 and DPT) and myofibroblast markers (COL11A1 and ACTA2) in fibroblast clusters in all disease states. D, Dot plot showing selected (among the top 30) significantly upregulated genes in each fibroblast subtype. E, Dot plot showing top enriched pathways for each fibroblast cluster using ReactomeDB pathways. F, UMAP visualization showing CXCL1, 2, 10, and 14 chemokines in fibroblast clusters in combined disease states. G, Bar chart showing distribution of fibroblast subtypes by patient, separated by disease state. H, UMAP visualization showing coarsely labeled myeloid subtypes in all disease states. I, Dot plot showing selected (among the top 30) significantly upregulated genes in each myeloid cell type. J, Bar chart showing the distribution of myeloid subtypes by patient separated by disease state. ER, endoplasmic reticulum; HIV, human immunodeficiency virus; CIT, citron rho-interacting kinase; PAK, p21-activated kinase.

Figure 4.

Associations between fibroblast, myeloid, and ductal clusters are identified in bulk RNA-seq. A, Correlation heatmaps showing hierarchical clustering–based groupings of fibroblast, myeloid, and ductal clusters in four PDAC bulk RNA-seq datasets (TCGA-PAAD, GSE71729, Puleo, and ICGC). B, Scatter plots depicting the correlation between the proportion of Ductal populations and the proportion of the CXCL10+ myFb population in the scRNA-seq samples with p-value and correlation coefficient reported. C, UMAP visualization showing expression of the CXCL10 receptor, CXCR3, in ductal, myeloid, and T and NK cells. Additionally, a violin plot showing expression of IFNγ in coarse labels from all cells, UMAP visualizations showing expression of IFNγR1 and IFNγR2 in fibroblasts, and violin plot showing the expression of CXCR3 in each of the nine ductal clusters. D, Grouped bar plot showing the distributions of normalized gene signature scores for ductal, myeloid, and fibroblast cells in the four PDAC bulk RNA-seq datasets. ICGC, International Cancer Genome Consortium.

Figure 5.

Associations between fibroblast, myeloid, and ductal clusters in ST. A, Deconvolved ST showing the existence and colocalization of fibroblast, myeloid, and ductal clusters in 22 ST samples. B, Higher magnification H&E of sample S2A_9 highlighting ductal 1 and 8 in two PanIN regions (R1 and R2). C, Higher magnification regions of S1A_16 highlighting ductal cell heterogeneity within sample, in which R1 is well differentiated with classic cell types whereas R2 and R3 are poorly differentiated with basal cell types. D, Bar chart showing the distribution of each ductal cluster in the ST by sample. E, Bar chart showing the distribution of the fibroblast clusters in the ST by sample. F, Bar chart showing the distribution of myeloid clusters in the ST by sample. G, Boxplots showing the average number of fibroblast neighbors for ductal 4 and ductal 1 in the ST. H, Ductal cell type deconvolution of six pancreatic cancer cell lines bulk RNA-seq from Bryant and colleagues (66). I, Ductal cell type deconvolution of 49 pancreas organoid bulk RNA-seq containing 11 normal and 38 PT samples.

Figure 6.

In situ validation of ductal 4 and ductal 7 in human PDAC tissues. A, Key showing the markers for expression. B, Bar chart showing relative proportions of ductal 4, ductal 7, transitory populations, and double-negative ductal cells (likely classic). C, Multiplex immunostaining of LAMC2 (green) for ductal 4 and S100A2 protein (red) for ductal 7 in situ. D, Multiplex immunostaining of LAMC2 (green) for ductal 4 and S100A2 protein (red) for ductal 7 in one tumor samples. E, Multiplex immunostaining of LAMC2 (green) for ductal 4 and S100A2 protein (red) for ductal 7 in situ. DAPI in blue for identification of nuclei and PANCK protein (white) to identify epithelial cells.