Upregulation of STAT1-CCL5 axis is a biomarker of colon cancer and promotes the proliferation of colon cancer cells.

BACKGROUND: Colorectal cancer (CRC) is the third most commonly diagnosed cancer in men and women globally. Investigating genetic ground differences between normal and CRC tissues would be significant for identifying some key oncogenic pathways and developing anti-cancer agents.
METHODS: Weighted gene co-expression network analysis (WGCNA) method was used to screen out core pathways related to the clinical traits of CRC patients. Then, multiple databases were utilized to further verify the hub genes obtained from data mining. Finally, to explore the role of hub genes in CRC, cell counting and EdU assays were performed.
RESULTS: The results of the WGCNA analysis showed that a module (turquoise module) was highly related with CRC differentiation grade (R =0.53, P<0.0001). Enrichment analysis indicated that genes of the turquoise module were remarkably enriched in multiple inflammatory processes and pathways. Among all hub genes of the turquoise module, the mRNA levels of STAT1 and CCL5 were significantly higher in CRC than in normal colon tissues. STAT1 expression was highly positively correlated with the level of CCL5. The results of the cell counting, EdU, CCK-8, and CFSE staining assays showed that interfering with STAT1 and CCL5 could inhibit the proliferation of CRC cells.
CONCLUSIONS: Our study indicated that the STAT1-CCL5 axis is an important modulator in the development of CRC through promoting cell proliferation. Moreover, the levels of STAT1 and CCL5 might be valuable biomarkers for CRC screening. 2020 Annals of Translational Medicine. All rights reserved.

Introduction

Colorectal cancer (CRC) was the third most commonly diagnosed cancer in men and women in 2019 (1). The statistical data of the American Cancer Society showed that CRC is the second leading cause of cancer-related death in the United States (1). Improvements in cancer prevention, early diagnosis, and treatment strategies leaded to a 35% decrease in CRC-associated mortality between 1990 and 2007 (2). Encouragingly, the advancements in precision medicine have further improved the outcomes of patients without compromising treatment effects. Nevertheless, a recent report presented an unsatisfactory 5-year survival rate for stage IV CRC patients of about 14% (3). With the developments of gene sequencing techniques, there is a growing body of evidence that indicates that some genetic alterations, such as those in KRAS, NRAS, epidermal growth factor receptor (EGFR), as well as APC, are closely related to carcinogenesis and poor clinical outcomes in CRC (4). For colon cancer patients, these gene expression statuses not only predict prognosis but also herald the efficacy of some regimens such as bevacizumab, cetuximab, and panitumumab (5). Therefore, investigating genic background differences between normal and colon cancer cells would be significant to identifying some key oncogenic pathways and understanding the malignant biological behaviors of cancer cells, and may also provide valuable information for the development of anti-cancer drugs.

In the present study, we calculated the correlation values of genes through a WGCNA-based systemic biology method and identified a novel signaling pathway associated with the development of colon cancer and the prognosis of patients. This study was based on dataset GSE29612. The raw gene expression data of GSE29621 were provided by Chen et al., and the pre-processed gene matrix was downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo) (6). GSE29612 contains the expression profiles of 65 colon cancer samples (tumor purity >90%, without genetic cancer syndromes and inflammatory bowel disease) (6).

Methods

Data cleaning and processing

shows the workflow followed in this study. The gene expression values were detected with the Affymetrix Genechip Human Genome U133 Plus 2.0 Array and were further transformed using the Robust Multi-array Average (RMA) algorithm. Then, based on the results of variance analysis, the top 50% of the variant genes (10,826 genes) were used in the following WGCNA analysis.

Figure 1

The work flow of the study.

Constructing co-expression network

To construct the co-expression network, the gene expression matrices of 10,826 genes were extracted by WGCNA package (version of R software: 3.6.0). The analysis was performed as previously described (7). The connection strength between two nodes was calculated by adjacency matrix aij using the following formula:

