Distinct gut virome profiles are associated with response to anti-PD-1 therapy in non-small cell lung cancer.

Introduction
Lung cancer is a leading malignancy worldwide, with approximately 2.48 million new cases in 2024, accounting for about one-eighth of all cancers [1]. Non–small cell lung cancer (NSCLC) constitutes 80–85% of all lung cancer cases and is characterized by poor prognosis, exhibiting one of the lowest five-year survival rates among solid tumors [2, 3]. Since early-stage NSCLC is typically asymptomatic, most patients are diagnosed at an advanced or metastatic stage, which limits opportunities for curative treatment [3]. Although conventional cytotoxic chemotherapy long served as the primary treatment, the introduction of targeted therapies and immune checkpoint inhibitors (ICIs) has substantially improved survival in selected patient subgroups [4–7]. In particular, ICIs targeting the programmed cell death protein 1/programmed cell death ligand 1 (PD-1/PD-L1) axis have reshaped the therapeutic landscape for advanced NSCLC [8, 9]. However, only about 20–30% of patients derive durable clinical benefit, while the majority develop primary or acquired resistance [10, 11], underscoring the need to identify biomarkers and mechanisms that influence treatment response.
In recent years, the gut microbiota has attracted considerable attention as a key modulator of host immunity and a determinant of ICI efficacy. Several studies have demonstrated that specific bacterial taxa, such as Bifidobacterium, can positively regulate antitumor immunity in vivo [12]. Patients enriched with Ruminococcaceae and Faecalibacterium exhibit enhanced systemic and intratumoral immune responses, characterized by more effective antigen presentation and more robust effector T-cell functions in both the peripheral circulation and the tumor microenvironment [13]. In NSCLC patients, Akkermansia muciniphila (Akk) has been consistently associated with favorable responses to PD-1 blockade, underscoring the significant impact of host–microbiome interactions on immunotherapy outcomes [14].
In contrast to the rapidly expanding literature on gut bacteria, research on the gut virome remains relatively limited, particularly with respect to bacteriophages, which constitute its dominant component. Accumulating evidence indicates that bacteriophages can shape bacterial community composition and metabolic activity, modulate mucosal immune homeostasis, and thereby potentially influence host antitumor immunity [15, 16]. However, whether and how the gut virome influences the response to ICIs in patients with NSCLC remains to be systematically elucidated.
To address this gap, we integrated fecal metagenomic sequencing data and clinical information from two cohorts of NSCLC patients treated with PD-1 inhibitors. We reanalyzed the primary cohort of 338 samples [14] and used an independent dataset of 30 samples [17] as an external validation set for the random forest model, aiming to characterize virome alterations associated with treatment response. Furthermore, we investigated the interactions between Akkermansia muciniphila and gut viruses to uncover potential virus–bacteria interplay that may influence immunotherapy outcomes. Our study provided novel insights into the underappreciated role of the gut virome in shaping the efficacy of immune checkpoint blockade.

Materials and methods

Sample information and preprocessing
This study analyzed a total of 368 fecal metagenomic samples collected from patients with advanced NSCLC receiving ICI therapy, sourced from two independent clinical cohorts (Table S1). The primary cohort consisted of 338 patients with histologically confirmed advanced NSCLC enrolled in a clinical study investigating the association between gut microbiota and ICI efficacy [14]. Based on RECIST 1.1 criteria [18], treatment responses were categorized as follows: 75 patients achieved complete or partial response (CR/PR), 102 had stable disease (SD), and 161 had progressive disease (PD). Accordingly, the cohort was stratified into 75 responders and 263 non-responders (comprising 102 SD and 161 PD cases). Patients received standard ICI therapies, including nivolumab or atezolizumab as second-line treatment following platinum-based chemotherapy failure, or pembrolizumab—either as monotherapy or combined with chemotherapy—for those with PD-L1 expression ≥ 1%.
The validation cohort comprised 30 patients with advanced NSCLC who received ICI treatment between 2019 and 2020. This cohort included 11 patients with CR/PR, 7 with SD, and 12 with PD. Treatment involved first-line pembrolizumab monotherapy or second-line nivolumab or atezolizumab.
Raw metagenomic sequencing data from both cohorts were obtained from the Sequence Read Archive (SRA) of the National Center for Biotechnology Information (NCBI). The corresponding BioProject accession numbers are PRJNA751792 and PRJNA1068493 respectively. Quality control of raw sequencing reads was conducted using fastp with parameters “-u 30 -q 20 -l 60 -y -trim_poly_g” [19]. High-quality reads were subsequently aligned to the human reference genome GRCh38 using Bowtie2 [20] to remove host-derived sequences.

