Single cell RNA sequencing reveals a complex immune landscape with heterogeneous expression of immune checkpoints and ligands in the pancreatic cancer microenvironment.

We performed single cell RNA sequencing (scRNA seq) on 16 PDA samples, including surgical (n=6) and fine needle biopsy specimens (n=10) (Extended Data Fig. 2A). All the patients were treatment-naïve at the time of sample acquisition. We also included three non-malignant pancreas samples (1196, 1258, 19732) from patients undergoing surgery for duodenal adenoma, ampullary carcinoma, or adjacent to PDA, respectively, where an uninvolved portion of pancreas was included in the resection (as determined by pathologic evaluation). To capture a window into the systemic immune response in PDA patients, we also collected peripheral blood mononuclear cells (PBMCs) from these patients and healthy subjects (Extended Data Fig. 2A, right panel).

In total, we sequenced 8541 cells from adjacent/normal samples and 46,244 from PDA, while from the blood samples we sequenced 14,240 cells from 4 healthy subjects, and 55,873 cells from 16 PDA patients. To define and visualize cell subpopulations, we batch corrected our tumor and blood sequencing samples (Extended Data Fig. 2B) and then used unbiased clustering and a dimensionality reduction through Uniform Manifold Approximation and Projection (UMAP) (Fig. 2A and Extended Data Fig 3A). We identified each subpopulation based on published lineage markers (Fig. 2B and Extended Data Fig. 3B). We observed variability in the total immune cell composition of individual tumors, and in the relative abundance of individual immune cell components similar to CYTOF and mfIHC data (Extended Data Fig. 2C and ​and2D2D).

Open in a separate window

Figure 2.

Single cell RNA sequencing reveals heterogenous expression of immune checkpoints in PDA tissue.

(A) UMAP on 3 adjacent/normal pancreas (left) and 16 PDA patient (right) tissues. Populations identified as follows: acinar (pink), epithelial (red), fibroblasts (dark and light teal), CD8⁺ T cells (dark green), CD4⁺ T cells, (neon green), Tregs (light green), NK cells (purple), B cells (light blue), plasma cells (dark blue), mast cells (yellow), macrophages (dark orange), granulocytes (light orange), dendritic cells (brown), endothelial cells (dark pink), and endocrine (dark red). (B) Dot plot of key markers used to define the identified cell populations. Color of dot represents average expression, while the size of the dot represents percent expression. Dot plot represents merged n=3 adj/norm patients and n=16 PDA patients gene expression of lineage markers. (C) Average expression of immune checkpoint ligands and receptors in the identified cell populations in n=16 tumor tissue samples. (D) Average expression of immune checkpoint receptors on CD8⁺ T cells in n=16 PDA patients and n=3 adj/norm patients merged tissues. (E) Pathway annotations from gene set enrichment analysis (GSEA) using the R package clusterprofiler in n=16 PDA samples compared to n=3 adj/norm samples. The color of the bar represents the p-value adjusted for multiple comparisons using the Benjamini Hochberg method. Enrichment score is plotted on the x-axis. (F) Unbiased differential expression between CD8⁺ T cells from adj/norm pancreas (black) and PDA (grey). Significantly up- and down-regulated genes are plotted as average expression per patient. This analysis was performed on all CD8⁺ T cells found in the adjacent/normal and PDA tissue. Disease stage is plotted on the left: resectable (green), locally advanced (light green), borderline resectable (blue), metastatic (pink), N/A (light blue).

We identified abundant pro-inflammatory cells in the PDA microenvironment including CD8⁺ T cells and natural killer cells (NK) (Fig. 2A). To gather insight into the possible mechanisms preventing anti-tumor immune responses, we profiled average expression of immune checkpoint receptors and ligands by cell type, both in tumors and PBMCs of PDA patients (Fig. 2C and Extended Data Fig. 3C). We observed expression of multiple immune checkpoint receptors in T and NK cell subsets, while myeloid populations were enriched for their corresponding ligands (Fig. 2C). CD8⁺ T cells had elevated ICOS, TIGIT, PDCD1 and LAG3, among others, but relatively low expression of CTLA4. CTLA4, as well as all other checkpoints except for LAG3, was high in CD4⁺ T cells. NK cells also had elevated CD47, TIGIT, TNFRSF18 and LAG3, and modest expression of PDCD1. The expression of checkpoint ligands was heterogeneous, with epithelial cells mainly expressing PVR and LGALS9 (encoding for PVR and GALECTIN 9, ligands for TIGIT and HAVCR2(TIM3) respectively). Myeloid and dendritic cells expressed several genes encoding checkpoint ligands, including SIRPA, LGALS9, PVR, and ICOSLG (Fig. 2C). Similarly, in PBMC samples, CD4⁺ T cells, CD8⁺ T cells and NK cells had elevated expression of multiple immune checkpoint receptors (TIGIT was elevated in all three cellular compartments), and granulocytes, monocytes and B cells, plasma cells and dendritic cells expressed the ligands (Extended Data Fig. 3C). We also detected expression of immune checkpoint ligands in other non-immune cell types, which included fibroblasts, endocrine, and endothelial cells (Fig. 2C). Our single cell data revealed a complex, patient-specific landscape of immune checkpoint ligand and receptor expression across multiple immune and non-immune cell types.

