Identification of neoadjuvant radiotherapy resistance genes in rectal cancer using integrated bioinformatics analysis

The incidence of rectal cancer is increasing year by year. And most patients have been diagnosed in the advanced stages largely because of the special anatomical location of the rectum, resulting in uneasy diagnosis in its early stage. The combination of radiotherapy and chemotherapy is widely accepted as effective treatment to the rectal cancer, significantly reduce the local recurrence rate and improve long-term survival, but the occurrence of radiotherapy resistance of rectal cancer strikingly weaken the responding of treatment. So finding reasons or factors leading to the resistance of neoadjuvant radiotherapy for rectal cancer become very urgent not just for the current clinic but also for development of new therapeutical approaches. To identify potential biomarkers associated with neoadjuvant radiotherapy resistance in rectal cancer using bioinformatics analysis. Firstly, we extracted relevant data from GEO database and analyzed the corresponding data by a series of analytic methods. Then, the candidate genes were identified by performing of the analysis of protein-protein interaction network, KEGG and GO, subsequently comfirmed at the levels of mRNA and Protein with Oncomine and The Human Protein Atlas respectively. Furthermore, these comfirmed genes were submitted to CMap analysis to identify candidate therapeutic compounds for neoadjuvant radiotherapy resistance. We found that the abnormal expression of some chemokines were tightly associated with the rectal cancer neoadjuvant radiotherapy resistance. Especially the upreuglation of Chemokines 9(CXCL9) and Chemokines 10(CXCL10) and downregulation of Chemokines 13(CXCL13) and Chemokines 5(CCL5) might be the potential biomarks of the radiotherapy resistance of rectal cancer. Through a series of analyses, the molecules that may be related to the radiotherapy resistance of rectal cancer were identified that may have a hint effect on the treatment of radiotherapy resistance of rectal cancer. proteolysis, positive regulation of osteoblast proliferation, cytoskeletal anchoring at magnesium ion transport, digestion, complement activation, classical pathway, magnesium ion transmembrane transport, chemokine-mediated signaling pathway, D-amino acid metabolic process, folic acid metabolic process, drug transmembrane transport, muscle cell cellular homeostasis, error-prone translesion synthesis, neuron projection morphogenesis, proteolysis involved in peroxisome, photoreceptor CXCR3 chemokine receptor binding, receptor binding, antigen binding, chemokine activity, PDZ domain binding, electron carrier activity, cofactor binding, magnesium ion transmembrane transporter activity, D-amino-acid oxidase activity, receptor signaling protein tyrosine kinase activator activity, repair; insulin secretion; cytokine-cytokine interactions; protein digestion and


Introduction
The advantages of preoperative neoadjuvant radiotherapy include a good blood supply of target cells, high sensitivity to radiation, and efficacy. Radiotherapy can shrink tumors, leading to significantly improving the curative tumor resection and anal preservation rates [1,2]. Therefore, neoadjuvant radiotherapy has become the standard of care for patients with resectable rectal cancer. But, because of radiotherapy resistance, a considerable proportion of patients were unable to be benefit from neoadjuvant radiotherapy. Hence, the clinical response of neoadjuvant radiotherapy needs further evaluation with respect to resistance genes and mechanisms in order to help physicians prepare more effective therapeutic approaches.
Until recently, some progresses were obtained in the field of neoadjuvant radiotherapy resistance. For example, MiR-17-5p has been reported to significantly enhance the sensitivity of X-rays to radioresistant cells as well as downregulate and promote antiradioactive phenotypes in esophageal adenocarcinoma cancer stem-like cells [3].
Overexpression of TCN1 is reportedly associated with adverse treatment responses and outcomes in rectal cancer patients receiving concurrent neoadjuvant chemoradiotherapy, indicating the potential prognostic value of this gene [4]. Elevated levels of β-catenin have been shown to inhibit proliferation of colorectal cancer cells in vitro and growth of mouse xenograft tumors [5], and heightened expression of Rhoa has been shown to be a good predictor of poor prognosis due to high treatment resistance rates [6]. Nonetheless, specific mechanism(s) associated with newly developed radiation resistance are not fully understood as it is a complex process involving a large number of genes.
Rapid development of high-throughput biotechnologies, such as microarray, has made it possible to analyze RNA expression profiles, thereby enhancing our understanding of various diseases and their molecular mechanisms [7][8][9]. However, few studies have systematically analyzed gene expression profiles of neoadjuvant radiotherapy resistance by integrating bioinformatics analysis. In the present study, neoadjuvant radiotherapy resistance gene expression datasets were obtained from Gene Expression Omnibus and systematic bioinformatics analysis performed, including identification of differentially expressed genes (DEGs), functional enrichment and co-expression network analyses, and identification of significantly differently expressed biomarkers of neoadjuvant radiotherapy resistance in colorectal cancer. The results of the present study will significantly improve our understanding of radiotherapy resistance mechanisms and efficacy.

