Abstract

   The intratumoral mycobiome plays a crucial role in the tumor
   microenvironment, but its impact on renal cell carcinoma (RCC) remains
   unclear. We collected and quantitatively profiled the intratumoral
   mycobiome data from 1044 patients with RCC across four international
   cohorts, of which 466 patients received immunotherapy. Patients were
   stratified into mycobiota ecology-depauperate and mycobiota
   ecology-flourishing (MEF) groups based on fungal abundance. The MEF
   group had worse prognosis, higher fungal diversity, down-regulated
   lipid catabolism, and exhausted CD8^+ T cells. We developed the
   intratumoral mycobiota signature and intratumoral mycobiota-related
   genes expression signature, which robustly predicted prognosis and
   immunotherapy outcomes in RCC and other cancers. Aspergillus tanneri
   was identified as a potential key fungal species influencing RCC
   prognosis. Our findings suggest that the intratumoral mycobiome
   suppresses lipid catabolism and induces T cell exhaustion in RCC.
     __________________________________________________________________

   Intratumoral mycobiome modulates tumor microenvironment, affecting
   immunotherapy outcomes in renal cancer.

INTRODUCTION

   Renal cell carcinoma (RCC) constitutes approximately 90% of all
   malignant kidney tumors, with clear cell renal cell carcinoma (ccRCC)
   being the most prevalent subtype, characterized by high metastatic
   potential and poor prognosis ([70]1–[71]4). The intratumoral
   microbiome, a crucial component of the tumor microenvironment, has been
   demonstrated to exist widely in various tumor tissues and is closely
   associated with cancer initiation, progression, metastasis, and
   treatment outcomes ([72]5–[73]11). Intratumoral microbiota can
   influence the host through multiple aspects, such as modulating
   immunogenicity and altering metabolic activities, which consequently
   affect the biological behavior of the tumor and patient prognosis
   ([74]12–[75]16). While previous studies have mainly focused on the
   bacterial microbiota, less attention has been given to the fungal
   microbiota (mycobiota), which is also an important element of the human
   microbiome. Existing evidence has confirmed that mycobiota may play a
   crucial role in the occurrence and progression of tumors ([76]17,
   [77]18). Therefore, further investigation into the role of intratumoral
   mycobiota in RCC is of great significance.

   Immune checkpoint inhibitors (ICIs) enhance the capacity of the immune
   system to attack cancer cells by blocking or inhibiting immune
   checkpoint molecules ([78]19, [79]20). However, only 10 to 20% of
   patients with RCC benefit from immunotherapy, and the complex
   mechanisms behind this are still being explored ([80]21, [81]22).
   Recent studies have suggested that the intratumoral microbiome may
   modulate the tumor immune microenvironment, subsequently influencing
   patients’ responses to immunotherapy ([82]14, [83]16, [84]23). In
   colorectal and gastrointestinal cancers, the intestinal mycobiota in
   patients has been demonstrated to be closely associated with prognosis,
   tumor progression, and response to immunotherapy ([85]24–[86]26).
   Studies of the intratumoral mycobiota have reported that it plays an
   important role in various types of tumors by modulating immune
   responses and regulating metabolic processes ([87]6, [88]27–[89]29).
   However, the impact of intratumoral mycobiota on patients with RCC and
   its role in immunotherapy remain to be elucidated.

   In this study, we integrated multiple clinical cohorts and multi-omics
   datasets to investigate the impact of the intratumoral mycobiome on the
   tumor microenvironment, prognosis, and immunotherapy outcomes in
   patients with RCC.

RESULTS

Prevalence of mycobiota in RCC and its association with prognosis

   To comprehensively investigate the ecology of intratumoral mycobiota in
   RCC, our study incorporated a total of 1084 patients with RCC from five
   distinct cohorts, encompassing 21 countries. The intratumoral mycobiome
   data of 522 patients with RCC were derived from the study by
   Narunsky-Haziza et al. ([90]18) [The Cancer Genome Atlas Kidney Renal
   Clear Cell Carcinoma (TCGA-KIRC)]. We collected data from 56 patients
   with primary RCC at West China Hospital of Sichuan University (WCH-RCC)
   ([91]30). Two large-scale, international, multicenter RCC immunotherapy
   clinical trials—IMmotion150 ([92]31) (n = 162) and IMmotion151 ([93]32)
   (n = 304)—were included. By processing raw RNA sequencing (RNA-seq)
   data through a microbial assignment pipeline, we obtained the
   intratumoral mycobiome profile of patients in the WCH-RCC, IMmotion150,
   and IMmotion151 cohorts (see Materials and Methods). Moreover, we
   retrospectively collected tumor tissues from a cohort of 40 patients
   with RCC at Shanghai Changhai Hospital (CHH-RCC). This cohort was used
   to corroborate the impact of intratumoral mycobiota on the host immune
   microenvironment. To further expand our research, we included a total
   of 2338 patients who received ICI therapy from 19 cohorts encompassing
   different tumor types, including non–small cell lung cancer (NSCLC),
   metastatic urothelial carcinoma (mUC), melanoma, and so on ([94]Fig. 1,
   A and B, and table S1) ([95]33–[96]50).

Fig. 1. Overall study design.

   [97]Fig. 1.
   [98]Open in a new tab

   (A) Basic information of all cohorts included in this study. The world
   map showed the countries or regions of origin for the RCC cohorts
   included in the relevant studies. (B) Study flowchart. NSCLC, non–small
   cell lung cancer; mUC, metastatic urothelial cancer; CRC, colorectal
   cancer; RCC, renal cell carcinoma; PDAC, pancreatic ductal
   adenocarcinoma.

   Among the 35 cancer types included in TCGA, RCC exhibited a higher
   proportion of intratumoral fungal reads relative to total reads
   compared to other cancers, with marked heterogeneity across patients
   ([99]18). In our study, we included 522 patients with RCC from TCGA,
   all of which exhibited the presence of at least one fungal species.
   Ascomycota and Basidiomycota constituted the dominant fungal phyla in
   RCC (fig. S1). We characterized the intratumoral mycobiota landscape of
   TCGA-KIRC patients at the class level ([100]Fig. 2A). Notably,
   Dothideomycetes and Tremellomycetes were predominant in more than
   three-fifths (327 of 522) of the samples, whereas Sordariomycetes
   dominated approximately one-sixth (88 of 522). In nearly one-sixth (79
   of 522) of the samples, Eurotiomycetes and Malasseziomycetes were more
   common. These findings suggested substantial heterogeneity in the
   intratumoral mycobiota composition among different patients with RCC.
   To validate the presence of mycobiota in RCC tissues, we performed
   fluorescence in situ hybridization (FISH) staining on tissue sections
   from 40 patients with RCC in the CHH-RCC cohort using fungal-specific
   28S rRNA probes. All these samples exhibited positive 28S probe
   fluorescence of varying intensities, further confirming the widespread
   presence of mycobiota in the tumor tissues of patients with RCC
   ([101]Fig. 2B).

Fig. 2. Intratumoral mycobiota landscape in RCC and group construction.

   [102]Fig. 2.
   [103]Open in a new tab

   (A) Composition of intratumoral fungi (class level) and clinical
   information of TCGA-KIRC patients. From inner to outer, the mapped
   information includes OS, gender, age, AJCC stage, T stage, N stage, M
   stage, and the relative abundance of intratumoral fungi at the class
   level in each sample. (B) Fluorescence in situ hybridization (FISH)
   showing the presence of fungi in tumor tissues of patients with RCC.
   Scale bar, 100 μm (10×) and 20 μm (40×). (C and D) Consensus clustering
   based on the intratumoral mycobiota abundance in TCGA-KIRC. (C)
   Consensus score matrix for all samples when k = 2. (D) Consensus matrix
   CDF curves for different k values. (E) Kaplan-Meier curves of OS for
   the two groups (named group A and group B) obtained based on consensus
   clustering [P = 2.1 × 10^−4, HR (95% CI) =1.78 (1.31 to 2.42)]. (F)
   Differences in various clinical characteristics between group A and
   group B, with significance tested using the chi-square test.

   We further explored the heterogeneity and clinical implications of the
   intratumoral mycobiome in patients with RCC. The optimal number of
   clusters (k = 2) was determined based on the consensus score matrix and
   cumulative distribution function (CDF) curve ([104]Fig. 2, C and D).
   Using consensus clustering based on the intratumoral mycobiota
   profiles, TCGA-KIRC patients were stratified into two clusters (group A
   and group B). Patients in group B displayed significantly poorer
   prognosis compared to group A [P = 2.1 × 10^−4, hazard ratio (HR)
   =1.78, 95% confidence interval (CI) =1.31 to 2.42, [105]Fig. 2E]. In
   addition, group B comprised a higher proportion of patients with
   advanced American Joint Committee on Cancer (AJCC) stages (P = 0.003)
   and T stages (P = 0.009), as well as a higher proportion of patients
   with metastasis (M stage: P = 0.033). However, there was no difference
   in age between two groups (P = 0.432; [106]Fig. 2F).

