The data, including gene sequencing data and corresponding clinical data, were obtained from the TCGA and GEO databases. TCGA-CRC dataset functioned as a discovery set for marker selection and model construction, consisting of 558 samples. Samples with survival times under 30 days were removed. GEO datasets GSE38832 and GSE17537 served as external validation sets to assess the robustness of the model. IMvigor210 immunotherapy cohort was obtained from the original article (18). Clinical information for all samples were detailed in Supplementary Table 1 .

The mutation landscapes of CRC patients were analyzed by the “maftool” R package (19). Unsupervised clustering, via the “ConsensusClusterPlus” R package, segmented patients into two subclusters (20). This segmentation was guided by maximizing intergroup variances and minimizing intragroup variances. The immune infiltration environment in each patient was assessed using CIBERSORT and xCell (21, 22). XCell was also utilized to compare the environment scores in different groups.

Using the “DESeq2” package, we obtained differentially expressed genes (DEGs) comparing two clusters, tumor and normal. The Lasso-Cox regression was applied to filter out DEGs that are associated with OS, reduce dimensions, and construct the prognostic model. Patients were classified into two risk groups by the risk score, which was calculated as: R i s k   S c o r e = ∑ i = 1 14 C o e f i ∗ E x p r i Coef_i denotes the estimated regression coefficient, while Expr_i denotes the gene’s expression level. The model’s predictive efficiency and accuracy were evaluated via the “survival”, “survminer”, “timeROC”, and “pheatmap” R packages. MSI status was obtained from the “cBioPoralData” R package. According to the results of multivariate Cox regression, we plotted the forest plot and nomogram.

The 14 model genes were ranked based on the feature importance calculated through “randomForest” package. All feature importance was normalized to range from 0 to 1. Genes with importance exceeding 0.75 were classified as TLS-related hub genes. To comprehensively analyze the biological pathways of TLS-related genes, a dual approach involving gene set enrichment analysis (GSEA) and gene set variation analysis (GSVA) were used. Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), and hallmark annotation were conducted with the “clusterProfiler” and “GSVA” packages (23, 24). The visualization of results using the “ggplot2”.

Human CRC tissue samples were collected from QuXian People’s Hospital, including 20 TLS+ tissue samples and 10 TLS- tissue samples (initial diagnosis of CRC, denial of other tumors, unaccepting of neoadjuvant therapy and complete clinical data). Clinical information for all samples was detailed in Supplementary Table 3 . The research received approval from the Ethics Committee of the QuXian People’s Hospital. Two pathologists independently evaluated the TLS existence in all CRC samples. Questionable results were then co-reviewed for an agreed annotation. TLS were classified into three stages based on maturity: (1) early TLS (compact clumps of lymphocytes without segregated T and B cell zones), (2) primary follicle-like TLS [(have follicular dendritic cell (FDC), lack of GC)] and (3) secondary follicle-like TLS (have an active GC) (25). Samples with TLS were categorized as TLS+ group, while others were deemed TLS- group.

Firstly, slides were deparaffinized and placed in the citrate-phosphate buffer (heat them in a microwave at low heat for 15 min). Circled the tissue with an immunohistochemical pen to block endogenous peroxidase. Secondly, incubated the slides in hydrogen peroxide water at room temperature and avoid light for 25 min. Added BSA buffer dropwise onto the slides and sealed for 30 min. After removing the BSA buffer, we added the prepared primary antibody (CD3: ab16669, Abcam, Britain, 1:1000) dropwise and incubated the slides overnight at 4 °C in a wet box. Next, the slides were washed, and horseradish peroxidase (HRP) labelled secondary antibody (Goat anti-mouse/rabbit IgG H&L, RCB054, RecordBio, China) were dripped onto the slides. Tyramide signal amplification (TSA) buffer was used to visualize each biomarker. Dripped corresponding TSA dye onto the slides and incubated them at room temperature in the dark for 10 min to stain the antibody. Repeated the previous steps using different primary antibodies (CD20: ab78237, Abcam, 1:1000; CD21: ab315160, Abcam, 1:1000) and switched different TSA dyes. Finally, stained nucleus with DAPI and coverslipped anti-fluorescence quenching to preserve fluorescence. The slides were digitized using OLYMPUS (OlyVIA, Japan).

