Integrated Genomic and Immunophenotypic Classification of Pancreatic Cancer Reveals Three Distinct Subtypes with Prognostic/Predictive Significance

Pancreatic ductal adenocarcinoma (PDAC) is considered a highly immunosuppressive and heterogeneous neoplasm (1, 2). Despite improved knowledge regarding the genetic background and an increasing understanding of the tumor microenvironment, targeted therapies and immunotherapy, although efficient for many solid malignancies, including metastatic melanoma and lung cancer, have not yielded any clinical benefit in PDAC (2, 3). Moreover, there are currently no approved treatments that target driver mutations in PDAC, such as KRAS, TP53, CDKN2A, and SMAD4 (3). Ongoing studies aiming to target the stroma and the immune microenvironment either alone or in combination are expected to lead to more endurable treatment outcomes (4–6). Specifically, inhibition of the bromodomain and extraterminal (BET) family of epigenetic readers was shown to block stroma-inducible transcriptional regulation in vitro and tumor progression in vivo (4). Moreover, expression of antitumor immunity genes to uncover biomarkers predictive of response to systemic therapies has been undertaken (5). Furthermore, mutations in the oncogenes PIK3CA, FGFR3, and RAS/RAF family members, as well as the tumor suppressor TP53, were linked to changes in immune infiltration, supporting the relevance of tumor genetics to questions of efficacy and resistance in checkpoint blockade therapies (6). A detailed characterization of the tumor microenvironment combining tumor- and host-associated factors may lead to the identification of novel biomarkers and pathways for a more targeted approach to the use of immunotherapy and/or other combinatorial treatments to overcome the immunosuppressive mechanisms in the PDAC microenvironment.

The immune microenvironment is also of importance for the local aggressiveness of tumor cells and PDAC progression, as specific immune expression signatures may render the tumor microenvironment permissible for single cancer cell invasion (7). Aggressive PDACs are characterized by increased numbers of dissociative growing tumor cells at the invasive front with epithelial–mesenchymal transition (EMT)–like features, coined “tumor buds,” that have been shown to represent an independent adverse prognostic factor in many gastrointestinal cancers, including PDAC (8–14). The microenvironment surrounding the tumor buds is therefore especially interesting, as it likely promotes their migration, angioinvasion, and metastatic potential. We recently reported that PDACs with EMT-like tumor budding display a tumor-favoring immune microenvironment that supports further tumor growth (15).

More recently, an integrated genomic analysis of PDACs identified four molecular subtypes—squamous, pancreatic progenitor, immunogenic, and aberrantly differentiated endocrine exocrine—indicating differences in the molecular evolution of PDAC subtypes and identifying opportunities for therapeutic development (16–18).

We hypothesized that the local immune phenotype will be of importance for more than the immunogenic “gene program.” We therefore analyzed the immune response patterns in the tumor microenvironment of PDAC from a well-characterized cohort and correlated the immune pattern to tumor-related factors, such as the mutational background and EMT-like tumor budding as well as clinicopathologic features, including survival of the patients, to classify PDAC into clinically/biologically relevant subgroups. We also attempted to accommodate these subgroups to the abovementioned molecular subtypes of PDAC according to Bailey and colleagues (16). Our specific aims were to identify valuable clues for the different immunosuppressive mechanisms of PDACs, outline a morphologic substrate for the recognition of the molecular subtypes in daily clinical practice, and provide a basis for a more successful and individualized therapeutic approach.

Histomorphologic data from 110 PDACs, stage I to III, surgically resected between 2002 and 2011, were reviewed from the corresponding hematoxylin and eosin (H&E)–stained slides while clinical data were obtained from corresponding reports (database). Clinicopathologic information for all patients included age; gender; tumor diameter; number of positive lymph nodes; total number of lymph nodes harvested; pathologic tumor–node–metastasis (pTNM) stage; perineural, blood vessel, and lymphatic invasion; and resection margin status. Preoperative CA19-9 values were available from 70 patients (Supplementary Table S1). Tumors treated neoadjuvantly, stage IV tumors, as well as carcinomas with histology other than ductal adenocarcinoma (acinar cell carcinomas, adenosquamous carcinomas, neuroendocrine carcinomas, mixed subtypes, and mucinous neoplasms with associated invasive carcinoma) were excluded from the study. Staging was performed using the AJCC Cancer Staging Manual 8th edition (19). This study was conducted in accordance with recognized ethical guidelines (Declaration of Helsinki), including written informed consent from the patients or their living relatives, and the use of the material was approved by the local ethics committee of the University of Bern (Bern, Switzerland; KEK Nr 200/14).