Gut virome analysis
We constructed a comprehensive gut virome reference catalog, termed the Chinese Gut Virome Catalog (cnGVC) [21], which was built from over 10,000 publicly available human fecal metagenomic datasets and comprises more than 90,000 non-redundant viral operational taxonomic units (vOTUs). High-quality reads from all samples were aligned against the cnGVC database using Bowtie2 [20], with viral species-level delineation defined at 95% nucleotide identity.
To generate vOTU abundance profiles for each fecal sample, reads mapped to each vOTU was counted and normalized by the total number of mapped reads per sample to obtain relative abundances. Relative abundances of vOTUs belonging to the same viral family were subsequently aggregated to determine family-level viral abundance.
Functional annotation of viral proteins was performed using DIAMOND [22], against the Kyoto Encyclopedia of Genes and Genomes (KEGG) database [23], with the parameters “–query-cover 50 –subject-cover 50 -e 1e-5 –min-score 50 –max-target-seqs 50”. Each protein was assigned to a KEGG Ortholog (KO) identifier based on the top-scoring hit in the database.
Host prediction was performed using both CRISPR spacer matching and homology-based approaches. Specifically, CRISPR-based host assignments were inferred from spacer–vOTU alignments identified by BLASTn (bit score ≥ 45 with full-length coverage). Homology-based predictions were derived from BLASTn alignments between vOTUs and prokaryotic genomes showing ≥ 90% nucleotide identity across ≥ 30% of the vOTU length. These curated host annotations were subsequently used for downstream functional and host-association analyses in this study.

Gut bacteriome analysis
Taxonomic profiling of bacterial communities was performed on all fecal metagenomes using MetaPhlAn4 [24]. The relative abundances of microbial tax were estimated and subsequently aggregated at the genus and species levels to generate comprehensive taxonomic profiles.

Functional comparison of responder- and non-responder-enriched vOtus
For each KOs, its prevalence within the responder (R) or non-responder (NR) group was calculated as the proportion of vOTUs encoding that KO relative to the total number of vOTUs in the group. Fisher’s exact test was performed using the fisher.test function in the R platform to assess whether the prevalence of each KO differed significantly between R and NR groups. KOs with p-values < 0.05 were considered significantly different prevalent.

