The recruitment of CC patients was conducted at the Instituto Jalisciense de Cancerología in Guadalajara, Jalisco, Mexico. Clinically healthy female donors (HD) from the community participated as the control group. This study was performed following the ethical principles outlined in the Declaration of Helsinki (2024 revision) and was approved by the Research and Ethics Committees of the health institution (PRO-72/23) as well as by the University Center (22-92-CI-00323). All participants were informed about the objectives of the study, and written informed consent was obtained before their inclusion.

This cross-sectional study included 77 female participants: 49 patients diagnosed with cervical cancer (CC) and 28 age-adjusted healthy donors (HD). From each participant, fecal samples were collected for gut microbiota profiling, and peripheral blood was obtained for NK cell phenotypic analysis. CC patients were stratified into two groups: 27 treatment-naïve individuals with newly diagnosed CC, assigned as CC pre-treatment group, and 22 patients who had completed standard-of-care radio-chemotherapy (RCT), which consisted of 50 Gy delivered in 25 fractions with concomitant platinum-based chemotherapy, within the previous two weeks, assigned as CC post-treatment. The control group comprised 28 healthy women with no history of malignancy, confirmed by a negative Papanicolaou test within the past year. Written informed consent was obtained from all participants prior to enrollment.

Inclusion criteria for the CC group were: (i) age ≥18 years, (ii) histopathologically confirmed cervical cancer (stages I–IV, classified by TNM system), and (iii) provision of informed consent. Exclusion criteria for all participants included: (i) use of prebiotics or probiotics within four weeks before recruitment, (ii) any chronic infection other than HPV, (iii) hospitalization within the past three months due to COVID-19-related illness, (iv) diagnosed gastrointestinal disorders, (v) autoimmune diseases and (vi) pregnancy.

Fecal samples were collected and immediately stored at −80°C. DNA was extracted from 150 mg of frozen feces with Quick-DNA™ Fecal/Soil Microbe Miniprep (Zymo Research, USA, cat: D6010) according to the manufacturer’s protocol. DNA was quantified using a NanoDrop™ OneC spectrophotometer (Thermo Scientific, Waltham, MA, USA). The 16S metagenomic sequencing library preparation was performed according to the Illumina MiSeq System protocol (Illumina, San Diego, CA, USA) (23). V3 and V4 regions from 16S were amplified with Platinum Taq DNA Polymerase High fidelity (Invitrogen, Waltham, MA, USA) using primers with adaptors. The sequence of the primers used was: Forward: (5’TCGTCGGCAGCGTCAGATGTGTATAAGAGACAGCCTACGGGNGGCWGCAG-3’), reverse: (5’GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAGGACTACHVGGGTATCTAATCC- 3’). PCR conditions were performed according to the Illumina protocol. Product purification was achieved using AMPure XP^® (Beckman Coulter, Indianapolis, IN, USA) magnetic beads and was quantified with the Qubit^® 3 dsDNA HS kit (Invitrogen, Waltham, MA, USA) according to product indications. Next, index incorporation was achieved with the Nextera XT Index Kit v2 Set A (No. Cat. FC-131-2001, Illumina, San Diego, CA, USA) by a second PCR amplification. Finally, amplicons were pooled to equimolar concentrations into a 4 nmol/L solution tube, which were then denatured, and was further diluted to the recommended loading concentration for the MiSeq Sample Loading (kit Miseq Reagent V3 600-cycle, Illumina, San Diego, CA, USA). The sequencing was performed according to the manufacturer´s protocol.

Analysis of 16S rRNA (V3-V4) sequences was performed using QIIME2 version 2024.2 amplicon distribution (24). Previously, the sequences whose Phred score were higher than 30 were processed. Then, raw reads (at least 100, 000 raw reads per sample) were further denoised using DADA2 via q2-dada2 (25) at default settings. Analyses were conducted with an average of 40, 000 denoised amplicon sequence variants (ASVs) per sample. Taxonomy assignation of our sequences was performed using a full-length 16S trained classifier (26), further employing Silva 138.1 as a reference taxonomic database (27, 28). ASVs identified as mitochondria and chloroplasts were removed. Then, filtered ASVs were aligned using multiple alignment using fast Fourier transform (MAFFT) via q2-alignment, and phylogeny was built with FastTree2 via q2-phylogeny (29). A-diversity indices (30–33) and β-diversity distances (Bray-Curtis, Jaccard, unweighted and weighted Unifrac) (34, 35) were computed via q2-diversity. Principal coordinate analysis (PCoA) plots were generated to visualize β-diversity distances using Emperor via q2-emperor (36) and the different distances were further analyzed using permutational multivariate analysis of variance (PERMANOVA) tests.