Identification of DEGs
The GSE60331 date set was obtained from the Gene Expression Omnibus database.
Analysis of radiotherapy-resistant group (8 radiotherapy-resistant rectal tumor group, 8 radiotherapy-resistant rectal normal group samples), and radiotherapy-sensitive group (6 radiotherapy-sensitive rectal tumor group, 6 radiotherapy-sensitive rectal normal group) obtained 1325 and 1247 related genes, respectively. A total of 1000 DEGs were common to the radiotherapy-resistant and radiotherapy-sensitive groups. These DEGs were thought to be common in rectal cancer and related to its occurrence, but not necessarily all specifically related to neoadjuvant radiotherapy resistance, so we ignore these genes. In the radiotherapy-resistant group, we found 325 unique DEGs thought to be related to neoadjuvant radiotherapy resistance in rectal cancer and 247 in the radiotherapysensitive group (Fig. 1A). Genes with corrected P < 0.05 and an absolute fold change > 2 were considered as DEGs after the original data was corrected by Limma (Fig. 1B). The specific resistance genes founded in the radiotherapy-resistant group showed that normal and cancerous tissues could be separated by cluster analysis (Fig. 1C); the left 8 samples were in the radiotherapy-resistant rectal tumor group, while the right 8 samples were in the radiotherapy-resistant rectal normal group. This result confirmed the reliability of the analysis. Based on this method, we found that 325 genes between the neoadjuvant radiotherapy-resistant and radiotherapy-sensitive groups were significantly differentially expressed [169 up-regulated, 156 downregulated] (Fig. 1D).

KEGG and GO Enrichment
To get deeper information about DEG function, GO enrichment on the above 325 DEGs was carried out with the aid of the Database for Annotation, Visualization, and Integrated Discovery with a P < 0.05. GO results showed that 20 pathways were enriched ( Fig. 2A), mainly including metabolic pathways; cellular response to zinc ion, negative regulation of growth, cellular response to cadmium ion, drug transport, endonucleolytic cleavage in ITS1 to separate SSU-rRNA from 5.8S rRNA and LSU-rRNA from tricistronic rRNA transcript (SSU-rRNA, 5.8S rRNA, LSU-rRNA), proteolysis, positive regulation of osteoblast proliferation, cytoskeletal anchoring at plasma membrane, magnesium ion transport, digestion, complement activation, classical pathway, magnesium ion transmembrane transport, chemokine-mediated signaling pathway, D-amino acid metabolic process, folic acid metabolic process, drug transmembrane transport, muscle cell cellular homeostasis, error-prone translesion synthesis, neuron projection morphogenesis, proteolysis involved in cellular protein catabolic process, extracellular exosome, peroxisomal matrix, peroxisome, extracellular space, cell, external side of plasma membrane, photoreceptor outer segment membrane, extracellular region, CXCR3 chemokine receptor binding, receptor binding, antigen binding, chemokine activity, PDZ domain binding, electron carrier activity, cofactor binding, magnesium ion transmembrane transporter activity, Damino-acid oxidase activity, receptor signaling protein tyrosine kinase activator activity, 1-alkylglycerophosphocholine O-acetyltransferase activity. The top 10 results of the GO enrichment analysis are shown in Table 1. Moreover, we also found DEGs enriched in different GO terms. Twenty pathways were enriched in biological processes, 8 were enriched in cellular components, and 11 were enriched in molecular unction (P < 0.05).
Then we classify these GO terms according to biological processes, cell component and  Table 2.

PPI network construction and analysis
The PPI relationships between DEGs were determined using the Search Tool for the Retrieval of Interacting Genes/Proteins, then a PPI network was constructed with Cytoscape. Finally, a network with 51 nodes and 116 edges was built (Fig. 3A). The top 18 degree hub nodes in the PPI network were DCAF13, BRIX1, NMU, RRS1, CIRH1A, PNO1, CXCL10, CXCL9, CCL5, NPY, POLRIE, DDX27, RCL1, GRM8, CXCL13, BMP2, DPP4, and POLR1D. These genes/proteins might play an important role in neoadjuvant radiotherapy resistance of rectal cancer (Fig. 3B). Importantly, 4 genes (CXCL9, CXCL10, CXCL13, and CCL5) were found to be enriched in each analysis. Therefore, it was postulated that these enrichment scores. As a result we identified top 7 negatively related small molecules drugs were listed in Table 3 and PCMap-value of < 0.0001 (Table 3).

