The landscape of tissue-resident microbiota across normal, polyp, and colorectal cancer tissues

Domain: chemistry
Published: 2026-02-24
Editorial outcome: major_revision

Abstract
Emerging evidence suggests that tissue-resident microbiota (TRM) is associated with tumor biology; however, their distribution and compositional characteristics across different colorectal tissue types remain incompletely defined. Here, we conducted a comprehensive cross-sectional analysis of TRM distribution and abundance across 1134 clinical specimens, including normal mucosa, precancerous polyps, and colorectal cancer (CRC) tissues. Our results reveal distinct microbial profiles among these diagnostic groups, with consistent differences in community composition between normal, polyps, and CRC samples. Integrative analyses further identified microbial signatures capable of distinguishing tissue categories and reflected differences in the local tumor microenvironment. In contrast, intratumoral microbiota composition showed subtle variation across established tumor stages and was not associated with clinical outcomes. These findings define diagnostic group-associated patterns of TRM in colorectal tissues and establish a foundation for future mechanistic investigations into the biological roles of TRM in colorectal cancer.

Introduction
Colorectal cancer (CRC) is a major global health burden, ranking as the third most diagnosed cancer and the second leading cause of cancer-related deaths according to GLOBOCAN 2020 estimates1. Despite therapeutic improvements, CRC continues to carry a high mortality burden, emphasizing the pressing need to understand the biological mechanisms that initiate and sustain tumor development.

Beyond conventional risk factors, including genetic predisposition, epigenetic modifications, and stromal microenvironment alterations, the gut microbiota has emerged as a critical mediator of CRC susceptibility and progression. Mounting evidence suggests that gut microbiota influence colorectal tumorigenesis through multifaceted mechanisms, such as microbial metabolite production, intricate immune system modulation, and direct epithelial interactions2–7. In addition to gut microbiota, recent discoveries of tissue-resident microbiota (TRM), defined as bacterial communities that establish stable colonization within host tissues, have opened new avenues for understanding colorectal tumorigenesis. Unlike luminal microbiota, TRM establishes direct physical interactions with host cells, suggesting they may play more immediate roles in epithelial transformation8–10. While TRM profiles have been linked to cancer progression in other malignancies, their specific contributions to colorectal tumor initiation and progression remain poorly understood11–14.

In this study, we conducted a comprehensive multi-cohort analysis of the TRM landscape across 1134 clinical specimens, encompassing normal mucosa, precancerous adenomatous polyps, and invasive carcinomas. Using pre-amplified 16S rRNA sequencing coupled with advanced bioinformatic pipelines, we identified distinct differences in microbial community composition among these tissue categories. Notably, TRM differences were most pronounced between diagnostic groups, with minimal variation across later clinical stages or patient survival. These findings define tissue-associated patterns of TRM composition in colorectal tissues and provide a basis for future mechanistic investigations aimed at understanding the biological relevance of TRM in colorectal cancer.

Results
Colorectal cancer is associated with distinct taxonomic differences in TRM
To systematically profile the TRM in CRC, we conducted a cross-sectional study utilizing CRC and normal colorectal tissue specimens from the Sir Run Run Shaw Hospital of Zhejiang University Biobank. The baseline demographic and clinical characteristics of the included participants are summarized in Tables S1 and S2. After quality control, 1,022 samples were retained, with 818 (182 normal, 68 polyps, 568 CRC) assigned to the discovery set and 204 (45 normal, 17 polyps, 142 CRC) to the validation set.

To characterize the TRM, we employed an optimized 16S rRNA gene sequencing protocol in which the V4 region was pre-amplified using biotin-labeled primers, followed by streptavidin bead-based enrichment to enhance specificity. High-throughput sequencing was subsequently performed with stringent contamination controls to ensure data reliability (Fig. 1A,B). Parallel culture-based assays confirmed viable bacteria in control, polyp, and CRC tissues (Fig. 1C). Quantitative analysis revealed a statistically significant increase in bacterial load in CRC tissues compared to both normal and polyp tissues, while no significant difference was observed between normal and polyp tissues (Fig. 1C).

Due to the relatively small number of polyp samples, we first focused on analyzing compositional differences in TRM between normal and CRC tissues. We conducted multivariate analyses using Bray-Curtis distances. Principal Coordinates Analysis (PCoA) and Detrended Correspondence Analysis (DAC) revealed distinct clustering patterns separating normal intestinal tissues from CRC samples at the genus level (Fig. 1D,E). Taxonomic profiling demonstrated marked variations in bacterial abundance between groups (Fig. 1F). Linear discriminant analysis effect size (LEfSe) identified specific bacterial taxa differentially abundant between normal and CRC tissues (Linear Discriminant Analysis (LDA) score > 2.0, q-value<0.05; Fig. 1G). Notably, normal tissues exhibited enrichment of Pseudomonadaceae, aligning with previous reports that highlight Pseudomonas aeruginosa’s tumor-suppressive capability15,16. In contrast, CRC tissues showed enrichment of several bacterial taxa with an established role in promoting CRC development. This included Fusobacterium specie17,18, Bacteroides species19, and Enterobacteroides species20,21. These findings collectively reinforce the established association between specific microbial communities and colorectal tumorigenesis while providing new insights into their tissue-specific distribution patterns.

To distinguish TRM from luminal microbiota, we compared tissue-derived microbial profiles with fecal microbiota from CRC patients aggregated across seven public datasets22–28. PCoA analysis revealed a clear separation between fecal and tissue microbiomes (PERMANOVA, p-value<0.001; Figure S1A), accompanied by significant differences in alpha diversity (Fig. S1B) and dominant genera composition (Fig. S1C). Parallel differential abundance analyses identified partially overlapping but weakly correlated taxa (ρ=0.106, p-value=0.451; Fig. S1D), indicating that TRM represent distinct microbial communities not adequately captured by fecal profiling.

TRM shows a limited association with clinical features in CRC
While substantial differences exist between TRM in normal versus CRC tissues, their potential to discriminate clinical features of CRC remains unexplored. To address this, we conducted PCoA analysis of intratumoral microbiota stratified by multiple clinical parameters at the genus level. Contrary to our expectations, microbial communities showed no significant segregation when analyzed by age, gender, differentiation status, tumor node metastasis classification (TNM) stages, anatomical, or survival status (Fig. 2A–F). We observed the same pattern after controlling inter-individual confounding factors by analyzing on samples of matched tumor and adjacent normal tissue pairs. (Figures S2A–F). These findings indicate that while TRM composition differs fundamentally between normal and cancerous tissues, these microbial signatures do not appear to reflect the clinical heterogeneity of established CRC.

To investigate potential associations between TRM and clinical features, we performed a hierarchical clustering analysis of bacterial genera abundance across the aforementioned clinical characteristics. The resulting heatmap revealed diverse abundance patterns of specific bacterial genera across different clinical features, with certain taxa showing relative enrichment in particular subgroups (Fig. 2G). However, these associations failed to reach statistical significance, confirming no robust correlation between TRM and clinical features in CRC.

To further explore the relationship between TRM and clinical features in CRC, we developed random forest classification models for key clinical parameters and assessed their performance through receiver operating characteristic (ROC) curve analysis. In the training cohort, area under the curve (AUC) values ranged from 0.565 to 0.714 (Fig. 2H). The validation cohort showed comparable limitations, with AUC values between 0.541 and 0.762 (Fig. 2I). Analyses on matched tumor-normal tissue pairs from the same individuals yielded similar results (Fig. S2G). This suggests that the observed microbial patterns are independent of intra-individual confounders.

To evaluate potential clinical confounders, we performed multivariate analyses incorporating age, gender, tumor stage, location, differentiation grade, and chemotherapy status. Cox proportional hazards analysis, adjusting for clinical confounders, revealed that no individual microbial taxon marginally associated with survival outcome (Fig. S3A). To further investigate the joint impact of the microbiota on survival outcome, the Lasso-Cox model was used to select important taxa, and the selected features were further analyzed using multivariable Cox regression. The results revealed that their Hazard Ratios (HR) were near 1.0, indicating a negligible contribution to survival compared to clinical covariates (Fig. S3B). Analysis on 5-year survival outcome by using logistic regression gave similar results (Fig. S4A,B). Furthermore, community-level models (RandomForest) showed that adding microbial data did not improve the predictive performance beyond clinical variables alone (Fig. S4C).

Overall, these results demonstrate that no strong associations were observed between TRM composition and survival outcomes or other clinical features under the applied statistical frameworks in our cohort. Consistent suboptimal model performance (AUC<0.7) indicates that TRM lacks sufficient discriminatory power to serve as a reliable biomarker for CRC tumor characteristics or clinical outcomes.