Differential abundance analyses at the genus level were performed using analysis of compositions of microbiomes with bias correction (ANCOM-BC) via q2-composition (37). Before analysis, a frequency filter was applied in which features that appeared more than 50 times in at least 10% of the samples were retained. A q ≤ 0.05 cut-off was used to assess significance, and a log fold change (LFC) ≥|1.0| to evaluate the effect size. To assess the potential metabolic profile of the gut microbiota, phylogenetic investigation of communities by reconstruction of unobserved states 2 (PICRUSt2) pipeline (38–42) was employed, coupled with the MetaCyc Database (43). The resulting pathways were further analyzed using ANCOM-BC, using the previously described parameters. Different taxonomic ratios were calculated following previously published methods (44).

Peripheral blood samples were collected using 10 mL K2-EDTA spray-coated tubes (Becton Dickinson, Franklin Lakes, New Jersey, USA; 366643) for the separation of peripheral blood mononuclear cells (PBMCs). PBMCs were isolated via density gradient centrifugation using Lymphoprep (Stemcell Technologies, Vancouver, British Columbia, Canada; 07851).

After isolation, PBMCs were stored in liquid nitrogen at −196°C until further analysis. Cell viability was assessed using trypan blue exclusion, and only samples exhibiting ≥90% viability were included in the following analyses.

A multi-parametric flow cytometry panel was employed to evaluate the expression of immune checkpoint molecules (PD-1, TIGIT, NKG2A, BTLA, TIM-3, and LAG-3) on NK cell subsets. Viability was evaluated in every sample in an independent staining using Zombie NIR™, and every sample showed a viability higher than 90%. The following monoclonal antibodies were utilized for staining 5 × 10⁵ PBMCs: Anti-CD3-FITC [fluorescein isothiocyanate, 300406], Anti-CD56-BV711 [Brilliant Violet 711, 362542], Anti-CD16-BV605 [Brilliant Violet 605, 302040], anti-CD45-AF700 [Alexa Fluor 700, 304024), Anti-TIGIT- PE [phycoerythrin, 372704], Anti-TIM-3-BV510 [Brilliant Violet 510, 345030], Anti-PD-1-BV421 [Brilliant Violet 421, 329920], Anti-LAG-3-PerCP/Cy5.5 [peridinin chlorophyll/cyanine 5.·5, 369312], Anti-BTLA-PE/Cy7 [phycoerythrin/cyanine 7, 344516], Anti-NKG2A-APC [allophycocyanin, 375108]. All antibodies were sourced from BioLegend (San Diego, CA, USA). Data acquisition was conducted on an Attune™ NxT Flow Cytometer (Thermo Fisher Scientific, Waltham, MA, USA). Compensation was made with compensation beads (Becton Dickinson, Franklin Lakes, NJ, USA; 55284). For each sample, singlet cells were identified using forward scatter-area (FSC-A) vs. forward scatter-height (FSC-H) dot plots, followed by gating lymphocytes using FSC-A vs. side scatter-area (SSC-A) plots. A total of 250, 000 events were recorded within the lymphocyte gate. Single and co-expression of immune checkpoint receptors was analyzed in peripheral NK cell populations using Kaluza software (version 2.1, Beckman Coulter, Brea, CA, USA). PBMCs were analyzed as CD3^-CD56^dim and CD3^-CD56^bright NK cell populations.

To standardize NK cell phenotyping and provide an integrative measure of exhaustion, we developed a composite Global NK Exhaustion Score. For each inhibitory receptor (PD-1, LAG-3, BTLA, TIM-3, TIGIT, NKG2A), Z-scores were calculated relative to the healthy donor distribution. The calculated Z-scores were then classified by specified percentiles, as explained in Supplementary Material 1 , to obtain a score (1 to 3). These values were summed separately for CD56^dim and CD56^bright NK subsets. The Global NK Exhaustion Score for each patient was defined as the mean of the two subset-specific scores (CD56^dim and CD56^bright), yielding a composite index of inhibitory receptor burden. This approach was inspired by previous work in T cells, where transcriptomic-based exhaustion scores were generated by integrating multiple exhaustion-related features into a unified index (45). Conceptually, our strategy also reflects the notion that NK dysfunction cannot be defined by a single marker but rather emerges from the combined expression of multiple inhibitory receptors and phenotypic alterations (46). Together, this percentile-based score provides a standardized framework to compare NK exhaustion across patients and to correlate it with microbiota alterations.

The microbiota dysbiosis score included α-diversity metrics (Shannon, Pielou, Simpson, and Strong indices) and the centered log-ratio (CLR)-transformed relative abundances of bacterial taxa previously identified as expanded or depleted in CC patients by ANCOM-BC analyses. After Z-score calculation, each parameter was scored using a similar percentile-based classification scheme, capturing deviations in either direction from control values, as both increases and decreases may reflect dysbiosis. Finally, correlation analysis between NK exhaustion and microbiota dysbiosis scores was performed. Detailed formulas and score calculation workflows are provided in Supplementary Material 1 .

