Comprehensive Analysis of PPMs in Pancreatic Adenocarcinoma Indicates the Value of PPM1K in the Tumor Microenvironment

Simple Summary Pancreatic adenocarcinoma is a devastating disease, with an extremely poor survival rate worldwide. Its poor responsiveness to chemotherapy and immunotherapy has a bearing on the unique tumor microenvironment. Our study demonstrates that Metal-dependent protein phosphatases (PPMs) be implicated in cell–cell adhesion and immune cell infiltration in pancreatic cancer. Among these, PPM1K was downregulated in the tissue and peripheral blood of pancreatic adenocarcinoma patients, negatively related to PD-L1 expression and poor prognosis. The knockdown of PPM1K markedly promoted the proliferation and migration of pancreatic cancer cells, confirming its role in tumor suppressor activity in pancreatic adenocarcinoma. This study reveals the potential clinical utility of PPM1K in tumor immunotherapy and brings about novel insights into the prognostic value of PPM1K in pancreatic adenocarcinoma. Abstract Early metastasis and resistance to traditional therapy are responsible for the poor prognosis of pancreatic adenocarcinoma patients. Metal-dependent protein phosphatases (PPMs) have been proven to play a crucial role in the initiation and progression of various tumors. Nevertheless, the expression and function of distinct PPMs in pancreatic adenocarcinoma have not been fully elucidated. In this study, we investigated the mRNA expression level, prognostic value, and the relationship between the expression of PPMs and the tumor microenvironment in pancreatic adenocarcinoma using Oncomine, TCGA and GTEx, GEO, Kaplan–Meier plotter, STRING, GeneMANIA, and HPA databases and R packages. GO and KEGG analysis revealed that PPMs and their differential co-expression genes are attributed to cell–cell adhesion and immune cell infiltration. Among these, PPM1K was downregulated in the tissue and peripheral blood of PAAD patients, whose expression level was negatively related to poor prognosis. Further to this, PPM1K was found to play a role in the epithelial–mesenchymal transition and immune infiltration. ROC curves showed that PPM1K had a good predictive value for pancreatic adenocarcinoma. The knockdown of PPM1K markedly promoted the proliferation and migration of pancreatic cancer cells, confirming its role in tumor suppressor activity in PAAD. This study demonstrates the potential clinical utility of PPM1K in tumor immunotherapy and brings about novel insights into the prognostic value of PPM1K in pancreatic adenocarcinoma.


Background
Metal-dependent protein phosphatases (PPMs), also known as the type 2C family of protein phosphatases (PP2Cs), belong to protein Ser/Thr phosphatases (PSPs). Members of the PPMs family possess an N-terminal catalytic lobe and a C-terminal 90-residue-lobe,

