Integrative Analysis of Gene Expression Based on Multi-Omics Databases Identied CXC Chemokines as Immune Prognostic Biomarkers of Rectal Adenocarcinoma

Background: Rectal adenocarcinoma (READ) is one of the most frequent malignancies with a high recurrence rate. CXC chemokines, as indispensable components of the immune system, are considered broadly involving in tumorigenesis by orchestrating the immune cell chemotaxis and thus affect the prognosis of READ patients. However, the values of CXC chemokines as prognostic biomarkers and potential regulatory mechanisms for READ remain unclear. Results: The expression levels of CXCL3, CXCL12, and CXCL13 were aberrant in both TCGA and ONCOMINE databases. Lower expression of CXCL3 and CXCL13 predicted poor survival of READ patients. Additionally, both CXCL3 and CXCL13 were associated with several clinicopathological features. CXCL3 and CXCL13 expressions were signicantly correlated with the tumor inltration levels of immune cells in READ tissue. CeRNA networks of mRNA-miRNA-lncRNA were constructed to reveal the potential mechanisms that regulated the expressions of CXCL3 and CXCL13. Furthermore, GSEA revealed the association between immune-related pathways and CXCL3 as well as CXCL13. Conclusions: CXCL3 and CXCL13 could be valuable prognostic biomarkers in READ. CXCL3/miR-425-5p/chr22-38_28785274–29006793.1 was identied as the most potential ceRNA network in READ. Our results might provide novel insights in READ immunotherapy. 0.05 and FDR q < 0.25 were considered signicant. The results were presented using R 3.6.2.


Background
Rectal adenocarcinoma (READ) is one of the most frequent malignancies with a high recurrence rate, with estimated 400,000 new cases every year around the world (1). Although great progress in standard treatments for locally advanced rectal cancer, there are still a considerable amount of patients with local or distant recurrence (2). Recently, immunotherapy, especially immune checkpoint blockade therapy has emerged as a novel strategy to treat different types of tumors, including colon cancer and rectal cancer (3). Immunotherapy has great potential for treating patients with READ, but the therapeutic targets and the corresponding mechanisms remain unclear.
Chemokines are subordinate to the family of low-molecular-weight heparin-binding chemotactic cytokines, which induce immunocyte migration (4). They are currently classi ed into four subfamilies based on the N-terminal cysteine (C) motifs: C, CC, CXC, and CX 3 C, in which X represents any amino acid (5,6). In-depth studies on cancer immunotherapy suggested that CXC chemokines can modulate anti-tumor immunological responses by regulating tumor immune in ltrates of speci c immunocytes such as CD8 + T cells and CD4 + T cells, which regulate the outcome of immunotherapy (7,8). A total of sixteen CXC chemokines (CXCL1 -CXCL17,excluding CXCL15) have been discovered, and they might be key circuits regulating tumor immunosuppression via mediating effector T cell tra cking (7,9). The expression of CXC chemokines may also predict the prognosis of patients with different kinds of tumors. For instance, higher CXCL9, CXCL10, and CXCL13 predict better overall survival of melanoma patients (10), while breast cancer patients with decreased expression of CXCL8 presented with better outcomes (11). Although several studies have suggested the prognostic and therapeutic values of CXC chemokines in rectal cancer (12)(13)(14), comprehensive analysis of CXC chemokines with su cient samples are yet to be performed.
In this study, we presented an integrative bioinformatics analysis of multi-omics data including expressions, gene alterations, immune in ltrates, and functional analyses of CXC chemokines in READ ( Fig. 1). We evaluated the potential of sixteen CXC chemokines as prognostic biomarkers for READ. We also explored the association between CXC chemokines and immune in ltrates as well as the underlying mechanisms. These data might illustrate the prognostic and immunotherapeutic values of CXC chemokines in READ treatment.

