Ethics approval was obtained from the institutional review boards of NanFang Hospital and Fujian Provincial Cancer Hospital (NFEC-2023-221). The requirement for informed consent was waived for this study. The study was conducted following the guidelines of the Declaration of Helsinki and the Transparent reporting of a multivariable prediction model for individual prognosis or diagnosis (TRIPOD) statement criteria.

The flow chart of patient recruitment in this study is shown in Supplementary Figure 1. The inclusion criteria were patients who underwent radical surgery with pathologically diagnosed stage I-III CRC, available follow-up data and clinicopathological characteristics, and hematoxylin and eosin (HE) slides with invasive tumor components. The exclusion criteria were patients with unavailable formalin-fixed paraffin-embedded (FFPE) specimens, a history of cancer, or received neoadjuvant treatment. As a result, a total of 327 consecutive patients were included in the training cohort between January 2011 and December 2013 from Nanfang Hospital. An independent validation cohort contained 327 consecutive patients from Fujian Provincial Cancer Hospital between October 2011 and December 2013. Two independent pathologists reassessed all samples based on the 8th edition AJCC staging criteria.

Clinicopathological characteristics included age, sex, primary tumor location, preoperative carcinoembryonic antigen (CEA) level, preoperative carbohydrate antigen 199 (CA199) level, tumor differentiation, tumor size, pT stage, and pN stage. Adjuvant chemotherapy after radical surgery is recommended for patients with high-risk stage II and stage III CRC according to NCCN guidelines.

A standardized follow-up protocol was implemented, including a serum CEA test every 3 months after surgery and every 6 months after 3 years; CT examination from chest to pelvis every 6 months in the first 5 years after surgery; and colonoscopy at 1 year after surgery.

FFPE samples were cut into 4-µm sections and stained with antibodies against CD3 and CD8 (Maixin Biotech. Co., Ltd., Fuzhou, China). Immunohistochemical staining was performed as previously described (36, 37). Whole slide images of stained slices were digitized by Aperio ImageScope (Leica Biosystems, CA, USA) at 20× magnification as.svs format files.

The Immunoscore was assessed in the following steps (Figure 1). First, two pathologists who were blinded to the prognostic information selected five representative regions at the TC and five representative regions at the IM. Second, CD3+ and CD8+ stained immune cells were quantified using QuPath software (version 0.2.3). Third, CD3+ and CD8+ density was used to divide the individual cases into “high” or “low” immune groups, and patients with a mean density ≥ 75th percentile were considered a “high” immune group. A high immune group score was set as 1, and a low immune group score was set as 0. The CD3_TC, CD3_IM, CD8_TC, and CD8_IM scores were added and converted into an Immunoscore (I0 - I4). Finally, patients were divided into two groups based on their Immunoscore: I0–I1 was classified as low Immunoscore (Lo IS), and I2-I4 was classified as intermediate-high Immunoscore (Int-Hi IS).

The regions at the TC and IM, which were used to calculate the density of CD3+ and CD8+, were used for multiphoton imaging. Image acquisition for multiphoton imaging was performed with a 200× original magnification objective on another unstained serial section and then compared with the HE image for histologic assessment (27). More information about the multiphoton imaging system can be found in the Supplementary Methods.

The framework we constructed for the quantitative extraction of collagen features is shown in the Supplementary Methods. In summary, 142 collagen features (Supplementary Table 1), including morphological features, histogram-based features, gray level concurrence matrix (GLCM) features, and Gabor wavelet transform features, were achieved automatically via MATLAB 2016b (Mathworks, Natick, MA, USA) (30–32). Finally, a total of 284 collagen features were obtained, including 142 from TC and 142 from IM, for further statistical analyses.

LASSO regression, which is suitable for the regression of high-dimensional data, was used to select the most useful predictive features (33–35). The LASSO regression used an L1 penalty to shrink the coefficients to zero. The penalty parameter λ, also called the tuning constant, controls the number of collagen features to enter the model. In this study, we applied 10-fold cross-validations to select the optimal value of λ via 1-standard error (SE) criteria in the training cohort, and the collage signature was calculated for each patient via a linear combination of selected features that were weighted by their respective coefficients in the training cohort. Then, the collage signature in the validation cohort was calculated by the selected features with their respective coefficients obtained from the training cohort. Details of the LASSO regression are provided in the Supplementary Methods.

The collagen signature and clinicopathologic characteristics were included in univariate analysis to investigate their association with Lo IS, and variables with P < 0.10 were included in multivariable analysis. A backward stepwise selection method with Akaike’s information criterion as the stopping rule was used to select the independent predictors of Lo IS (38). To facilitate clinical application, we developed a collagen nomogram according to the independent predictors in the training cohort (39).