mRNA and Protein Expression of PPMs in PAAD
The transcriptional levels of different PPM family members in various cancer tissues were analyzed in the Oncomine gene expression array dataset (www.oncomine.org, accessed on 18 January 2021) [23]; a p-value < 0.05, fold change of 1.5, and a gene rank in the top 10% were set as the significant thresholds. The mRNA expression profiles and clinical data of pancreatic adenocarcinoma (PAAD) patients (n = 179) and normal controls (n = 4) in TCGA (the Cancer Genome Atlas, https://portal.gdc.cancer.gov/, accessed on 31 December 2020) [24], and the transcription expression of PPMs in normal pancreas tissue (n = 167) in the GTEx project (Genotype-Tissue Expression) [25] were obtained from UCSC XENA (https://xenabrowser.net/datapages/, accessed on 31 December 2020). In addition, Gene expression profiling data sets (GSE28735, GSE15471, GSE16515, and GSE71989) involved 130 PAAD patient tissues and 104 normal controls were downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/gds, accessed on 18 January 2021) [26]. Moreover, the expression profile of peripheral blood (GSE 74629 involved 36 PAAD patients and 14 healthy controls) was also analyzed to evaluate the application of the target PPMs. Finally, GSE23952 (dataset inclusion of pancreatic cancer cell line PANC-1 with (n = 3) or without (n = 3) TGF-β for EMT induction) was used to verify the role of PPMs in EMT. The protein expression data of PPMs in PAAD and normal tissues was explored using HPA (Human Protein Atlas) (https://www.proteinatlas.org/, accessed on 29 June 2021) [27].

Prognostic Value of PPMs in PAAD
The Kaplan-Meier plotter database (http://kmplot.com/analysis/, accessed on 29 June 2021) was used to conduct a survival and prognosis analysis of PPMs in PAAD (n = 177) [28]. Patient samples were divided into high and low-expression groups due to the median mRNA expression of PPMs. Univariate and multivariate Cox regression analyses of PPMs were also conducted to assess the impact of expression levels of PPMs on pathologic features using the "survival" package [24]. ROC (receiver operating characteristic) curve analysis was conducted to describe the predictive value of PPMs using the "pROC" and "ggplot2" packages [29,30]. The RNA expression data and corresponding clinical characteristics of PAAD patients were collected from the TCGA and GTEx datasets (normal = 171, tumor = 179) and the GSE21501 (including 102 patients) dataset. Survival analysis was performed using Cox regression models with "survival" and "survminer" packages [24].

Functional Enrichment Analysis
To further explore the underlying biological function of PPMs molecules in PAAD, functional enrichment analysis was applied based on the transcription data from TCGA. Differential expression genes related to PPMs were analyzed with the "DeSeq2" package [31]. All PPMs expressions were divided into two groups (the high expression and the low-expression group). The differentially expressed genes associated with PPMs were distinguished by |log 2 (FoldChange)| > 1 and adjusted p < 0.05. The "ClusterProfiler" package and "ggplot2" were used to analyze Gene Ontology (GO) and the Kyoto encyclopedia of genes and genomes (KEGG) [32].

Correlations between PPMs Expression and Tumor Environment
The "GSVA" and "ESTIMATE" packages were employed to analyze the correlation between PPMs and immune cell infiltrates using the data from TCGA. p < 0.05 was defined as the threshold of significant difference [33][34][35]. The co-expression relationship between PPMs and target genes (including epithelial-mesenchymal-associated genes and immune cell biomarkers and immune checkpoints) was analyzed using TCGA data and visualized using the "ggplot2" package with Spearman's correlation. The protein-protein interactions (PPI) of PPMs binding proteins were investigated using the GeneMania [36] (http:// genemania.org/, accessed on 30 June 2021) and STRING websites [37] (https://string-db. org/, accessed on 30 June 2021) and were visualized through Cytoscape (3.8.2) [38].
The human siRNAs targeting PPM1K (siRNA-1, siRNA-2, and siRNA-3) and negative control (NC) were purchased from Kidan Biosciences (Guangzhou, China). For stable transfection, cells were plated in six-well plates 24 h in advance, then transfected with siRNAs using lipofectamine 2000 reagents (Invitrogen, Waltham, MA, USA) according to the manufacturer's instructions. Two days after transfection, these cells were harvested for the following experiments. The siRNA sequences are listed in Supplementary Table S1.

Quantitative PCR (qPCR) Detection
Total RNA was isolated using an RNA Purification Kit (EZBioscience, Roseville, MN, USA) following the manufacturer's protocol. PrimeScript TM RT Master Mix (Takara, Tokyo, Japan, RR036A) was used to carry out the reverse transcription. qPCR was performed with a ChamQ SYBR qPCR Master Mix (Vazyme, CN). The cycling conditions were 95 • C for 30 s, 95 • C for 10 s, and 60 • C for 30 s for 40 cycles using CFX Connect (Bio-Rad, Hercules, CA, USA). β-actin was used as an endogenous control. The primer sequences are listed in Supplementary Table S1.

Cell Proliferation Detection
The Cell Counting Kit-8 (CCK-8; APExBIO, Houston, TX, USA) assay was used to estimate the cell growth capacity. A total of 2 × 10 3 cells were seeded into each well of 96-well plates. Subsequently, the culture medium and CCK-8 were mixed at a ratio of 10:1, and 100 µL of the resulting mixture was added to each well at different time points (24,48,72, and 96 h) after seeding. Absorbance was measured at 450 nm with a Synergy™ H1m Microplate Reader (BioTek, Winooski, VT, USA).

Transwell Assays
Transwell inserts (8-µm pore, Costar, Washington, DC, USA) were performed to evaluate cell migration ability. DMEM with 10% FBS (600 µL in total) was added to the lower chamber. After transfection, cells were re-suspended, harvested, and seeded in serum-free DMEM. A total of 1 × 10 5 Cells in 200 µL serum-free medium were added to the upper chamber. After incubation for 24 h at 37 • C, the non-migrated cells were removed from the upper surface of the membrane. Cells on the bottom surface of the membrane were fixed with paraformaldehyde and stained with crystal violet. Migration potential was assessed by calculating the number of stained cell nuclei from three random fields using a NI-U (Nikon, Tokyo, Japan) microscope system.

Statistical Analysis
Data are presented as mean ± standard deviation (SD). GraphPad Prism 8 software was used to assess the statistical significance between groups. Student's t-test and the Wilcoxon rank-sum test were used to compare the transcription levels of PPMs and PD-L1 in PAAD. Fisher's exact and Chi-square tests were used for the analysis of contingency tables. Survival analysis of patients was conducted using Kaplan-Meier curves (Log-rank test and Cox regression). All R packages were deployed using R software version (v3.6.3), and p < 0.05 was defined as statistical significance (ns, p ≥ 0.05; *, p < 0.05; **, p < 0.01; ***, p < 0.001).

Transcriptional Levels of PPMs and Clinicopathological Parameters of Patients in PAAD
To explore the mRNA expression of distinct PPMs members in various cancer and normal tissue specimens, we analyzed the Oncomine database ( Figure 1A). Furthermore, data from TCGA and GTEx were also analyzed ( Figure 1B). Finally, we analyzed the transcriptional levels of PPMs in PAAD according to the data from the GEO datasets ( Figure 2A, Supplementary Figures S1 and S2). PPM1K was found to be downregulated in tumor tissue in Oncomine, TCGA, and GSE 16515 ( Figure 2A). Moreover, we analyzed the mRNA expression level of PPMs in peripheral blood and found that PPM1K and ILKAP decreased in patients with PAAD ( Figure 2B).    Figure S3) was carried out. The expression of PPM1K in T 1+2 was markedly higher than that in the T 3+4 subgroup in TCGA; however, there were no significant differences in GSE21501. Then, we conducted univariate and multivariate Cox regression analyses to assess the prognostic value of PPMs in PAAD. Univariate Cox regression analyses suggested that a higher expression of PPM1E/K, PHLPP2, and ILKAP had a protective effect on overall survival (OS) in patients with PAAD. However, patients with higher PDP1 had a poorer OS rate. Then, variables with p-value ≤ 0.1 were included for multivariate regression analysis, and the conclusion that PHLPP2 and ILKAP favored PAAD patients' survival was confirmed by multivariate Cox regression analysis (p < 0.05, Figure 3E,F). The protein level of these PPMs in PAAD and normal tissues were compared using immunohistochemistry (IHC) staining via the HPA database. Consistent with transcription expression, protein levels of PPM1E/K were lower and PPM1G/L, PDP1, and PPTC7 were higher in PAAD ( Figure 3G).  Figure S3) was carried out. The expression of PPM1K in T1+2 was markedly higher than that in the T3+4 subgroup in TCGA; however, there were no significant differences in GSE21501. Then, we conducted univariate and multivariate Cox regression analyses to assess the prognostic value of PPMs in PAAD. Univariate Cox regression analyses suggested that a higher expression of PPM1E/K, PHLPP2, and ILKAP had a protective effect on overall survival (OS) in patients with PAAD. However, patients with higher PDP1 had a poorer OS rate. Then, variables with p-value ≤ 0.1 were included for multivariate regression analysis, and the conclusion that PHLPP2 and ILKAP favored PAAD patients' survival was confirmed by multivariate Cox regression analysis (p < 0.05, Figure  3E,F). The protein level of these PPMs in PAAD and normal tissues were compared using immunohistochemistry (IHC) staining via the HPA database. Consistent with transcription expression, protein levels of PPM1E/K were lower and PPM1G/L, PDP1, and PPTC7 were higher in PAAD ( Figure 3G).

Prognostic Value of PPMs in PAAD Patients
Associations between distinct PPMs expression and patients' clinical outcomes were analyzed using the Kaplan-Meier plotter database and the TCGA and GEO datasets. Kaplan-Meier plotter analysis showed that patients with higher PPM1K expression levels were significantly associated with better OS (Figure 4A), and similar results could be concluded from T 1+2 , N 1 , Stage I, and G1 subgroup analysis (Supplementary Figure S4, the sample number in T1, T4, Stage III, Stage IV, and G4 subgroups was too low for meaningful analysis). Furthermore, we conducted an analysis based on data from GSE21501. Similar to the result in TCGA, a higher PPM1K expression was associated with a better OS (Supplementary Figure S5). As presented later, PPM1K might play a role in tumor immune infiltration. We further conducted a prognostic analysis, using the Kaplan-Meier plotter database, on the correlation between PMM1K expression and immune cell infiltration. It was shown that a higher PPM1K expression correlated to the prolonged OS in patients with enriched CD4 + T cells, mesenchymal stem cells, Treg cells, and patients with decreased macrophages, and Th1 cells (Supplementary Figure S4D). In summary, PPM1K might be a novel and promising biomarker in PAAD.

9
PPM1E, PDP1, PHLPP2, and ILKAP were correlated to PAAD patients' overall survival (OS). (F) Multivariate Cox regression analysis showed that higher PHLPP2 and ILKAP favored the OS of PAAD patients. (G) PPM1E/K were lower and PPM1G/L, PDP1, and PPTC7 were higher in PAAD protein levels using IHC staining via the HPA database (Protein levels of other PPMs showed no significant differences and are not presented).

Prognostic Value of PPMs in PAAD Patients
Associations between distinct PPMs expression and patients' clinical outcomes were analyzed using the Kaplan-Meier plotter database and the TCGA and GEO datasets. Kaplan-Meier plotter analysis showed that patients with higher PPM1K expression levels were significantly associated with better OS (Figure 4A), and similar results could be concluded from T1+2, N1, Stage I, and G1 subgroup analysis (Supplementary Figure S4, the sample number in T1, T4, Stage III, Stage IV, and G4 subgroups was too low for meaningful analysis). Furthermore, we conducted an analysis based on data from GSE21501. Similar to the result in TCGA, a higher PPM1K expression was associated with a better OS (Supplementary Figure S5). As presented later, PPM1K might play a role in tumor immune infiltration. We further conducted a prognostic analysis, using the Kaplan-Meier plotter database, on the correlation between PMM1K expression and immune cell infiltration. It was shown that a higher PPM1K expression correlated to the prolonged OS in patients with enriched CD4 + T cells, mesenchymal stem cells, Treg cells, and patients with decreased macrophages, and Th1 cells (Supplementary Figure S4D). In summary, PPM1K might be a novel and promising biomarker in PAAD.

Co-Expression, PPI, and Functional Enrichment Analysis of PPMs in PAAD Patients
GO and KEGG functional enrichment analysis revealed that differential expression PPMs were associated with cell adhesion molecules, cell growth and proliferation, immune infiltration and response, the cAMP signal pathway, the JAK-STAT signal pathway, and the PI3k-AKT signal pathway (Supplementary Figure S6). Protein-protein interaction analysis of PPMs was performed using STRING and GeneMANIA. As presented in Supplementary Figure S7A, the PPI network resulting from STRING showed a correlation with amino  Figure S7B). Therefore, we conducted a correlation analysis between PPMs and EMT-associated genes and immune infiltration. It revealed that there were close relationships between PPMs and EMT-associated genes in PAAD, especially PPM1F/K/M and PDP1 ( Figure 5A, Supplementary Figure S7C). The association was confirmed by the analysis of GSE23952 (Supplementary Figure S8). After TGF-β treatment for EMT induction, PPM1K expression decreased in pancreatic cancer cell line PANC-1, which suggested that PPM1K was closely related to EMT.