Data distribution was assessed for normality using the D’Agostino–Pearson normality test. For comparisons between two groups, parametric data were analyzed using the Student’s T-test, while non-parametric data were analyzed using the Mann–Whitney U test. For comparisons involving three or more groups, a one-way analysis of variance (ANOVA) was applied for parametric data, with post-hoc comparisons adjusted using the Benjamini-Hochberg false discovery rate (FDR) method. The Kruskal-Wallis test with the FDR method of Benjamini-Hochberg for multiple comparisons was used for non-parametric data.

Microbiota composition and diversity analyses were conducted using non-parametric tests (Kruskal-Wallis) within the QIIME2 package. Correlation between microbiota and NK cell exhaustion markers was assessed using Spearman’s rank correlation coefficient.

Sociodemographic data were analyzed using the Chi-square test and are presented as frequencies and percentages. Statistical analyses were performed using GraphPad Prism, version 10.4. and R Studio version 4.4.2. P-values ≤ 0.05 were considered statistically significant.

A Random Forest (RF)-based classification model was developed to predict CC using microbiota composition, NK cell exhaustion markers, and patient demographic data. The dataset consisted of 104 variables, including microbial genera abundances, diversity indices, NK-cell receptor expression profiles, along with clinical variables such as age and BMI. Only HD and newly diagnosed, pre-treatment CC patients were considered for analysis. This approach was adapted from Tsakmaklis et al., 2023 (47). Briefly, feature selection was performed using Recursive Feature Elimination (RFE) with leave-group-out (Monte-Carlo) cross-validation (1000 iterations) to minimize overfitting and retain the 20 most predictive variables. The final RF model was trained with 200 decision trees (ntree=200) and an empirically determined mtry of 2. Model performance was assessed through repeated random partitioning (1000 iterations) into training (70%) and test (30%) subsets. In each iteration, the model was independently trained and tested, and its predictive ability was evaluated using the area under the receiver operating characteristic curve (ROC-AUC), reported as the mean and standard deviation across iterations reported as the mean and standard deviation (SD) across iterations. The SD represents the variability of model performance across these repeated cross-validation runs, providing an estimate of the model’s stability and generalizability. Classification accuracy was assessed using a confusion matrix, calculating sensitivity, specificity, and overall accuracy. To further validate the model and reduce predictor variables, stepwise logistic regression was applied to the most important features identified through RF analysis. The final logistic regression model was developed with the selected variables and assessed for its predictive accuracy.

To assess the impact of microbiota and NK cell exhaustion markers on patient survival, an independent mortality prediction model was developed, which was restricted to pre-treatment CC patients using a 106-variable data set. The XGBoost (Extreme Gradient Boosting) model demonstrated superior performance to RF and was selected as the final model, which was trained using 10-fold cross-validation with hyperparameter optimization through grid search. Performance was assessed through 1000 iterations of random partitioning (70% training, 30% testing). Similarly, stepwise logistic regression was applied to the most important predictors identified by XGBoost to refine the model and assess its predictive power.

For each selected variable, Kaplan-Meier survival curves were generated to assess the differences in survival between groups. The cutoff value for categorizing patients into “high” and “low” groups was determined by the point of maximum sensitivity and specificity on a ROC curve for each variable. This approach allowed a distinction between patients who survived and those who did not. Subsequently, Kaplan-Meier curves were constructed for survival analysis of the CC pre-treatment patients over the 15-month follow-up period.

The Clinicopathological characteristics of the study participants are presented in Table 1 . No significant differences were observed in the age across groups (p = 0.9658). However, BMI was significantly higher in CC pre-treatment patients compared to HD (p = 0.0238), with a greater prevalence of patients classified as overweight or with obesity. The most common histological type was squamous cell carcinoma, accounting for most cases in both the pre- and post-treatment groups. Stool consistency, assessed by the Bristol Scale, shifted from a median of 4 (IQR: 3–4) in HD to 3 (IQR: 2–6) in CC pre-treatment and 5 (IQR: 4–6) in post-treatment patients, consistent with bowel habit alterations linked to disease and treatment. All post-treatment patients underwent a standard regimen of fractionated radio-chemotherapy (RCT), receiving a total dose of 50 Grays in 25 sessions, concomitant with platinum-based chemotherapy.