The Hosmer−Lemeshow test was applied to estimate the goodness of fit of the model (40). The multicollinearity of the collagen nomogram was evaluated through the variance inflation factor (VIF) (41). The area under the curve (AUC) and the calibration curve were applied to assess the discrimination and calibration of the collagen nomogram. Then, the collagen nomogram was performed in the validation cohort, and its AUC and calibration curve were acquired. More information on the nomogram is shown in the Supplementary Methods.

To assess the clinical application value of the collagen nomogram, a traditional model was developed for comparison with the collagen nomogram. In our study, the traditional model was constructed based on clinicopathological predictors after univariate and multivariable logistic regression in the training cohort. The clinical usefulness of the collagen nomogram was evaluated by decision curve analysis (DCA) and clinical impact curves (CICs) (42). The maximum Youden index value of the ROC curve of the two models was measured to estimate the specificity, sensitivity, accuracy, positive predictive value (PPV), and negative predictive value (NPV). Moreover, the net reclassification improvement (NRI) and integrated discrimination improvement (IDI) were used to show the improvement of the collagen nomogram compared with the traditional model (43, 44). Details of the NRI and IDI are provided in the Supplementary Methods.

Baseline characteristics were compared between the training and validation cohorts by t test, U test, Fisher’s exact test, and χ2 test when applicable. The odds ratio (OR) and 95% confidence interval (CI) of the predictors were calculated using multivariable logistic regression. Survival curves were generated by using the Kaplan–Meier method and compared by log-rank tests. Univariate and multivariable analyses with Cox proportional hazards regression determined the hazard ratio (HR) of predictors for disease-free survival (DFS) and overall survival (OS). All statistical analyses were performed with SPSS version 22.0 software (IBM, Armonk, New York USA) and R version 4.0.3 (http://www.r-project.org/). All P values were two-sided, and statistical significance was defined as P < 0.05.

The baseline characteristics of the patients in the training and validation cohorts are summarized in Table 1. A total of 421 (64.3%) patients were < 65 years old, with 405 (61.9%) men. The clinicopathological characteristics of the two cohorts were similar (Supplementary Table 2).

The density of CD3+ and CD8+ TILs in the TC and IM is shown in Supplementary Figure 2, with a higher density of TILs in the IM than in the TC for both CD3+ and CD8+ cells. The cutoff values of CD3+ and CD8+ cells were 593 and 382 cells/mm² in the TC and 1382 and 714 cells/mm² in the IM, respectively (Supplementary Table 3). Finally, the proportions of patients with Lo IS and Int-Hi IS were 34.3% and 65.7% in the training cohort and 35.8% and 64.2% in the validation cohorts, respectively.

The median follow-up duration [interquartile range (IQR)] was 72 (42–85) months in the training cohort and 71 (40–83) months in the validation cohort. The 5-year DFS and OS rates were 67.6% and 74.6%, respectively, in the training cohort. Similarly, the DFS and OS rates were 67.0% and 74.3%, respectively, in the validation cohort (Supplementary Figure 3). Patients with Int-Hi IS from the training cohort had a significantly better 5-year DFS (76.7% vs. 49.1%; P < 0.001) and 5-year OS (84.2% vs. 55.5%; P < 0.001) than patients with Lo IS (Figure 2A). Likewise, patients with Int-Hi and Lo IS had significant differences in 5-year DFS (80.5% vs. 42.7%; P < 0.001) and 5-year OS (83.8% vs. 56.4%; P < 0.001) in the validation cohort (Figure 2B).

The framework of the collagen signature is presented in Figure 3. As a result, a collagen signature was constructed based on sixteen collagen predictors from 284 collagen features by LASSO regression. (Supplementary Figure 4). The calculation formula for the collagen signature is proposed in the Supplementary Results. The distributions of the 16 collagen predictors and Immunoscore for each patient in the training and validation cohorts are shown in Supplementary Figure 5. The patients with a high collagen signature were more likely to show Lo IS in both cohorts (Figure 4). The collagen signature yielded an AUC of 0.896 (95% CI, 0.854-0.936) in the training cohort and 0.903 (95% CI, 0.863-0.944) in the validation cohort. A significant association between the collagen signature and Lo IS was found when stratified analysis was performed (Supplementary Tables 4, 5).

We also assessed the performance of the collagen signature and the single selected collagen feature in predicting Immunoscore. The results indicated that the collagen signature was more powerful than any individual parameter, demonstrating the added predictive value of the collagen signature (Figure 5).

Univariate and multivariable logistic regression was performed to identify independent predictors of Lo IS. The results showed that the collagen signature (OR: 4.632, 95% CI: 3.068-6.993; P < 0.001), tumor differentiation (OR: 2.537, 95% CI: 1.121-5.741; P = 0.026), pT stage (OR: 2.602, 95% CI: 1.106-6.121; P = 0.028), and pN stage (OR: 2.550, 95% CI: 1.197-5.433; P = 0.015) were independent predictors of Lo IS (Table 2). Then, the collagen nomogram was developed, integrating the above four predictors (Figure 6A). ROC curve analysis indicated that the collagen signature had the most discrimination ability compared with the other predictors (Supplementary Figure 6). Alluvial diagrams were employed to intuitively illustrate the relationship between the four predictors and Immunoscore (Supplementary Figure 7). The variance inflation factor (VIF) values of each predictor were < 10, indicating that there was no multicollinearity among the four predictors (Supplementary Table 6). The Hosmer−Lemeshow test yielded a nonsignificant statistic (P = 0.299), demonstrating that there was no departure from a perfect fit.

In the training cohort, the collagen nomogram yielded satisfactory discrimination with an AUC of 0.925 (95% CI: 0.895-0.956). The calibration curve showed good agreement between the predicted and the actual Lo IS probability (Figure 6B). Similar results were observed in the validation cohort (AUC: 0.911, 95% CI: 0.872-0.949) and all patients (AUC: 0.918, 95% CI: 0.893-0.942) (Figure 6C).

A traditional model was developed based on tumor differentiation, pT stage, and pN stage in the training cohort (Supplementary Table S7). The traditional model yielded AUCs of 0.683 (95% CI, 0.622-0.744) in the training cohort, 0.680 (95% CI, 0.619-0.742) in the validation cohort, and 0.685 (95% CI, 0.642-0.728) in all patients. The collagen nomogram exhibited better discrimination ability than the traditional model (training cohort: 0.925 vs. 0.683; validation cohort: 0.911 vs. 0.680; all patients: 0.918 vs. 0.685; all P < 0.001) (Figure 6B). Moreover, the stratified analysis showed that the collagen nomogram was still superior to the traditional model among the subgroups in the training cohort, the validation cohort, and all patients (Supplementary Figures S8–S10). DCA revealed that the collagen nomogram could add more benefits than the traditional model (Figure 7A). CICs were generated to intuitively recognize the application value of the collagen nomogram to more accurately identify patients with Lo IS (Figure 7B).

Furthermore, the collagen nomogram exhibited better sensitivity (97.6% vs. 82.1%), specificity (87.3% vs. 43.7%), accuracy (89.9% vs. 56.6%), PPV (72.3% vs. 43.2%), and NPV (99.1% vs. 82.5%) in the training cohort. Similar results were observed in the validation cohort and all patients (Table 3). The corresponding NRI and IDI both showed significantly improved classification accuracy of the collagen nomogram compared with the traditional model in the training cohort, validation cohort and all patients (Table 4).

Patients were divided into high- and low-probability Lo IS groups based on the ROC curve of the collagen nomogram. We found that patients with a low-probability Lo IS subgroup showed a better prognosis than patients with a high-probability Lo IS subgroup in the training cohort (Supplementary Figure S11A), the validation cohort (Supplementary Figure S11B) and all patients (Supplementary Figure S11C). This result was also observed in stage I-II (Supplementary Figure 12) and III patients (Supplementary Figure S13). Cox regression analysis demonstrated that the probability of Lo IS was an independent prognostic factor after adjusting for other variables in the training cohort [DFS: HR 2.475 (95% CI, 1.667-3.675), P < 0.001; OS: HR 2.179 (95% CI: 1.409-3.370), P < 0.001] (Supplementary Table 8), the validation cohort [DFS: HR 2.211 (95% CI, 1.510-3.239), P < 0.001; OS: HR 2.111 (95% CI: 1.366-3.262), P < 0.001] (Supplementary Table 9), and all patients [DFS: HR 2.350 (95% CI, 1.787-3.091), P < 0.001; OS: HR 2.119 (95% CI: 1.559-2.881), P < 0.001] (Supplementary Table 10). The collagen signature and clinicopathological predictors with the corresponding DFS and OS status are presented in Figure 8.

In addition, we investigated the chemotherapy benefits of high-risk stage II and stage III CRC patients in the high- and low-probability Lo IS subgroups. The results of the survival analysis showed that chemotherapy was associated with high-risk II and stage III CRC patients (Supplementary Figure S14). A test for an interaction between the probability of Lo IS and chemotherapy demonstrated that in either high-risk stage II or stage III, the benefit observed in the low-probability Lo IS patients [high-risk stage II (Figure 9): DFS, HR: 0.486 (95% CI: 0.280-0.842), P = 0.010; OS, HR: 0.441 (95% CI: 0.229-0.852), P = 0.015; stage III (Figure 10): DFS, HR: 0.464 (95% CI: 0.284-0.758), P = 0.002; OS, HR: 0.452 (95% CI: 0.266-0.770), P = 0.003; all P < 0.05 for interaction; Table 5] was superior to that observed in the high-probability Lo IS patients. The results indicated that chemotherapy significantly improved survival outcomes in the low-probability Lo IS group (high-risk stage II: P = 0.010 and P = 0.015; stage III: P = 0.002 and P = 0.003, respectively) but had no significant influence in the high-probability Lo IS group (high-risk stage II: P = 0.459 and P = 0.319; stage III: P = 0.535 and P = 0.449, respectively).