RNA-Seq transcriptomic data of colon adenocarcinoma patients were retrieved from TCGA-COAD project via the R TCGAbiolinks package (version: 1.12.0) (Colaprico et al., 2015). The corresponding clinical information of the TCGA-COAD samples was obtained from the UCSC Xena database (Goldman et al., 2020). Microarray datasets (GSE39582, GSE81653, GSE87211, GSE28702 and GSE17536) for colorectal cancer samples and matching clinical characteristics were downloaded from the NCBI Gene Expression Omnibus (GEO) database. Microarray datasets for 5-Fluorouracil and oxaliplatin resistant cellular models (GSE76489 and GSE83131) were also obtained from GEO. Gene expression profiles of colorectal cancer cell lines and corresponding half-maximal inhibitory concentration (IC50) of 5-Fluorouracil and oxaliplatin were obtained from the Genomics of Drug Sensitivity in Cancer (GDSC) database (Yang et al., 2013).

Colorectal cancer (CRC) samples (four with lymph node metastasis and three without lymph node metastasis) were obtained from patients who underwent curative surgery at the First Affiliated Hospital of Soochow University (Suzhou, Jiangsu, P.R. China). All patients included in this study provided written informed consent for participation. The study was approved by the Biomedical Research Ethics Committee of Soochow University (Suzhou 215123, Jiangsu Province, China).

Total RNA of each sample was extracted using TRIzol reagent. The Agilent Human lncRNA-mRNA Microarray V2.0 4 × 180 K (Agilent Technologies, Inc.) was used to compare transcriptomic profiles between subgroups with and without lymph node metastasis. RNA labeling, microarray hybridization, and data acquisition were performed by Shanghai OE Biotech. Co., Ltd. The raw microarray data reported in this paper has been deposited in the OMIX, China National Center for Bioinformation/Beijing Institute of Genomics, Chinese Academy of Sciences, under the accession number OMIX505 and can be accessed at https://ngdc.cncb.ac.cn/omix.

For TCGA-COAD RNA-Seq data, the raw reads count matrix was converted to the counts per million mapped reads (CPM) format. Trimmed mean of M values (TMM) normalization was carried out to correct the variation of gene expression abundance across different samples. Both log2-transformed or raw formats were generated for further bioinformatics analyses. Samples were categorized according to the median expression level of target genes. Differentially expressed genes were filtered using the Benjamini-Hochberg false discovery rate approach (adjusted p-value < 0.05, and absolute value of log2 fold-change > 1). All preprocessing and differential gene expression analysis was conducted using an R pipeline based on “limma” (version: 3.48.1) (Ritchie et al., 2015) and “edgeR” (version: 3.34.0) (Robinson et al., 2010) R packages.

For Affymetrix microarray data, data was normalized using the Robust Multi-array Average (RMA) method. For Agilent microarray data, data were normalized using the “quantile” method. Gene level expression estimates were log2-transformed and summarized into a single matrix for subsequent analysis. Raw p-values < 0.01 were used as the cutoff to identify differentially expressed genes. All microarray analyses were performed using the “limma” R package (version: 3.48.1).

Genes were ranked according to the median absolute division (MAD) values across samples in the GSE81653 and GSE28702 dataset, and the top 5,000 genes were subjected to the subsequent weighted gene co-expression network analysis (WGCNA). Briefly, a minimum power with R ² > 0.85 was set to optimize the power (β) for automatic network and further topological overlap matrix construction. After that, minModuleSize = 30 was utilized to initially classify co-expression modules. In one data set with less than 100 samples (GSE28702), the cutHeight was set to 0.2. Conversely, when the sample size exceeded 100 (GSE81653), this parameter was set to 0. A p-value less than 0.05 indicated statistical significance for the association between co-expression modules and clinical traits. Module membership (MM) values >0.7 served as the cutoff for key gene screening. The intersection between the results of GSE81653 and GSE28702 was taken for further analysis. All computational analyses were performed using the “WGCNA” R package (version: 1.70-3) (Langfelder and Horvath, 2008).

The STRING database (version 11.0) (Szklarczyk et al., 2019) was used to identify a protein-protein interaction network linked by the target gene connecting the significant gene modules yielded by WGCNA and the up-regulated genes in the colorectal cancer samples with lymph node metastasis (filtered by p-value < 0.05 in the “limma” workflow). The “clusterProfiler” R package (version: 4.0.2) (Wu et al., 2021) was used to annotate the biological function of the genes in the network. A raw p-value less than 0.05 indicated statistically significant enrichment.