Discussion
Radiation tolerance is a major factor related to failure of radiotherapy and poor prognosis for cancer patients. Therefore, radiotherapy resistance has become a major obstacle to overcome [10,11]. Radiotherapy resistance gene chip analysis of rectal cancer has the advantage of being able to provide a much larger amount of data than other previous biotechnology methods, allowing for more extensive experimental design and analytical scope [12,13]. In the present study, an integrated bioinformatics analysis approach was used to identify potential biomarkers of neoadjuvant radiotherapy resistance in rectal cancer. Ultimately, the 4 differentially expressed chemokine genes (CXCL 9, CXCL10, CXL13, and CCL5) were identified as being related to radiotherapy resistance in colorectal cancer.
These results are in agreement with previous studies reporting a close relationship between chemokines and neoadjuvant radiochemotherapeutic efficacy [14]. CXCR4, CXCL12, and CXCR7 have been suggested to be able predict neoadjuvant radiochemotherapy and radiotherapy resistance in T-and N-state rectal tumors, with high CXCR4 expression indicating a poor prognosis [15]. Moreover, upregulation of CXCL12 (Sdf-1a) has been shown to induce neoadjuvant radiotherapy resistance [16]. CXCL9, CXCL10, and CXCL11 proteins have also been reported to inhibit the formation of tumor blood vessels through their coligands; hence, overexpression of CXCL9, CXCL10, and CXCL11 can inhibit the growth and metastasis of colorectal cancer [17]. Thus, we concluded that abnormal expression of chemokines maybe the key factors that leading to occurrence of radiotherapy resistance.
Herein, high CXCL9 and CXCL10 expression was found to be related to radiotherapy resistance in rectal cancer cells. This is in agreement with previous studies showing that high CXCL10 expression inhibits growth of certain cancer cells [18], and high CXCL10, CXCL8, and CXCL11 expression lead to a decrease in the survival rate of breast cancer patients [19]. Moreover, overexpression of CXCL9 and CXCL10 has been related to a high degree of tumor invasion [20]. CXCR3directed chemokines CXCL9 and CXCL10 induce Cathepsin B in human breast cancer cells, thereby supporting tumor invasion and metastasis [21].
The results of the present study also predicted that low CXCL13 and CCL5 expression is related to radiotherapy resistance in rectal cancer cells. A previous study reported that transforming growth factor-β signaling and Oct4 synergistically induce expression of the epithelial-to-mesenchymal transition-related genes Snail, Slug, and CXCL13, which accelerates disease progression, especially in late stages, and may indicate a poor prognosis in breast cancer patients [22]. Overexpression of CCL5 has also been shown to lead to migration and proliferation of breast cancer cells [23]. CXCL13 has also been shown to participate in the migration and invasion of androgen-induced prostate cancer due to upregulation by the androgen receptor axis at mRNA and protein levels. Thus, CXCL13 is considered to be an androgen-responsive gene in prostate cancer patients with higher than normal androgen levels [24]. In contrast, PKCε overexpression and PTEN loss have also been reported to cooperatively upregulate production of CXCL13, inhibiting the migration and tumorigenicity of prostate cancer [25]. At the same time, low CXCL13 expression in gastric cancer patients has been associated with a better prognosis [26].
However, blocking the CCR5/CCL5 axis may be prevent metastasis and provide more therapeutic strategies to control the progression of pancreatic cancer [27]. Nevertheless, it is obvious that these chemokines (CXCL9, CXCL10, CXCL13, and CCL5) have certain influences on the occurrence and development of tumors.
Neoadjuvant radiotherapy kills tumor cells through the cytotoxic action of radiation and drugs and its effects on the surrounding microenvironment. The tumor microenvironment is composed of a variety of components, including chemokines, plays an important role in the tumor progression and pathophysiological processes [28]. Chemokine receptors and ligands are abundantly expressed in inflammatory environments and increasing evidence indicates that chemokines are involved in the regulation of tumor proliferation, angiogenesis, apoptosis, and metastasis in cancer-related inflammation and can be used as therapeutic targets. CXCR4 has been reported to be upregulated by the microenvironment in colon cancer, and CXCR4 inhibitors have the potential to inhibit metastatic growth of tumors [29]. CXCR2/CXCL7 have also shown prognostic value in colorectal cancer. Colorectal cancer cells secrete interleukin-8, which mediates the homing of mesenchymal stem cells and induces CXCL7 production. The disease-free and overall survival of patients with high CXCR2/CXCL7 expression in colon cancer metastases are both shortened [30]. Moreover, overexpression of CCR6 has been shown to be a key factor in promoting the metastasis of colorectal cancer cells in vivo, especially to the liver [14]. The concentration gradients of CCR6 and CCL20 in colorectal cancer provide a pathway for the specific transfer (metastasis) of colorectal cancer cells to the liver. In breast cancer cells, transfection of CXCL7 stimulates Matrigel invasion and secretion of the lymphangiogenic factors vascular endothelial growth factor-C and -D [31,32].
Although the results of our analysis are not completely consistent with the outcomes of the above studies, we suspect this may be due to a decrease in the expression of CXCL13 and CCL5, which leads to changes in the tumor microenvironment. However, such mechanisms of action remain to be explored. microenvironment. Therefore, Therefore, we speculate that chemokines may have a certain effect on radiotherapy for rectal cancer. This will be of some significance for us to explore their molecular mechanisms.

