This meta-analysis was designed and reported in accordance with the Preferred Reporting Items for Systematic Reviews and Meta-Analyses (PRISMA 2020) guidelines (19). A detailed PRISMA checklist is provided in Supplementary Table 1. The study protocol was prospectively registered in the International Prospective Register of Systematic Reviews (PROSPERO; registration ID CRD42023493604).

A comprehensive literature search was performed in PubMed, Embase, Web of Science, and the Cochrane Library to identify studies published from database inception through January 1, 2026. The search strategy was developed according to the Population–Intervention–Comparator–Outcome (PICO) framework, as summarized in Supplementary Table 2. Detailed search strategies for each database are presented in Supplementary Table 3. All retrieved records were screened in two stages: (i) title and abstract review and (ii) full-text assessment against the predefined inclusion and exclusion criteria.

Studies were eligible for inclusion if they met all of the following criteria: (i) enrolled patients with histopathologically confirmed CRC; (ii) assessed TAN infiltration in primary tumor tissue using histopathological evaluation or transcriptomic estimation; (iii) reported the association between TAN levels and at least one survival outcome, including cancer-specific survival (CSS), overall survival (OS), or disease-free survival (DFS); studies reporting recurrence-free survival (RFS) were also included and analytically grouped with DFS; (iv) provided sufficient information to extract or calculate a hazard ratio (HR) with its 95% confidence interval (CI); and (v) were published as full-text articles in English. Exclusion criteria were as follows: (i) non-original publications, including reviews, editorials, letters, conference abstracts, and case reports; (ii) non-human studies; (iii) studies with overlapping or duplicate patient populations (with the largest or most complete dataset retained); (iv) studies with a sample size < 50 patients; or (v) studies lacking adequate data for effect size estimation.

A standardized data extraction form was used to collect key information from each eligible study. The extracted variables included: first author; year of publication; country; study period; sample size (and the number of patients evaluated for TANs); follow-up duration; tumor stage reported using either the Union for International Cancer Control (UICC) TNM classification or the Dukes staging system (with Dukes A–D broadly corresponding to UICC stages I–IV); receipt of preoperative treatment; TAN assessment method; TAN-related markers; cut-off criteria or methods for defining high versus low TAN density; cut-off values; TAN localization; survival outcomes (CSS, OS, and DFS); and data source (univariate and/or multivariate analysis). Studies reporting RFS were grouped with DFS as a unified time-to-event outcome based on overlapping endpoint definitions (Supplementary Table 4).

TAN localization was classified into four predefined groups. Whole-tumor section (WTS) referred to a global, section-level assessment of TAN density across the entire tumor tissue without regional subdivision into tumor nest (TN), tumor stroma (TS), or invasive margin (IM); TN was defined as intraepithelial tumor-infiltrating neutrophils within tumor nests; TS represented the stromal compartment surrounding tumor nests; and IM referred to the tumor–host interface at the invasive margin. TAN density was dichotomized into high and low levels according to the definitions used in the original studies. For meta-analytic modeling, “overall” effects were obtained by jointly modeling all available region-specific effect sizes (WTS, TN, TS, and IM), whereas region-specific analyses were conducted separately for each anatomical category.

For studies in which survival outcomes were reported only as Kaplan–Meier curves, survival probabilities at multiple time points were digitized using Engauge Digitizer (version 4.1) to enable reconstruction of time-to-event estimates. Hazard ratios (HRs) and corresponding 95% confidence intervals were then reconstructed according to the methods described by Tierney et al. (20), based on the estimated numbers at risk, events, and log-rank statistics where available. When essential information was incomplete, attempts were made to contact the corresponding authors; studies were excluded if no response or usable data were obtained. Data extraction was independently performed by two authors (R. Yang and C. Zhu), with discrepancies resolved by discussion with additional researchers.

The methodological quality of the included studies was evaluated using the Newcastle–Ottawa Scale (NOS) (21). This scale comprises three domains: selection (0–4 points), comparability (0–2 points), and outcome assessment (0–3 points). Two reviewers independently assessed study quality, and disagreements were resolved by discussion. Studies with NOS scores of ≥7 were considered high quality.

All statistical analyses were performed in R (version 4.5.0). For each study, HRs and corresponding 95% CIs were extracted or derived. HRs were inverted when necessary to ensure that HR < 1 uniformly indicated a favorable prognosis associated with higher TAN infiltration. The main survival outcomes were CSS, OS, and DFS; studies reporting RFS were analytically grouped with DFS. Statistical significance was defined as a 95% CI not crossing 1.

Univariate and multivariate HRs were synthesized separately. Because many studies reported multiple correlated effect sizes from different tumor regions (WTS, TN, TS, IM) or from multiple cohorts within the same publication, a multilevel random-effects model (rma.mv, REML) was applied, with effect sizes nested within studies (Study_ID/Effect_ID). Between-study and within-study variance components were estimated simultaneously, and heterogeneity was quantified using an I²-like statistic. For the main analysis, all available effect sizes from WTS, TN, TS, and IM were combined to obtain an overall pooled estimate, and region-specific analyses were conducted separately for each anatomical category. Pooled HRs were reported with 95% CIs and 95% prediction intervals (PIs).

