Introduction
Colorectal cancer (CRC) ranks as the second most common cause of cancer-related mortality globally [
1]. Advances in the understanding of CRC pathophysiology have expanded the available therapeutic options, such as immunotherapy [
2,
3]. However, the 5-year survival rate of people with CRC metastasis is approximately 14% [
4]. Uncontrolled carcinogenic events promote genetic mutations and epigenetic modifications, ultimately resulting in the development of CRC. Genomic instability plays a crucial role in carcinogenesis [
5]. Due to molecular heterogeneity, the relapse and mortality rates of CRC might significantly differ among individuals with identical clinicopathological characteristics [
6,
7]. Recent genomics advancements have enabled the TCGA Research Network to characterize primary subtypes of CRC [
8].
The utilization of single-cell RNA sequencing (scRNA-seq) has shown the potential for investigating the heterogeneity of different malignancies.
Zhang et al
. used scRNA-seq to analyze differentially expressed genes (DEGs), distribution, and T cell receptor profiles of T cells from patients with CRC [
9]. scRNA-seq facilitates analyses of cellular diversity at the individual cell level and investigations of the functions and disease processes that are associated with distinct cell clusters [
10‐
12]. Via scRNA-seq analysis,
Wang et al
. confirmed that monocytes and macrophages significantly influence the tumor microenvironment (TME) of CRC [
13]. However, none of these studies were verified via bulk RNA-seq. Bulk RNA-seq provides a comprehensive overview of transcription, although this method has limitations in accurately identifying individual cell types and determining intratumoral heterogeneity [
14]. Due to the inconsistencies that have been observed between subtypes that are classified based on a single molecule and clinical outcomes [
15‐
17], combinations of key molecular markers have been proposed to elucidate relationships between molecular events and clinical outcomes and to increase the accuracy of prognosis [
18‐
21].
Luo et al
. used scRNA-seq and bulk RNA-seq to molecularly categorize CRC based on necroptosis and to make prognostic predictions [
22].
Joanito et al
. identified two distinct states of epithelial tumor cells and improved the consensus molecular categorization of CRC via scRNA-seq and bulk RNA-seq analyses [
23]. In addition, tumor mutational processes continuously influence the somatic genome, resulting in immunodeficiency, aging, and other diseases.
Thyroid hormone receptor-interacting protein 6 (TRIP6) is an adapter protein that utilizes its LIM domain to engage in interactions with numerous proteins [
24,
25]. Previous studies have demonstrated that activation of TRIP6 may be induced by the transcriptional activating factor v-rel and acts as a promising target for inhibiting v-rel-induced tumor formation and transcriptional activity [
26]. TRIP6, a multifunctional protein, regulates a multitude of biological processes associated with diverse types of cancer; for example, TRIP6 promotes carcinogenesis by enhancing the malignant proliferation and invasion of cancer cells [
27,
28]. The overexpression of TRIP6, which is upregulated in glioma cells and tissues, is associated with unfavorable clinical outcomes among glioma patients [
29]. However, there is still not enough known about TRIP6, especially what role it plays in CRC.
In this study, we revealed main cell subtypes of CRC that are associated with the surrounding tissues. Furthermore, we established a CRC model that is associated with prognostic genes. By conducting in vitro cell experiments and bioinformatics analyses of data from multiple public databases, we have identified the expression patterns and function of TRIP6 and proposed its potential role as an oncogene in CRC.
Materials and methods
Data sources and processing
From the TCGA COAD and READ cohorts, bulk RNA-seq data and clinical information of patients with CRC were downloaded, which included 51 normal colorectal mucosa samples and 647 CRC tissue samples. The scRNA-seq dataset (GSE161277) [
30] was obtained from the GEO database (
https://www.ncbi.nlm.nih.gov/geo/) and contained data of three paracancerous tissues and eight tumor tissues from patients with CRC. The GSE17536 and GSE39582 [
31,
32] datasets (containing data of 177 and 556 patients, respectively, with complete expression profiles and survival information) were also downloaded from the GEO and used as external validation datasets.
scRNA-seq data analysis
The samples were merged utilizing the R package "Seurat" [
33], and cells were isolated by filtration using scRNA-seq with the following exclusion criteria: low expression levels of genes were identified only with nFeature_RNA values > 50 and percent.mt values < 5.
Gene expression levels in isolated cells were standardized using a linear regression model, and the 10 genes with the greatest variability in expression were identified by analysis of variance (ANOVA). The dimensions of the scRNA-seq data were reduced by performing principal component analysis (PCA), and 18 PCs were selected for subsequent analysis using ElbowPlot [
34]. The relative location of each cluster was established by conducting t-distributed stochastic neighbor embedding analysis. We annotated the clusters using cell marker genes in the CellMarker database [
35] to determine possible positional relationships between cells. We used HumanPrimaryCellAtlasData [
36] as a reference for supplementary annotation. Subsequently, we identified marker genes for each cell subtype by analyzing the single-cell expression profiles by adjusting the logfc.threshold argument of the FindAllMarkers function [
37] to a value of 1. Finally, genes with a log
2-fold change > 0.585 and a
P < 0.05 were selected as unique marker genes of each cell subtype based on our data and on other existing literature [
38,
39].
Construction and validation of a prognostic model
Univariate Cox regression analysis was conducted on the candidate genes in the training set to identify genes that are specifically linked to prognosis. Variables with
P values < 0.05 were analyzed by LASSO regression, which was conducted using the "glmnet" R package [
40] to minimize the number of genes in the risk model. The prognostic model was constructed using the subsequent equation: risk score = gene exp1 × β1 + gene exp2 × β2 + … + gene expression n × βn (gene expression represents the numerical number of the expression levels, and β represents the coefficient obtained via LASSO regression analysis). The R package "survminer" [
41] was used to construct survival curves. The R package "survROC" [
42] was utilized to generate ROC curves for assessing the risk scores in predicting OS at 1, 3, and 5 years. The prognostic model's validity was validated through internal and external datasets.
Analysis of immune cell infiltration
The RNA-seq data were processed utilizing the R package “CIBERSORT” [
43] to determine the proportions of 22 different infiltrating immune cell types. Associations between risk scores and the presence of tumor-infiltrating immune cells were analyzed using Pearson correlation coefficients.
Gene set enrichment analysis (GSEA)
GSEA [
44] (
http://www.broadinstitute.org/gsea) was conducted on all the DEGs in the TCGA dataset using the clusterProfiler package [
45]. After performing 1000 permutations, we identified enriched gene sets with a
P value < 0.05 and a false discovery rate of 0.25. In conclusion, we conducted enrichment analyses using the Kyoto Encyclopedia of Genes and Genomes (KEGG) [
46] and Gene Ontology (GO) [
47] to illustrate the functional pathways that differentiate high-risk and low-risk groups.
Gene set-variant analysis (GSVA)
GSVA [
48] was utilized to assess gene set enrichment in the transcriptome data. The gene sets were acquired from the Molecular Signatures Database (MSigDB) [
49] (
http://www.gsea-msigdb.org/gsea/index.jsp), and limma software was utilized to perform GSVA, which can accurately assess possible differences in biological function across various samples.
Analysis of chemotherapy sensitivity
We predicted the sensitivity of each CRC sample to chemotherapy using the Genomics of Drug Sensitivity in Cancer (GDSC) database [
50] (
https://www.cancerrxgene.org/) and the R package “pRRophetic” [
51]. We screened drugs using Spearman’s analysis to evaluate correlations between the IC
50 values and risk scores of various medicines.
Tumor mutational burden (TMB) analysis
We analyzed TCGA mutation annotation data using the "Maftool" R package [
52]. The difference in TMB between the two risk groups was assessed. In addition, waterfall plots were generated to visualize mutations in the 30 genes that were most frequently different between the low- and high-risk groups. After performing SubMap analysis, the results were visualized with the “complexHeatmap” R package [
53].
Cell communication analysis
The CellPhoneDB database [
54] (version 4.0) was used to examine ligand‒receptor interactions by analyzing the single-cell expression patterns of relevant molecules. Additionally, separate analyses of differences in cell communication between the high-risk and low-risk groups were conducted.
Random survival forest (RSF) analysis
The RSF algorithm was applied to assess the significance of associated genes via the R software “randomForestSRC” [
55]. This algorithm involved 1000 iterations with Monte Carlo simulation (
n rep = 1000). Genes that had a relative relevance greater than 0.3 were selected for inclusion in the final signature.
Clinical specimens and immunohistochemistry (IHC)
The harvested specimens were stored in a solution specifically designed for preserving tissues (Miltenyi Biotec, catalog number 130–100-008) from February 2020 to February 2023 from CRC patients in the First People's Hospital of Foshan. All the specimens were cut into pieces that were approximately 1 mm × 1 mm in size and incubated in preheated digestion solution with an enzyme mixture (Sigma‒Aldrich) at 37 °C for 30 min. The tissue samples were fixed, paraffin-embedded, dewaxed, rehydrated, and subjected to antigen retrieval. The samples were stained with primary antibodies at 4 °C overnight and then incubated with a suitable secondary antibody for 30 min at 37 °C. Next, the samples were visualized using DAB solution. Finally, images were captured using an optical microscope. Quantitative analysis was performed on five representative images at a magnification of 40 × using ImageJ software.
Cell culture and transfection
The NCM460, LoVo, HCT-116, and SW620 cell lines were acquired from the Shanghai Cell Bank of the Chinese Academy of Sciences (Shanghai, China). The cells were cultivated in DMEM medium with 10% FBS (Life Technologies, USA) and penicillin (100 U/mL) and streptomycin (100 U/mL) in a humidified environment at 37 °C with 5% CO
2. A scrambled shRNA and TRIP6 shRNA were purchased from GeneChem (Shanghai, China). The transfections were performed using Lipofectamine 3000 (Invitrogen, Carlsbad, CA) according to the manufacturer’s instructions [
56].
Western blotting analysis
Western blotting assays were conducted using a previously published protocol [
57]. Total proteins were extracted from CRC samples with RIPA buffer (Thermo Fisher Scientific, USA). The protein concentrations were measured with a BCA assay kit (Thermo Fisher Scientific, USA). An equivalent amount of protein from each sample was separated using SDS‒PAGE. Then, the proteins were transferred to polyvinylidene membranes (Millipore; Burlington, MA, USA). The membranes were blocked with Tris-buffered saline-0.1% Tween-20 (TTBS) supplemented with 5% skim milk for 2 h at 25 °C. The membranes were then incubated at 4 °C overnight with the following primary antibodies at the indicated dilutions: anti-CYP2W1 (1:500, PA5-101315, Invitrogen), anti-GDE1 (1:1000, PA5-43012, Invitrogen), anti-PTPN6 (1:1000, ab124942, Abcam), anti-PTTG1IP (1:500, ab128040, Abcam), anti-SEC61G (1:500, PA5-21384, Invitrogen), and anti-TRIP6 (1:500, ab137478, Abcam) antibodies. Next, the membranes were incubated with a secondary antibody, anti-rabbit IgG-horseradish peroxidase conjugate (1:4000; Cell Signaling Technology) for 1 h at 25 °C. The expression levels of the target proteins were detected using an ECL kit (Bio-Rad, USA) and analyzed via ImageJ software ((NIH V1.8.0.112, USA).
Reverse transcription-qPCR (RT-qPCR) analysis
TRIzol reagent (Roche) was used to extract total RNA from cultivated cells. Next, reverse transcription was performed on each sample using the PrimeScript™ RT Reagent Kit with gDNA Eraser (Takara). Gene expression was quantified with TB Green® Premix Ex Taq™ II. The primer sequences that were utilized are provided in Table
S1. Relative gene expression levels were determined utilizing the comparative 2
−ΔΔCT method [
58].
The cells were evenly distributed in 96-well plates. After a 96-h incubation period, the viability of cells was assessed according to the manufacturer’s instructions using a CCK-8 assay. For colony formation assays, a total of 1 × 103 LoVo cells were incubated in 6-well plates for 2 weeks. Subsequently, the cells were immobilized using a 4% paraformaldehyde (PFA) solution for 15 min, followed by staining with 0.2% crystal violet for 15 min. The images were acquired with a digital camera and processed using ImageJ software.
Wound-healing and Transwell assays
The cell monolayer was linearly wounded using a pipette tip, and cells were cultured in serum-free medium using a pipette tip in wound-healing assays. The wound-healing rate was imaged and examined after 48 h. In a 24-well transwell cell culture apparatus containing 8-μm pore size multiporous polycarbonate membrane inserts, migration assays were performed. Briefly, the upper chamber was filled with a suspension of cells, while the bottom chamber was filled with medium that included 10% FBS. The filters were removed, rinsed twice with PBS, fixed with methanol, and stained with 0.5% crystal violet reagent following a 24-h incubation at 37°C in 5% CO2. The cells located on the upper side of the filter were removed using cotton swabs. In order to identify migrated cells on the lower side of the filter, particular cross-sectional fields on the filters were quantified. Transwell invasion assays were conducted using Matrigel-coated transwells under the same conditions.
Flow cytometry assay
LoVo cells were treated with prechilled 70% ethanol for 12 h, incubated with RNase in a 37 °C water bath for 30 min, and prepared with a propidium iodide (PI) detection kit (KeyGen, China). The cell nuclear DNA was labeled with PI at 4 °C. The cells were examined using flow cytometry within 1 h following the manufacturer's recommendations to analyze cell cycle distribution. To analyze apoptosis, the LoVo cells were trypsinized and collected by centrifugation at 500 × g for 3 min at 4 °C. Afterward, the cells were suspended again using a FITC Annexin V Apoptosis Detection Kit (BD Biosciences). The data were analyzed using ModFit-LT software (Verity Software, Topsham, ME) and FlowJo V10 software (Treestar, Inc. Ashland OR).
5‐Ethynyl‐2′‐deoxyuridine (EdU) assay
After transfection, LoVo cells were seeded in 12-well plates with 14 µl of slippers and cultivated with 50 µm of EdU reagent (diluted 1:1000 in DMEM with 10% FBS) for 2 h at 37 °C. The cells then were fixed with 4% paraformaldehyde (PFA) and stained with Hoechst solution.
Discussion
The prognostic value of the heterogeneous and aggressive characteristics of CRC remains limited [
59]. It is crucial to evaluate innovative risk models or subtype-specific risk variables to assist in the development of patient-specific treatments and improve patient prognosis [
60]. scRNA-seq is used to characterize distinct cell populations and identify unique biomarkers and variations within different cell types in many cancers, including CRC [
61,
62]. Using scRNA-seq, Poonpanichakul et al. discovered that cancer subpopulations with significant heterogeneity exhibit distinct responses to a specific chemotherapy [
63]. In contrast, bulk RNA-seq reveals mean gene expression across all cells [
64]. To construct a novel risk model, we performed an exhaustive analysis of bulk RNA-seq and scRNA-seq data in this study that demonstrates good efficacy in predicting prognosis and determining the response to immunotherapy in patients with CRC.
Zheng et al. verified a nine-gene profile associated with cancer-associated fibroblasts (CAFs) that may be used as a standalone predictive marker for patients with CRC who might not benefit from immunotherapy, but this study did not consider intercellular communication [
65]. The heterogeneity of tumor samples also results in increased variation in cellular communication, and heterogeneity of interactions with the TME is essential for tumorigenesis and resistance to therapy [
66]. Our results suggest that intercellular communication is an important mechanism among cell subtypes. GO and KEGG analyses of DEGs that were obtained from TCGA revealed significant enrichment of genes related to the cell cycle, EMT, and KRAS signaling pathways, which contribute to the development of CRC. Alterations in cyclins and EMT-related genes are common in CRC, and therapeutic intervention that targets aberrant cell cycle regulators may be advantageous in the treatment of CRC [
67]. Multiple studies have confirmed that KRAS-related signaling pathways play crucial roles in CRC development [
68].
The high-risk group showed significant enrichment in immunological processes. Thus, we hypothesized that the risk score could function as a possible prognostic marker for patients with CRC who are receiving immunotherapy. We evaluated several aspects of TMB, immune infiltration, and immune checkpoints. The results showed that the expression of various immune factors, including immunomodulators, chemokines, and cell receptors, significantly differed between the high- and low-risk groups, suggesting that patients with different risk scores may exhibit differences in their responses to and ability to benefit from targeted immunotherapy. Furthermore, the results of correlations between the model and clinicopathological features demonstrated that the model was strongly correlated with lymph node metastasis and tumor stage, suggesting that the predictive capability of the model also predicts overall survival.
Furthermore, druggable targets in patients with CRC and their corresponding drugs in the GDSC database were discovered using our prognostic models. The aforementioned findings supported the DEG enrichment outcomes. The high-risk group demonstrated greater drug resistance to gemcitabine, doxorubicin, mitomycin C, and bleomycin, whereas the low-risk group demonstrated significantly reduced drug resistance to docetaxel and metformin, indicating that people with low risk scores may be more likely to benefit from chemotherapy.
The six prognosis-related genes,
TRIP6, PTTG1IP, PTPN6, GDE1, SEC61G, and
CYP2W1, were screened using RSF analysis. We found that
TRIP6, PTPN6, and
GDE1 were independent prognostic factors for OS of CRC patients. Zhang et al. established subcategories for CRC based on the cGAS-STING pathway. These subcategories may be utilized to predict patient prognosis using 27 DEGs but without identifying DEG-associated ceRNAs [
69]. The regulatory networks involved play crucial roles in CRC progression. Therefore, we developed a ceRNA regulatory network to elucidate the potential mechanisms underlying CRC using six crucial DEGs. This approach can provide insights into the unidentified regulatory networks involved in CRC. In addition, based on the TCGA analysis, TRIP6 was highly expressed in CRC, and its high expression reduced the overall survival of CRC patients. Therefore, we performed further experiments focusing on
TRIP6.
The patterns of gene expression and heterogeneity identified through bulk RNA-seq and scRNA-seq have often not been validated experimentally [
70]. We first examined and confirmed to be highly expressed in CRC tissues by immunohistochemistry, RT-qPCR, and western blotting.
TRIP6 encodes a protein with three LIM (Lin-11, Isl-1, and Mec-3) zinc-binding domains [
71]. TRIP6 activates Akt signaling to facilitate CRC drug resistance by directly interacting with PARD3 [
72]. Moreover, by targeting TRIP6, miR-7 inhibits the proliferation and migration of CRC cells [
73]. However, the investigation of TRIP6's role in CRC remains unexplored. In this study, the inhibition of TRIP6 in LoVo cells led to a substantial inhibition of cell proliferation and growth, induction of cell cycle arrest, and facilitation of tumor cell apoptosis. This suggests that TRIP6 induces apoptosis and inhibits the cell cycle, thereby promoting cell survival. Molecular network analysis suggested that TRIP6 expression levels were correlated with that of some oncogenes and pathways related to carcinogenesis in CRC (Fig.
6C, D). For example, the inhibition of liver metastasis in colon cancer is significantly observed upon suppression of the oncogene AKT1. Additionally, epithelial-to-mesenchymal transition (EMT) is a pivotal factor in facilitating CRC invasion and metastasis. Tumor metastasis is the result of cancer cell migration and invasion. Our results indicated the inhibitory effects of TRIP6 silencing on LoVo cell invasion and migration in vitro. Therefore, TRIP6 may be a possible therapeutic target for CRC treatment.
Conclusions
In conclusion, our research exhaustively characterized the various cell subpopulations found in the colorectal cancer tissues. Key molecular mechanisms and potential therapeutic implications of targeting immune modifications in colorectal cancer were identified in our study. Moreover, TRIP6 expression is upregulated in CRC. TRIP6 promotes cell proliferation, migration, invasion, cell cycle dysregulation, inhibition of apoptosis, which may serve as a prognostic indicator for CRC. Clearly, further exploration of the use of TRIP6 as a therapeutic target for CRC is warranted.
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.