Pan-cancer analyses reveal cancer-type-specific fungal ecologies and bacteriome interactions

Abstract

Cancer-microbe associations have been explored for centuries, but cancer-associated fungi have rarely been examined. Here, we comprehensively characterize the cancer mycobiome within 17,401 patient tissue, blood, and plasma samples across 35 cancer types in four independent cohorts. We report fungal DNA and cells at low abundances across many major human cancers, with differences in community compositions that differ among cancer types, even when accounting for technical background. Fungal histological staining of tissue microarrays supported intratumoral presence and frequent spatial association with cancer cells and macrophages. Comparing intratumoral fungal communities with matched bacteriomes and immunomes revealed co-occurring bi-domain ecologies, often with permissive, rather than competitive, microenvironments and distinct immune responses. Clinically focused assessments suggested prognostic and diagnostic capacities of the tissue and plasma mycobiomes, even in stage I cancers, and synergistic predictive performance with bacteriomes.

Keywords: biomarkers; cancer; fungi; liquid biopsy; metagenomics; metatranscriptomics; microbial interactions; tumor microbiome; tumor mycobiome.

Graphical abstract

Figure 1

Fungal nucleic acids exist in human cancers (A) Table of all studied samples. (B) Fungal DNA abundance in WIS cohort quantified by 5.8S qPCR. Blue bars show medians. Values clipped at 1,000. One-sided t tests between tumor types and extraction controls (n = 89, far left): paraffin controls (n = 48), p = 5.8 × 10⁻⁴; GBM (n = 25), p = 3.1 × 10⁻⁵; melanoma (n = 31), p = 2.1 × 10⁻⁷; colon (n = 19), p = 6.6 × 10⁻⁵; lung (n = 56), p = 2.1 × 10⁻⁶; ovary (n = 26), p = 4.2 × 10⁻⁶; pancreas (n = 25), p = 4.5 × 10⁻¹⁰; bone (n = 25), p = 0.014; and breast (n = 54), p = 1.5 × 10⁻⁵. All p values have an FDR of ≤0.2. (C) Percentage of fungal or bacterial reads in TCGA primary tumors versus total reads. Sample sizes inset in blue, and vary slightly when samples had only bacterial counts. Two-sided Wilcoxon tests for each cancer type; ∗∗∗∗p ≤ 0.0001; ∗∗∗p ≤ 0.001; ns, not significant. Boxplots show median, 25^th and 75^th percentiles, and 1.5 × interquartile range (IQR). See Data S1.2F for paired analysis. (D) Phylogenomics of TCGA-derived fungal bins >85 kbp using Benchmarking Universal Single-Copy Orthologs (BUSCO) against NCBI fungal genomes.

Figure 2

Visualization of fungi in human cancer tissue (A) Table summarizing percent of tumor microarray cores from five cancer types with positive fungal staining of α-β glucan, α-Aspergillus, and 28S rRNA FISH probes, and their localization. (B) Representative stained tumor microarrays from five cancer types using hematoxylin and eosin (H&E), antibodies against β-glucan, Aspergillus, CD45, CD68, CD8, and by FISH probes against fungal 28S rRNA sequences. Negative controls for all these cores are in Data S2.3. Scale bars shown. Squares in H&E images demarcate areas presented at higher magnification. PDAC, pancreatic adenocarcinoma; FISH, fluorescence in situ hybridization.

Figure 3

Different cancer types exhibit distinct mycobiomes (A) Fungal and bacterial species richness for WIS and TCGA cohorts. NC, negative controls. Boxplots show median, 25^th and 75^th percentiles, and 1.5 × IQR. t test p values inset on plots. (B) Scatter plot demonstrating significant Spearman correlations (ρ) and p values between fungal and bacterial richness in four tumor types shared between WIS and TCGA cohorts and no correlation in negative extraction controls. Linear regression lines and 95% confidence intervals shown. (C) Rarefaction plot of the number of species detected in the WIS cohort per tumor type with 100 random subsamples per number of samples. Mean and SD shown. Extraction and paraffin controls were grouped together. (D) Fungal beta diversity analyses using robust Aitchison PCA (Martino et al., 2019) on decontaminated mycobiome data from TCGA MD Anderson primary tumor (WGS) samples (n = 259, 8 cancer types). Permutational multivariate analysis of variance (PERMANOVA) statistics (999 permutations) shown on plot. (E) Mean relative abundance bar plots at class-level phylotypes across WIS tumor types. Colors correspond to fungal class. (F) Unsupervised hierarchical clustering of fungal prevalence in the WIS cohort using species that appear in ≥10% of samples in ≥1 tumor/NAT/normal tissue types. Values represent Z scores per row. Amplicon sequence variants (ASVs) without species level classification were aggregated by the lowest classification they received. (G) Principal coordinate analysis (PCoA) of Jaccard dissimilarities between composite fungal species profiles across tissues.

Figure 4