Statistical analysis
All statistical analyses and visualizations were performed in the R platform unless otherwise specified. Alpha-diversity indices, including the Shannon index and Simpson index, were computed using the diversity function in the vegan package [25], while observed vOTU numbers were calculated using the specnumber function. Differences in alpha-diversity between groups were assessed using two-tailed Wilcoxon rank-sum tests. Beta-diversity was assessed based on Bray–Curtis dissimilarities, calculated with the vegdist function. Principal coordinate analysis (PCoA) was conducted on the Bray–Curtis distances using the cmdscale function in vegan. Permutational multivariate analysis of variance (PERMANOVA) with 999 permutations was performed using the adonis2 function (vegan package) to test for compositional differences between groups.
Differential abundance analysis of microbial features (viruses and bacteria) was performed using microbiome multivariate association with linear models (MaAsLin2), adjusting for potential confounders (sex, age, and antibiotic exposure) [26]. Using the R group as the reference, comparisons were made between R vs. SD (stable disease) and R vs. PD (progressive disease). Only features with a minimum relative abundance > 0.01% and a prevalence > 10% were retained for analysis, and the MaAsLin2 model was specified as “CPLM”. For each microbial feature, the two p-values from the R vs. SD and R vs. PD comparisons were combined using Fisher’s method. Features present in both comparisons with a combined p-value < 0.05 were considered differentially abundant between R and NR group.
Correlations between significantly altered gut viruses and bacteria were evaluated using Spearman’s rank correlation. Correlation pairs with an absolute coefficient |ρ| > 0.6 and p-value < 0.05 were retained for network construction and visualized using the ggraph package.
Random forest models were built using viral markers, bacterial markers, or a combination of both. Model training was followed by 10-fold cross-validation. Model performance was evaluated based on the area under the receiver operating characteristic curve (AUC), calculated with the roc function. Feature importance was ranked using the importance function. The robustness of the optimal model was further validated using the independent validation cohort.

Results

Altered gut virome diversity and composition across clinical response groups
To investigate differences in the gut virome among advanced NSCLC patients with distinct clinical outcomes (R, SD, and PD), we first compared α-diversity at the vOTU level using the Shannon index, Simpson index, and observed vOTU numbers (Fig. 1A). Both the SD and PD groups exhibited significantly lower Shannon diversity than the R group (p < 0.05), with the lowest values observed in PD and the highest in R, indicating stepwise reduction in viral community diversity associated with poorer clinical outcomes. We next conducted PCoA and PERMANOVA based on Bray–Curtis dissimilarities. The first two principal coordinates explained 6.2% and 10% of the total variance, respectively, and PERMANOVA confirmed that clinical response status had a significant effect on overall virome composition (adonis, p < 0.05; Fig. 1B). These findings suggest a marked disruption of the gut virome in individuals with less favorable treatment responses. At the viral family level, stacked bar plots showed that Microviridae and Winoviridae were the most dominant families across all three groups (Fig. 1C).

Identification of viral signatures of different clinical response groups
To identify virome signatures associated with differential clinical responses to PD-1 inhibitor therapy in advanced NSCLC, we performed differential abundance analyses using MaAsLin2 with the R group as the reference, adjusting for sex, age, and antibiotic exposure. At the family level, Crevaviridae was the only viral family that differed significantly in both the R vs. SD and R vs. PD comparisons (p < 0.05; Fig. 2A). Fisher’s combined test further confirmed its significant association when comparing the R group to the combined NR groups (combined p < 0.05; Fig. 2B).
We then conducted MaAsLin2 at the vOTU level using the same covariate adjustments. After combining the p-values from the R vs. SD and R vs. PD comparisons, a total of 194 vOTUs were found to be significantly enriched in the NR group, whereas 594 vOTUs were enriched in the R group (combined p < 0.05; Fig. 2C; Table S2). Taxonomically, NR-enriched vOTUs predominantly belonged to the families Peduoviridae and Inoviridae, while R-enriched vOTUs were largely affiliated with Herelleviridae and Microviridae (Fig. 2D).
Host prediction analyses revealed that 86.7% of R-enriched vOTUs were bacteriophages, with predicted bacterial hosts mainly from Ruminococcus_D (n = 45), Faecalibacterium (n = 34), Roseburia (n = 19), Bacteroides (n = 16), and Blautia_A (n = 14) (Table S3). In contrast, 75.6% of NR-enriched vOTUs were predicted to be bacteriophages targeting bacterial genera such as Clostridium_M (n = 26), Bacteroides (n = 20), and Escherichia (n = 11) (Fig. 2E). Notably, no R-enriched vOTUs were predicted to infect Escherichia, whereas no NR-enriched vOTUs were predicted to target Ruminococcus_D, suggesting fundamentally different virus–host ecological structures between R and NR.