When substantial heterogeneity was present (I²-like > 50%), prespecified subgroup analyses and meta-regression were performed to explore potential moderators. Moderators included study region, sample size, tumor stage, preoperative treatment status, TAN detection method, and TAN marker type. Meta-regression was conducted only when ≥10 studies were available, using multilevel random-effects models with cluster-robust variance estimators and CR2 small-sample correction to obtain robust standard errors (robust SEs). A P value < 0.05 for subgroup comparisons or regression coefficients was considered statistically significant.

Sensitivity analyses included leave-one-out procedures and influence diagnostics. Robustness was further evaluated by applying CR2-based robust SEs and by repeating analyses after excluding studies with lower quality, smaller sample sizes, or differing TAN assessment methods or markers. Publication bias was assessed using funnel plots, Egger’s regression, and Begg’s rank correlation tests; these tests and the trim-and-fill method were applied when ≥10 studies were available for a given outcome. All tests were two-sided, with P < 0.05 considered statistically significant.

A total of 6,380 records were identified through the database search. After removing 4,661 duplicates and screening titles and abstracts, 89 records remained, of which 21 were excluded as reviews, case reports, or letters. The remaining 68 full-text articles were assessed against the predefined eligibility criteria, and 43 were excluded. Ultimately, 25 studies fulfilled the inclusion criteria and were included in the meta-analysis (Figure 1).

This meta-analysis included 25 studies involving 10,356 patients with CRC, published between 2005 and 2025 (22–46). Geographically, 12 studies were conducted in Euro–American populations, 12 in East Asia (China and Japan), and one study analyzed publicly available transcriptomic datasets. Patient cohorts covered tumor stages I–IV, reported using either the UICC TNM classification or the Dukes staging system, with Dukes A–D broadly corresponding to UICC stages I–IV. Regarding preoperative treatment, 11 studies exclusively included treatment-naïve patients, 3 included patients who had received preoperative radiotherapy or chemotherapy, and 11 did not report preoperative treatment status. Detailed clinical and methodological characteristics of the included studies are summarized in Table 1.

A range of methodological approaches was used to quantify TANs. Immunohistochemistry (IHC) was the predominant technique (18 studies), employing markers such as CD66b (CEACAM8), MPO, CD15, CD177, and CD10. Five studies used hematoxylin–eosin (H&E)–based morphological assessment, one used manual cell counting, and one utilized transcriptomic deconvolution. Cut-off definitions for high versus low TAN density varied across studies and were based on median values, ROC curve analysis (Youden index), X-tile software, regression tree analysis, semi-quantitative scoring, or predefined absolute cell-count thresholds.

TAN localization was assessed in four predefined anatomical compartments: WTS, TN, TS, and IM. Nine studies evaluated WTS, 12 focused on TN, 4 on TS, and 5 on IM. Survival outcomes included CSS, OS, and DFS; studies reporting RFS were analytically grouped with DFS.

All included studies demonstrated high methodological quality, with NOS scores ranging from 7 to 9, indicating an overall low risk of bias.

In assessing potential publication bias and small-study effects, Egger’s regression test, Begg’s rank correlation test, and the trim-and-fill method were applied across six analytical sets corresponding to CSS, OS, and DFS in both univariate and multivariate models (Table 4). Overall, most tests yielded non-significant results (P > 0.10), indicating no clear evidence of publication bias. For example, in the CSS–uni model, the pooled HR was 0.61 (95% CI, 0.48–0.77), and the trim-and-fill procedure did not impute missing studies (k₀; = 0), yielding an identical adjusted estimate and a symmetrical funnel plot.

In the OS–uni model, Egger’s test suggested mild asymmetry (P = 0.034), whereas Begg’s test remained non-significant (P = 0.454), and trim-and-fill still imputed no additional studies (k₀ = 0), implying that any potential small-study effects had limited influence on the overall estimate. A similar pattern was observed in the DFS–uni model, where Egger’s test indicated modest asymmetry (P = 0.011), but no studies were added by trim-and-fill and the adjusted HR (0.75; 95% CI, 0.44–1.27) remained consistent with the original estimate. For multivariate models (CSS–multi, OS–multi, DFS–multi), Begg’s tests were non-significant and trim-and-fill did not impute additional studies (all k₀ = 0); although Egger’s test yielded a nominally significant P value for DFS–multi (P = 0.008), the absence of imputed studies and the similarity between original and adjusted estimates suggest that any small-study effects are unlikely to materially bias the results. Notably, because several models included fewer than ten studies (k < 10), the statistical power of Egger’s and Begg’s tests was limited, and these findings should therefore be interpreted with caution.

Taken together, no substantial publication bias or small-study effects were detected, supporting the credibility and stability of the pooled effect estimates.

