CXCL10 is a Tumor Microenvironment and Immune Infiltration Related Prognostic Biomarker in Pancreatic Adenocarcinoma

Pancreatic adenocarcinoma (PAAD) is the 10th most common cancer worldwide and the outcomes for patients with the disease remain extremely poor. Precision biomarkers are urgently needed to increase the efficiency of early diagnosis and to improve the prognosis of patients. The tumor microenvironment (TME) and tumor immune infiltration are thought to impact the occurrence, progression, and prognosis of PAAD. Novel biomarkers excavated originating from the TME and immune infiltration may be effective in predicting the prognosis of PAAD patients. In the current study, the ESTIMATE and CIBERSORT algorithms were applied to estimate the division of immune and stromal components and the proportion of tumor-infiltrating immune cells in 182 PAAD cases downloaded from The Cancer Genome Atlas database. Intersection analyses of the Protein-Protein Interaction networks and Cox regression analysis identified the chemokine (CXC-motif) ligand 10 (CXCL10) as a predictive biomarker. We verified that CXCL10 in the TME negatively correlates with prognosis in PAAD and positively correlates with tumor cell differentiation. GSE62452 from the GEO database and cumulative survival analysis were performed to validate CXCL10 expression as an independent prognostic indicator. We also found that memory B cells, regulatory T cells, and macrophages M0 and M1 were correlated with the expression of CXCL10 indicating that expression of CXCL10 influenced the immune activity of the TME. Our data suggest that CXCL10 is beneficial as a prognostic indicator in PAAD patients and highlights the potential for immune targeted therapy in the treatment of PAAD.


INTRODUCTION
Pancreatic adenocarcinoma (PAAD) is one of the most deadly malignant tumors and is ranked the seventh leading cause of cancer death (n 432,000) (Bray et al., 2018). PAAD is commonly diagnosed in the advanced stage and currently no effective therapies are available (Yan et al., 2017). The prognosis in PAAD is highly unsatisfactory with 1-year survival less than 25% and 5-year survival of no more than 5% (Ottenhof et al., 2009). There is an urgent clinical need to identify novel biomarkers for the diagnosis and prognosis of PAAD.
Recently, research has focused on the tumor microenvironment (TME) and its role in cancer. The TME consists of the extracellular matrix, soluble molecules, and tumor stromal cells (Thomas and Radhakrishnan, 2019). Multiple studies have demonstrated that the TME can improve the invasiveness of tumor cells by modulating proliferation, chemotherapy resistance, immune escape, and metastasis (Chen et al., 2015;J. H. Park et al., 2017). Immune and stromal cells are the two main types of nontumor components in TME which have been considered to be valuable for tumor diagnosis and prognosis evaluation (Wörmann et al., 2014;Xu et al., 2014). T cells, B cells, and macrophages are chemotactic constituents within the TME. Several studies have demonstrated the role of stromal cells in extracellular matrix remodeling and tumor angiogenesis. Also, PAAD is characterized by an intense stromal desmoplastic reaction around the cancer cells (Bussard et al., 2016;Ligorio et al., 2019). Previous studies have demonstrated that infiltrating immune cells can be isolated from tumors suggesting that tumor immune infiltration is a crucial biological process that occurs as a result of immune cell migration from the blood into the TME (Ge et al., 2019;R. Liu et al., 2020). Multiple reports have suggested that the proportion and functional alterations of different tumor infiltrated immune cells contribute to the initiation and progression of PAAD (Meng et al., 2020;Tian et al., 2020). However, the regulatory mechanisms of immune cells, stromal components of TME, and tumor immune infiltration in PAAD remain largely unknown.
Based on the essential roles of the immune and stromal components within the TME, and tumor immune infiltration in carcinogenesis and progression, molecules extracted from the TME and tumor immune infiltration cells could have major potential as biomarkers in PAAD. Considering the pivotal roles of immune, stromal cells and tumor immune infiltration in PAAD, we first calculated immune and stromal scores and quantified the proportion of tumor immune cells using the Estimation of Stromal and Immune cells in Malignant Tumors using Expression data (ESTIMATE) algorithm. We used the CIBERSORT computational method to identify novel and potential therapeutic targets in PAAD. From differential gene expression analysis generated by comparing the immune and stromal components, we identified chemokine (CXCmotif) ligand 10 (CXCL10) as a potential predictive biomarker. CXCL10 belongs to the subfamily of the CXC chemokines identified as the main chemokine family in humans . Originally, CXCL10 is shown to be proinflammatory and proliferative and is associated with advanced human cancer (Lunardi et al., 2015;Dai et al., 2020). In this study, using bioinformatics analyses, we demonstrated that CXCL10 was a potentially valuable biomarker for alterations within the TME in PAAD.