Gene set enrichment analysis (GSEA) was performed to determine whether certain biological pathways or a priori defined gene sets showed statistically significant differences between two subgroups with high and low expression of the target gene. The “clusterProfiler” R package (version: 4.0.2) was used for enrichment analyses (Wu et al., 2021). Gene sets including Chemical and Genetic Perturbations, KEGG pathway, and Gene Ontology Biological Process, were accessed using the Molecular Signatures Database (MSigDB). Gene markers for lymph node metastasis were defined as the top 300 significantly up-regulated genes (p-value less than 0.01) ranked by log2-transformed fold change obtained from our own microarray data. A raw p-value less than 0.05 indicated statistically significant enrichment.

GSEA was also used to determine whether the gene signatures that changed upon the expression of the target gene showed statistically significant, concordant differences in 5-Fluorouracil or Oxaliplatin resistant and naïve cell lines. Significant differentially expressed genes between the TCGA-COAD subgroups with high and low expression of the target gene were ranked according to the log2-transformed fold change. The top 300 up- and down-regulated genes were defined as the gene sets. GSEA was then applied to examine the enrichment of these gene sets in the GSE76489 (5-Fluorouracil resistance model) and GSE83131 (Oxaliplatin resistance model) datasets, respectively. A raw p-value less than 0.05 indicated statistically significant enrichment.

The Pearson correlation between log2-transformed levels of the target gene expression and ln-transformed IC50 from the GDSC database were calculated. A p-value less than 0.05 indicated statistical significance.

We used the GSE39582 dataset for evaluation of recurrence risk. Samples were categorized according to the median expression level of the target gene. Odds ratios (ORs) and corresponding 95% confidence intervals (95% CIs) were calculated to assess the risk of colorectal cancer recurrence. Stratified analyses according to TNM stage, size of primary tumor (T), regional lymph node involvement (N), presence of distant metastatic spread (M), and mutation status of TP53, KRAS, and BRAF were performed to assess whether the effect of the target gene is influenced by the above-mentioned confounding factors.

TCGA-COAD, GSE17536, GSE87211 and GSE39582 datasets were used for prognostic analyses to assess the impact of the target gene on colorectal cancer survival. The optimal cutoff values for each dataset were determined by the algorithms embedded in the “survminer” (version: 0.4.9) R package. The logrank test was performed to compare the survival distributions between the assigned subgroups. A p-value less than 0.05 indicated statistical significance.

CPM reads of the TCGA-COAD dataset were assessed for immune cell infiltration analysis. Samples were categorized according to the median expression level of the target gene. The cell-type identification by estimating relative subsets of RNA transcripts (CIBERSORT) method (Newman et al., 2015) was used to profile tumor infiltrating immune cells and compare the differences between the two subgroups. A p-value less than 0.05 indicated statistical significance.

All bioinformatics analyses were carried out in the R programming environment (version: 4.1.0).

To identify key genes associated with specific clinical traits of colorectal cancer, we employed the WGCNA method to identify gene module-trait correlations. We categorized the top 5,000 genes with the highest MAD in the GSE28702 dataset into 22 co-expression modules (Figure 1A). The MEgrey60 and Melightgreen modules, containing 195 genes, were positively associated with FOLFOX chemotherapy resistance (Figure 1B) and we selected 44 key genes with a MM value greater than 0.7. Next, we performed WGCNA analysis for GSE81653 using the same methods and parameters as previously to explore the relevant factors of recurrence risk after FOLFOX therapy which resulted in clustering of 7 co-expression modules (Figure 1C). The Mebrown and Meblue modules consisting of a total of 216 genes were positively associated with recurrence risk (Figure 1D) and we again selected 114 hub genes with a MM value higher than 0.7. We then intersected the key gene members from both WGCNA analyses corresponding to FOLFOX resistance and tumor recurrence which identified the complement component 3 (C3) gene as the only common element (Figure 1E).

We then constructed a protein-protein interaction network by integrating the interaction relationships of the significant modules from the WGCNA analyses above and the differentially expressed genes between samples with and without lymph node metastasis to delineate the core genes connecting subnetworks responsible for FOLFOX chemotherapy resistance, tumor recurrence, and spread to lymph nodes. As highlighted in the interaction network (Figure 2A), the C3 gene represents one of the bridge hubs linking the three subnetworks. Further functional annotation analysis identified a significant enrichment of gene sets related to cancer stemness, invasiveness, migration, proliferation, angiogenesis, and inflammation response. In addition, cellular responses to chemical stimulus, a pathway associated with chemotherapy resistance, was also significantly enriched (Figure 2B). Based on these results, we propose that the C3 gene may act as a potential regulator of the prognostic CRC network and serves as a shared hub gene, rather than as a node with scattered functional edges participating in a single signaling pathway.