Tissue microarrays were constructed using the next-generation tissue microarray (ngTMA) approach (20). For each patient, one H&E-stained whole tissue slide containing representative regions of tumor center and invasion front was scanned (Pannoramic P250; 3DHISTECH). Using a tissue microarray annotation tool of 0.6 mm in diameter, slides were digitally annotated as follows: four regions of tumor center (green) and four regions of invasion front (red) including regions of highest tumor budding if available. Next, corresponding formalin-fixed (10% buffered formalin) paraffin-embedded tissue blocks were loaded into an automated tissue microarrayer (TMA Grand Master; 3DHISTECH). The digital slides were aligned with the corresponding donor block. Annotated regions were cored from the donor block and transferred to the recipient ngTMA. Three different ngTMA blocks were produced, resulting in 960 spots for evaluation.

Tumor budding has been evaluated for the purposes of another study according to the International Tumor Budding Consensus Conference (ITBCC) method (21) and was defined according to ITBCC as single tumor cells or tumor cell clusters of up to four cells. Whole tissue sections of the PDACs, stained with H&E as in the routine diagnostics, were used. Two experienced pathologists (E. Karamitopoulou and M. Wartenberg) independently searched all tumor slides throughout at low magnification. Densest budding area (hot spot) was selected by visual estimation regardless of the topographic area within the tumor (intratumoral or at the invasive front). Tumor buds in this area were counted at ×20 magnification (field area, 0.785 mm²). Density of tumor buds was assigned into four groups: no budding (BD-0): 0 buds; low budding (BD-1): 1–4 buds; intermediate budding (BD-2): 5–9 buds; and high budding (BD-3): ≥10 buds. For simplicity reasons, categories 0 and 1 (no and low budding) were grouped together, as only very few patients had 0 buds. The interobserver agreement was assessed using the Kappa statistic (k) for categorical counts and the intraclass correlation coefficient (ICC) for the number of buds.

ngTMAs were sectioned at 3 μm, dewaxed, and rehydrated in dH2O and stained immunohistochemically for CD3 (1:400, clone SP7, Abcam ab16669), CD4 (1:100, clone CD4/4B12, Dako M7310), CD8 (1:100, clone C8/144B, Dako M7103), CD20 (1:100, clone L26, Dako M0755), hyaluronan-mediated motility receptor (RHAMMl 1:100, clone 2D6, 1:50; Leica Biosystems), p63 (1:40, clone 7JUL, Novocastra NCL-L-p63), PD-L1 (1:400, clone E1L3N, Cell Signaling Technology 13684), MLH1 (1:200, clone ES05, Novocastra NCL-L-MLH1), PMS2 (1:100, clone A16-4, BD Pharmingen 556415), MSH2 (1:500, clone G219-1129, Cell Marque 286M-15), and MSH6 (1:50, clone PU29, Novocastra NCL-L-MSH6). Antigen retrieval was performed with Tris-HCl pH 9 for 30 minutes at 95°C. Antibody testing and staining protocols have been established, and staining was performed by an automated Leica BOND RX system (Leica BOND RX, Leica Biosystems) with the Bond Polymer Refine Detection Kit (with DAB as chromogen) and Bond Polymer Refine Red Detection Kit for the double staining (Leica Biosystems).

Protein expression was evaluated by estimating the percentage of positive cells per TMA punch. Because each patient had multiple tumor cores taken from different regions within the tumor, the percentage of positive tumor cells across all cores was averaged.

Proportions of PD-L1–positive carcinoma cells were scored according to a six-step system (22). Accordingly, carcinoma cells were considered positive if the cell membrane was partially or completely stained. Cytoplasmic staining was disregarded, and necrotic areas were excluded. The following categories were used: score 0: 0% to 1%, score 1: ≥1%, score 2: ≥5%, score 3: ≥10%, score 4: ≥25%, and score 5: ≥50%.

Scoring of RHAMM IHC was adapted from Rüschoff and colleagues (23). In more detail, score 0 was assigned when no positivity was observed; score 1 was assigned when high-power magnification (×20–×40) was needed to detect expression; score 2 was recorded when staining was observed under medium power (×10); and score 3 was assigned when positive staining was observed at low magnification (×5). In addition, the total number of tumor buds and proportion of buds expressing RHAMM (score ≥1) were recorded for each tumor. Cases with more than 20% RHAMM-positive buds were classified as “high RHAMM–expressing buds,” and cases with less than 20% RHAMM-positive buds were classified as “low RHAMM–expressing buds” based on ROC analysis (24).

