Opposite Roles of Tumor Cell Proliferation and Immune Cell Infiltration in Postoperative Liver Metastasis of PDAC

Background Recurrence of liver metastasis after pancreatectomy is often a predictor of poor prognosis. Comprehensive genomic analysis may contribute to a better understanding of the molecular mechanisms of postoperative liver metastasis and provide new therapeutic targets. Methods A total of 67 patients from The Cancer Genome Atlas (TCGA) were included in this study. We analyzed differentially expressed genes (DEGs) by R package “DESeq2.” Weighted gene co-expression network analysis (WGCNA) was applied to investigate the key modules and hub genes. Immunohistochemistry was used to analyze tumor cell proliferation index and CD4+ T cells infiltration. Results Functional analysis of DEGs between the liver metastatic and recurrence-free groups was mainly concentrated in the immune response. The liver metastasis group had lower immune and stroma scores and a higher TP53 mutation rate. WGCNA showed that the genes in key modules related to disease-free survival (DFS) and overall survival (OS) were mainly enriched in the cell proliferation process and tumor immune response. Immunohistochemical analysis showed that the pancreatic cancer cells of patients with early postoperative liver metastasis had higher proliferative activity, while the infiltration of CD4+ T cells in tumor specimens was less. Conclusion Our study suggested that increased immune cell infiltration (especially CD4+ T cells) and tumor cell proliferation may play an opposite role in liver metastasis recurrence after pancreatic cancer.


INTRODUCTION
Pancreatic ductal adenocarcinoma (PDAC) is the fourth leading cause of cancer-related death in the United States (Siegel et al., 2021). Surgical resection is the only curative treatment for PDAC. However, the poor postoperative prognosis is usually common for most patients because of the locoregional recurrence and distant metastasis (Mizrahi et al., 2020). Several studies have investigated recurrence after curative resection of PDAC. And more than half of patients with pancreatic cancer develop hepatic recurrences after surgery, according to the literature (Sperti et al., 1997;Van den Broeck et al., 2009). In another study, which evaluated 692 patients with comprehensive and detailed followup data, more than a quarter of patients had only liver metastasis after surgery, with a median duration of 6.9 months. Meanwhile, the author also found that the lung-only recurrence occurred later, with a median time of 18.6 months (Groot et al., 2018b). Similar results were confirmed in another study, suggesting that early liver metastasis after surgery is a risk factor for shorter survival in pancreatic cancer patients (Kolbeinsson et al., 2020). Thus, it indicated that postoperative liver metastasis (especially liver-only recurrence) is one of the reasons for the poor prognosis of pancreatic cancer patients after surgery. Moreover, the molecular phenotype and biological behavior of tumors may vary significantly among different types of metastases.
As a result, accurate prediction of the risk of postoperative liver metastasis and identification of therapeutic targets are of great significance for guiding the postoperative treatment of pancreatic cancer, and gene-level analysis can be the most profound interpretation of the molecular characteristics associated with liver metastasis. Therefore, the purpose of this study was to investigate the possible molecular mechanisms of postoperative liver metastasis of pancreatic cancer to find potential molecular targets.

Acquisition of Datasets
Clinical information (n = 67) of PDAC samples and corresponding mRNA sequencing data (n = 65) were downloaded from The Cancer Genome Atlas (TCGA). 1 In addition, normalized reverse phase protein array (RPPA) data were obtained from cBioPortal, 2 which measured 218 proteins in 49 out of 67 PDAC samples (Figure 1).

Patients and Samples
Twenty-seven patients who underwent pancreatectomy for PDAC in our center between June 2015 and December 2016 were considered suitable for this study. By the last followup, 13 patients had liver metastases recurrence (disease-free survival (DFS) of 10 patients was less than 12 months, and the other 3 patients had a DFS of 13, 15, and 22 months, respectively), while 14 patients had no recurrence (follow-up time of these patients was no less than 36 months). There were no significant differences in age, sex, smoking, alcohol consumption, diabetes mellitus, TNM stage, tumor weight, positive lymph nodes, maximum tumor dimension, pathological stage, and postoperative chemotherapy between the two groups ( Table 1). Surgically resected specimens were fixed with 10% formalin and cut into 5 mm thick continuous sections for subsequent studies. This study was approved by the Ethics Committee of the First Affiliated Hospital of Nanjing Medical University.