Higher fungal diversity and complexity in the poorer prognosis group

   The multiple taxonomic levels of fungi exhibited heterogeneity between
   group A and group B ([107]Fig. 3, A and B, and fig. S2, A and B). The
   alpha diversity analysis revealed significantly higher intratumoral
   fungal richness and diversity in group B, compared to group A
   (P < 0.001, [108]Fig. 3, C to E, and fig. S2C). Beta diversity analysis
   revealed distinct fungal communities between group A and group B
   [permutational multivariate analysis of variance (PERMANOVA) P < 0.001,
   [109]Fig. 3F]. Linear discriminant analysis effect size (LEfSe) was
   used to further identify the significantly differed fungi between group
   A and group B [linear discriminant analysis (LDA) score > 2, P.adj <
   0.05]. At the phylum level, Mucoromycota and Microsporidia were found
   to be significantly enriched in group B. Regarding the class level, 16
   taxa, including Eurotiomycetes, Sordariomycetes, and Malasseziomycetes,
   exhibited higher abundance in group B, while Xylonomycetes were more
   enriched in group A. In terms of the order level, a total of 32 taxa,
   including Hypocreales and Eurotiales, revealed significantly higher
   abundance in group B. In contrast, only Xylonales, Eremomycetales, and
   Dothideales were enriched in group A ([110]Fig. 3G). As for the family
   and genus levels, 53 and 88 taxa, respectively, were significantly
   enriched in group B, while only 2 and 2 taxa were enriched in group A
   (table S2). At each taxonomic level, group B exhibited a greater number
   of enriched taxa. Furthermore, fungal co-occurrence network analysis
   revealed that the network complexity and the degree of fungal
   interactions were markedly higher in group B, with higher number of
   nodes, edges, and average degree ([111]Fig. 3, H and I, and table S3).
   The dominant fungal genera showed differences between the co-occurrence
   networks of group A and group B. In group A, the predominant taxa were
   Fusarium (27.2%) and Malassezia (16.3%), while in group B, the top taxa
   were Aspergillus (20.7%) and Fusarium (6.8%). Notably, Aspergillus and
   Fusarium occupied relatively critical positions in the fungal
   communities of both groups. Together, these results suggested that the
   intratumoral fungal community in group B manifested a more intricate
   and interconnected structure compared to that of group A.

Fig. 3. Heterogeneity of fungal communities between groups.

   [112]Fig. 3.
   [113]Open in a new tab

   (A and B) Relative abundance of fungi at different taxonomic levels in
   group A and group B, showing only the top eight taxa. (A) Family level;
   (B) genus level. (C to E) Comparison of fungal alpha diversity between
   group A and group B using three indices: Shannon, Simpson, and Chao 1.
   (F) Comparison of fungal community differences between group A and
   group B based on beta diversity using principal coordinate analysis
   (PCoA) with Bray-Curtis distance, statistically tested using PERMANOVA.
   (G) Identification of fungi with the abundance differences between
   group A and group B using LEfSe analysis. The bar plot showed the
   results at the phylum, class, and order levels (LDA score > 2, FDR
   <0.05). (H and I) Fungal community co-occurrence networks and network
   topology parameters for group A and group B. The nodes implied
   individual fungi, and edges represented significant Spearman
   correlations between fungi (|ρ| > 0.80, FDR < 0.05). The nodes and
   edges of the networks were colored according to fungi (genus level),
   with the top eight most abundant fungi assigned different colors and
   the remaining fungi colored gray. The blue edges of the network
   represented the positive correlations between fungi. The thickness of
   the edges was the weight of the correlation. (H) Interaction network
   for group A (node: 21, edges: 46, modularity: 0.23, average degree:
   4.38). (I) Interaction network for group B (node: 120, edges: 1396,
   modularity: 0.28, average degree: 23.27). (J) Group definition based on
   the features of fungal communities. LEfSe, linear discriminant analysis
   effect size; LDA, linear discriminant analysis; PERMANOVA,
   permutational multivariate analysis of variance. ***P < 0.001.

   Collectively, the fungal communities demonstrated marked distinctions
   between group A and group B. Multiple community feature analyses
   revealed that the intratumoral fungal community in group B had greater
   diversity and complexity than that of group A. To more precisely
   describe the distinct community attributes of the two groups, we
   designated group A as the mycobiota ecology-depauperate (MED) group,
   while group B was termed the mycobiota ecology-flourishing (MEF) group
   ([114]Fig. 3J).

Down-regulation of lipid catabolism pathway in the MEF group

   To elucidate the mechanisms by which the intratumoral mycobiota
   influences the host microenvironment, clinical outcomes, and tumor
   characteristics, we investigated the differences in host biology
   between the MEF and MED groups. Differential expression analysis
   identified 110 up-regulated and 276 down-regulated genes in the MEF
   group compared to the MED group [|logFC| > 1, false discovery rate
   (FDR) < 0.05, fig. S3A]. These differentially expressed genes (DEGs)
   were primarily enriched in metabolic and immunological pathways,
   including lipid metabolism, amino acid metabolism, and
   interleukin-related signaling ([115]Fig. 4, A and B). In addition, we
   assessed the somatic mutation landscape in the MEF and MED groups,
   identifying 24 genes with significantly different mutation frequencies
   (P < 0.05, fig. S3B). The three genes exhibiting the most significant
   differences were GRIN2B, SPHKAP, and LAMA2, all of which have been
   previously implicated in tumor progression ([116]51–[117]53).

