Prognostic and therapeutic potential of gene profiles related to tertiary lymphoid structures in colorectal cancer

Formulation of the prognostic TLSRGs signature

We initiated our analysis by amalgamating the bulk RNA-seq datasets, rectifying batch effects using the “limma” and “sva” packages. Subsequently, a log2 transformation was applied to standardize the data. The relationship between TLSRGs expression levels and patient overall survival (OS) was scrutinized using univariate Cox regression analysis, utilizing the “survival” package in R (version 3.2–13) (https://CRAN.Rproject.org/package=survival).

Consensus clustering for prognostic TLSRGs

We used the R package ‘ConsensusClusterPlus’ to perform unsupervised clustering analysis on COADREAD patients, stratifying them based on the expression of prognostic-TLSRGs. By employing the K-means algorithm, we identified the optimal number of clusters ranging from 2 to 10. This process was reiterated 1,000 times to affirm the stability and reliability of our results.

Clinicopathological traits

Investigating the clinical implications, we analyzed the correlation between the molecular subtypes and clinicopathological features, subsequently visualized in a heatmap. The clinicopathological features that were examined included the patient’s age, sex, tumor grade, and stage of tumor node metastasis (TNM). The “survminer” packages in R facilitated our survival analysis.

Pathway enrichment analysis

To elucidate the TIME characteristics in various molecular subtypes, GSVA was employed utilizing gene sets from hallmark, KEGG, and Reactome obtained from the MSigDB database.

Immune infiltration

The stromal score, immune score, and ESTIMATE score in tumor tissues were quantified using the “estimate” R package (Scire et al., 2023) by employing the ESTIMATE algorithm from the ESTIMATE project (https://sourceforge.net/projects/estimateproject/). Furthermore, we utilized ssGSEA to delineate the immune infiltration pattern in COADREAD by employing marker genes for 23 immune cell types (Hänzelmann, Castelo & Guinney, 2013).

Analysis of differentially expressed genes

The “limma” package was utilized to identify differentially expressed genes (DEGs) among clusters, specifically targeting those with an adjusted p-value < 0.05 and |log2 (Fold Change)| >1. The “clusterProfiler” R package was used to perform Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses. To discern DEGs with prognostic significance, univariate Cox regression analysis was executed, with results presented in a forest plot. Additionally, the GSCA datasets (https://guolab.wchscu.cn/GSCA/#/) (Liu et al., 2023) were used to perform mutation, SNV, CNV, and Methylation analyses on these prognostic genes.

Refinement of the TLS score system

We employed the “ConsensusClusterPlus” package in R for unsupervised consensus clustering to categorize COADREAD patients into distinct molecular subtypes. The classification was determined by analyzing the expression patterns of previously identified prognostic TLS-DEGs. The clustering process involved the K-means algorithm, assessing optimal cluster numbers from k = 2–10. This procedure was iteratively performed 1,000 times for robust stability of results. Kaplan-Meier plots elucidated significant prognostic variations between the clusters. Additionally, a heatmap delineated the interplay between clinical attributes, gene expression profiles, and clustering outcomes. The TLS score, computed via the GSVA package using 16 pivotal prognostic TLS-DEGs, facilitated subsequent survival analysis. This analysis discerned the prognostic predictive capability of the TLS score in COADREAD cases. A Sankey diagram illustrated the intricate associations among signature genes, scoring, and prognostic categorization.

Comprehensive immune landscape analysis utilizing the TLS score

To further understand the TLS score’s implications, we created a correlation heatmap showcasing its association with immune cell infiltration levels. Moreover, the interrelations between the TLS score and various cytokines, chemokines, and their receptors were visualized through additional heatmaps. The GSVA method was pivotal in correlating the TLS score with 50 hallmark pathways, with the findings depicted using R’s “ggplot2” package. We also gathered and analyzed mRNA expression data of immune checkpoints and gene mutation information from the TCGA-COADREAD cohorts. This led to the calculation of mean normalized expression levels of immune checkpoints per sample, subsequently normalized to convey their relative expression magnitudes.

Analysis of MSI, TMB, CSCs, and somatic mutations

We conducted an in-depth analysis of MSI, CSC markers, and TMB across varying TLS score groups. For this, somatic mutation data specific to COADREAD patients were extracted from the TCGA database. To discern somatic mutations between different TLS score groups, we utilized the “maftools” package to construct detailed waterfall plots. A comparative analysis of the somatic mutation variances across the TLS score groups was then conducted using forest plots.

Evaluation of immunotherapy response and chemotherapy sensitivity

The predictive relevance of the TLS score for immunotherapy response was investigated using transcriptomic data from two distinct cohorts. The dataset GSE135222 includes information from patients with advanced non-small cell lung carcinoma who received treatments targeting anti-PD-1/PD-L1. Meanwhile, GSE176307 includes transcriptomic information from patients with metastatic urothelial cancer. The prognostic significance of the TLS score in COADREAD was evaluated using Kaplan-Meier analysis in the R package. Moreover, by utilizing the CancerRxGene database (https://www.cancerrxgene.org/) (Yang et al., 2013), we calculated the IC50 values for different target-therapeutic medications in relation to the TLS score. This estimation was performed using the “pRRophetic” package in R, offering insights into the chemotherapy sensitivity linked to the TLS score.

RNA extraction and quantitative real-time PCR

Total RNA from colorectal cells was isolated employing the Trizol method (Invitrogen, Waltham, MA, USA) as per the manufacturer’s instructions. Subsequent reverse transcription to synthesize cDNA was conducted using Prime Script RT reagent Kit (RR047A, Takara, Shiga, Japan). Gene expression was measured using TB Green Premix Ex Taq (RR420A, Takara, Shiga, Japan), and GAPDH was used as the reference for normalization. Details of the primers used are provided in Table S1.

Cultivation of cell lines and reagent preparation

The Chinese Academy of Sciences Cell Bank (Shanghai, China) provided several colorectal cancer cell lines, namely NCM460, HCT116, HT29, and LOVO. Each cell line underwent STR profiling for authentication and was confirmed to be free of mycoplasma contamination (https://www.cellbank.org.cn/). NCM460 and HCT116 were cultured in RPMI 1640 medium (Gibco, Waltham, MA, USA) with the addition of 10% FBS, whereas HT29 and LOVO were maintained in DMEM (Gibco, Waltham, MA, USA) with a 10% FBS enrichment. Cultures were incubated at 37 °C in a humidified atmosphere with 5% CO2.

Targeted knockdown of TLS-related genes

To conduct gene silencing experiments, HCT116 and HT29 cells were plated at a density of 3 × 10⁵ cells per 60 mm dish. Following 24 h of incubation, the medium was replaced with fresh growth medium. Next, the cells were transfected with siRNAs that targeted C5AR1, APOE, CYP1B1, and SPP1, or a control siRNA obtained from Genepharma in Shanghai, China. This transfection was performed using Lipofectamine RNAiMax (13778075; Life Technologies, Carlsbad, CA, USA). After transfection, the cells were cultured in RPMI 1640 or DMEM, with the addition of 10% FBS, for at least 24 h. The siRNA sequences are listed in Table S1.

Cell proliferation assessment

Cells were seeded in 96-well plates at a concentration of 3,000 cells per well. The growth of cells was observed for a period of 24 h using a Cell Counting Kit-8 (CCK-8) from Dojindo, following the manufacturer’s provided instructions. A microplate reader (Thermo Fisher, Waltham, MA, USA) was used to acquire absorbance measurements at 450 nm. GraphPad Prism version 9.5.0 was utilized for conducting data analysis.

Transwell migration assay

To perform migration experiments, a total of 100,000 cells were placed in the top compartment of a Transwell device on a Matrigel layer (Corning, Corning, NY, USA), utilizing either RPMI 1640 or DMEM excluding FBS. The bottom section contained 600 μL of the corresponding solution enriched with 20% FBS. After a 24-h incubation period, cells were fixed and stained with crystal violet. Cells remaining in the upper chamber were carefully removed, and the migrated cells were imaged for further analysis.

Statistical methodology

The statistical analysis was conducted using R software (version 4.1.2). Correlation assessments were conducted using Pearson or Spearman correlation coefficients. The Wilcoxon test was applied for comparisons between two groups. The overall survival (OS) across various subgroups was analyzed using Kaplan-Meier survival curves and log-rank tests. Univariate Cox regression analyses were utilized to estimate prognostic values. The data from the qRT-PCR analysis were examined utilizing the t-test of Student. A p-value less than 0.05 was deemed to be statistically significant, with an asterisk (*) representing p-values less than 0.05, two asterisks (**) representing p-values less than 0.01, and three asterisks (***) representing p-values less than 0.001.

Ethical consideration

The authors conducted this study without the participation of any humans or animals.

Dimensionality reduction

In our initial step, the single-cell sequencing data from the CRC_EMTAB8107 dataset, were subjected to an integration analysis to consolidate information across different samples, Fig. S1 exhibited the minimal batch effects and 19 distinct clusters categorized by the UMAP algorithm. Further, we examined the expression of specific surface marker genes across these clusters. This analysis enabled the identification of 10 unique cell types, encompassing malignant cells, CD4+ T cells, plasma cells, fibroblasts, CD8+ T cells, monocytes-macrophages, B cells, endothelial cells, mast cells, and dendritic cells, as depicted in Fig. 2A.

Landscape of transcriptome

In the realm of gene expression, Fig. 2B highlighted the top five genes making the most significant contributions. Specifically, in malignant cells, the genes exhibiting the highest expression levels included C19orf33, KRT19, TSPAN8, AGR2, and LGALS4. In stark contrast, genes such as SRGN, RGS1, TSC22D3, CXCR4, and CD52 were among the least expressed. Within CD8+ T cells, genes like CCL5, NKG7, GZMA, CCL4, and GNLY were highly expressed, while others including CALD1, CST3, SPARCL1, GSN, and IFITM3 showed minimal expression. For monocytes-macrophages, the top expressed genes were AIF1, TYROBP, C1QA, FCER1G, and C1QB, whereas genes like TRAC, CALD1, IL32, CD69, and FKBP11 were the least expressed.

Pathway enrichment analysis

The GSVA brought to light hallmark pathways that were predominantly enriched in each cell type (see Fig. 2C). In malignant cells, pathways such as mTORC1 signaling, glycolysis, xenobiotic metabolism, fatty acid metabolism, were notably prevalent. Fibroblasts showed a significant enrichment of the epithelial mesenchymal transition pathway. The monocytes-macrophages were enriched in pathways such as complement, IL2-STAT5 signaling, inflammatory response, and interferon response pathways.

Cell-type specific expression of TLS insights

Employing the “FindAllMarkers” function, a total of 39 TLSRGs were analyzed. This detailed examination culminated in deciphering the expression patterns of TLSRGs across various cell types (Fig. 2D). Focusing on cell-specific expression, MS4A1 emerged as a gene predominantly expressed in B cells, while its expression remained minimal in other cell types. In the realm of CD8+ T cells, genes like CCL4 and CCL5 showcased high expression levels, contrasting with their subdued expression in other cells. Additionally, TNFRSF17 was notably expressed in plasma cells, yet exhibited low expression in other cell types.

GSVA-based TLS score calculation

Utilizing the UMAP cell annotation map, we calculated the TLS score for each cell through the GSVA method, as illustrated in Fig. 2E. Subsequently, cells were categorized into clusters with low and high TLS scores based on their median expression levels, depicted in Fig. 2F. Predominantly, malignant cells and fibroblasts were found in the low TLS score cluster, whereas T cells and monocyte-macrophages were more prevalent in the high TLS score cluster.

Cell cluster distribution and TLS score visualization

This distribution was further detailed through box plots representing TLS scores across different cell clusters (Fig. 2G), where CD8+ T cells received the highest scores, and mast cells the lowest. The distribution patterns of cell types within these clusters were effectively represented in a composition chart based on cell count and proportion (Figs. 2H and 2I), highlighting the tendency of malignant cells and fibroblasts to fall into the low scoring group, in contrast to monocyte-macrophages and CD8+ T cells which were more commonly found in the high scoring group.

TLS score correlation with hallmark pathways

Further analysis revealed a significant correlation between the TLS score and specific hallmark pathway scores (Fig. 2J). In malignant cells, the TLS score showed a positive correlation with pathways such as apoptosis and KRAS signaling up. In monocyte-macrophages, a positive correlation was observed with the Inflammatory response, interferon gamma response, IL2-STAT5 signaling pathways. For dendritic cells, the TLS score positively correlated with pathways including Inflammatory response, Interferon gamma response, Interferon alpha response. We identified hallmark pathways that were differentially expressed between high- and low-TLS score clusters, as shown in Fig. 2K. These enriched pathways in high-score group included inflammatory response, interferon gamma response, interferon alpha response, markedly distinct compared to the low-TLS score group.

TLS molecular subtypes: identification and analysis

Our exploration focused on discerning the TLS patterns in COADREAD tumorigenesis. To this end, we examined 1,184 COADREAD patient samples sourced from the TCGA and GEO databases (including GSE39582 and GSE17538). Survival analysis pinpointed 12 TLSRGs (CCL2, CCL20, CXCR3, IL1R1, TIGIT, SGPP2, ICOS, CCL8, CXCL9, CXCL11, CXCL10, CXCL13) significantly correlated with the overall survival (OS) in COADREAD patients (Fig. 3B, p < 0.05). We then constructed a TLS network to delineate the interrelationships among TLSRGs and their prognostic values (Fig. 3A). Univariate Cox regression analysis on these genes revealed positive associations of CCL20, CXCR3, TIGIT, SGPP2, ICOS, CXCL9, CXCL11, CXCL10, and CXCL13 with patient survival, while CCL2, CCL8, and IL1R1 showed negative associations.

Based on the 12 prognostic TLSRGs, we employed a consensus clustering algorithm to divide COADREAD patients into two distinct molecular subtypes, namely Cluster A and Cluster B, with the analysis suggesting k = 2 as the optimal cluster number (Fig. 3D). Patients in subtype B exhibited significantly higher survival rates compared to those in subtype A, as demonstrated by the log-rank test (p = 0.001; Fig. 3C). A notable divergence in the expression patterns of TLSRGs was observed between clusters A and B (Fig. 3E). Particularly, a subset of TLSRGs, including CXCL9, CXCL10, CXCL11, CXCL13, ICOS, TIGIT, among others, were found to be predominantly upregulated in cluster B (p < 0.05). Further, we profiled the associations of TLSRG expression with various clinical characteristics such as age, disease stage, gender, recurrence, metastasis, fustat, and futime across the different subtypes (Fig. 3F). This profiling provided insights into the clinical relevance of TLSRG expression variations in distinct patient subgroups. To enhance the distinct biological traits of the subtypes, our study incorporated a comprehensive GSVA enrichment analysis using Hallmark, KEGG, and Reactome pathway databases. Collectively, these comprehensive pathway analyses consistently underscored the prominence of interferon gamma response, interferon alpha response, chemokine signaling pathway, and cytokine-cytokine receptor interaction pathways in subtype B. These pathways were critically linked to immune activation (Yenyuwadee et al., 2022), further delineating the distinct immunological landscape of subtype B (Figs. 3G–3I).

The PCA analysis reinforced the significant transcriptional distinctions between subtype A and B (Fig. 3J). To delve deeper into the role of TLSRGs in shaping the TIME, we quantified human immune cell subsets in each COADREAD sample. As depicted in Fig. 3K, cluster B was characterized by elevated stromal, immune, and ESTIMATE scores. A more granular examination through ssGSEA revealed marked disparities in immune cell infiltration between subtypes A and B (Fig. 3L). Specifically, subtype B demonstrated notably higher infiltration levels of activated B cells, activated CD4+ T cells, activated CD8+ T cells, activated dendritic cells, macrophages, natural killer T cells, natural killer cells, T follicular helper cells, Type 1 T-helper cells, and Type 2 T-helper cells compared to subtype A.

From these observations, it can be inferred that subtype B was significantly enriched in immune activation pathways, which appears to be intricately linked with the dynamics of TIME. In the quest to decode the biological intricacies of TLS subtypes, we scrutinized 189 DEGs between subtypes A and B. These were visually represented through a volcano plot (Figs. 4A, |logFC| > 1, p < 0.05). To identify pertinent biological pathways, we engaged in comprehensive GO and KEGG enrichment analyses. The GO and KEGG enrichment analysis pinpointed the DEGs’ significant involvement in pathways such as cytokine-cytokine receptor interaction, upregulation of cytokine production, lymphocyte-mediated immunity, and adaptive immune responses founded on somatic leukocyte proliferation (Figs. 4B, 4C). A cnetplot was constructed to visualize the interplay within these pathways, highlighting the top five pathways (Fig. 4D).

Identification of TLS gene subtypes

In the continuum of our research, a univariate Cox regression analysis was conducted to evaluate the prognostic significance of the 189 identified DEGs. This rigorous analysis led to the isolation of 16 key signature genes closely associated with overall survival (OS) in COADREAD patients (p < 0.001). These genes, namely HOXC6, SFRP2, CYP1B1, GAS1, SPOCK1, SPP1, POSTN, MARCO, VSIG4, THBS2, OLR1, COL10A1, FAP, GZMB, C5AR1, and APOE, were depicted in a forest plot for clarity (Fig. S2A).

Subsequently, COADREAD patients were stratified into two distinct groups based on the expression profiles of these 16 prognostic genes, utilizing a consensus clustering algorithm. The partitioning suggested that k = 2 was an optimal division, forming two gene clusters, namely A and B (Fig. 4E). Survival analysis revealed a noteworthy finding: patients in subtype B exhibited a significantly enhanced survival probability compared to those in subtype A (log-rank test, p < 0.001; Fig. 4F). The box-map analysis highlighted that most of the signature genes, including APOE, VSIG4, C5AR1, CYP1B1, and SPP1, were predominantly 11 upregulated in cluster B (Fig. 4G, p < 0.05). Moreover, we extended our analysis to explore the correlation of these gene expressions with various clinical features such as age, stage, gender, recurrence, metastasis, fustat, and futime, stratifying them across the molecular subtypes (Fig. 4H).

Construction of TLS score

In advancing our study, a TLS score was formulated based on GSVA of the 16 identified signature genes. This scoring system facilitated the stratification of patients into two distinct genomic subtypes: those with low and high TLS scores. Survival analysis revealed a critical insight; patients categorized in the low-score group exhibited a markedly higher survival probability compared to their high-score counterparts (Fig. 4I, p < 0.001). A Sankey diagram was employed to eloquently depict the patient distribution across the two TLS molecular subtypes, gene subtypes, and score groups, alongside their respective prognostic implications (Fig. 4J).

Association of TLS score with immune infiltration

Spearman correlation analysis established a positive correlation between the TLS score and the infiltration of 23 different immune cells. This correlation suggested an augmented immune component presence in the TIME of patients within the high-score cluster, potentially indicating a more favorable immune prognosis (Fig. 4K, p < 0.05). Building upon these findings, we delved deeper into the association of TLS with cytokine-chemokine interaction networks. Notably, the high-score TLS cluster demonstrated significant enrichment in chemokines and their receptors, interleukins, and interferons along with their respective receptors (Fig. 4L). The hallmark pathway analysis revealed a significant positive correlation 12 of the TLS score with pathways indicative of immune activation. These pathways included interferon gamma response, interferon alpha response, inflammatory response, and IL-2 STAT5 signaling. Conversely, an inverse correlation was observed between the TLS score and pathways typically associated with malignancy, such as peroxisome proliferation, MYC target activation, oxidative phosphorylation, fatty acid metabolism, and DNA repair mechanisms (Fig. 4M). These findings underscored the strong association of the TLS score with predominantly immune-activated pathways.

Association of TLS score with MSI, TMB, CSC, and somatic mutations

Cancer stem cells (CSCs), which were acknowledged for their significant roles due to their selfrenewal capacity and differentiation potential (Nassar & Blanpain, 2016). We delved into the relationship between the TLS score and CSC indices, including mDNAsi and mRNAsi. Our analysis, as depicted in Fig. 5A, revealed a notable negative linear correlation between the TLS score and the mDNAsi index (R = −0.18, p = 6e−04). This trend was further demonstrated by the lower mDNAsi observed in the high TLS score cluster compared to the low TLS score cluster (p < 0.01; Fig. 5B). Similarly, Fig. 5D presented a negative linear correlation between the TLS score and the mRNAsi index (R = −0.63, p = 0), with the high TLS score cluster exhibiting significantly lower mRNAsi than its low-score counterpart (p < 0.0001; Fig. 5E). Heat maps further corroborate these findings, illustrating the inverse relationship between both mDNAsi and mRNAsi index values and the TLS score (Figs. 5C, 5F). These insights suggested that COADREAD patients with a higher TLS score may exhibit reduced stem cell characteristics and diminished cell differentiation capabilities.

In the context of microsatellite instability (MSI), a promising genomic biomarker for gauging 13 patient responsiveness to immunotherapy (Vilar & Gruber, 2010), our findings revealed a significantly higher MSI in the high TLS score cluster compared to the low (p < 0.05; Fig. 5H). Spearman correlation analysis further substantiates this, indicating a positive association between MSI and TLS score (R = 0.18, p = 7e−04; Fig. 5G). The positive correlation between MSI expression and TLS score was also visually represented in a heat map (Fig. 5I). This pattern suggested a potential link with enhanced immunotherapy efficacy.

Regarding tumor mutational burden (TMB), which served as an indicator of cancer mutation volume and has clinical relevance with ICIs outcomes (Chan et al., 2019), our data analysis from TCGACOADREAD cohorts indicated a higher TMB in the high TLS score cluster than in the low score cluster (p < 0.001; Fig. 5K). Spearman’s analysis corroborated a positive relationship between TMB and TLS score (R = 0.2, p = 1e−04; Fig. 5J), with a heat map further illustrating this correlation (Fig. 5L). These findings were particularly noteworthy for patients with higher TMB who typically show more benefit from ICIs.

Our analysis extended to the somatic mutation landscape within the two distinct TLS score clusters. In the high TLS score cluster, the predominant mutated genes were TP53, TTN, APC, KRAS, SYNE1, MUC16, PIK3CA, OBSCN, FAT4, and RYR2. Conversely, the low TLS score cluster was characterized by mutations primarily in APC, TP53, KRAS, TTN, SYNE1, PIK3CA, FAT4, MUC16, OBSCN, and MUC4. Notably, the mutation frequency across these genes was more pronounced in the high TLS score cluster, as illustrated in Figs. 5M and 5N. Additionally, a forest plot detailed the top 12 genes exhibiting the most significant differences in mutation frequency between the high and low score clusters (Figs. 5O). This led to the inference that gene mutation frequency in the high TLS score cluster was not only higher but also more14 diverse. Detailed SNV, CNV, and methylation profiles of prognostic TLS-DEGs in COADREAD patients were depicted in Figs. S2 and S3. Analysis revealed significant disparities in recurrence/metastasis, stage, and fustat between the low- and high-TLS score clusters (Figs. 5P–5S).

Interplay of TLS score and immune checkpoints expression

Subsequent research focused on the expression levels of immune checkpoints in relation to TLS scores. Elevated expression of key immune checkpoints, including PDCD1 (PD-1), CD274 (PD-L1), CTLA4, TIGIT, and LAG3, was predominantly observed in the high-score cluster, as depicted in Fig. 5T. This trend suggested that TLS score was potentially served as a predictive marker for immunotherapy efficacy.

Evaluating TLS score as a predictor for PD-L1 blockade immunotherapy

In the realm of T-cell immunotherapy, a promising approach for cancer treatment, we explored the significance of the TLS score in COADREAD, building on prior findings and our current results (Waldman, Fritz & Lenardo, 2020). Within the GSE135222 cohort, a notable observation was the higher progression rate (100%) in patients with lower TLS scores compared to those with higher scores (75%) (Fig. 5V). Furthermore, a markedly prolonged progression-free survival (PFS) was evident in patients with higher TLS scores (Fig. 5U, p = 0.0061). Analyzing the GSE176307 cohort, which comprised 90 patients treated with anti-PD-L1 receptor blockers, revealed a spectrum of responses ranging from complete response (CR) and partial response (PR) to stable disease (SD) and progressive disease (PD). Significantly, patients exhibiting CR/PR had elevated TLS scores compared to those with SD/PD (Fig. 5X). In this cohort, patients in the high-TLS cluster demonstrated substantial clinical benefits and notably longer overall survival (OS) than 15 their low-TLS counterparts (Fig. 5W, p = 0.018).

Analysis of drug susceptibility in relation to TLS score

Subsequent to establishing the TLS score as an effective predictor for immunotherapy response, our study delved into the variances in IC50 values for various targeted anticancer drugs across low- and high-TLS clusters. Intriguingly, patients within the high-TLS cluster exhibited notably lower IC50 values for a range of anticancer agents, namely A.770041, AP.24534, AG.014699, AICAR, ABT.263, and AMG.706 (Fig. 5Y). These findings underscored the potential utility of the TLS score as a guiding metric in selecting appropriate anticancer drugs.

Impact of TLSRGs on colorectal cancer cell behavior in vitro

In our investigation, RT-q PCR was utilized to ascertain the expression levels of C5AR1, APOE, CYP1B1, and SPP1 in a standard colon cell line and three colorectal cancer cell lines (Figs. 6A–6D). These four genes were markedly upregulated in cancer cells compared to normal cells. We employed siRNA techniques targeting C5AR1, APOE, CYP1B1, and SPP1. This was conducted in HCT116 and HT29 cell lines, selecting siRNA-1 and siRNA-2 for further analysis based on their transfection efficiency exceeding 70% (Figs. S4A–S4H). A CCK8 assay was subsequently performed to evaluate cell proliferation. Results demonstrated that the knockdown of APOE, C5AR1, CYP1B1, and SPP1 significantly impeded the proliferation of both HCT116 and HT29 cells (Figs. 6E–6L). The Transwell assay was conducted to assess cellular invasion capabilities. Notably, the suppression of APOE, C5AR1, CYP1B1, and SPP1 led to a marked reduction in the invasion capacity of HCT116 and HT29 cells (Figs. 6M and 6N). Furthermore, an unexpected observation was the considerable decrease in PD-L1 expression in both HCT116 and HT29 cells following the knockdown of these genes. This reduction in PD-L1 expression could hint at the potential of combining TLS inducers with PD-L1 inhibitors as a therapeutic strategy (Figs. 6P and 6Q). Additionally, the immunohistochemical data from the HPA database revealed an elevated protein expression of APOE, C5AR1, CYP1B1, and SPP1 in the stroma of COADREAD tissue (Fig. 6O).