Next, we performed IHC to examine the PRRX1 expression level in TLS+/- CRC formalin-fixed paraffin-embedded (FFPE) samples, using PRRX1 antibody (ab211292, Abcam, 1:100). Tissues were scored according to the product of stained area percentage (0 = 0%, 1 = 1%–25%, 2 = 26%–50%, 3= 51%–75%, and 4 = more than 75%) and the staining intensity scale [ranged from 0 (no staining) to 3 (strong staining), measured by Halo]. The IHC score for all samples was detailed in Supplementary Table 3 .

The data analysis was performed using R software version 4.2.2. For categorical data, the Chi-square test and Fisher’s exact test were used to assess differences between groups. For continuous data, if the data are normally distributed and have equal variances, an independent samples T-test was used; if the data do not follow a normal distribution, the Mann-Whitney U test was used. Spearman’s correlation analysis was used to identify the relationship between two variables. The threshold for statistical significance was set at P< 0.05 (two-tailed).

By integrating current studies and reviews, we identified 42 TRGs, of which 25 were expressed in TCGA-CRC dataset ( Supplementary Table 2 ). Among 517 CRC patients, 94.39% (488) experienced somatic mutations. APC and TP53 had the highest mutation frequencies ( Figure 1A ). We further analyzed the somatic mutation pattern of TRGs ( Figure 1B ). 11.41% (59) of patients experienced somatic mutations, predominantly missense mutations.IL1R1 and IL1R2 exhibited the highest mutation frequencies, while CD4 and LAMP3 primarily underwent frame shift deletions. Using STRING (V12.0), a protein-protein interaction (PPI) network was constructed ( Figure 1C ). GISTIC2.0 was used to calculate the copy number variations (CNV). The results showed that copy number amplification (including homozygous and heterozygous mutations) was the major mutation ( Figure 1D ).

Unsupervised consensus clustering based on TRGs categorized CRC patients into two molecular subclusters. The consensus clustering results indicated that the optimal clustering number was 2 ( Supplementary Figures 1A–C ). Principal component analysis (PCA) proved the optimal division into two clusters, C1 and C2 ( Supplementary Figure 1D ). The Kaplan-Meier analysis demonstrated a statistically difference in C1 and C2 (P<0.01, Figure 2A ). The prognosis for C1 was more favorable compared to C2. Meanwhile, we examined the clinicopathological characteristics of two subclusters. C1 had a higher proportion of males, N0, M0, and stage I/II ( Supplementary Figures 1D–I ).

Next, we analyzed the TRGs expression profiles of two subclusters in the TCGA-CRC dataset ( Figure 2B ). We found C1 had higher expression levels of CCL18, CCL21, CXCL10 and CXCL11. C2 was characterized by the remarkable enrichment in SDC1 and CCL20. In general, TRGs were more highly expressed in C1 than C2. Immune infiltrating cells are a crucial component of TME. We explored the differences in TME infiltration between two subclusters. ( Figures 2C, D ). B cells, CD8+ T cells, follicular helper T (T_FH) cells, activated dendritic cells (DC), neutrophils and granulocyte-macrophage progenitor (GMP) cells were higher in C1, while naive CD4+ T cells and natural killer T (NKT) cells were higher in C2. Both macrophage M1 and macrophage M2 had higher infiltration levels in C1. Besides, we further noticed that C1 correlated positively with higher immune and microenvironment scores, while there was no difference between two subclusters in stroma scores ( Figure 2E ).