Altered bacterium–virus co-abundance networks between R and NR
To identify bacterial features for network construction, we first performed differential abundance analysis. Using MaAsLin2 with the R group as reference and adjusting for sex, age, and antibiotic use, we identified 71 bacterial species that were significantly differentially abundant between R and NR groups (Table S4).
To investigate differences in virus–bacteria interactions differ between patients with distinct clinical response, we constructed co-occurrence networks for the R and NR groups based on Spearman correlation (|ρ| > 0.6, p < 0.05). Topological analysis showed that the R network consisted of 228 nodes and 204 edges, with an average degree of 1.964 and an average path length of 1.964. The NR network contained 212 nodes and 190 edges, with both the average degree and average path length equal to 1.960 (Fig. 3A–B; Table S5).
Comparison of node connectivity revealed that 20 features displayed higher degrees in the NR network, including bacterial nodes such as Enterocloster aldensis and Enterocloster bolteae, and viral nodes such as v10923 and v05571, suggesting distinct virus–bacterium ecological organization between R and NR (Fig. 3C–D).
We next compared the shared and group-specific vOTUs associated with bacterial families across the two networks (Fig. 3E). Overall, most vOTUs were shared between the R and NR groups. For Lachnospiraceae and Erysipelotrichaceae, vOTUs were largely shared but also included network-specific subsets in both groups, with more Lachnospiraceae-associated vOTUs uniquely observed in the R network than in the NR network. Notably, among vOTUs related to Ruminococcaceae, a distinct subset of associated vOTUs was unique to the R network, whereas all Ruminococcaceae-associated vOTUs observed in the NR network were shared between the two groups.
Regarding enrichment patterns, vOTUs related to Tannerellaceae, Enterobacteriaceae, and Erysipelotrichaceae tend to be enriched in the NR groups, while vOTUs related to most other bacterial families were primarily enriched in the R group. Consistently, the NR network also contained a higher proportion of NR-enriched vOTUs, with Lachnospiraceae-associated vOTUs showing a greater enrichment ratio in NR than in R.

Functional annotation of differential vOTUs
To explore the potential functional mechanisms of the gut virome associated with ICIs outcomes in advanced NSCLC, we performed KEGG functional annotation on the 788 differentially abundant vOTUs identified in the previous analyses. A total of 17 KOs showed significant differences in their prevalence between the R and NR groups (Fisher’s exact test, p < 0.05; Fig. 4A; Table S6).
Among these, only one KO term (K01449) occurred more frequently in the R group and was therefore classified as R-enriched. All remaining KO terms were detected more frequently in the NR group and were classified as NR-enriched (Fig. 4B). Most NR-enriched KOs encoded metabolic functions, including enzymes involved in amino acid metabolism (e.g., asparagine synthase, K01953), central carbohydrate metabolism (e.g., enolase, K01689), and cofactor/vitamin biosynthesis (e.g., a folate biosynthesis enzyme, K01737, and a vitamin B12/porphyrin-related enzyme, K09883). Notably, the NR group also showed enrichment of UDP-glucuronic acid 4-epimerase (K08679), a key enzyme involved in O-antigen nucleotide-sugar biosynthesis. In contrast, the only R-enriched function was a cell-wall hydrolase (K01449), suggesting enhanced phage-encoded lytic activity that may preferentially target clinically detrimental bacteria and thereby contribute to favorable treatment responses.