We next compared the effect of high and low expression of C3 on tumor recurrence in patients with CRC in the GSE39582 dataset. Patients with high C3 expression exhibited an increased risk of recurrence (OR = 1.52, 95% CI: 1.05–2.22; p = 0.023) (Figure 3). Further subgroup analysis showed that the levels of C3 expression had significantly positive effect on colorectal cancer recurrence in patients during early TNM stages (p = 0.002) and those without spread to lymph nodes (p = 0.009). In addition, we found a statistically significant susceptibility in subgroups without mutations in TP53 (p = 0.010), KRAS (p = 0.004), and BRAF (p = 0.019).

Subsequent survival analysis was performed to validate the prognostic value of C3 in four independent transcriptome datasets with a total of 1,358 samples. Kaplan-Meier plots revealed that high C3 expression was associated with a statistically significant unfavorable progression−free survival (Figure 4A, p = 0.023), disease−free survival (Figures 4B,C; both p-values < 0.05), and recurrence−free survival (Figure 4D, p = 0.022). Therefore, we propose that C3 represents a candidate biomarker for poor prognosis in colorectal cancer patients.

To explore the impact of C3 gene expression on resistance to FOLFOX chemotherapy, we assessed the correlation between C3 gene expression in different colorectal cancer cell lines and the corresponding IC50 of 5-Fluorouracil and Oxaliplatin. C3 gene expression was positively associated with Oxaliplatin IC50 (Figure 5 A, p = 0.019), but not with 5-Fluorouracil (Figure 5B, p = 0.120). Next, we compared the transcriptome concordance between chemotherapy reagent resistance and C3 overexpression. The top 300 up-regulated genes were obtained by comparing the samples with a high and low abundance of C3 expression (no significantly down-regulated gene were identified because of insufficient log2-transformed fold changes). The GSEA results found these genes were significantly enriched in the comparison between Oxaliplatin-resistant cells and naïve cells (Figure 5C, p = 0.003). However, the enrichment was not statistically significant in the cell model of 5-Fluorouracil resistance (Figure 5D, p = 0.832). This was consistent with the results of in silico IC50 tests. Based on the observations above, it could be reasonably speculated that a high C3 gene expression in colorectal cancer might promote resistance to Oxaliplatin.

To explore the impact of C3 gene expression on the tumor immune microenvironment, we employed the CIBERSORT algorithm to compare the abundance of 22 human infiltrating immune cell types between colorectal cancers with C3 gene expression above and below median level (Figure 6). This revealed a large heterogeneity of infiltrating immune cells between the two different C3 gene expression groups. There were significantly more infiltrating dendritic cells activated, mast cells activated, monocytes, and T cells CD4 memory resting in colorectal cancers with C3 lower expression compared to those with higher C3 expression (all p-values < 0.001). Conversely, tumors with high C3 expression showed significantly higher fractions of M0 macrophages, M1 macrophages, and mast cells resting (all p-values < 0.01). No significant difference of other immune cells fraction was found (all p-values > 0.05).

We then performed GSEA to elucidate the biological pathways regulated by C3 that are involved in the progression of CRC. This revealed that clinical CRC samples with higher C3 expression were enriched for the general pathway of colorectal cancer development and progression (Figure 7A). Furthermore, lymphatic metastasis signatures identified by our microarray data were also significantly overrepresented (Figure 7B) and likewise, genes associated with invasion of lymphatic vessels during metastasis also showed marked enrichment (Figure 7C). In addition, further gene set mining suggested that higher C3 expression could promote cancer proliferation (Figure 7D), invasion (Figure 7E), and adhesion (Figure 7F). Taken together, these observations corroborated the hypothesis that a higher C3 gene expression could drive poor prognosis of colorectal cancer. To further uncover the underlying mechanisms, we next investigated alterations of biological pathways upon C3 hyperexpression. Colorectal cancers with high C3 gene level exhibited a hyperactivation of several pathways including pathways in cancer (Figure 7G), TGF-BETA signaling (Figure 7H), JAK-STAT signaling (Figure 7I), MAPK signaling (Figure 7J), VEGF signaling (Figure 7K), and WNT signaling (Figure 7L).

To summarize, overexpression C3 gene is likely to activate several signaling pathways and plays a prominent role in initiating cancerous deterioration, driving poor prognosis.