Fig. 4. Unique metabolic characteristics and immune microenvironment between
different groups.

   [118]Fig. 4.
   [119]Open in a new tab

   (A and B) Metabolic and immune-related pathways obtained by ORA based
   on DEGs between groups (different ribbons represented different sources
   of gene sets). Pink, GO database; light blue, KEGG database; lime
   green, Reactome database; orange, WikiPathways; purple, Hallmark
   database. (C) Lipid catabolism pathways significantly down-regulated in
   the MEF group compared to the MED group in GSEA. (D) Spearman
   correlation between lipid catabolism pathways and differentially
   abundant fungi (genus level) between MED and MEF groups. Lipid
   catabolism pathways were obtained from GSEA results (FDR < 0.05), and
   differentially abundant fungi were obtained from LEfSe results (LDA
   score > 2, FDR < 0.05). The size of the bubbles represented the level
   of correlation. The color of the bubbles indicated the FDR level: blue,
   negative correlation; red, positive correlation. (E) Proportions of
   major cells in the immune microenvironment of MED and MEF groups
   obtained using the xCell. (F) Proportions of CD8^+ T cells in MED and
   MEF groups based on CIBERSORT, EPIC, TIMER, and quanTIseq. (G)
   Proportions of CD8^+ T cell subtypes within MED and MEF groups obtained
   using gene sets provided by Zheng et al. ([120]59). Tn, naïve T cells;
   Tem, effector memory T cells; Temra, terminally differentiated effector
   memory or effector; KIR, killer cell immunoglobulin-like receptors;
   Trm, tissue-resident memory T cells; Tex, exhausted T cells; ISG,
   interferon-stimulated genes. (H) Scores of terminal Tex in MED and MEF
   groups. (I) Expression levels of exhausted CD8^+ T cell–related markers
   (CTLA-4, LAG-3, and PDCD1) in MED and MEF groups. (J) Expression levels
   of immunosuppressive genes IDO1, SLAMF7, and VEGFA in MED and MEF
   groups. ORA, overrepresentation analysis; DEGs, differentially
   expressed genes; GSEA, gene set enrichment analysis; NES: normalized
   enrichment score; TME: tumor microenvironment. *P < 0.05, **P < 0.01.

   Gene set enrichment analysis (GSEA) showed significant down-regulation
   of multiple lipid metabolism pathways in the MEF group compared to the
   MED group, including lipid catabolism, neutral lipid catabolism, and
   glycerolipid catabolism ([121]Fig. 4C and fig. S3C). Furthermore, the
   peroxisome proliferator-activated receptor (PPAR) pathway, closely
   linked to lipid metabolic activity ([122]54), was also down-regulated
   in the MEF group (fig. S3D). These findings denoted that the metabolic
   activity of patients in the MEF group underwent substantial changes,
   potentially further leading to the accumulation of lipid accumulation
   within the tumor microenvironment. The association between lipid
   accumulation in RCC and poor prognosis has been well documented in
   previous studies ([123]55, [124]56). To gain further insights, we
   evaluated the specific changes in lipid catabolism pathways among the
   different groups. The expression patterns of key metabolic genes were
   found to be in concordance with the observed pathway changes. Of
   particular note, CPT1A, members of the acyl-CoA dehydrogenase family,
   and PPARs exhibited lower expression levels in the MEF group (fig.
   S3E). To explore the potential relationship between intratumoral fungi
   and lipid metabolism pathways, we performed a correlation analysis
   between the differentially abundant fungal genera (LDA score > 2, FDR <
   0.05) and lipid metabolism pathways. We observed that some fungal
   genera exhibited significant correlations with multiple lipid
   metabolism pathways (FDR < 0.05). Among these, Coccidioides and
   Trichoderma exhibited robust negative correlations with several lipid
   metabolism pathways. In contrast, Anthracocystis and Coniophora
   displayed positive correlations with multiple lipid metabolism pathways
   ([125]Fig. 4D).

Elevated exhausted CD8^+ T cell infiltration in the MEF group

   Lipids play an important role in modulating immune responses to tumor
   and their accumulation within the microenvironment inducing dysfunction
   of CD8^+ tumor-infiltrating lymphocytes ([126]57, [127]58). The immune
   infiltration analysis using xCell showed a higher proportion of CD8^+ T
   cells in the MEF group compared to the MED group ([128]Fig. 4E). To
   verify the composition of CD8^+ T cells across different groups, we
   integrated four immune infiltration algorithms: Cell-type
   Identification By Estimating Relative Subsets Of RNA Transcripts
   (CIBERSORT), Estimating the Proportion of Immune and Cancer cells
   (EPIC), Tumor IMmune Estimation Resource (TIMER), and quantify the
   fractions of ten immune cell types from bulk RNA-sequencing data
   (quanTIseq). All algorithms confirmed a relatively higher proportion of
   CD8^+ T cells in the MEF group, with quanTIseq showing statistical
   significance (P < 0.05, [129]Fig. 4F). To further investigate the
   immune microenvironment characteristics, we calculated the composition
   of different CD8^+ T cell subtypes in the two groups ([130]59). The
   results exhibited significantly higher proportions of terminal and
   TCF7^+ exhausted CD8^+ T (Tex) cells in the MEF group ([131]Fig. 4G).
   We also observed higher proportions of GZMK^+ effector memory T (Tem)
   cells and ZNF683^+CXCR6^+ tissue-resident memory T (Trm) cells in the
   MEF group (P < 0.01). They participate in and dominate the two pathways
   from naive T cells to terminal Tex cells, playing a critical role in
   the development trajectory of T cell exhaustion ([132]59). We further
   evaluated the exhaustion scores and Tex subtypes of different groups
   ([133]60, [134]61). The results suggested a significantly higher
   exhaustion score in the MEF group, with high infiltration of all Tex
   subtypes ([135]Fig. 4H and fig. S3F). Although the overall regulatory T
   (Treg) cell levels did not present substantial differences between
   groups (fig. S3G), both TNFRSF9^+ Treg and ISG^+ Treg subtypes had
   higher proportions in the MEF group (fig. S3H). Among the Tex-related
   markers, the expression levels of cytotoxic T-lymphocyte antigen 4
   (CTLA-4) and Lymphocyte-activation gene 3 (LAG-3) were significantly
   higher in the MEF group (P < 0.05, [136]Fig. 4I). Meanwhile, the PDCD1
   was also relatively up-regulated, although not reaching statistical
   significance (P = 0.070, [137]Fig. 4I). Furthermore, we found
   significantly higher expression of indoleamine 2,3-dioxygenase 1
   (IDO1), signaling lymphocytic activation molecule family member 7
   (SLAMF7), and vascular endothelial growth factor A (VEGFA) in the MEF
   group (P < 0.05, [138]Fig. 4J). These genes have been reported to play
   a crucial role in inducing T cell exhaustion ([139]62–[140]64).

Intratumoral mycobiota signature predicted prognosis and immunotherapy
response in patients with RCC

   To evaluate the association between the fungal abundance
   characteristics and patient prognosis, we developed an integrated
   machine learning (ML) framework. This framework was used to develop an
   intratumoral mycobiota signature (IMS) model associated with overall
   survival (OS). In TCGA-KIRC patients, a total of 224 fungal species
   were detected. We eliminated low-abundance species (abundance of 0 in
   more than 80% of samples) to reduce excessive noise in the model.
   Subsequently, fungi significantly associated with OS were screened
   using a univariate Cox regression model, and the abundances of the 40
   screened fungi were input into the ML framework ([141]Fig. 5A). We
   constructed 100 different ML models and calculated the mean C-index of
   each model across all datasets (fig. S4A). The model with the highest
   mean C-index was selected for further investigation. The results showed
   that the combination of random survival forest (RSF) and bidirectional
   stepwise Cox regression [Stepcox (both)] has the optimal performance
   (mean C-index = 0.613, [142]Fig. 5B). In this ML algorithm combination,
   the RSF model evaluated and selected important variables based on the
   principle of minimum tree depth, obtaining 23 important species.
   Further application of Stepcox regression constructed the IMS model,
   which included seven fungal species ([143]Fig. 5C).