Tumor-infiltrating CD8⁺ T cells have a distinct gene expression profile, with progressive dysfunction in advanced disease.

Cytotoxic T cells are a fundamental component of anti-tumor immune responses and the target of immunotherapy (for review see [30]). To gather deeper insight into the functional status of tumor-infiltrating CD8⁺ T cells, we investigated their transcriptional profile. To investigate patient-specific variability, we mapped the average expression of immune checkpoint receptors in CD8⁺ T cells in each individual patient’s tumor and blood samples (Fig. 2D and Extended Data Fig. 3D). Infiltrating CD8⁺ T cells expressed markedly distinct immune checkpoint profiles in individual patient samples, both in tumors and blood (Fig. 2D and Extended Data Fig. 3D). In tumor samples, LAG3 was elevated in patients 1141, 1294, 1261 and 1229. Patient 1229 (locally advanced) also had high expression of ICOS, CTLA4, TIGIT and CD47. Conversely, patient 1261 (also locally advanced), had elevated CD27, LAG3, PDCD1, HAVCR2, TNFRSF18, CSF1, TIGIT, CD40LG, CD47 and CD28, but not CTLA4. The immune checkpoint landscape did not cluster by disease stage, and while some metastatic patients had high expression of multiple immune checkpoints (3210) others expressed only a limited subset (1253). Analysis of circulating CD8⁺ T cells revealed a similarly complex landscape, but no clear overall correlation in the expression of individual checkpoints between patient tumor and blood T cells at a gene expression level (Extended Data Fig. 3D).

We then investigated tumor-infiltrating CD8⁺ T cells compared to CD8⁺ T cells in normal/adjacent tissue. By functional annotation, we observed pathways relating to cell cytotoxicity, chemokine signaling, T cell receptor signaling, and antigen processing were significantly enriched in PDA samples, compared to adjacent/normal tissue CD8⁺ T cells, an indication that immune responses had been elicited in the tumors (Fig. 2E). We performed unbiased clustering of the CD8⁺ T cell gene expression signatures in patients. Circulating CD8⁺ T cells gene expression patterns did not clearly segregate by disease stage (Extended Data Fig. 3E). We then performed an unbiased differential expression analysis on the tissue-infiltrating CD8⁺ T cells in tumors versus adjacent/normal tissue. This analysis revealed distinct expression patterns in healthy versus tumor infiltrating total CD8⁺ T cells (Fig 2F). We noted that the CD8⁺ T cells in adjacent/normal samples clustered together, while tumor-infiltrating CD8⁺ T cell signatures spanned a spectrum in individual samples. Some tumor signatures partially resembled non-malignant tissue and others were greatly diverging. Interestingly, upon clinical annotation, we discovered that the divergence from the normal signature was more pronounced in samples from advanced disease stage. When we considered which genes were differentially expressed across the groups, we observed an increase in T cell activation and trafficking markers (GZMB, GZMA, KLF2) [3, 4, 31] in tumor CD8⁺ T cells, compared to adjacent/normal tissues. Further, T cell exhaustion markers such as EOMES and GZMK were low in healthy CD8⁺ T cells, and elevated in the majority of tumor CD8⁺ T cell samples [32, 33]. The only immune checkpoint receptor identified as differentially overexpressed in tumor-infiltrating CD8⁺ T cells, compared to CD8⁺ T cells in adjacent/normal tissue, was TIGIT, a gene encoding a receptor belonging to the Ig superfamily [14]. Overall, comparison of the transcriptional profile of tumor and adjacent/normal infiltrating CD8⁺ T cells revealed unique profiles of expression of immune checkpoint genes in individual patients. However, some common features emerged, such as expression of activation markers, as well as an exhausted gene expression signature which progressively increased with advanced disease stage.

Tumor infiltrating CD8⁺ T cells include an expanded exhausted population characterized by TIGIT expression.

Distinct populations of tumor infiltrating CD8⁺ T cells have been described [34]. Given the progressive dysfunction of tumor-infiltrating CD8⁺ T cells, we hypothesized that the transcriptional profile shift might be caused by changes in CD8⁺ populations. By unbiased clustering, we distinguished 6 populations of CD8A-expressing T cells in both adjacent/normal and PDA samples (Fig. 3A). To identify sub-populations, we plotted the top expressed genes per cluster (Fig. 3B), and compared them with published signatures of CD8⁺ T cells subtypes [32]. We identified two populations of effector CD8⁺ T cells (T_eff), expressing PRF1 and GZMB; a population of likely memory (mem)/ precursor effector (pec) CD8⁺ T cells (T_mem/pec) expressing CCR6 [35]; and two populations of exhausted CD8⁺ T cells (T_ex) expressing EOMES, GZMK, and TIGIT (Fig. 3C, Extended Data Fig. 4A and average expression heatmap in Extended Data Fig. 4C). Interestingly, comparison of tumor infiltrating versus adjacent/normal CD8⁺ T cells revealed a relative increase in exhausted T cells and memory T cells, with converse reduction of effector T cell levels (Fig. 3D and Extended Data Fig. 4B). We then examined expression of immune checkpoints and immune activation markers across CD8⁺ T cell subsets, and observed uniform expression of PDCD1, HAVCR2 and LAG3 (Fig. 3C), while TIGIT was enriched in exhausted and memory CD8⁺ T cells. We found that 13.3% of effector CD8⁺ T cells, compared to 44.2% of exhausted CD8⁺ T cells expressed TIGIT. To compare gene expression changes within distinct clusters of CD8⁺ T cells, we performed two separate differential expression analyses. We first compared PDA effector CD8⁺ T cells to adjacent/normal effector CD8⁺ T cells, and we found that GZMA and GZMB were higher in tumor infiltrating effector CD8⁺ T cells, suggesting T cell activation (Fig. 3E). RORA expression, a marker associated with effector T cells, was also upregulated [32] (Fig. 3E), consistent with an ongoing immune response. We then compared PDA exhausted CD8⁺ T cells to adjacent/normal exhausted CD8⁺ T cells and tumor infiltrating exhausted CD8⁺ T cells had higher expression of EOMES and KLF2, markers of exhaustion [32] (Fig. 3F).