The α-diversity of the gut microbiota was evaluated using Shannon, Simpson, Pielou evenness, Simpson evenness, and Strong metrics ( Figure 1A ). All comparisons showed significant differences between groups ( Supplementary Table 2A ). CC patients exhibited a marked reduction in overall diversity compared with HD, which was further aggravated after RCT. This reduction reflected both a loss of diversity (Shannon, Simpson) and evenness (Pielou, Simpson evenness), alongside increased bacterial dominance (Strong). This shift indicates a profound imbalance in gut microbiota, characterized by the overrepresentation of specific taxa that may contribute to chronic inflammation and dysbiosis. RCT further exacerbated this disruption, amplifying bacterial dominance, intensifying the loss of microbial diversity and homeostasis in CC patients.

β-diversity was assessed using Bray–Curtis ( Figure 1B ), Jaccard, weighted, and unweighted UniFrac metrics ( Supplementary Figure 2B ), and visualized through Principal Coordinate Analysis (PCoA). HD samples formed a relatively tight cluster, while pre-treatment CC samples displayed a broader distribution, indicating greater inter-individual variability. In contrast, post-treatment CC samples were more distantly separated from HD and clustered more compactly, reflecting a pronounced shift in microbial community structure. PERMANOVA confirmed significant differences between groups ( Supplementary Table 2C ), underscoring the strong disruption of gut microbiota composition in CC patients, which was further exacerbated following RCT.

Additionally, we assessed gut microbiota alpha and beta diversity in CC patients stratified by clinical stage, both before and after treatment. However, no significant compositional differences were observed across stages, and further analyses of microbiota-related parameters were not pursued (not shown).

At the phylum level, pre-treatment CC patients displayed a reduction in the abundance of Firmicutes (80.35% vs. 61.66%), an enrichment of Bacteroidota (16.53% vs. 31.72%) and Proteobacteria (0.90% vs. 3.04%) in comparison to HD. This shift was even more pronounced in post-treatment patients (Firmicutes 52.92%, Bacteroidota 34.48%, Proteobacteria 8.52%), with a further expansion of Actinobacteriota (1.08% vs. 2.91%) ( Figure 2A ).

At the genus level, pre-treatment CC patients showed an expansion of pathobionts, including Prevotella (3.72% vs. 10.55%) and Escherichia-Shigella (0.29% vs. 3.29%), alongside a decrease in SCFA-producing bacteria such as Faecalibacterium (15.53% vs. 12.43%), Ruminococcus (1.35% vs. 0.53%), and Roseburia (2.31% vs. 1.64%) compared to HD. Post-treatment patients exhibited even greater enrichment of Escherichia-Shigella (4.68%), accompanied by an increase in Bacteroides (8.56% vs. 21.67%) and Anaerostipes (0.85% vs. 1.73%), and pronounced reduction in Faecalibacterium (7.55%), Ruminococcus (0.30%) Roseburia (0.78%), Akkermansia (0.72% vs. 0.09%), and Subdolingranulum (2.50% vs. 1.42%), in contrast to HD ( Figure 2B ).

The ANCOM-BC analysis identified statistically significant microbial taxonomic shifts between pre-treatment, post-treatment, and HD groups, providing insights into the key bacteria associated with the disease and treatment. In pre-treatment patients, when compared to HD, Prevotella 9 showed the highest enrichment (LFC = 2.49) and emerged as a hallmark species of this group. This was accompanied by an increased abundance of several other microbial taxa, notably within the Proteobacteria phylum, including Escherichia-Shigella (p = 0.0319) and Bilophila (p = 0.0019). Additionally, significant enrichment was observed in several bacteria within the Firmicutes phylum, including Streptococcus, Enterococcus, Ruminococcaceae UBA 1819 (p = 0.0020), Parabacteroides, Phascolarctobacterium, and Anaerovoracaceae XII AD3011 (p = 0.0360) ( Figure 3A ).

Post-treatment patients exhibited a marked alteration in their microbial communities. While some similarities to the pre-treatment microbiota were observed, post-treatment individuals demonstrated a more pronounced enrichment of inflammation-associated bacteria. Compared to HD, key bacteria with significant increase included Phascolarctobacterium (p = 0.0039), Escherichia-Shigella (p = 0.0400), Streptococcus (p = 0.0168), Bilophila (p = 0.0009), Parabacteroides (p = 0.0351), and Anaerovoracaceae XII AD3011 (p = 0.0073). Compared to pre-treatment, there was a particular enrichment of Phascolarctobacterium (LFC 4.20), along with Oscillibacter, Collinsella, and Intestinibacter, suggesting potential microbial signatures associated with the therapeutic intervention ( Figure 3B ). Additionally, post-treatment patients exhibited a pronounced depletion of short-chain fatty acid (SCFA)-producing bacteria, including Ruminococcus (p = 0.0206), Christensenellaceae R-7 group (p = 0.0015), Akkermansia, Oscillospiraceae, Lachnospiraceae, Eubacterium, and Clostridium ( Figure 3C ). These results suggest a shift towards a pro-inflammatory microbial profile in the post-treatment state.

