Identification of Transcriptomic Signatures of Pancreatic Ductal Adenocarcinoma–Derived Exosomes That Promote Macrophage M2 Polarization and Predict Prognosis: S100A9 Reveals Tumor Progression

Abstract

Background:

Exosomes play a role in intercellular communication and participate in the interaction between pancreatic ductal adenocarcinoma (PDAC) cells and immune cells. Macrophages can receive tumor cell–derived exosomes to polarize into M2-type macrophages, which can enhance the invasion and metastasis of pancreatic cancer, leading to poor prognosis. However, the mechanism by which pancreatic cancer cell–derived exosomes promote M2-type macrophages is still unclear.

Methods:

M2 macrophage–associated exosome-derived key module genes were identified by differentially expressed genes (DEGs) and weighted gene co-expression network analysis (WGCNA) analysis using exoRbase 2.0, The Cancer Genome Atlas (TCGA), and The International Cancer Genome Consortium (ICGC) databases. Multivariate Cox regression analysis was used to identify key prognostic genes and obtain regression coefficients to establish prognostic signature. Immune infiltration, tumor mutations, and GSEA among different risk groups were compared. exoRbase 2.0, Gene Expression Profiling Interactive Analysis 2 (GEPIA2), HPA, and TISCH2 databases were used to further analyze the expression pattern of S100A9 in pancreatic cancer. In vitro experiments, cell-derived exosome isolation, quantitative polymerase chain reaction (qPCR), western blot, flow cytometry analysis, cell transfection, transwell assay, and CCK-8 assay were applied to investigate the roles of S100A9 in macrophage M2 polarization and tumor progression.

Results:

The key genes of PDAC-derived exosomes promoting M2-type macrophage polarization were identified, and a risk score model was established. The risk score is related to the expression of common immune checkpoints, immune score, and stromal score, and the tumor mutational burden and biological function of high- and low-risk groups were also different. S100A9 was positively correlated with M2-type macrophage marker. In addition, scRNA-seq data from the TISCH2 database revealed that S100A9 is predominantly expressed in pancreatic cancer cells and mono/macrophage cells, suggesting that S100A9 in pancreatic cancer cells could be received by macrophages, thereby inducing macrophage polarization. In vitro, we used exosomes from BxPC-3 cell lines to coculture macrophages and found that macrophages were mainly polarized toward M2 type, which further promoted the proliferation and metastasis of PDAC.

Conclusions:

Our study established a reliable risk score model for PDAC-derived exosomes and M2 macrophages, identified the important role of S100A9 in macrophage M2 polarization, which provides a new strategy for the diagnosis and treatment of PDAC, and strengthened the understanding of the mechanism of tumor development and metastasis.

Introduction

Pancreatic ductal adenocarcinoma (PDAC) is one of the most lethal solid tumors with a poor prognosis (5-year survival rate less than 8%). ¹ Carbohydrate antigen 19-9 (CA19-9) is the only serum biomarker approved by the Food and Drug Administration (FDA), but due to its low sensitivity and specificity, it does not meet the clinical needs. ² Even with the advent of several new diagnostic markers (eg, macrophage inhibitory cytokine-1 and mucin 5AC) and liquid biopsies (eg, ctDNA, Circulating tumor cells, and exosomes), early diagnosis of pancreatic cancer remains difficult.^3-5 Therefore, the diagnosis of PDAC is often confirmed in the advanced stage of the tumor, and the patient loses the opportunity for surgery. A part of patients with PDAC could benefit from immune checkpoint inhibitors (PD-L1 and CTLA4), chimeric antigen receptor T (CAR) cells, immunomodulators, and vaccines. ⁶ In addition, the tumor is surrounded by a dense stroma and a tumor microenvironment (TME) with persistent mild inflammation, which leads to the evasion of immune attack, the ineffectiveness of chemotherapy, and poor prognosis. Among them, tumor-associated macrophages (TAMs) are one of the main components of the TME, accounting for about 15% to 20% of the total cellular tumor mass. ⁷ In general, macrophages can be polarized into M1 or M2 macrophages. M1 macrophages are induced by proinflammatory cytokines such as interferon-gamma (IFN-γ), interleukin (IL)-12, and IL-23. In contrast, TAMs generally exhibit M2-type macrophages activated by T helper 2 (Th2) cytokines such as IL-4, IL-10, and IL-13. ⁸

Tumor-associated macrophages are considered to contribute to tumor cell proliferation, invasion, and metastasis, and previous studies have indicated that high levels of M2 macrophage (CD163+) infiltration are significantly associated with poorer prognosis in patients.^9,10 Hence, targeting TAMs as potential therapeutic targets for tumors has also attracted people’s attention and interest. In addition, in recent years, extracellular vesicles, especially exosomes, have become increasingly important as involved in communication with the tumor stroma. Exosomes are vesicles released by fusion of multivesicular bodies (MVB) with the plasma membrane. The exosomes are 30 to 150 nm in diameter, with a circular bilayer structure and a density between 1.13 and 1.19 g/mL. ¹¹ Under electron microscopy, the exosomes are disk-shaped or cup-shaped. The exosomes contain proteins, lipids, mRNA, miRNA, lncRNA, and other biologically substances, which can participate in intercellular communication, regulation of immune response, cell differentiation, tissue repair and regeneration, tumor growth, and other physiological or pathological processes.^12-14 Exosomes derived from PDAC have been proven to play an important role in this cellular communication. ¹⁵ Exosomes derived from tumor cells can induce macrophage polarization, which in turn further affects the survival, development, and metastasis of tumors. Therefore, it is important to identify key molecular targets in the promotion of M2 macrophage polarization by PDAC-derived exosomes and to develop novel and reliable signatures that could predict prognosis.