For the immune cell counts, both intraepithelial and stromal counts were assessed by visual estimation on a numerical scale. Because a single tissue microarray spot may contain various degrees of stroma versus tumor epithelial content, the percentage of stroma and tumor tissue per spot was recorded. Cell counts were normalized for tumor cell and stromal content of each spot. Evaluation was performed blinded to clinical endpoints. To verify the robustness of the data by visual estimation, stained TMA slides were digitally analyzed using the Visiopharm Integrator System (VIS) for Windows 8.1, version 6.5.0.2303 (Visiopharm A/S). A custom Analysis Protocol Package (APP) was developed with the VIS author module. The thus obtained area measurements and cell counts were exported. To estimate the interrater reliability for the counts of stromal and intraepithelial lymphocytes, and tumor buds between visual estimation and digital image analysis, Pearson ρ and Krippendorff α values were calculated using R, including the "irr" package.

Analysis of FOXP3 and PTEN was performed previously (15, 25).

Library preparation using the Ion AmpliSeq Cancer Hotspot Panel v2 (Thermo Fisher Scientific) was performed following the manufacturer's instructions and as described previously (26–28). The Ion AmpliSeq Cancer Hotspot Panel v2 was designed to allow amplification-based capture and sequencing of coding regions of 50 cancer-related genes (https://tools.thermofisher.com/content/sfs/brochures/Ion-AmpliSeq-Cancer-Hotspot-Panel-Flyer.pdf). This panel includes 207 primer pairs and requires a minimum initial amount of 10 ng of DNA as input. After DNA extraction and quantification, multiplex PCR for target enrichment was performed using genomic DNA mixed with primer pools and the Ion AmpliSeq HiFi Master Mix (Ion AmpliSeq Library Kit 2.0; Thermo Fisher Scientific) for 2 minutes at 99°C, followed by 22 cycles of 99°C for 15 seconds and 60°C for 4 minutes and holding at 10°C. The obtained PCR products were subsequentially treated with 2 μL of FuPa reagent (Ion AmpliSeq Library Kit 2.0; Thermo Fisher Scientific) to partially digest the primer sequences and then phosphorylated at 50°C for 20 minutes, followed by 55°C for 20 minutes, and finally at 60°C for 20 minutes. The generated amplicons were ligated to Ion beads adapters (Ion AmpliSeq Library Kit 2.0; Thermo Fisher Scientific) and sample-specific barcodes from the Ion Xpress Barcode Adapters Kit (Thermo Fisher Scientific) for 30 minutes at 22°C and then 72°C for 10 minutes. Afterward, the generated barcoded libraries were purified using Agencourt AMPure XP reagents (Beckman Coulter). The library concentration was determined using the Ion Universal Library Quantitation Kit (Thermo Fisher Scientific) and the amplicons size controlled with the 2100 Bioanalyzer (Agilent Technologies). A total of 40 pmol/L for each library were loaded into the Ion Chef System using the Ion PGM IC 200 Sequencing Kit (Thermo Fisher Scientific) for fully automated emulsion PCR and chip loading (37 samples were multiplexed on one 530 chip). Finally, loaded chips were sequenced (500 flows) using the Ion S5 Sequencing Kit using the S5XL Sequencer (Thermo Fisher Scientific).

Raw data (FASTQ files) for each sample were processed for the alignment of sequencing reads with the human genome reference (hg19) using the Torrent Suite software v5.2 (Thermo Fisher Scientific). The alignment pipeline also included signaling processing, base calling, quality score assignment, adapter trimming, and control of mapping quality. Local realignment and quality base score recalibration were also carried out using the Genome Analysis Toolkit (GATK). Coverage metrics for each amplicon (minimal acceptable coverage threshold was set at 500×) were obtained by running the Coverage Analysis plugin software v5.2.2 (Thermo Fisher Scientific). Base calling was performed using the Ion Reporter v5.2 (Thermo Fisher Scientific). In addition, MuTect4 and VarScan25 algorithms were used to call somatic single-nucleotide variants and VarScan25, Strelka6, and Scalpel7 for insertions and deletions (indels). We chose these specific mutations caller strategy as they allow identification of somatic mutations down to 5% mutant allele fraction (i.e., 5% of the reads harboring a given mutation). Finally, all candidate mutations were manually reviewed using the Integrative Genomics Viewer8 (29).