Establishing pan-cancer mycotypes through mycobiome-bacteriome-immunome interactions (A) Co-occurrence analyses of WIS-overlapping TCGA fungal and bacteria genera (Table S7.5), and TCGA immune cell compositions (Thorsson et al., 2018) using MMvec (Morton et al., 2019a). Hierarchical clustering linkage information identified three distinct clusters (“mycotypes”) associated with groups of fungal genera: F1, F2, and F3. (B) Log-ratios of fungal mycotype abundances across TCGA cancer types, revealing significantly differing values (one-way ANOVAs). (C) Varying mycotype immune log-ratios across pan-cancer immune subtypes (Thorsson et al., 2018). C1, wound healing; C2, IFN-γ dominant; C3, inflammatory; C4, lymphocyte depleted (but with second most macrophages); C5, immunologically quiet (but with most macrophages); C6, TGF-β dominant. Table S7.3 shows pairwise log-ratio comparisons across all immune subtypes. (D) Significant associations with overall survival in 20 cancer types based on the F1/F2 fungal mycotype log-ratio (left) or mycotype immune log-ratios (F1/F3, middle; F2/F3, right). Table S7.6 shows the sample sizes above and below the medians. Note: the C3 immune subtype has the best prognosis. (B and C) Boxplots show median, 25^th, and 75^th percentiles and 1.5 × IQR.

Figure 5

Machine learning (ML) analyses reveal cancer-type-specific tumor and blood mycobiomes (A) One-cancer-type-versus-all-others predictions on Harvard Medical School tumors (HMS, n = 876). (B) Negative control analyses for (A) using scrambled metadata or shuffled samples. All one-cancer-type-versus-all-others performances are aggregated. ^∗∗∗∗ q < 0.001; ns, not significant. (C) Multi-class pan-cancer discrimination among TCGA WGS tumor samples using WIS-overlapping features across 500 independent folds (50 iterations of 10-fold CV). (D) Aggregated one-cancer-type-versus-all-others ML performance in WIS cohort tumors. (E) One-cancer-type-versus-all-others predictions using batch-corrected, TCGA primary tumor data (n = 10,998). (F) One-cancer-type-versus-all-others predictions using HMS blood samples (n = 835). (G) Multi-class pan-cancer discrimination among TCGA WGS blood samples using WIS-overlapping features across 500 independent folds (50 iterations of 10-fold CV). (H) One-cancer-type-versus-all-others predictions using batch-corrected, TCGA blood data (n = 1,771). (A, E, F, and H) Area under ROC curve (AUROC) and area under precision-recall curve (AUPR) measured on independent holdout folds (10-fold cross-validation [CV]) to estimate averages (dots) and 95% confidence intervals (brackets). “High coverage,” 31 fungal species with ≥1% aggregate genome coverage; “∩ Weizmann,” 34 WIS-overlapping fungal species; “decontaminated,” 224 decontaminated fungal species. Horizontal lines denote null AUROC or AUPR. (B, C, D, and G) Two-sided Wilcoxon tests with Benjamini-Hochberg correction. Boxplots show median, 25^th, and 75^th percentiles and 1.5 × IQR.

Figure 6

Clinical utility of cancer mycobiomes (A and B) Differential prevalence of fungal taxa in WIS breast tumors by (A) age or (B) human epidermal growth factor receptor 2 (HER2) status. (C) Kaplan-Meier survival probability of WIS breast cancer patients positive (n = 11) or negative (n = 69) for Malassezia globosa. p value from log-rank test. (D) Fungal richness in WIS lung tumors by smoking status. Boxplots: median, 25^th and 75^th percentiles, and 1.5 × IQR. (E) Differential prevalence of fungi in WIS lung tumors by smoking status. (F) Kaplan-Meier plot demonstrating progression free survival (PFS) probability in WIS ovarian patients positive (n = 9) or negative (n = 45) for Phaeosphaeriaceae family. p value from log-rank test. (G) Differential prevalence of fungi in WIS melanoma tumors by immune checkpoint inhibitor response. (H) Treatment-naive pan-cancer versus healthy discrimination in the Hopkins plasma cohort across all database hits (red, 7,418 features), WIS-overlapping fungi and bacteria (blue, 287 species), or decontaminated fungi (orange, 209 species). (I) Per cancer type versus healthy discrimination in the Hopkins cohort with 10-fold cross-validation. The “top 20 fungi” (green) are derived from pan-cancer versus healthy machine learning model. Dots and brackets represent average performance and 95% confidence intervals, respectively. Horizontal lines (gray or colored) denote null AUROCs and AUPRs. (J) Stage I pan-cancer versus healthy discrimination in the Hopkins cohort with equivalent feature sets and colors as (H). (K) Pan-cancer versus healthy controls discrimination in the Hopkins (purple) and UCSD (teal) plasma cohorts using the “top 20 fungi” features. (A, B, E, and G) p values calculated by Fisher’s exact test. (H, J, and K) 10-fold cross-validation repeated ten times. Mean performance with 99% confidence intervals (colored ribbons) and gray or lightly colored lines each denoting single repeats.