TRM exhibits distinct profiles across normal, polyp, and CRC tissues
Our observation of stable TRM composition across established tumor stages, contrasted with marked differences between normal and CRC tissues, indicates cross-sectional differences in TRM profiles between tissue types rather than stage-dependent variation within colorectal cancer. To investigate this hypothesis, we expanded our analysis to include precancerous polyps (68 in the discovery cohort, 17 in the validation cohort), encompassing samples from normal mucosa, adenomatous polyps, and invasive carcinoma.

Through PCoA and DCA analyses based on Bray-Curtis distances, we revealed three distinct microbial ecologies corresponding to each histological stage. The polyp-associated microbiota exhibited a distinct profile relative to both normal and colorectal cancer tissues, demonstrating clear microbial differences among normal, polyp, and CRC samples (Fig. 3A,B). Taxonomic profiling identified specific bacterial families with distinct abundance patterns across histological groups (Fig. 3C). Fusobacteriaceae, Streptococcaceae, Prevotellaceae, Leptotrichiaceae, and Burkholderiaceae showed the lowest relative abundance in normal tissue, intermediate abundance in polyp tissue, and the highest abundance in CRC tissue. Conversely, Coriobacteriaceae, Bifidobacteriaceae, and Tannerellaceae were most abundant in normal tissue. Genus-level analysis confirmed these patterns, with Fusobacterium, Streptococcus, Peptostreptococcus, Leptotrichia, Massilia, and Prevotella displaying the lowest abundance in normal tissue, intermediate levels in polyps, and the highest abundance in CRC, while Phascolarctobacterium, Collinsella, Bifidobacterium, Pseudomonas, and Akkermansia were enriched in normal tissue, reduced in polyps, and further depleted in CRC samples (Cochran-Armitage trend test, q-value<0.05). Notably, several genera, including Bacteroides, Porphyromonas, and Campylobacter exhibited markedly higher abundance in CRC compared with both normal and polyp tissues (Fig. 3D,E).

These findings demonstrate distinct TRM profiles across normal, polyp, and colorectal cancer tissues. While each histological category is characterized by a distinct microbial signature, TRM composition appears relatively stable between polyps and established carcinoma, highlighting the relevance of TRM profiles in distinguishing diagnostic tissue categories, rather than across clinical strata within colorectal cancer.

TRM signatures distinguish colorectal tissue states
Building upon our observation of distinct TRM profiles in normal, polyp, and CRC tissues, alongside a relative stability of TRM composition across established tumor stages, we next sought to determine whether these microbial signatures could reliably differentiate tissue states. We initially developed a three-class random forest model (normal mucosa, polyps, CRC), which achieved acceptable performance (AUC>0.8; Fig. S5). However, the limited number of polyp samples (n=68) constrained classification precision. To address this and better capture microbial differences between tissue types, we employed pairwise models. Specifically, we compared polyps with Stage I CRC (n=61). These groups were selected because of comparable sample sizes, allowing a balanced cross-sectional comparison of microbial profiles between polyp tissue and early-stage colorectal cancer tissue.

To identify core microbial features, we first ranked taxa by their importance scores derived from the random forest algorithm and selected the top 10 candidates for further evaluation. These taxa were then filtered based on prevalence, excluding those present in only a small subset of samples. Subsequent refinement incorporated pairwise correlation analyses to prioritize taxa with low inter-correlation, thereby minimizing redundancy. Finally, we excluded taxa with ambiguous taxonomic classification, retaining only those with clear biological relevance for modeling (highlighted in boldface; Fig. 4A,B).

To optimize feature parsimony, we monitored both training and validation performance (measured by mean of AUCs) using 5-fold cross-validation, and selected the number of variables that achieved the highest validation performance while controlling the risk of overfitting based on training performance (Fig. 4C; See Methods). The same selection protocol was applied to the normal-vs-polyp (Fig. S6) and normal-vs-CRC (Figure S7) models. Notably, we found that four key taxa: Pseudomonas, Pigmentiphaga, Parvimonas, and Fusobacterium were sufficient to achieve robust discrimination between Polyp and CRC_I tissues (highlighted in red; Fig. 4C). Similarly, six taxa: Pseudomonas, Parabacteroides, Bacteroides, Pigmentiphaga, Fusobacterium, and Collinsella provided effective discrimination between normal mucosa and polyps (highlighted in red; Fig. 4D).

The resulting models consistently classified tissue states with high performance within our cohort, achieving AUCs of 0.911 (training) and 0.954 (validation) for normal-versus-polyp, 0.986 and 0.996 for polyp-versus-CRC, and 0.944 and 0.917 for normal-versus-CRC comparisons (Fig. S8A,B). To ensure the robustness of these findings and reduce the risk of overfitting, we implemented a 5-fold cross-validation strategy across the discovery cohort for model evaluation. Additionally, up-sampling was applied during model training to address class imbalances, particularly the relatively small number of polyp samples, which produced more conservative and reliable estimates. Across all classification tasks, this rigorous validation yielded mean AUCs of 0.87-0.90, with standard deviations reflecting consistency across folds (Fig. 4E), demonstrating that the models capture reproducible differences in tissue-associated microbial composition. Notably, 55 patients in our dataset contributed both tumor and matched normal samples. To further assess model generalizability, we performed Leave-One-Patient-Out (LOPO) validation, which ensures that samples from the same patient are never included in both the training and validation sets. This analysis yielded results consistent with our original cross-validation findings (Fig. S8C).

Intersection analysis of discriminatory taxa across models (Fig. S9A) identified key overlapping features: Pseudomonas, Pigmentiphaga, and Fusobacterium were shared between normal-vs-polyp and polyp-vs-CRC_I models, while Pseudomonas and Fusobacterium appeared consistently across all three models. Differential abundance analyses across various clinical features showed that most microbial taxa did not differ significantly, indicating that these signatures are largely independent of common clinical variables (Fig. S9B–G). This stability across diverse patient characteristics highlights the robustness of TRM composition as a discriminator of colorectal tissue states.

To evaluate generalizability, we applied our models to an independent external cohort comprising 112 samples collected from a separate clinical center. The model demonstrated robust performance, achieving AUCs of 0.779 for Normal vs. Polyp, 0.796 for Polyp vs. CRC, and 0.767 for Normal vs. CRC (Fig. 4F). These AUCs provide a realistic estimate of model performance in an external clinical setting. Although modestly lower than the internal validation metrics, such reductions are expected in rigorous external testing and are consistent with patterns reported in prior large-scale microbiome studies27. These findings indicate that the microbial patterns identified in our study are reproducible across distinct populations and clinical contexts.

To further examine the influence of clinical heterogeneity, we stratified CRC samples by tumor location (right- vs. left-sided) and microsatellite instability (MSI) status. Due to the limited number of MSI-High samples (n=28), analyses focused on the predominant MSS subgroup. Differentially abundant taxa were largely consistent across Right-sided, Left-sided, and MSS subgroups, and models rebuilt using stratum-specific microbial signatures maintained similar performance (Fig. S10). These findings indicate that TRM compositional differences are robust and largely independent of common clinical variables such as tumor location or MSI status, reinforcing the generalizability of the observed microbial patterns across diverse CRC subtypes.

Taken together, these findings show that TRM composition is distinct across the normal, polyp, and CRC status and remains stable across clinical subgroups and cohorts. The recurrent presence of taxa such as Pseudomonas, Pigmentiphaga, and Fusobacterium points to reproducible microbial patterns associated with colorectal tumorigenesis.

Validation of model-identified microbial signatures in clinical specimens
To experimentally validate the microbial signatures identified by our machine learning models, we performed fluorescence in situ hybridization (FISH) on clinical specimens representing normal mucosa, polyps, and CRC tissues. This approach enabled direct visualization and quantification of the signature bacterial taxa at their native tissue microenvironments.

For taxa shared between normal-vs-polyp and polyp-vs-CRC_I models, we observed distinct colonization patterns consistent with our taxonomic profiling results. Pseudomonas fluorescence intensity declined from normal tissue to polyps and CRC (Fig. 5A,B), whereas Pigmentiphaga was least abundant in polyp specimens (Fig. 5C, D). Notably, Fusobacterium intensity increased across tissue types, reaching its highest levels in CRC samples (Fig. 5E,F), in line with its established role in CRC development. FISH analysis of model-specific bacterial signatures further validated our machine learning models. In normal versus polyp comparisons, Parabacteroides, Collinsella, and Bacteroides all maintained significantly higher abundance in normal tissues, while Parvimonas was exclusively enriched in CRC specimens (Fig. S11).

To assess the consistency between FISH and 16S rRNA sequencing data, we conducted a group-level comparative trend analysis. For each diagnostic group (normal, polyp, and CRC), we compared mean relative abundance values from sequencing with corresponding FISH quantification results (Fig. S12A). Both methods revealed highly concordant patterns across tissue types. For example, Fusobacterium abundance increased from normal mucosa to tumor tissues in both datasets (Fig. S12B), demonstrating the robustness and reproducibility of this microbial trajectory across distinct analytical platforms.