Immunohistochemical Analysis
Paraffin-embedded tissue from pancreatic cancer patients was processed using standard methodology. Tissue samples were deparaffinized and rehydrated, and antigen retrieval was performed by citrate buffer (pH 6) for 15 min at 95 • C. Slides were incubated overnight with primary antibody anti-PCNA (1:200,

Identification of Differentially Expressed Genes
The R package "DESeq2" was used to identify the differentially expressed genes (DEGs). Adjusted p-value <0.05 and | Log2fold change| >1 were set as cut-off criteria.

Construction of Weighted Gene Co-expression Networks
The mRNA sequencing data from 30 postoperative liver metastasis patients were used to construct the weighted gene co-expression network by WGCNA package in R software (version 4.0.2). First, we set the soft threshold as 6 (scale-free R 2 = 0.85) and establish an adjacent matrix according to the connectivity degree so that our genes can satisfy the precondition of scale-free network distribution. Then, the adjacency matrix was transformed into a topological overlap matrix (TOM), which was used to describe the similarity of gene expression, and 1-TOM represented the inter-gene heterogeneity. The TOM is clustered by dissimilarity between genes, and then the tree is cut into different modules by a dynamic tree cutting algorithm. Here, we set the cut height at 0.25 and the minimum module size to 30. Module eigengene (ME) of each module is calculated as the first principal component of a specific module.

Identification of Clinically Relevant Modules and Key Module Genes
We calculated the correlation between the clinical traits and the MEs through the Pearson correlation coefficient, and p < 0.05 was considered to be significant. Key module genes were identified in specific modules under the threshold of module membership | MM| >0.80 and gene significance | GS| > 0.2.

Functional Enrichment Analysis
We used the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis to reveal the key pathways associated with these DEGs. In addition, Gene Ontology (GO) terms enrichment analysis was used to identify the functions of these DEGs in terms of biological process (BP). The clusterProfiler package in R software conducted the analysis, and p < 0.05 was considered to have statistical significance. ggplot2 packages were utilized for the visualization of the identified DEGs.

Construction of Protein-Protein Interaction Networks
In this study, the STRING database 3 was applied to construct the PPI network of DEGs and hub genes identified in specific modules. Confidence score >0.9 was set as the minimum required interaction score. Using the cytoHubba plugin based on Cytoscape (version 3.8.0), the top 10 ranked genes were considered as hub genes in the PPI network.

Identification of Hub Genes and Construction of a Prognostic Signature
We used DFS to represent the time of postoperative recurrence of liver metastasis. A univariate Cox proportional hazards regression analysis was performed to investigate the association between the top 10 ranked genes and DFS and overall survival (OS). p < 0.05 was considered as prognostic genes. The prognostic genes were further applied to stepwise multivariate Cox proportional hazards regression analysis to predict the risk value of each patient. The patients were divided into highand low-risk groups, and the predictive value of these genes was estimated by mapping the receiver operating characteristic (ROC) curve. The sensitivity and specificity of the DFS risk model were evaluated by calculating the area under the ROC curve using the "timeROC" R package.

Statistical Analysis
The Student's t-test was applied to compare the diverse expressions of PCNA and numbers of CD4 + T cells in the recurrence-free and liver metastases recurrence groups in our specimens. Correlation between the immune infiltration level and the DFS and OS was analyzed by using a Pearson correlation test.
In addition, Kaplan-Meier curves were conducted to compare the DFS and OS between the hub genes based on the different gene expressions. The statistical analysis was performed using R software (version 4.0.2) and GraphPad Prism (version 8.4.0).

Patient Baseline Characteristics
A total of 67 patients were enrolled in this study, including 34 recurrence-free patients and 33 liver metastasis patients. There were no significant differences in age, sex, smoking, alcohol consumption, diabetes mellitus, TNM stage, tumor stage, and postoperative chemotherapy between the two groups. However, patients with postoperative liver metastasis had larger tumor maximum dimension trend (3.945 ± 1.012 vs. 3.448 ± 1.065; p = 0.0689), higher histopathological grade (p = 0.0047), and more positive lymph node number (4.061 ± 3.905 vs. 1.794 ± 2.042; p = 0.0039) ( Table 2). We also found that patients with postoperative liver metastasis had a higher TP53 mutation rate (p = 0.0125) ( Table 3).

Identification of DEGs and Functional Analysis
Sixty-five patients, including 34 recurrence-free patients and 31 liver metastasis patients with mRNA sequencing data, were subsequently examined, and a total of 857 DGEs were identified (364 downregulated and 493 upregulated). All DGEs were illustrated by volcano plots and heatmap (Figures 2A,B). Then, we annotated functions of DEGs by KEGG pathway and GO analysis for biological processes (GO_BP). KEGG pathway analysis suggested that the DEGs were significantly enriched in "cytokine-cytokine receptor interaction, primary immunodeficiency, hematopoietic cell lineage, systemic lupus erythematosus, and pathway in cancer." The results of GO_BP analysis showed that DEGs were significantly enriched in "immune response, adaptive immune response, epidermis development, B cell receptor signaling pathway, and inflammatory response" (Figures 2C,D). This enriched KEGG pathway and GO_BP terms could further help us understand

Construction of PPI Network
The PPI network of DEGs was constructed for further analysis, including 460 nodes, 329 edges, and PPI enrichment Bolded values: statistically significant or trend. *P < 0.05.
p-value <1.0e −16 . Using the cytoHubba plugin based on Cytoscape to check the first stage nodes, the top 10 ranked proteins were identified as hub genes, which may play a significant role in the development of postoperative liver metastasis. Moreover, the shortest path and the expanded subnetwork associated with the 10 hub genes were identified ( Figure 3A). KEGG pathway and GO_BP analysis were carried out to annotate the 28 functions of genes. KEGG pathway analysis suggested that the function was significantly enriched in the "chemokine signaling pathway, cytokinecytokine receptor interaction, and neuroactive ligand-receptor interaction." GO_BP analysis showed significant enrichment in "G-protein coupled receptor signaling pathway, inflammatory response, immune response, and chemokine-mediated signaling pathway" (Figure 3B).

RPPA Analysis
The differences in expression levels of specific proteins between the two groups were further analyzed in 49 specimens: 25 in the recurrence-free group and 24 in the liver metastasis group. The results showed that CCNB1, CCNE1, GSK3, PAI-1, and p-Myosin IIa S1943 were significantly increased in the liver metastasis group, and TIGAR, HER3, PTEN, and G6PD were increased considerably in the recurrence-free group ( Figure 3C). Proteins with a p-value less than 0.1 were further included in the KEGG pathway and GO_BP analysis. KEGG pathway analysis suggested that the function was significantly enriched in "central carbon metabolism in cancer, prostate cancer, p53 signaling pathway, and pathways in cancer and insulin resistance." GO_BP analysis showed significant enrichment in "cellular response in  hypoxia, positive regulation of neuron death, and regulation of cell cycle" (Figure 3D).

Infiltration Level of Immune Cells
Functional analysis of DEGs suggested that there might be differences in tumor immune microenvironment between the two groups. With the estimation of XCELL, the infiltration levels of B cells, CD4 + naive T cells, myeloid dendritic cells, and natural killer T (NKT) cells were higher in the recurrence-free group, and Th2 CD4 + T cells, common myeloid progenitor cells, were higher in the liver metastasis group (Figure 4A).
With the estimation of TIMER, CD4 + T cells were significantly higher in the liver metastasis group. There is no significant difference in the B cells ( Figure 4B). We further measured the ESTIMATE score of these two groups, and the results suggested that the recurrence-free group had significantly higher immune and stroma scores ( Figure 4C). We did not find a correlation between the ESTIMATE scores with the DFS and OS for 65 patients (Supplementary Figure 1). However, there is a significant correlation between the ESTIMATE scores with the DFS and OS in the liver metastasis group (Figures 4D,E).

Construction of WGCNA and Identification of Key Modules
Based on the abovementioned studies, we found that the immune response within the tumor and cell cycle may be related to the occurrence of postoperative liver metastasis in pancreatic cancer. To further explore its mechanisms, we constructed the weighted gene co-expression network analysis (WGCNA).
WGCNA is a method to find out the significant cluster of genes associated with specific clinical traits. The WGCNA network comprised 30 samples with a specific time for the occurrence of postoperative liver metastasis. Clinical sample information included DFS, time of liver metastasis recurrence, and OS. All samples were clustered, and 28 samples were chosen for further analysis ( Figure 5A). The soft threshold was "6" (Figure 5B), and the cluster dendrograms of samples matching the clinical traits were obtained ( Figure 5C). Finally, 54 modules were obtained. We found that there were three modules positively correlated with DFS and OS. Also, there were two modules negatively correlated with DFS and OS ( Figure 5D). In this study, three modules (turquoise, brown, and purple) were extracted for further analysis because they had the most significant correlation with DFS and OS. Based on the eigengene adjacent heatmap, the three modules showed independent validation to each other ( Figure 5E).

Functional Analysis of Key Module Genes
From the module-trait relationship heatmap, the brown module was highly correlated with clinical traits (correlation coefficient = 0.54, p = 6.3E −92 ; Figure 6A). The brown module contained 1,202 genes and was correlated with DFS and OS positively. The key module genes were identified in a brown module under the threshold of module membership | MM| >0.85 and gene significance | GS| >0.3. To reveal the functions of these genes, we conduct the KEGG pathway and GO_BP analysis. KEGG pathway analysis suggested that the function was significantly enriched in the "chemokine signaling pathway, osteoclast differentiation, and hematopoietic cell lineage" (Figure 6B). GO_BP analysis exhibited significant , and ESTIMATE score of the two groups (C). Positive correlation between stroma score (D), immune score (E), and the DFS and OS in the liver metastasis group. * ,# P < 0.05, **P < 0.01. enrichment in "inflammatory response, immune response, innate immune response, regulation of immune response, cell surface receptor signaling pathway, and T cell activation" (Figure 6B). The turquoise module contained 3,205 genes and was correlated with DFS and OS positively (correlation coefficient = 0.44, p = 7E −152 ; Figure 6C). The KEGG pathway analysis showed that the function was significantly enriched in "renal cell carcinoma, proteoglycans in cancer, and pathways in cancer" (Figure 6D). GO_BP analysis showed significant enrichment in "transcription, DNA-templated, positive regulation of transcription, DNA-templated, and protein ubiquitination" (Figure 6D).
The purple module contained 407 genes and was correlated with DFS and OS negatively (correlation coefficient = 0.56, p = 5.7E −35 ; Figure 6E). The KEGG pathway analysis showed that the function was significantly enriched in the "cell cycle, oocyte meiosis, and Fanconi anemia pathway" (Figure 6F). GO_BP analysis was significantly enriched in "mitotic nuclear division, cell division, cell proliferation, and sister chromatid cohesion" (Figure 6F).
We further constructed the PPI network of the three hub modules (Figures 7A,C,E). Moreover, the outcomes of the KEGG pathway and GO_BP analysis of the three module genes are shown in Figures 7B,D,F

Hub Genes in the Purple Module
Finally, we found that the purple module was negatively correlated with DFS and OS, and genes CDK1, PLK1, AURKB, CENPN, KIF2C, etc., played the most important role in this module. This suggests that the increased expression of these genes might play an essential role in promoting the occurrence of postoperative liver metastasis of pancreatic cancer. Survival analysis of the top five genes in the purple module showed an inverse relationship between gene expression and prognosis (Figures 8A,B). Meanwhile, we found that the expression of these five genes was elevated in pancreatic cancer (Supplementary Figure 2A). Further analysis in the TCGA database showed that the expression of these five genes was closely related to postoperative recurrence and OS (Supplementary Figures 2B,C). (C) Clustering dendrograms of identified co-expressed genes in modules in the liver metastasis group. Each colored row represents a color-coded module that contains a group of highly connected genes. (D) Heatmaps of the correlation between eigengene and DFS and OS of the liver metastasis group. Each row corresponds to an ME, and each cell contains the corresponding correlation (red, positively correlated; blue, negatively correlated) and p-value. (E) Dendrogram of MEs obtained by WGCNA and heatmap plot of the adjacencies of modules. Red represents high adjacency (positive correlation), and blue represents low adjacency (negative correlation).

Tumor Cell Proliferation, CD4 + T Cells Infiltration, and Postoperative Liver Metastasis
The above analysis indicated that increased cell proliferation could promote the occurrence of postoperative liver metastasis of pancreatic cancer and the increase of CD4 + T cells infiltration was more obvious in the group without liver metastasis recurrence. In this study, immunohistochemical analysis showed that the tumor cells of patients with early postoperative liver metastasis of pancreatic cancer had higher proliferative activity ( Figure 9C) and less CD4 + T cells infiltration (Figure 9D). Moreover, the proliferative activity of tumor cells decreased with the prolongation of liver metastasis (Figure 9A). In addition, we found that most of the tumor cells were quiescent in recurrence-free patients, and that there was a significant infiltration of CD4 + T cells around the lesions of the pancreatic duct (follow-up time (FT) >3 years) (Figures 9A,B).

DISCUSSION
Pancreatic cancer is a malignant tumor of the digestive system. Due to local recurrence and distant metastasis, the postoperative prognosis of most patients is poor. Distinct locations of the first recurrence have unique independent predictive values for OS, which is helpful for stratification of prognosis and treatment decision after diagnosis of recurrence (Groot et al., 2018a). Several studies have reported the recurrence patterns after pancreatectomy for PDAC. Most postoperative patients with pancreatic cancer recurred at isolated distant sites, and liver recurrence tended to occur earlier (Groot et al., 2018b). The patterns of postoperative liver recurrence could be used as a basis for specific surveillance and treatment of patients with pancreatic cancer (Tanaka et al., 2019), which was verified by a metaanalysis including 89 studies with 17,313 patients undergoing PDAC resection. In this study, we presented an integrated view of hepatic recurrence after pancreatectomy at the molecular level.
We examined TCGA gene expression profile to identify potential biomarkers for postoperative liver metastasis of pancreatic cancer. Based on the KEGG analysis of DEGs, we found that DGEs were enriched highly in cytokine-cytokine receptor interaction. GO enrichment analysis showed that DEGs were significantly enriched in immune response and adaptive immune response. Furthermore, the KEGG analysis results revealed that key genes in the PPI network regarding DGEs were mostly enriched in the chemokine signaling pathway. In GO analysis, the enrichment of these genes mainly functioned to G-protein coupled receptor signaling pathway, inflammatory response, and immune response. Muller et al. (2001) validated that chemokines and their receptors play a vital role in determining the direction of tumor cell metastasis. G-protein coupled receptor kinases and receptors are also involved in tumorigenesis and cancer aggressiveness (O'Hayre et al., 2013;Wang et al., 2017;Nogues et al., 2018;Brown and Ganapathy, 2020). The immune response is implicated in tumor progression and metastasis in various ways (Berraondo et al., 2016;O'Sullivan et al., 2019;Riera-Domingo et al., 2020). It can protect the CD4 + naive T cells, myeloid dendritic cells, and NKT cells. However, Th2 CD4 + T cells, common myeloid progenitor cells, were higher in the liver metastasis group. These results imply that the occurrence of liver metastasis after pancreatectomy for PDAC is closely correlated with immune cell infiltration. Solid tumors, namely, breast cancer, renal cell carcinoma, melanoma, ovarian cancer, and gastrointestinal stromal tumors, have been shown to have spontaneous T cell infiltration (Zhang et al., 2003;Kreike et al., 2007;Mahmoud et al., 2011;Azimi et al., 2012;Rusakiewicz et al., 2013). In the TME of ovarian cancer, the ratio of CD8 + T cells to Foxp3 + regulatory T cells was significantly high, which was associated with a particularly favorable prognosis (Curiel et al., 2004).
Our further WGCNA showed that the brown module was correlated with DFS and OS positively. In addition, we found that the function of key module genes was highly enriched in the chemokine signaling pathway by KEGG pathway analysis and in the inflammatory response and immune response by GO analysis. However, the turquoise module and the purple module were correlated with DFS and OS negatively. The function analysis in the purple module possessed cell cycle in KEGG pathway analysis and mitotic nuclear division, cell division, and cell proliferation in GO analysis. Furthermore, in the purple module, the top five genes (CDK1, PLK1, AURKB, CENPN, and KIF2C) showed an inverse relationship between gene expression and prognosis, and their expression was closely related to postoperative recurrence and OS.
Among the top five genes, PLK1 was reported to be closely correlated with pancreatic cancer. Polo-like kinase 1 (Plk1), a serine/threonine kinase of polo-like kinases family, plays an essential role in cell mitosis, spindle assemblies, DNA damage, and so on (Lens et al., 2010;Strebhardt, 2010). Weichert et al. (2005) have reported that PLK1 was found to be overexpressed in pancreatic neoplasia as early as in pancreatic intraepithelial neoplasia III lesions and invasive pancreatic adenocarcinomas were PLK1, strongly positive in 47.7% of cases. Additionally, a study suggested that PLK1 showed >50% downregulation in gemcitabine sensitive cases and no change in the resistant cases (Jimeno et al., 2010). Correspondingly, inhibition of PLK1 also synergized with gemcitabine in gemcitabine-refractory in vitro (Jimeno et al., 2010) and in vivo models (Li et al., 2016). Moreover, the combination of PI3K/Akt pathway inhibition and PLK1 depletion can enhance chemosensitivity to gemcitabine in pancreatic cancer (Mao et al., 2018). According to the above results, PLK1 was proposed as a promising target for the treatment of gemcitabine-resistant pancreatic cancer. In a phase I trial, rigosertib, a PLK1 inhibitor, was combined with gemcitabine to treat patients with advanced PDAC (Ma et al., 2012). This phase I study determined the recommended phase II dose (RPTD), and under this, RPTD rigosertib is well tolerated with a toxicity profile similar to the gemcitabine alone. After demonstrating the safety of rigosertib, a phase II/III randomized study was then conducted (O'Neil et al., 2015). This phase II/III randomized study was designed to compare the efficacy and safety of rigosertib plus gemcitabine vs. gemcitabine alone in patients with metastatic pancreatic cancer. However, the combination of rigosertib plus gemcitabine failed to improve prognosis or response compared with gemcitabine alone in patients with metastatic pancreatic adenocarcinoma. The failure may be owed to the inherent heterogeneity of pancreatic cancer, and more comprehensive research regarding this topic is urgently needed. Hence, if overexpression of PLK1 for patients with pancreatic cancer after pancreatectomy was detected, non-gemcitabine chemotherapy would be preferred to avoid potential chemotherapeutic resistance. If gemcitabinebased chemotherapy has to be applied, PLK1 inhibitor could be considered to combine with it, and the resulting combination is safe and theoretically helpful.

CONCLUSION
Although we found that reduced immune cell infiltration was strongly associated with postoperative liver metastasis of pancreatic cancer, the internal mechanism still needs further investigation. The direct immune killing effect on tumor cells may be the main reason. Furthermore, previous studies have shown that tumor metastasis is closely related to the epithelial mesenchymal transition (EMT) process of tumor cells rather than tumor cell proliferation. However, our study showed that the proliferation of primary tumor cells was significantly increased in patients with postoperative liver metastasis. Therefore, we hypothesized that some patients with pancreatic cancer already had preclinical liver micro-metastases, which could not be detected by current detection methods. Such micro-metastases showed liver metastasis recurrence at different times after surgery due to the different proliferation abilities of tumor cells. According to our clinical tissue samples, patients with early postoperative liver metastasis had significantly more proliferative primary tumors than patients without recurrence.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/ Supplementary Material.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by the Ethics Committee of the First Affiliated Hospital with Nanjing Medical University. The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
GW designed and drafted the manuscript. SD and HG conducted the analyses. YG, LY, KZ, and XH conducted the literature search. ZL and YM designed, led, supervised the study, and critically revised the manuscript. All authors approved the final draft submitted, read and agreed to the published version of the manuscript.