The HD group was predominantly characterized by members of the Firmicutes phylum such as Ruminococcus, Christensenellaceae R-7 group, Eubacterium ruminantium, Roseburia, and genera such as Bifidobacterium and Akkermansia, which are associated with intestinal health.

Overall, CC patients are characterized by Prevotella and Escherichia-Shigella, while treatment led to the enrichment of specific bacteria, particularly Phascolarctobacterium, Oscillibacter, Collinsella, and Intestinibacter, while SCFA-associated taxa, including Ruminococcus, Akkermansia, Oscillospiraceae, and Lachnospiraceae, were depleted.

Following these results, a ratio analysis was performed using CLR-transfomed relative abundance of representative taxa ( Figure 3D ) ( Supplementary Material 3 , for ratio calculation), such as Proteobacteria/Firmicutes (HD vs. CC pre- treatment: p = 0.0428; HD vs. CC post- treatment: p = 0.050; CC pre- treatment vs. CC post- treatment: p = 0.9332), Bacteroidota/Firmicutes (HD vs. CC pre-treatment: p = 0.0781; HD vs. CC post-treatment: p = 0.0112; CC pre-treatment vs. CC post-treatment: p = 0.3501), Escherichia/Ruminococcus (HD vs. CC pre-treatment: p = 0.002; HD vs. CC post-treatment: p = 0.0017; CC pre-treatment vs. CC post-treatment: p = 0.7831), Prevotella/Ruminococcus (HD vs. CC pre-treatment: p = 0.0215; HD vs. CC post-treatment: p = 0.1446; CC pre-treatment vs. CC post-treatment: p = 0.7777), Escherichia/Erysipelotrichaceae UCG-003 (HD vs. CC pre-treatment: p = 0.0048; HD vs. CC post-treatment: p = 0.0615; CC pre-treatment vs. CC post-treatment: p = 0.4137).

PICRUSt2 analysis showed that pre-treatment patients exhibit an enrichment of amino acid biosynthesis and metabolism (L-tryptophan biosynthesis, ornithine degradation, L-arginine degradation, L-ornithine degradation, chorismate metabolism, arginine degradation II, and chorismite biosynthesis II), degradation of aromatic compounds and xenobiotics (4-hydroxyphenylacetate degradation, toluene degradation I and II, 3-phenylpropanoate degradation, protocatechuate degradation II, and 4-methylcatechol degradation), inflammation process (enterobactin biosynthesis, enterobacterial common antigen, Kdo2-lipid A and LPS biosynthesis) bacterial stress (ppGpp biosynthesis), antibiotic resistance (polymyxin resistance) and pathways involved in energy metabolism (Glycolysis-TCA-GLYOX-Bypass, TCA-GLYOX-Bypass, sulfoglycolysis, glyoxylate cycle). In contrast, the HD group present an enrichment of adenosylcobalamin biosynthesis and octane oxidation ( Figure 4A ).

Post-treatment patients, in comparison with HD exhibit enrichment of antioxidant-related pathways (ubiquinol-8, -7, -9, and -10 biosynthesis), pathways involved in energy metabolism (TCA cycle IV, Glycolysis-TCA-GLYOX-Bypass, Fatty acid and beta-oxidation), amino acid metabolism (L-tryptophan biosynthesis, and L-histidine degradation II), degradation of aromatic compounds and xenobiotics (toluene degradation I and II, 3-phenylpropanoate degradation, and 4-hydroxyphenylacetate degradation), inflammation-related pathways (enterobactin biosynthesis, enterobacterial common antigen biosynthesis, and Kdo2-lipid A biosynthesis), bacterial stress response (ppGpp biosynthesis), and antibiotic resistance (polymyxin resistance). In contrast, the HD group presented an enrichment of adenosylcobalamin biosynthesis I, ethylmalonyl-CoA pathway, methyl aspartate cycle, formaldehyde assimilation I, taurine degradation, isopropanol biosynthesis ( Figure 4B ).

Comparing pre- and post-treatment patients, pre-treatment showed an enrichment of pathways involved in xenobiotic degradation and detoxification (formaldehyde assimilation I, isopropanol biosynthesis, and sulfolactate degradation), and pathways related to amino acid metabolism (L-threonine metabolism, taurine degradation, and phenylwthylamine degradation), other pathways such as phenylacetate degradation I, along with methylcatechol degradation II, were also observed. In contrast, in the post-treatment group, enrichment of degradation of nucleotide derivatives (allantoin degradation IV and pyrimidine ribonucleosides degradation) and L-histidine degradation II were observed ( Figure 4C ).