Fig. 5. Construction of IMS to predict prognosis and immunotherapy response
in patients with RCC through integrated machine learning.

   [144]Fig. 5.
   [145]Open in a new tab

   (A) Screening process for selecting fungi as input for machine learning
   model training. (B) Machine learning based on 40 fungi, establishing
   100 predictive models and calculating the C-index of each model in all
   datasets (only the first 10 models were shown here; the complete result
   can be found in fig. S4A). (C) Information on the seven fungi in IMS.
   The left side showed the model coefficients of the fungi; the right
   side showed the HR and 95% CI of each fungus in TCGA-KIRC. (D) KM curve
   of OS obtained by applying IMS scores in TCGA-KIRC [P = 1.6 × 10^−10;
   HR (95% CI) = 2.73 (1.98 to 3.76)]. (E) KM curve of OS obtained by
   applying IMS scores in WCH-RCC [P = 0.016; HR (95% CI) =4.25 (1.18 to
   15.3)]. (F) ROC curve for OS prediction based on IMS scores in
   TCGA-KIRC (AUC = 0.682, 95% CI: 0.633 to 0.731). (G) ROC curve for OS
   prediction based on IMS scores in WCH-RCC (AUC = 0.671, 95% CI: 0.500
   to 0.842). (H) PFS curve obtained by applying IMS scores in the
   IMmotion150 cohort of patients receiving immunotherapy [P = 0.0095, HR
   (95% CI) =1.96 (1.17 to 3.29)]. (I) PFS curve obtained by applying IMS
   scores in the IMmotion151 cohort of patients receiving immunotherapy
   [P = 0.018, HR (95% CI) =1.42 (1.06 to 1.90)]. (J) Clinical response
   rates to immunotherapy [complete response (CR)/partial response (PR)
   and stable disease (SD)/progressive disease (PD)] in the high and low
   IMS score groups of the IMmotion150 cohort, with significance tested
   using the chi-square test. (K) Clinical response rates to immunotherapy
   (CR/PR/SD/PD) in the high and low IMS score groups of the IMmotion151
   cohort, with significance tested using the chi-square test.

   We divided patients with RCC into high-score (IMS-High) and low-score
   (IMS-Low) groups based on their IMS scores. Patients in the IMS-High
   group had significantly worse OS compared to the IMS-Low group in both
   the TCGA-KIRC cohort (training dataset) and the WCH-RCC cohort
   (validation dataset) [TCGA-KIRC: P = 1.6 × 10^−10; HR (95% CI) = 2.73
   (1.98 to 3.76); WCH-RCC: P = 0.016; HR (95% CI) =4.25 (1.18 to 15.3);
   [146]Fig. 5, D and E]. The receiver operating characteristic (ROC)
   curves showed that the area under the curve (AUC) of IMS scores was
   0.682 (95% CI: 0.633 to 0.731; [147]Fig. 5F) in TCGA-KIRC and 0.671
   (95% CI: 0.500 to 0.842; [148]Fig. 5G) in WCH-RCC. Moreover, the
   multivariate Cox regression model indicated that IMS remained a
   significant prognostic factor after adjusting for clinical
   characteristics including age, AJCC stage, and TNM stage (P < 0.05,
   fig. S4, B and C). The above evidence substantiated the reliability and
   robustness of IMS in predicting the prognosis of patients with RCC.

   Previous studies have suggested that the intratumoral mycobiome is
   closely associated with the tumor immune microenvironment and may
   influence the efficacy of ICIs in treated patients ([149]28, [150]65).
   Accordingly, we further assessed the relationship between IMS and the
   prognosis of patients receiving immunotherapy in two international
   multicenter RCC cohorts (IMmotion150 and IMmotion151). The IMS-High
   group manifested significantly worse progression-free survival (PFS)
   compared to the IMS-Low group [IMmotion150: P = 0.0095, HR (95% CI)
   =1.96 (1.17 to 3.29); IMmotion151: P = 0.018, HR (95% CI) =1.42 (1.06
   to 1.90), [151]Fig. 5, H and I]. The IMS-High group had a significantly
   lower response rate to immunotherapy (IMmotion150: P = 0.044;
   IMmotion151: P = 0.003, [152]Fig. 5, J and K).

Aspergillus tanneri as the key fungus influencing RCC prognosis

   We carried out an in-depth analysis of the fungi in IMS. To identify
   the key species, we integrated the results from LEfSe, survival
   analysis, and ML through the framework of “rationality-congruence
   assessment” (RCA) (see Materials and Methods, fig. S5A). We performed a
   multidimensional interpretation of the fungi to select the key species
   with consistent performance through the framework. Among the seven
   fungi involved in the IMS, only Aspergillus tanneri achieved an RCA
   score of 6 (fig. S5B). The feature importance analysis also revealed
   that the importance of A. tanneri in the IMS model was higher than the
   other species ([153]Fig. 5C). Consequently, we considered A. tanneri as
   the key fungus and executed further analyses.

   In the TCGA-KIRC cohort, A. tanneri was detected in nearly 80%
   (415/522) of tumor tissues. It showed a significant positive
   correlation with the alpha diversity of the fungal community (Shannon,
   Simpson, ACE, and Chao1: P < 0.001, fig. S6A). Survival analysis showed
   that the high-abundance group of A. tanneri had a significantly worse
   OS [P = 2.2 × 10^−8, HR (95% CI) = 2.31 (1.71 to 3.13), [154]Fig. 6A].
   In the high–A. tanneri abundance group, the tumor progression was more
   advanced as demonstrated by higher AJCC stage, T stage, and M stage
   (fig. S6B). Moreover, we noticed that the patients in the high–A.
   tanneri abundance group, the IMS-High group, and the MEF group had a
   high degree of overlap ([155]Fig. 6B), all showing poorer prognosis. We
   further confirmed the presence of A. tanneri in tumor tissues within
   the CHH-RCC cohort through quantitative polymerase chain reaction
   (qPCR) analysis (fig. S6C). The above results indicated that A. tanneri
   had a notable impact on the prognosis of patients with RCC.

Fig. 6. A. tanneri as a key fungus influencing prognosis in patients with
RCC.

   [156]Fig. 6.
   [157]Open in a new tab

   (A) KM curve of OS obtained by grouping TCGA-KIRC patients into high-
   and low-abundance groups based on A. tanneri abundance
   [P = 2.2 × 10^−8, HR (95% CI) = 2.31 (1.71 to 3.13)]. (B) Sankey
   diagram showed the distribution of samples in TCGA-KIRC from consensus
   clustering groups based on mycobiota abundance to high and low IMS
   score groups obtained from machine learning, and then to high and low
   A. tanneri abundance groups. MED and MEF groups obtained from consensus
   clustering; IMS: High and low score groups obtained based on IMS from
   machine learning; A. tanneri abundance: Groups obtained based on high
   and low A. tanneri abundance. (C) Spearman correlation between the
   seven fungi in IMS and lipid catabolism pathways. Lipid catabolism
   pathways were selected from GSEA results (FDR < 0.05). The size of the
   bubbles represented the level of correlation. The color of the bubbles
   indicated the FDR level: blue, negative correlation; red, positive
   correlation. (D) Proportions of CD8^+ T cells between high and low A.
   tanneri abundance groups obtained from five algorithms (TCGA-KIRC). (E)
   Expression levels of exhausted CD8^+ T cell–related markers (CTLA-4,
   LAG-3, and PDCD1) in high and low A. tanneri abundance groups. (F and
   G) Representative images of A. tanneri FISH staining and multiplex
   immunofluorescence in the CHH-RCC cohort. FISH: A. tanneri (red); DAPI
   (blue). Scale bar, 80 μm (10×) and 20 μm (40×). Multiplex
   immunofluorescence: positive cells detected, CD8 (red), PD-1 (green),
   LAG-3 (purple), CTLA-4 (orange), and DAPI (blue). Scale bar, 80 μm
   (10×) and 20 μm (40×). (H to J) Spearman correlation scatter plot
   between A. tanneri and PD-1^+ CD8^+ T cells (H), LAG-3^+ CD8^+ T cells
   (I), and CTLA-4^+ CD8^+ T cells (J) in the CHH-RCC cohort. *P < 0.05;
   **P < 0.01.

   The correlation analysis showed that A. tanneri was significantly
   negatively correlated with lipid catabolism pathways (FDR < 0.01,
   [158]Fig. 6C). In addition, A. tanneri had a negative correlation with
   the PPAR pathway (P = 0.002, ρ = −0.136, fig. S6D). Notably, among all
   the fungi found so far, the A. tanneri NIH1004 strain has been reported
   to have the strongest specialized metabolic capacity ([159]66).
   Together, there may be a complex and close cross-talk between A.
   tanneri and host metabolic activities.

   To further explore the potential impact mechanism of A. tanneri on the
   host, we compared the composition of the immune microenvironment
   between the high- and low-abundance group of A. tanneri. Similar to the
   MEF group, the proportion of CD8^+ T cells was higher in the high–A.
   tanneri abundance group with worse prognosis (P < 0.05, fig. S6E). The
   multiple immune infiltration algorithms provided consistent evidence
   ([160]Fig. 6D). Meanwhile, we found that Treg cells showed
   significantly higher infiltration in the high–A. tanneri abundance
   group (fig. S6E). In terms of the expression of Tex-related markers,
   CTLA-4, LAG-3, and PDCD1 were up-regulated in the high–A. tanneri
   abundance group (CTLA-4: P = 0.038, LAG-3: P = 0.004, PDCD1: P = 0.050;
   [161]Fig. 6E).

   In patients with CHH-RCC, FISH revealed the widespread presence of A.
   tanneri in tumor tissues ([162]Fig. 6, F and G). To confirm the
   association between A. tanneri and T cell exhaustion, we performed a
   correlation analysis on the results of FISH and multiplex
   immunofluorescence. We found that A. tanneri abundance exhibited a
   significant positive correlation with PD-1 and LAG-3 (PD-1: ρ = 0.423,
   P = 0.007; LAG-3: ρ = 0.351, P = 0.027; [163]Fig. 6, H and I). However,
   no notable correlation was observed between A. tanneri and CTLA-4
   (P = 0.945, [164]Fig. 6J). These results suggested that A. tanneri may
   influence patient prognosis by modulating the immune microenvironment.