Volcano map was used to demonstrate DEGs between C1 and C2, with 639 genes upregulated in C1 and 86 in C2 ( Supplementary Figure 2A ). Firstly, we performed over-representation analysis (ORA) to determine the biological behavior behind these DEGs. The GO enrichment analysis revealed that immune receptor activity enrichment in molecular function (MF); plasma membrane signaling receptor complex and immunological synapse enrichment in cellular component (CC); T cell proliferation and activation enrichment in biological process (BP) ( Figure 3A ). The KEGG analysis indicated predominantly enrichment in cytokine-cytokine receptor interaction, calcium signaling pathway, and cell adhesion molecules ( Figure 3B ). GSEA showed that C1 was enriched in immune response, activation, and antigen presentation, including immunoglobulin mediated immune response, Th1 and Th2 cells differentiation, and Toll-like receptor signaling pathway ( Figures 3C, D ). Metabolic processes and biosynthesis related pathways were enriched in C2, such as glucuronate metabolic process, and biosynthesis of amino acids and steroid hormone ( Figures 3E, F ). GSVA showed differentially active signaling pathways and immune responses between two clusters ( Figure 3G ). Taken together, this suggested that TLS-enriched subcluster, C1, may inhibit tumor progression by activating non-specific anti-tumor immune responses. While TLS-depleted subcluster, C2, was for hypermetabolic and immunosuppression state.

To further analyze and quantify the TLS characteristics, a TLS-related prognostic model was established. TCGA-CRC dataset served as the training set, while GSE38832 and GSE17537 were utilized for the external validation. DEGs between tumor and normal samples were measured with the “DESeq2” ( Supplementary Figure 2B ). There were 508 intersection genes left ( Supplementary Figure 2C , Supplementary Table 2 ). After univariable Cox regression analysis, we found 43 OS-related DEGs. Ultimately, 14 genes were selected for the risk model: CAVIN2, CCL19/21, CD8A, CXCL5, FHL1, IGHG1, MMP1, MRC1, NEXN, MOTUM, PRRX1, SELL and VSIG4 ( Supplementary Figures 2D, E ). The model genes displayed a distinct expression pattern between two groups ( Supplementary Figure 2F ).

Based on the optimum cutoff value, patients were divided into two risk groups. Kaplan-Meier curves indicated lower mortality and better prognosis in the low-risk group (P<0.001, Figures 4A, B ). Furthermore, the AUC revealed the 1-, 3-, and 5-year survival rates were all above 0.7, which proved the accuracy of this model was wonderful ( Figure 4C ). External validation datasets GSE38832 and GSE17537 also demonstrated the model’s excellent prognostic performance ( Figures 4D–I ).

Different kinds of immune cell infiltration were estimated using CIBERSORT ( Figure 5A ). Plasma cells, macrophages and follicular helper T cells were higher in high-risk group. Additionally, we estimated the proportion of each cell in TME ( Figure 5B ). There were remarkably positive associations between risk score and CD8+ T cells, CD4+ naive T cells, macrophages and activated mast cells (P<0.01, Figure 5C ). Further analysis revealed that there were more patients with KRAS mutations and fewer patients with BRAF mutations in the high-risk group (no statistically significant, Figures 5D, E ). Next, we estimated the association between immune checkpoint genes and 14 model genes ( Figure 5F ). We found that the vast majority of model genes were positively correlated with immune related genes, except NOTUM.

Our above findings demonstrated a potential correlation between risk score and immunotherapy response. The low-risk group had higher stroma, immune, and TIME scores ( Supplementary Figures 3A–C ). TIDE score has been extensively employed to assess the likelihood of immune evasion and resistance to immunotherapy. We discovered the high-risk group had a much higher TIDE score ( Supplementary Figure 3D ). A higher TIDE score was linked to worse ICB effectiveness as well as shorter OS with ICB therapies. The IMvigor210 dataset, involving patients treated with atezolizumab, was utilized to validate these observations. In IMvigor210, low-risk patients had a greater prognosis and higher disease control rate (p<0.05, Supplementary Figures 3E, F ).