Open in a separate window

Figure 3.

Single cell RNA sequencing reveals exhausted CD8⁺ T cell phenotype in PDA patients is defined by the immune checkpoint TIGIT.

(A) UMAP analysis of CD8⁺ T cells from n=3 adjacent/normal pancreas samples (left) and n=16 PDA tumors (right). The 6 identified subsets of CD8⁺ T cells were collapsed into potential memory (blue), effector (pink) and exhausted (green). (B) Single cell resolution heatmap analysis of top 10 genes for each identified CD8⁺ T cell subset. (C) Violin plots of normalized expression for selected markers mapped across the CD8⁺ T cell subsets. (D) Quantitation of potential exhausted (p=9.11E-6), effector (p=2.209E-5) and memory (p=0.0031) T cells in adjacent/normal pancreas and PDA patients, plotted as % total CD8⁺ T cells. Plots represent n=3 adj/norm and n=16 PDA patients. Two-sided Student’s t-test was performed to compare between groups and a p value of 0.05 or less was considered statistically significant. Panel of genes differentially expressed in (E) effector and (F) exhausted CD8⁺ T cells in PDA (red) compared to adjacent/normal pancreas (blue). Plots represent n=3 adj/norm and n=16 PDA patients. Violin plots are shown as normalized expression. All violin plots in (E) and (F) have an adjusted p-value of p<0.01 and are considered statistically significant.

Our data suggest that exhausted CD8⁺ T cells are abundant in pancreatic tumors, and that their exhausted phenotype is more profound than the equivalent population in adjacent/normal tissue. Further, TIGIT was the sole immune checkpoint receptor that specifically defined exhausted CD8⁺ T cells.

A complex landscape of NK and CD4⁺ T cells cell subsets in pancreatic cancer.

Similar to effector CD8⁺ T cells, NK cells display cytotoxic activity and express immune checkpoint receptor (Fig. 2C), however these cells are not well defined in human PDA [36]. Unsupervised sub-clustering of NK cells revealed three populations (Fig. 4A and ​and4B).4B). Along with the three NK cell subsets, we found two additional populations that we labeled cluster 1 and 2, which expressed NKG7, CD3E, and CD8A. While these cells did not cluster with the CD8 T cells based on their transcriptional profile, they could potentially be highly cytolytic CD8 T cells expressing some NK markers [37, 38]. A “NK” cell population expressing CD8A was identified in a recent PDA scRNA-seq paper [12]. Highly variable gene expression analysis highlighted differences among subpopulations (Fig. 4B). In NK subsets we detected expression of markers of antigen presentation (HLA-DRA), cytolytic activity (PRF1, GZMB), and chemokines/chemokine receptors (CCL3, IL7R), in agreement with a recently published characterization of human NK cells by single cell RNA sequencing [39]. NK cluster 1 was enriched for immune checkpoint HAVCR2, while NK cluster 3 expressed high levels of immune checkpoint TNFRSF4 (Fig. 4B and ​and4C).4C). We then performed differential gene expression analysis and unbiased clustering of individual patient samples to compare tumor infiltrating NK cells with NK cells in non-tumor tissue. The signatures of tumor-infiltrating and non-tumor NK cells were not as divergent as was the case for CD8⁺ T cells, and, with one exception (1324), resectable tumor samples clustered closely with the non-tumor samples (Fig. 4D). However, advanced disease samples had a different expression signature than healthy counterparts, with increased expression of activation markers such as GZMA and elevated expression of two immune checkpoint genes, TIGIT and HAVCR2. The role of NK cells in PDA is not well understood; our findings set the stage for future functional studies on the role of NK cells in this disease.

Open in a separate window

Figure 4.

Single cell RNA sequencing of pancreatic tissues reveals TIGIT is differentially expressed in NK cells from PDA patients and is a defining marker of Tregs.