Co-Expression, PPI, and Functional Enrichment Analysis of PPMs in PAAD Patients
GO and KEGG functional enrichment analysis revealed that differential expression PPMs were associated with cell adhesion molecules, cell growth and proliferation, immune infiltration and response, the cAMP signal pathway, the JAK-STAT signal pathway, and the PI3k-AKT signal pathway (Supplementary Figure S6). Protein-protein interaction analysis of PPMs was performed using STRING and GeneMANIA. As presented in Supplementary Figure S7A, the PPI network resulting from STRING showed a correlation with amino acid metabolism and glycolysis. GeneMANIA analysis revealed the relationship between protein dephosphorylation and the integral component of the postsynaptic membrane (Supplementary Figure S7B). Therefore, we conducted a correlation analysis between PPMs and EMT-associated genes and immune infiltration. It revealed that there were close relationships between PPMs and EMT-associated genes in PAAD, especially PPM1F/K/M and PDP1 ( Figure 5A, Supplementary Figure S7C). The association was confirmed by the analysis of GSE23952 (Supplementary Figure S8). After TGF-β treatment for EMT induction, PPM1K expression decreased in pancreatic cancer cell line PANC-1, which suggested that PPM1K was closely related to EMT.

Immune Cell Infiltration of PPMs in PAAD
As is known, immune cell infiltration has a close bearing on the progression and therapy outcome of cancers. The ESTIMATE algorithm based on tumor immune and stromal features is believed to be a good prediction model for patient prognosis. PPMs were closely correlated to immune activation and response based on previous functional enrichment analyses (Supplementary Figure S6). Therefore, we investigated the correlation between PPMs family members and immune infiltrates (Supplementary Table S2), and it was suggested that PPM1K expression was positively associated with the infiltration of B cells, Mast cells, T cells including CD8 + cytotoxic cells, T helper cells, T follicular helper cells (TFH), and Th1 cells ( Figure 5C). We also found that the immune score and stromal score were higher in PAAD tissue with higher PPM1K expression (Supplementary Figure  S4C). The correlation between PPM1K and markers of immune cells was also analyzed. Consistent with the cell infiltration analysis, PPM1K had a close relationship with immune cell markers and immune checkpoints ( Figure 5B,D).