Gene Expression Data
A total of 182 cases of PAAD transcriptome FPKM data and clinicopathological characteristics were downloaded from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer. gov). These cases included 4 normal cases and 178 tumor cases. Due to the small scale of the normal samples in TCGA database, we used the UCSC Xena database (https://xenabrowser.net/datapages/) to obtain gene expression RNAseq FPKM data from the GTEx cohort. Gene expression data from these two databases were combined with the 'normalize Between Arrays' function in R language (version 4.0. 2). To verify the precision of gene expression, mRNA expression data in PAAD were searched and downloaded from the GEO DataSets in NCBI (https://www.ncbi.nlm.nih.gov/gds/) using the keywords "pancreatic cancer," "PAAD," "microarray," and "adjacent." The GSE62452 cohort contained 69 tumor samples and 61 adjacent normal tissue samples that were selected for further analysis.

ESTIMATE Algorithms
The ESTIMATE algorithms included three kinds of scores, specifically an immune score, a stromal score, and an ESTIMATE score. The immune and stromal score were calculated based on the relative proportion of the immune and stromal elements. ESTIMATE scores were the sum of the immune and stromal scores. Three scores were calculated from the features of "limma" and "estimate" packages in R language loaded from gene expression files.

Differential Expression Analysis and Enrichment of Genes
Based on the results of ESTIMATE algorithms, all PAAD cases were divided into high and low groups according to the immune and stromal scores, respectively. The package named "pheatmap" was run to screen the heat-maps of the immune and stromal scores. The "limma" package from R language was performed to analyze the differentially expressed genes in the immune and stromal score groups with a threshold value of p < 0.05 and |log fold change|>1 by a Wilcoxon rank-sum test.

GO and KEGG Enrichment Analyses
Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of the 772 differentially expressed genes were selected by the "clusterProfiler," "ggplot2," "org.Hs.eg.db," and "enrichplot" packages. Both the p values and q values less than 0.05 of the samples were considered to be significantly enriched.

PPI Network Construction and Cox Regression Analysis
The Search Tool for the Retrieval of Interacting Gene (STRING) database was used to establish the Protein-Protein Interaction (PPI) network. Nodes where the interactive relationships were greater than 0.99 were selected to build the network. R language with the aid of package "survival" was set up for univariate Cox regression analysis which listed the top 16 genes ordered by the p values (p < 0.05).
"survival." Based on the Kruskal-Wallis rank-sum test, the "ggpubr" package in R language was run to explore the relationships between the stromal and immune scores and the clinicopathological parameters. Multivariate independent prognostic analysis was employed with the "survival" package in R language. A difference of p < 0.05 indicated statistical significance.

Gene Set Enrichment Analysis
C7.all.v.7.1symbols.gmt from MSigDB was used as the target set to carry out Gene Set Enrichment Analysis (GSEA) analysis of all tumor cases with the software GSEA (version 4.1.0) downloaded from the Broad Institute. Only gene sets with NOM p value < 0.05 were considered as statistically significant.

Evaluation of Tumor Infiltrated Immune Cells
The fractions of the tumor infiltrated immune cells of all tumor cases were calculated using the CIBERSORT algorithm. When the p value of CIBERSORT was less than 0.05, the data were filtered and selected for further analysis with the "limma" package in R language. Difference and correlation analyses were applied using the "limma," "vioplot," "gglot2," "ggpubr," and "ggExtra" packages in R language. p values less than 0.05 were regarded as statistically significant using the Wilcoxon rank-sum and Pearson coefficient tests. Also, cumulative survival analysis was based on the Tumor Immune Estimation Resource (Timer) 2.0 database (http://timer.cistrome.org/).

Flowchart of the Data Analysis Procedure
A flowchart of the data analysis procedure used in this study is summarized in Figure 1. The fractions of tumor infiltrated immune cells and proportions of the immune and stromal component in the 182 cases of PAAD downloaded from TCGA were evaluated using the ESTIMATE and CIBERSORT algorithms. The tumor samples were evaluated based on the immune, stromal, and ESTIMATE scores, and each cohort included its own set of differentially expressed genes. A total of 772 overlapping differentially expressed genes were found in the groups from the up & downregulated immune and stromal scores. Using PPI network construction and Cox regression analysis to calculate the differentially expressed genes, an intersection analysis was applied to the datasets. Based on these analyses, we found that CXCL10 was closely associated with the prognosis of PAAD patients. Subsequently, analysis of survival, clinical correlation, validation, and correlation were performed to evaluate the prognostic value of CXCL10.

Remodeling of Tumor Immune Infiltration and TME Indices Are Closely Related to Clinicopathological Characteristics of PAAD
A previous study has reported that the ESTIMATE algorithms, either the stromal or the immune score, are indicators of patient survival, relapse, metastasis, and chemotherapeutic drug resistance (W. Liu et al., 2018). Higher immune or stromal scores represented larger amounts of immune or stromal components in the TME. The feasible associations of the scores and clinicopathological characteristics were then explored. The clinicopathological characteristics of the PAAD patients are summarized in Supplementary Table S1, which were analyzed after detecting the data. The data showed the immune, stromal, and ESTIMATE scores were all significantly associated with gender as well as G1 and G2 histologically graded disease ( Figure 2). The p values of gender however, no significance was found between the scores and the indicated clinicopathological features (Supplementary Figures  S1B-S3A). These data suggested that the immune and stromal scores were closely related to the progression of PAAD. Heat-maps of top 20 differentially expressed genes between high score group and low score group based on immune score and stromal score with the threshold value of p < 0.05 and |log fold change|>1 by Wilcoxon ranksum test, respectively.
Frontiers in Molecular Biosciences | www.frontiersin.org February 2021 | Volume 8 | Article 611508 5 Identification of Differentially Expressed Genes in Tumor-Infiltrating Immune Cells and the TME According to the median score of the stromal and immune scores, the PAAD patients were divided into two groups. Differential expression analysis identified differentially expressed genes in the immune (both low and high scores) and stromal score groups (both low and high scores). Heat-maps displayed the differential expression of these genes. The upregulated (logFC > 0) genes represent genes that were significantly elevated compared to the downregulated genes (logFC < 0). Based on the immune scores, a Frontiers in Molecular Biosciences | www.frontiersin.org February 2021 | Volume 8 | Article 611508 6 total of 901 differentially expressed genes were acquired, among which 822 genes were upregulated and 79 genes were downregulated ( Figure 3A). Similarly, 1,372 differentially expressed genes were obtained in the two stromal score groups, among which 1,103 genes were upregulated and 169 genes were downregulated ( Figure 3B). Also, intersections between the high stromal and immune score group or the low stromal and immune score group were visualized in a Venn diagram. Data showed 715 common differentially expressed genes that were upregulated in both the immune and stromal score groups and 57 common differentially expressed genes that were downregulated in both the immune and stromal score groups (Supplementary Figure S3B). A total of 772 differentially expressed genes may be crucial components of the TME.

Differentially Expressed Genes Have Principal Roles in Mediating Immune-Related Functions
GO and KEGG pathway enrichment analyses were conducted to identify the latent molecular mechanisms of 772 differentially Frontiers in Molecular Biosciences | www.frontiersin.org February 2021 | Volume 8 | Article 611508 7 expressed genes. The results of the GO functional enrichment analysis indicated that the differentially expressed genes were enriched in immune-related functions including T-cell activation, regulation of lymphocyte activation, and leukocyte migration ( Figure 4A). The results of the KEGG functional enrichment analysis showed differentially expressed genes mainly enriched in cytokine-cytokine receptor interactions, hematopoietic cell lineages, and chemokine signaling pathways ( Figure 4B). These data highlighted both functional enrichment analyses appeared to focus on immune-related functions suggesting that the participation of immune components has principal roles in the TME in PAAD.

CXCL10 Is the Only Intersectional Gene in the PPI Network Construction and Univariate Cox Regression Analysis
To visualize potential interactions between the differentially expressed genes, PPI network construction based on the STRING database was performed that covered 222 nodes and 476 edges (Supplementary Figure S4). The top 30 differentially expressed genes were sorted and summarized by the number of nodes ( Figure 5A). Univariate Cox regression analysis was performed to further investigate the significant elements of these differentially expressed genes in the survival of PAAD and the top 16 genes sorted by p value were identified ( Figure 5B). According to the intersection analysis of the PPI network and the univariate Cox regression analysis, CXCL10 was shown to be the only overlapping element across these two analyses ( Figure 5C).

CXCL10 Is Strongly Correlated With Prognosis and the Clinicopathological Characteristics of PAAD Patients
A previous study demonstrated that CXCL10 is a proinflammatory chemokine and chemoattractant for T cells. Overexpression of CXCL10 has been shown to promote tumor growth, migration, and invasion via targeting of the cognate receptor chemokine (CXC-motif) receptor (CXCR3) (Lunardi et al., 2015). In this study, based on the median expression of CXCL10, all PAAD cases were divided into high or low CXCL10 expression groups. The Wilcoxon rank-sum test was carried out and showed the correlation of CXCL10 expression with clinical characteristics. No significant difference in the expression of CXCL10 was found between normal and tumor cases from the TCGA database ( Figure 6A). Interestingly, when the normal cases were added from the GTEx database, the expression of CXCL10 in the tumor cases was higher compared to the normal cases ( Figure 6B, p < 0.001). Survival analysis indicated that CXCL10 expression was negatively related to outcomes in PAAD patients ( Figure 6C). Also, the expression of CXCL10 correlated with gender, G1, G2, and G3 of histologic grades with p values of 0.015, 0.029, and 0.034, respectively ( Figures 7A,B). All of the above results show that overexpression of CXCL10 in the TME had a negative relationship with prognosis in PAAD and was positively related to tumor cell differentiation. No significance was found in the indicated clinical parameters and CXCL10 expression ( Figures 7C-G).

Validation of CXCL10 Expression and Survival in PAAD
GSE62452 from the GEO databases was analyzed to verify the expression of CXCL10 in PAAD. The expression of CXCL10 was higher in the tumor samples compared to the normal cases ( Figure 6D, p < 0.001). Higher CXCL10 expression was associated with poorer survival compared to low CXCL10 expression ( Figure 6E, p 0.016). Moreover, the multivariate independent prognostic analysis revealed that the expression of CXCL10 was an independent prognostic indicator for PAAD (p 0.035) ( Figure 6F). Based on these data, we suggest that CXCL10 may be a potential biomarker that reflects the status of TME in PAAD.

CXCL10 Impacts the Immune Activity of the TME
To further confirm the association of CXCL10 expression and the immune microenvironment, C7 collection from MSigDB in GSEA was used to analyze data between the high and low CXCL10 expression groups. Most immune functional gene sets were enriched in the high CXCL10 expression group while few gene sets were enriched in the low CXCL10 expression group (Figure 8; Supplementary Table S2). The fraction of tumor infiltrated immune cells was analyzed using CIBERSORT and 22 kinds of immune cells in PAAD were constructed (Figure 9). Through the intersection of difference and correlation analyses, 3 kinds of tumor-infiltrating immune cells were negatively associated with the expression of CXCL10 which included memory B cells, regulatory T cells (Tregs), and macrophages M0. Macrophages M1 were shown to be positively correlated with CXCL10 expression (Figure 10; Supplementary Table S3).
To validate the correlation of CXCL10 expression in immune cells, cumulative survival analysis was performed. Indeed, the infiltration levels of Tregs, B cell, and macrophage M0-CIBERSORT subtype were negatively correlated with CXCL10 expression, while the infiltration level of macrophage M1-CIBERSORT subtype showed positively associated with CXCL10 expression (Supplementary Figure S5). Collectively, all the data prove that aberrant expression of CXCL10 impacts the immune activity of the TME which is a promising biomarker for the prognosis of PAAD patients.

DISCUSSION
Recently, remodeling of the TME has been shown to play a crucial role in the initiation and progression of multiple cancer types (Cheng et al., 2019;Jiang et al., 2020;Park and Surh, 2017). Also, several studies have been performed to dissect the features of the TME in PAAD (Karamitopoulou, 2019;Hessmann et al., 2020;Meng et al., 2020). As surgical resection is the only option for a limited number of eligible patients (Petrusel et al., 2020), Frontiers in Molecular Biosciences | www.frontiersin.org February 2021 | Volume 8 | Article 611508 specific signaling factors and molecules from TME involved in the carcinogenesis and progression of PAAD remain to be fully determined. Exploring the potential mechanisms of immune and stromal cells within the TME (Tsai et al., 2018) may be a promising strategy that could potentially inform the development of novel therapeutic options in PAAD. To date, the ESTIMATE algorithms could be applied to calculate the components of immune and stromal cells (H. Wang et al., 2019) to obtain corresponding immune and stromal scores. Our data from transcriptomic analysis based on PAAD samples revealed the proportion of immune and stromal components in TME contributed to the prognosis of PAAD. These consequences highlighted the significance of exploring the association of tumor cells and immune cells which provided constructive insight for developing much more effective treatment regimen. Though immunotherapy has become support of cancer treatment and brings new hope for PAAD (Fritz and Lenardo, 2019), the developments of individual patient confront with tough challenges due to side effects (Sahin et al., 2017). It is necessary to investigate some novel candidates for the immunotherapy of PAAD. Here, we set out from transcriptomic analysis of PAAD, which implied the increased expression of CXCL10 was significantly associated with poor prognosis and advanced clinicopathological characteristics.
CXCL10, known as interferon-(IFN-) γ-induced protein 10 (IP-10), belongs to the CXC chemokine subfamily and contains a single and variable amino acid between two of the four highly conserved cysteine residues (Müller et al., 2010;Antonelli et al.,  2014; Qian et al., 2019). Several studies have demonstrated the significant role of CXCL10 in tumor, including PAAD. The effects of CXCL10 are mediated by binding to CXCR3 receptor. Significantly higher levels of CXCL10 and CXCR3 are observed in PAAD specimens compared to those from chronic pancreatitis (Singh et al., 2007;Moin et al., 2018). Recent studies have confirmed a stromal origin of CXCL10 in human tumors which is overexpressed in human pancreatic cancer and associated with poor survival of PAAD patients (Delitto et al., 2015). Recently, CXCL10 and CCL21 are demonstrated to promote migration of pancreatic cancer cells toward sensory neurons and increase the frequency of cancerassociated pain (Hirth et al., 2020). Moreover, High expression of CXCL10 is also contributed to chemotherapeutic gemcitabine (Delitto et al., 2015). Consistently, in our study, we also validated the bioinformatic characteristics of CXCL10, whereby CXCL10 expression was an independent predictor of poor survival in PAAD. Further work needs to investigate molecular mechanism of CXCL10 in regulating malignant phenotypes of PAAD. Tumor infiltrated immune cells can regulate the proportion and distribution of immune cells in tumors and tumor inflammation is abnormally active especially in the early stage of carcinogenesis (Neviani et al., 2019). We found CXCL10 was correlated with early stage of PAAD which indicated potential role of CXCL10 in PAAD carcinogenesis. Recently, macrophage M0 is demonstrated to harbor antitumorigenic activities to suppress PAAD cell growth by TNF-α secretion, but not M1 or M2 (Tekin et al., 2020). This could be an explanation why CXCL10 is negatively related to macrophage M0 in our study. Another issue is that we found CXCL10 was positively correlated with macrophage M1. Macrophages M1, that act as antitumoral immune components, are clarified to prevent immune escape in PAAD (Young et al., 2020). Tumor expressed galectin 9 activates dectin 1 on macrophage which subsequently promotes PAAD and enhances immune tolerance (Daley et al., 2017). Combined with the proinflammatory properties of macrophage M1, we speculated that aberrant upregulation of CXCL10 and increased macrophage M1 may be an early event in the process of carcinogenesis.
Tumor-infiltrating Treg cells can suppress tumor-associated DC immunogenicity to promote immune tolerance in pancreatic cancer (Jang et al., 2017). Moreover, FOXP3 + Tregs recruited in TME contribute to immunosuppressive microenvironment which subsequently promote immune escape and metastasis of PAAD (X. Wang et al., 2017). Also, circulating Tregs are identified as prognostic factors and can predict chemotherapeutic response in unresectable PDAC patients (C. Liu et al., 2017). These works collectively emphasize infiltration of Tregs in TME of PAAD. However, in current work, we found CXCL10 was negatively correlated with Tregs. Previous study identified systemic dysfunction and plasticity of the immune macroenvironment and found CD103 + Tregs were recruited at day 7 but decreased with tumor progression, indicating dynamic distribution of Tregs in TME during the process of tumor initiation and progression (Allen et al., 2020). Additionally, pancreatic stellate cells (PSCs) secrete CXCL10 induced by pancreatic cells to recruit CXCR3 + and FOXP3 + Tregs, which subsequently suppress adaptive immune responses and promote immunosuppression (Lunardi et al., 2015). Thus, we speculate that, during the progression of PAAD, the correlation between the expression of CXCL10 in tumor cells and the abundance of Tregs cells will be weakened. Further work needs to focus on exploring the correlation between CXCL10 and different subtypes of Tregs and then systematically understand the dynamic changes of CXCL10 in the progress of PAAD.
Collectively, based on the ESTIMATE algorithms, functional enrichment analysis, PPI network construction, and Cox regression analysis, we suggest CXCL10 is a promising diagnostic and prognostic indicator in PAAD patients which has potential to provide novel immunotherapy insights for PAAD. Further studies are required to explore CXCL10 as a therapeutic target in the treatment of PAAD.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
LL conceived and designed the study. HH and WZ collected and analyzed and interpreted the data. HH, RC, BX, and SZ participated in interpreting and revising the manuscript. LL and HH wrote the manuscript. All authors contributed to manuscript revision, read, and approved the submitted version.