Spearman correlation coefficients (r_s) were used to determine the strength of the relationship between immune cell counts and clinicopathologic tumor characteristics. For consistency, scores were then dichotomized according to the mean into low and high groups. The χ² test was used to determine the association of low and high cell counts with categorical variables and t test for tumor size and age. The presence of association between categorical variables or between normally distributed variables, such as age, was calculated with the Wilcoxon–Mann–Whitney test or t test, respectively. Classification trees were used to calculate threshold values and predictive algorithms for survival. Univariate Cox regression analysis was performed to determine the presence of association of continuous variables with survival prior to dichotomization. Kaplan–Meier survival curves were plotted and survival time differences analyzed using the log-rank test. Multivariable Cox regression analysis was performed to determine the independent prognostic effect of the feature after adjustment for potential confounders. HRs and 95% confidence intervals (CI) were used to determine effect size. Pearson correlation coefficient was used to evaluate the strength of the linear relationship between CD68 counts by a human observer and by software. All P values were two sided and considered significant when P < 0.05. Analyses were conducted on SPSS (V21), SAS (V9.3), and Statistical Package Software R [Version 3.4.1 (2017-06-30), www.r-project.org].

The study design is shown in Supplementary Fig. S1. The study was designed to comply with the REporting recommendations for tumor MARKer (REMARK) guidelines for tumor marker prognostic studies (30).

After comparing the budding scores, the ICC for budding density between the two pathologists was 0.85, and the weighted kappa value for categorical assessments was 0.62 (0.5–0.73), suggesting substantial agreement.

We evaluated both intraepithelial and stromal immune cell counts (CD3, CD4, CD8, and CD20) per tumor. Normalized intraepithelial counts were generally very low (0–14.7/0.274 mm²) and showed only few correlations with clinicopathologic features. On the contrary, stromal counts were more abundant (0–347, 5/0.274 mm²) and showed an inverse correlation to many adverse clinicopathologic features (Supplementary Table S2). Representative findings of the immune cell infiltrates with respect to PDAC phenotypes are shown in Figs. 1–3 and Supplementary Fig. S2.

Presence of tertiary lymphoid tissue (TLT; 21/120 patients; 17.5%) was correlated with low tumor budding (P = 0.0011) and increased overall survival (OS, P < 0.0001), as well as disease-free survival (DFS, P = 0.0067). The impact of TLT on survival was found to be independent after adjusting for TNM and budding itself (Supplementary Table S3; Fig. 3).

An increased stromal CD8/FOXP3 [effector T cells (Teff)/regulatory T cells (Treg)] ratio correlated with a lower rate of lymph nodes (P = 0.0391) and DFS (P = 0.0449), as well as better OS in univariate (P = 0.019), but not in multivariate analysis. Overall, patients with low tumor budding (BD-1), presence of TLT, and low FOXP3 stromal counts had the best prognosis (Table 1). Stromal CD8/FOXP3 ratio showed an inverse correlation with tumor budding (Supplementary Fig. S3). Correlation between digitally analyzed TMA slides and visual estimation is shown in Supplementary Fig. S4. Substantial correlations were found between stromal counts by visual estimation and by software (Pearson ρ 0.87). The density of CD3, CD8, and FOXP3 counts for each subtype is shown in Supplementary Fig. S5.

Nuclear p63 protein expression, a marker of a squamous phenotype (31), was focally observed in 18% of the PDACs (Fig. 1), frequently also in the tumor buds. Eight percent of the p63-positive cases belonged to the high budding category (BD-3, median: 26 buds).

Positive RHAMM staining (score ≥1) was observed in 49% of the tumors and high RHAMM–expressing buds were seen in 31% of the RHAMM-positive cases (15% of the cohort; Figs. 1 and 2). All cases with high RHAMM–expressing buds belonged to the high budding category (BD-3; see also PDAC subtypes below).

Overall, 27% of the cohort showed a membranous PD-L1 positivity in the tumor cells, with most of them showing only focal positivity: 40% of the positive cases exhibited positivity in <5% of the tumor cells, and 70% of the cases exhibited positivity in ≤10% of the tumor cells. Cases with scores 4 and 5 (staining in ≥25% of the tumor cells) constituted 30% of the PD-L1–positive cases, and 7% of the whole cohort belonged to the high budding category (BD-3, mean number of buds: 24.4) and had reduced OS (median: 10 months). Moreover, they showed similar immunophenotypical features, being characterized by a high CD8/FOXP3 stromal ratio (median: 19.32) and a generally immune cell–rich microenvironment (Table 1; Supplementary Fig. S2).

