A flowchart of this study is presented in Figure 1 . Gene expression profiles of Datasets (GSE3964, GSE19860, GSE104645, GSE106584, GSE69657) were downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi). A total of 230 CRC samples were retrieved for subsequent analysis from the five datasets. The first four datasets contained expression profiling of 200 clinical samples collected from CRC patients before the exposure to 5-FU-based chemotherapy and the last one contained expression profiling of 30 clinical samples collected from CRC patients after the exposure to 5-FU-based chemotherapy. Analysis of gene expression profiles between chemotherapy-resistant patients and chemotherapy-sensitive patients in datasets GSE3964, GSE19860, GSE104645, and GSE106584 were used to identify biomarkers associated with innate drug responses. Gene expression profiles of dataset GSE69657 were used to explore the tumor microenvironment expression patterns of the patients who received 5-FU-based chemotherapy. Another gene expression profile (FPKM) of 584 CRC patients (426 COAD, 158 READ) was downloaded from UCSC (https://gdc.xenahubs.net/download/TCGA-COAD.htseq_fpkm.tsv.gz; https://gdc.xenahubs.net/download/TCGA-READ.htseq_fpkm.tsv.gz). Samples with a follow-up time of fewer than 30 days were excluded.

Included in our study were 46 CRC patients of GSE3964, 29 CRC patients of GSE19860, 54 CRC patients of GSE104645, and 71 CRC patients of GSE106584. Differentially expressed genes (DEGs) in datasets GSE3964, GSE19860, GSE104645, and GSE106584 were analyzed by R packages, respectively. Genes with p < 0.05 was considered a significantly differentially expressed gene between the chemotherapy-resistant and chemotherapy-sensitive patients.

The DEGs of datasets GSE3964, GSE19860, GSE104645, and GSE106584 were separately constructed to gene co-expression networks using the WGCNA package in R to seek the modules of highly associated genes among samples for relating modules to external sample traits (19). A weighted adjacency was constructed by calculating Pearson correlations of all gene pairs and a soft power value was selected to construct a standard scale-free network in each dataset. The similarity matrix completed by Pearson correlation of all gene pairs was transformed into topological overlap matrix (TOM) and corresponding dissimilarity (1-TOM). Similar gene expression was classified into different gene co-expression modules using a hierarchical clustering dendrogram of the 1-TOM matrix. Module-resistant associations were calculated to locate functional modules in the co-expression network. Modules with high correlation coefficients (the absolute value of correlation coefficient ≥0.4, P<0.05) were considered to be related to resistance and were extracted for further analysis.

The immune-related genes (IRGs) list was retrieved from the ImmPort database (https://www.immport.org/shared/genelists). Overlapping immune-related genes, from the results of WGCNA and the IRGs, were selected for further analysis.

Both related genes from the GSE106584 were verified on survival data from itself. The genes screened from the other GEO datasets were validated on survival data from the TCGA cohort. Univariate Cox regression analysis was used to filter all candidate genes and multivariate Cox regression analysis was used to filter the genes from univariate Cox regression analysis. The meaning genes from multivariate Cox regression analysis were selected and compared to the expression of these genes of 5-FU resistant and 5-FU sensitive patients in each dataset.

The ESTIMATE algorithm was used to determine the scores of CRC patients in GSE69657 to explore the role of the tumor microenvironment in 5-FU resistance. We compared the differentially expressed genes of the high-score and low-score groups to acquire 5-FU resistance-related genes in the tumor microenvironment.

The R package of pRRophetic was used to predict IC50 of 5-FU in GSE19860. IC50 implies the efficiency of a substance in restraining certain biological or biochemical functions. “cgp2016” was selected to be the predicted library (20). The quartile of RBP7 expression as the cut-off, patients were divided into RBP7-high expressed and RBP7-low expressed groups.

Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses for genes in the tumor microenvironment differential genes of 5-FU resistance were performed using an R package “clusterProfiler”. GO annotation was based on three categories including, biological processes (BP), cellular compartments (CC), and molecular functions (MF). Terms in GO and KEGG with adj.P.value < 0.05 were considered significantly enriched and were visualized by R package “enrichplot” and “ggplot2”.

Excluding incomplete information samples, there were 401 of 584 CRC patients (286 COAD, 115 READ) to be explored the RBP7 clinical meaning in the TCGA cohort. The median of RBP7 expression as the cut-off, patients were divided into RBP7-high expressed and RBP7-low expressed groups. Chi-square or Fisher’s exact tests were used to compare the differences in clinical parameters between RBP7-high expressed and RBP7-low expressed groups.

To identify signaling pathways that are differentially activated in 5-FU resistance CRC with different expression RBP7, we selected an ordered list of genes through the limma R package and conducted Gene Set Enrichment Analysis (GSEA) with adjusted p < 0.05 using the clusterfiler R package in the GSE19860 dataset, in which BBP7 expression was statistically different between the 5-Fu sensitive and the resistant groups.

To identify downstream signaling pathways that are differentially activated in 5-FU resistance CRC with different expression RBP7, we performed GSVA with adjusted p < 0.05 to explore correlated pathways of RBP7 in the GSE19860. Hallmark gene sets “h.all.v7.5.1.symbols.gmt” were downloaded from Molecular Signatures Database (https://software.broadinstitute.org/gsea/downloads.jsp).

The cell plots of single-cell sequencing of colorectal tumors and adjacent non-tumor colon tissue were downloaded from ArrayExpress databases (https://www.ebi.ac.uk/gxa/sc/home). We searched for colorectal cancer, obtained the result of single-cell sequencing for colorectal tumors and adjacent non-malignant colon tissue, drew pictures online, and downloaded. Drawing path: plot type: UMAP, plot options: n_neibors:100, color plot by: ontology labels, gene name: RBP7.

Lovo and Lovo/5-FU cells (induced by increasing continuous exposure concentration of parental cells Lovo) were cultured in DMEM/F-12 medium (BasalMedia) containing 10% fetal calf serum (BI), and 100 U/ml each of penicillin and strepcomycin (BI) at 37°C with 5% CO2. The maintenance concentration of Lovo/5-FU cells was 0.01mg/mL.

To establish 5-FU resistant cells, Lovo cells were induced by increasing continuous exposure concentration. Lovo cells cells were treated in a culture medium containing 0.1μg/mL 5-FU. When the cell survival rate was greater than 90% and the cell growth could be maintained, the dose was increased by 1.5 times, and repeated until the cells could grow stably in the culture medium with 5-FU concentration of 0.01mg/ml. The 5-FU resistant Lovo cells were obtained by continuous exposure to gradually increased concentrations of 5-FU for eight months.

Proteins were separated by 10% SDS-PAGE and then transferred onto a nitrocellulose membrane. The following primary antibodies were used: anti-GAPDH (Cell Signaling Technology) and anti-RBP7 (ABclonal Technology). The following secondary antibody was used: goat anti-rabbit IgG-HRP antibody. The proteins were visualized using an ECL detection kit (FDBIO).

Full-length RBP7 DNA was amplified by PCR with the primer 5’-CTTTGCCACTCGTAAAATAGCCA-3’ and 5’-CGTGTGGATGGTAAAAGAATCCC-3’ (Tsingke Biotechnology, China). Full-length GAPDH DNA was amplified by PCR with the primer 5’- GGAGCGAGATCCCTCCAAAAT-3’ and 5’-GGCTGTTGTCATACTTCTCATGG-3’ (Tsingke Biotechnology).

All statistical analyses were performed using Prism (version 8) and R (version 4.1.2). The proportional composition of the two variables was compared using the Chi-square or Fisher’s exact tests. The comparison of RBP7 expression and the ESTIMATE scores between the two datasets were performed using the Mann-Whitney U tests. Predicted 5-FU sensitivity difference between groups was tested by unpaired t test.The median value was set as the cut-off. Survival analysis, WGCNA, ESTIMATE algorithm, and the COX regression analysis were carried out by R version 4.1.2 and corresponding packages. Statistical significance was set at P.value <0.05 or adj.P.value < 0.05.

Firstly, datasets GSE3964, GSE19860, GSE104645, and GSE106584 all contained chemotherapy-sensitive and chemotherapy-resistant samples. To explore the differentially expressed genes related to chemotherapy response, we used the limma package for differential analysis. There were 1,068 DEGs in the GSE19860 mRNA profile, 651 DEGs in the GSE3964 mRNA profile, 1,392 DEGs in the GSE 104645 mRNA profile, and 1,237 DEGs in the GSE 106584 mRNA profile ( Figure 2 ).

To further explore the functional clusters related to 5-FU based chemotherapy resistance, the weighted gene co-expression network was constructed from the DEGs sets which filter from GSE3964, GSE19860, GSE104645, and GSE106584. The soft-thresholding power in WGCNA was determined based on a scale-free R2 (R2 = 0.9). Hub modules were identified based on the average linkage hierarchical clustering and the soft-thresholding power. To evaluate the link between modules and clinical traits (5-FU resistant and 5-FU sensitive), the heatmap of the module-trait relationship was plotted ( Figure 3 ). All modules showed a high correlation with 5-FU resistance of CRC in GSE3964 and contained 651 5-FU resistance-related genes. All modules showed a high correlation with 5-FU resistance of CRC in GSE19680 and contained 1,068 5-FU resistance-related genes. Only one module showed a high correlation with 5-FU resistance in GSE104645 and contained 381 5-FU resistance-related genes. All modules showed a high correlation with 5-FU resistance of CRC in GSE106584 and contained 1,237 5-FU resistance-related genes.

After differentially expressed genes analysis and further screening by WGCNA, we identified gene clusters associated with 5-FU resistance. We obtained 2,483 immune-related genes from the Immport database. The genes related to 5-FU resistance, screened out from GEO datasets, were intersected with immune-related genes to obtain immune-related genes to 5-FU resistance ( Figure 2 ). There were 23 IRGs to 5-FU resistance for CRC patients in GSE3964. There were 75 IRGs to 5-FU resistance for CRC patients in GSE19860. There were 28 IRGs to 5-FU resistance for CRC patients in GSE104645. There were 39 IRGs to 5-FU resistance for CRC patients in GSE106584. HLA-DQA1 and MX2 were both in GSE104645 and GSE19860. TUBB3 was also found both in GSE19860 and GSE3964. One research suggested that patients with high TUBB3 had a statistically significant poorer OS when undergoing docetaxel-based versus 5-FU/LV chemotherapeutic (21).

The IRGs from GSE106584 were verified on its survival data and the IRGs from the other three datasets were verified on survival data in the TCGA cohort. In GSE106584, univariate Cox regression analysis suggested that 15 genes were related to DFS, and 9 genes mean to OS as well. Meanwhile, in the TCGA cohort, 19 genes were mean to OS and 9 genes were correlated to DFS, but none of the genes correlated to OS in the GSE106584 resulted from the multivariate Cox regression. In The TCGA cohort, 4 genes were associated with OS by the multivariate Cox regression analysis. Finally, we found 13 IRGs linked to 5-FU resistance.

Multivariate Cox regression analysis was used to develop immune-related prognostic gene signatures. After obtaining the coefficients of each gene, we calculated the risk score of each CRC patient with the computational equation. In GSE106584, the DFS-prognostic model included nine genes: HSPA8, RARB, RABEP2, ICAM2, CHGB, GALP, ICOS, RELA, and CSH2 ( Figure 4A ). In the TCGA cohort, the OS-prognostic model contained four genes: CCL22, FABP7, LTBR, and RBP7 ( Figure 4B ). Risk scores were constructed with the regression coefficients from these models and thresholds were chosen manually at the median. In the DFS-prognostic model, high-risk patients had significantly worse DFS (P <0.0001) ( Figure 5A ). In the OS-prognostic model, high-risk patients had statistically notable worse OS (P <0.0001) ( Figure 5B ).

The ESTIMATE algorithm was used to determine the scores of each sample in dataset GSE69657 by R software. The scores showed a remarkable difference between the 5-FU resistant and the 5-FU sensitive groups. 5-FU resistant patients showed lower scores and were statistically significant in StromalScore and ESTIMATEScore but not in ImmuneScore ( Figures 6A–C ). Furthermore, we compared the differentially expressed genes between high StromalScore and low StromalScore groups to obtain tumor microenvironment-related genes. We found two gene sets which contained 988 genes: 355 down-expressed genes and 633 up-expressed genes ( Figure 6G ).

The potential function of the tumor microenvironment-related genes to 5-FU resistance was performed by GO and KEGG enrichment analysis. For BP enrichment, the genes most enriched in organelle fission, extracellular matrix organization, and extracellular structure organization are shown in Figure 7A . For CC enrichment, the genes most enriched in the collagen-containing extracellular matrix and spindle are shown in Figure 7A . For MF enrichment, the genes most enriched in the extracellular matrix, tubulin binding, and actin-binding are shown in Figure 7A . For KEGG enrichment, large portion of the genes concentrated on proteoglycans in cancer, focal adhesion, cell adhesion molecules, staphylococcus aureus infection, phagosome, complement and coagulation cascades, hematopoietic cell lineage, viral myocarditis, viral protein interaction with cytokine, cytokine receptor, and cell cycle are shown in Figure 7B ).

We obtained immune-related genes of 5-FU resistance genes in the tumor microenvironment by intersecting previously screened, immune-related drug resistance genes with tumor microenvironment-related genes ( Figure 6I ). We gained two overlapping genes RARB and RBP7 and further verified the expression of these two genes in 5-FU-sensitive and -resistant patients in the other three GEO datasets. RBP7 was up-regulated in patients with 5-FU- resistance and was statistically significant in GSE19860 ( Figure 6D ). However, RBP7 was down-regulated in CRC tissue compared with normal tissue in the TCGA cohort ( Figures 6E, F ). In addition, in GSE19860, the IC50 of 5-FU predicted by "pRRophetic" in patients with high expression of RBP7 was higher than patients with low RBP7 expression, but the results were not statistically significant(p-value=0.33)( Figures 6H ). We defined the RBP7 as a tumor immune microenvironment-related gene for CRC to 5-FU resistance.

We selected the median as the cut-off to divide the groups into the RBP7-high expressed and RBP7-low expressed groups. We analyzed the relationship between the expression of RBP7 and the clinical characteristics of CRC patients and the results showed that the expression of RBP7 was closely correlated with T stage, but not with age, N stage, M stage, or mismatch repair gene deletion ( Table 1 ). Although RBP7 was not associated with lymphatic invasion, it did mean for the number of lymph nodes. Also, RBP7 expression was linked to survival status.

GSEA analysis showed 36 significant KEGG pathways were associated with RBP7 expression. The top 10 percent of pathways included Herpes simplex virus 1 infection, PI3K-Akt signaling pathway, and cytokine-cytokine receptor interaction. GO enrichment for 593 gene sets was performed. 480 gene sets were enriched in the BP process, 65 gene sets accessed in the CC enrichment, and only 48 gene sets were acquired in the MF enrichment. The top 30 significant GO term and KEGG pathways are shown in Figures 8A, B .

In order to further obtain information related to the function of RBP7 gene, we selected GSE19860 dataset for GSVA analysis to obtain enrichment of downstream pathways related to RBP7. We found that two pathways gene sets were significantly altered in chemotherapy-sensitive (RBP7^Low) and chemotherapy-resistant (RBP7^High) patients, and those were Hypoxia and TNFα signaling via NFκB gene sets ( Figure 8C ).

To understand how the tumor cells influence the immune microenvironment, Lee et al. analyzed the transcriptome of 91,103 single cells from 23 Korean and 6 Belgian patients. Their result showed that intercellular network reconstruction supported the link between cancer cell signatures and specific matrix or immune cell populations (22). Based on the Lee et al. study results, we preliminarily understood the distribution of RBP7 in colorectal cancer. RBP7 was centrally distributed in endothelial stalk cells, endothelial tip cells, and myeloid cells ( Figures 8D, E ).

To evaluate whether RBP7 mediates the response of CRC cells to 5-FU treatment, we preliminarily verified in cell experiments that both transcription level and protein level of RBP7 were increased in 5-FU resistant cells ( Figure 9 ), further suggesting that RBP7 might induce 5-FU resistance in CRC.