Predictive Value of PPMs in Clinical Applications
In order to explore the clinical utility of PPMs, ROC curves were used to assess the predictive value of PPMs in PAAD. As was shown in Figure 5E, most PPMs had good predictive values in PAAD tissue. Moreover, PPM1K was differentially expressed in peripheral blood between PAAD patients and normal controls, which has promising value for clinical applications ( Figure 5F).

PPM1K Acts as a Tumor Suppressor and Participates in PD-L1 Regulation in PAAD
Owing to its significantly differential expression and favorable prognostic value in PAAD, we selected PPM1K uniquely in our study. According to the analysis mentioned above, PPM1K expression is concerned with immune cell infiltrates and immune checkpoints. For further validation, we explored the expression of PPM1K in pancreatic cancer cell lines and the immortal pancreatic duct cell line hTERT-HPNE. It was suggested that, compared to hTERT-HPNE, PPM1K expression was downregulated in pancreatic cancer cell lines ( Figure 6A). PPM1K expression was relatively high in PANC-1 and HPAF-II cell lines. Thus, we conducted knockdown experiments in PANC-1 and HPAF-II cells. SiRNA-3 was selected due to its high knockdown efficiency ( Figure 6B). As described previously, PPM1K had an impact on immune infiltration. It is shown that the knockdown of PPM1K can upregulate PD-L1 expression ( Figure 6C). Furthermore, PPM1K knockdown can enhance the capacity of proliferation and migration in PAAD cells, indicating that PPM1K acts as a tumor suppressor in PAAD ( Figure 6D,E).