An initial frequency analysis of NK cell populations ( Figure 5A ) revealed a significant increase in total NK cell frequency in both patient groups compared to healthy women (HD vs. CC pre-treatment: p = 0.014, HD vs. CC post-treatment: p = 0.015). In the pre-treatment group, there was a trend toward an expansion of the CD56^dim NK cell subset, although this increase did not reach statistical significance (HD vs. CC pre-treatment: p = 0.1312), accompanied by a decrease in the CD56^bright NK cell population (HD vs. CC pre-treatment: p = 0.1299). Conversely, the post-treatment group exhibited the opposite trend, with a reduction in the CD56^dim population (HD vs. CC post-treatment: p = 0.6263) and a concomitant increase in the CD56^bright NK cell subset (HD vs. CC post-treatment: p = 0.7816). Notably, a direct comparison between pre- and post-treatment patients revealed a near-significant difference in the frequency of both CD56^dim (p = 0.0508) and CD56^bright NK cell subsets (p = 0.0809) ( Figure 5B ). NK cells were further subdivided into four subsets based on CD56 and CD16 expression: CD56^dimCD16⁺, CD56^dimCD16^-, CD56^brightCD16⁺, and CD56^brightCD16^- (48). Among these, we observed a significant expansion of the CD56^brightCD16⁺ population in post-treatment patients compared with both HD (p = 0.0239) and CC pre-treatment groups (p = 0.0370). In contrast, the frequencies of the other subsets did not show significant differences across groups ( Supplementary Figure 4 ). Analyses of immune checkpoint expression were conducted only on the classical CD56^dim and CD56^bright NK cell subsets.

The CD56^dim NK cell population revealed a significant increase of PD-1 expression in treatment-naive patients compared to HD (HD vs. CC pre-treatment: p = 0.0457), which was further amplified post-treatment (HD vs. CC post-treatment: p < 0.0001; CC pre-treatment vs. CC post-treatment: p = 0.0261). A similar trend was noted for LAG-3 (HD vs. CC pre-treatment: p = 0.0313; HD vs. CC post-treatment: p < 0.0001; CC pre-treatment vs. CC post-treatment: p = 0.0306) and BTLA (HD vs. CC pre-treatment: p = 0.0051; HD vs. CC post-treatment: p < 0.0007). However, no significant difference was observed in BTLA expression between pre- and post-treatment patient groups (CC pre-treatment vs. CC post-treatment: p = 0.4516). No statistically significant differences were found for the markers TIM-3, TIGIT, and NKG2A between the groups ( Figure 6A ).

In the CD56^bright NK cell population, a similar pattern of immune checkpoint expression was observed, with an increase in the percentage of positive cells in the CC pre-treatment group, followed by an exacerbation post-treatment, compared to the HD group. However, these differences did not reach statistical significance between the two CC patient subgroups. Notably, this trend was particularly evident for LAG-3 (HD vs. CC pre-treatment: p < 0.0001; HD vs. CC post-treatment: p < 0.0001; CC pre-treatment vs. CC post-treatment: p = 0.5857) and TIM-3 (HD vs. CC pre-treatment: p = 0.0018; HD vs. CC post-treatment: p = 0.0008; CC pre-treatment vs. CC post-treatment: p = 0.6651). Though TIGIT followed a similar trend, statistical significance was only reached in the post-treatment group compared to HD (HD vs. CC pre-treatment: p = 0.1038; HD vs. CC post-treatment: p = 0.0148; CC pre-treatment vs. CC post-treatment: p = 0.3632). Interestingly, PD-1 expression increased significantly in both CC groups compared to HD. Post-treatment patients exhibited a decrease in PD-1 expression compared to pre-treatment, although this did not reach statistical significance (HD vs. CC pre-treatment: p < 0.0001; HD vs. CC post-treatment: p = 0.0002; CC pre-treatment vs. CC post-treatment: p = 0.2414). No significant changes were observed in the percentage of cells positive for NKG2A and BTLA between the groups ( Figure 6B ).

The co-expression analysis revealed consistent evidence of a putative exhausted phenotype in NK cells across both patient groups, with a notable exacerbation following treatment, especially in the CD56^dim population. Significant increases were detected in the co-expression of several inhibitory receptors: PD-1⁺BTLA⁺ (HD vs. CC pre-treatment: p = 0.3513; HD vs. CC post-treatment: p = 0.0029; CC pre-treatment vs. CC post-treatment: p = 0.0329), PD-1⁺LAG-3⁺ (HD vs. CC pre-treatment: p = 0.3683; HD vs. CC post-treatment: p = 0.0105; CC pre-treatment vs. CC post-treatment: p = 0.0829), PD-1⁺TIM-3⁺ (HD vs. CC pre-treatment: p = 0.3798; HD vs. CC post-treatment: p = 0.0063; CC pre-treatment vs. CC post-treatment: p = 0.0533), PD-1⁺TIGIT⁺ (HD vs. CC pre-treatment: p = 0.7292; HD vs. CC post-treatment: p = 0.0025; CC pre-treatment vs. CC post-treatment: p = 0.0061), and TIGIT⁺TIM-3⁺ (HD vs. CC pre-treatment: p = 0.9917; HD vs. CC post-treatment: p = 0.0488; CC pre-treatment vs. CC post-treatment: p = 0.0438). In contrast, the co-expression of NKG2A⁺TIGIT⁺ showed a different trend, with a significant decrease in both patient groups compared to healthy controls (HD vs. CC pre-treatment: p = 0.0017; HD vs. CC post-treatment: p = 0.0054; CC pre-treatment vs. CC post-treatment: p = 0.8195). No significant differences were found in the co-expressions of PD-1⁺NKG2A⁺ and NKG2A⁺TIM-3⁺ ( Figure 7A ).