Predictive performance of gut viral and bacterial profiles for immunotherapy response
We next built random forest classification models with 10-fold cross-validation to assess the predictive performance of viral, bacterial, and combined virus–bacteria profiles in distinguishing clinical outcomes (R vs. NR) among advanced NSCLC patients treated with PD-1 blockade. In the discovery cohort, the model based solely on viral features demonstrated the highest discriminative power (AUC = 76.8%, 95% CI: 74.3%–79.2%), followed by the combined model (AUC = 76.2%, 95% CI: 73.7%–78.7%), whereas the bacterium-only model showed the lowest performance (AUC = 66.4%, 95% CI: 63.6%–69.1%) (Fig. 5A–B). Feature importance analysis revealed that the top 20 predictors in both the virus-only and combined models consisted predominantly of viral taxa enriched in the R group. In the combined model, most of the important features were viral, with only three bacterial features ranking among the top predictors (Fig. 5C).
To evaluate model generalizability, we validated all three models in an independent external cohort. Consistent with the discovery cohort, the virus-only model again yielded the best predictive performance (AUC = 74.2%, 95% CI: 53.6%–94.7%), followed by the combined model (AUC = 70.8%, 95% CI: 52.6%–89.1%) and the bacterium-only model (AUC = 67.0%, 95% CI: 46.7%–87.3%) (Fig. 5D).

Association between Akk colonization status and gut virome structure
To further investigate the potential influence of Akk on gut virome structure, we stratified samples based on its presence (Akk-positive) or absence of Akk (Akk-negative) and compared their virome profiles. At the vOTU level, α-diversity indices—including the Shannon, Simpson, and observed vOTUs numbers)—were significantly higher in Akk-positive individuals, indicating greater viral richness and diversity compared to Akk-negative subjects (Fig. 6A; Table S7). We next assessed overall community structure using PCoA based on Bray–Curtis distance. The first two principal coordinates (PCo1 and PCo2) explained 10% and 6.2% of the variance, respectively. PERMANOVA analysis confirmed a significant separation of virome composition between the Akk-positive and Akk-negative groups (R2 = 0.008, p = 0.002; Fig. 6B).
At the viral family level, both groups were dominated by Microviridae and Winoviridae. However, Akk-positive individuals exhibited higher relative abundances of Herelleviridae and Suoliviridae, whereas Peduoviridae and Inoviridae were enriched in Akk-negative samples (Fig. 6C).
To further characterize virome features associated with Akk status, we restricted the analysis the 788 response-associated vOTUs previously identified as differentially abundant between clinical outcomes groups. Afte reanalyzing these vOTUs using MaAsLin2 while adjusting for age, sex, and antibiotic use, Akk-positive and Akk-negative samples remained clearly separated based on this subset of vOTUs (Fig. 6D), suggesting that Akk colonization is independently associated with gut virome composition.

Discussion
Accumulating evidence indicates that the gut microbiome as a key modulator of antitumor immunity and clinical response to ICIs in advanced NSCLC [27–29]. However, although bacterial components have been extensively investigated, the gut virome remains comparatively understudied, despite being dominated by bacteriophages that can regulate bacterial community structure, metabolic outputs, and mucosal immune homeostasis [30]. Here, we conducted a large-scale virome-focused analysis by integrating fecal metagenomic data from 338 NSCLC patients treated with PD-1 blockade, together with an independent external validation cohort. Using a multi-step framework that included differential abundance testing, phage–host prediction, cross-kingdom co-occurrence network analysis, functional annotation, and random forest modeling, we systematically characterized virome alterations associated with treatment response. Notably, virome features showed stronger discriminatory power for predicting clinical response than bacterial features, and they exhibited distinct compositional, ecological, and functional patterns linked to clinical outcomes. These findings establish the gut virome as a pivotal and previously underappreciated determinant of immunotherapy efficacy in NSCLC.