Data sources
The GSE60331 gene expression profile was obtained from the Gene Expression Omnibus database (http://www.ncbi.nlm.nih.gov/geo/) of the National Center for Biotechnology Information; this platform was applied to expression arrays. Data were analyzed from a total of 28 tissue samples include 8 radiotherapy-resistant rectal tumor group, 8 radiotherapy-resistant rectal normal group, 6 radiotherapy-sensitive rectal tumor group, 6 radiotherapy-sensitive rectal normal group, the criteria were data from tumor samples and normal tissue; biopsies were taken at two different timepoints: the first time is when tumor and normal mucosa before treatment and the second time is only tumor (week 3) that loading dose of bevacizumab but before chemoradiation. CEL forms and annotation files were downloaded for further analysis.

Identification of DEGs
Microarray data was preprocessed in Bioconductor with the Affy software package [33] of R (version 3.3.2) software. The robust multivariate average algorithm [34] was used for background correction, quartile normalization, and probe summarization to obtain gene expression matrices. The data were normalized using the Limma package of R software and then by the Bayes method [35]. DEGs between cancerous and normal tissues in radiotherapy-resistant and -sensitive groups were screened separately using a P ≤ 0.05 and |log(Fold Change)| > 1 as filter conditions. Then, DEGs in both groups were merge, and the duplicates were removed. Finally, specific DEGs associated with colon cancer radiotherapy resistance or sensitivity were obtained. The original expression datasets under all conditions were normalized with the robust multiarray average method [36] using the ClusterProfiler package in R software to obtain the gene expression matrix. The t-test method in the R software Limma package was used to identify DEGs between radiotherapy-resistant and -sensitive group genes. Values with a |log(Fold Change)| > 2.0 and P < 0.05 were selected as cutoff criteria.

Functional enrichment analysis
In order to better understand the functions of the selected DEGs, Gene Ontology [GO] [37] functional and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed for radiotherapy-resistant DEGs to identify signaling pathways potentially involved in the development of radiotherapy resistance in rectal cancer. The Database for Annotation, Visualization, and Integrated Discovery (https://david.ncifcrf.gov) was used for GO enrichment of DEGs, and the KEGG Automatic Annotation Server (http://www.genome.jp/tools/kaas/) was used for KEGG pathway enrichment [37,38] with a P < 0.05 and number of genes > 2. The ggplot2 package in R software was used to create GO and KEGG histograms. Then we draw the column diagram through barplot.

Analysis of DEGs by CMap
The DEGs were converted into probes. These probes were collated into tag sets of "up" or "down" genes, and analyzed using the CMap website (https://portals.broadinstitute.org/cmap/) [41]. The similarity was reported as connectivity score ranging from + 1 to -1. Positive connectivity score denotes increased similarity while negative score denotes inverse similarity between the expression patterns induced by the compounds in the CMap database and the query DEGs.     CXCL9: X axis label 1 is normal of rectal tissue. X axis label 2 is tumor of reatal tissue. (B right) CXCL10: X axis label 1 is normal of rectal tissue. X axis label 2 is tumor of reatal tissue. (C right) CXCL13: X axis label 1 is normal of ascending colon. X axis label 2 is normal of descend colon. X axis label 3 is normal of reatal.

Declarations
X axis label 4 is normal of transverse colon, X axis label 5 is normal of sigmoid colon, X axis label 6 is tumor of colorectal. (D right) CCL5: X axis label 1 is normal of rectal tissue. X axis label 2 is tumor of reatal tissue.