Very few cases (4% of the cohort) showed an IHC loss of the DNA mismatch repair proteins (MLH1, PMS2, MSH2, and MSH6) in different combinations [microsatellite instability-high (MSI-H); Supplementary Fig. S6]. The patients in this category were generally older (median: 69 years), and the tumors were characterized by a microenvironment rich in immune cells with a very high stromal CD8/FOXP3 ratio (55.2%).

The genes included in the Ion AmpliSeq Cancer Hotspot Panel v2 are depicted in Supplementary Table S4A. KRAS mutation (mostly in exon 2) was detectable in 92.8% of all cases (allelic frequency between 7% and 71%) and was the sole mutation detected in 15% of the PDACs. The KRAS wild-type cases (7.2%) showed mutations in TP53, CDKN2A, ATM, and/or BRAF. TP53 was the second most common mutation (65.4%, allelic frequency between 6% and 67%), in 98% as co-mutation with KRAS and/or other genes. SMAD4 was mutated in 25% of the PDACs, and CDKN2A mutations were found in 15.3% of the cohort (allelic frequency between 6% and 33%). Mutations in other genes were less frequent (Supplementary Table S4B; Fig. 4A). Seven percent of the cases also showed many subclonal populations with minor allelic frequency (<5%). Twenty-eight percent of the cases showed mutations in genes additional to the four main drivers (KRAS, TP53, SMAD4, and CDKN2A). The number of mutations did not impact survival (P = 0.0819). All detected mutations are depicted in Supplementary Fig. S7 (OncoPrint).

Preoperative carbohydrate antigen 19-9 (CA 19-9) values were available from 70 patients. They ranged between 10 and 6,835 U/mL (median: 270 U/mL) and correlated with OS and DFS (P < 0.0001). Patients with CA19-9 values <100 U/mL (21.4% of the cohort) survived longer, whereas patients with CA19-9 values >800 U/mL (28.6% of the cohort) had the worst prognosis (Supplementary Fig. S8).

On the basis of the immune cell composition of the microenvironment, we were able to define three distinct groups of PDAC with clinicopathologic and tumor cell–specific implications:

• The “immune-escape” (54%) subtype has a microenvironment poor in T (CD3, CD4, and CD8) and B cells (CD20) but enriched in FOXP3⁺ Tregs. These PDACs show morphologically high tumor budding, focal p63 expression (especially in the tumor buds), and adverse clinicopathologic features (Table 1; Fig. 1). This group is molecularly characterized by KRAS mutations in 94%, along with a high TP53, CDKN2A, SMAD4, and PIK3CA mutation rate (Fig. 4A and B). Although the prevalence of MET and ERBB4 mutations was generally low, all detectable mutations in these genes belong to this category. All this correlates with a poor outcome (median OS: 10 months). Moreover, immune-escape cases often showed high preoperative CA19-9 levels (median: 800 U/mL; Supplementary Fig. S8), while the majority of the cases with high RHAMM–expressing buds belonged to this category (69% of the cases with high RHAMM–expressing buds, representing 22% of the immune-escape cases; Fig. 1; Supplementary Fig. S9).

• The “immune-rich” (35%) subtype has a microenvironment rich in T (CD8, CD3, and CD4) and B cells (CD20) and poorer in FOXP3⁺ Tregs. These tumors exhibit infrequent tumor budding and favorable clinicopathologic features (Table 1; Fig. 2). Molecularly, they are characterized (apart from KRAS mutations in 94.5%) by a lower CDKN2A, SMAD4, and PIK3CA mutation rate as compared with the previous subtype. The few GNAS and IDH2 mutations identified in our cohort are present in this category (Fig. 4A and B). All these characteristics contribute to a better outcome (median OS: 19 months). The presence of TLT defined a subpopulation of the immune-rich subtype (44% of the immune-rich cases, 17.5% of the whole cohort) with an additional high prevalence of STK11, ATM, and SMARCB1 mutations, even more favorable features (Supplementary Table S3; Table 1; Fig. 3), and the best outcome (median OS: 23 months). Only 12% of the cases with high RHAMM–expressing buds belonged to this category (5% of the immune-rich cases; Fig. 2; Supplementary Fig. S9). Patients of the immune-rich subtype had the lowest preoperative CA19-9 levels (median: 109 U/mL; Supplementary Fig. S8).

