scRNAseq data integration reveals diverse tumor microenvironmentTo conduct this study, primary tumor samples (n = 51), metastatic samples (n = 6), and adjacent normal (control) tissue (n = 6) from a small subset of patients were obtained from five public scRNAseq datasets. After filtering for low quality cells or genes, a total population of 98,749 cells and 35,677 genes remained for downstream analysis, revealing an atlas of the TME (Fig. 1a). The major clusters were manually annotated according to differentially expressed (DE) canonical marker genes (Fig. 1b, c): T cells (n = 27,862) were positive for CD3E, CD4, CD8A, and GZMB; B-cells (n = 5071) expressed MS4A1 and CD79A; plasma B (n = 1516) had high expression of CD79A and MZB1; acinar cells (n = 1694) exhibited PRSS1 and REG1A, epithelial cells (n = 29,781) expressed EPCAM and KRT8; endothelial (n = 1803) showed PECAM1 and CDH5; mast cells (n = 2010) were positive for TPSAB1, myeloid cells (n = 21,645) expressed CD16b, CD68, and CD14; and stromal cells (n = 5564) were defined by COL1A1 and LUM. We then combined epithelial with acinar cells, mast cells with myeloid cells, B cells with plasma B cells, and endothelial cells with stromal cells to form five major groups for downstream analysis, with calculated total counts (Fig. 1d) and patient type composition (Fig. 1e), revealing a diverse tumor population composed primarily of epithelial, T-cells, and myeloid cells. These preliminary results show the abundance of immune cells in the TME; however, the poor prognosis and failures of current immunotherapies suggest mechanisms that prevent an effective anti-tumor immune response.Figure 1Canonical cell markers in PDAC TME. (a) Uniform Manifold Approximation Projection (UMAP) of the comprehensive tumor microenvironment, with a total of 63 samples, 98,749 cells and 35,677 genes. (b) Dot plot of canonical cell marker genes used to annotate major cells. (c) UMAP of canonical marker genes. (d) Bar plot of number of cells per cell type. (e) Bar plot of percentages of each condition per cell type.Characterization of T-cells in the PDAC TMEFirst, we investigated the composition of T cells in the TME, as they are key parts of the adaptive immunity response against tumors. We labeled 11 groups under the general classification of CD4+/CD8+ T-cells and GZMB+ natural killer (NK) cells, and used reported genes in literature to identify specific T cell subtypes (Fig. 2a, b). CD4+ T cells were comprehensively found in greater proportions in tumor populations (Fig. 2c), with CD4+ central memory (Tcm) (Diff = 7.2%, p = 0.031), FOXP3+ T regulatory (Treg) (Diff = 7.6%, p = 0.035), Naive CD4+ (Diff = 13.2%, p = 0.0021), and T helper 2 (Th2) (Diff = 3.6%, p = 0.0059) all enriched. Populations of FOXP3+ Treg have been confirmed mediators in angiogenesis and immune-suppressive functions, contributing to the pro-tumor TME25, although a recent study stated the role of Th2 cells in anti-tumorigenic responses26. Conversely, the majority of CD8+ T cells were found in lower levels in tumor samples, with activated CD8+ (Diff = 26.7%, p = 0.0026) and CD8+ resident memory (Trm) (Diff = 4.3%, p = 0.0095) being significant in adjacent normal tissue (Fig. 2c). CD8+ effector memory (Tem) (Diff = 8.43%, p = 0.012) was the only CD8+ enriched group in the tumor. Although studies have logically shown the increased presence of exhausted CD8+ T cells (Tex) in tumors due to the immunosuppressive environment27,28, we found no significant difference in Tex proportions between healthy adjacent and tumor (Diff = 3.2%, p = 0.088), with even activated CD8+ cells in the tumor surpassing the Tex proportion (Fig. 2c). This suggests that the infiltration and effectiveness of CD8+ T cells in the tumor is extremely low, hence its classification as a “cold tumor” owing to the failures and challenges of anti-PD1/PDL1 therapies in PDAC29. The increased number of CD4+ Treg and reduced CD8+ anti-tumor responses indicate the prevalent pro-tumor mechanisms in the TME, and more research is needed to uncover these biological processes.Figure 2Proportions of T cells in TME underscore low CD8+ T cell infiltration and activation. (a) UMAP characterization of T cell populations in the TME. (b) UMAP of defining marker genes in T-cells populations. (c) Bar plot of percentage per group per patient. CD4+ Tcm (p = 0.031), FOXP3+ Treg (p = 0.035), Naive CD4+ (p = 0.0021), Th2 (p = 0.0059), activated CD8+ (p = 0.0026), and CD8+ Trm (p = 0.0095) are statistically significant using Mann–Whitney U Test.Characterization of stromal cells in the TMEAs stromal cells have been known to contribute to an immunosuppressive TME while playing a critical role in the production of extracellular matrix in PDAC30, we decided to investigate this cell type further to classify specific identities of these populations. Altogether, 12 populations of stromal cells were identified based on previously described gene markers (Fig. 3a, Supplementary Fig. 2A), including PECAM1+ endothelial cells (n = 1803). Smooth muscle cells (n = 1274) and pericytes (n = 523) were classified by high expression of MYL9 and α-SMA (ACTA2), which suggests that α-SMA may not be suitable as a hallmark for myCAF given its high expression across multiple groups (Supplementary Fig. 2B). We therefore distinguished our large subset of myCAF based on their expression of matrix-associated genes (collagens, proteoglycans, and matrix metalloproteinases) and CTHRC1, which were then further divided into two clusters: canonical myCAF (n = 444) with high CTHRC1 expression, and CTHRC1+GREM1+ myCAF (n = 1413). CLU+ fibroblasts (n = 515) resembled a smooth muscle phenotype with expression of CLU and ADIRF. A large population of fibroblasts expressed the complements C3 and C7, suggesting their proinflammatory nature; we differentiated them into three subtypes, C3+RARRES1+ CAF (n = 873), C3+SFRP1+ CAF (n = 264), and a group that exhibited both an inflammatory and myofibroblast signature, which we designated C3+CTHRC1+ CAF (n = 1413). Mesothelial cells acquired high expression of KRT8 and KRT18 (n = 90). A small subset of fibroblasts were shown to be antigen-presenting fibroblasts based on expression of MHC-II and RGS5 (n = 78), while another group expressed MKI67, represented as proliferative CAF (n = 30).Figure 3CTHRC1+GREM1+ myCAF contribute to fibrosis and pro-tumorigenic ECM production and remodeling. (a) UMAP characterization of stromal cell populations in the TME. (b) Bar-plot of percentage per group per patient. CTHRC1+GREM1+ myCAF (p = 0.0014) and CLU+ CAF (p = 0.025) are statistically significant using Mann–Whitney U-Test. (c) UMAP of DE genes in CTHRC1+GREM1+ myCAF. (d) Violin plots of proteoglycan and TGF-β signaling enrichment signatures. (e) GSEA plot of EMT pathway in CTHRC1+GREM1+ myCAF. (f) Bar plot of upregulated pathways in CTHRC1+GREM1+ myCAF using Enrichr. (g) Kaplan–Meier survival curve analysis of CTHRC1+GREM1+ myCAF (TCGA cohort).
CTHRC1
+
GREM1
+ myCAF contribute to fibrosis, epithelial mesenchymal transition and are linked to poor survivalTwo stromal cell subtypes, CTHRC1+GREM1+ myCAF (Diff = 23.7%, p = 0.0014) and CLU+ CAF (Diff = 5.9%, p = 0.025) were found in significantly higher proportions in tumor patients (Fig. 3b), while other groups such as smooth muscle (Diff = 10.4%), C3+RARRES1+ CAF (Diff = 9.7%), and canonical myCAF (Diff = 3.5%) were trending towards tumor-enriched. In particular, CTHRC1+GREM1+ myCAF expressed high levels of collagens (Fig. 3c, Supplementary Fig. 2B, C), including type-1 (COL1A1, COL1A6), type-3 (COL3A1), type-5 (COL5A2), and type-6 (COL6A1), which contribute to tumor cell migration and proliferation through collagen/integrin interactions31,32. These fibroblasts also expressed matrix metalloproteinases MMP2, MMP11 and MMP14 (Supplementary Fig. 2B), which are not only key components of matrix remodeling33, but lead to cancer progression and worse survival in patients34,35. We also observed increased expression of proteoglycans such as CTHRC1 and THBS2 that were unique to myCAF (Fig. 3c, d), which have both been determined to contribute to EMT through the Wnt/b-catenin and NF-kB pathways, leading to tumor invasiveness and poor survival36,37,38. In line with these findings, expression of CTHRC1 was upregulated in tumors from all cancers in the TCGA cohort (Supplementary Fig. 2D), suggesting these ECM protein networks are not limited to PDAC. FN1 was also an important gene in CTHRC1+GREM1+ myCAF (Fig. 3c), which has been shown to promote angiogenesis and metastasis of tumor cells through integrin signaling, leading to the activation of the FAK pathway and also contributing to EMT39. Interestingly, ITGB1 and ITGB5 were discovered to be present on these myCAF (Supplementary Fig. 2B), implying the role of these two integrins as a vital form of crosstalk between myCAF and ECM proteins (Supplementary Fig. 2E)40,41.To elucidate the biological pathways of this fibroblast phenotype, gene set enrichment analysis revealed that CTHRC1+GREM1+ myCAF contributed towards increased levels of TGF-β signaling (Fig. 3d). Further analysis revealed that EMT was an extremely significant mechanism regulated by CTHRC1+GREM1+ myCAF (Fig. 3e, f, Supplementary Fig. 2F), indicating these fibroblasts support cancerous cell differentiation and proliferation leading to rapid tumor growth and potential metastasis. In addition, these cells were also found to help promote hypoxia (ANXA2, SDC2, LOX) and angiogenesis (POSTN, VCAN, LUM) (Fig. 3f). To evaluate the clinical effect on patients, we extended our study to include the TCGA cohort, finding that patients with a high signature of the top 200 DE genes of CTHRC1+GREM1+ myCAF were shown to have significantly worse survival, implicating these cells as not only pro-fibrotic but pro-tumorigenic (Fig. 3g).Characterization of myeloid cells in the TMEAlthough fibroblasts are the main contributor to ECM deposition and remodeling, the function of myeloid cells as crucial components of the pro-tumorigenic TME in PDAC reveals the need to investigate these complex interactions further42. Altogether, we classified our myeloid subset into 15 clusters (Fig. 4a, Supplementary Fig. 3A), including a known population of TPSAB1+ mast cells (n = 2010). Neutrophils were categorized into three groups, with high expression of CD16b and CD62L throughout all cells; the first expressed high levels of interferons (IFIT2 and IFIT3) which were labeled as IFN+ neutrophils (n = 1647). Another group expressed high levels of MMP9 and MMP25, which we labeled as MMP9+ neutrophils (n = 356), while the last group was called GMFG+ neutrophils (n = 2610). Four subtypes of monocytes were also found, which were labeled as CD14+CD16- monocytes (n = 1337), CD16+ monocytes (n = 440), ITGB2+ monocytes (n = 471), and IL1B+ monocytes (n = 2880). A main group of classical dendritic cells (cDC) expressed high levels of MHC-II markers and CD74 (n = 1856), and interestingly a small population of plasmacytoid DC were also noted by their expression of CLEC4C and IL3RA (n = 148). Macrophage populations characterized by CD68 expression composed of the largest proportion of total myeloid cells, which were labeled into four groups: monocyte-like macrophages expressed CD16a and MHC-II (n = 2256), C1Q-high macrophages expressed C1QA and C1QB (n = 2748), SPP1+APOE+ macrophages (n = 3226), and SPP1+VEGFA+ macrophages (n = 1320). A small proportion of proliferating macrophages expressing MKI67 (n = 176) were also identified.Figure 4SPP1+APOE+ TAM exhibit pro-fibrosis and pro-tumor properties. (a) UMAP characterization of myeloid cell populations in the TME. (b) Bar plot of percentage per group per patient. SPP1+APOE+ TAM (p = 0.011) and IFN+ (p = 0.046) TAN are statistically significant using Welch’s T-Test. (c) UMAP plots of significant DE genes in SPP1+APOE+ TAM. (d) UMAP plots of M1 and M2 macrophage signatures. (e) Violin plot of EMT enrichment signature. (f) GSEA plot of EMT pathway. (g) Bar plot of upregulated pathways in SPP1+APOE+ TAM using Enrichr. (h) Kaplan–Meier survival curve for SPP1+APOE+ TAM (TCGA cohort).
SPP1
+ TAM enriched in tumor and contribute to pro-tumorigenic functions and lead to worse survivalWe compared proportions of each myeloid subtype in each patient, finding that SPP1+APOE+ TAM (Diff = 9.5%, p = 0.011) and CD62L+IFN-high neutrophils (Diff = 5.0%, p = 0.046) were significantly enriched in primary tumor patients versus adjacent normal (Fig. 4b). SPP1+VEGFA+ TAM were also trending towards increased proportions in tumor samples (Diff = 3.4%, p-value = 0.21). All SPP1+ TAM exhibited the M2-marker CD68 with other marker genes such as CSTB and NUPR1 (Fig. 4C), as well as being enriched for M2 signature under the M1 and M2 classification (Fig. 4d)43; additionally, SPP1+APOE+ TAM exhibited genes associated with ECM remodeling, such as FN1, MMP14, and LGALS1 (Supplementary Fig. 3B) that are also indicative of a pro-fibrotic M2 polarization (Fig. 4d, Supplementary Fig. 3C). As the role of FN1 in ECM signaling and EMT is paramount, we noticed that SPP1+ TAM also exhibited a high EMT signature (Fig. 4e–g), with FN1, LRP1, PLAUR, and TIMP1 highly expressed (Supplementary Fig. 3B) and identified as important components leading to EMT44,45. Moreover, pathway analysis revealed that SPP1+APOE+ TAM were regulated by the HIF-1 signaling pathway (Supplementary Fig. 3C, 3D), inducing expression of VEGFA, LDHA, and ALDOA46,47, which not only promotes hypoxia (Fig. 4g) but is another approach to induce EMT in tumor cells48. We also identified SPP1+APOE+ TAM as major downstream components of mTORC1 signaling (Fig. 4g, Supplementary Fig. 3C), which has important roles in macrophage polarization, tumor metabolism, and protein synthesis49, possibly hinting that mTORC1 is a viable candidate for inducing M2 macrophages in PDAC. Survival analysis conducted in the TCGA cohort showed that the top 200 DE gene signature of these TAM led to worse survival for patients (Fig. 4h). Thus, we have established the monumental impact of these macrophages in the immuno-suppresive, pro-fibrotic niche, which we hypothesize may help support the role of CTHRC1+GREM1+ myCAF and tumor cells further through EMT and other pro-tumor mechanisms. Together, these interactions reveal the multi-faceted roles of SPP1+ TAM in the PDAC TME.Expressions of CTHRC1
+
GREM1
+ myCAF and SPP1
+
APOE
+ TAM correspond with worse survivalBecause of the significant contributions of CTHCR1+GREM1+ myCAF and SPP1+APOE+ TAM to a pro-tumorigenic TME, we investigated the potential of a synergistic relationship. To address this hypothesis, we first conducted expression plots of CTHRC1, GREM1, SPP1, and APOE in the TCGA cohort (Fig. 5a, Supplementary Fig. 4A), which were all significantly enriched in tumor samples versus control. To determine if these two groups were related, Spearman correlation (Fig. 5b) revealed a remarkable correlation between the myofibroblast and macrophage signatures (r = 0.87), while the relationship between CTHRC1 and SPP1 alone was also deemed significant (r = 0.25). Survival analysis of CTHRC1 and SPP1 gene expression in PAAD demonstrated that patients with high expression of these genes led to worse survival (Fig. 5c), demonstrating the pro-tumor functions of these genes in the clinical context. Furthermore, the combined gene signature of both groups also resulted in worse prognosis for patients, as well as being enriched in tumor populations (Fig. 5c, d). These results highlight the pro-tumorigenic and positive correlation of CTHRC1+GREM1+ myCAF and SPP1+APOE+ TAM in the TME in single-cell resolution. Next, we aimed to demonstrate their relationship using spatial transcriptomics, amplifying research significance if found to be consistent with single cell transcriptomics.Figure 5SPP1+APOE+ TAM and CTHRC1+GREM1+ myCAF are positively correlated and contribute to worse prognosis in cancer patients. (a) Box plots of CTHRC1 and SPP1 gene expression in normal versus tumor patients (TCGA cohort). (b) Spearman correlation graphs of CTHRC1 and SPP1, and of the CTHRC1+GREM1+ myCAF signature and SPP1+APOE+ TAM signature. (c) Overall survival (OS) curves of CTHRC1+GREM1+ myCAF, SPP1+APOE+ TAM, and the combined signatures of both. (d) Expression boxplot of the combined signature.Spatial transcriptomics reveals co-localization of CTHRC1
+
GREM1
+ myCAF and SPP1
+
APOE
+ TAMTo assess cell interactions in the spatial landscape, we conducted spatial transcriptomics on three public datasets acquired from PDAC patients (Fig. 6a, Supplementary Fig. 5A). A large population of ductal cells were identified based on expression of PRSS1 and REG1A; populations of alpha cells (GCG) and beta cells were also observed (INS), but were all categorized under the ductal population. We also identified mixed cellular compartments, including a fibroblast/ductal population characterized by high expression of COL1A1, ACTA2, and REG1A; a fibroblast/malignant population expressed GREM1, TIMP1, and TFF1, and a mixed immune population had markers CD3E and MS4A1. In addition to normal cells of the pancreas, tumor cells were determined based on their expression of EPCAM, CDH1, and ID1, while an epithelial cluster was established based on the unique expression of MUC4, MUC6, and MUC5AC. We also noticed a substantial population of CTHRC1+GREM1+ myofibroblasts which expressed high levels of ECM genes (CTHRC1, GREM1, FN1, POSTN), aligned with scRNAseq results (Fig. 6b, Supplementary Fig. 5B). Additionally, we determined the aforementioned SPP1+ macrophages, with expression of SPP1, CD68, MARCO, and FN1 (Fig. 6b). We verified this by plotting the CTHRC1+GREM1+ myCAF, SPP1+APOE+ TAM, and EMT signatures from previous parts of our study (Fig. 6c, Supplementary Fig. 5C), which showed clear overlays, supporting the co-localization of fibroblasts and macrophages while promoting the pro-tumor TME. Furthermore, the myCAF and EMT signature were almost identical, spatially confirming the presence of CTHRC1+GREM1+ myCAF in supporting EMT. Spearman correlation between the myCAF and TAM signatures (r = 0.64) further confirmed the synergistic relationship between these two groups in the spatial landscape (Fig. 6d). Since all three groups of fibroblasts, macrophages, and tumor cells were identified in close proximity to each other (Fig. 6a), we then quantitatively visualized these spatial localizations with a proximity enrichment heatmap (Fig. 6e), showing that CTHRC1+GREM1+ myCAF were only in relationship to SPP1+ TAM and tumor cells, while SPP1 + TAM were in close proximity to tumor cells and also to immune cells.Figure 6Spatial analysis confirms co-localization of CTHRC1+GREM1+ myCAF and SPP1+ TAM. (a) Representative histopathological H&E (FFPE) staining of pancreas tissue cross-section. UMAP of Leiden clusters shows colocalization of CTHRC1+GREM1+ myCAF and SPP1+ TAM. (b) Spatial gene expression of important marker genes. (c) Spatial enrichment of SPP1+ TAM, CTHRC1+GREM1+ myCAF, and EMT pathway signatures based on scRNAseq. (d) Spearman correlation graph of SPP1+ TAM and CTHRC1+GREM1+ myCAF signatures. (e) Neighborhood enrichment correlation matrix to quantitatively visualize spatial cellular cluster proximity. (f) Interaction weights of SPP1+ macrophage signaling to other cell types. (g) Dot plot highlights SPP1 ligand-receptor interactions from macrophages to other cells.Discovering the intricate crosstalk between CTHRC1
+
GREM1
+ myCAF, SPP1
+ TAM and tumor cellsTo understand how these three populations interact, cellular ligand-receptor analysis revealed that SPP1 expressed by macrophages interacted with tumor cells, fibroblasts, and endothelial cells the most (Fig. 6f). SPP1 bound with integrins such as ITGB1 and ITGB5 (Fig. 6g), implying the role of SPP1 in integrin signaling between fibroblasts, macrophages, and tumor cells that lead to EMT and cell-adhesion to the ECM. The SPP1-CD44 axis was also significantly enriched in macrophage crosstalk (Fig. 6g), corresponding to an important pathway for cell surface adhesion and metastasis through the activation of PI3K/Akt and MAPK signaling50. SPP1+ macrophages also exemplified high expression of TGFB1 (Supplementary Fig. 5B), suggesting an increase in TGFβ signaling towards tumor cells and fibroblasts that directly contribute towards fibroblast recruitment, cell proliferation and EMT, although TGFβ1 was high across all groups. Similarly, CTHRC1+GREM1+ myCAF communicated with SPP1+ macrophages and tumor cells through integrin signaling (Supplementary Fig. 5D), where collagen/integrin pairs and FN1/integrin pairs were prevalent. These interactions not only enhance cell proliferation and transformation51, but could serve as a potential factor of polarization and recruitment of macrophages to the TME. Together, these interactions reveal the pathways between fibroblasts, macrophages, and tumor cells that contribute to tumor differentiation, proliferation, and worse prognosis.In-silico knockout of osteopontin receptors in CTHRC1
+
GREM1
+ myCAFIn line with our previous findings from spatial ligand-receptor analysis, we then conducted an in-silico knockout of three potential osteopontin receptors—CD44, ITGB5, and ITGB1 in the CTHRC1+GREM1+ myCAF population to validate the impact of SPP1-CD44, SPP1-ITGB5, and SPP1-ITGB1 crosstalks (Fig. 6g, Supplementary Fig. 5D). Given the pro-fibrotic properties of SPP1, we hypothesized observational changes related to ECM associated genes following the virtual KO. Indeed, the top 50 perturbed genes resulting from CD44 KO were closely associated with ECM organization and interaction, including genes such as CTHRC1, FN1, THBS2, SPARC, and various collagens (Fig. 7a). These outcomes suggest that the SPP1-CD44 pair significantly affects the fibrotic functions of myCAF, as the absence of CD44 introduces major perturbations in genes related to fibrosis and ECM. We then repeated our analysis with ITGB5 and ITGB1, with similar results (Fig. 7b, Supplementary Fig. 6A). The loss of ITGB5 affected FN1, MMP11, SDC1, and MMP14, although CTHRC1 was not among the top 50 perturbed genes. ITGB1 KO resulted in changes to CTHRC1, VCAN, and LGALS1. Nevertheless, the convergence of virtual gene perturbations resulting from the knockouts of all three genes underscores the multiple mechanisms of osteopontin secreted by SPP1+APOE+ TAM in inducing fibrotic CTHRC1+GREM1+ myCAF.Figure 7In-silico KO of CD44 and ITGB5 in CTHRC1+GREM1+ myCAF. (a) Bar plot of top GSEA terms from CD44 KO perturbed genes. Egocentric plot of CD44 KO perturbed genes, along with corresponding GSEA terms. (b) Bar plot of top GSEA terms from ITGB5 KO perturbed genes. Egocentric plot of ITGB5 KO perturbed genes, along with corresponding GSEA terms.
CTHRC1+ fibroblasts and SPP1+ macrophages synergistically contribute to pro-tumorigenic tumor microenvironment in pancreatic ductal adenocarcinoma