MSI status was used to be an effective predictor for patients receiving ICB treatment. Typically, patients with high level of MSI were classified as MSI-High, indicating genomic instability resulting from repeated microsatellite expansions or contractions. Conversely, patients with low or no MSI were classified as microsatellite instability-low (MSI-L)/Microsatellite Stability (MSS). We analyzed the relationship between patients with different MSI status and risk score, and found no significant correlation ( Figure 6A ). Then, we combined the risk score with MSI status for stratified survival analysis ( Figure 6B ). Investigations indicated that MSI-H and low-risk scores patients had the best prognosis, followed by those with MSI-L and low-risk scores. Univariate and multivariate Cox regression confirmed the risk score as an independent prognostic factor for CRC (HR=3.706, 95%CI=2.552–5.380, Figures 6C, D ). Exploiting these Clinicopathological characteristics, a new-type nomogram was constructed to assess the clinical prognostic more accurately ( Figure 6E ).

Random forest analysis was used to determine the feature importance of 14 model genes, among which three genes—VSIG4, SELL and PRRX1—stood out with importance scores exceeding 0.75 ( Supplementary Figure 4A ). In the high-risk group, the expression level of VSIG4, SELL and PRRX1 were higher ( Supplementary Figure 4B ). Correlation analysis revealed that the vast majority of model genes were significantly positively correlated with hub genes ( Supplementary Figure 4C ). We conducted a further study to explore the correlation between the expression of TLS-related hub genes and the infiltration of immune cells in TME ( Supplementary Figure 4D ). VSIG4, SELL and PRRX1 exhibited a positive relationship with macrophages M2 and neutrophils.

GSEA and GSVA were conducted to investigate the molecular mechanisms and immune functions of the three hub genes. Results indicated that VSIG4 was positively correlated with signaling pathways such as NF-kappa B signaling pathways, PI3K-AKT signaling pathways, and cytokine-cytokine receptor interaction ( Figures 7A, B ). High expression of SELL primarily enriched in the intestinal immune network for IGA production, oxidative phosphorylation, arginine and proline metabolism, Fc epsilon signaling pathway, and others ( Figures 7C, D ). GSEA of PRRX1 subsequently uncovered that PRRX1 was significantly linked to epithelial-mesenchymal transition (ECM) -receptor interaction, proteasome and TGF-beta signaling pathway ( Figure 7E ). High expression of PRRX1 was mainly enriched in riboflavin metabolism, regulation of actin cytoskeleton, and others ( Figure 7F ). The low expression of PRRX1 was enriched in pyruvate metabolism, purine metabolism, peroxisome, etc. In addition, PRRX1 may serve as a prominent regulatory factor in tumorigenesis and progression.

Increasing evidence has shown that the mutual regulatory patterns between mRNA, lncRNA, miRNA, and downstream target genes are intimately related to tumors. Our study exploited the Human MicroRNA Disease Database (HMDD) and miRNA Walkthrough (miRWalk) databases to identify the miRNA and lncRNA associated with CRC. Excluding interactions unrelated to disease, only two mRNA (SELL, PRRX1) and eight miRNA were retained ( Figure 8A ). Secondly, those mRNA-miRNA interactions were validated in the ENCORI database, hsa-miR-485–5p and hsa-miR-20a-5p regulated PRRX1 were discovered. Utilizing Cytoscape, we constructed a ceRNA network, including 1 mRNA, 2 miRNA and 127 lncRNA (56 for hsa-miR-20a-5p, 71 for hsa-miR-485–5p) interaction pairs ( Figure 8B ).

Using IHC to assess the PRRX1 expression in the TLS+/- groups. Firstly, we observed TLS images through hematoxylin-eosin (H&E) and mIF ( Figure 9A ). The lymphocyte aggregates containing CD3+T cells, CD20+B cells, and CD21+FDCs that form corresponding regions are considered secondary follicle-like TLS, also known as mature TLS. To further validate the expression of PRRX1, we performed IHC on TLS+/- CRC tissues ( Figures 9B, C ). PRRX1 was significantly higher in the TLS- group ( Figure 9D ). Through these above findings, we speculate that PRRX1 may be a negative regulator of TLS and play a crucial role in tumor metastasis and immunocompetence.