M2‐phenotype tumour‐associated macrophages upregulate the expression of prognostic predictors MMP14 and INHBA in pancreatic cancer

Abstract Pancreatic cancer is one of the most lethal gastrointestinal tumours, the most common pathological type is pancreatic adenocarcinoma (PAAD). In recent year, immune imbalanced in tumour microenvironment has been shown to play an important role in the evolution of tumours progression, and the efficacy of immunotherapy has been gradually demonstrated in clinical practice. In this study, we propose to construct an immune‐related prognostic risk model based on immune‐related genes MMP14 and INHBA expression that can assess the prognosis of pancreatic cancer patients and identify potential therapeutic targets for pancreatic cancer, to provide new ideas for the treatment of pancreatic cancer. We also investigate the correlation between macrophage infiltration and MMP14 and INHBA expression. First, the gene expression data of pancreatic cancer and normal pancreatic tissue were obtained from The Cancer Genome Atlas Program (TCGA) and The Genotype‐Tissue Expression public database (GTEx). The differentially expressed immune‐related genes between pancreatic cancer samples and normal sample were screened by R software. Secondly, univariate Cox regression analysis were used to evaluate the relationship between immune‐related genes and the prognosis of pancreatic cancer patients. A polygenic risk score model was constructed by Cox regression analysis. The prognostic nomogram was constructed, and its performance was evaluated comprehensively by internal calibration curve and C‐index. Using the risk model, each patient gets a risk score, and was divided into high‐ or low‐ risk groups. The proportion of 22 types of immune cells infiltration in pancreatic cancer samples was inferred by CIBERSOFT algorithm, correlation analysis (Pearson method) was used to analyse the correlation between the immune‐related genes and immunes cells. Then, we applied macrophage conditioned medium to culture pancreatic cancer cell line PANC1, detected the expression of MMP14 and INHBA by qRT‐PCR and Western blot methods. Knock‐down MMP14 and INHBA in PANC1 cells by transfected with shRNA lentiviruses. Detection of migration ability of pancreatic cells was done by trans‐well cell migration assay. A subcutaneous xenograft tumour model of human pancreatic cancer in nude mice was constructed. In conclusion, an immune‐related gene prognostic model was constructed, patients with high‐risk scores have poorer survival status, M2‐phenotype tumour‐associated macrophages (TAMs) up‐regulate two immune‐related genes, MMP14 and INHBA, which were used to establish the prognostic model. Knock‐down of MMP14 and INHBA inhibited invasion of pancreatic cancer.


| INTRODUC TI ON
Pancreatic cancer is one of the most lethal solid organ tumour types. As early-stage pancreatic cancer patients have mild or no clinical symptom, most pancreatic cancer patients are diagnosed at an advanced stage of cancer after obvious symptoms appear or during physical examination and hardly undergo radical resection.
Pancreatic cancer patients always result in poor prognosis, with a 5-year survival rate of <20%. 1 It is especially important to identify high-risk pancreatic cancer patients, predict the prognosis of cancer and adjust the treatment option accordingly at the first opportunity.
With the rapid development and cross-pollination of immunology and oncology in recent decades, a large number of studies have shown that the immune response is closely related to the progression of tumours. One sign of cancer is that tumour cells would adopt the ability to evade the host immune response and maximize their continuous growth potential. 2 Under normal conditions, tumour cell antigens stimulate an immune response that can lead to the elimination of cancer cells. 3 In some cases, this response might be inhibited, leading to a suitable immune microenvironment 4 where inflammatory cells are believed to promote development and tumour formation. 5,6 Macrophages are the major population of tumour-associated inflammatory cells and divided into two different phenotypes: macrophage M1 and macrophage M2. 7 Macrophage M1 cells have anti-tumour behaviours, while macrophage M2 cells can promote tumour growth, angiogenesis, invasion and metastasis. 8,9 A previous study illustrated that breast cancer cells separated and progressed into an epithelial-mesenchymal-transition (EMT) morphology in the presence of macrophage M2-conditioned medium (MCM). 9 Currently, classical TMN stage is the main method to determine the prognosis of cancer; however, due to the high heterogeneity of pancreatic cancer, patients with pancreatic cancer at the same stage often behave differently in terms of tumour recurrence and response to antitumour therapy. Therefore, there is an urgent need to find more reliable parameters to guide prognostic stratification and personalized treatment. With the rapid development of gene expression profiling technology, represented by second-generation sequencing, we are able to predict patient prognosis by analysing gene expression profiles to further screen for molecular markers with different characteristics. However, since the development of pancreatic cancer is usually a complex process involving multiple genes, a single molecular marker is not sufficient to predict the prognosis of pancreatic cancer patients and guide individual-