(A) Merged UMAP of 5 identified subsets of NK cells from adjacent/normal pancreas (left) and PDA (right). Plots represent n=3 adj/norm and n=16 PDA patients. (B) Single cell resolution heatmap of each NK cell subset identified. Immune checkpoints (HAVCR2, TNFRSF4) are bolded. (C) Violin plots of normalized average expression within NK cell subsets demonstrating specific lineage markers for NK cells (such as NCAM1/FCGR3A) and immune checkpoint receptors. (D) Unbiased differential average expression of merged NK cells from adjacent/normal pancreas (black) and PDA (grey). Disease stage is plotted on the left. (E) Merged UMAP of all CD4⁺ T cells with 13 identified cell subsets. Naïve CD4⁺ T cells are denoted as T_h0 (CCR7⁺) and Tregs as Tregs (FOXP3⁺). All other subsets are denoted as CD4 T cells. (F) Single cell resolution heatmap of each CD4⁺ T cell subset. Boxes on the left designate naïve CD4⁺ T cells (T_h0) and the CD4⁺ T cell subsets that are defined by immune checkpoint expression (TIGIT, TNFRSF18, PDCD1). (G) Feature plots of CTLA4 and TIGIT in regulatory CD4⁺ T cells (outlined). In all panels, plots represent n=3 adj/norm and n=16 PDA patients.

We then investigated CD4⁺ T cells, a complex population that includes regulatory T cells, and plays a fundamental role in regulating pancreatic carcinogenesis [40–42]. Unlike the CD8⁺ T and NK cells, we could not perform differential expression on CD4⁺ T cells as we did not capture enough cells from adjacent/normal tissue samples; instead we focused on studying the tumor-infiltrating component. We identified 13 transcriptionally distinct populations of CD4⁺ T cells, although many of the different clusters tended to merge together (Fig. 4E). A highly enriched gene analysis of CD4⁺ subsets revealed expression of naïve T cell markers CCR7 and SELL [43] in cluster 0 and expression of regulatory T cell marker FOXP3 in cluster 3 (Fig. 4F). As expected, based on previous studies, CTLA4 was highly expressed within Tregs [44] (Fig. 4G). Specific immune checkpoint genes were highly expressed in individual clusters, such as TNFRSF18 and PDCD1. TIGIT appeared as a top expressed gene in Tregs (Fig. 4F), although it was also expressed in other clusters more sparsely (Fig. 4G). The other clusters did not correspond to known subsets of CD4⁺ T cells, and were relatively similar to one another, a possible reflection of low transcriptional activity of non-Treg CD4⁺ T cells. Taken together, these analyses suggest multiple immune checkpoint receptors are expressed in CD4⁺ and NK cell subsets. In particular, the immune checkpoint TIGIT was differentially overexpressed on tumor-infiltrating PDA CD8⁺ T cells, regulatory T cells and NK cells.

Myeloid and dendritic cells are an important source of immune checkpoint ligands in human PDA.

We next analyzed myeloid cells, which are an important source of immune checkpoint ligands in PDA. By unbiased clustering, we identified 6 transcriptionally distinct populations of myeloid cells. (Fig. 5A). We observed an abundant granulocyte population expressing CXCR1 and CXCR2, FCGR3B and S100A8 (Fig. 5B). Consistent with previous studies [12, 25], we detected resident macrophages and alternatively activated macrophages (MARCO⁺), and classical monocytes (Fig. 5A and ​andB).B). We also observed an additional myeloid population, denoted as alternatively activated macrophages 2, that resembled alternatively activated macrophages and was uniquely defined by abundant expression of CHIT1 and multiple immune checkpoint ligands (Fig. 5A, ​,5B5B and Extended Data Fig. 5A). We next mapped immune checkpoint ligand expression within specific myeloid compartments and observed heterogeneous expression of immune checkpoints in specific clusters (Fig. 5C). Differential expression analysis between adjacent/normal and tumor-infiltrating myeloid sub-clusters revealed multiple upregulated checkpoint ligands. LGALS9, the ligand for HAVCR2 (encoding for TIM3), was significantly increased within alternatively activated macrophages, while SIRPA, the ligand for CD47, was higher in PDA granulocytes compared to granulocytes in adjacent/normal tissue (Fig. 5D). PVR, the ligand for TIGIT was enriched in total macrophages (Fig. 5D). Average expression heatmaps of macrophages (Extended Data Fig. 5B) and granulocytes (Extended Data Fig. 5C) demonstrated that the expression of immune checkpoint ligands was highly variable in individual patients.

Open in a separate window

Figure 5.

Single cell RNA sequencing reveals distinct myeloid and dendritic cell subsets.

(A) Merged UMAP of 6 identified myeloid cell subsets in adjacent/normal pancreas (left) and PDA (right). (B) Single cell resolution heatmap of each myeloid cell subset identified. Boxes on the left designate the top expressing genes for each myeloid subset. (C) Selected feature plots of the immune checkpoints, LGALS9, CD274, PVR, CSF1R, SIRPA, HLA-DQA1 in myeloid cells. (D) Selected feature plots of markers that define alternatively activated macrophages, granulocytes, and total macrophage subsets (left) and violin plots of immune checkpoint ligands that are upregulated in PDA patients (right). (E) UMAP analysis of dendritic cells in merged normal/adjacent pancreas and PDA. (F) Top ten highly enriched gene signature analysis of dendritic cell subclusters identifying potential DC subsets, including plasmocytoid DCs (pDCs), Langerhans-like DCs (Lang_DCs), conventional DCs (cDCs), and activated DCs (Act_DCs). In all panels, plots represent n=3 adj/norm and n=16 PDA patients.