Intratumoral mycobiota-related gene expression signature predicted prognosis
and immunotherapy response

   In this study, we found a substantial association between the
   intratumoral mycobiota and patient prognosis in patients with RCC.
   Considering that the interaction between fungi and the host is
   generally reflected in gene expression patterns, we investigated the
   gene expression signatures related to intratumoral mycobiota. An ML
   prognostic model was constructed based on these signatures.

   Three distinct gene sets were defined to screen for gene expression
   features associated with intratumoral mycobiota. The first one included
   genes significantly correlated with at least one fungal species that
   had differential abundance between the MEF and MED groups (FDR < 0.05).
   The second gene set consisted of DEGs between the MEF and MED groups
   (|logFC| > 1, FDR < 0.05). The third set encompassed genes
   significantly associated with prognosis (univariate Cox regression
   model: P < 0.05). By taking the intersection of the three sets, we
   obtained a gene collection associated with the characteristics of the
   intratumoral mycobiome (n = 69, [165]Fig. 7A). On the basis of these 69
   genes and through the integration of ML algorithms, we constructed the
   intratumoral mycobiota-related genes expression signature (IMRGES).
   Each model was fitted in the training set TCGA-KIRC and the validation
   set WCH-RCC, and the mean C-index of the models was calculated
   ([166]Fig. 7B and fig. S7A). Among all ML models, the combination of
   bidirectional stepwise Cox regression and RSF exhibited the best
   performance with the highest mean C-index (mean C-index = 0.925).
   Eventually, this combination was used to construct the IMRGES model,
   which included 23 genes ([167]Fig. 7C).

Fig. 7. Construction of IMRGES through integrated machine learning based on
intratumoral fungi-related gene sets.

   [168]Fig. 7.
   [169]Open in a new tab

   (A) Venn diagram showing the intersection of three conditional gene
   sets (n = 69). Fungi-related genes: genes correlated with at least one
   fungal species that had differential abundance between the MEF and MED
   groups (FDR < 0.05); prognostic genes: genes that can independently
   predict patient prognosis (Cox proportional hazards test P < 0.05);
   DEGs: differentially expressed genes between MED and MEF groups
   (|logFC| > 1, FDR < 0.05). (B) Machine learning based on 69
   fungi-related genes, establishing 92 predictive models and calculating
   the C-index of each model in all datasets (only the first 15 models
   were shown here; the complete result can be found in fig. S7A). (C)
   Importance index of the 23 genes in the selected IMRGES. (D) KM curve
   of OS for IMRGES scores in the TCGA-KIRC cohort [P = 2.4 × 10^−120; HR
   (95% CI) =52 (31 to 86)]. (E) KM curve of OS for IMRGES scores in the
   WCH-RCC cohort [P = 2.8 × 10^−7; HR (95% CI) =16 (4 to 66)]. (F) KM
   curve of PFS for IMRGES scores in the IMmotion150 cohort of patients
   receiving immunotherapy [P = 0.039, HR (95% CI) =1.55 (1.02 to 2.34)].
   (G) KM curve of PFS for IMRGES scores in the IMmotion151 cohort of
   patients receiving immunotherapy [P = 0.012, HR (95% CI) =1.52 (1.10 to
   2.11)]. (H) Clinical response rates to immunotherapy (CR/PR/SD/PD) in
   the high and low IMRGES score groups of the IMmotion150 and IMmotion151
   cohorts, with significance tested using the chi-square test. (I) Forest
   plot showing the results of a meta-analysis of the impact of IMRGES on
   OS in patients with different cancers at the pan-cancer level. The
   squares represented the study-specific HRs. The horizontal lines
   implied the 95% CI. The diamond indicated the summary HR and its
   corresponding 95% CI. The dashed vertical line represented the summary
   HR.

   According to the IMRGES, patients were divided into a high–IMRGES score
   group (IMRGES-High) and a low–IMRGES score group (IMRGES-Low). In the
   training set TCGA-KIRC and the validation set WCH-RCC, the IMRGES-High
   group demonstrated significantly worse OS [TGCA-KIRC:
   P = 2.4 × 10^−120; HR (95% CI) =52 (31 to 87); WCH-RCC:
   P = 2.8 × 10^−7; HR (95% CI) =16 (4 to 66); [170]Fig. 7, D and E]. The
   ROC curve showed that the AUC of the IMRGES score was 0.952 (95% CI =
   0.935 to 0.970; fig. S7B) in TCGA-KIRC and 0.869 (95% CI = 0.744 to
   0.994; fig. S7C) in WCH-RCC. In addition, multivariate Cox regression
   analysis revealed that IMRGES remained significant (P < 0.05, fig. S7,
   D and E) after adjusting for clinical characteristics such as age, AJCC
   stage, and TNM stage, indicating that IMRGES has an independent
   predictive effect on the OS of patients with RCC.

   In the RCC cohort receiving immunotherapy, IMRGES also demonstrated the
   capability to differentiate patients’ PFS [IMmotion150: P = 0.039, HR
   (95% CI) =1.55 (1.02 to 2.34), [171]Fig. 7F; IMmotion151: P = 0.012, HR
   (95% CI) =1.52 (1.10 to 2.11), [172]Fig. 7G]. Furthermore, it exhibited
   promising predictive performance for the response of patients with RCC
   to immunotherapy (IMmotion150: P = 0.026; IMmotion151: P = 0.201;
   [173]Fig. 7H).

   To further validate the predictive potential of IMRGES, we collected
   transcriptome data from 2338 patients across six cancer types for
   extensive validation. All patients underwent ICI therapy. We conducted
   a meta-analysis to examine the relationship between IMRGES and OS in
   these cohorts, using a random-effects model to account for intercohort
   heterogeneity (P < 0.01, I^2 = 59%, [174]Fig. 7I). Our findings
   indicated that, at the pan-cancer level, higher IMRGES scores were
   frequently associated with poor patient OS (HR = 1.506, 95% CI =1.195
   to 1.897, P < 0.01, [175]Fig. 7I).