To verify the robustness of the main findings, sensitivity analyses were performed across all six outcome models (Supplementary Table 6). Leave-one-out analyses showed minimal changes in pooled HRs (absolute difference < 0.15) with consistent effect directions, indicating that no single study materially influenced the overall estimates. Influence diagnostics (Cook’s D and DFBETAS) did not identify high-leverage or influential outlier studies, and all models exhibited stable convergence during weighted least-squares refitting, suggesting that the results were not driven by aberrant observations.

Robust standard error (robust SE) adjustments produced central HRs that were highly consistent with the primary analyses (differences < 5%), with only modest widening of confidence intervals, indicating low sensitivity to variance assumptions. Directed exclusion strategies based on study quality, sample size, detection method, and TAN markers further confirmed the stability of the results: pooled HRs remained within 0.46–0.66 after removing low-quality or small-sample studies, and stronger protective effects were observed in analyses restricted to IHC-based assessments or CD66b markers (HR ≈ 0.27–0.50).

Overall, the sensitivity analyses demonstrated consistent robustness across all outcome models. High TAN infiltration showed stable and significant improvement in CSS (HR ≈ 0.58–0.66; I²-like ≈ 0%) in both univariate and multivariate analyses, while OS and DFS exhibited concordant protective trends (HR ≈ 0.46–0.73) without effect reversal. These results collectively support the reliability and statistical stability of the association between high TAN infiltration and favorable survival outcomes.

Given the low heterogeneity observed in CSS models, our exploration of heterogeneity focused on OS and DFS. In the OS univariate model, overall between-study heterogeneity was substantial (τ² = 0.592; I²-like = 81.1%). Subgroup analyses indicated that study region and tumor stage were major drivers of this variability. Among non-Asian (Euro–American) studies, heterogeneity markedly decreased (τ² = 0.0064; I²-like = 7.2%), and in the multivariate OS model, heterogeneity in Euro–American cohorts was essentially abolished (τ² = 0; I²-like = 0%). Similarly, in the multivariate OS model stratified by UICC stage, the stage I–III subgroup showed no residual heterogeneity (τ² = 0), whereas studies including stage I–IV disease retained high heterogeneity (τ² = 0.4099; I²-like = 88.7%), suggesting that differences in stage composition may underlie inconsistent effect estimates. For DFS, overall heterogeneity was likewise high (τ² = 1.0533; I²-like = 92.1%), and subgroup patterns by region and stage paralleled those observed for OS. Detection method (IHC vs H&E) modestly attenuated heterogeneity in selected models, whereas sample size and preoperative treatment had relatively minor impact.

On this basis, multilevel meta-regression analyses were conducted to quantify the contribution of key moderators to between-study variance (Supplementary Table 7). In the OS univariate model, none of the covariates reached conventional statistical significance, although UICC stage showed a positive trend (β = 0.62; P = 0.142), which remained directionally consistent after CR2-robust variance correction (p_CR2 = 0.176). In contrast, the multivariate OS model showed a marked reduction of heterogeneity after inclusion of covariates (τ² decreasing from 0.933 to approximately 0), and the overall regression was highly significant (QM = 75.67; P < 0.001). UICC stage (β = 2.31; 95% CI, 1.70–2.91; P < 0.001) and TAN detection marker (β = −1.25; 95% CI, −1.77 to −0.73; P < 0.001) emerged as the principal explanatory variables, while smaller sample size (≤200 patients) was also associated with lower HRs (β = −1.38; P = 0.001). After CR2 correction, stage and marker remained statistically significant (p_CR2 < 0.01), indicating that these two factors robustly account for nearly all between-study heterogeneity in OS; the extremely small CR2 p-value for sample size likely reflects numerical instability rather than a genuinely stronger effect.

For DFS in the univariate setting, several covariates significantly explained heterogeneity: UICC stage (β = −1.23; P = 0.0025), detection marker (β = −0.96; P = 0.021), and detection method (β = −0.61; P = 0.048), suggesting that differences in staging and assessment strategies may yield divergent effect directions. After CR2 adjustment, the associations for UICC stage and detection marker were attenuated and became non-significant (p_CR2 ≥ 0.27), although the directions of effect remained unchanged; CR2-based tests were not available for detection method owing to convergence limitations. In the multivariate DFS model, only UICC stage remained significant (β = 2.01; 95% CI, 1.14–2.87; P < 10^-5), explaining approximately 90% of the between-study variance (τ² reduced from 1.053 to 0.045). However, because the comprehensive multivariate model failed to converge stably in the presence of multiple correlated moderators, CR2-robust testing could not be reliably applied.

Taken together, subgroup and meta-regression analyses yielded highly consistent patterns: tumor stage and TAN detection markers are the predominant sources of heterogeneity in OS and DFS models, whereas detection method, sample size, and preoperative treatment exert comparatively modest influence. The staging variable remained consistently significant or directionally stable across multiple models, suggesting that differences in tumor stage distribution are a key determinant of variability in the prognostic impact of TAN infiltration.