Dendritic cells (DCs) cells are professional antigen presenting myeloid cells, and support anti-tumor activity by stimulating T cells; their relative rarity in PDA is one of the potential causes for ineffective immune responses in this disease [45]. Clustering analysis revealed multiple populations of tumor-infiltrating DCs (Fig. 5E). We found two populations of plasmacytoid DCs (IRF8, GZMB) and two populations of Langerhans-like DCs (CD207, CD1A) as previously described [12] (Fig. 5F and ​and6A).6A). Langerhans-like DC2 had robust expression of IL22RA2, also known as IL22BP, and the IL-22-IL22BP axis is known to be a crucial mediator of tumorigenesis in the colon [46] and pancreas [47]. We also detected a population of conventional DC1s (CLEC9A, IRF8) [12, 48], and two additional populations of potential conventional DC2s that expressed immune checkpoint ligand SIRPA (Fig. 5F and ​and6A).6A). We also detected two unique populations of activated DCs (LAMP3, CCL22) [12] (Fig. 5F and ​and6A).6A). We then plotted the average expression of known immune checkpoint ligands in the different DC subsets. We discovered that activated DC1 had elevated expression of nearly all the immune checkpoint ligands, including PVR (Fig. 6B), suggesting that in pancreatic tumors some subsets of DCs may be immunosuppressive [49].

Open in a separate window

Figure 6.

Predicted ligand receptor mapping in PDA patients demonstrate myeloid and non-immune cell types as sources of immune checkpoint ligands.

(A) Violin plots, where each dot represent a single cell, of select dendritic cell lineage markers across all 9 identified subsets. (B) Immune checkpoint ligand expression heatmap within dendritic subclusters. (C) Circos plot map of all putative ligand receptor interactions that are upregulated in PDA macrophages, (D) granulocytes, (E) dendritic cells, (F) endothelial cells (G) epithelial cells compared to adjacent/normal pancreas visualized by circos plot using the Circos software V0.69-9 (circos.ca). The heatmap within the circos plots is the scaled average expression of each gene within PDA tissue cell populations. The interactions plotted are those in which the expression level of either the ligand, the receptor, or both are increased in expression in PDA samples compared to adjacent/normal tissue. (H) Violin plots for the normalized expression of TIGIT, CD96, and CD226 in CD8⁺ T cells in PDA (red) compared to adjacent/normal pancreas (blue). Between adj/norm and PDA groups, the asterisk indicates P<0.0001, and exact P=4.8E-32. For Figure 6 panels A through H, n=3 adj/norm samples were examined and n=16 PDA patients were analyzed. (I) Dot plot analysis of TIGIT family members within PDAC tissue. Color of dot represents average expression, while the size of the dot represents percent expression. Dot plot represent n=16 PDA patients gene expression.

Mapping predicted interactions and tissue heterogeneity in pancreatic cancer samples by single cell sequencing.

To explore potential cross-talk between T/NK and myeloid populations, we applied a predicted interaction algorithm [15] based on known ligand-receptor (LR) pairs interacting with high affinity [50]. We curated the list to specifically add immune checkpoints and limit the receptor-ligand pairs to cytokines, chemokines and specific signaling pathways [for a comprehensive list, see Supplemental Table 5]. We first plotted all the receptor-ligand interactions that were statistically higher in tumor versus non-malignant samples, based on the level of ligand expression, and observed a complex landscape of potential interactions involving multiple cell types (Extended Data Fig. 5D). We then visualized upregulated ligands in macrophages (Fig. 6C), granulocytes (Fig. 6D), dendritic cells (Fig. 6E), endothelial cells (Fig. 6F) and epithelial cells (Fig. 6G) and mapped the predicted binding partners in CD4, CD8, and NK cells. Among interactions upregulated in cancer compared to adjacent/normal, we detected known putative immune suppressive interactions, such as those mediated by the chemokine receptors CXCR2 in granulocytes and CCR2 in macrophages [25, 26]. Predicted interactions mediated by IL1A and IL1B with their receptor encoding genes IL1R1 and IL1R2 were also upregulated, consistent with their known roles in pancreatic cancer [12, 51, 52]. Multiple putative interactions linked T and NK cells to myeloid immune checkpoint ligands, which is consistent with a key role for myeloid cells in establishing immune suppression in pancreatic cancer (for review see [8]). Predicted immune checkpoint-mediated interactions such as ICOS/ICOSLG and SIRPA/CD47 were among those upregulated in pancreatic cancer compared to healthy/adjacent tissue. TIGIT/PVR interactions were elevated between macrophages and CD4⁺ T cells, CD8⁺ T cells and NK cells (Fig. 6C). Interestingly, the putative TIGIT/PVR interaction was also elevated between tumor endothelial and epithelial cells, T cell and NK cell subsets (Fig. 6F, ​,G).G). We then endeavored to investigate the expression of other genes involved in the TIGIT pathway. We investigated the expression of TIGIT’s costimulatory counter receptor, DNAM1 (CD226), which competes for PVR and PVRL2 and promotes T cell activation [53]. We found that while TIGIT was significantly increased on PDA CD8⁺ T cells (P= 4.8E-32), CD226 expression was not altered in CD8⁺ T cells between adjacent/normal and PDA cells (Fig. 6H). CD96 and PVRIG act similarly to TIGIT, inhibiting T cell activation. Expression analysis showed that mRNA levels of these receptors were not altered between adjacent/normal and PDA CD8⁺ T cells (Fig. 6H). PVRL2 encodes a second ligand for TIGIT, although it binds with a lower affinity compared to PVR [54]. We detected expression of PVRL2 in epithelial, myeloid, and endothelial cells (Fig. 6I). TIGIT, CD96, PVRIG, and CD226 were mainly expressed by T and NK cells in PDA tissue (Fig. 6I). We then investigated the Adenosine pathway, because of its immune suppressive role [55]. We profiled this pathway in tumor samples and found expression of the adenosine receptor ADORA1 in epithelial and mast cells, ADORA2B in epithelial, mast cells, and dendritic cells, ADORA3 in dendritic, mast, and myeloid cells (Extended Data Fig. 5E).