• The “immune-exhausted” (11% of all cases) subtype consists of two subpopulations showing considerable similarities concerning their microenvironmental features: one with high PD-L1 expression (65% of the immune-exhausted cases, 7% of the whole cohort), frequent tumor budding, and a high PIK3CA mutation rate, along with an immunogenic microenvironment with a high CD8/FOXP3 ratio and unfavorable outcome (median OS: 10 months; Table 1; Supplementary Fig. S2), and a small, microsatellite-unstable subpopulation (35% of the immune-exhausted cases, 4% of the whole cohort) with loss of DNA mismatch repair proteins in various combinations, an antitumoral immune microenvironment with the highest CD8/FOXP3 ratio, and a high prevalence of JAK3 in addition to PIK3CA mutations (Supplementary Fig. S5; Fig. 4A). Cases with high RHAMM–expressing buds were also highly represented in the PD-L1 category (19% of the cases with high RHAMM–expressing buds, representing 43% of the PD-L1 cases; Supplementary Fig. S9).

KRAS mutations were present in 94% to 100% of the cases in all subtypes.

The prognostic impact of the three main PDAC subtypes is depicted in Fig. 5. The immune-rich subtype has better prognosis, whereas the immune-escape and immune-exhausted subtypes have worse prognosis (P < 0.001 for OS and P = 0.004 for DFS). Moreover, the various combinations of tumor budding, stromal FOXP3 counts, and CDKN2A mutations seem to significantly affect survival in PDAC. Stromal FOXP3 counts can substratify low-grade budders into two further categories, with low budding/low stromal FOXP3 having the best prognosis, whereas the presence of CDKN2 mutations seems to add independent prognostic information and leads to substratification of the high-grade budders (Supplementary Fig. S10). Thus, PDACs with high tumor budding, presence of CDKN2A mutations, and high stromal FOXP3 counts have the worst prognosis in our cohort (Supplementary Fig. S10; Supplementary Table S5A and S5B).

Here, we show that the immune response patterns in the microenvironment of PDAC, in correlation with the mutational background and the clinicopathologic features, define three distinct subtypes with clinical and prognostic relevance. The immune-exhausted subtype might be helpful in stratifying patients for different immunomodulatory therapies.

One of our most interesting findings are the differences between the immune-escape and the immune-rich subtype, which seem to have reverse microenvironmental response patterns and almost the opposite mutational background. The immune-escape category is characterized by a high CDKN2A, SMAD4, and PIK3CA mutation rate, high-grade tumor budding, and a microenvironment rich in FOXP3⁺ Tregs and poor in T and B cells, whereas the immune-rich subtype (with or without TLT) is characterized by a lower CDKN2A, SMAD4, and PIK3CA mutation rate and a high prevalence of STK11, ATM, and SMARCB1 mutations, along with infrequent tumor budding and a microenvironment rich in T and B cells and poorer in FOXP3⁺ Tregs. This leads to different morphologic and clinicopathologic features and impacts clinical outcome. Furthermore, it suggests that tumor budding cells, which are significantly more frequent in the immune-escape subtype, are embedded in and probably interact with the tumor-permissive microenvironment, while their close interaction with the FOXP3-Tregs opens possibilities for therapeutic intervention (31). Although this interplay with the immune cells guarantees the survival of the prognostically unfavorable tumor buds, there is evidence that other components of the tumor microenvironment, like surrounding stromal cells, also support EMT-like tumor budding by expressing high levels of E-cadherin suppressors and/or by contributing to miRNA dysregulation (32, 33). In addition, increasing evidence has linked EMT to cancer stem cell features, thus supporting the suggestion that EMT-like tumor buds may represent a population of migrating cancer stem cells (34). For example, the WNT pathway, one of the common pathways undergoing genetic alterations in PDAC (35) and strongly associated with the promotion of a stem cell–like phenotype (36), is involved in the development of EMT-like tumor budding (9, 32). Furthermore, it has been shown that EMT-like cells (like tumor budding cells) share similar molecular characteristics with cancer stem cells, for example, they are both drug resistant and have higher metastatic potential (11, 37). This may further contribute to the worse prognosis of the immune-escape PDAC subtype.

