We retrospectively assessed fresh frozen (FF) tumor tissue samples from patients who underwent resection for histologically proven PDAC between December 2004 and December 2016, at the Erasmus MC University Medical Center (EMC) in Rotterdam, the Netherlands. Patients were screened for eligibility based on their survival time. We selected treatment naïve tumors of patients following surgical resection for their PDAC that remained recurrence-free and survived for at least 3 years (36 months, long-term survivors). We compared these tumors to those of patients who survived less than 6 months (short-term survivors) due to tumor recurrence. We excluded patients who died from postoperative complications and other causes. Furthermore, patients with neuroendocrine, duodenal, distal-bile duct, and ampullary carcinoma were also excluded.

Clinical, histopathological, and laboratory data were retrieved from the electronic medical records. Information obtained from pathology reports included: tumor grade (well, moderate, or poor), lymph node status (positive and negative), tumor location (head, body, or tail), tumor stage (according to the AJCC 8^th edition), and margin status (radical (R0) vs non-radical (R1; ¾ 1mm)). From the laboratory data, we collected baseline cancer antigen (CA) 19-9 (kU/L), carcinoembryonic antigen (CEA) (ng/ml) and, calculated the systemic immune inflammation index (SIII) (29). The study was approved by the Medical Ethical Committee of EMC (MEC-2020-0252).

Intergroup differences in baseline characteristics (short-term vs. long-term survivors) on continuous variables were determined using the non-parametric Mann-Whitney U test, and categorical data were compared using Fisher’s exact test. Cancer-specific survival and recurrence-free survival were calculated from the date of surgery to the date of the event (death from cancer or recurrence of cancer, respectively). In the case of no event, the information of the patients was censored at the date of the last follow-up. Patients were followed according to the standard of care guidelines (30). Follow-up information was retrieved through the electronic medical records and by contacting patients’ general practitioners in the event of missing information. Significance for statistical tests was inferred at a p-value of <0.05. All statistical analyses, if not mentioned otherwise, were performed using SPSS (version 26.0).

To validate our results, an independent cohort of samples was collected from La Paz University hospital of Madrid. Formalin-Fixed, Paraffin-Embedded (FFPE) samples of patients who had the same characteristics as the discovery cohort, were selected. The study was approved by the Medical Ethical Committee of the University of Madrid (PI3468). The overview of the study is presented in ( Figure 1 ).

Fresh Frozen (FF) tumor samples of the included PDAC patients were sectioned with 5 µM thickness, stained with hematoxylin and eosin (H&E), and examined by a pathologist who highlighted the location of the tumor-rich areas. Each sample was then further sectioned, and a total of 60 µM was used for RNA isolation using the Fresh Frozen Micro Kit (QIAGEN, Hilden, Germany) according to the manufacturer’s protocol. The quality and quantity of RNA samples were examined using the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, California, USA). RNA concentration was corrected to include the proportion of the sample that was >300 base pairs.

To measure the differences in the immune infiltration, the nCounter^® PanCancer Immune Profiling Panel (NanoString Technologies, Seattle, WA) was used. This panel consisted of 730 cancer immune-related genes, 40 housekeeping genes, 6 positive and 8 negative controls. RNA samples of 200 ng/5-7 uL were hybridized to target probes for 17 hours at 65°C. RNA expression levels were measured using the nCounter^® FLEX Instrument (31). The counting of the genes was done by scanning 490 Fields Of Views (FOV).

Gene count data were analyzed using the nSolver™ Analysis Software (version 4.0) and its nSolver™ Advanced Analysis module (version 2.0) (www.nanostring.com). The quality control of the measurements was done according to the general workflow used in the nSolver™ Analysis Software (32). The expression levels of the negative controls were used to determine the detection limit by calculating the average of the negative control expression plus double the standard deviation. The expression data were normalized using the most stable housekeeping genes using the geNorm algorithm (33). The differentially expressed genes were identified using the recommended workflow of the nSolver™ Advanced Analysis, which includes a mixture of negative binomial models, simplified negative binomial models, or log-linear models based on the convergence of each gene. Genes were regarded differentially expressed if the p-value was < 0.05 and the |Fold-Of-Change (FOC)| > 2.

