The Cancer Genome Atlas (TCGA) database (https://www.cancer.gov/tcga) and the Gene Expression Omnibus (GEO) repository (https://www.ncbi.nlm.nih.gov/geo/) were used to obtain gene expression data. To obtain comprehensive gene expression data, the present study used two major public databases: The Cancer Genome Atlas (TCGA; http://www.cancer.gov/tcga) and the Gene Expression Omnibus (GEO; http://www.ncbi.nlm.nih.gov/geo/). TCGA is a landmark cancer genomics program that has characterized over 20,000 primary cancer and matched normal samples across 33 cancer types. The present study specifically used the TCGA dataset for colorectal adenocarcinoma (COAD), which provides rich genomic and transcriptomic information, enabling the exploration of the molecular underpinnings of colorectal cancer. Additionally, the present study focused on the GSE41258 dataset from GEO, which includes samples from patients with colonic neoplasms treated at Memorial Sloan-Kettering Cancer Center between 1992 and 2004. This dataset comprises 390 expression arrays from primary colon adenocarcinomas, adenomas, metastases, and corresponding normal mucosae. Data processing was performed using R (version 4.4.1, http://www.R-project.org/).

Least absolute shrinkage and selection operator (LASSO) regression was used to create a prognostic signature using the samples from the TCGA cohort. The regression coefficients (β) from the CXCL family model were then combined with the relevant gene expression levels to generate a multi-gene marker-based predictive risk score. The TCGA data was used to train the risk score model, which was built as follows: Risk score=expression level of CXCL1 × (−0.097330812) + expression level of CXCL6 × (0.08298018) + expression level of CXCL8 × (0.009198228) + expression level of CXCL9 × (−0.053426341) + expression level of CXCL11 × (−0.10920036) + expression level of CXCL12 × (0.051584401) + expression level of CXCL14 × (−0.055678695). LASSO regression was implemented with the glmnet package (version 4.1.2; http://CRAN.R-project.org/package=glmnet). Cross-validation and model evaluation were carried out using the ROCR package (version 1.0.11; http://CRAN.R-project.org/package=ROCR). ggplot2 (version 3.3.6; http://CRAN.R-project.org/package=ggplot2) was used for data visualization.

Using the STRING database, a PPI network map was constructed encompassing the CXCL family and their interacting proteins. The network was subsequently refined and visualized within the Cytoscape (version 3.9.1; http://cytoscape.org/download) (31). Hub genes were identified employing the Degree algorithm (32), where the size, color, and spatial arrangement of the nodes within the network corresponded to their respective Degree scores.

The Kaplan-Meier Plotter (http://kmplot.com/) (33) was used to analyze the correlation between gene expression levels and overall survival (OS), Recurrence-Free Survival (RFS) and Post-progression Survival (PPS) in colorectal cancer cohorts.

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed on the related genes and co-expressed genes using the clusterProfiler package (version 4.4.4; http://CRAN.R-project.org/package=clusterProfiler) to identify overrepresented biological processes and pathways.

For the immune infiltration analysis, the single-sample gene set enrichment analysis (ssGSEA) algorithm provided in the R package GSVA (version 1.46.0; http://bioconductor.org/packages/release/bioc/html/GSVA.html) was used to calculate the correlation between CXCL9 and immune cells. The Tumor Immune System Interaction Database (TISIDB) (http://cis.hku.hk/TISIDB/) was employed to analyze the chemokines and chemokine receptors related to CXCL9. TISIDB is a comprehensive database that integrates multi-omics data related to the interaction between tumors and the immune system. It contains information from various sources, including genomics, transcriptomics, proteomics and immunology. This database enables researchers to explore the complex relationships between tumor-associated genes, immune cell infiltration, immune-related pathways, and therapeutic responses (34).

Target RNA sequences were obtained from the National Center for Biotechnology Information (NCBI; http://www.ncbi.nlm.nih.gov). The NCBI is a leading bioinformatics hub under the U.S. National Library of Medicine, offering extensive biological databases, such as GenBank, and powerful analysis tools, such as BLAST, to support global scientific research in genetics, genomics and molecular biology. The prediction of m6A methylation sites was performed using the sequence-based RNA adenosine methylation site predictor (SRAMP; http://www.cuilab.cn/sramp) database. SRAMP is a mammalian m6A sites predictor, which was constructed by extracting and integrating sequence and predicted structural features around m6A sites within a machine learning framework. It shows promising performance in cross-validation tests on its training dataset and in rigorous independent tests (35).

The human colorectal cancer cell line SW480, RKO and the human kidney epithelial cell line 293T were procured from the American Type Culture Collection. Cells were routinely cultured in Dulbecco's Modified Eagle's Medium (DMEM; Gibco; Thermo Fisher Scientific, Inc.) supplemented with 10% fetal bovine serum (FBS; Zhejiang Tianhang Biotechnology Co., Ltd.), 1% penicillin-streptomycin solution. Cells were maintained at 37°C in a humidified incubator with 5% CO₂.

To initiate the present study, three distinct small interfering (si)RNA oligonucleotides were designed specifically targeting the coding sequence of ALKBH5, as well as a scrambled negative control sequence (Table SI). These siRNA sequences were synthesized by Beijing Tsingke Biotech Co., Ltd. Subsequently, each siRNA (at a concentration of 50 nM) was transiently transfected into RKO cells cultured in 6-well plates using Lipofectamine 2000^® (Invitrogen; Thermo Fisher Scientific, Inc.; cat. no. 11668019) at 37°C for 24 h. This process was conducted to evaluate the knockdown efficiency of each siRNA (Fig. S1A). Based on the results, the siRNA sequence with the highest knockdown efficiency (GAAAGGCTGTTGGCATCAATA) was selected to design shRNA primers for further construction of the knockdown vector.

For the overexpression vector, the pCDH plasmid was used as the backbone. For the knockdown vector, the pLKO.1 plasmid was employed. The primers used for amplification and annealing were designed based on the sequences provided in Table SI. These primers, along with the empty vectors, were synthesized or purchased by Beijing Tsingke Biotech Co., Ltd. The detailed procedure involved PCR amplification of the target sequences or annealing of primers at high temperatures, followed by restriction enzyme digestion of the vectors. For the overexpression vector, the pCDH vector was digested using EcoRI (NEB, R3101VVIAL) and BamHI (NEB, R3136VVIAL) restriction enzymes, For the knockdown vector, the pLKO.1 vector was digested using AgeI (NEB, R3552SVIAL) and EcoRI (NEB, R3101VVIAL) restriction enzymes. Homologous recombination (2X MultiF Seamless Assembly Mix, ABclonal, cat. no. RM20523) was performed to integrate the target sequences into the vectors. The ligation products were transformed into competent DH5α cells, which were subsequently cultured overnight at 37°C. Single colonies were selected and further verified by sequencing to ensure the successful construction of the vectors.

To achieve efficient gene knockdown and overexpression, a lentivirus vector system was employed. For lentivirus packaging, 293T cells were co-transfected with a mixture of the recombinant vector or the control vector, the pCMV–VSV-G envelope vector, and the pCMV-Gag-Pol packaging vector (both from Promega Corporation) at a 4:3:1 ratio (totaling 10 µg of DNA) using the PEI reagent (Polyplus-transfection SA). The cells were then incubated at 37°C for 72 h. Following incubation, the viral supernatant was collected and filtered through a 0.22-µm filter to remove cellular debris. The viral titer was determined using the Lenti-Pac^™ HIV qRT-PCR Titration Kit (GeneCopoeia, Inc.). Subsequently, the lentivirus was used to infect RKO cells at a multiplicity of infection of 5. After infection, puromycin was applied to select stable RKO cell lines.

3-Methyladenine (3-MA, MCE) was prepared in DMSO and stored at −80°C until use. Actinomycin D was purchased from MilliporeSigma. Stock solutions were prepared in sterile distilled water at a concentration of 1 mg/ml and stored at −20°C until use. 3-MA was used at a final concentration of 1 mM. Actinomycin D was used at a final concentration of 5 µg/ml.

RNA extraction, cDNA synthesis and qPCR performed according to the manufacturer's protocols. All experiments were repeated at least three times. Total RNA was isolated from cellular samples using the Ultrapure RNA Kit, a product of CWBio. Subsequent to RNA extraction, the Reverse Hifair^® III 1st Strand cDNA Synthesis SuperMix for qPCR, manufactured by Shanghai Yeasen Biotechnology Co., Ltd., was used to perform reverse transcription, converting the isolated total RNA into complementary DNA (cDNA). Following this, the 2X Universal Blue SYBR Green qPCR Master Mix (with UDG), provided by Wuhan Servicebio Technology Co., Ltd., was employed for qPCR, which was conducted on the CFX Connect Real-Time PCR Detection System (Bio-Rad Laboratories, Inc.), under the following thermal cycling conditions: Initial denaturation step at 95°C for 10 min, followed by 40 cycles consisting of denaturation at 95°C for 10 sec, annealing at 60°C for 20 sec, and extension at 72°C for 15 sec. Relative gene expression was quantified using the comparative Cq (2^−ΔΔCq) method (36), with GAPDH serving as an endogenous reference gene for normalization of the data. The primers used for qPCR were shown in Table SII.

The Cell Counting Kit-8 (CCK-8) assay was employed to quantitatively assess cell proliferation. Cells were seeded at a density of 1×10⁴ cells per well in 96-well plates and incubate at 37°C with 5% CO₂ for 12 h to allow for cell adhesion. After cells attached to 96-well plate for 0–120 h, 10 µl of the Cell Counting Kit-8 (Dalian Meilun Biology Technology Co., Ltd.) was added to each well and incubated for 1 h at 37°C. Subsequently, the optical density of the samples was measured using a spectrophotometer, with the wavelength set to 450 nm which is indicative of viable cell count.

The upper aspect of the Transwell membrane was uniformly layered with Matrigel, followed by a 30-min incubation period at 37°C to facilitate solidification. Subsequently, a cellular suspension consisting of 1×10⁵ cells in 100 µl of serum-free medium was dispensed into the upper chamber of each Transwell insert. Concurrently, the lower chamber was supplemented with 600 µl of medium enriched with 10% serum to act as a chemotactic agent. After being incubated for 20 h at 37°C, the cells remaining at the upper surface of the membrane were removed by a cotton swab. The cells that had successfully traversed to the inferior surface of the membrane were then immobilized using a 4% paraformaldehyde solution for a 15-min fixation at room temperature and stained with 1X modified Giemsa stain (Beyotime Institute of Biotechnology) for 45 min at room temperature. The migratory cells on the underside of the membrane were subsequently captured using a fluorescence microscope for documentation.

Once a confluent monolayer was established, the cell monolayers were subjected to mechanical wounding via a sterile 200 µl micropipette tip to create a standardized wound. To eliminate cellular debris and ensure accurate wound assessment, the wounded monolayers were copiously rinsed with PBS multiple times. Subsequently, the cells were maintained in a serum-free medium for a duration of 72 h to simulate wound healing conditions. The progression of in vitro wound closure was monitored using a fluorescent microscope at designated time intervals. The in vitro wound healing was recorded and assessed using the closure rate, which was calculated as follows: [(Wound area at 0 h-Wound area at × h)/Wound area at 0 h] ×100%.

Cells were subjected to lysis in a RIPA buffer, which was complemented with Phenylmethanesulfonyl Fluoride (PMSF, Beyotime Institute of Biotechnology). The resultant lysate was centrifuged at 12,000 × g for a duration of 15 min at 4°C to achieve clarification. Subsequently, the protein concentration within the supernatant was ascertained using a BCA (Beyotime Institute of Biotechnology; cat. no. P0012) assay. The Omni-Easy^™ One-step Color PAGE Gel Rapid Preparation Kit (Epizym, Inc.; PG211 and PG213) was used to prepare 7.5 and 12.5% gels. An aliquot of the protein lysate, containing 20 µg of total protein, was resolved on SDS-PAGE and subsequently transferred onto PVDF membrane (MilliporeSigma; cat. no. IPVH00010). The membranes were pre-incubated with a blocking solution consisting of 5% non-fat milk for 1 h at room temperature. Thereafter, the membranes were incubated with the respective primary antibodies at 4°C overnight. Following this incubation, the membranes were washed thoroughly and then exposed to anti-rabbit/mouse/rat secondary antibodies for 1 h at room temperature. The immunoreactive bands were visualized using an ECL substrate (Biosharp Life Sciences; cat. no. BL520A), and the membrane was detected using a chemiluminescence imaging system (eBlot). Primary antibodies used in this study were: Mouse anti-GAPDH (1:5,000; Proteintech Group, Inc.; cat. no. 60004-1-Ig), rabbit anti-LC3B (1:1,000; ABclonal Biotech Co., Ltd.; cat. no. A11282), rabbit anti-P62 (1:2,000; ABclonal Biotech Co., Ltd.; cat. no. A19700), rabbit anti-CXCL9/MIG (1:1,000; ABclonal Biotech Co., Ltd.; cat. no. A25183); and the secondary antibodies were: HRP Goat Anti-Rat IgG (H+L; 1:10,000; ABclonal Biotech Co., Ltd.; cat. no. AS028), HRP Goat Anti-Mouse IgG (H+L; 1:10,000; ABclonal Biotech Co., Ltd.; cat. no. AS003) and HRP Goat Anti-Rabbit IgG (H+L; 1:10,000; ABclonal Biotech Co., Ltd.; cat. no. AS014). The gray value was measured by ImageJ software (version 1.54f; National Institutes of Health).

Total RNA was extracted from cells, for both the IgG group and the m6A group, 100 µg of total RNA was added, along with an appropriate quantity of RNase inhibitor. Subsequently, 1 µl of either IgG antibody (Proteintech Group, Inc.; cat. no. B900620; for the IgG group) or m6A antibody (Proteintech Group, Inc.; cat. no. 68055-1-Ig; for the m6A group) was introduced. IP buffer was then supplemented to adjust the volume to 200 µl. After thorough mixing, the samples in both groups were incubated overnight at 4°C. Protein A/G Magnetic beads are washed with IP buffer, and the premix is added and incubated at 4°C for 3 h. The beads with bound RNA are washed 4–5 times with IP buffer. After adding TRIzol^® (Thermo Fisher Scientific, Inc.) an RNA purification kit (CWBio.) to was used to extract the RNA on the magnetic beads. The purified m6A-enriched RNA were analyzed by RT-qPCR as aforementioned to quantify specific m6A/modified RNAs.

All data analyses were conducted using GraphPad Prism software, version 8.0.1 (Dotmatics). The experimental groups were compared using the non-parametric two-tailed Student's t-test to assess the statistical significance between two independent samples. For multiple group comparisons, a one-way analysis of variance was employed, followed by Tukey's post hoc test to correct for multiple comparisons and to identify the source of significant variance within the dataset. Each treatment condition for the cellular samples was replicated a minimum of three times to ensure the reproducibility and reliability of the experimental outcomes. Data are presented as the mean ± standard error of the mean to reflect the variability within the sample set. P-values are reported for all statistical tests. P<0.05 was considered to indicate a statistically significant difference. Specific P-values are indicated in the figure legends and results sections (*P<0.05, **P<0.01, ***P<0.001).

Variance analysis indicated that the CXCL family displayed significant differential expression in COAD, as presented in Fig. 1A. A prognostic model of CXCL family in COAD was constructed through LASSO regression, which has good diagnostic value (Fig. 1B-E). COX regression analysis was performed on seven genes in the prognostic model of CXCL family (Fig. 1F), and the association between molecular and methylation genes with P<0.2 in univariate regression was analyzed. The results showed that CXCL8, CXCL9, and CXCL11 have an association with methylation genes (Fig. 1G). The PPI network with CXCL family and their interacting genes shows that CXCL9 is the hub gene with the highest degree score in the prognostic model genes of the CXCL family (Fig. 1H).

Considering the significant role of CXCL9 in the prognostic model of the CXCL family, the expression of CXCL9 in cancer was further analyzed, particularly in COAD, as well as its effect on the prognosis of COAD patients. Results showed that CXCL9 is highly expressed in most types of cancer (Fig. 2A) and is markedly associated with patient PFI in cancers such as breast cancer (BRCA), COAD, KIRP, LGG, SARC and UCEC (Fig. 2B). The TCGA and GSE41258 dataset from GEO indicated that CXCL9 was highly expressed in COAD (Fig. 2C and D). An analysis of the relationship between CXCL9 expression and various clinical indicators in the TCGA dataset revealed that CXCL9 expression exhibited a trend of first increasing and then decreasing with the progression of COAD (Fig. 2E). Slight differences in CXCL9 expression were observed among different ethnicities and in different residual tumor classifications (Fig. 2F and G). KMplot showed that patients with elevated CXCL9 expression exhibited markedly improved OS, RFS and PPS compared to those with lower expression levels (Fig. 2H-J).

To further clarify the molecular mechanisms underlying the role of CXCL9 in COAD, the present study conducted enrichment analyses using the GO, KEGG and GSEA. GO and KEGG enrichment analysis was performed using the top 200 most-associated genes. The EMAP plot revealed significant enrichment of CXCL9 and its associated genes in pathways related to T cell activation and leukocyte adhesion (Fig. 3A). In GSEA enrichment analysis using co-expressed genes, numerous pathways related to the functions of T cells and B cells were identified, such as Co-stimulation by the CD28 family and CD22 mediated BCR regulation (Fig. 3B). In addition, pathways related to Fc epsilon Receptor I (FcεRI) and mitochondria related pathways as well as DNA methylation, and the Wnt Beta catenin signaling pathway, cell surface interactions at the vascular wall and others were also enriched (Fig. 3C-H). GO and KEGG combined with FC bubble chart demonstrated a correlation with the cell cycle related pathway (Fig. 3I).

Immune infiltration analysis indicated that CXCL9 in CRC was positively associated with the majority of immune infiltration cells, predominantly cytotoxic cells, activated dendritic cells (aDC), Th1 cells, and T cells (Fig. 4A). The effect of CXCL9 on CRC immune infiltration varies across different stages of cancer progression. In Stage I, CXCL9 shows a positive correlation with cytotoxic cells, aDC, and Th1 cells (Fig. 4B). In Stage III, it maintains a positive correlation with cytotoxic cells, macrophages, and T cells (Fig. 4D). In Stage IV, the positive correlation was observed with cytotoxic cells, Th1 cells and CD8 T cells (Fig. 4E). However, a notable alteration in immune infiltration is observed in Stage II, where CXCL9 exhibited a negative correlation with numerous immune cells, including cytotoxic cells, aDC, macrophages, Th1 cells, T cells, and neutrophils (Fig. 4C). The TISIDB analysis suggested that CXCL9 in COAD was most closely associated with chemokines CCL4, CCL5, CCL8, CCL18, CXCL10, CXCL11 and CXCL13, as well as chemokine receptors CCR1, CCR2, CCR5, CCR8 and CXCR6 (Fig. S2A and B). Therefore, the expression of these chemokine receptors in different stages of COAD was analyzed and it was found that the expression patterns of CCR5, CCR8 and CXCR6 were similar to CXCL9 (Fig. S2C).

Subsequent to bioinformatics analysis of CXCL9, in vitro assays were performed to explore its potential role in promoting tumor growth in colorectal cancer cells. In colon cancer cell lines, the expression of the CXCL9 gene was examined and it was found that its expression level was markedly higher in RKO and SW480 cells compared to NCM460 (Fig. 5A). Through qPCR and western blotting validation, cell lines with overexpression of CXCL9 were successfully constructed (Fig. 5B-D).

Wound healing assays were performed to assess the effect of CXCL9 overexpression on cell migratory capacity. The results demonstrated that cells overexpressing CXCL9 closed the wound area more rapidly than control cells within 48 h (Fig. 5E and F). The invasive potential of colorectal cancer cells with CXCL9 overexpression was further evaluated through Transwell assays. Notably, the number of cells that traversed the chamber was markedly higher when compared with the control cohort (Fig. 5G and H), indicating a significant increase in invasiveness (P<0.01) attributed to the overexpression of CXCL9. CCK-8 assays were utilized to measure the impact of CXCL9 overexpression on cell proliferation. The absorbance at OD450, which reflects cell viability, showed that cells with overexpressed CXCL9 exhibited a higher proliferation rate during 0–120 h (Fig. 5I). These results indicated a potential role of CXCL9 in the progression of colorectal cancer.

The present study treated SW480 and RKO cells with 3-MA and found changes in the expression of the CXCL9 gene within 0–8 h (Fig. 6A and B). These observations suggested a potential role of CXCL9 in autophagy regulation Following CXCL9 overexpression, an increase in the RNA levels of autophagy-related genes was observed, indicating a potential involvement of CXCL9 in the regulation of autophagic flux (Fig. 6C). Western blotting results indicated that overexpression of CXCL9 can inhibit the degradation of p62 and suppress the conversion of LC3-I to LC3-II, thereby blocking autophagy flux (Fig. 6D).

SRAMP identified a high-confidence m6A methylation site within the CXCL9 gene (Fig. 7A and B). This implied that CXCL9 may be subject to post-transcriptional regulation through m6A modification, which could have significant implications for its expression and function in cellular processes. Through bioinformatics analysis, a strong correlation between CXCL9 and the m6A eraser ALKBH5 was predicted (Fig. 7C). Consequently, the expression levels of CXCL9 in cell lines with ALKBH5 knockdown and overexpression were examined. The results showed that the expression of CXCL9 increased in response to ALKBH5 overexpression and decreased following ALKBH5 knockdown (Fig. 7D and E). Actinomycin D chase experiments revealed that the stability of CXCL9 was markedly reduced in ALKBH5 knockdown cells compared to the control group, while overexpression of ALKBH5 markedly enhanced the stability of CXCL9 (Fig. 7F and G). The MeRIP experiment demonstrated that the m6A level of CXCL9 decreased upon overexpression of ALKBH5, while it increased upon the knockdown of ALKBH5 (Fig. 7H and I). The present study stably expressed a knockdown of ALKBH5 in RKO cells and transiently transfected them with the OE-CXCL9 plasmid. The results from wound healing, CCK8, and Transwell assays demonstrated that CXCL9 can partially restore the effects of ALKBH5 knockdown on colorectal cancer cells (Fig. 8).