The expression pattern as well as the molecular and clinical characteristics of the immune-escape phenotype largely overlap with the squamous subtype described by Bailey and colleagues (16) or the quasi-mesenchymal subtype described by Collisson and colleagues (17), whereas the pattern that characterizes the immune-rich phenotype is more compatible with the pancreatic progenitor subtype described by Bailey and colleagues (16). We could thus assume that the presence of high-grade budding and increased FOXP3 counts can be regarded as the morphologic signature of the squamous subtype, whereas the picture associated with low-grade budding and low stromal FOXP3 can be regarded as the morphologic substrate of the pancreatic progenitor type. Our results are also compatible with Siemers and colleagues (6), who observed a strong link between CDKN2A mutation (as in our immune-escape subtype) and reduced immune cell infiltrates across 40 tumor cohorts. Most importantly, the molecular differences create opportunities for various treatment modalities toward a more individualized therapy approach. For example, MET pathway–targeted anticancer therapies may be more effective in PDACs of the immune-escape subtype, to additionally target EMT-like tumor budding, as mutations in the receptor tyrosine kinases MET and ERBB4 were more frequent in the immune-escape cases (38, 39). In a recent study, dual MET and ERBB inhibition was shown to overcome intratumor plasticity in osimertinib-resistant advanced non–small cell lung cancer (40). Given the limited therapy options for PDAC, administration of such combinatorial treatment can be exploited in PDACs of the immune-escape subtype.

SMAD4 acts as the central mediator of TGFβ signaling, and its inactivation is frequently found in PDACs (3, 41). The role of the TGFβ pathway as a tumor promoter or suppressor is still a matter of debate. It is generally accepted that TGFβ functions as a tumor suppressor in the early phase of tumorigenesis but can be converted to a tumor promoter during cancer progression (42), thereby contributing to a permissive microenvironment for tumor growth and metastasis (43). Our own results show that inactivating SMAD4 mutations are more prevalent in the unfavorable immune-escape subtype, compatible with a tumor suppressor activity.

The CD20⁺ and CD3⁺ stromal infiltrates can give rise to two different, sometimes interrelated patterns, including the loose infiltrates and the formation of TLT (44). Presence of TLT (a subtype of the immune-rich phenotype) conveyed a strong antitumor impact, resulting in favorable features, with very low tumor budding and a strong and independent survival advantage in our cohort. This is in agreement with recently published data associating rich B-lymphocyte infiltrates and especially TLT with better prognosis in PDAC (44). Moreover, PDACs of the “immune rich with TLT” subpopulation frequently exhibit mutations in DNA damage repair genes (ATM and STK11), potentially leading to an increased number of tumor antigens and/or to defective DNA repair that may have direct therapeutic implications (45). Mutations in STK11 have been found to promote activation of the PI3K/AKT/mTOR pathway, pointing to new drug targets (46).

A further interesting finding is the identification of an immune-exhausted subtype with two subpopulations, showing similarities to the immunogenic subtype by Bailey and colleagues (16). The one characterized by high PD-L1 expression exhibits unusual characteristics combining unfavorable clinicopathologic features like frequent tumor budding, along with an immune-rich microenvironment and increased Teff/Treg ratio. Molecularly, these tumors are characterized by a higher frequency of PIK3CA mutations as well as a higher rate of PTEN loss (80%) compared with the other subtypes. Consistent with this, PIK3CA mutations have been associated with higher estimates of CD8⁺ T cells and natural killer (NK) cells, along with decreased Treg/CD8 ratios, across multiple tumor types (47). Moreover, in an orthotopic mouse pancreatic cancer model, overexpression of PD-L1 correlated with lack of PTEN expression, which aberrantly activated the PI3K/Akt/mTOR signaling pathway (48). Our findings further indicate that even though these PDACs show molecular and immunophenotypic similarities to the immune-rich subtype, the antitumor effect of the immune response is probably cancelled by the PD-1/PD-L1 blockade, rendering the microenvironment tumor permissive and allowing for the formation of tumor buds, provoking a biologic behavior more similar to the immune-escape phenotype.

The second subpopulation of the immune-exhausted subtype consists of a small number of MSI-H PDACs characterized by an immune-rich microenvironment while showing prevalence of PIK3CA and JAK3 mutations. Consistent with this, a recent study showed that the gene encoding mutLhomolog-1 (MLH1), a key component of the DNA mismatch repair system, is silenced by promoter methylation in cells harboring JAK mutations (49). Interestingly, both subpopulations of the immune-exhausted subtype (PD-L1 and MSI-H) show remarkable microenvironmental similarities. Although these patients may represent good candidates for the administration of checkpoint inhibitors, the unusual characteristics and low frequency of these tumors may explain the limited success of immunotherapy in PDAC (50).