Identi cation of differentially expressed CXC chemokines
Multiple databases were utilized to demonstrate the expression patterns of CXC chemokines. We rst explored the expression levels of sixteen CXC chemokines using 131 rectal adenocarcinoma (READ) samples and normal rectum tissues from the TCGA-READ dataset with complete clinicopathological information. CXCL3 and CXCL16 were upregulated, while CXCL12 and CXCL13 were downregulated in READ tissues compared with normal tissues (Figs. 2A-2D). We then assessed the expression of CXC chemokines in READ and normal tissues using the ONCOMINE database (Table 1 and Additional le 1:   Table S1). Based on the data from ONCOMINE, nine CXC chemokines (CXCL1, CXCL2, CXCL3, CXCL5, CXCL6, CXCL8, CXCL10, CXCL11, and CXCL17) were upregulated, while three CXC chemokines (CXCL9, CXCL12, and CXCL13) were downregulated in READ samples. The expression patterns of CXCL3, CXCL12, and CXCL13 were consistent in both TCGA and ONCOMINE (Fig. 2E). We also compared the relative expression levels of sixteen CXC chemokines (CXCL15 was not included as the data were not available) using GEPIA (Additional le 2: Figure S1).

2.2.
Genetic alteration analysis, co-expression analysis, and protein-protein interaction (PPI) network construction

The prognostic value of CXC chemokines in READ patients
To assess the prognostic value of CXCL3, CXCL12, and CXCL13 in READ patients, we explored the correlation between the expression of these CXC chemokines and the overall survival of READ patients using the TCGA dataset. Kaplan-Meier curves showed that at median cutoff, lower expression levels of CXCL13 (P = 0.021) were associated with poor survival (Fig. 6). On the other hand, lower CXCL3 (P = 0.006) and CXCL13 (P = 0.001) levels correlated with worse outcome at best cutoff. Hence, we analyzed the correlation between clinicopathological features and expression levels of CXCL3 and CXCL13, respectively ( Table 2). Lower CXCL3 expression was signi cantly correlated with a more advanced pathological stage (P = 0.007), T stage (P = 0.036), and N stage (P = 0.035). Lower CXCL13 expression was associated with an advanced pathological stage (P < 0.001).

CeRNA networks construction & gene set enrichment analysis of CXCL3 and CXCL13
To further explore the potential mechanism underlying the dysregulation of CXCL3 and CXCL13, we constructed mRNA-miRNA-lncRNA ceRNA networks using DIANA tools based on the theory of ceRNA. For the CXCL3 network, a total of 12 miRNAs and 14 lncRNAs were predicted (Fig. 8A). As for the CXCL13 network, a total of eight miRNAs and 22 lncRNAs were involved (Fig. 8B). As for these miRNAs, miR-425-5p was the unique oncogene reported in colorectal cancer (CRC) (30,31). LncRNA chr22-38_28785274-29006793.1 was predicted as the ceRNA of miR-425-5p. We ascertained that the CXCL3/miR-425-5p/chr22-38_28785274-29006793.1 network was the hub potential regulator in the pathogenesis of READ.
GSEA on CXCL3 and CXCL13 using the TCGA-READ dataset was also presented for downstream functional analysis (Additional le 3: Figure S2). CXCL3 was positively related to immune-related GO terms such as antimicrobial humoral immune response mediated by antimicrobial peptide and antimicrobial humoral response. CXCL13 was also positively associated with several immune-related GO terms, such as B cell activation, B cell differentiation, and T cell activation. As for KEGG pathway analysis, most of the pathways predicted in the higher CXCL3 group are not related to the immune system. On the other hand, terms such as antigen processing and presentation, B cell receptor signaling pathway, and T cell receptor signaling pathway are related to the higher CXCL13 group.