Together, these experimental results demonstrate strong concordance between computational predictions and experimental observations at the single-microbe level. Our findings reinforce the biological validity of TRM-based microbial signatures and provide experimentally supported candidates for future mechanistic studies on microbiota-host interactions in colorectal cancer development.

Methods
This study was conducted in accordance with the ethical principles of the Declaration of Helsinki and was approved by the Institutional Review Board (IRB) of Sir Run Run Shaw Hospital (SRRSH) (Approval No. IRB-2022-0423). Written informed consent was obtained from all patients prior to the collection and analysis of their samples.

Humans samples collection
This retrospective cohort study utilized samples and data identified by systematically querying the SRRSH Biobank and clinical archives for patients who underwent colonoscopy or surgical resection between January 2017 and January 2018. Participants were categorized into three groups based on archived histopathological reports: CRC, polyp, and normal controls. The polyp group encompassed various pathological types, including adenomatous and non-adenomatous polyps.

Inclusion and exclusion criteria
The following criteria were applied to select the final study cohort:
Inclusion Criteria: Availability of a definitive histopathological diagnosis of CRC, polyp (any type), or normal mucosa; age between 18 and 90 years at the time of procedure; no documented history of anticancer therapy prior to sample acquisition; availability of high-quality, frozen tissue specimens in the Biobank; and completeness of essential clinical and pathological data.
Exclusion Criteria: Documented use of antibiotics, prebiotics, or probiotics within one month before the procedure; co-diagnosis of other gastrointestinal diseases (e.g., inflammatory bowel disease); prior history of colorectal resection; and insufficient tissue quality or missing critical clinical data.

Sample Collection, Quality Control and Processing
Normal, polyp, and CRC samples were collected from different individuals, meeting the above criteria. Among these samples, there are 55 paired normal and CRC tissues obtained from the same patients. Fresh human tissues were collected in the sterile surgery room from the medical centers of SRRSH. Fresh normal intestinal tissues (located at least 10 cm away from the colorectal cancer lesion), polyp tissues, and colorectal tumor tissues were immediately transferred to germ-free 15 ml centrifuge tubes containing sterile DMEM culture medium. The −80°C frozen human tissue was obtained from the biobank of SRRSH. Each tissue sample was rinsed three times with PBS, each rinse lasting 5 minutes. Then, using autoclaved forceps and sterile blades, the surface tissue was dissected off, and the central core (about 200 mg) was retained for subsequent DNA extraction.

To minimize environmental contamination, the core regions of each tissue sample were carefully dissected prior to DNA extraction. Bacterial load across different tissue types was quantified by qPCR targeting the 16S rRNA gene. While all environmental control samples fell well below the threshold of 10³ CFU/g (as established by standard curve analysis) (Fig. S13), all examined tissue samples significantly exceeded this level. Consequently, samples below this threshold were excluded from subsequent analysis.