In CD56^bright NK cells, a significant increase in populations exhibiting signs of immune exhaustion was observed, including PD-1⁺BTLA⁺ (HD vs. CC pre-treatment: p = 0.0001; HD vs. CC post-treatment: p = 0.0048), PD-1⁺LAG-3⁺ (HD vs. CC pre-treatment: p < 0.0001; HD vs. CC post-treatment: p = 0.0191), PD-1⁺TIM-3⁺ (HD vs. CC pre-treatment: p = 0.0064; HD vs. CC post-treatment: p = 0.0261), PD-1⁺TIGIT⁺ (HD vs. CC pre-treatment: p = 0.0001; HD vs. CC post-treatment: p = 0.0052), PD-1⁺NKG2A⁺ (HD vs. CC pre-treatment: p = 0.0013; HD vs. CC post-treatment: p = 0.0108), TIGIT⁺TIM-3⁺ (HD vs. CC pre-treatment: p = 0.0170; HD vs. CC post-treatment: p = 0.0003), and NKG2A⁺TIM-3⁺ (HD vs. CC pre-treatment: p = 0.0216; HD vs. CC post-treatment: p = 0.0073), compared to the healthy control group. No significant differences were found between the pre- and post-treatment patient subgroups. Consistent with the individual PD-1 expression patterns observed in this cell population, post-treatment patients tended to decrease the percentage of PD-1⁺ cells co-expressing other immune checkpoints compared to pre-treatment patients. Nevertheless, these changes did not reach statistical significance ( Figure 7B ).

To evaluate the association between microbiota dysbiosis and NK cell exhaustion, a Spearman correlation analysis was performed between the dysbiosis score and NK cell exhaustion scores derived from CD56^dim and CD56^bright subpopulations, as well as their global exhaustion score. A positive correlation was observed between the dysbiosis score and the global NK cell exhaustion score (R = 0.50, p < 0.0001), suggesting that higher levels of dysbiosis are associated with increased NK cell exhaustion when using this approach. When analyzing individual NK cell subsets, the dysbiosis score showed a significant positive correlation with both CD56^dim (R = 0.44, p < 0.0001) and CD56^bright NK cell exhaustion scores (R = 0.42, p < 0.00), indicating that alterations in microbiota composition could be linked to exhaustion across different NK cell populations ( Supplementary Figure 5 ).

The random forest (RF) classification model identified the most important predictors for CC as the expression of PD-1, LAG-3, and their co-expression on CD56^bright NK cells, along with the Escherichia/Ruminococcus ratio ( Figure 8A ). The model demonstrated robust discriminative power with an average ROC-AUC of 0.950 (SD = 0.0545). Repeated 10-fold cross-validation (5 repetitions) confirmed the stability and reliability of the model, with a mean ROC-AUC of 0.935 ( Figure 8B ). This RF model effectively distinguished pre-treatment CC patients from HD, as shown by the confusion matrix ( Supplementary Figure 6 ), achieving an accuracy of 92.86% (p = 0.00455, 95% CI: 0.6613-0.9982), with a kappa coefficient of 0.8571. The model yielded a sensitivity of 87.5%, specificity of 100%, positive predictive value of 100%, and negative predictive value of 85.71%, resulting in a balanced accuracy of 93.75%. Individual ROC curve analyses were performed for the five top-ranked variables identified by the RF model. These features exhibited moderate-to-high discriminative power, with AUCs ranging from 0.778 to 0.891 ( Supplementary Figure 7 ).

Subsequent simplified logistic regression using the top predictors stepwise selected revealed that increased PD-1 expression on CD56^bright NK cells (OR = 1.81, 95% CI: 1.31-2.84, p = 0.002) and a higher Escherichia/Ruminococcus ratio (OR = 14.0, 95% CI: 1.64-176, p = 0.023) were significantly associated with increased CC risk ( Figure 8C ). This logistic regression model yielded an AUC of 0.905 ( Figure 8D ). The final predictive equation was:

An XGBoost classification model was applied to predict survival outcomes in CC patients. Feature importance analysis identified the most relevant predictors as the expression of TIGIT⁺ CD56^bright NK cells, CD56^dim NK cells, and the Proteobacteria/Firmicutes ratio, followed by other immune and microbial features ( Figure 9A ). The model achieved an AUC of 0.875 (SD = 0.236). Repeated cross-validation yielded a mean ROC-AUC of 0.806 (SD = 0.119) ( Figure 9B ). Despite a limited sample size, this model effectively distinguished survival from death in CC patients, as shown by the confusion matrix ( Supplementary Figure 8 ). The model achieved an accuracy of 83.3% (95% CI: 0.3588-0.9958), with a sensitivity of 75.0%, specificity of 100%, balanced accuracy of 87.5%, and a kappa value of 0.6667.

The stepwise algorithm revealed TIGIT⁺TIM-3⁺ CD56^bright and CD56^dim NK cell frequency as the most important variables. Posterior logistic regression showed that elevated TIGIT⁺TIM-3⁺ CD56^bright (OR = 1.23, 95% CI: 1.04-1.62, p = 0.050) and CD56^dim NK cells (OR = 2.89, 95% CI: 1.38-13.22, p = 0.0499) were associated with higher mortality ( Figure 9C ). The two-variable model yielded an AUC of 0.885 ( Figure 9D ). The final predictive equation was:

A second logistic regression using TIGIT⁺TIM-3⁺ CD56^bright NK cells and Proteobacteria/Firmicutes (log) ratio did not reveal any statistically significant associations with mortality ( Supplementary Figure 9A ). The combined model yielded an AUC of 0.823 ( Supplementary Figure 9B ), but these findings should be considered preliminary.

Kaplan–Meier survival analysis was conducted over a 15-month follow-up period to evaluate the prognostic significance of the identified markers. Analyses were restricted to pre-treatment CC patients, as RCT profoundly alters both the gut microbiota and NK cell phenotypes; therefore, untreated samples provide the most appropriate setting to assess predictors of survival. Of the 27 pre-treatment CC patients, 23 were evaluable for survival analysis. 4 patients were lost to follow-up and, therefore, excluded from survival analyses. Patients alive at the end of the 15-month follow-up were censored at that time. Patients were stratified by cut-off values determined via ROC curve analysis for each variable. Patients were stratified into high and low groups using cut-off values derived from ROC curves based on the survival status of pre-treatment CC patients, selecting the point with maximum sensitivity and specificity for each variable.

Patients with elevated levels of TIGIT⁺TIM-3⁺ expression on CD56^bright NK cells were associated with worse prognosis, with a median survival of 12 months compared to higher survival in the low-expression group (log-rank p=0.01), and a log-rank hazard ratio (HR) of 5.789 (95% CI: 1.503-22.29). In this analysis, 7 deaths occurred in the high-expression group and 2 in the low-expression group, while 3 and 11 patients, respectively, were censored at 15 months ( Figure 10A ). Similarly, a high Proteobacteria/Firmicutes ratio exhibited significantly reduced overall survival (log-rank p = 0.05), with a median survival of 12 months and a log-rank HR of 3.921 (95% CI: 1.059-14.52). For this variable, 7 deaths occurred in the high-ratio group and 2 in the low-ratio group, while 5 and 9 patients, respectively, were censored at 15 months ( Figure 10B ). An analysis for the CD56^dim NK cell population was performed; patients with a high percentage showed a trend toward lower survival probability (log-rank HR = 2.436), but the difference was not statistically significant (log-rank p = 0.1831; 95% CI: 0.6553-9.055) (not shown).

Combined survival analyses were conducted to assess synergistic effects; as these analyses were restricted to extreme categories (high/high vs. low/low variables), the effective sample size was reduced. Patients with high levels of both TIGIT⁺TIM-3⁺ expression on CD56^bright NK cells and a high frequency of CD56^dim NK cells had a significantly reduced survival with a mean of 11.5 months (Mantel-Haenszel HR = 17.83, 95% CI: 2.742-115.9, log-rank p = 0.0026) compared to those with low levels of both variables. In this comparison, 5 deaths occurred in the high/high group, whereas no deaths were observed in the low/low group; 7 patients in the low/low group and 1 in the high/high group were censored at 15 months ( Figure 10C ). Similarly, patients presenting both high TIGIT⁺TIM-3⁺ expression on CD56^bright NK cells and a high Proteobacteria/Firmicutes ratio exhibited significantly reduced survival, with a mean of 11 months (HR = 10.68, 95% CI: 1.918-59.44, log-rank p = 0.0054), compared to patients with low values for both markers. For this variable combination, 5 deaths occurred in the high/high group and 1 in the low/low group, while 6 patients in the low/low group and none in the high/high group were censored at 15 months ( Figure 10D ).