DISCUSSION

   In this study, we systematically investigated the intratumoral
   mycobiota in patients with RCC, revealing the heterogeneity of fungal
   communities among different patients and their correlation with
   prognosis and immune response. The widespread presence of fungi was
   detected in RCC tumor tissues. On the basis of the intratumoral
   mycobiome abundance, patients could be divided into two groups, MED and
   MEF groups. In the MEF group, the prognosis of patients was
   significantly worse compared to the MED group. Prior researches have
   demonstrated that the diversity of intratumoral bacterial microbiota
   was frequently associated with patient prognosis ([176]7, [177]11,
   [178]67). These findings underscore the intimate relationship between
   the intratumoral microbiota and patients. Nevertheless, the connection
   between the intratumoral mycobiota and the prognosis and immunotherapy
   outcomes of patients with RCC remained unexplored. Our in-depth
   analysis in this study revealed that the MEF group harbored a
   significantly higher fungal community diversity and complexity,
   implying that the heterogeneity of the intratumoral mycobiota was also
   closely related with patient prognosis. The MEF group displayed
   significant down-regulation of lipid catabolism pathways and an immune
   microenvironment characterized by high infiltration of Tex cells.
   Moreover, we constructed the IMS based on the fungal community, which
   could robustly discriminate the prognosis and response to immunotherapy
   of patients with RCC. We also pinpointed A. tanneri as a potential key
   species influencing the prognosis of patients with RCC. Ultimately, by
   developing the IMRGES, we extended the predictive capability of fungi
   to the transcriptome level and further applied it to other cancer
   types, thereby expanding the application value of our research.

   Lipid metabolism reprogramming is a hallmark of cancer that frequently
   shapes a unique tumor microenvironment and cancer cell phenotype,
   playing a pivotal role in tumor development and progression ([179]56,
   [180]68). In the MEF group, several pathways related to lipid
   catabolism were significantly down-regulated, and the expression of
   certain key genes involved in lipid metabolism was markedly reduced.
   This change potentially leads to the continuous accumulation of lipids
   in the tumor microenvironment. Lipid accumulation is a prominent
   feature in RCC, especially in ccRCC, signifying malignant tumor
   progression and poor prognosis ([181]69, [182]70). Moreover, we found
   correlations between certain intratumoral fungal genus and lipid
   metabolism pathways. These results suggested that the intratumoral
   mycobiota may influence tumor progression by modulating host metabolic
   activities. Prior studies have shown that intratumoral microbes and
   their derived metabolites can affect the host’s metabolic patterns.
   They played a role in regulating the tumor microenvironment and
   affecting long-term prognosis and immunotherapy efficacy ([183]5,
   [184]15, [185]27, [186]71–[187]73). The intratumoral mycobiota
   stimulates cancer cells to secrete interleukin-33 (IL-33), which, in
   turn, recruits and activates T[H]2 and innate lymphoid cells 2 within
   the pancreatic ductal adenocarcinoma (PDAC) tumor microenvironment and
   further stimulates the tumor growth ([188]27). Liu et al. ([189]29)
   have found that Aspergillus sydowii could accelerate the lung tumor
   progression by proliferation and stimulation myeloid-derived suppressor
   cells through IL-1β. This process leads to the suppression of cytotoxic
   T lymphocyte cell activity and the accumulation of PD-1^+ CD8^+ T
   cells. The intratumoral microbiota could also mediate the function of
   immune effector cells in the tumor microenvironment, influencing
   antitumor immunity and tumor biological behavior by altering tumor
   immunogenicity ([190]9, [191]71, [192]74–[193]77). In our study, the
   MEF group had a highly exhausted state of CD8^+ T cells despite their
   high infiltration. Simultaneously, this group showed high expression of
   T cell inhibitory receptors in the tumor tissue, and certain Treg
   subtypes were significantly up-regulated. These Treg cells played a
   crucial role in tumor immune escape and could induce the generation of
   Tex cells ([194]59, [195]78–[196]80). The immune microenvironment in
   the MEF was in a highly suppressed state according to the above
   analyses. Evidence from various studies have pointed that lipid
   metabolism imbalance in tumor tissues can induce CD8^+ T cell
   dysfunction ([197]57, [198]58, [199]81). In RCC, patients with worse
   prognosis always have tumor immune microenvironments enriched with
   CD8^+ T cells exhibiting reduced cytotoxic activity ([200]41, [201]82).
   Therefore, in the MEF group, the significant down-regulation of lipid
   catabolism pathways may lead to lipid accumulation in the tumor
   microenvironment, thereby inhibiting the antitumor activity of CD8^+ T
   cells and causing them to exhibit an exhausted state.

   Compared to the MED group, we observed significantly higher expression
   levels of IDO1, SLAMF7, and VEGFA in the MEF group. These genes have
   been reported to play a role in promoting or sustaining T cell
   exhaustion within the tumor microenvironment ([202]62–[203]64). IDO1
   has been shown to drive immunosuppressive effects through the
   metabolism of tryptophan into kynurenine. It would ultimately suppress
   effector T cell responses and foster an environment conducive to T cell
   dysfunction. Similarly, SLAMF7 activation can induce multiple
   inhibitory receptors on T cells and reprogram them toward a more
   exhausted phenotype. VEGFA, a potent mediator of angiogenesis, also
   exerts immunomodulatory functions by impairing T cell infiltration,
   enhancing the presence of immunosuppressive cells, and contributing to
   the up-regulation of T cell exhaustion markers. Therefore, beyond the
   notable up-regulation of immune exhaustion markers such as LAG3, the
   higher expression of genes capable of indirectly suppressing T cell
   function may also play a crucial role in molding the T cell exhaustion
   state. Notably, IDO1 has been identified as a crucial immunoregulatory
   enzyme connecting host metabolic pathways and microbial signals in
   previous studies ([204]83, [205]84). Zhang et al. ([206]85) found that
   the gut commensal bacterium Dubosiella newyorkensis produces l-lysine
   (Lys), which activates IDO1 and thereby induces the immunotolerogenic
   capacity of dendritic cells. This metabolic circuitry ultimately gives
   rise to an immunosuppressive microenvironment. While much attention has
   focused on bacterial contributions, fungal constituents may similarly
   alter or amplify IDO1-related pathways. In conclusion, evaluating IDO1
   expression alongside fungal community changes may provide insights into
   the broader role of host-mycobiota metabolic cross-talk in shaping
   antitumor responses.

   Recent studies have highlighted the intimate association between
   intratumoral microbiome, particularly bacterial microbiota, and the
   prognosis or outcome of immunotherapy ([207]5, [208]14, [209]86).
   However, the relationship between intratumoral mycobiota and
   immunotherapy remains unclear. In the context of the gut microbiome,
   gut fungi have exhibited superior predictive capabilities compared to
   gut bacteria in terms of patient response to ICIs ([210]25, [211]87).
   In our study, we observed that intratumoral mycobiota in RCC
   potentially played an essential role in reshaping the immune
   microenvironment. By using integrated ML algorithms, we identified the
   most representative and influential fungi combinations from the fungal
   microecological community to construct the IMS. This signature could be
   used to predict patient survival and differentiate patients’ responses
   to immunotherapy. Furthermore, applying the framework of RCA, we
   identified A. tanneri as a potentially key fungal species influencing
   the RCC fungal community. Medema et al. ([212]66) have indicated that
   A. tanneri had the most robust metabolic capabilities among currently
   found fungal species. Our findings demonstrated that A. tanneri was
   intimately associated with the host immune microenvironment.
   Consequently, there existed a potential nexus between intratumoral
   fungi and immunotherapy efficacy. Further studies are warranted to
   elucidate the intricate interaction mechanisms between fungi and the
   tumor immune microenvironment ([213]18, [214]27, [215]88).

   Our research had certain limitations. Although we have demonstrated the
   critical role of A. tanneri in RCC through various analyses, the
   difficulty in obtaining and culturing this species posed technical
   challenges that currently prevented experimental verification of its
   effects in RCC. We have not investigated the potential impact of the
   metabolic products derived from the intratumoral mycobiota and A.
   tanneri on the tumor microenvironment. Moreover, our study only deeply
   explored A. tanneri, while other important components of the
   intratumoral mycobiota required further attention and investigation.

   In conclusion, by integrating data from multiple large-scale clinical
   cohorts and multi-omics datasets, this study systematically unveiled
   the heterogeneity of the RCC intratumoral mycobiota and its intricate
   associations with host metabolism and immune status. These findings
   offered different perspectives for future research on the role of
   intratumoral mycobiota in tumor development and immunotherapy. Further
   comprehensive exploration into the fungal communities within tumors may
   provide unique insights into the mechanisms underlying tumor occurrence
   and progression. In addition, such research could guide the development
   of precise immunotherapeutic strategies and advance the application of
   precision medicine approaches in the field of oncology.

MATERIALS AND METHODS