Discussion
Since the imbalance between normal cells and the immune system is considered to induce tumorigenesis, several studies have demonstrated signi cant roles for CXC chemokines in tumor progression (32)(33)(34), tumor microenvironment (35), and prognostic outcomes(36). Previous studies suggested that CXC chemokines, especially CXCL3 and CXCL13, presented as the prognosis markers in rectal adenocarcinoma (READ). For instance, lower expression of CXCL3 exhibited poor overall survival in READ (13). In the rectal cancer liver metastases, M2 macrophage polarization promoted tumor colonization by secreting CXCL13 and activated the CXCL13/CXCR5 axis in rectal cancer cells (37). These ndings indicated that both CXCL3 and CXCL13 possibly had dual biological functions (tumor promoter/tumor suppressor) in READ, which consisted of the previous researches about chemokines (38). In detail, chemokines might promote metastasis by acting directly on the invasion and migration of tumor cells, while chemokines also recruited immune cells towards tumors (39,40). In different diseases, each chemokine might have its unique function (38). However, the role of the immune microenvironment in READ remained unclear, and studies with large samples were needed. Therefore, we aimed to investigate whether CXC chemokines played crucial roles in the regulation of the immune microenvironment in READ.
We rst characterized the CXC chemokine expression patterns in READ. The aberrant expression of three CXC chemokines (CXCL3, CXCL12, and CXCL13) was consistent in READ tissues compared with normal ones in both TCGA and ONCOMINE databases. Signi cantly increased CXCL3 level was observed in READ tissues compared with normal ones, while decreased CXCL12 and CXCL13 levels were observed in READ tissues compared with normal ones. Further survival analysis revealed that both lower expression of CXCL3 and CXCL13 predicted a poorer prognosis of READ. As CXC chemokines mediate immunocyte migration(4), we also found that the expression of CXCL3 and CXCL13 correlated with the in ltration of several immune cells in READ. The expression of CXCL3 was negatively and signi cantly correlated with macrophage in ltration, while positively correlated with neutrophil in ltration. On the other hand, CXCL13 level was negatively and signi cantly correlated with tumor purity, while positively correlated with the in ltration of B cell, CD8 + T cell, neutrophil, and dendritic cell. Furthermore, there was a signi cant positive correlation between CNV of CXCL13 and macrophage in ltration.
CXCL3 is reported as a chemoattractant of neutrophils (41,42), which is consistent with our TIMER results. Multiple studies indicate that CXCL3 is an oncogene that promotes the progression of various malignancies, such as prostate cancer, breast cancer, and colorectal cancer (20,43,44). Our results show that CXCL3 is upregulated in READ, indicating that CXCL3 could also be oncogenic in READ. Surprisingly, upregulated CXCL3 is associated with a better prognosis. Similarly, CXCL3 is also upregulated in colon cancer, in which upregulated CXCL3 indicates better clinical outcomes (45). In esophageal squamous cell carcinoma(46), upregulated CXCL3 secreted from cancer cells unexpectedly increased the recruitment and the antitumor effect of neutrophils. Collectively, these ndings indicate that in READ, CXCL3 may present dual character, which promotes tumor progression and upregulates neutrophil recruitment simultaneously.
On the other hand, CXCL13 was originally identi ed in the stromal cells from B cell follicles which could regulate homing of B cells and follicular T cells (Tfh) (47,48). It is constitutively secreted by stromal cells in secondary lymphoid tissues, including the spleen and lymph nodes. This is consistent with our results as CXCL13 expression positively correlated with the in ltration level of B cells in READ. In several solid tumors, CXCL13 can also promote tumor cell proliferation and metastasis through CXCL13/CXCR5 axis, such as colorectal cancer (49), breast cancer (50,51), and prostate cancer (52,53). In colorectal cancer(54), lung cancer (55), and gastric cancer(56), CXCL13 is upregulated and serves as a prognostic biomarker as well as a therapeutic target, in which lower CXCL13 expression associates with the better prognostic outcome. Reversely, our results show that CXCL13 is downregulated in READ, and lower CXCL13 expression is linked to a worse prognosis. In terms of mechanism, downregulated CXCL13 may reduce B cells and Tfh cells in ltration, and thus lead to READ progression and poor prognosis, which is similar to breast cancer [57] and colorectal cancer [58]. These ndings indicate the dual character of CXCL13, similar to CXCL3, in different tumors [59], and CXCL13 mainly performs anti-tumor effect on READ.
Long non-coding RNAs (lncRNAs) function effect as miRNAs sponges to interfere the functions of miRNAs targeting to mRNAs in various tumors[60, 61]. In our study, 12 miRNAs and 14 lncRNAs were identi ed as ceRNA networks of CXCL3, while eight miRNAs and 22 lncRNAs were identi ed targeting CXCL13. Among all these miRNAs, miR-425-5p, which downregulated CXCL3[62], was reported as an oncogene [20,21] and inducing chemo-resistance[63] in colorectal cancer (CRC). It was consistent with our identi cation that lower expression level of miR-425-5p and higher expression level of CXCL3 were presented in CRC. Additionally, it was reported that overexpression of miR-425-5p could induce the M2 polarization of macrophages [31], which enhanced tumor metastasis by secreting vascular endothelial growth factors (VEGF). It suggests that miR-425-5p may regulate M2 polarization of macrophages via CXCL3. According to prediction of miRNA-lncRNA from database DIANA, we found that miR-425-5p bound to lncRNA chr22-38_28785274-29006793.1. Regarding to lncRNA chr22-38_28785274-29006793.1, it was considered as the member of mRNA/miRNA/lncRNA ceRNA networks and involved in the in ltration of CD4 + and CD8 + T cell in colon cancer by bioinformatic analysis without biological validation [64]. Thus, we assume that CXCL3/miR-425-5p/chr22-38_28785274-29006793.1 network may upregulate CXCL3 in READ and promote immune in ltration to exert the anti-tumor immunological effect.
Several limitations should be addressed. First, the bioinformatics analysis only focused on the transcriptional level of CXC chemokines, which could only re ect speci c changes but not global changes of READ immune status. Second, the data were acquired from public resources, so the bias in the study cannot be completely eliminated. Last, ceRNA networks were predicted based on bioinformatics algorithms instead of experiments. To overcome these issues, further experiments in vitro and in vivo are warranted.