Virome signatures and putative mechanisms linking to clinical response
We observed a stepwise decline in the viral Shannon index with worsening therapeutic outcomes, mirroring the reduced bacterial diversity reported in NR [31]. This concordant pattern indicates a coordinated dysbiosis of the intestinal ecosystem that may accompany resistance to PD-1 blockade. Moreover, NR-enriched vOTUs were predominantly assigned to Peduoviridae and Inoviridae. Notably, members of Inoviridae have been reported to carry virulence-associated genes and can enhance bacterial pathogenicity by converting otherwise non-virulent hosts into strains with increased virulence [32].
Host prediction analysis further revealed that these NR-enriched vOTUs frequently targeted bacterial genera including Clostridium_M, Bacteroides, and Escherichia. The Clostridium_M has been reported to be enriched in patients with non-gastrointestinal cancers, consistent with clinical observations [33]. In addition, Bacteroides has been consistently associated with poor prognosis during ICI therapy across multiple cancer types, including NSCLC [34], melanoma [35], and gastric cancer [36]. The enrichment of these taxa may not merely reflect a dysbiosis state but could indicate an active, virus-mediated modulation of their functional activity.
Supporting this notion, functional annotation of the differentially enriched vOTUs provides crucial molecular insights. KOs significantly overrepresented in the NR group includeK01448 (amiABC; N−acetylmuramoyl−L−alanine amidase) and K01689 (ENO1_2_3; enolase). K01448 is involved in peptidoglycan hydrolysis—a key step in viral progeny release—and also regulates bacterial cell division and wall metabolism [37]. K01689 is a central glycolytic enzyme whose human homologues (ENO1 and ENO3) are upregulated in colorectal and pancreatic cancers, with high expression correlating positively with poor patient prognosis and advanced disease [38, 39]. The enrichment of these viral-encoded genes points to a potential two-pronged mechanism fostering a resistance phenotype in NR patients: 1) direct disruption of bacterial community integrity via enzymes like K01448 that interfere with cell wall dynamics, and 2) viral reprogramming of core host bacterial metabolism (e.g., glycolysis via K01689), which may subsequently alter the gut metabolite landscape to favor an immunosuppressive tumor microenvironment.
In stark contrast, phages enriched in the R group predominantly targeted bacterial genera with established beneficial roles in host immunity, including Ruminococcus_D, Faecalibacterium, and Roseburia. Consistent with our observations, Ruminococcus was significantly enriched in patients with stage III/IV NSCLC who responded to ICI treatment [40], and the Ruminococcaceae has also been associated with response to anti-PD-1 treatment in melanoma [41], and higher baseline Faecalibacterium abundance predicts better outcomes in multiple cancers [40, 41]. Roseburia, a key butyrate producer, has been shown to enhance anti-PD-1 efficacy in preclinical models by promoting butyrate-mediated activation of CD8+ T cells [42, 43]. This is particularly relevant as NSCLC responders exhibit higher PD-1 expression on peripheral CD8+ T cells, indicative of a more activated state [44].
Our co-occurrence network analysis further substantiated these findings, revealing distinct ecological architectures. Beyond differences in abundance, vOTUs linked to Ruminococcaceae (encompassing Faecalibacterium and Ruminococcus) and Lachnospiraceae (encompassing Roseburia) formed interaction modules that were unique to the R network and absent or less prominent in the NR network. This finding suggests that, in responders, these beneficial bacteria are embedded within a specialized phage–bacteria ecological framework, rather than merely being present at higher abundance.
Integrating this evidence, we propose that R-enriched phages infecting genera like Roseburia may modulate the abundance and/or metabolic output of their hosts, thereby influencing the production of immunoregulatory metabolites such as butyrate. This phage-mediated tuning could enhance CD8+ T-cell activation and potentially affect PD-1 expression dynamics, collectively contributing to a more effective anti-tumor immune response. Consequently, butyrate-related readouts and specific Roseburia-phage interaction patterns emerge as promising candidate biomarkers for predicting clinical outcomes under anti–PD-1 outcomes.