The expression of multiple immune checkpoints has been previously described in pancreatic cancer [4]. Since a deconvolution approach was used, the specific immune cell types expressing receptors and ligands could not be assessed. In contrast, our analysis provides a comprehensive view of the multiple, redundant potential immune suppressive interactions within the pancreatic cancer microenvironment.

Multiplex fluorescent immunohistochemistry (mfIHC) imaging, cell segmentation, and basic phenotyping

Images were taken using the Mantra^™ Quantitative Pathology Work Station (Akoya Biosciences) as described in the Online Methods.

Cytometry Time-of-Flight (CyTOF) Immune Phenotyping and Data Analysis

For CyTOF, we collected 8 samples from non-malignant pancreas specimens, including non-involved pancreas tissue adjacent to a duodenal adenoma (1196, (patient ID)), ampullary carcinoma (1258), insulinoma (19–700), and non-involved pancreas tissue adjacent to PDA (1172, 19–262, 19–561, 19–732, 1252), all obtained surgically (Fig. S1A). PDA samples were collected from either surgical (n=7) or fine needle biopsy (FNB) (n=3) procedures, and clinical annotation is shown in Fig. S1A. Human patient tissues from FNB or surgery were immediately placed into DMEM media supplemented with Y27632 (Rho-Kinase inhibitor) for transport to the laboratory. Sample Preparation, Analysis and Data analysis are described in detail in the Online Methods.

Immunofluorescent staining