Here, we first compared the mRNA transcriptome expression in exosomes of PDAC patients with those of healthy donors from the exoRbase database and obtained 111 differentially expressed genes (DEGs). Then, we downloaded and merged the PDAC sample from The International Cancer Genome Consortium (ICGC) and The Cancer Genome Atlas (TCGA) database to identify key gene modules highly relevant to M2 macrophages. By intermingled DEGs and key module genes, we obtained 36 genes that were highly expressed in exosomes of PDAC and could induce M2 macrophage polarization. Finally, we developed a risk score model based on PDAC-derived exosomes and M2 macrophages and found that S100A9 was highly expressed in PDAC tissues and exosomes, and was significantly positively correlated with M2 macrophages (CD163, CD206, CD301, and CD115), which may improve the understanding of pathogenesis and provide new anti-cancer therapeutic strategies.

Materials and Methods

Data download and processing

ExoRBase 2.0 database

The transcripts per kilobase million (TPM) type expression data from the exosomes of 164 patients with PDAC and 118 healthy donors were downloaded through the exoRBase 2.0 database (http://www.exorbase.org/) and then transformed into log2 (TPM + 1) for subsequent DEG analysis. ¹⁶

TCGA and ICGC database

Transcript sequencing data of PDAC patients with clinical and survival information were downloaded from the TCGA database (https://cancergenome.nih.gov/), and eventually, 178 tumor samples were obtained, and somatic mutation data of patients with PDAC were also obtained. In addition, the Australian cohort (ICGC-PAAD-AU: 90 tumor samples) and the Canadian cohort (ICGC-PAAD-CA: 186 tumor samples) with complete prognostic information were also downloaded. The 3 independent cohorts covered gene transcripts from diverse populations in Europe and the United States. The TPM values were converted to log2 (TPM + 1) for analysis, and the above 3 data sets were merged and normalized using the “SVA” R package, and finally, the normalization results were verified by principal component analysis (PCA). ¹⁷

Gene Expression Omnibus database

The gene expression profiles of PDAC patients with prognostic information were searched in the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). Finally, 4 data sets, GSE28735 (45 tumor samples), GSE57495 (63 tumor samples), GSE62452 (69 tumor samples), and GSE85916 (80 tumor samples), were selected and standardized to be merged in this study. Likewise, the PCA was used to validate the normalized result.

Identification of exosome-derived differentially expressed genes

To identify exosome-derived DEGs between 164 patients with PDAC and 118 healthy donors, we used the “Limma” R package and set the filter condition to P < .05 and|log2FC| > 0.5.

Functional annotation analysis

We used the “clusterProfiler” R package to perform functional enrichment analysis of the above exosome-derived DEGs, focusing on Gene Ontology (GO), biological processes (BP), cell components (CC), molecular functions (MF), and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway. To screen the enrichment information, P < .05 was considered statistically significant.

Protein-protein interaction analysis

Furthermore, we performed protein-protein interaction (PPI) analysis of the exosome-derived DEGs through the Retrieval of Interacting Genes Database (STRING) database (http://string-db.org/) and then calculated the degree of nodes. Finally, the “igraph” R package was used to identify and visualize the central nodes. ¹⁸

Identification of modules associated with M2 macrophages using weighted gene co-expression network analysis

First, 3 data sets (TCGA-PDAC, ICGC-AU, and ICGC-PA) were combined after debating through the “SVA” R package to obtain 454 samples with complete prognostic and gene expression profiles. The proportion of tumor-infiltrating immune cells (TIICs) was counted based on gene expression data using the “Xcell” R package. The 454 samples were then divided into high and low M2 macrophage infiltration according to the median value of the M2 macrophage score, which was used as trait data. Then, the gene expression profiles in the top 25% of variance were screened as input for weighted gene co-expression network analysis (WGCNA). The “WGCNA” R package was used to construct XCell-defined M2 macrophage–related gene coexpression networks. The connectivity of the network was estimated by transforming the adjacency matrix into the topological overlap matrix (TOM). We used the “sft$powerEstimate” function to automatically calculate the power of β = 12 and set the cut height to 0.25 and the minimum module size to 10. After recalculating the module eigengenes (MEs), a module dendrogram describing the relationship between the eigengenes and M2 macrophage traits was drawn. According to the Pearson correlation between MEs and M2 macrophage traits, the module with the highest correlation with M2 macrophage traits was determined as the key module. ¹⁹

Construction of the prognostic signature associated with PDAC-derived exosomes and M2 macrophages

The key module genes related to M2 macrophages identified by WGCNA and the exosome-derived DEGs were interleaved to obtain candidate genes, and these genes were subjected to univariate Cox regression analysis in the 454 PDAC cohort. Only genes with P < .05 were subjected to multivariate Cox regression analysis to identify key prognostic genes and obtain regression coefficients. Risk score = a × gene A × expression level of gene A + b × gene B × expression level of gene B + ····· + x × gene X × expression level of gene X. Here, a, b, and x were the regression coefficients in the multivariate Cox regression analysis. ²⁰

Validation of the prognostic signature

All 454 PDAC samples obtained a risk score according to the previous risk formula. The median of the risk score was set as the cutoff point, and all samples were divided into low- and high-risk subgroups. First, Kaplan-Meier (K-M) survival curves were plotted using the R package “Survival” to determine prognostic differences. The combined 4 GEO data sets (GSE28735, GSE57495, GSE62452, and GSE85916: 257 PDAC samples) were then used for external validation of the risk score model. The receiver operating characteristic (ROC) curves were also performed on the training cohort (454 PDAC samples) and the validation cohort (257 PDAC samples) to verify the prognostic value of the risk score model.

Development and validation of nomograms

In the training cohort, the risk score, age, sex, grade, and stage were used for Cox multivariate regression analysis, and a nomogram integrating risk score and other clinicopathological characteristics were constructed to predict 1-, 3-, and 5-year overall survival (OS) of patients with PDAC. Moreover, calibration curves and ROC curves at 1-, 3-, and 5-year OS were plotted to evaluate the accuracy of the nomogram.

Correlations between risk score model and immunization

In the training cohort (454 PDAC samples), we used the single-sample gene set enrichment analysis (ssGSEA) method of the R software Gene Set Variation Analysis (GSVA) package to calculate 29 immune-related functional analyses between the high- and low-risk groups. In addition, we used the ESTIMATE algorithm to estimate stromal and immune cells in tumor tissues to determine the stromal score, immune cell infiltration score, and tumor purity in 454 PDAC samples, and used the “pheatmap” R package to create a heat map. Finally, the “corrplot” R package was used to analyze the expression level of immune checkpoints (eg, PD-L1) between the high- and low-risk groups.

Analysis of tumor mutations

Somatic mutation data were obtained from the TCGA GDC portal (https://portal.gdc.cancer.gov/). Waterfall plots were then drawn using the “maftools” R package to chart the mutational landscape of patients between high- and low-risk groups.

Gene set enrichment analysis

Gene set enrichment analysis (GSEA) was used to determine the statistical significance of priori specified gene set and consistent heterogeneity between high- and low-risk groups. We set “h.all.v7.4.symbols.gmt” as the reference gene set. An adjusted P < .05 was considered statistically significant in the enriched information screen.

Expression of genes derived from the risk score model

Through the exoRbase 2.0 database, we compared the expression of S100A9, RDX, and AFF3 in PDAC-derived exosomes and healthy donor-derived exosomes, and also compared the expression differences in tissues between patients with PDAC and healthy donors through The Gene Expression Profiling Interactive Analysis 2 (GEPIA2) (http://gepia.cancer-pku.cn/index.html) database. The correlation between S100A9 and M2 macrophage–specific markers (CD163, CD206, CD301, and CD115) was analyzed according to the GEPIA2 database. In addition, The scRNA-seq data of S100A9 in pancreatic cancer were analyzed using the TISCH2 database (http://tisch.comp-genomics.org/home/).

Cell culture

Human pancreatic cancer cell lines Panc-1, Capna-2, BxPC-3, Miapaca-2, CFPAC-1, human pancreatic ductal cell lines HPNE, and human monocytic cell line THP-1 were purchased from the Cell Bank of Chinese Academy of Sciences, Shanghai, China. THP-1 cell line was cultured in 1640 medium (Gibco, Invitrogen, USA). Other cell lines were cultured in DMEM (Gibco, Invitrogen, USA) supplemented with 10% fetal bovine serum (FBS, Lonsera) and 1% penicillin/streptomycin, and incubated in a constant temperature cell culture incubator (CO₂ concentration 5%) at 37°C. THP-1 cells (1 × 10⁶) were induced into macrophages using 100 ng/mL phorbol 12-myristate 13-acetate (PMA; Abcam, UK) treatment for 48 hours.

Cell-derived exosome isolation

Cells were cultured in DMEM medium containing 10% FBS, and when cell density reached about 80%, the cells were changed to FBS-free DMEM medium and cultured for 24 hours. The cell culture supernatant was collected by centrifugation at 500 g for 10 minutes at 4°C. The supernatant was collected and centrifuged at 4°C 12000 g for 20 minutes, and 200 000 g and resuspend the precipitate with 20 mL of phosphate-buffered saline (PBS). Finally, exosomes were obtained by ultracentrifugation at 200 000 g for 2 hours. The collected exosomes were used for subsequent experiments. ²¹

Transmission electron microscopy

Samples dissolved in PBS buffer were dropped into a carbon-coated copper grid and then stained with 2% uranyl-oxalate solution for 5 minutes, and 2% methylcellulose and 4% uranyl acetate for 10 minutes. The image was captured by transmission electron microscopy (TEM) (Hitachi, Japan). ²²

Nanoparticle tracking analysis

Exosome size and concentration were assessed by nanoparticle tracking analysis (NTA) using the Nanosight NS300 (Malvern Panalytical). Exosomes were homogenized by vortex and then diluted 1:500 in filtered phosphate saline buffer and analyzed by NS300. Each sample analysis was performed for 60 seconds. Data were analyzed by Nanosight NTA 3.2 analysis software, with detection thresholds optimized for each sample and a screen gain of 10 to track as many particles as possible with minimal background. Blank 0.1 μm filtered 1× PBS was also used as a negative control. Each sample was analyzed at least 3 times.

Quantitative reverse transcription polymerase chain reaction

Total RNA was extracted using the MiPure Cell/Tissue miRNA Kit (Vazyme, Nanjing) according to the manufacturer’s instructions. RNA sample concentration and purity were measured by Colibri Spectrometer (Germany). The extracted RNA was reverse transcribed to obtain cDNA using SweScript RT I First Strand cDNA Synthesis Kit (Servicebio, Wuhan). Quantitative expression was performed using 2× SYBR Green qPCR Master Mix (High ROX) (Servicebio, Wuhan) on a LineGene 9600 Plus PCR Assay System (Bioer Technology, Hangzhou). GAPDH was used as an internal reference for mRNA. The primer sequences are shown in Table 1.

Table 1.

Primers were used for this study.

Primers	Sequence (5′-3′) forward	Sequence (3′-5′) reverse
S100A9	ACACATCATGGAGGACCTGGACAC	GGTTAGCCTCGCCATCAGCATG
GAPDH	TGACATCAAGAAGGTGGTGAAGCAG	GTGTCGCTGTTGAAGTCAGAGGAG

Western blot

RIPA lysate (Servicebio, Wuhan) was added with 50× Cocktail protease inhibitor, PMSF, and phosphorylated protease inhibitor for extraction of total cellular protein. Protein concentration was measured using the BCA Protein Assay Kit (Beyotime, Shanghai). 5× Loading Buffer (Servicebio, Wuhan) was added 1:4 to the total protein solution and boiled at 95°C for 10 minutes according to the reagent instruction manual. In total, 20 μg of total protein was added to the loading wells of SDS-PAGE gels (Epizyme, Shanghai, China) and transferred to PVDF membranes (Millipore, USA) after completion of electrophoresis. About 5% skim milk was used to seal the whole PVDF. TBST was washed 3 times and incubated with primary antibody at 4°C overnight. The PVDF membrane was removed the next day, washed 3 times with TBST, and incubated with secondary antibody for 1 hour at room temperature before exposure. The antibodies are listed in Table 2.

Table 2.

Antibodies used in western blot.

Antibody name	Manufacturer (item number)
TSG101	Abcam (ab125011)
CD63	Abcam (ab134045)
CD81	Abcam (ab109201)
Calnexin	Abcam (ab92573)
S100A9	Abcam (ab92507)
GAPDH	Proteintech (60004-1-Ig)

Flow cytometry

PMA-treated THP-1 cells and exosome-treated macrophages were suspended in ice-cold PBS (with 10% FCS and 1% sodium azide) and then incubated with CD68 and CD163 (See Table 3) for 30 minutes at 4°C in darkness. Wash the cells by centrifugation at 1000 r for 5 minutes and incubated with FITC-labeled secondary antibody (Servicebio, Wuhan) for 30 minutes at 4°C in the dark. Cells were washed after incubation and analyzed using the flow cytometer Invitrogen AttuneNxt (Thermo Fisher Scientific, USA).

Table 3.

Antibodies used in flow cytometry.

Antibody name	Manufacturer (item number)
CD68	Abcam (ab201340)
CD163	Abcam (ab182422)

Cell transfection

The siRNA (See Table 4) to inhibit S100A9 expression was purchased from GenePharma (Shanghai, China). Following the instructions of siRNA-mate (GenePharma, Shanghai, China), cells were seeded in 6-well plates (2 × 105 well) for 24 hours and then transfected with 100 pmol siRNA for 48 to 72 hours.

Table 4.

siRNA sequences were used for this study.

Name	Sequence
siRNA-S100A9-1	5′-CCTGGACACAAATGCAGACAA-3′
siRNA-S100A9-2	5′-GTCATAGAACACATCATGGAG-3′

Transwell assay

Following manufacturer’s protocol, cell migration and invasion assays were performed using a Transwell 24-well Boyden chamber (Corning, USA) with an 8.0 μm polycarbonate membrane. Matrigel is required for the invasion assay, while the migration assay does not. 3 × 10⁴ tumor cells were suspended in 200 µL of serum-free medium in the upper chamber, and 600 µL of macrophage culture supernatant was added to the lower chamber. After 24 hours in the incubator, the cells were fixed in methanol and stained with crystalline violet.

CCK-8 assay

Cell proliferation was detected using the CCK-8 assay. Based on CCK-8 kit protocol (KeyGEN BioTECH, Nanjing, China), the 96-well plate was added with 1000 cells/well and incubated in the cell culture incubator for 24 hours. Macrophage supernatant was added 100 µL/well. Remove the supernatant at 0, 24, 48, 72, and 96 hours, respectively, add 10 µL/well CCK-8 assay solution, incubate at 37°C for 2 hours, and detect the absorbance at 450 nm by zymography.

Statistical analysis

The Wilcoxon test was performed to compare whether there were statistically significant differences between 2 groups. The differences in OS between 2 groups were analyzed by the log-rank test. Correlation coefficients were calculated by Spearman analysis. All statistical analyses were performed using the R software (version 4.1.0). P < .05 indicates statistical significance.

All experiments were statistically analyzed using SPSS (version 26.0). Plotting was done using GraphPad Prism 9. Data were expressed as mean ± standard deviation (SD). Statistical significance was compared using t test. One-way analysis of variance (ANOVA) was used to compare data from 2 or more groups. The significance level was set as P < .05.

Results

Identifying exosome-derived DEGs and functional annotation

We obtained the transcriptome data of 164 PDAC-derived exosomes and 118 healthy donor-derived exosomes from exoRbase 2.0 database (Figure 1A). Differentially expressed gene analysis between groups obtained 111 DEGs, of which 35 were upregulated and 76 were down-regulated in patients with PDAC compared with healthy donors (Figure 1B). Subsequently, we performed the pathway enrichment analysis of 111 DEGs, and we identified pathways associated with immunity (eg, GO:1903131: mononuclear cell differentiation, GO:0019882: antigen processing and presentation, GO:MHC class II protein complex, and has:Th1 and Th2 cell differentiation) (Figure 1C). Moreover, we constructed the PPI network of DEGs through the STRING database, and the gene networks with threshold weights greater than 0.4 were displayed in Figure 1D.

Figure 1.

Identifying exosome-derived DEGs and functional annotation. (A) The transcriptome data of 164 PDAC-derived exosomes and 118 healthy donor-derived exosomes were downloaded through the exoRBase 2.0 database. (B) Volcano plot of exosome DEGs between patients with PDAC and healthy donors. (C) Functional annotation of 111 exosome DEGs (GO: BP, CC, MF; KEGG). (D) The network of the genes in 111 exosome DEGs.

Identification of M2 macrophage–associated exosome-derived key module genes

A total of 454 PDAC samples from the TCGA-PDAC, PACA-AU, and PACA-CA cohorts were merged through the “SVA” R package, and PCA indicated that the cross-platform batch effect was successfully removed (Figure 2A). We then used the Xcell algorithm to calculate M2 macrophage proportion for 454 patients with PDAC and divided patients with PDAC into high and low macrophage infiltration groups according to a median M2 macrophage proportion score (0.025) (Figure 2B and C). To assess the potential value of M2 macrophages, we performed WGCNA. First, no outlier samples were found through the “Hclust” function (Supplemental Figure 1). The soft threshold β was set to 12 according to the “sft$powerEstimate” function. After constructing the TOM network and setting the height to 0.25, 19 modules were identified by clustering dendrograms, among which the brown module had close correlation strength with M2 macrophage traits (r = 0.42, P = 3e−20), and this module included 2309 genes (Figure 2D and E). In addition, there was a significant relationship between module membership (MM) and gene significance (GS) in brown module (Cor = 0.72, P < 1e−200) (Figure 2F). To explore the mechanism of macrophage M2 polarization mediated by PDAC-derived exosomes, we integrated 111 exosome-derived DEGs and 2309 genes in brown module and finally obtained 36 genes (Figure 2G).

Figure 2.

Identification of M2 macrophage–associated exosome-derived key module genes. (A) PCA of the transcript sequencing data from TCGA-PDAC, PACA-AU, and PACA-CA cohorts. (B) Box plot of M2 macrophage proportion. (C) Heat map of macrophage proportion between the 2 groups. (D) Dendrogram of M2 macrophage–related genes clustered. (E) Heat map of correlation between the MEs and M2 macrophages. (F) Scatterplot between gene salience and module members in brown. (G) Venn plot of 111 exosome-derived DEGs and 2309 brown module genes.

Establishment of prognostic characteristics related to PDAC exosomes and M2 macrophages

The 36 genes above were analyzed by Cox univariate regression analysis, and we identified 4 genes that were significantly associated with prognosis of 454 patients with PDAC (CD44: HR: 1.198, P = .004; RDX: HR: 1.190, P = .001; AFF3: HR: 0.826, P = .009; S100A9: HR: 1.14, P < .001) (Figure 3A and B). Then, Cox multivariate regression analysis was performed on these 4 genes, and RDX, AFF3, and S100A9 were identified as pivotal genes involved in the construction of prognostic signature, and risk score was calculated: risk score = (0.1742 × expression value of RDX) + (−0.2544 × expression value of AFF3) + (0.1028 × expression value of S100A9). Patients with PDAC were divided into high- and low-risk groups according to the median cutoff value of the risk score, and the K-M survival curve showed that the OS of patients in the high-risk group was significantly lower than that in the low-risk group (P < .001; Figure 3C). The point distribution of survival status and risk scores indicated that OS time was also lower in the high-risk group compared with the low-risk group (Figure 3D).

Figure 3.

Establishment of prognostic characteristics related to PDAC exosomes and M2 macrophages. (A) Forest plot of univariate Cox regression analysis of prognosis in 454 patients with PDAC. (B) Kaplan-Meier survival curves of high and low gene expression groups (CD44, RDX, AFF3, and S100A9) (C) Kaplan-Meier survival curves between high-and low-risk groups. (D) Distribution of risk score and OS of 454 patients with PDAC.

External cohort validation of risk score model

We then downloaded and merged 4 PDAC gene expression data sets (GSE28735, GSE57495, GSE62452, and GSE85916: 257 patients with PDAC) with prognostic information from the GEO database, also using the “SVA” R package to remove batch effect of different data sets (Supplemental Figure 1). We applied the above risk score model to 257 PDAC GEO cohort to validate the external prognostic prediction performance, which also found that the OS of the high-risk group was significantly lower than that of low-risk group (Figure 4A). The point distribution of survival status and risk scores indicated that patients with high-risk score also had lower OS (Figure 4B). The area under the curves (AUCs) of the risk score model for 1-, 3-, and 5-year OS were 0.641, 0.638, and 0.650 in training cohort (454 patients with PDAC) (Figure 4C). In the validation cohort (257 patients with PDAC), the AUCs of the risk score model for 1-, 3-, and 5-year OS were 0.626, 0.661, and 0.556 (Figure 4D). To further explore the prognostic value of the risk score model, age, sex, grade, stage, and risk score model were included in the ROC analysis of 1-, 3-, and 5-year OS, respectively, and we found that risk score model had the highest AUC value (Figure 4E-G). Furthermore, we included these above variables (age, sex, grade, stage, and risk score model) in Cox univariate and multivariate regression analysis, and found that risk score model was an independent prognostic factor for OS (HR = 2.606, P < .001; Figure 4H). Finally, we included these variables together to establish nomogram for prognostic prediction (Figure 4I). The calibration curves indicated that the nomogram has excellent prognostic performance (Figure 4J).

Figure 4.

External cohort validation of risk score model. (A) OS curves of GEO cohort (257 patients with PDAC). (B) Distribution of risk score and OS of 257 patients with PDAC. (C) ROC analysis to estimate the prediction value of the risk score model in ICGC and TCGA cohort (454 patients with PDAC). (D) ROC analysis to estimate the prediction value of the risk score model in GEO cohort (257 patients with PDAC). (E-G) The prognostic value of the risk score model, age, sex, grade, stage, and risk score model was included in the ROC analysis of 1-, 3-, and 5-year OS, respectively. (H) Forest plot of univariate and multivariate Cox regression analysis of prognosis in 454 patients with PDAC. (I) Nomogram predicting OS for patients with PDAC. (J) Calibration curve of nomogram for OS.

Correlations between risk score model and immunization

The ESTIMATE scores demonstrated higher stromal and immune scores in the high-risk group than in the low-risk group, while tumor purity was lower than in the low-risk group (all P < .05) (Figure 5A-E). According to 29 immune-related functional analyses, some parameters were significantly different between high- and low-risk groups, such as aDCs, Parainflammation, MHC_class_I, DCs, APC_co_stimulation, Check-point, T_cell_co-inhibition, Th1_cells, CCR, Treg, Macrophages, T_helper_cells, Th2_cells, and inflammation promoting (all P < .05; Figure 5A). Finally, the heat map of expression of common immune checkpoints between high- and low-risk groups indicated significant differences in the expression of CD274, POLE2, FEN1, MCM6, MSH6, MSH2, FAP, and LOXL2 between high- and low-risk groups (all P < .05; Figure 5F).

Figure 5.

Correlations between risk score model and immunization. (A) Heat map of 29 immune-related functional and ESTIMATE score between high- and low-risk groups. (B-E) Differences in the distribution of tumor purity, immune score, stromal score, and ESTIMATE score between high- and low-risk groups. (F) The heat map of expression of immune checkpoints between high- and low-risk groups.

Analysis of tumor mutations and GESA between high- and low-risk groups

An oncoplot revealed the top 20 genes with the highest mutation frequency in the high- and low-risk groups (Figure 6A and B). The tumor mutational burden (TMB) was significantly higher in the high-risk group compared with the low-risk score (87.64% vs 76.71%). Furthermore, PDAC-related mutated genes KRAS, TP53, and CDKN2A were more frequently mutated in the high-risk group, while KRAS, TP53, and SMAD4 were more frequently mutated in the low-risk group (Figure 6A and B). Gene set enrichment analysis was carried out in training cohort (454 patients with PDAC) between high- and low-risk groups. The gene sets of the high-risk group were enriched in tumor proliferation, metabolism, and metastasis-related pathways (eg, E2F_TARGETS, EPITHELIAL_MESENCHYMAL_TRANSITION, G2M_CHECKPOINT, HYPOXIA, and INFLAMMATORY_RESPONSE) (Figure 6C), while the gene sets of the low-risk group were enriched in pancreas beta cells PANCREAS_BETA_CELLS (Figure 6D).

Figure 6.

Analysis of tumor mutations and GESA between high- and low-risk groups. Correlation between the high-risk group (A) and the low-risk group (B) with tumor mutations. GESA between the high-risk group (C) and the low-risk group (D).

Expression analysis of genes for the construction of risk score model

In the exoRbase 2.0 database, expression of S100A9 and RDX was significantly higher in exosomes from patients with PDAC than in exosomes from healthy donors, while AFF3 was significantly lower than in exosomes from healthy donors (all P < .05; Figure 7A). Through the GEPIA2 database, only S100A9 was found to be significantly higher in PDAC tissues than in healthy tissues (AFF3: P > .05; RDX: P > .05; and S100A9: P < .05) (Figure 7B-D). Furthermore, the protein expression level of S100A9 in pancreatic cancer tissues was also significantly higher than that in non-tumor tissues using the HPA database (Figure 7E). Finally, we found a significant positive correlation between S100A9 and the specific markers (CD163, CD206, CD301, and CD115) of M2 macrophages through the GEPIA2 database (Figure 7F). Finally, through scRNA-seq data sets (PAAD_CRA001160), we discovered that S100A9 is highly expressed in both cancer tissues and mono/macrophage cells (Figure 7G-I). Moreover, compared with adjacent non-cancerous tissues and corresponding macrophages, cancer tissues and macrophages within the cancer tissue exhibit significantly elevated expression of S100A9, respectively (Figure 7J).

Figure 7.

Expression analysis of genes for the construction of risk score model. (A) The transcriptome expression levels of S100A9, RDX, and AFF3 in exosomes between patient with PDAC and healthy donors from exoRbase 2.0 database. (B-D) The transcriptome expression levels of AFF3, RDX, and S100A9 in tissues between patient with PDAC and healthy donors from the GEPIA2 database. (E) The protein expression levels of S100A9 in tissues from the HPA database. (F) The correlation between S100A9 and M2 macrophage markers (CD163, CD206, CD301, and CD115). (G, H) Cell clustering annotations of UMAP plot in PAAD_CRA001160 data set. (I, J) UMAP and violin plots of S100A9 in PAAD_CRA001160 data set.

Identification of pancreatic cancer cell–derived exosome extraction and S100A9 expression

The cell culture supernatant without FBS was collected. Exosomes were extracted from the cell supernatant by differential and ultra-high speed centrifugation, and typical cup-shaped (Figure 8A) vesicles of approximately 50 to 150 nm in size (Figure 8B) were visible under TEM. Meanwhile, the exosome markers TSG101, CD81, and CD63 were positively expressed (Figure 8C and Supplemental Figure 2), and calnexin was negatively expressed in exosomes and positively expressed in cells (Figure 8C and D). Western blot (Wb) results also showed that S100A9 protein expression in pancreatic cancer cells and their secreted exosomes was higher than that in normal pancreatic ductal epithelial cells (Figure 8C and D). The results of quantitative reverse transcription polymerase chain reaction (qRT-PCR) also showed increased expression of S100A9 in pancreatic cancer cells and their secreted exosomes (Figure 8E and F).

Figure 8.

Extraction and identification of exosomes of pancreatic cancer cell and pancreatic ductal epithelial cell origin with S100A9 expression analysis. (A) TEM showing pictures of exosomes (scale bar: 100 nm). (B) NTA analysis of the particle size of exosomes. (C) WB analysis of exosome marker proteins and S100A9 proteins. (D) WB analysis of S100A9 proteins in cellular secretory exosomes. (E) qRT-PCR detection of S100A9 expression in pancreatic cancer cell lines and normal pancreatic ductal epithelial cells (**P < .05). (F) qRT-PCR detection of S100A9 expression in cellular secretory exosomes (**P < .01; *P < .05).

Pancreatic cancer cell–derived exosome S100A9 induces macrophage polarization and promotes pancreatic cancer cell proliferation, migration, and invasion

To further verify whether exosome could induce macrophage polarization and the effect of polarization on pancreatic cancer cell function. We co-cultured the extracted exosomes with macrophages and then cultured pancreatic cancer cells using macrophage culture supernatants from the co-cultured cells (Figure 9A). Macrophages were obtained by PMA-induced THP-1 and detected by flow cytometry for CD68 (Figure 9B). Co-culture of exosomes with macrophages resulted in increased expression of CD163 (Figure 9C), indicating that exosomes from pancreatic cancer cells induced M2 polarization of macrophages. M2 macrophages further promoted pancreatic cancer cell proliferation, migration, and invasion (Figure 9D and E).

Figure 9.

Pancreatic cancer cell–derived exosome S100A9 induces macrophage polarization and promotes pancreatic cancer cell proliferation, migration, and invasion. (A) Protocol of tumor cell–derived exosomes induces macrophage polarization to promote a malignant tumor phenotype. (B) Induction of macrophages. (C) Detection of macrophage M2 polarization by flow cytometry. (D) CCK-8 assay detects cell proliferation (****P < .0001; ***P < .001). (E) Transwell assay detects cell migration, invasion (*P < .05). (F, G) qRT-PCR and WB detection of mRNA and protein expression after S100A9 knockdown. (H) Effect on macrophage polarization after knockdown of S100A9 in exosomes. (I, J) Effect of S100A9 knockdown on cell proliferation and migration invasion (**P < .001; ****P < .0001).

To further explore the role of exosomal S100A9, we used siRNA to interfere with the expression of S100A9 in the highest-expressing BxPC-3 cells (Figure 8C and F). Quantitative reverse transcription polymerase chain reaction and WB results showed that the levels of S100A9 were down-regulated in both the cells and their secreted exosomes (Figure 9F and G). Macrophage M2 polarization was reduced after knockdown of exosome S100A9 (Figure 9H). In addition, the proliferation, migration, and invasion ability of pancreatic cancer cells were diminished (Figure 9I and J). In conclusion, exosome S100A9 induces macrophages to undergo M2 polarization promoting a malignant phenotype in pancreatic cancer cells.

Discussion

Pancreatic ductal adenocarcinoma is highly aggressive and metastatic and has immune escape and chemo-resistant properties, in which the tumor stroma plays a key role. Macrophages are an important component of the stroma in PDAC, and the M2 phenotype macrophages are more associated with poor clinical outcomes in PDAC, which contribute to tumor progression and metastasis. ²³ The mechanism regarding the polarization and infiltration of M2 macrophages into pancreatic cancer is unclear. However, the role of exosomes in intercellular communication in cancer development is being extensively studied in recent years. Exosomes released by tumor cells have the effects of anti-apoptosis, promoting blood vessel growth, changing the blood-brain barrier, and changing cell metabolism, which is conducive to tumor growth and metastasis.^24-27 Previous studies have found that macrophages phagocytizing exosomes released by pancreatic cancer cells can be polarized into M2 macrophages, thereby promoting tumor cell growth and metastasis. ²⁸ As exosomes belong to a class of extracellular vesicles, which can carry a variety of mRNA, proteins, and nucleic acids. ²⁹ It is of great significance to identify the key targets in PDAC exosomes that could induce M2 polarization of macrophages, which could provide a novel and reliable therapeutic target.

In this study, we identified 111 DEGS in serum exosomes from patients with PDAC and health donors by exoRbase database. Then, by combining mRNA expression data of 454 PDAC from TCGA and ICGC databases, we classified patients with PDAC into high and low M2 macrophage groups by Xcell algorithm, followed by WGCNA to identify M2 macrophage phenotypes with highly associated brown modules. To integrate the expression information of PDAC-derived exosomes and M2 macrophages, we intersected the DEGs and brown module and finally obtained 36 genes. Then, we performed Cox regression analysis based on TCGA and ICGC cohorts and identified 3 signature genes (AFF3, RDX, and S100A9) which were used to construct a prognostic risk score model associated with PDAC exosomes and M2 macrophages, and the high- and low-risk score groups could well differentiate the OS of PDAC. In addition, external verification was performed through 4 GEO data sets. Then, it was also found that the risk score model had a better AUC value than other clinical variables, and the risk score was included with the clinical variables to establish a nomogram, which had a good ability to predict 1-, 3-, and 5-year OS, which may enable clinicians to determine the prognosis of individual patients. As far as we know, this is the first risk score model for prognostic features constructed from the PDAC exosome, M2 macrophage polarization-associated gene set. In addition, in our risk score model, we found that the immune score and stromal score of the high-risk group were significantly higher than those of the low-risk group, and there were also significant differences between the 2 groups in the analysis of immune-related functions and the expression of immune checkpoints (CD274), which is helpful in identifying immunotherapy-sensitive patients, guiding patient treatment regimens, and improving the therapeutic benefits of clinical PD-L1 inhibitors (atalizumab and divalirudin). There are also differences in tumor mutation and functional enrichment between the high- and low-risk groups. Although our risk score model ultimately did not include CD44, CD44 plays an important role in tumor invasion and metastasis. CD44 has been shown to increase the activity of c-Met and inhibit the Hippo signaling pathway, which is associated with proliferation, invasion, metastasis, differentiation, and drug resistance in pancreatic cancer.^30-32 Then, we compared the expression levels of 3 signature genes (AFF3, RDX, and S100A9) in exosomes and tissues, and interestingly, only S100A9 was found to be highly expressed in exosomes and tumor tissues, and the HPA data also showed that S100A9 was also higher in cancer tissues than in paracancerous tissues at the protein level. Finally, we compared the mRNA expression levels of 3 signature genes (AFF3, RDX, and S100A9) in exosomes and tissues. Interestingly, only S100A9 was found to be highly expressed in exosomes and tumor tissues, and the protein level of S100A9 was also higher in cancer tissues than in adjacent tissues according to HPA database. Furthermore, we analyzed the correlation between S100A9 and M2 macrophage marker (CD306 and CD163) expression based on the TCGA cohort and explored the expression pattern of S100A9 in scRNA-seq data. The results showed a significant positive correlation between S100A9 and M2 macrophage markers, and S100A9 was mainly expressed in pancreatic cancer cells and mono/macrophage cells, suggesting that S100A9 is an important factor in promoting macrophage M2 polarization. In vitro, we found that S100A9 expression was significantly higher in pancreatic cancer cells (PANC-1, CAPAN-2, BXPC-3, MIPACA-2, and CFPAC-1) than in normal pancreatic cells (HPNE), both in exosomes and cells.

Finally, we co-cultured macrophages using exosomes derived from the highest S100A9-expressing BxPC-3 cell line and found that macrophages were predominantly polarized toward the M2 type. M2 macrophages have been shown to promote tumor progression. ³³ This is consistent with our results, where cell proliferation, migration, and invasion assays revealed that macrophage M2 polarization promoted a malignant phenotype in PDAC cells. However, this procedure was attenuated with the knockdown of S100A9 in exosomes. Therefore, we believe that pancreatic cancer exosomes as carriers carrying S100A9 were phagocytosed by macrophages, which induced M2 polarization of macrophages and accelerated tumor progression. Similarly, Wang et al ³⁴ had found that exosome S100A9 could promote stemness and growth of CRC cells. Zhu et al ³⁵ found that exosome S100A9 could be used as a diagnostic marker for hepatocellular carcinoma.

S100A9 is a calcium-binding protein, which could exist in the form of homodimer or heterodimer with S100A8.^36,37 S100A9 is also upregulated in many other cancers, including breast, colon, hepatocellular, gastric, colorectal, non-small-cell lung, and cervical cancers. S100A9 can bind to cell surface receptors to trigger signaling pathways associated with many cellular processes, including proliferation, differentiation, and migration. In addition, S100A9 expression can be detected in infiltrating macrophages in rheumatoid arthritis and other inflammatory diseases. ³⁸ The secreted S100A9 protein has been reported to play a role in establishing an environment conducive to cancer growth.^39,40 This is in line with our finding that this study also identified that S100A9 in tumor cell–derived exosomes can be involved in M2 macrophage polarization as cellular communication, ultimately leading to poor prognosis of PADC. Our study for the first time suggested that S100A9 contributes to M2 macrophage polarization and invasion through tumor cell–derived exosomes pathway, which provided a new strategy for the treatment of PDAC and strengthened the understanding of the mechanism of PDAC developmental metastasis.

Conclusions

In this study, 3 signature genes (AFF3, RDX, and S100A9) were identified to establish a risk score model associated with PDAC exosomes and M2 macrophages, which has a good ability to predict OS. In addition, the potential role of tumor cell–derived S100A9 on macrophage M2 polarization was explored. Nevertheless, further experimental and clinical validation in different centers and larger cohorts is needed.

Supplemental Material

sj-docx-1-onc-10.1177_11795549241239042 – Supplemental material for Identification of Transcriptomic Signatures of Pancreatic Ductal Adenocarcinoma–Derived Exosomes That Promote Macrophage M2 Polarization and Predict Prognosis: S100A9 Reveals Tumor Progression

Supplemental material, sj-docx-1-onc-10.1177_11795549241239042 for Identification of Transcriptomic Signatures of Pancreatic Ductal Adenocarcinoma–Derived Exosomes That Promote Macrophage M2 Polarization and Predict Prognosis: S100A9 Reveals Tumor Progression by Siyuan Tan, Haodong Tang, Zheng Zhang, Yang Wang, Haifeng Li, Wenyuan Shi, Hao Ye, Peng Xie and Jiahua Zhou in Clinical Medicine Insights: Oncology