Conclusions
CXCL3 and CXCL13 could be valuable prognostic biomarkers in READ. CXCL3/miR-425-5p/chr22-38_28785274-29006793.1 was identi ed as the most potential ceRNA network in READ. Our results might provide novel insights in READ immunotherapy.

The Cancer Genome Atlas (TCGA)
The Cancer Genome Atlas (TCGA) (https://tcga-data.nci.nih.gov/tcga/) is a publicly funded project that aims to catalog and discover major cancer-causing genomic alterations to create a comprehensive "atlas" of cancer genomic pro les [65]. The transcriptome expression pro les and copy number variation (CNV) of CXC chemokines, along with clinical data from the "TCGA-READ" cohort with 159 READ patients were extracted from the TCGA database. A total of 131 READ patients with full clinicopathological data were applied for further analyses. Co-expression analysis was presented. Prognostic analysis was performed using a Kaplan-Meier method.. The association between expression levels of CXC chemokines and the clinicopathological characteristics of READ patients was also analyzed by the Chi-square test. P < 0.05 was considered signi cant.

ONCOMINE
ONCOMINE (https://www.oncomine.org/) is an online microarray database for genome-wide expression analysis [66]. We compared the mRNA levels of CXC between READ tissues and corresponding normal tissues. The criteria for data ltration were set as a P < 0.05, a |fold change| > 2, and a gene rank in the top 10%. Student's t-test was used to analyze the difference in mRNA levels of CXC chemokines in READ samples.

GEPIA
GEPIA (http://gepia.cancer-pku.cn/index.html) is an interactive web application for gene expression analysis based on 9736 tumor samples and 8587 normal samples from the TCGA and the GTEx databases[67]. We used the "Multiple gene comparison" module in the "Expression DIY" to generate a relative expression heatmap of CXC chemokines in READ.

cBioPortal
The cBio Cancer Genomics Portal (cBioPortal, http://cbioportal.org) is an open-access resource for interactive exploration and visualization of multidimensional cancer genomics data, which provides access to data of more than 5,000 tumor samples from different cancer studies [68]. Genetic alteration analysis of CXCL3, CXCL12, and CXCL13 were presented using 155 rectal adenocarcinoma (READ) samples extracted from 594 colorectal adenocarcinoma samples (TCGA, PanCancer Atlas). The criterium of z-score for data ltration was set as ± 2.0. The top 50 frequently altered neighbor genes associated with CXCL3, CXCL12, and CXCL13 were also identi ed, in which 41 genes were applied for protein-protein interaction (PPI) and functional analyses as the other nine genes did not functionally relate to others.

Metascape
Metascape (http://metascape.org/gp/index.html) is an effective and e cient tool to analyze and interpret OMICs-based studies comprehensively [69]. In this study, the "Express Analysis" module was chosen. Metascape subdivided the groups and built a network according to the enriched pathways of CXCL3, CXCL12, CXCL13, and the 41 neighbor genes.

STRING
The STRING database (https://string-db.org/) aims to collect, score, and integrate all publicly available sources of PPI information, and to complement these with computational prediction [70]. The PPI networks of upregulated CXC chemokines, as well as frequently altered neighbor genes, were constructed using STRING. The con dence score of ≥ 0.4 was set as the threshold. All the networks were visualized by Cytoscape 3.7.0.
We presented the interactive functional association network predicted by GeneMANIA among CXCL3, CXCL12, CXCL13, and the 41 neighbor genes [72].

Competing endogenous RNA (ceRNA) network construction
The ceRNA network construction of CXCL3 and CXCL13 was presented based on the theory of ceRNA [73]. The construction process included two steps. First, miRNA-mRNA interactions of CXCL3 and CXCL13 were predicted respectively by DIANA-MicroT-CDS (http://www.microrna.gr/microT-CDS) [74]. The threshold was set as the miTG-score > 0.9. After that, we used the selected miRNA to predict miRNA-lncRNA interactions by DIANA-LncBase Predicted v.2 (www.microrna.gr/LncBase) [75]. The criteria were set as the interaction score = 1.000. For each network, we only selected the long non-coding RNAs(lncRNAs) that interacted with at least two of the predicted miRNAs above for construction. Finally, the two ceRNA networks were visualized by R 3.6.2.

Gene Set Enrichment Analysis
Gene Set Enrichment Analysis (GSEA) that evaluates microarray data at the level of gene sets can be used to link prior knowledge to newly generated data and thereby help uncover the collective behavior of genes in states of health and disease [76]. We used GSEA for Gene Ontology (GO) enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways analysis of CXCL3 and CXCL13 utilizing the TCGA-READ dataset. P < 0.05 and FDR q < 0.25 were considered signi cant. The results were presented using R 3.6.2.

Declarations
Ethics approval and consent to participate Not applicable.

Consent for publication
Not applicable.

Availability of data and materials
The datasets used and analyzed during the current study are available from the corresponding author on reasonable request.

Competing interests
The authors declare that they have no competing interests.    The transcription levels of CXC chemokines in READ. The expression levels of (A) CXCL3, (B) CXCL12, (C) CXCL13, and (D) CXCL16 in READ tissues and normal rectum tissues generated from the TCGA database were shown. Data were analyzed by Student's t-test and P < 0.05 was considered signi cant. * P < 0.05, ** Page 17/22 P < 0.01, *** P < 0.001. (E) Expression patterns of CXC chemokines generated from TCGA database and ONCOMINE database, respectively. READ: Rectal adenocarcinoma.