In the formula, x and x represent the vectors of the expression data of genes i and j. S represents the Pearson’s correlation coefficients of x and x. a reflects the network connection strength between two genes. In this study, β =5 was served as the soft-threshold value for constructing the scale-free network (scale-free R2 >0.85). In WGCNA analysis, the correlations between two connected genes as well as the correlations between multiple associated genes are vital parameters. A series of modules were generated by hierarchical clustering of weighted coefficient matrix. Then, topological overlap measure (TOM) was calculated using the following formula:

In this formula, A represents weighted adjacency matrix calculated by A. TOM-based dissimilarity method (minimized gene group =50) was employed to perform average linkage hierarchical clustering. The DynamicTreeCut algorithm was used to distribute genes with similar expression patterns into the same module.

Identifying clinically significant modules

Two parameters, module gene significance (MS) and module eigengenes (MEs), were adopted to evaluate the correlations between the clinical traits of patients and clustered modules. As the first principal component, the expression levels of all the genes in one clustered module could be profiled with ME. Gene significance (GS) refers to the absolute values of P in the linear regression analysis. MS was calculated based on the average GS values of all the genes in each module.

KEGG pathway and GO term enrichment analysis

As previously described, the enrichment analysis was conducted using the online analysis tool G: Profiler (https://biit.cs.ut.ee/gprofiler/) (version: e98_eg45_p14_ce5b097). The significantly enriched KEGG pathways and GO terms were identified with <0.05 serving as the cutoff for determining the false positive rate.

Hub genes

Hub genes were screened out by calculating their connectivity with other genes (|cor.geneModuleMembership| ≥0.8). After that, the hub genes had to meet the criterion |cor.geneTraitSignificance| ≥0.2. To further investigate the role of hub genes in CRC, the gene expression matrix was downloaded from the GEO, the Cancer Cell Line Encyclopedia, The Cancer Genome Atlas (TCGA), and the Genotype-Tissue Expression (GTEx) databases. Comparisons of hub gene expression levels between tumor and normal tissues were performed using the online Gene Expression Profiling Interactive Analysis (GEPIA) tool (http://gepia.cancer-pku.cn/).

Cell counting

To evaluate the effect of two hub genes, STAT1 and CCL5, on colon cancer cells, colon cancer cell lines HCT-116, SW480, and SW620 were cultured in high-glucose Dulbecco’s Modified Eagle Medium (DMEM) supplemented with 10% or 1% fetal bovine serum (FBS) at 37 °C with 5% CO2. After the cells were incubated with STAT1 inhibitor Fludarabine (10 µM, Catalog No: HY-B0069, MCE) (8) and CCL5 neutralizing antibody (1 µg/mL, Catalog No: 12000-1-AP, Proteintech), the number of cells of the three cancer cell lines were counted with a cell counting chamber.

EdU assay

An EdU assay was conducted according to the manufacturer’s instructions (Cell-LightTM EdU Apollo 567 Kit, Catalog No: C10310-1, RIBOBIO) (9). Firstly, cells were seeded in 24-well plates at a density of 3×103 per well and incubated for 2 days. Next, the cells were stimulated for 12 h with or without 10% FBS. Then, 250 µL EdU medium was added to each well. After incubation lasting 2 h, the cells were fixed with 4% paraformaldehyde for 20 min. Subsequently, glycine solution was added to neutralize aldehydes. After the permeability of the cell membrane had been increased with 0.5% TritonX-100 in PBS, the cells were incubated with Apollo dyeing reaction solution in the dark for 30 min. Nuclei were counterstained with Hoechst 33342. Finally, fluorescent images were taken at 200× and the positive cells were counted in 5 fields per well.

Small-interfering RNA (Si-RNA) transfection

Si-RNA of STAT1 was purchased from RiboBio Co., Ltd. (Guangzhou). Tumor cells were seeded at a density of 2×106 cells in a 6-well plate on the day prior to transfection. The transfection was performed using Hieff Trans™ in vitro siRNA/miRNA Transfection Reagent (Yeasen Biotech Co., Ltd.) according to the manufacturer’s recommendations. Twenty picomole Si-RNA was required for each well.

Real-time PCR analysis

The primers of STAT1 were designed as follows: forward primer, 5'-CAGCTTGACTCAAAATTCCTGGA-3'; reverse primer, 5'- TGAAGATTACGCTTGCTTTTCCT-3'. The primers of CCL5 were designed as follows: forward primer, 5'-CCAGCAGTCGTCTTTGTCAC-3'; reverse primer, 5'-CTCTGGGTTGGCACACACTT-3'. Total mRNA was extracted using FastPure Cell/Tissue Total RNA Isolation Kit (Vazyme) according to the manufacturer’s protocol. Approximately 500 ng total mRNA was reversely transcribed to complementary DNA (cDNA) using HiScript II Q Select RT SuperMix for qPCR (Vazyme). Quantification reactions were performed with ABI 7500 platform. Samples were amplified using the following conditions: 95 °C for 30 s, 45 cycles of 95 °C for 15 s, 60 °C for 60 s. The expression fold changes were calculated with the formula: 2−∆∆CT (10).

CCK-8 assay

Cell viability was measured using CCK-8 cell proliferation kit (Dojindo). Ten µL CCK-8 reagent was added into 100 µL cell culture medium. After incubation for 2.5 hours in 37 °C, the absorbance values were measured at 450 nm with the microplate reader (Molecular Devices).

CFSE staining and flow cytometry

SW480 and SW620 were labeled with 2.5 µM CFSE (Invitrogen). Then, cells were seeded at the density of 1 × 106 per well in 6-well plates and cultured with RPMI-1640 media containing 10% FBS. After culture for 2 days, we performed flow cytometry assay to measure the ratio of proliferating cells (CFSE-low) with BD Celesta platform. The data of flow cytometry were analyzed with Flow Jo software (Flow Jo LLC, Ashland, USA).

Results

Constructing the co-expression network and identifying clinically significant modules

Patients’ clinical traits were summarized in . After Pearson’s correlation coefficients and average strength were calculated, the 65 samples of GSE29621 were clustered (). Then, with the soft-thresholding power as β =5 (R2 >0.85), the co-expression analysis was conducted by generating a scale-free network. Based on an average hierarchical clustering algorithm, 26 modules were generated (). Notably, the turquoise module possessed a strong correlation with CRC differentiation grade (R =0.53, P<0.0001) (, ). Thus, the turquoise module was identified as having clinical significance.

Table 1

The clinical traits of patients in GSE29621

Clinical traits	Group	Values
Gender	Male	40
	Female	25
Tumor size	T2	8
	T3	52
	T4	5
Lymph node metastasis	N0	32
	N+	32
	NA	1
Distant metastasis	M0	46
	M1	18
	NA	1
AJCC stage	Stage I–II	29
	Stage III–IV	36
Vital status	Alive	40
	Dead	25
Overall survival time		45.87 (median value)

Figure 2

The clustering dendrogram of 65 CRC samples. CRC, colorectal cancer.

Figure 3

Calculating soft-thresholding power. (A) Calculating scale-free fit indices using different soft-thresholding powers (β). (B) Calculating mean connectivity using different soft-thresholding powers.

Figure 4

Identifying clinically significant modules of GSE29621. (A) The dendrogram of the most variant genes clustered by Dynamic Tree Cut algorithm. (B) The heat map profiling the correlations between all module eigengenes and the clinical traits of patients including gender, adjuvant chemotherapy history, tumor differentiation grade, tumor size, tumor stage, overall survival, and death event.

Table 2

The correlations between the turquoise module and patients’ clinical traits

Module	Traits	R	P values
turquoise	Gender	−0.19824	0.113414
turquoise	Adjuvant chemotherapy	−0.15069	0.23084
turquoise	Histologic grade	0.519138	9.39E-06
turquoise	Tumor size	0.155133	0.217226
turquoise	AJCC stage	−0.17554	0.1619
turquoise	Overall survival time	0.127409	0.311825
turquoise	Death event	−0.15132	0.228881

The results of KEGG pathway and GO term enrichment analyses

The turquoise module was a highly complicated gene panel comprising 1,723 genes. Therefore, the biological functions of the turquoise module were profiled through KEGG pathway analysis and GO term enrichment analysis. The results of enrichment analysis showed that the genes in the turquoise module were significantly enriched in inflammation-associated pathways such as the B-cell receptor signaling pathway, antigen processing and presentation, the chemokine signaling pathway, cytokine-cytokine receptor interaction, the NF-κB pathway, PD-L1 expression, and the PD-1 pathway in cancer (). Furthermore, this panel of genes was remarkably enriched in extracellular cytokine release processes including in the cytokine-mediated signaling pathway, extracellular regions, and chemokine receptor activity. Cytoscape software (software version: 3.6.0) was employed to construct the interaction networks of GO terms ().

Figure 5

KEGG pathway enrichment analysis of genes in the turquoise module. (A) The histogram showing the significantly enriched representative KEGG pathways. (B) The bubble diagram showing the significantly enriched representative KEGG pathways.

Figure 6

The interaction networks of significantly enriched GO biological process and GO cellular component terms. (A) The interaction networks of significantly enriched GO BP terms. (B) The interaction networks of significantly enriched GO CC terms. BP, biological process; CC, cellular component.

Figure 7

The interaction networks of enriched GO molecular function terms.

The identification and verification of hub genes

A total of 82 genes in the turquoise module met the criteria of |cor.geneModuleMembership| ≥0.8 and |cor.geneTraitSignificance| ≥0.2 and were subsequently screened out as hub genes. Among these 82 genes, two genes, STAT1 and CCL5, stood out. The results from the TCGA and GTEx databases showed that the mRNA levels of STAT1 and CCL5 were significantly elevated in CRC tissues compared with normal colon tissues (). In addition, the correlation analysis of the TGCA, GEO, and Cancer Cell Line Encyclopedia databases indicated that the expression of STAT1 was closely correlated with that of CCL5 (P<0.0001) (). Moreover, in SW489 and SW620 cells, when STAT1 expression was inhibited by Si-STAT1, the level of CCL5 was decreased simultaneously ().

Figure 8

The mRNA levels of STAT1 and CCL5 in CRC and normal colon tissues. (A) The mRNA levels of STAT1 in CRC and normal colon tissues. (B) The mRNA levels of CCL5 in CRC and normal colon tissues. This analysis was performed using the online tool GEPIA. CRC, colorectal cancer.

Figure 9

The correlations analysis between STAT1 and CCL5 mRNA levels. (A) The correlations analysis between STAT1 and CCL5 mRNA levels in CRC tissues based on the GSE29621 dataset. (B) The correlations analysis between STAT1 and CCL5 mRNA levels in CRC tissues based on the TCGA database. (C) The correlations analysis between STAT1 and CCL5 mRNA levels in CRC tissues based on the GSE28722 dataset. (D) The correlations analysis between STAT1 and CCL5 mRNA levels in colon cancer cell lines based on the Cancer Cell Line Encyclopedia database. (E) The histogram reflecting the effect of Si-STAT1 on CCL5 expression in SW480 cells. (F) The histogram reflecting the effect of Si-STAT1 on CCL5 expression in SW620 cells. *, P<0.05; ****, P<0.0001. CRC, colorectal cancer.

STAT1-CCL5 pathway promoting the proliferation of colon cancer cells

For the cancer cells cultured with medium containing 10% FBS, STAT1 inhibitor and CCL5-neutralizing antibody exhibited anti-proliferative activity. After the HCT-116, SW480, and SW620 cells underwent STAT1 inhibitor treatment at 10 µM, a significant decrease was observed in their numbers. Meanwhile, CCL5-neutralizing antibody at 1 µg/mL also decreased the number of cancer cells in the culture medium with 10% FBS (a decreasing trend was observed, but statistical significance was not reached in the HCT-116 cells) (). The results of the EdU assay showed that treatment with STAT1 inhibitor remarkably reduced the number of EdU-positive cancer cells (). Moreover, we used medium with 1% FBS to minimize the potential influence of growth factors from serum. Similar to the findings under the conditions of 10% FBS, blocking STAT1 and CCL5 inhibited the proliferation of three cancer cell types.

Figure 10

Cell counting assays. (A) The histogram reflecting the effect of STAT1 inhibitor on the proliferation of colon cancer cell lines cultured with 10% FBS. (B) The histogram reflecting the effect of STAT1 inhibitor on the proliferation of colon cancer cell lines cultured with 1% FBS. (C) The histogram reflecting the effect of CCL neutralizing antibody on the proliferation of colon cancer cell lines cultured with 10% FBS. (D) The histogram reflecting the effect of CCL neutralizing antibody on the proliferation of colon cancer cell lines cultured with 1% FBS. *, P<0.05; **, P<0.01; ****, P<0.0001.

Figure 11

EdU assays. (A) The representative images and relative quantities of EdU positive colon cancer cells treated with STAT1 inhibitor or DMSO cultured with 10% FBS. (B) The representative images and relative quantities of EdU-positive colon cancer cells treated with STAT1 inhibitor or DMSO cultured with 1% FBS. In this assay, proliferative cells could be dyed red after the reaction of EdU reagent and fluorochrome. The cell nuclei were counterstained with Hoechst 33342 (blue). FBS, fetal bovine serum. ****, P<0.0001.

In addition, we used Si-STAT1 to explore the effect of STAT1 on colon cancer cell’s proliferation in SW480 and SW620 cells. The results of CFSE staining assay indicated that Si-STAT1 treatment significantly reduced the ratio of proliferating (CFSE-low) colon cancer cells (). The results of CCK-8 and cell number counting assays also demonstrated that Si-STAT1 impaired the proliferation of colon cancer cells ().

Figure 12

The effect of Si-STAT1 on the proliferation of colon cancer cell lines. (A) The representative images CFSE-labeled colon cancer cells treated with Si-STAT1 or Si-NC in flow cytometry assay. (B) The ratio of proliferating (CFSE-low) cancer cells treated with Si-STAT1 or Si-NC. (C) The histogram reflecting the relative cell counts of colon cancer cells treated with Si-STAT1 or Si-NC. (D) The results of CCK-8 assay reflecting the effect of Si-STAT1 on the proliferation of colon cancer cell lines. *, P<0.05; **, P<0.01; ***, P<0.001; ****, P<0.0001.

Discussion

Both sporadic CRC and colitis-associated cancer are well known to be related to inflammation. At some stages of cancer, the localized inflammation not only counteracts the anti-tumor immune response through an immunosuppressive effect, but it also participates in the proliferation, survival, and metastasis of CRC cells (11). Before the tumor lesion forms, cytokines from inflammatory cells may upregulate the generation of reactive nitrogen intermediates (RNI) and reactive oxygen species (ROS) (12). Terzic et al. found that intracellular RNI and ROS in pre-malignant cells silenced tumor suppressors through epigenetic alterations and promoted tumor initiation (11). Their study also showed that, apart from inducing the initiation of cancer cells, multiple pro-tumor cytokines could activate the receptors on colon cancer cells, which subsequently activated intracellular oncogenic signaling pathways or transcriptional factors such as NF-κB, AP-1, and STAT3 (11). These oncogenic signals regulate a wide range of biological processes including cell proliferation, survival, and apoptosis (13,14). Furthermore, Li et al. also found that the activation of these oncogenic signals enhanced chemokine secretion, recruited more inflammatory cells, and self-propagated localized inflammation (15). In advanced colon cancer, some proinflammatory cytokines including IL-1, TNF-α, and IL-6 activate transcriptional factors such as ZEB, Twist, and Snail, which promote epithelial-mesenchymal transition and increase the migratory properties of cancer cells (11). In this study, we comprehensively analyzed gene expression data from the TCGA, GTEx, GEO, and Cancer Cell Line Encyclopedia (CCLE) databases. We found that several tumor immune microenvironment-related signals played a key role in CRC.

The present WGCNA analysis was performed based on the GSE29621 dataset with the aim of identifying the co-expression modules associated with tumorigenesis or outcomes of CRC. To save computational memory, the most differentially expressed genes (in the top 50%) were used to construct a weighted co-expression network. Subsequently, the turquoise module was identified as a clinically significant module. Further enrichment analysis showed that the genes of the turquoise module were remarkably enriched in multiple inflammatory processes and pathways, such as cytokine pathways, chemokine signaling, and immunosuppressive pathways. These results were consistent with previous studies. Chronically localized inflammation not only induces the initiation of CRC, but also promotes the progression of CRC. Among all hub genes in the turquoise module, the expressions of two genes, STAT1 and CCL5, were highly positively correlated. Moreover, the mRNA levels of STAT1 and CCL5 were both significantly elevated in CRC compared with normal colon tissues. This finding indicated that STAT1 and CCL5 overexpression might be valuable biomarkers for predicting the progression of CRC.

STAT1 is a vital participator of innate immunity, protecting the host from infection by pathogens. As a downstream mediator of interferon receptors, STAT1 regulates the expression of multiple anti-pathogens and immunoregulatory genes (16). STAT1 is generally believed to have a bi-directional role in different cancer cells. In some cancers, STAT1 performs its anti-tumor role by enhancing immune surveillance and increasing the tumor-killing effect of natural killer (NK) cells (17). STAT1 has also been shown to suppress the proliferation of cancer cells via the upregulation of cyclin-dependent kinase inhibitor expression or the downregulation of c-myc expression (18,19). In other cancers, such as leukemia, STAT1 promotes leukemogenesis by inhibiting immune clearance (20). Moreover, in some solid tumors, STAT1 increases the expression of several interferon-inducible genes and repairs treatment-mediated DNA damage, and is highly related to chemotherapy resistance (21). Wang et al. reported that STAT1 exhibited pro-survival effect in CRC by increasing PI3K class IB signaling and decreasing the expression of programmed cell death protein 4 (22). Furthermore, when STAT1 signaling was blocked, they also observed that the tumorigenic potency was impaired and there was reduced sensitivity to the inhibition of translation initiation factor eIF4A (22). For some oncogenes, such as CUG2, oncogenic signaling pathway-mediated STAT1 activation significantly enhances the metastatic and drug-resistance properties of colon cancer cells (23). In our study, STAT1 inhibitor markedly suppressed the proliferation of multiple colon cancer cell lines, which suggested STAT1 might contribute to cancer progression.

CCL5, which is also known as RANTES, is a C-C chemokine secreted by immune cells, fibroblasts, and epithelial cells (24). CCL5 was originally found to be an inducer, recruiting leukocytes to inflammatory sites (24). Further investigations showed it to be overexpressed in multiple cancers (25,26). Elevated CCL5 promotes cancer development through regulating cell proliferation, the generation of matrix metalloproteinase, and angiogenesis (27,28). Recent studies demonstrated that increased CCL5 had a substantial influence during the development of CRC. In CRC, CCL5 not only promotes growth and metastasis, but also undermines immune surveillance by increasing the ratio of Tregs and myeloid-derived suppressor cells (29-31). Furthermore, CCL5-β-catenin-Slug has been shown to be a core pathway by which epithelial-mesenchymal transition of CRC is induced (32). In the present study, neutralizing CCL5 significantly decreased the proliferation of CRC, indicating that CCL5 might serve as a potential target for CRC treatment.

Previous studies have shown that STAT1-CCL5 axis is the downstream pathway of cancer-associated molecules such as interleukin 35 and interferon-induced protein 35 (33,34). In pancreatic ductal adenocarcinoma, the activated STAT1-CCL5 pathway promoted micro-vessel formation and monocyte infiltration (34). However, to our knowledge, the role of the STAT1-CCL5 axis in CRC has not been systemically studied. The results of our study show that increased levels of STAT1 and CCL5 might be potential biomarkers for CRC screening. The STAT1-CCL5 axis may be a vital regulator in the development of CRC through its promotion of cell proliferation.

Conclusions

In conclusion, using WGCNA big data mining, we screened out a clinically significant module. The genes of this module were remarkably enriched in inflammation and cytokine-associated pathways. Among all of the hub genes of the clinically significant module, the mRNA levels of STAT1 and CCL5 were significantly higher in CRC than in normal colon tissues. STAT1 expression was highly positively correlated with CCL5. The results of in vitro study showed that interfering with STAT1 and CCL5 could inhibit the proliferation of CRC cells. Our study indicated that the STAT1-CCL5 axis is an important modulator of CRC development through its regulation of cell proliferation. The levels of STAT1 and CCL5 may become valuable biomarkers for CRC screening.

The article’s supplementary files as