TRM sample preparation
There are three major challenges that affect bacterial DNA extraction: the presence of an enormous host genome, release of bacterial DNA, and contamination in reagents. To obtain high-quality bacterial DNA, we tested different DNA extraction kits, optimized the amount of tissue, and employed bead shearing techniques. However, contamination in reagents is difficult to eliminate completely. Therefore, we established stringent environmental controls that underwent the exact same procedure as the samples, allowing us to generate a comprehensive contamination landscape. To determine whether to eliminate the host cell genome, we tested two DNA extraction kits: the QIAamp DNA Microbiome kit (QIAGEN #51704), which can deplete host DNA, and the QIAamp Power Fecal (pro) DNA kit (QIAGEN #51804). All tested kits included a bead-beating step, which was performed in a TissueLyser II (QIAGEN, 85300) for 10 minutes at 30 Hz. The bacterial DNA quantity decreased without the bead-beating step. The QIAamp DNA Microbiome kit contains steps to lyse host cells with buffer AHL and Benzonase to digest host DNA. However, these additional steps consistently decreased the sample yield for tissues, resulting in lower sensitivity of qPCR quantification and frequent ambiguous results. Therefore, we decided to use the QIAamp Power Fecal (pro) DNA kit, which does not include a host genome elimination step, for subsequent bacterial DNA extraction.

For tissue processing, we determined that 150 mg of tissue should be used for DNA extraction, as larger amounts would decrease extraction efficiency due to the limited binding capacity of the spin column. Tumor tissue samples were ground in 1 ml sterile PBS. After centrifugation at full speed for 10 minutes, the supernatant was discarded, and the pellet was either extracted or frozen at −80°C.

To set up controls, as contamination in reagents is difficult to eliminate completely, we established appropriate negative controls. We used PBS as tissue surrogate to undergo the same processing steps as tissue samples. The resulting DNA product served as the environmental bacteria control (EC). Additionally, we used DEPC-treated water as the no-template control (NTC) when running the qPCR quantification assay.

qPCR quantification
YEASEN Hieff qPCR reaction mixture was used for bacterial quantification. The V9 region of the 16S ribosomal RNA was amplified using the following primers: 1369 Forward 5’-CGGTGAATACGTTCYCGG-3’, 1492 Reverse 5’-GGWTACCTTGTTACGACTT-3’. To calculate bacterial DNA concentration in the samples, we generated a standard curve using Escherichia coli DNA. No-template controls (NTCs) were included in all reactions. The qPCR program consisted of an initial denaturation at 94 °C for 10 minutes, followed by 40 cycles of 94 °C for 15 seconds and 60 °C for 60 seconds. Raw threshold cycle (Ct) values were normalized against the E. coli DNA standard curve for accurate quantification.

16S library construction
To increase the recovery efficiency of PCR products in the first step during library preparation, the overhang of the amplification primer was modified by 5’ biotin. After the first-round biotin PCR, Dynabeads MyOneTM Streptavidin C1 Beads (Thermo Fisher Scientific, #5002) were used to pull down the biotin PCR products following the manufacturer’s instructions. The Streptavidin beads containing the first round PCR product were then directly subjected to the second round index PCR. Specifically, dual indices from TruePrep® Index Kit V3 for Illumina® (5 ul for each index, Vazyme #TD203) were added to target amplicons in a second round of PCR using 23 Kapa HiFi Hotstart ReadyMix DNA polymerase. Then the Ampure XP beads (Beckman, A63881 AMPure XP) were added to clean up the final libraries, and the quality of these purified libraries was checked by Fragment Analyzer-12/96 (GENE-QC006). Finally, samples were sequenced on the Novaseq with 2× 250 bp paired-end reads (Novogene, Tianjin).

For tumor-resident microbiota library construction, we used PBS that underwent all the sample preparation procedures as an environmental negative control and DEPC-treated water as no template control, as well as E. coli as the positive control. The V4 hypervariable region of 16S was amplified with biotinylated primer set 515F (0.4 uM, 5¢-TCGTCGGCAGCGTCAGATGTGTATAAGAGACAGGTGYCAGCMGCCGCGGTAA -3¢) and 805R (0.4 uM, 5¢- GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAGGGACTACNVGGGTWTCTAAT-3¢), each with overhang adapter sequences (IDT), using 23 Kapa HiFi Hotstart ReadyMix DNA polymerase (Kapa Biosystems). Amplification was performed with an initial heating step of 95 °C for 3 minutes, 25 cycles of 30 seconds at 95 °C, 30 seconds at 55 °C, 30 seconds at 72 °C, and final extension of 72 °C for 5 minutes. The PCR cycles were optimized to 30 cycles. PCR products were then purified and enriched by Dynabeads MyOne Streptavidin C1 Beads (Thermo Fisher Scientific, #65002) following the manufacturer’s instructions. Biotin beads containing the first round PCR products were directly subjected to the second round index PCR. Dual indices from TruePrep® Index Kit V3 for Illumina® (5 ul for each index, Vazyme #TD203) were added to target amplicons in a second PCR using 23 Kapa HiFi Hotstart ReadyMix DNA polymerase (25 ul). PCR conditions were 3 minutes at 95 °C, with 8 cycles of 30 seconds at 95 °C, 30 seconds at 55 °C, 30 seconds at 72 °C, and final extension of 5 minutes at 72 °C. Libraries were purified by the AMPure XP (0.8X) bead before quantification by using the Qubit DNA assay (Thermo Fisher Scientific) and qPCR, or quality check by Fragment Analyzer-12/96 (GENE-QC006). Negative controls, including no template controls and environmental controls, were included in all sequencing runs. Samples were sequenced on the Novaseq with 2× 250 bp paired-end reads (Novogene, Tianjin).

Microbiome data analysis
Processing amplicon sequencing data. The sequencing data were split using de-multiplexing tools bcl2fastq2 (double check with the sequencing service company for the de-multiplexing method they used) to generate fastq format files for every sample using barcode sequences. All reads were analyzed using a standardized amplicon bioinformatics pipeline using vsearch v2.14.243 and usearch v1044. First, paired-end reads were merged using “vsearch -fastq_mergepairs” with default parameters. Then, adaptor trimming and low-quality read removal were carried out by “vsearch -fastx_filter” with a maximum expected error rate of 0.01. After quality control filtering, reads were processed. Amplicon templates were constructed using the quality control-passed reads. Reads were dereplicated using “vsearch -derep_fulllength” with singleton reads discarded. All remaining unique sequences were denoised by the unoise3 method to get candidate sequence features, which balances between resolution and specificity45. Lastly, chimeric features were removed using “vsearch -uchime_ref” against the SILVA reference release 123.

The abundance table was generated by mapping all quality control passing reads onto the templates sequence features using “vsearch -usearch_global” with an identity threshold of 0.97. Sequence features were classified by “vsearch -sintax” with a cut offset to 0.6. Non-specific sequence features from eukaryote, and mitochondria were filtered out. Only samples with enough prokaryote-specific reads were kept using the cutoff of 1000 for human samples, respectively. All kept samples were normalized into the same number of reads for a fair comparison (1000 reads for human) using the “rrarefy” function of R package vegan (v2.6-4)46.

Contamination correction. The contamination effect in each human sample was numerically corrected using the data from qPCR quantification experiment. First, the relative abundance of each bacterial species was obtained from the 16S rRNA sequencing data. Then, the absolute bacterial amount per species per sample was calculated according to the qPCR quantification value and relative abundance within the sample. The contamination bacteria from background environment was measured as the median absolute amount of a species among all negative control samples, formulated as
C = c * Q, C = Median C
ij ij i j ij
where c is the relative abundance of species j from sample i, i.e. the percentage of species j data among the data of all species in sample i; Q is the qPCR quantification result of sample i in unit of CFU/g. C is the overall contamination effect of species j, as the median relative abundance of species j among all samples.

The correction simply took off the overall contamination effect of species j from the measured species j in each sample: A = a * Q - C
ij ij i j
where a is the relative abundance of species j from sample i; Q is the quantification of total bacteria of sample i; A is the corrected amount for species j of sample i.

Contamination correction was based on binomial tests of every taxon’s abundance between samples and negative controls. The probability p of the binomial distribution was estimated as the taxon’s occurrence frequency in negative control samples. The number of occurrence samples and the total number of samples were assigned as x and n in the binomial test, respectively. Taxon with p-value < 0.05 were kept in the following analysis.

Beta diversity analysis
Beta diversity analysis was performed to assess compositional differences between microbial communities at the genus level. PCoA and Non-Metric Multidimensional Scaling (NMDS) were implemented using the vegan package (v2.6-4) in R (v4.3.1) to visualize multivariate patterns in community structure. Both methods utilized Bray-Curtis dissimilarity matrices calculated from genus-level abundance profiles. While PCoA preserves original distance relationships through eigen decomposition47, NMDS optimizes the stress function to represent rank-order distances in reduced dimensions48, providing complementary perspectives on community dissimilarity.

Statistical significance of group differences was quantified using PERMANOVA through the adonis2 function (vegan package) with 999 permutations. This non-parametric multivariate analysis of variance tests centroid dispersion and compositional differences between predefined groups while accounting for potential covariates49. Model construction followed sequential (Type I) sums of squares to control for confounding variables. Both visualization approaches and PERMANOVA outcomes were interpreted synergistically to characterize microbial community variation across experimental conditions.

Alpha diversity analysis
Alpha diversity analysis was conducted to evaluate within-sample microbial diversity across taxonomic hierarchies (family, genus, and species levels). The Shannon index (quantifying community richness and evenness) and Simpson index (emphasizing dominant taxa proportions) were computed using the diversity function in the vegan package (v2.6-4)50. Statistical comparisons between experimental groups were performed through pairwise Student’s t-tests with Benjamini-Hochberg (BH) correction for multiple comparisons, implemented via the pairwise.t.test function in the base R stats package51. This parametric approach was selected after confirming normality assumptions using Shapiro-Wilk tests and homogeneity of variance via Levene’s test. Results were visualized through boxplots annotated with significance levels (p<0.05, p<0.01, p<0.001) to highlight inter-group diversity differences across taxonomic resolutions.

Discovery of microbial taxa differentially abundant among the three groups
Microbial biomarker discovery was performed using Linear Discriminant Analysis Effect Size LEfSe through the microeco package (v0.11.0) in R52. The analysis identified phylogenetically discriminative features between normal tissues and CRC tumor microbiomes using a three-stage approach: (1) Kruskal-Wallis tests53 detected genus-level abundance differences across groups (α=0.25), (2) pairwise Wilcoxon rank-sum tests54 confirmed subgroup consistency, and (3) LDA quantified effect sizes (log10 LDA score > 2.0). Low-abundance features were filtered (prevalence threshold: 5% across samples via filter_thres=0.05) prior to analysis. The top 20 discriminative taxa with the highest LDA scores were selected for visualization, annotated with their relative abundance distributions (Fig. 1F,G). Multiple testing correction was applied using the BH procedure55 during feature selection.

We next sought to identify microbial genera that exhibit a significant trend in abundance across the normal, polyp, and tumor. To achieve this, we performed a Cochran-Armitage trend test56. First, for each genus, we converted its relative abundance in every sample into a binary variable (e.g., high vs. low) based on a quantile threshold calculated from its distribution across all samples. The Cochran-Armitage test was then applied to this binary data to assess the statistical significance of any monotonic trend across the three groups by the R package DescTools (v0.99.60). The resulting p-values were corrected for multiple hypothesis testing using the BH method.

Composition and distribution of dominant bacterial communities
We constructed stacked bar plots to illustrate the relative mean abundance of microbial communities across different groups, including the environmental group, fecal group, normal tissue group, polyp tissue group, and CRC patient group (Fig. 3C,D). The analysis focused on the top 19 most abundant taxa at the family and genus levels, while the remaining taxa were collectively categorized as “others.” For each group, we calculated the proportion of the mean abundance contributed by these selected taxa. Similarly, we presented the proportional distribution of dominant genera in our TRM dataset and in fecal data from public databases at the genus level.

Microbial signature profiling was conducted to identify genus-level abundance patterns associated with clinical phenotypes. The top 30 most abundant genera across all samples were selected based on cumulative relative abundance. Group-specific mean values of log10(x+1)-transformed abundances were calculated for comparative groups: 1) normal tissues vs. CRC tumors, and 2) CRC subgroups stratified by clinical characteristics (e.g., TNM stage, histological grade). Visualization was achieved using the pheatmap package (v1.0.12)57 with hierarchical clustering (Euclidean distance, complete linkage) applied to both rows (genera) and columns (clinical groups). Color gradients were scaled by Z-score normalization within genus to emphasize relative abundance variations.

Differential analysis between the TRM and the fecal microbiota dataset
A total of seven publicly available cohorts from the NCBI database were aggregated as our reference fecal microbiota dataset. The potential batch effects across diverse public datasets and our own cohort were reduced by the adjust_batch function implemented in the R package MMUPHin (v1.4.2)58. We performed a detailed analysis to investigate the relationship between the cancer-associated microbial signals in the fecal and tissue datasets. We conducted parallel differential abundance analyses using ANCOM (R package ANCOMBC, v2.4.0)59 on both datasets to identify taxa significantly associated with CRC in each environment. To formally test the relationship between these signals, we generated a scatter plot of the p-values from the tissue analysis versus the p-values from the fecal analysis for these taxa. The Spearman correlation test was used to confirm this observation.

Multivariate analysis of microbiota and patient survival
To control for potential clinical confounders that may influence patient survival outcomes, we conducted a multivariate analysis. This analysis was adjusted for the following covariates: patient age, gender, tumor stage (I–IV), tumor location (left or right), differentiation grade (MA, PD, MD, HD), and chemotherapy status (Yes or No). The patients’ ages were divided into Young (≤60) and Senior (>60).

To evaluate the association between microbial features and survival outcomes, we primarily employed Cox proportional hazards models for each individual bacterial on genus-level, adjusting for all listed clinical covariates. As a complementary analysis, we also evaluated the 5-year survival outcome using multivariate logistic regression. The resulting p-values were adjusted by the BH method.

We then proceeded to a community-level analysis. First, a Lasso-Cox proportional hazards model (R package glmnet, v4.1.10) was employed for feature selection to identify a combination of microbial and clinical variables most strongly associated with survival. These selected features were then used to build a Cox proportional hazards model. The p-values and Hazard Ratios with their corresponding 95% confidence intervals were shown in a forest plot.

We also built a multivariate logistic regression model based on the features selected by Lasso for the 5-year survival outcome. The regression coefficients, standard errors, p-values, and Odds Ratios (OR) with their corresponding 95% confidence intervals were shown in a forest plot. Finally, a random forest model was employed for potential non-linear and high-dimensional effects. After partitioning the data into training (70%) and testing (30%) sets, we first built a baseline model using only demographic, clinical, and treatment variables. Another model was then constructed by adding the full microbial data to demonstrate its impact on the model.

Machine learning modeling workflow
Microbiome data preprocessing prior to machine learning modeling was performed using the metagenomeSeq package (v1.44.0)60. Low-prevalence features were first filtered by removing genera with occurrence frequency <5% across samples (prevalence threshold = 0.05) to reduce sparsity-induced noise. The retained microbial features were then normalized using cumulative sum scaling (CSS) normalization through three-stage processing:
1. Data structure conversion with the function newMRexperiment to create a metagenomic sequencing object.
2. Quantile normalization via the function cumNorm with scaling percentile parameter p=0.5 (median scaling).
3. Extraction of normalized counts using function MRcounts (norm=TRUE, log=FALSE).

This normalization strategy addresses compositional bias by scaling counts to the median cumulative sum of features up to the specified quantile60. Compared to rarefaction, CSS preserves quantitative relationships while mitigating sequencing depth artifacts, as demonstrated in sparse microbial datasets61. The final normalized matrix was used as input for downstream machine learning pipelines.

Random forest modeling based on the microbiome sequencing data at the family and genus levels was performed to address various classification tasks using default parameters. Initially, the cohorts were randomly divided into a discovery set (80%) and a validation set (20%). The discovery cohort comprised 182 normal, 68 polyp, and 568 CRC samples, whereas the validation cohort included 45 normal, 17 polyp, and 142 CRC samples. The discovery data were used for model training, the validation data for model validation, and an independent external cohort from another hospital served as the test set for final model evaluation.

For each classification problem, models were built using the caret62 and randomForest63 packages in R. The training set underwent three repetitions of 5-fold cross-validation to optimize the model performance, and the resultant models were subsequently validated on the test set. In the binary classification scenario, the AUC was computed from both the cross-validation results on the training set and from the independent testing set, with ROC curves being plotted for visual assessment. For the three-class classification problem involving normal tissue, polyp tissue, and tumor tissue, class-specific metrics—including accuracy, recall, and F1 score—were calculated on both training and testing sets using the MLmetrics package. Additionally, a confusion matrix detailing the model’s prediction outcomes on the test set was generated.

Minimal feature set selection
For pairwise modeling among normal tissue, polyp tissue, and cancer tissue, we aimed to select a limited number of key microbial genera that are both critical for the classification task and exhibit low inter-correlation while also being sufficiently prevalent in the samples. The sequential forward selection strategy for minimal feature identification was implemented as follows:
- Step 1 (Initial Ranking): An initial Random Forest model was trained on the entire training dataset, and all microbial features were ranked based on their Mean Decrease Accuracy (MDA) scores.
- Step 2 (Candidate Pool): The top K highest-ranking features were selected as a candidate pool for the subsequent selection process. K=10 works well in our study.
- Step 3 (Filtered Forward Selection): Starting from a null model, features were iteratively added one at a time from the candidate pool in descending order of importance, subject to the following quality filters: (i) Features classified as ‘unclassified’ were systematically excluded. (ii) To avoid collinearity and information redundancy, pairwise correlation (absolute value) with all selected features must be < 0.7. (iii) The feature must exhibit a prevalence (relative abundance > 0) in at least 25% of all samples64.
- Step 4 (Performance Evaluation): After adding each feature, a new model was trained. Model performance was evaluated using 5-fold cross-validation. Note that when the number of features included in the model is less than 2, we train a linear model. Otherwise, an RF model is trained. In summary, for the kth iteration (k=1, 2, ..., K), there are k features included in the model. We evaluated the performance of the model using the mean of training and validation AUC from 5-fold cross-validation, denoted by AUCtr(k) and AUCval(k), respectively.
- Step 5 (Criteria for Variable Selection): If there exists some k such that AUCtr(k)=1, define K*=smallest k such that AUCtr(k)=1, else define K*=K+1. To avoid overfitting, we use the 1-SE Rule65 to define the maximum number of variables that can be included in the model as
k*=smallest k such that AUCtr(k) ≥ AUCm tr - SE for k < K*
where AUCm tr is the maximum of {AUCtr(k) : k=1, 2..., K*-1} and SE is its standard error: With this information, we select the first k0 variables to build our prediction model with k0 defined by k0 = argmax k≤k* AUCval(k). That is, we select variables giving the best validation performance in cross-validation while controlling the risk of overfitting.

16S RNA FISH
Different gene probes (Table S3) labeled with different fluorophores were used to detect the bacterial abundance within human tissues by FISH. The specific steps for FISH staining are as follows:
1. Tissue fixation: After the tissue was rinsed with PBS, it was immediately placed in the in situ hybridization fixative for more than 12 h.
2. Dehydration wax embedding: After the tissue was fixed, the tissue block in the target area was cut into about 3 mm thick in the fume hood, dehydrated by low to high gradient alcohol, transparent in xylene, and embedded in wax.
3. Paraffin section: Paraffin blocks were sliced 5 μm thick by a slicing machine, sliced by a spreader, and sliced in a 60 °C oven for 2 h.
4. Paraffin sections dewaxed to water: The sections were placed in dewaxing transparent liquid I, 15 min - dewaxing transparent liquid II, 15 min - anhydrous ethanol I, 5 min - anhydrous ethanol II 5 min - 85% alcohol, 5 min - 75% alcohol, 5 min, DEPC water immersion.
5. Repair and digestion: According to the type of tissue, the length of fixed time, and the positioning of the index, different repair fluids were selected for repair. After natural cooling, the stroke circle was organized. The samples were incubated with 1 mg/mL lysozyme (Sigma #62970) at room temperature for 10 minutes. According to the characteristics of different tissues, protease K (20 μg/ml) (Servicebio #G1234) was added dropwise to digest at 37 °C. The digestion time is shown in the table above. After washing with pure water, PBS was washed 3 times, 5 min each time.
6. Pre-hybridization: The pre-hybridization solution (Servicebio #G3045) was added dropwise and incubated at 37 °C for 1 h.
7. Hybridization: Dump the pre-hybridization solution, drop the probe-containing hybridization solution, and mix overnight in the incubator at 40 °C.
8. Washing after hybridization: Wash the hybridization solution, 2× SSC (Servicebio #G3015), 37 °C for 10 min, 1× SSC, 37 °C for 2× 5 min, 0.5× SSC for 10 min at room temperature.
9. DAPI double staining nucleus: DAPI staining solution was added to the slices, incubated in the dark for 8 min, and the anti-fluorescence quenching sealing agent was added after washing.
10. Microscopic photography: Sections were observed under a Nikon upright fluorescence microscope, and images were collected. (UV excitation wavelength 330-380 nm, emission wavelength 420 nm, blue light. The green light excitation wavelength of FAM (488) is 465-495 nm, and the emission wavelength is 515-555 nm. Cy5 pink light excitation wavelength of 608-648 nm, emission wavelength of 672-712 nm, pink light).

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Dr. Xiang et al. analyzed 1,022 specimens—including normal mucosa, polyps, and colorectal cancer (CRC) tissues—to characterize tumor-resident microbiota (TRM) composition and abundance. They integrated microbial signatures into a machine learning framework to develop a diagnostic model that distinguishes normal, precancerous, and cancerous tissues. Their findings suggest that TRM alterations may be involved in the early stages of colorectal tumorigenesis rather than later disease progression.\nThis is a comprehensive microbiome analysis in CRC, using a large sample size and incorporating tumor tissue, normal tissue, and stool samples. The inclusion of multiple biospecimen types expands upon prior other fecal-based studies and could provide some additional insight into CRC pathogenesis.\n• While this study's integration of multiple sample types is potentially valuable, the manuscript does not clearly demonstrate any added clear impact and clinical utility in predictive modeling for normal, polyps and CRC, beyond current clinical methods that distinguish tumor, normal, and adenoma tissues.\n• There are several methodological concerns:\n- The use of 16S rRNA sequencing is outdated, with limited resolution for species-level identification and limited capacity for functional inference.\n- The analysis appears to be taxon-based and does not address the compositional nature of microbiome data. Compositionally aware methods, such as ANCOM, should be used to minimize statistical bias.\n- Multiple comparisons should be corrected using appropriate methods (e.g., false discovery rate adjustment).\n- For the machine learning model, further detail is needed on feature importance and the minimal feature set required for accurate multiclass classification (normal, polyp, tumor). Validation in an independent external cohort is also necessary.\n• Key methodological details are lacking. There is not enough detail provided in the methods for the work to be reproduced.\n- The manuscript does not sufficiently describe the study design, participant characteristics, inclusion and exclusion criteria, or treatment information.\n- It is unclear whether “normal” tissue refers to adjacent tissue from CRC patients or from individuals without CRC or adenomas.\n- Additionally, it is not specified whether fecal, normal, polyp, and tumor samples were collected from the same participants or from different individuals.\n- The process for defining discovery and validation cohorts is also unclear.\n• Due to the limited clinical data and the lack of clarity in predictive model performance, the conclusions are only partially supported. More evidence is needed to clarify the utility of each sample type and any association with clinical outcomes. The manuscript refers to a “fundamental temporal pattern in TRM dynamics,” but since the study is cross-sectional (not longitudinal), this language may be misleading. Terms like “temporal” and “dynamics” should be avoided unless causal inference over time is supported by the data.\n• The null associations with CRC survival are not well supported without demographic, clinical, treatment, and follow-up data. The methodological description—especially for sequencing processing and statistical analysis—falls short of current standards. In addition, references to temporal progression are misleading given the cross-sectional design. These issues require revision before the manuscript is suitable for publication.", "reviewer_2": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.\nI was able to clone the repository, install all dependencies, and execute the main analysis pipeline without issues.", "reviewer_3": "In this paper Haoyi Xiang and colleagues examine the distribution of tissue-resident microbiota in a large cohort of CRC, adenomas and healthy controls. The authors argue that there is a signal associated with case vs. control but not with tumor stage or clinical outcome. This is an important question, and the paper makes a solid contribution advancing our knowledge of this difficult to study problem.\nOverall, the paper is well written, and the results are compelling. This is a challenging study as the low microbial biomass in these samples makes it difficult to determine signal from noise. The authors have put a good deal of thought into quality control of these difficult to characterize samples. But the manuscript could provide additional details of the decisions that were made and better establish that their results cannot be explained by kit contaminant or other low abundance nuisance signals.\nThe authors used qPCR against universal 16S rRNA primers and state that samples containing 10^3 CFU/g (or presumably less) were excluded from downstream analysis (line 81). But why was this value chosen? Are there samples that are just above this threshold that may contain contaminant DNA that were still included? Were all the negative controls below this value? Why do the authors believe that samples above this threshold are not dominated by contaminant?\nLikewise, the authors state that “High-throughput sequencing was then performed with stringent contamination controls to ensure data reliability”. But what were the controls? How many? How did the authors convince themselves that the signals in their samples were real?\nI would encourage the authors to look at the presentation of negative controls in another cancer microbiome study (The human tumor microbiome is composed of tumor type–specific intracellular bacteria; https://www.science.org/doi/10.1126/science.aay9189). In Fig 1B of that paper, the data from negative controls are explicitly shown and are clearly different from the distributions from the cancer samples. Such a presentation that allows readers to see directly these differences is much more compelling than simply stating that controls have ensured adequate reliability.\nI appreciate the culture controls seen in Fig. 1C, but again more detail would make these results more compelling. Were these results only seen in a single sample representing each sample type? If there are multiple samples, can summary data from all samples be shown? Are there statistically significant differences between the different sample types? Were the negative controls in Fig. 1C and elsewhere in the paper taken at the beginning of the pipeline and exposed to every step the other samples were?\nThe fecal samples in Fig S1 seem like they are under-analyzed. I appreciate that these samples are distinct from tissue samples, but can they also be used to discriminate case from control? Are the fecal samples more sensitive to stage or patient outcome than the tissue samples? Is there any relationship between the fecal samples and tissue samples despite the PCOA convergence. If for example, the authors calculated a p-value vs. p-value plot for case control inference, would there be a relationship between the taxa significantly associated with cancer in the fecal samples and the tissue samples? A connection between the fecal samples and the tissue samples is potentially important and could be used to argue against the signal in the tissue samples arising from contaminant.\nOverall, the FISH data are very compelling and make a strong case that there are taxa associated with cancer that cannot just be explained by sequencing error or artifact.\nCould the authors directly correlate the FISH data shown in Figs. 5B, 5D and 5E with the relative abundance of the corresponding taxa by reporting an r-squared and p-values for the correlation of the two different types of data? This would be a way to determine how quantitative the different methods are and how much they are (or are not) in agreement with each other?\nI found the section on contamination correction somewhat confusing. Although this is presumably a copy and paste error, lines 558 and 562 refer to mouse samples so it is not clear if the contamination correction section applies to the data in this paper. There is a deep literature on contamination correction (with algorithms such as decontam and SCRuB) and it appears as if the approach the authors take in lines 562-582 mirror these other algorithms. The authors should explicitly state why they have not used these other algorithms and how their approach differs. All decontamination methods have serious limitations and the authors should consider presenting data with and without decontamination to establish that their results are not an artifact of a decontamination algorithm making inappropriate assumptions about the data. If on the other hand, the decontamination algorithm was not used on these data and this section is a copy and paste error, that should be clarified.", "reviewer_4": "This manuscript investigates tissue-resident microbiota (TRM) in colorectal cancer (CRC), polyps, and normal mucosa across a large cohort. The authors report progressive shifts in TRM during the adenoma–carcinoma sequence, stability across tumor stages, and strong classification accuracy using machine learning models. They conclude that TRM alterations may drive early tumorigenesis and provide diagnostic potential.\nThe study is ambitious, includes a large sample size, and employs rigorous contamination controls. However, I have significant concerns regarding the clinical annotation of samples, translational framing, and mechanistic depth. These issues limit both interpretability and the claimed clinical relevance.\n1. Clinical Annotation and Patient-Level Information: The manuscript provides minimal clinical data beyond histologic category (normal, polyp, CRC). This is a major deficiency. Not all polyps are biologically equivalent. Were all polyps adenomatous, or did the cohort include serrated/hyperplastic lesions? These follow distinct molecular pathways. Tumor location is highly relevant, as right- and left-sided CRC differ in microbiome associations, molecular profiles, and clinical outcomes. Similarly, MSI vs. MSS status has profound implications. Without this information, it is difficult to know whether observed TRM signatures reflect true biology or unmeasured heterogeneity. The authors should provide detailed patient metadata (location, MSI status, stage distribution, molecular alterations). At minimum, stratified analyses by right vs. left colon and MSI status are needed.\n2. Interpretation of TRM as Diagnostic Biomarkers: TRM profiling requires tissue acquisition via polypectomy, biopsy, or surgery. At that point, diagnosis is already established by pathology, limiting the utility of TRM as diagnostic biomarkers. While differences between TRM and luminal microbiota are scientifically interesting, their translational application is overstated. TRM cannot realistically serve as a screening tool in their current form due to the invasive nature of tissue acquisition. The authors should clarify whether they view TRM analysis as a mechanistic research tool, an adjunct for risk stratification in already-resected polyps, or a step toward identifying non-invasive correlates (stool, blood) of tissue signatures. Without this clarification, the proposed “diagnostic framework” appears conceptually flawed.\n3. Reported AUC values are extremely high, in some cases approaching 1.0. This raises concerns about model overfitting, particularly given the relatively small number of polyp samples. Validation was performed within a single institution. External, geographically distinct cohorts are needed to establish generalizability.\n4. The study is largely descriptive. While the FISH validation is a strength, there is little investigation of functional pathways or host-microbe interactions that could explain how TRM influence tumor initiation. Integration with host transcriptomic, immune, or metabolomic data would strengthen the manuscript substantially. At minimum, the discussion should acknowledge this gap and frame the work as hypothesis-generating.\n5. The observation that TRM differences are associated with tumor initiation but not progression is important. However, the manuscript repeats this point extensively across sections without expanding on its biological implications. A more balanced discussion should highlight both strengths and limitations, including the ongoing controversy regarding the existence and significance of intratumoral microbiota in low-biomass samples."}, "author_response": "Response to reviewers for NCOMMS-25-46758\nDear Reviewers,\nWe sincerely appreciate your insightful and constructive comments, which have helped us significantly improve the manuscript. In response to your suggestions, we have carefully revised the paper, incorporating additional experiments, new data, and enhanced data analysis.\nBelow, we provide a summary of the major revisions:\n1. Clarified study focus and narrative framework: We reframed the overall narrative to emphasize the primary aim of our study is to delineate community-level shifts and compositional remodeling of TRM along the normal-polyp-CRC sequence. Our focus lies in providing a mechanistic foundation for future studies rather than to develop or propose TRM profiling as a diagnostic tool.\n2. External cohort validation and enhanced clinical context: We incorporated an independent external cohort to further assess the robustness and reproducibility of our findings (Figure 4E-F), and we now provide demographic and clinical metadata for both cohorts (Tables S1 and S2). We additionally performed multivariate and stratified analyses across clinical subgroups, which consistently supported our original conclusions (Figures S3, S9).\n3. Comparison of TRM and fecal microbiota: We included comparative analyses between fecal microbiota and TRM, which provided additional support for the concept of TRM as a distinct biological system.\n4. Refined sample screening and decontamination criteria: We expanded methodological clarity regarding threshold selection for sample quality control and decontamination procedures (Figure S12), enhancing transparency and rigor.\n5. Strengthened experimental validation: We added further validation analyses, including FISH-sequencing concordance (Figure S11), additional tissue plating experiments, and inter-group comparisons, that reinforce the reliability of the observed TRM patterns.\n6. Methodological enhancements: We refined the Methods section to provide clearer procedural details and adopted more appropriate statistical approaches in response to your comments.\nIn summary, we added six new supplementary figures (Figures S1, S3, S7, S9, S11, S12), revised Figures 1 and 4, expanded the analyses of fecal microbiota, decontamination workflows, tissue plating, and FISH validation, and incorporated external cohort testing as well as modeling methodology. We have provided all the code for analyses at https://github.com/Lingning927/microbioInCRC.\nPlease find the point-to-point responses to the comments below. For clarity, our responses are highlighted in blue. We sincerely hope that the revisions and improvements we have made address your concerns and contribute to the overall quality of the manuscript. We eagerly await your feedback and further suggestions.\nThank you once again for your invaluable input and guidance.\nBest regards,\nZhangfa Song, M.D.\nReviewer #1: ... [Author response continues for all reviewers] ..."}, {"round": 2, "reviewer_comments": {"reviewer_1": "Dr. Xiang et al. analyzed 1,134 specimens—including normal mucosa, polyps, and CRC tissues—to characterize tumor-resident microbiota (TRM) composition and developed a machine learning model to distinguish normal, precancerous, and cancerous tissues. Their findings suggest that TRM shifts may occur early in colorectal tumorigenesis. The study benefits from a large sample size and inclusion of multiple biospecimen types, expanding upon prior fecal-only studies. \nDespite some revisions, several major issues remain unaddressed.\n1. Unclear Clinical Impact of the Predictive Model:\nThe clinical utility of distinguishing normal tissue, adenomas, and cancer using TRM profiles is unclear, as these distinctions are routinely made during colonoscopy. The authors state that their TRM-based signature could serve as an adjunct to histopathology, yet the predictive model does not incorporate current histopathological factors. The practical clinical implications and potential added value of the model should be explicitly stated in the manuscript.\n2. Unclear Source of Normal, Polyp, and Tumor Samples:\nThe authors state in the revision (Line 427) that “normal, polyp, and tumor samples were collected from different individuals.” However, the Methods indicate that normal intestinal tissue was collected at least 5 cm from a CRC lesion, implying that some “normal” samples were taken from cancer patients. This directly contradicts the claim that samples were obtained from different individuals. It remains unclear:\n• Whether “normal” tissues were from cancer patients or cancer-free individuals;\n• Whether polyp samples were taken from the same or different individuals.\nThis ambiguity raises concerns about sample classification and the validity of downstream analyses.\n3. Incomplete Application of False Discovery Rate (FDR) Correction:\nThe authors note that FDR-adjusted results are shown in Figure 1F–G; however, all taxon-level analyses require FDR correction given the large number of comparisons. It is unclear whether FDR adjustment was applied to other figures or analyses. In the Results, none of the taxon-based findings report FDR-adjusted values. This is a substantial issue, as the reported differences among tumor, adenoma, and normal tissues do not remain significant after FDR correction.\nThe authors should:\n• Apply FDR correction to all taxon-level comparisons;\n• Clearly distinguish p-values from q-values in every figure and legend;\n• Report which findings remain significant after correction.\n4. Insufficient Detail on Prediction Model Development:\nThe revised description of the predictive model remains inadequate. The authors should clarify:\n• How cutoffs were selected;\n• Which baseline clinical and demographic characteristics were included;\n• The criteria used for variable selection.\nWithout this information, the rigor and reproducibility of the model cannot be evaluated.\n5. Stool Samples Not Addressed\nThe origin and characteristics of stool samples are not described. It is unclear who provided them, how they were collected, and how they relate to the tissue samples. These details are essential for interpreting the microbiome results.\nOverall Assessment: Although the study includes a large number of samples, concerns about methodological detail, sample classification, and overall study quality remain.", "reviewer_2": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.", "reviewer_3": "The authors have done a thorough job of addressing reviewer concerns and I have no further comments.", "reviewer_4": "The authors have been highly responsive and have substantially strengthened the manuscript. The revisions directly and fully address the concerns raised in my initial review. In particular, the expanded methodological detail, clearer cohort definitions, improved description of contamination controls, incorporation of an external validation cohort, clarification of statistical procedures, and removal of misleading temporal language all materially improve the rigor and clarity of the study. The additional FISH validation further increases confidence in the biological relevance of the identified TRM signatures.\nOverall, the manuscript is now significantly improved. I find the responses and revisions fully satisfactory, and I have no remaining concerns."}, "author_response": "Response to reviewers for NCOMMS-25-46758A\nDear Reviewers,\nWe sincerely appreciate your insightful and constructive comments, which have helped us significantly improve the manuscript. In response to your suggestions, we have carefully revised the paper, and enhanced data analysis. Please find the point-to-point responses to the comments below. For clarity, our responses are highlighted in blue. We sincerely hope that the revisions and improvements we have made address your concerns and contribute to the overall quality of the manuscript. We eagerly await your feedback and further suggestions.\nThank you once again for your invaluable input and guidance.\nBest regards,\nZhangfa Song, M.D.\nReviewer #1: ... [Author response continues] ..."}, {"round": 3, "reviewer_comments": {"reviewer_1": "The authors have addressed reviewer concerns and I have no further comments.", "reviewer_2": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.", "reviewer_3": "Overall Assessment\nWhile the authors have addressed some earlier comments, the revised manuscript introduces\nseveral new methodological, interpretational, and conceptual issues that materially affect the\nrobustness of the results. Many of the central claims—particularly those related to microbial\nprogression, “transitional states,” predictive modeling, and survival associations—remain\nunsupported given the cross-sectional design, heterogeneous sampling contexts, and analytical\nlimitations. Minor revisions to the language are required to align the statistical methods,\ninterpretation, and wording with what can be validly inferred from the data.\n1. Study design terminology (Abstract vs. Results)\nThe Abstract correctly refers to the study as cross-sectional, but the Results (line 74) describe it\nas a retrospective cohort study, which is misleading. The wording should be made consistent\nand should reflect the cross-sectional nature of the analyses. All microbial measurements in this\nwork are cross-sectional and obtained at a single time point per patient. The normal, polyp, and\nCRC samples come from different individuals, collected by different procedures (endoscopic vs\nsurgical), and do not represent sequential stages within the same subjects. Therefore, the\ndataset does not capture tumor initiation, temporal progression, or within-person transitions.\nEven if survival analyses were performed on the same patients, these relate only to post-\ndiagnosis outcomes and do not provide information about when microbial alterations occurred.\nSurvival associations cannot be used to infer “early” (line 38) versus “late” microbial changes.\nIn the absence of longitudinal data, repeated sampling, pre-diagnostic follow-up, the authors\ncannot ascribe their findings to “early tumorigenesis,” (line 39) “initial phases,” (line 37)\n“progression,” “stage-specific transitions,” (line 156) or any other temporal or causal process. All\nsuch language constitutes overinterpretation and should be removed or replaced with strictly\ncross-sectional wording (e.g., “group differences”).\n2. Fig S1A: Batch effects in external fecal datasets\nFor the fecal datasets analyzed in Figure S1A, it is unclear how batch effects were addressed or\ndistinguished from biological differences. Public metagenomic data vary substantially in\nextraction protocols, sequencing platforms, and QC standards. Please clarify how batch\nvariation was controlled or whether the observed differences may simply reflect technical\nartifacts.\n3. Misinterpretation of “paired analyses” (Lines 115–118)\nThe statement that “this pattern persisted in paired analyses of tumor and adjacent normal\ntissues” is unclear. Stratifying paired tumor–normal comparisons by age, gender, TNM stage, or\nother patient-level covariates is not possible, as both samples in a pair share identical metadata.\nPlease clarify what analysis was actually performed or remove language implying stratified\npaired analyses that are not statistically meaningful.\n4. Machine learning analyses and interpretation (Line 142)\nBecause no microbial features were significantly associated with clinical parameters in\nunivariate analyses, forcing these same features into Random Forest models does not “evaluate\npredictive value,” but instead reanalyzes noise with a more complex method. The claim that\npaired analyses demonstrate that “inter-individual variation did not drive the results” is also\nunsupported by the modeling framework. The ML section should be framed explicitly as\nexploratory, with appropriate discussion of its limitations.\n5. Survival analysis is methodologically inappropriate (Lines 151–155; Fig S3)\nThe survival analysis uses logistic regression, which does not model time-to-event outcomes,\nignores censoring, and produces odds ratios instead of hazard ratios. This approach is not valid\nfor survival data. Additionally, the covariates included in the models are not specified. A proper\nCox proportional hazards model (or alternative survival model) with reporting of events,\ncovariates, effect sizes, and hazard ratios is required.\n6. Overstated conclusions regarding survival (Line 156)\nThe statement that “TRM composition is largely independent of survival outcomes and other\nclinical features” is too strong. A lack of statistically significant findings in underpowered or\ninappropriate models does not demonstrate independence. This conclusion should be\nweakened to reflect the analytical limitations and absence of robust survival modeling.\n7. Overinterpretation of cross-sectional data as progression (Lines 168–182)\nThe interpretation in this section overstates what can be inferred from cross-sectional, between-\nindividual comparisons. Terms such as “three distinct microbial ecologies,” “clear transitional\nstate,” and “stage-specific microbial changes” imply longitudinal processes or causal\nprogression that the study design cannot support. The use of a Cochran–Armitage trend test\nfurther reinforces an unsupported progression model. These results should be reframed as\ncross-sectional group differences, not developmental transitions or ordered biological stages.\n8. Confounding from heterogeneous sampling methods\nNormal and CRC tissues are surgical specimens, whereas polyps are typically collected via\ncolonoscopy. Classification models may therefore reflect sampling modality or batch effects\nrather than biological signal. This must be explicitly tested (e.g., including sampling modality as\na covariate or performing batch-corrected analyses).\n9. Misinterpretation of polyp vs Stage I CRC comparison\nThe polyp vs Stage I CRC comparison does not isolate “early tumorigenesis.” Both groups\ncome from different individuals and different sampling contexts, and no longitudinal data support\na developmental sequence. Claims suggesting temporal or mechanistic insight should be\nremoved.", "reviewer_4": "Thanks to the authors for addressing my main concerns with evidence from external validation and the mouse model. My remaining concerns are as follows: The manuscript still narrates a “progressive transition” of the microbiome from normal → polyp → cancer, but all samples are cross-sectional and mostly from different people in your response. Please consistently rephrase throughout Results, not just Discussion, to avoid implying longitudinal evolution or causality. It seems there is an internal inconsistency between the description of the study design and the sample source. In line 427, the revised text states that “Normal, polyp, and tumor samples were collected from different individuals,” implying these are independent samples from distinct participants. However, in line 74 of the Results section, the manuscript says “we conducted a retrospective cohort study utilizing matched CRC and normal colorectal tissue,” which indicates paired samples (tumor and normal) from the same individuals. You need to explicitly show that “normal,” “polyp,” and “CRC” groups are not just clinically different populations sampled in different contexts, and please clarify the sample design. Report and adjust for sampling method, anatomic site, and other clinical/technical covariates when claiming a “progressive” gradient. Saying “TRM is critical in initiation but not progression” is stronger than what the data can prove. You can say “TRM differences are most apparent between diagnostic classes (normal / polyp / CRC), not across later-stage clinical strata or survival.” Determine whether any per-patient paired samples are present in both the training and validation folds, and rerun the models with leave-one-patient-out if necessary to show generalizability. Highlight the external cohort AUC (~0.78–0.80) as the realistic performance."}, "author_response": "Response to reviewer for NCOMMS-25-46758B\nDear Reviewer,\nWe sincerely appreciate your insightful and constructive comments, which have helped us significantly improve the manuscript. In response to your suggestions, we have carefully revised the paper, and enhanced data analysis. Please find the point-to-point responses to the comments below. For clarity, our responses are highlighted in blue. We sincerely hope that the revisions and improvements we have made address your concerns and contribute to the overall quality of the manuscript. We eagerly await your feedback and further suggestions.\nThank you once again for your invaluable input and guidance.\nBest regards,\nZhangfa Song, M.D.\nWhile the authors have addressed some earlier comments, the revised manuscript introduces several new methodological, interpretational, and conceptual issues that materially affect the robustness of the results. Many of the central claims—particularly those related to microbial progression, “transitional states,” predictive modeling, and survival associations—remain unsupported given the cross-sectional design, heterogeneous sampling contexts, and analytical limitations. Minor revisions to the language are required to align the statistical methods, interpretation, and wording with what can be validly inferred from the data. ... [Author response continues]"}]

Decision letter
Overall Assessment
While the authors have addressed some earlier comments, the revised manuscript introduces
several new methodological, interpretational, and conceptual issues that materially affect the
robustness of the results. Many of the central claims—particularly those related to microbial
progression, “transitional states,” predictive modeling, and survival associations—remain
unsupported given the cross-sectional design, heterogeneous sampling contexts, and analytical
limitations. Minor revisions to the language are required to align the statistical methods,
interpretation, and wording with what can be validly inferred from the data.
1. Study design terminology (Abstract vs. Results)
The Abstract correctly refers to the study as cross-sectional, but the Results (line 74) describe it
as a retrospective cohort study, which is misleading. The wording should be made consistent
and should reflect the cross-sectional nature of the analyses. All microbial measurements in this
work are cross-sectional and obtained at a single time point per patient. The normal, polyp, and
CRC samples come from different individuals, collected by different procedures (endoscopic vs
surgical), and do not represent sequential stages within the same subjects. Therefore, the
dataset does not capture tumor initiation, temporal progression, or within-person transitions.
Even if survival analyses were performed on the same patients, these relate only to post-
diagnosis outcomes and do not provide information about when microbial alterations occurred.
Survival associations cannot be used to infer “early” (line 38) versus “late” microbial changes.
In the absence of longitudinal data, repeated sampling, pre-diagnostic follow-up, the authors
cannot ascribe their findings to “early tumorigenesis,” (line 39) “initial phases,” (line 37)
“progression,” “stage-specific transitions,” (line 156) or any other temporal or causal process. All
such language constitutes overinterpretation and should be removed or replaced with strictly
cross-sectional wording (e.g., “group differences”).
2. Fig S1A: Batch effects in external fecal datasets
For the fecal datasets analyzed in Figure S1A, it is unclear how batch effects were addressed or
distinguished from biological differences. Public metagenomic data vary substantially in
extraction protocols, sequencing platforms, and QC standards. Please clarify how batch
variation was controlled or whether the observed differences may simply reflect technical
artifacts.
3. Misinterpretation of “paired analyses” (Lines 115–118)
The statement that “this pattern persisted in paired analyses of tumor and adjacent normal
tissues” is unclear. Stratifying paired tumor–normal comparisons by age, gender, TNM stage, or
other patient-level covariates is not possible, as both samples in a pair share identical metadata.
Please clarify what analysis was actually performed or remove language implying stratified
paired analyses that are not statistically meaningful.
4. Machine learning analyses and interpretation (Line 142)
Because no microbial features were significantly associated with clinical parameters in
univariate analyses, forcing these same features into Random Forest models does not
“evaluate predictive value,” but instead reanalyzes noise with a more complex method. The claim that
paired analyses demonstrate that “inter-individual variation did not drive the results” is also
unsupported by the modeling framework. The ML section should be framed explicitly as
exploratory, with appropriate discussion of its limitations.
5. Survival analysis is methodologically inappropriate (Lines 151–155; Fig S3)
The survival analysis uses logistic regression, which does not model time-to-event outcomes,
ignores censoring, and produces odds ratios instead of hazard ratios. This approach is not valid
for survival data. Additionally, the covariates included in the models are not specified. A proper
Cox proportional hazards model (or alternative survival model) with reporting of events,
covariates, effect sizes, and hazard ratios is required.
6. Overstated conclusions regarding survival (Line 156)
The statement that “TRM composition is largely independent of survival outcomes and other
clinical features” is too strong. A lack of statistically significant findings in underpowered or
inappropriate models does not demonstrate independence. This conclusion should be weakened to
reflect the analytical limitations and absence of robust survival modeling.
7. Overinterpretation of cross-sectional data as progression (Lines 168–182)
The interpretation in this section overstates what can be inferred from cross-sectional,
between-individual comparisons. Terms such as “three distinct microbial ecologies,”
“clear transitional state,” and “stage-specific microbial changes” imply longitudinal processes or causal
progression that the study design cannot support. The use of a Cochran–Armitage trend test
further reinforces an unsupported progression model. These results should be reframed as
cross-sectional group differences, not developmental transitions or ordered biological stages.
8. Confounding from heterogeneous sampling methods
Normal and CRC tissues are surgical specimens, whereas polyps are typically collected via
colonoscopy. Classification models may therefore reflect sampling modality or batch effects
rather than biological signal. This must be explicitly tested (e.g., including sampling modality as
a covariate or performing batch-corrected analyses).
9. Misinterpretation of polyp vs Stage I CRC comparison
The polyp vs Stage I CRC comparison does not isolate “early tumorigenesis.” Both groups
come from different individuals and different sampling contexts, and no longitudinal data support
a developmental sequence. Claims suggesting temporal or mechanistic insight should be removed.