Discussion
Our study shows that most PPMs differentially express between normal tissue and tumors. However, the expression status of PPMs in different databases does not coincide exactly. In this study, we first explored the expression and prognosis of PPMs and their correlation with EMT and immune infiltrates in PAAD. Survival analysis revealed that PPM1D/E/G/H/J/K/L, PDP1, PP2D1, PHLPP2, PPTC7, and ILKAP are associated with prognosis in TCGA. Among these, PPM1E/G/K in GSE21501 was consistent with the results in TCGA. PPM1E/G expression is higher in tumor tissue in TCGA+GTEx but lower in Oncomine. While PPM1K is downregulated in most of the datasets including TCGA+GTEx, Oncomine, GSE 16515, and GSE 74629 (PPM1K is not detected in GSE 28735 and is lower

Discussion
Our study shows that most PPMs differentially express between normal tissue and tumors. However, the expression status of PPMs in different databases does not coincide exactly. In this study, we first explored the expression and prognosis of PPMs and their correlation with EMT and immune infiltrates in PAAD. Survival analysis revealed that PPM1D/E/G/H/J/K/L, PDP1, PP2D1, PHLPP2, PPTC7, and ILKAP are associated with prognosis in TCGA. Among these, PPM1E/G/K in GSE21501 was consistent with the results in TCGA. PPM1E/G expression is higher in tumor tissue in TCGA+GTEx but lower in Oncomine. While PPM1K is downregulated in most of the datasets including TCGA+GTEx, Oncomine, GSE 16515, and GSE 74629 (PPM1K is not detected in GSE 28735 and is lower in tumor tissue but not significantly different in GSE15471 and GSE 91989). We also explored the mRNA expression of PPM1K, and the results were aligned with the conclusions drawn from the public databases that PPM1K expression is lower in cancer cell lines than in pancreatic duct cell lines. Hence, we chose PPM1K as the target molecular in our study. Further knockdown experiments in vitro suggest that the down-regulation of PPM1K can promote the proliferation and migration of pancreatic cancer cells.
The functional enrichment and co-expression analyses in our study showed that PPM1A/B/D/F/K/M, PDP1/2, and PHLPP2 had a strong correlation with EMT (only PPM1K is presented). PPM1K expression was positively correlated with mesenchymal-associated genes [21,22,39] (SNAI1/2, Twist1/2, ZEB1/2, VIM, and CDH2) and negatively with the epithelial-associated gene CDH1. However, PPM1K expression decreases in PANC-1 cells after TGF-β treatment in GSE23952. This inconsistent result may be due to tumor heterogeneity. In addition, as has been presented before, PPM1K upregulation is positively correlated with the PAAD stromal score. Positive relationships between PPM1K and mesenchymal-associated genes may be due to stromal and immune components in tissue.
Moreover, our study demonstrates that PPM1K participates in tumor immune infiltrates. It is suggested that PAAD tissue with a high PPM1K expression has a higher immune score, showing that PPM1K expression is associated with immune cell infiltration. Correlation analysis reveals that PPM1K is positively correlated with B cells, Th1 cells, and CD8 + T cells, which always predicts a good prognosis in PAAD [20]. Further, the knockdown of PPM1K can increase the expression of PD-L1 in pancreatic cancer cells. It is well known that EMT and enriched mesenchymal cells always give clues to tumor progression and bad prognoses [22]. Moreover, enriched Treg and decreased Th1 and macrophage infiltration may contribute to tumor immune suppression and lead to bad prognoses [20,40]. However, KM plot analysis showed that the overexpression of PPM1K might reverse this condition and prolong overall survival. EMT status is found to be concerned with checkpoint therapy response [22]. Therefore, combination therapy both targeting EMT and checkpoints may be more efficient. As previously indicated, PPM1K may play a valuable role in both EMT and the tumor immune microenvironment. The ROC model showed that PPM1K has a good predictive value not only in tissue but also in peripheral blood. The area under the ROC curve (AUC) for PAAD diagnosis is 0.823 in peripheral blood, which is more suitable for a wide range of clinical applications.

Conclusions
Though inhibitors of some PPMs including PPM1A/B/D/E/F and PHLPP1/2 are already in pharmaceutical research and development, there is still a long way to go before a full understanding of PPMs is achieved. Our study provides a guiding light in the search for novel tumor biomarkers and therapeutic targets. For its correlation with the tumor microenvironment and differential expression in both tissue and peripheral blood between tumor and normal control, PPM1K may serve as a prognostic marker or potential target for PAAD therapy.

Conflicts of Interest:
The authors declare that they have no competing interest.

EMT epithelial-mesenchymal transition GTEx
Genotype-Tissue Expression IHC immunohistochemistry JNK c-Jun N-terminal kinase KEGG Kyoto Encyclopedia of Genes and Genomes PAAD pancreatic adenocarcinoma PP2Cs type 2C family of protein phosphatases PPMs metal-dependent protein phosphatases PSPs protein Ser/Thr phosphatases ROC receiver operating characteristic TCGA The Cancer Genome Atlas TGF-β transforming growth factor-β TME tumor microenvironment HPA Human Protein Atlas