Study population and sample collection

   The intratumoral mycobiome data for TCGA-KIRC were derived from a
   comprehensive pan-cancer study investigating the fungal ecologies
   ([216]18). The RNA-seq data, mutation data, and clinical information
   for TCGA-KIRC were obtained using the TCGAbiolink R package ([217]89).
   We collected data from 56 patients with WCH-RCC ([218]30). In addition,
   we included patients who received ICI therapy from two large-scale,
   international, multicenter RCC clinical trials, IMmotion150 ([219]31)
   (n = 162) and IMmotion151 ([220]32) (n = 304). In both cohorts, all
   patient samples were collected before treatment initiation, and the
   patients had not received any prior systemic therapy for RCC before
   undergoing immunotherapy. In IMmotion150 cohort, patients received
   either atezolizumab monotherapy or atezolizumab plus bevacizumab. In
   IMmotion151 cohort, patients received atezolizumab plus bevacizumab. A
   summary table of baseline information for patients in different cohorts
   is shown in table S4.

   Moreover, we included transcriptome data from 2338 patients receiving
   immunotherapy across six types of cancers, including melanoma, NSCLC,
   mUC, colorectal cancer (CRC), RCC, and PDAC.

   For multiplex immunofluorescence and FISH, we retrospectively collected
   40 patients who were clinically diagnosed with ccRCC and underwent
   renal cancer surgery at Changhai Hospital between 2019 and 2021 in our
   study (CHH-RCC). The baseline characteristics of patients were
   described in table S4. The inclusion criteria were as follows: (i)
   preoperative imaging and clinical evaluations indicated resectable
   renal tumors; (ii) postoperative pathological diagnosis confirmed
   ccRCC; and (iii) comprehensive initial medical records were available,
   including detailed medical history, imaging data, and preoperative
   laboratory test results. The exclusion criteria were as follows: (i)
   receipt of targeted therapy, immunotherapy, interventional therapy, or
   other antitumor treatments before surgery; (ii) antibiotic use within 1
   month before surgery; (iii) postoperative pathological findings
   indicating non-ccRCC subtypes; and (iv) presence of congenital
   autoimmune diseases or other malignancies documented in the medical
   records. This study was approved by the Ethics Committee of Changhai
   Hospital (ethical approval no. CHEC2021-091). Following surgical
   resection, tumor tissues were fixed in 10% neutral buffered formalin
   and embedded in paraffin blocks for preservation. We sectioned these
   samples into 4- to 5-μm-thick slices for the study.

   For qPCR, tumor tissues from 16 patients who underwent radical
   nephrectomy for RCC were collected from Changhai Hospital. Patient
   recruitment and sample collection were approved by the Ethics Committee
   of Changhai Hospital (no. CHEC2021-091). Before the start of this
   study, written informed consent was obtained from each participant. All
   cases were pathologically confirmed as ccRCC post-surgery. To ensure
   the integrity and sterility of the samples, we performed strict aseptic
   procedures. Tumor tissues were collected from patients under sterile
   conditions in the operating room. Tissue samples were excised within 30
   min of removal from the body using sterile surgical instruments. The
   tumor tissues were carefully dissected, rinsed with sterile saline, and
   immediately stored in sterile cryovials. Samples were transferred to
   liquid nitrogen tanks within 15 min of collection. To control for
   environmental contamination, four 1.5-ml sterile saline tubes were used
   as negative controls. During the sampling process, the cryovials
   containing sterile saline were left open and exposed to the operating
   room environment for 5 min. Negative controls were transported and
   stored under the same conditions as the tissue samples. These controls
   were subsequently used to identify and exclude potential contamination
   introduced during tissue collection and downstream processing.

Microbial assignments

   A microbial assignments workflow was established based on the research
   by Narunsky-Haziza et al. ([221]18). (i) Host depletion: We constructed
   a host depletion pipeline using a combination of fastp ([222]90),
   Minimap2 ([223]91), and SAMtools ([224]92). Initially, quality control
   was conducted on the raw RNA-seq data using fastp, and low-quality
   reads (based on default fastp parameters) were eliminated.
   Quality-controlled reads were subsequently aligned to the human genome
   GRCh38.p14 and the Phi X 174 viral genome using Minimap2. Successfully
   aligned reads were removed using SAMtools. (ii) Microbial assignments:
   The host-depleted sequences were aligned to a standard microbial
   reference genome constructed based on the RefSeq release 220 (rep220)
   database using Bowtie2 ([225]93). Last, the alignment results were
   inputted into Woltka ([226]94) to obtain a microbial abundance matrix
   based on operational genomic units. All reference genomes were
   downloaded in September 2023.

Transcriptomic quantification

   The raw reads were subjected to quality control using fastp (under the
   same conditions as previously mentioned) and then the expression of
   transcripts was quantified using Salmon v1.10.2 ([227]95) based on the
   human reference genome (GRCh38.p14). The gene-level expression matrix
   was obtained using the tximport package. With the exception of the
   TCGA-KIRC, WCH-RCC, IMmotion150, and IMmotion151 cohorts, the source of
   transcriptomic data for the remaining datasets is shown in table S1.

Consensus clustering

   Consensus clustering was performed based on the intratumoral mycobiome
   abundance in RCC. The optimal number of clusters was identified using
   the consensus score matrix and CDF curve. This process was executed
   using the ConsensuClusterPlus ([228]96) R package.

Microbial signature analysis

   The alpha diversity indices, including ACE, Chao1, Shannon, and
   Simpson, were calculated using the vegan package ([229]97). The
   Wilcoxon rank-sum test was used to compare the statistical differences
   in alpha diversity between two groups. Principal coordinate analysis
   (PCoA) based on Bray-Curtis distances was used to assess the beta
   diversity of the fungal communities. PERMANOVA was performed to compare
   the differences in fungal community structure between groups, and
   visualization was conducted using the microeco ([230]98) package. LEfSe
   was applied to identify differentially abundant fungal taxa from the
   phylum to genus level (LDA score > 2, FDR < 0.05). To compare the
   ecological complexity of the fungal communities in different groups, we
   calculated the Spearman correlation between species to construct a
   co-occurrence network. Species with a prevalence below 10% in the
   sample were excluded from the analysis. The P values for the
   correlations were corrected using the FDR method. The threshold for the
   existence of interactions between species was set as a Spearman
   correlation >0.8 and an FDR <0.05. Subsequently, we visualized the
   network structure and computed their topological characteristics,
   including the number of nodes and edges, average degree, and other
   indicators.

Differential gene expression analysis and pathway enrichment analysis

   To explore the potential differences in biological characteristics
   between the MED and MEF groups, the edgeR ([231]99) package was used to
   identify DEGs between the groups. The thresholds for |fold change| and
   FDR were set to >2 and <0.05, respectively. Subsequently, we performed
   overrepresentation analysis (ORA) on the DEGs, comprehensively
   considering the Kyoto Encyclopedia of Genes and Genomes (KEGG), Gene
   Ontology (GO), Reactome, WikiPathways, and Hallmark databases for
   pathway annotation. ORA was conducted using the BioEnricher package.
   GSEA ([232]100) was performed using the clusterProfiler ([233]101)
   package to identify up-regulated and down-regulated pathways in the MED
   group. The gene sets were obtained from the MSigDB database. The
   enrichment P values were corrected for multiple testing using the FDR
   method. Visualization was performed using the GseaVis ([234]102)
   package. Single-sample GSEA (ssGSEA) was conducted using the “GSVA”
   ([235]103) package.

Immune infiltration

   We performed immune infiltration analysis using the IOBR ([236]104)
   package. To ensure the robustness of the assessment, we used five
   immune infiltration algorithms, including CIBERSORT ([237]105), EPIC
   ([238]106), TIMER ([239]107), xCell ([240]108), and quanTIseq
   ([241]109). In addition, we further evaluated cell subtypes and
   calculated exhaustion scores based on ssGSEA. The composition of CD8^+
   T cell subtypes and Treg subtypes was computed using gene sets provided
   by Zheng et al. ([242]59) and Li et al. ([243]61). The exhaustion
   scores were calculated based on the exhaustion-related gene sets
   provided by Sade-Feldman et al. ([244]60).