Clinical translation: predictive modeling and ecological insights
Accurate patient stratification prior to immune checkpoint inhibitor therapy remains an unmet clinical need. Our machine learning analysis directly addresses this challenge by demonstrating the superior predictive value of the gut virome. A random forest classifier based solely on viral features achieved the highest prediction performance (AUC = 76.8%), significantly outperforming a model using only bacterial data (AUC = 66.4%). Crucially, this advantage was robustly confirmed in an independent external validation cohort, where the virus-only model maintained the highest AUC. Notably, integrating bacterial features did not enhance model performance and even led to a slight decrease compared to the virus-only model. This suggests that bacterial signatures may contain more unrelated variants, thereby weakening the stronger and more specific signal provided by viral markers. These findings highlight the potential of virome based classifiers as clinically valuable tools for pre-treatment stratification and highlight the need to incorporate virome analysis into precision immuno-oncology research.
Previous studies have established the relative abundance of Akk as a reliable biomarker for favorable prognosis in NSCLC patients receiving PD-1 blockade [14], However, its relationship with the gut virome in this context was undefined. Our analysis revealed that Akk-negative individuals harbored a gut virome with significantly lower alpha diversity compared to Akk-positive subjects. Furthermore, the virome composition in Akk-negative patients was dominated by viral families such as Peduoviridae and Inoviridae, a pattern that mirrors the signatures enriched in non-responders. These findings suggest that Akk may not act in isolation but instead that its colonization may depend on, or contribute to, a healthier and more diverse gut virome ecosystem.

Limitations and future perspectives
Our study has several limitations should be acknowledged. First, although the primary cohort is sizable for metagenomic analysis, the validation cohort is relatively small; future multi-center prospective studies with larger sample sizes are needed to confirm the generalizability of our findings. Second, both viral host prediction and functional annotation depend on existing reference databases, which remain incomplete and may not fully capture novel or uncharacterized phage taxa. Third, the observational nature of our study precludes causal inference. Future mechanistic investigations using gnotobiotic animal models or targeted phage manipulation will be essential to establish whether the identified viral signatures actively modulate antitumor immunity. Finally, profiling fecal samples may not fully capture the mucosa-associated microbial communities or localized immune interactions within the tumor microenvironment.
Despite these limitations, our findings collectively highlight the gut virome as a critical and previously underappreciated determinant of immunotherapy response in NSCLC. By integrating assessments of viral diversity, taxonomic, ecological networks, and functional potential, we show that phage communities can distinguish between immune-favorable and immune-dysregulated gut ecosystems with a resolution superior to bacterial markers alone. These insights advocate for the incorporating virome profiling into the framework of microbiome-guided precision immunotherapy and suggest that targeted modulation of phage–bacteria interactions could represent a novel strategy to enhance antitumor efficacy. Future research aimed at dissecting the mechanistic basis of phage–host–immune crosstalk will be crucial for translating virome signatures into clinically actionable interventions.

Conclusion
In conclusion, this large-scale, integrative analysis establishes the gut virome as a critical determinant of clinical response to anti-PD-1 therapy in patients with advanced NSCLC. Clinical outcomes were associated with coordinated remodeling of the gut viral community, including changes in diversity, taxonomic composition, ecological interaction networks, and functional potential. Notably, the virus-only classifier consistently outperformed the bacteria-based model in both the discovery and validation cohorts, highlighting the distinct and robust predictive signal carried by the viral component of the microbiome. We further delineated a poor-response virome signature characterized by Peduoviridae and Inoviridae phages that preferentially target potential pathobionts and are enriched in auxiliary metabolic genes, in contrast to a favorable-response signature marked by Herelleviridae phages targeting beneficial butyrate-producing genera and forming more stable ecological networks. Together, these findings support a model in which the gut virome contributes to shaping host responses to ICIs, potentially through phage-mediated modulation of bacterial communities and their immunoregulatory outputs. Our results underscore the value of incorporating virome profiling into future microbiome-informed strategies for precision immunotherapy and suggest that targeting specific phage–bacteria interactions may represent a promising therapeutic avenue.

Electronic supplementary material
Below is the link to the electronic supplementary material.