Importantly, we find an increased number of cases with high RHAMM–expressing tumor buds in our immune-escape and immune-exhausted phenotypes as compared with the immune-rich subtype. High RHAMM expression has been associated with adverse prognosis in many carcinomas (51–53) and hematologic malignancies (54); experimental data suggest that targeting RHAMM could inhibit metastatic dissemination and cancer progression (55). In addition, it has been suggested that high RHAMM expression in malignant cells may evoke an antigen-specific antitumoral immune response, making RHAMM a promising target for cancer immunotherapy (56). Our results are consistent with these findings and also with Casalegno-Garduño and colleagues (57) who have shown that persistent RHAMM expression and decreasing CD8⁺ T-cell responses (like in our immune-escape subtype) might indicate the immune escape of leukemia cells.

Interestingly, we also find a good correlation between CA19-9 values and PDAC subtypes as well as with OS and DFS in our cohort. CA19-9, a modified Lewis(a) blood group antigen that is a component of glycoproteins and mucins, is the only established pancreatic cancer biomarker (58). Indeed, cases of the immune-escape phenotype showed significantly higher CA-19-9 values, followed by the immune-exhausted cases, whereas the immune-rich phenotype exhibited the lowest values. Our findings are in accordance with previous studies, where an association between high CA-19-9 values and worse prognosis has been demonstrated (59).

Antitumor immunity is currently believed to be conditioned by the number of neoantigens that represent mutated peptides shed from tumor cells and considered as nonself (2, 60). Tumors with high mutational burdens can thus elicit an antitumor immune response (2). Although the level of neoepitopes is lower in PDAC compared with highly immunogenic tumors, a subset of PDACs harbor a significant neoantigen number (2). Consistent with this, we show that PDACs with an immune-rich microenvironment generally exhibit a higher mutational burden, even in low-range mutations and/or a high prevalence of PIK3CA mutations, which favor higher immune cell infiltrates. Moreover, the balance of the immune cell infiltrates seems to correlate with the molecular and morphologic characteristics, influencing the clinical outcome. In a previous study (15), we have shown that FOXP3 counts have an independent prognostic effect in PDAC. Here, we additionally show that the prognostic impact of FOXP3 is stronger in cases with low-grade tumor budding, resulting in a subset of PDACs (low budding/low FOXP3) with especially good prognosis. On the contrary, tumors with high FOXP3 infiltrates show high-grade tumor budding and an adverse prognosis, although in these cases, the presence of CDKN2A mutations has an additive adverse prognostic impact.

This study may be hampered by the small number of cases and the use of ngTMAs for immune cell and protein evaluation. It has, however, several strong points: The study cohort is well characterized with full clinicopathologic information, multiple tumor punches from center and front were included to account for heterogeneity, and the new AJCC TNM staging system (8th ed.) has been implemented. Moreover, the frequency and type of the identified mutations of the examined genes are comparable with recently published The Cancer Genome Atlas (TCGA) data for pancreatic cancer (18). In addition, the study was designed to conform with the REMARK guidelines (30).

In conclusion, our findings indicate the presence of three different PDAC subtypes with distinct microenvironmental, morphologic, clinicopathologic, and molecular characteristics. The recognition of these subtypes in daily practice has prognostic and predictive implications and would help clinical decision-making toward a more individualized treatment approach for patients with PDAC.

No potential conflicts of interest were disclosed.

The funders had no involvement in the study design; in the collection, analysis, and interpretation of the data; in the writing of the report; and in the decision to submit the article for publication.

Conception and design: B. Gloor, A. Perren, E. Karamitopoulou

Development of methodology: M. Wartenberg, E. Karamitopoulou

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): M. Wartenberg, S. Cibin, I. Zlobec, E. Vassella, L. Terracciano, M.D. Eichmann, M. Worni, B. Gloor, E. Karamitopoulou

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): M. Wartenberg, I. Zlobec, E. Vassella, S. Eppenberger-Castori, L. Terracciano, M.D. Eichmann, M. Worni, B. Gloor

Writing, review, and/or revision of the manuscript: I. Zlobec, E. Vassella, M. Worni, B. Gloor, A. Perren, E. Karamitopoulou

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): M. Wartenberg

Study supervision: E. Karamitopoulou

This project was supported by the Werner and Hedy Berger-Janser Foundation and the Foundation for Clinical-Experimental Tumor-Research.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.