Patient tissue slides were rehydrated in xylene, 100% ethanol, 95% ethanol, then running deionized water sequentially. Antigen retrieval was performed with sodium citrate, pH 6.0. Tissue was blocked in 10% donkey serum overnight at 4°C. Primary antibodies (PVR/CD155 (1:100, Cell Signaling Technology), VE-Cadherin/CD144 (1:250, R&D), FOXP3 (1:100, Cell Signaling), Vimentin (1:100, Cell Signaling), CD163 (1:100, Novus Biologicals), or Pan Cytokeratin-488 (1:250, Thermo Fisher Scientific) were diluted in 5% donkey serum in PBST (DS/PBST) and incubated overnight at 4°C. For tissue co-stained for TIGIT-FOXP3 and PVR-Vimentin: the tyramide signal amplification kit with Alexa Fluor 488 (Thermo Fisher Scientific) was used following the manufacturer’s recommendation to enhance signaling for PVR and FOXP3. Samples underwent a second citrate antigen retrieval and were then multiplexed with TIGIT and Vimentin following the aforementioned standard IFC protocol. A 1% BSA block was used throughout the TSA protocol and subsequent multiplex staining. Tissue was mounted with DAPI ProLong^™ Gold Antifade Mountant (Thermo Fisher Scientific) and subsequently imaged by confocal microscopy on a Leica SP5. Supplementary Table 4 lists the antibodies used.

Statistical Analysis and Reproducibility

Significance was evaluated by the following statistical analyses: two-tailed, parametric, unpaired Student’s t-test, Student’s t-test with Welch’s correction, Wilcoxon rank-sum test, or a Mann–Whitney U-test in GraphPad Prism (version 7) or JMP Pro software (version 14). The data were presented as means ± standard error (SEM) or means ± standard deviation (STDV). A p value of p<0.05 was considered statistically significant. Pearson correlation coefficients were used to measure R and R². Intergroup comparisons (differential expression) of scRNA seq was performed using Wilcoxon ranked test and p-values were adjusted for multiple comparisons with the Bonferroni correction method. For the interactome analysis, differences of the ligand and receptors between groups were determined using Wilcoxon ranked test, and p-values were adjusted for multiple comparisons with the Bonferroni correction method. Ligand/Receptor pairs were considered significantly different if the p < 1.0×10⁻⁴. Ligand/Receptor pairs were then sorted by the adjusted ligand expression p-value. No statistical method was used to predetermine sample sizes, experiments were not randomized and mass cytometric analysis of samples was not blinded. No data were excluded from the analyses. Each patient is considered an independent biological sample in the analyses. For comparison of differential abundance analysis of mass cytometry data, the edegR package (version 3.11) was used.

Additional Experimental procedures are included in the Online Methods.

Data Availability

All raw data are publicly available without restrictions. All mass cytometry data used for this publication have been deposited in the FlowRepository. All fcs files of tissue (tumor and adjacent normal) have been uploaded to FlowRepository Experiment FR-FCM-Z2S4 and PBMC files have been uploaded to FlowRepository Experiment FR-FCM-Z2S3. The corresponding file key for the FCS files is in supplementary table 6 and refers to Figures 1 and ​and7,7, and Extended Data 6 and ​and7.7. Single cell RNA sequencing data with clinical metadata are available at NIH dbGAP database under the accession phs002071.v1.p1. Deidentified single cell RNA sequencing data are available at NIH GEO database under the accession GSE155698. Source data are available for this study. All other data supporting the findings of this study are available from the corresponding author on reasonable request.

Cytometry Time-of-Flight (CyTOF) Immune Phenotyping

Tissues were mechanically minced and enzymatically digested with collagenase P (1mg/mL DMEM) at 37 degrees Celsius with gentle shaking and subsequently filtered through a 40μm mesh to obtain single cells. Whole blood was collected pre-operatively into two 10mL EDTA tubes. EDTA tubes were inverted 10 times before centrifugation at room temperature (RT), 1700 × g for 20 minutes. Serum was removed and using a P1000 tip, the white layer of PBMCs at the interface between serum and RBCs was removed and placed into 15mL falcon tube. PBMCs were washed in 3X volume PBS centrifuged at RT, 300 × g for 15 minutes. Following centrifugation, the supernatant was removed, and 10ml ACK lysis buffer was added to lyse RBCs for 10 minutes at RT. Following this, PBMCs were centrifuged at 300 × g for 5 minutes. PBMC and tissue samples were washed twice with MaxPar® PBS (Fluidigm) prior to Cell-ID^™ Cisplatin (Live/Dead staining). Cell-ID^™ Cisplatin reagent (1.67μM) was incubated with tissue and PBMCs single cell suspensions for 5 minutes at RT. To quench this reaction, 4mL of Cell Staining Buffer (Fluidigm) was added to each sample and samples were centrifuged at 300 g for 5 minutes. The supernatant was removed, and cells were washed with 2 mL of MaxPar® Cell Staining Buffer. Cell fixation was achieved by removing the supernatant, and re-suspending the cell pellet in residual volume, prior to the addition of freshly prepared cell fixation buffer (1.6% Methanol-free Formaldehyde; Thermo Fisher 28906 in MaxPar® PBS) for 10 minutes at RT. After fixation, samples were washed twice with 2mL of MaxPar® Cell Staining Buffer and centrifuged at 300 g for 5 minutes. Samples were re-suspended in 1mL MaxPar® Cell Staining Buffer and stored at 4 degrees Celsius for up to one week prior to staining. Up to 3 million cells per sample were stained with cell surface antibody cocktail (All antibodies were purchased from Fluidigm and used at the following dilutions: (CD3 (1:200); CD19 (1:300); CD15 (1:400); CD163 (1:100); CD64 (1:100); CD16 (1:400); LAMP1 (1:100); CD66b (1:200); CCR2 (1:200); TIGIT (1:100); PD-1 (1:100); PD-L1 (1:100); CD8a (1:200); CD33 (1:200); CD45RO (1:200); CD34 (1:100); CD45RA (1:100); CD206 (1:100); CD25 (1:100); CTLA-4 (1:100); CD68 (1:100); PD-L2 (1:100); HLA-DR (1:400); CD14 (1:100); CD4 (1:100); CD11b (1:200); CD45 (1:200); LAG3 (1:100); CD23 (1:100); CD56 (1:100)). (see Supplemental Table 1 for additional antibody information) in 100μl volume of MaxPar® Cell Staining Buffer for 30 minutes at room temperature. After being washed once in 1mL MaxPar® Cell Staining Buffer cells were re-suspended in 2mL cell intercalation solution (125 nM Cell-ID Intercalator-Ir in MaxPar® Fix and Perm Buffer) and shipped to either the Flow Cytometry core at the University of Rochester Medical Center or the Indiana University Simon Cancer Center Flow Cytometry where CyTOF2 Mass Cytometer cell acquisition was performed.

CyTOF Data Preprocessing

Normalized FCS files were analyzed using the Premium CytoBank Software V7.3.0 (cytobank.org). Data were checked for quality of staining and normalized by the use of internal bead standards. Live singlet cells were identified using the combination of Ir191 DNA Intercalator, Event Length, and Pt195 Cisplatin staining intensity channels. Filtered live single cells were exported as new FCS files for downstream analysis.

CyTOF Analysis

Unbiased identification of cellular subpopulations was performed in parallel using multiple approaches – visualization through FlowSOM-viSNE in R, where an initial FlowSOM clustered cells into 100 initial nodes, followed by the ConsensusClusterPlus package which, along with manual annotation helped to further consolidate the clusters based on cell surface marker expression [18], or Astrolabe Cytometry Platform (Astrolabe Diagnostics, Inc.), where single-cell data was clustered using the FlowSOM R package [61] and labeled using the Ek’Balam algorithm [19]. The hierarchical clustering for all heatmaps uses the Pearson’s correlation as a distance metric. Differential abundance analysis was performed using the edgeR V3.11 R package [62, 63]. We used a combination of manual gating validation and unbiased approaches to analyze our datasets and included samples with >3000 live singlets in clustering algorithms.

Treatment of Batch Effects

In order to avoid batch effects within the data analysis, the Astrolabe Cytometry Platform did not compare numerical intensities between samples [19]. Each sample was analyzed separately, and then comparisons were done using either cell frequencies (such as comparing T Cell counts) or qualitative values (“CD3 high” versus “CD3 low”). The underlying assumption was that a given subset was the same regardless of if underlying marker intensity has shifted; in other words, a T Cell was defined as a T Cell whether the CD3⁺ peak was centered around a transformed intensity of 4, or a transformed intensity of 6. This mirrors the approach utilized in manual gating analysis.

t-SNE Visualization

For the t-SNE maps, each sample was uniformly downsampled into at most 10,000 cells. Samples were then concatenated, and the complete data set was uniformly downsampled into at most 500,000 cells. t-SNE algorithm was run using the Rt-SNE package: https://github.com/jkrijthe/Rt-SN.

Single-cell RNA sequencing

Tissues were mechanically minced and enzymatically digested with collagenase P (1mg/mL DMEM) and subsequently filtered through a 40μm mesh to obtain single cells. Dead cells were removed using MACS^® Dead Cell Removal Kit (Miltenyi Biotec Inc.). Single-cell cDNA libraries were prepared and sequenced at the University of Michigan Sequencing Core using the 10x Genomics Platform. Samples were run using paired end 50 cycle reads on HiSeq 4000 or the NovaSeq 6000 (Illumina) to a depth of 100,000 reads. The raw data were processed and aligned by the University of Michigan DNA Sequencing Core. Cellranger count version 3.0.0 with default settings was used, with an initial expected cell count of 10,000. In all cases the hg19 reference supplied with the cellranger software was used for alignment. R Studio V3.5.1 and R package Seurat version 3.0 was used for single cell RNA-seq data analysis similarly as previous described. Data were initially filtered to only include all cells with at least 200 genes and all genes in greater than 3 cells. Data were initially normalized using the NormalizeData function with a scale factor of 10,000 and the LogNormalize normalization method. Variable genes were identified using the FindVariableFeatures function. Data were assigned a cell cycle score using the CellCycleScoring function and a cell cycle difference was calculated by subtracting the S phase score from the G2M score. Data were scaled and centered using linear regression on the counts and the cell cycle score difference. PCA was run with the RunPCA function using the previously defined variable genes. Violin plots were then used to filter data according to user-defined criteria. All tissue samples were batch corrected through the R package Harmony V1.0 (https://github.com/immunogenomics/harmony). Harmony is a flexible multi-dataset integration algorithm for scRNA-seq by correcting the low-dimensional embedding of cells from principal component analysis (PCA). It first uses soft clustering to find potential clusters, and then uses a soft k-means clustering algorithm to find clusters that favors the cells from multiple datasets and penalizes for any specified unwanted technical or biological factors. It then learns a simple linear adjustment function by computing cluster-specific linear correction factors, such as individual cell-types and cell state, from the cluster-specific centroids from each dataset. Each cell is weighted and corrected by its cell-specific linear factor. It then iterates the clustering and correction until the cell cluster assignments are stable. We used Harmony V1.0 to integrate our scRNA-seq patient data, correcting for individual scRNA-seq Run IDs (as each individual patient was each their own Run ID). Cell clusters were identified via the FindNeighbors and FindClusters function using a resolution of 1.2–2 for all samples and Uniform Manifold Approximation and Projection (UMAP) clustering algorithms were performed. FindAllMarkers table was created and clusters were defined by user-defined criteria. Code is publicly available on GitHub.com (https://github.com/PascaDiMagliano-Lab/MultimodalMappingPDA-scRNASeq).

We thank Matthew Cochran and Terry Wightman at the Flow Cytometry core at the University of Rochester Medical Center and Andrea Michelle Gunawan at the Indiana University Simon Cancer Center Flow Cytometry for their support in cell CyTOF acquisition. We thank Vinicius Motta and Kevin Brown from Fluidigm for their assistance with panel design. We would like to thank Patricia Schnepp and Aquila Ahmed for their assistance with CyTOF experimental design. We would also like to thank Tricia Tamsen and Judy Opp from the University of Michigan Advanced Genomics Core. We would like to thank David Hill and Michael Czerwinski for their input on designing single cell analysis pipelines. We would very much like to thank Amir Gilado and Ido Amit for their expertise in building our pancreatic interactome network. We would like to thank the Tissue Procurement Center at the University of Michigan. Thanks to Ed Stack formerly with Perkin Elmer for assistance with initial R introduction and basic training using inForm 2.3.0 and earlier versions. and staining strategies. We thank Philip Turncliff for the excellent graphics. We also thank Jason Spence for the VE-cadherin antibody gift.