Integrated ML framework

   To further investigate the impact of intratumoral fungi on patient
   prognosis and immunotherapy, we integrated 10 ML algorithms to predict
   OS, including CoxBoost, generalized boosted regression modeling, Lasso,
   Ridge, elastic net (Enet) ([245]110), RSF ([246]111), partial least
   squares regression for Cox ([247]112), survival support vector machine
   ([248]113), supervised principal components ([249]114), and stepwise
   Cox. A total of 100 algorithm combinations were obtained to identify
   the most representative and stable species or genes in the complex
   intratumoral fungal microecological environment. Some combinations were
   skipped due to their inability to effectively reduce the dimensionality
   of the feature space. We constructed the IMS and IMRGES from the
   mycobiome and transcriptome levels, respectively. For IMS, the sparsity
   and high discreteness of species abundance in the fungal community
   could introduce noise and affect the accuracy of IMS. Therefore, we
   first removed low-abundance species (fungi with zero abundance in more
   than 80% of the samples). Subsequently, fungi significantly associated
   with OS were identified through univariate Cox analysis and served as
   candidate inputs for the model. For IMRGES, to maximize the transfer of
   fungal characteristics to the genome and enhance its broad
   applicability, we defined three major gene sets to serve as “bridges.”
   Set 1: Fungal-related genes. We calculated the Spearman correlation
   between each gene and all fungal species that had differential
   abundance between the MEF and MED groups, applying FDR to correct the
   correlation P values. Genes significantly correlated with at least one
   fungus were retained (FDR < 0.05). Set 2: DEGs. We identified DEGs
   between the MED and MEF groups (|logFC| > 1, FDR < 0.05), which were
   closely related to the fungal community ecology. Set 3: Prognosis
   genes. We ensured that each gene could predict the OS of patients (Cox
   P < 0.05). The intersection of genes from the three gene sets was
   considered to effectively characterize the fungal features and served
   as input for our model.

   We calculated the C-index of each model in the training and validation
   sets. The highest mean C-index was used as the basis for model
   selection. Subsequently, on the basis of the lastly selected model, we
   calculated the IMS score or IMRGES score for each patient. We then
   determined the optimal cutoff value using the survminer package,
   dividing all patients into high-score and low-score groups. We
   conducted a meta-analysis using Onlinemeta to investigate the impact of
   IMRGES on OS in patients receiving immunotherapy across six cancer
   types. Because of significant heterogeneity observed among the cohorts
   (P < 0.65, I^2 = 59%), a random-effects model was used.

Rationality-congruence assessment

   For the seven fungi in IMS, we constructed the framework of RCA by
   integrating the results of LEfSe, survival analysis, pathway analysis,
   and ML (fig. S5A). This framework allowed for the comprehensive
   screening of key species using a cumulative scoring method based on the
   consistency of multiple evidence. We combined the four analysis results
   in pairs, yielding a total of six combinations for testing. To evaluate
   the performance of fungi in the paired test combinations, we proposed
   the principle of “directionality” consistency. The directionality
   referred to the inference of the potential direct or indirect impact of
   fungi on prognosis based on the findings of previous studies and their
   performance in the four analyses (tests). In any test, if a fungus
   showed an unfavorable impact on survival, its directionality in that
   test is considered unfavorable and scored as −1. Conversely, if it has
   a favorable impact on prognosis, it is scored as +1.

   Specifically, in LEfSe, patients in group B had significantly worse
   prognosis compared to those in group A. If a fungus was enriched in
   group B, it was considered to have a potential unfavorable impact on
   prognosis, and the score for this test was recorded as −1. Conversely,
   fungi enriched in group A were scored as +1. If not enriched in either
   group, it was assigned a score of 0. For survival analysis, we directly
   referred to the HR of the fungi, where HR > 1 could be defined as a
   risk factor and scored as −1. In pathway analysis, we mainly considered
   the impact of lipid catabolism pathways on survival. In RCC, if any
   lipid catabolism pathways are significantly down-regulated, it may lead
   to lipid metabolism imbalance and lipid accumulation, which often
   indicates poor prognosis ([250]55). Therefore, if a fungus was
   negatively correlated with any lipid catabolism pathway, we inferred
   that the higher its abundance, the more likely lipid accumulation
   occurred in the tumor microenvironment, and the more unfavorable it is
   for patient survival, scoring −1. Conversely, if a fungus was
   positively correlated with lipid catabolism pathways, it was scored +1,
   and if it had no significant correlation, it was scored 0. In ML,
   patients in the high–IMS score group showed significantly worse
   prognosis compared to the low–IMS score group. We made judgments based
   on the coefficient of the fungus in the model. If the coefficient > 0,
   it indicated a positive contribution to the model and may be associated
   with worse prognosis, scoring −1.

   The score of a single paired test combination was the product of the
   scores of the two tests (+1, 0, or −1). A score of +1 indicates that
   the two tests have consistent directionality, meaning the paired test
   is passed. Species that passed all paired test combinations,
   accumulating a score of +6, were considered to have passed the RCA.

DNA extraction and qPCR quantification

   Approximately 200 mg of tissue was collected from each ccRCC tumor mass
   for DNA extraction. Total DNA was extracted using the QIAamp PowerFecal
   (Pro) DNA Kit according to the manufacturer’s instructions. Briefly,
   tumor tissues were thoroughly ground in liquid nitrogen using a sterile
   mortar and pestle. Subsequently, C1 buffer was added, and the samples
   were vortexed at maximum speed for 10 min. The samples were then
   centrifuged at 15,000g at 4°C for 10 min, after which the supernatant
   was discarded. DNA was further purified according to the manufacturer’s
   instructions. The quality of the extracted DNA was assessed using a
   NanoDrop spectrophotometer, and all samples had the ratio of the
   absorbance values of 260 nm versus 280 nm (A[260]/A[280] ratios)
   between 1.8 and 2.0. The study was designed with three experimental
   groups: ccRCC tumor blocks obtained under sterile postsurgical
   conditions, environmental controls, and nontemplate controls (detection
   of contamination of the PCR reagents). All procedures were performed
   under strict sterile conditions.

   For qPCR, we used the Takara Premix Ex Taq (Takara-#RR390A) kit to
   detect the presence of A. tanneri within ccRCC samples. The qPCR
   reaction program was as follows: initial denaturation at 94°C for 10
   min, followed by 40 cycles of denaturation at 94°C for 15 s and
   annealing/extension at 60°C for 60 s. For each sample, three replicate
   wells were prepared with 15 ng of DNA per well. The mean CT value of
   the triplicates was used as the CT value for A. tanneri in that sample.

Multiplex immunofluorescence

   Tissue sections were treated with EDTA antigen retrieval buffer. After
   circling with an immunohistochemistry pen, the sections were immersed
   in a 3% hydrogen peroxide solution for 25 min at room temperature in
   the dark. Subsequently, bovine serum albumin was added and incubated
   for 30 min, followed by the addition of the primary antibody and
   overnight incubation at 4°C. After washing with PBS, the corresponding
   fluorescent-labeled secondary antibody was incubated for 50 min. This
   process was repeated for the other three antibodies.
   4′,6-Diamidino-2-phenylindole (DAPI) staining was added at room
   temperature and incubated for 10 min. Analysis was performed using
   ImageJ (v 1.54g). A nuclear segmentation algorithm was applied on the
   DAPI channel to identify and count cells. The number of positive and
   negative cells for each marker was quantified.

Fluorescence in situ hybridization

   Formalin-fixed paraffin-embedded sections were immersed in
   diethylpryrocarbonate water and hydrolyzed with lysozyme (10 mg/ml) for
   30 min to disrupt the fungal cell walls. Subsequently, probes targeting
   the fungal 28S rRNA gene (AGGTGCCCCGAAGGGCGCTA) and A. tanneri–specific
   probe (TAGCGCCCTTCGGGGCACCT) were labeled with Cy3 fluorophore
   (excitation wavelength: 555 nm, emission wavelength: 590 nm) ([251]28).
   The intensity of the fluorescent signal from A. tanneri was quantified
   using ImageJ (v 1.54g).

Statistical analysis

   Comparisons of continuous variables between two groups were performed
   using the Wilcoxon rank-sum test. The chi-square test was used to
   compare categorical variables. The correlation between two continuous
   variables was evaluated using the Spearman correlation test. Cox
   regression and Kaplan-Meier analyses were conducted using the survival
   package. The ROC curves for the ML models were obtained using the pROC
   package. The ggplot2 package was used for graphical visualization. All
   statistical analyses and plotting were performed in R (version 4.2.3).
   All statistical tests were two-sided, and a P value less than 0.05 was
   considered statistically significant.

Acknowledgments