| Acquisition of data for pancreatic normal samples and tumour samples
The Cancer Genome Atlas (TCGA) data sets and The Genotype-Tissue Expression (GTEx) data sets were obtained from the UCSC Xena resource (http://xena.ucsc.edu/). 10 For gene expression, the TCGA RNA-Seq (polyA + Illumina HiSeq) data were downloaded as log2(norm_count + 1) values, and the GTEx RNA-Seq data were downloaded as log2(norm_count + 0.001) but was transformed to log2(norm_count + 1) values. Clinical parameters, including patient age, gender, clinic stage and pathological grade, were also acquired. All statistical analyses and bioinformatics analysis were performed using R software (version3.6.0) and its programme package.

| Identification of mRNA expressed differentially between pancreatic cancer samples and normal samples
In order to identify differentially expressed genes (DEGs) between pancreatic normal samples and PAAD, we used SVA R package migration ability of pancreatic cells was done by trans-well cell migration assay. A subcutaneous xenograft tumour model of human pancreatic cancer in nude mice was constructed. In conclusion, an immune-related gene prognostic model was constructed,   patients with high-risk scores have poorer survival status, M2-phenotype tumourassociated macrophages (TAMs) up-regulate two immune-related genes, MMP14 and   INHBA, which were used to establish the prognostic model. Knock-down of MMP14 and INHBA inhibited invasion of pancreatic cancer.

| Pathway enrichment and visualization
The gene set enrichment analysis (GSEA) was used to figure out if genes may be statistically critical in the studied gene set. 13 To investigate the potential function of DEGs in pancreatic samples, we performed pathway enrichment analysis based on GO Biological Process Ontology gene sets with GSEA, to explore the immune function between normal samples and tumour samples and to identify whether the corresponding immune genes were significantly different. Metascape (http://metas cape.org) is an online tool that was used for gene enrichment analysis while we tried to compare the different biological pathway in high-and low-risk group of patients, and it also found relevant signalling pathways. 14

| Establishment an immune prognostic model by estimating risk score
Risk scores for each patient were estimated according to genes expression that was multiplied a linear combination of regression coefficient that was acquired from the multivariate Cox regression.
Patients then were assigned to high-risk or low-risk subgroups based on the most suitable risk score, which were called the cut-off value; here, the cut-off value was 2.62. To validate the risk score model, a Kaplan-Meier OS curve for high-risk and low-risk groups was plotted and estimated the area under the curve (AUC) of the receiver that operates the characteristic (ROC) curve.

| Validation for prognostic risk prediction and Nomogram development
The nomogram that integrated a variety of immune prognostic model and clinical risk factors was assembled. To assess whether risk score can be used independently to predict prognosis of patients with pancreatic cancer, we performed univariate and multivariate Cox regression analyses on risk score and other clinical pathological factors (eg age and gender), and factors meeting HR >1 or <1

| Cell culture and gene knock-down techniques
The human monocytic leukaemia cell line THP-1 and the human pancreatic cancer cell line PANC1 were purchased from the American Type Culture Collection (ATCC). THP-1 cells were cultured in PRMI-1640 medium (Gibco). PANC1 cells were cultured in DMEM medium (Gibco). Medium was added with 10% FBS (foetal calf serum; Gibco), 100 mg/ml streptomycin and 100 U/ml penicillin. Cells were cultured at 37°C in a 5% CO 2 incubator with a humidified atmosphere. The cells were passaged every 2-3 days to sustain growth.

| Western blot analysis
Total protein was extracted using a lysis buffer (Cell lysis buffer for Western and IP; P0013; Beyotime). The protein extracts were separated by 10% SDS-polyacrylamide gel electrophoresis (PAGE) and were transferred to nitrocellulose filter (NC) membranes. After 2 h blocking in 5% non-fat milk, the membranes were incubated overnight at 4°C with rabbit anti-MMP14 antibody (13130S; Cell Signaling Technology), rabbit anti-INHBA antibody (ab233083, Abcam) or rabbit anti-GAPDH antibody (Santa Cruz Biotechnologies). The protein expression was determined using horseradish peroxidaseconjugated antibodies followed by enhanced chemiluminescence (ECL; Millipore) detection. GAPDH was used as internal control.

| Stromal landscape and immune landscape are different between pancreatic normal sample and PAAD
The flowchart of this study is shown in Figure 1

| A total 4620 genes were differentially expressed between normal pancreatic samples and PAAD
The genome wide mRNA expression profiles in pancreatic normal samples and tumour samples were combined from GTEx and TCGA databases (including 171 pancreatic normal samples and 177 tumour samples). |log2 FC (fold-change) | >1 and p < 0.05 were the cut-off criteria. Compared with normal pancreatic samples, we identified 2304 and 2316 genes that were significantly up-and downregulated respectively. The distribution of differentially expressed genes (DEGs) between normal samples and pancreatic cancer samples is shown by volcano plot ( Figure 2B).

| A weakened immune phenotype in pancreatic cancer
Considering the different immune landscape between normal pancreatic samples and PAAD samples, we evaluated the distinct im-  (Table S1). In these 9-biological immune-related processes, 487 immune-related genes were obtained from the 4,620 DEGs (Table S2). To further verify the association between these 487 immune-related genes and the prognosis of patients with pancreatic cancer, univariate Cox regression was used to analyse the 487 immune-related genes in the TCGA data sets, the results showed that 60 of the immune-related genes met the condition of HR >1 or <1 with p < 0.01, confirming that these 60 immune-related genes were indeed strongly associated with the prognosis of patients with pancreatic cancer (Table S3).

| An immune prognostic model based on MMP14 and INHBA expression was established and evaluated in TCGA data set
To confirm the robust of these 20 immune-related genes were indeed strongly associated with the prognosis of patients with pancreatic cancer, the GEO data set GSE62452 was utilized for further validation analysis. We first calculated stromal scores, immune scores and estimate scores for GEO data set GSE62452 ( Figure 3A).
Like the TCGA cohort, scores of cancerous microenvironments landscape were different from normal samples. We also analysed differentially expressed genes (DEGs) between pancreatic normal samples and tumorous samples in the GSE62452 cohort ( Figure 3B).
Next, we used a Venn plot to identify up-regulated genes in the GEO cohort and 20 prognosis-immune-related genes obtained from the TCGA cohort ( Figure 3C) to confirm which immune-related genes were strongly associated with prognosis; only two genes were iden- we divided patients with pancreatic cancer into high-and low-risk as stated by the optimum risk score cut-off value 2.62. Figure 3F presents the results of Kaplan-Meier survival analysis of the two groups of patients with pancreatic cancer in TCGA data set, and the F I G U R E 1 Workflow chart of this study prognosis of patients with pancreatic cancer was worse in the highrisk group compared with the low-risk group (p < 0.05). Figure 3D

| The immune prognostic model has the ability to accurately predict the prognosis of patients with pancreatic cancer in GEO data set
We utilized the GEO data set GSE62452 for validation and application of our immune prognostic model. Using the same method and cut-off value to calculate risk value, patients were categorized as high-risk or low-risk. The K-M survival curve was used to assess the ability of the risk score model to predict the prognosis of pancreatic cancer patients in GEO data set, and the results show that the overall survival of patients in the high-risk group was significantly shorter (p = 0.038) as shown in Figure 4C. Figure 4A sorts the patients into high-and low-risk groups based on the cut-off value (2.62) by ranking them from low-to high-risk scores and depicts the difference in survival

| A nomogram was developed and could facilitate the use the immune prognostic model
To investigate whether risk score could independently predict patient prognosis, we performed independent prognostic analysis on risk score and other clinical pathological factors in TCGA data set including age, gender, stage and grade. Univariate Cox regression analysis indicated that age (p < 0.01), grade (p < 0.05) and risk score (p < 0.01) calculated from the two immune-related genes were independent prognostic factors for OS. While multivariate Cox regression analysis indicated that age (p < 0.05) and risk score (p < 0.05) were independent prognostic factors for OS ( Figure 5A). In order to predict the prognosis of pancreatic cancer patients more conveniently and accurately, a nomogram was then constructed to predict 1 year, 3 years and 5 years OS based on two independent prognostic factors analysed by multivariate Cox regression including age and risk score. The first row is the score scale of individual factors, the second and the third rows are the age and risk score. Accordingly, the scores corresponding to the ages and risk scores are derived, and the total scores obtained are added and compared with the total score scale in the fourth row to derive the 1 year, 3 years and 5 years survival rates of patients and to visually analyse the prognosis of pancreatic patients ( Figure 5B). The nomogram was suitable to estimate the mortality showed by the calibration plots ( Figure 5C). Decision curve analysis showed that the nomogram demonstrated the most excellent net benefit for 1 year, 3 years and 5 years OS ( Figure 5D), further shows that the nomogram has high predictive consistency.

| The risk score could partially reflect the state of the immune microenvironment
The immune prognostic risk model is composed of immune-related genes, we hypothesized that the risk score could partially reflect the state of the immune microenvironment. First, CIBERSORTx was used to evaluate 22 immune cells proportions, Figure 6A (Figure 6K-N).  Figure 7E).

| Macrophage M2 could upregulate the prognostic genes MMP14 and INHBA in pancreatic cancer
To investigate the different clinicopathological features of the different risk groups of the pancreatic cancer, the pathway was identified by us underlying the risk score, differential expression analysis was performed and got 251 DEGs between the high-risk group and lowrisk group (|log FC| > 1, p < 0.05) ( Figure 8A). We used Metascape to annotate the function of DEGs and charted the top 20 significant biological processes, as shown in Figure 8B. GO analysis of the DEGs revealed that 'extracellular matrix organization' were enriched in biological processes and pathways, which might be correlated with pancreatic cancer's malignant progression. The 'extracellular matrix organization' is the most significant biological process, which reveals that in the high-risk subgroup patients, cancer cells may have more progression ability than the low-risk subgroup. Furthermore, only 1.1% patients were containing genetic alterations in MMP14 and 1.6% genetic alterations in INHBA ( Figure 8C) shown in cBioportal for TCGA database, demonstrated that these two genes contained only little genetic alteration. Moreover, the mRNA expression of MMP14 and INHBA was significantly up-regulated compared with non-tumour samples in the Oncomine Pei pancreas cohort ( Figure 8D). MMP14's protein expression was further explored by us in the human protein profiles and shown in Figure 8E. However, we did not find INHBA protein expression in the database.
We then used Metascape to annotate the function of DEGs through KEGG analyse, analysis of the DEGs revealed that 'ECM receptor interaction' was the most significant pathway, which might be highly correlated with malignant progression of pancreatic cancer.

| DISCUSS ION
Although the treatment of pancreatic cancer has made great progress, the incidence and mortality rate of pancreatic cancer are still in the forefront and on the rise. 17,18 The prognosis of pancreatic cancer is relatively poor worldwide, which leads to a fairly similar morbidity and mortality rate. 19 The number of pancreatic cancer patients in China is a major economic and social burden. 20 Tumours are highly heterogeneous, which prose a great challenge for patient prognosis prediction and individualized treatment.  cells may support cancer progression and tumour growth. 24,25 Matrix metalloproteinases (MMP) are a category of enzymes mediating alterations in the TME that occurs throughout progression and tumour development. 26,27 In humans, 6 distinct membrane-type MMPs (MT-MMPs) have been discovered. MMP14 was first depicted by Sato et al. 28,29 as a transmembrane protein activating pro-MMP2 to induce tumour cell invasion. 30,31 In several cancer types, MMP14 is up-regulated and contributes to advancing inflammation, angiogenesis, metastasis and cancer cell invasion. MMP14's potential clinical relevance in pancreatic cancer has been addressed applying patient-derived tumour substance. 32-35 A rather consistent picture was shown by studies with respect to MMP14 overexpression in tumours compared with control sections. 36 Tumour-associated macrophages M2 expressing MMP14 are involved in complicated matrix remodelling, as demonstrated in a colorectal cancer orthotopic mouse model. 37 We discovered that MMP14 expression is higher in pancreatic cancer samples compared with corresponding normal samples, and the 38 MMP14 expression is correlated with macrophage M2.
INHBA, which encodes a member of the TGFβ (transforming growth factor-beta) superfamily of proteins, has been discovered to play a vital function in distinct kinds of cancer. 39

CO N FLI C T O F I NTE R E S T
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.