To identify the immune cell types, candidate marker genes that define the cell types were tested for their specificity and stable expression, following our previously described method (34). The assumption for appropriate marker genes is that all selected genes for the same cell type decrease and increase in the same direction. This was tested by calculating the correlation of determination (R² ) and the slope for each pair of genes within a cell type definition. The gene pairs with a correlation of determination between 0.4 and 0.6 were checked separately using scatter plots to calculate their slope. The slope should be above 0.75 and below 1.25 which is equivalent to being between 33.75 and 56.25 degrees to indicate high pairwise similarity (35) between the candidate genes and to be included in the analysis ( Supplementary Table 1 ). Genes that did not pass the above criteria were excluded from cell type profiling analysis by changing the Probe Annotation and cell types contrast matrix in nSolver™ Advanced Analysis. Independent two-sample t-tests were conducted between the short-term and long-term survivors using the cell type scores calculated in the cell type profiling analysis.

FFPE samples of the same PDAC samples used in the targeted gene expression profile were used for validation in the GeoMx™ Digital Spatial Profiling analysis. A total of 10 samples were retrieved (5 samples from long-term and 5 samples from the short-term survivors). From each sample, one section of 5µM was used to be stained with GeoMx Solid Tumor TME Morphology Kit Human Protein Compatible (NanoString Item # 121300301) that contains: Pan-Cytokeratin (PanCK, 647 nm), alpha-smooth muscle actin (αSMA, 488 nm), CD45 (594 nm). Nuclear stain CYTO-13 (NanoString Item # 121300303) was used to detect all cells (532 nm). At the same time, the sections were incubated with a cocktail of 53 photo-cleavable, oligo-labeled primary antibodies (the Onco-Immune protein panel: 4-1BB, αSMA, B7H3, Bcl2, beta-2-microglobulin, CD11c, CD20, CD127, CD14, CD163, CD25, CD27, CD3, CD34, CD4, CD40, CD45, CD45RO, CD56, CD66b, CD68, CD8, CD80, CD86, CTLA4, EpCAM, ER-alpha, FAPa, FoxP3, GITR, GZMB, fibronectin, Her2/ErbB2, Ox40L, HLA-DR, ICOS, IDO1, Ki-67, LAG3, MART-1, NY-ESO-1, OX40-L, Pan-Cytokeratin, PD-1, PD-L1, PD-L2, PR, PTEN, S100b, STING/TMEM173, TIM-3, and VISTA). The incubation with antibodies was done overnight according to the manufacturer’s protocol. Afterward, slides were loaded onto the DSP instrument. Each slide was first scanned to produce a digital image of tissue morphology based on the fluorescent markers. Next, the regions of interest (ROI) were selected which was guided by the morphological markers. Finally, the photocleaved oligos from the spatially resolved ROIs were hybridized to the corresponding probes for 17 hours at 65°C and scanned using 280 fields of views (FOV). Microplates were analyzed using DSP data center software (www.nanostring.com).

Multiple ROIs were selected based on morphological markers including tumor-rich ROIs with a high morphological expression of PanCK, desmoplastic with a high morphological expression of αSMA, and immune-rich containing a high expression of CD45. However, some ROIs contained more than one feature and were therefore relabeled/classified as “tumor plus desmoplasia” or “desmoplasia plus CD45” ROIs. In addition, some ROIs contained all features, so they were categorized as “tumor plus desmoplasia plus CD45” ( Supplementary Figure 1 ). DNA staining was used to ensure selecting ROIs that contained cells. The size and the number of cells included in selected ROIs can vary.

The protein expression data generated by the GeoMx™ DSP were first normalized for technical variation during hybridization by using the positive controls included in the experiment (The External RNA Control Consortium (ERCC)). Subsequently, a normalization based on housekeeping genes that showed a stable expression across all samples (Histone 3 and S6 Ribosomal Protein) was conducted. GAPDH, a third housekeeping gene, was excluded from being used in normalization because it was not stably expressed in all ROIs and all samples ( Supplementary Figure 1G ). These two subsequent normalizations enables comparing ROIs of different size and number of cells to each other. The first comparison was performed between ROIs that contained tumor (i.e. tumor, tumor plus desmoplasia, and tumor plus desmoplasia plus CD45) and ROIs that contained stroma (i.e. desmoplasia, desmoplasia plus CD45, and CD45) in long-, and short-survivors separately ( Supplementary Figure 1 and Supplementary Table 2 ). The non-parametric Mann-Whitney U test was used to determine differentially expressed proteins. After that, a correlation analysis using Spearman’s rank-order correlation was run for all the significant proteins in the tumor-rich and stromal ROIs aiming to identify immune cell types that infiltrate together in various areas of PDAC tissue samples based on the survival of the patients. The correlations were regarded to be positively significant when the correlation coefficient rho (r) > 0.5, and negatively significant when r < -0.5 (36). Finally, to identify the cell types that were infiltrated with B cells, a correlation analysis was performed on all detected proteins in ROIs that showed high CD20 expression.

A total of 22 FFPE samples were used to confirm our results (n=12 long-term survival, n=10 short-term survival). Samples were sectioned (5µm thickness) and stained with Pan-Cytokeratin antibody that is available in the GeoMx Solid Tumor TME Morphology Kit Human Protein Compatible (NanoString Item # 121300301) (PanCK, 647 nm), Nuclear stain CYTO-13 (NanoString Item # 121300303) was used to detect all cells (532 nm). Anti-CD20 antibody (Abcam Recombinant Anti-CD20 antibody [EP459Y] (ab78237) was labeled with antibody labeling kit from Invitrogen Thermo Fisher called Alexa Fluor™ 594 Antibody Labeling Kit Item number A20185 (594 nm). Sectioned samples were hybridized with the immune cell profiling panel, consisting of 20 antibodies (PD-1, CD68, HLA-DR, Ki-67, Beta-2-microglobulin, CD11c, CD20, CD3, CD4, CD45, CD56, CD8, CTLA4, GZMB, PD-L1, PanCk, SMA, Fibronectin). The panel contains three positive controls (Ribosomal protein S6, Histone 3, and GAPDH) and three negative controls (two mouse and one rabbits IgGs) that were used in each measurement. Samples were prepared according to the manufacturer’s protocol (NanoString, GeoMx DSP, Seattle, USA). The GeoMx DSP instrument was used to perform the measurements and the software (v2.2).

A total of twenty samples of PDAC patients were included in the discovery cohort, 10 in the long-term survival group, and 10 in the short-term survival group. Median cancer-specific survival was 53.1 and 4.7 months, respectively, whereas the median time-to-recurrence was 46.3 and 3.7 months, respectively. The main characteristics of the patients were matched to the best of our ability. There were no significant differences in age, gender, tumor location, operation procedure, tumor differentiation, lymph node status (positive vs. negative), margin status (R0 vs. R1), T-stage (T1-T2 vs. T3), adjuvant systemic therapy, SIII, CEA, or CA19-9. Baseline clinicopathologic characteristics are summarized in Supplementary Table 3 .

In addition, a total of twenty-two samples of PDAC patients were included in the validation cohorts, 12 in the long-term survival group, and 10 in the short-term survival group. Also for these samples, the main characteristics of the patients were matched to the best of our ability. The clinicopathologic characteristics of the independent validation cohort were similar, except for T-stage, to the discovery cohort of samples and were summarized in Supplementary Table 4 . However, when comparing the short-term survival group to the short-term survival group, and the long-term survival group to the long-term survival in the discovery cohort and the validation cohort, respectively, tumor grade was found to be significantly different in the short-term survival groups.

Of the 730 immune-related genes, 631 genes were above the detection threshold (i.e. 21 counts), and 99 of these genes were found to be differentially expressed between the two groups (short-term vs. long-term survivors) with a p-value < 0.05. Thirty-nine genes had |FOC|>2, among which 14 genes were over-expressed, and 25 genes were under-expressed in long-term survivors compared to short-term survivors ( Supplementary Table 5 ; Figure 2A ). Four of the over-expressed genes in long-term survivors were related to the function of B cells, namely, CR2, CD19, CD79B, and BLNK ( Supplementary Table 5 ). Only three genes were found to be differentially expressed genes after multiple corrections, therefore, we presented the differences between long- and short-term survivals based on immune cell types.

Based on the cell type profiling, the total infiltration of immune cells (PTPRC (CD45⁺)) was comparable between long- and short-term survivors (p=0.309, Figure 2B ). A significantly higher B cell score (CD19, CD22, CR2, and MS4A1 (CD20)) was found in long-term survivors compared to short-term survivors (p=0.018, Figure 2C ). None of the other cell types were found to be significantly different between the two groups ( Supplementary Figure 2 ). The higher infiltration of B cells in long-term survivors was confirmed by using the GeoMx™ DSP analysis using the FFPE samples corresponding to the discovery cohort ( Table 1 ). B cells were found to be significantly higher in the stroma of the long-term survival samples (p=0.038, Figure 2D ). They also were found to be higher within tumor cells of the same group (p=0.049, Figure 2E ). The comprehensive analysis that included the desmoplastic and the tumor-rich area showed a significantly higher infiltration of B cells in the long-term survival group of samples (p=0.002, Figure 2F ).

Tumor ROIs that expressed high levels of CD20 in long-term survival showed a high expression of various T cell markers like CD3, CD4 and CD8 and of HLADR marker that might be expressed by memory B cells or dendritic cells (DCs) or antigen-presenting cells (APCs) ( Supplementary Table 6 ). Moreover, the expression of CD20 was associated with the high expression of CD34 which is a bone marrow stem cell marker that is expressed by progenitor cells of blood vessels and stromal tissue and may reflect the increased number and quality of blood vessel formation in long-term survivors (37, 38). Furthermore, a strong correlation between the expressions of CD20 and CD27 (r= 0.752, p<0.001) and Bcl2 (r= 0.734, p<0.001) was found in tumor areas of long-term survivors, which may reflect a specific subtype of B cells like memory B cells. All proteins which were correlated with CD20 expression are summarized in Supplementary Table 6 ; Supplementary Figures 3 , 4 .

Since we found higher infiltration of CD20-positive cells in the tumor tissue of the long-term survival group compared to the short-term survival group, and high correlations were found between CD20-positive cells and CD3, CD4, and CD8-positive cells, we examined high CD20-positive ROIs between short and long-term survival groups. We found that the expression of CD20 correlates negatively with the expression of αSMA (r=-0.656, p=0.002) and with fibronectin (r=-0.540, p=0.014) in long-term survivors. Positive correlations were found between CD20 and Bcl2 (r=0.695, p=0.001), CD3 (r=0.598, p=0.550), CD4 (r=0.550, p=0.012), CD8 (r=0.567, p=0.009) and Ki67 (r=0.558, p=0.011), which confirm our previous findings in the tumor areas.

Tumor ROIs of short-term survivors revealed a high expression of molecules correlated with immune suppression such as PD-1, PD-L1, CTLA4, LAG3, FoxP3, and CD25 compared to long-term survivors ( Table 1 ). The same suppressive molecules were found to be highly expressed in desmoplastic areas ( Table 2 ). In addition, the level of αSMA and FAPa expression was significantly higher in desmoplastic areas of short-term survivors compared to long-term survivors (p=0.001, Table 2 ; Figures 3A, C ). Importantly, our results do not show differences in the percentage of the desmoplastic areas between the two groups, it shows differences in the desmoplastic protein expression levels in very similar sizes of desmoplastic ROIs in short- and long-term survivors. The comprehensive analysis that included the desmoplastic and the tumor-rich area showed the same results ( Figures 3B, D ). All differentially expressed and correlated proteins in tumor and stromal areas are summarized in Supplementary Tables 6 , 7 .

Interestingly, we found a strong negative correlation between αSMA and CD20-positive cells (r=-0.749, p<0.001, Supplementary Table 7 ) in the stromal ROI of the short-term survivors. Negative correlations were also found between αSMA and other immune effector cells such as CD8-positive cells (r=-0.505, p=0.014, Supplementary Table 7 ), suggesting that these effector immune cells are unable to penetrate the stroma of short-term survivors. On the other hand, immunosuppressive cells, and receptors such as FoxP3, CTLA4, and PD-L1 were all positively associated with αSMA. Analyzing CD20-rich areas, we found similar results in long-term survivors. That is, CD20 was negatively associated with αSMA (r=-0.762, p=0.028), CTLA4 (r=-0.762, p=0.028), Fibronectin (r=-0.738, p=0.037), and GZMB (r=-0.810, p=0.015).