Gut Mycobiota Dysbiosis Is Associated with Melanoma and Response to Anti–PD-1 Therapy Open Access

Melanoma is a malignant skin cancer that arises from melanocytes, the pigment-producing cells in the skin. The incidence of melanoma has been increasing globally over the past few decades, and it is now one of the most common cancers worldwide (1). Although the major risk factors for melanoma are well-established, such as exposure to UV radiation, genetic predisposition, and skin pigmentation type (2), there are still some gaps in understanding the whole picture of melanoma progression and the efficacy of its treatments.

Recent studies have suggested that the microbiome plays a role in melanoma progression (3). Changes in microbiome diversity have been observed in melanoma patients (4). It has been shown that the microbiome can influence the development of melanoma by affecting the metabolism of host cells. For example, some microorganisms have been found to produce compounds that increase melanocyte proliferation, promoting melanoma development (5, 6). Additionally, research has shown that the microbiome can modulate the immune response to melanoma cells, influencing the progression of the disease, and the presence of certain bacteria has been found to be positively or negatively correlated with melanoma progression (7, 8). Although a recent study suggested that the mycobiota present in the gut may also have an impact on melanoma progression (4), the role of gut mycobiota in cancer development and particularly in melanoma still is severely understudied compared with the role of bacteria.

Gut fungi may have a profound effect on the human body. Many of them are believed to have a beneficiary effect, such as Saccharomyces cerevisiae var. boulardii (9), whereas others, such as certain species of Candida and Aspergillus, are considered harmful. The impact of fungi is associated with the production of a wide range of metabolites that strongly affect metabolic pathways (10). Also, components of the fungal cell wall (including chitin, |$\beta $|-glucan, and mannan) can interact with the host immune system to induce proinflammatory responses (11, 12). As inflammation is considered a key contributor to cancer development (13, 14), it can be hypothesized that also in the case of cutaneous melanoma, intestinal inflammation associated with altered microbial profiles may contribute to the maintenance of a systemic proinflammatory, and potentially pro-oncogenic, state (15, 16).

This study aimed to explore the relationship between gut mycobiota and melanoma, examining changes in gut mycobiome structure in response to anti–PD-1 therapy and mycobiota association with response to such treatment. We believe that our research will bring new insight into the intricate connection between gut fungi and melanoma progression, particularly in the context of anti–PD-1 therapy and its efficiency.

This study involved human participants and was approved by the Bioethical Commission of the Karol Marcinkowski University of Medical Sciences, Poznan, Poland (resolution No 316/22, passed on 14 April 2022, as an extension of resolution No. 485/19, passed on 11 April 2019). The study was conducted in accordance with the Declaration of Helsinki. Participants gave written informed consent for the collection of stool samples and study information before taking part. All methods were carried out in accordance with relevant guidelines and regulations.

Patients or the public were not involved in the design, conduct, reporting, or dissemination of the research.

Stool samples were collected from 86 melanoma patients (referred to as melanoma or oncological) and 134 healthy participants (referred to as controls) selected from the Polish Microbiome Map project (ClinicalTrials.gov study identifier: NCT04169867).

Among the 86 melanoma patients were 62 who had histologically confirmed unresectable stage III or stage IV cutaneous melanoma and were enrolled in treatment with anti–PD-1 therapy (nivolumab or pembrolizumab) as a part of the Ministry of Health (Poland) drug program (17), and 24 patients who had clinically stable long-term melanoma (stages III and IV) and were receiving genetic whole-cell therapeutic melanoma vaccine AGI-101H injections every month (ETAM2-5 Trial, EudraCT Number 2008-003373-40; refs. 18, 19). All patients were recruited to the study at the Department of Medical and Experimental Oncology, Heliodor Swiecicki Clinical Hospital, Poznan University of Medical Sciences (Poznan, Poland), from June 2018 to December 2021.

Clinical information, including tumor stage and serum lactate dehydrogenase (LDH) concentration, was collected from the medical records. Response to PD-1 therapy was assessed according to the response evaluation criteria in solid tumors (RECIST) v.1.1 criteria. Responders (R) were defined as patients with complete (CR) or partial (PR) responses, whereas nonresponders (NR) were defined as patients with stable disease (SD) or progression (PD). Twenty-seven R and 34 NR patients were analyzed during this study. Additional, independent classification categorized patients into those who exhibited clinical benefit (CB) from immune-checkpoint inhibitor (ICI) therapy (CR, PR, and SD) and those who did not (NB, consisting of PD). 37 CB and 24 NB patients were analyzed during this study.

Healthy participants were broadly representative of the overall Polish adult (>18 years old) population concerning age, health, education, and marital status. Healthy samples were taken from 2020 to 2021.

Detailed information about the patients and healthy controls, including information about the number of samples obtained from each individual, is presented in Supplementary Table S1.

All samples belonging to healthy participants who marked autoimmune diseases in the questionnaire were filtered out (74 samples). To account for differences in the number of samples taken from patients and healthy study participants, thus ensuring that the compared data sets are properly balanced while keeping the demographic characteristics within the compared groups as similar as possible, particularly in terms of age and body mass index (BMI), only selected samples taken from healthy patients were included in the analysis. The selection of healthy samples was as follows. First, healthy and melanoma data sets were stratified according to age (≤25, 26–30, 31–40, 41–50, 51–60, 61–70, 71–80, >80) and BMI (BMI < 18.5 underweight, 18.5 ≤ BMI < 25 normal, 25 ≤ BMI < 30 overweight, 30 ≤ BMI obesity). Then, the age and BMI structures of the melanoma data set were visually assessed. As a result, all healthy samples belonging to subjects younger than 31 were removed, as none of the samples in the patient data set belonged to patients below the range of 31 to 40. Then, in order to match the age and BMI structure of the melanoma data set, an appropriate number of samples were randomly drawn from the healthy data set taking into account data stratification: four samples from the healthy, obese, or overweight participants in a range of 31 to 40 years of age, and 40 samples from the healthy, obese, or overweight participants in a range of 41 to 50 years of age. Due to a scarcity of such samples, all samples from participants older than 50 were included in the analysis in order to ensure the closest possible similarity of demographic characteristics of healthy and melanoma subjects. As a result, samples from 134 healthy participants were included for further analysis, thus constituting the control group.

The sampling, DNA isolation, sequencing libraries preparation, and whole-genome high-throughput sequencing were performed according to the protocol described in our previous studies (20, 21). The details are provided below.

Fecal samples (approximately 1 g) were self-collected by all donors into vials containing 3 mL of RNAlater Stabilization Solution (Invitrogen, Thermo Fisher Scientific) and delivered by courier within 24 hours to the laboratory, where the samples were anonymized and stored at 4°C for up to one week. Each sample was homogenized. Tubes were centrifuged at 14,000 × g for 5 minutes, the supernatant was discarded, and residues were transferred to a −20°C freezer for storage until DNA extraction (usually a few weeks).

The frozen stool samples were thawed on ice, and DNA was extracted from them using a DNAeasy PowerSoil Pro Kit (Qiagen) according to the manufacturer's instructions. The following protocol adjustments were applied. The liquid phase of stabilized stool samples was thoroughly discarded to remove high salt content that may interfere with a subsequent DNA purification step. Next, the stabilized stool samples (250 mg) were bead-beaten in PowerBead Pro tubes (Qiagen, cat. no. 19301) containing a mix of zirconium beads of different diameters using a Mixer Mill MM400 (Retsch) for 15 minutes at 25 Hz. To remove RNA and increase DNA yield, each sample was incubated with 5 μL RNase (10 mg/mL concentration; A&A Biotechnology) at 60°C for 10 minutes. The DNA quality was verified with agarose gel electrophoresis. The final DNA concentration was measured by a Nanodrop ND-1000 spectrophotometer (Thermo Fisher Scientific).

Libraries were constructed with the TruePrep DNA library prep kit V2 for Illumina TD501 (Vazyme Biotech, China) according to the manufacturer's protocol. 500 ng of stool-extracted DNA was used for each library preparation. In the library amplification step, six PCR cycles were applied. Library concentration was measured using a Qubit fluorometer and Qubit DNA HS Assay Kit (Thermo Fisher Scientific). The quality of libraries and fragment distribution were analyzed using a Bioanalyzer 2100 and DNA 1000 Kit or High Sensitivity DNA Kit (Agilent Technologies), depending on the obtained library quantity.

Purified libraries were stored for a few weeks at −20°C until sequencing.

Before sequencing, all libraries were thawed on ice and normalized to a final concentration of 10 nmol/L. Libraries with distinctive index combinations were pooled together and diluted with EB Buffer (Qiagen) to obtain a mix of 2 nmol/L libraries, according to Protocol A: Standard Normalization Method for the NextSeq system (Illumina). Whole-genome sequencing was performed with NextSeq 550 (Illumina) using High Output Kit v2.5 reagents (Illumina); approximately 10 million 150 bp paired-end reads were generated per library. Neither human DNA sequence depletion nor microbial or viral DNA enrichment was performed.

Demultiplexing was run on the raw BCL intensity file with the bcl2fastq tool (22) for base calling and separating the reads from different samples. To assess the quality of the sequencing procedure, we generated quality control reports with FastQC (23) and MultiQC (24). We preprocessed the raw fastq reads with cutadapt (25) using the following procedure: we trimmed the adapter sequences (based on TruSeq adapter sequences) and poly-G tails observed in the data characteristic of the two-channel sequencing technology of NextSeq. We also filtered out reads shorter than 140 bases to remove the bias in taxonomy profiling that could emerge from the shorter sequences. The remaining reads were subjected to further analysis.

Profiling of the fungal microbiome was performed via Kraken2 (26, 27), followed by Bracken (28).

Statistical analysis of the data was performed in R language environment v. 4.1.1. The Kruskal–Wallis test (for continuous data), Chi-square test, and Fisher exact test (for categorical data) were used to study whether there were significant differences between analyzed groups (Oncological vs. Controls, R vs. NR, CB vs. NB) in terms of variables, accepting a confidence level of 95%.

Survival analysis was performed using the Kaplan–Meier method from R library survival. The Kaplan–Meier survival curves were visualized using the ggsurvplot function from the survminer package. Survival curves were compared using the log-rank (Mantel–Cox) test using the survdiff available in the survival package. Cox regression was used to examine associations between covariates of interest and survival outcomes using the coxph function from the survival R package. The level of fungi abundance (low or high) was determined using the maximally selected rank statistics using the surv_cutpoint function from the survminer R package.

For α-diversity, the following metrics were calculated: richness, evenness (Pilou evenness), and effective diversity (Shannon's Index, diversity function, R package vegan; refs. 29, 30). Fungal richness (species count) was calculated by counting the number of species detected at least once in a given sample. Pilou evenness was calculated by dividing the diversity (Simpson's Index, diversity function, R package vegan) by |$log( S )$|⁠, where |$S$| is the total number of species. Data were log-ratio centered to analyze the structure of mycobial communities across samples (β-diversity, clr transformation, clr function, compositions R package). Pseudocount of 0.005 was added to the species raw count prior to distance calculation to avoid zeros, as clr transformation cannot be used with non-positive data.

Statistical testing of α-diversity was performed using the Kruskal–Wallis test. Multiple pairwise comparisons between groups were performed with the pairwise Wilcoxon rank-sum tests (pairwise.wilcox.test function, stats R package, false discovery rate [FDR] <5%). For β-diversity significance testing, permutational multivariate analysis of variance using distance matrices was conducted (adonis (31), Euclidean distance, vegan R package, 999 permutations, FDR < 5%), accompanied by the Wilcoxon rank-sum test with continuity correction. For the data with more than two levels, β-diversity significance testing was performed with the pairwise implementation of the aforementioned methods, pairwise adonis (custom implementation), and pairwise Wilcoxon rank-sum tests (pairwise.wilcox.test function, stats library). The principal component analysis (PCA, R package stats) was used to illustrate how the community structure of individual samples across the study population differs.

To identify fungal species contributing to community compositional dissimilarities data were adjusted by sex, BMI groups, and age groups (as potential confounding factors; ref.32), and five different methods were utilized: Maaslin2 (ref. 33; Microbiome Multivariable Associations with Linear Models, Maaslin2 R package from Bioconductor), linDa (ref. 34; Linear Discriminant Analysis, MicrobiomeStat R package with custom changes), Linear Decomposition Model (ref. 35; LDM, LDM R library), pairwise Wilcoxon rank-sum test (pairwise.wilcox.test function, stats library), and Welch two-sample t test (ggbetweenstats function from ggstatplot R package; ref. 36). The result was considered confirmed when at least two out of five methods agreed. Correlation analysis was performed using Pearson correlation coefficient (cor.test function from stats R package).

Data generated in this study are available within the article and its Supplementary Data Files or from the corresponding author upon reasonable request. The raw data of the healthy human cohort are available in the NCBI SRA repository (https://www.ncbi.nlm.nih.gov/sra) under accession number PRJNA850529. The raw data of the melanoma cohort are available in the NCBI SRA repository (https://www.ncbi.nlm.nih.gov/sra) under accession number PRJNA972625.

The study cohort consisted of melanoma and healthy subjects. Eighty-six patients with unresectable stage III or stage IV cutaneous melanoma were recruited for the study. Stool samples from the enrolled melanoma patients treated with anti–PD-1 inhibitors were analyzed at three time points: before treatment, after 3 months of treatment, and after 12 months of therapy (61, 36, and 17 samples, respectively). For one patient, we have an additional fourth sample, which was obtained in the sixth month after treatment. Stool samples from patients receiving vaccine therapy were analyzed at two time points: before the monthly vaccine administration and 6 to 10 days after the next vaccine dose (23 and 23 samples, respectively). In total, we obtained 164 samples from melanoma patients. Among the melanoma patients treated with the anti–PD-1, 37 patients exhibited CB from ICI therapy, and 24 did not (NB); additionally, of these, 27 patients responded to therapy (R), whereas 34 did not (NR). For reference, out of 859 nonmelanoma samples described in our previous study (32), we selected 134 healthy samples that served as control (details of healthy sample selection are provided in the Supplementary Materials). The baseline clinical data and characteristics of patients and controls are shown in Table 1; Supplementary Fig. S1; Supplementary Table S1.

The average serum LDH level was significantly higher in NR compared with R (and similarly in NB compared with CB, P < 0.05, Fig. 1A). The number of patients with elevated serum LDH level (LDH > 1.5×, which corresponds to >375 U/L) was also significantly higher in R and CB (P < 0.05). Moreover, patients with an elevated level of serum LDH had a lower chance of surviving longer without melanoma progression (LDH > 1.5×, P < 0.05, Fig. 1B). The median progression time was much lower for patients with an increased versus normal level of serum LDH (∼2.25 vs. ∼16 months, respectively; Fig. 1C). The results from Cox's proportional hazards (CPH) model showed that patients with serum LDH > 1.5× (considering age group, sex, and BMI group) were 5.87 more likely to experience melanoma progression than patients with the normal LDH level (P < 0.05, Fig. 1D).

Taking a look at changes in the fungal species, gut fungal diversity (Shannon diversity) of patients with an increased serum LDH level was significantly higher than in patients with normal serum LDH levels (P = 0.02, Supplementary Fig. S2A). When analyzing individual fungal species, we observed that anti–PD-1 patients with elevated serum LDH had a higher abundance of Sugiyamaella lignohabitans and Neurospora crassa (Supplementary Table S2; Supplementary Fig. S2B). We also observed a positive correlation between increased N. crassa abundance and a rising serum LDH level (P < 0.01, R = 0.37).

The age, sex, and BMI of the patient had no significant effect (P > 0.05) on melanoma progression probability. However, when we divided patients into groups based on their BMI (BMI < 18.5 underweight, 18.5 ≤ BMI < 25 normal, 25 ≤ BMI < 30 overweight, 30≤BMI obesity), we found that being overweight affects melanoma progression probability (P < 0.05). Overweight patients were 2.79 more likely to experience melanoma progression than patients with normal BMI (Fig. 1D). Moreover, when analyzing differences in the gut mycobiota between normal and overweight patients, overweight subjects had a higher abundance of Nakaseomyces genus and Candida glabrata species (Supplementary Table S3; Supplementary Fig. S3).

We characterized the taxonomical structure of each sample. Fungi were detected in 148 (90.24%) melanoma samples and 115 (85.82%) healthy samples. Fifty-four species were detected in the melanoma and healthy samples (53 and 47 species, respectively; Fig. 2A and B), belonging to 36 genera (36 and 32, respectively, Fig. 2C; Supplementary Figs. S4 and S5). 93.62% of the species identified in the healthy samples and 92.45% in the melanoma samples belonged to the Ascomycota phylum. Ascomycota was also the dominant phylum regarding occurrence in particular samples, present in 78.52% of healthy and 83.03% of melanoma samples. Basidiomycota was the next most frequently occurring phylum, present in 42.96% and 38.79% of healthy and melanoma samples, respectively (Fig. 2D). The three most frequently occurring genera in healthy samples were Saccharomyces, Candida, and Sporisorium (present in 40.00%, 27.41%, and 25.19% of the samples, respectively, Fig. 2E). In melanoma samples, the most frequently occurring genera was Candida (present in 44.85% of the samples), followed by Saccharomyces (31.52% of melanoma samples) and Sporisorium (25.19% of melanoma samples). At the species level, the most prevalent species in healthy subjects was Saccharomyces cerevisiae (present in 39.55% of samples), followed by Sporisorium graminicola and Candida albicans (both present in 24.63% of samples; Fig. 2F). In melanoma patients, the occurrence of C. albicans was much higher (40.85% of samples) than in healthy samples, whereas the occurrence of S. cerevisiae and S. graminicola was lower (31.10% and 18.90% of samples, respectively) than in the healthy cohort.

To assess the importance of differences between melanoma and healthy samples, we first analyzed the differences in richness and α-diversity (Shannon and Pilou diversities) and β-diversity. Although the difference in fungal richness between melanoma and healthy samples was not statistically significant, the maximum number of species identified in a melanoma sample was much higher than in any sample in the healthy cohort (11 vs. 22 fungal species in a single sample; Supplementary Fig. S6). Neither richness nor Shannon nor Pilou diversity significantly differed between healthy and melanoma samples. However, when analyzing β-diversity, melanoma samples differed considerably from the control cohort (P < 0.05, Adonis, Fig. 3A). PCA showed that melanoma samples were skewed toward three fungi species from the Candida, whereas healthy samples were skewed toward S. cerevisiae (Fig. 3A).

To assess differences in the abundance of detected fungal species between healthy and melanoma patients, we used five statistical methods, and only results significant in at least two of the five methods were considered valid. Consistent with the PCA view, when analyzing the differences at the individual genus and species level, gut mycobiota of melanoma subjects proved to be depleted in Saccharomyces (P < 0.05, log₂ fold change (LFC) = −1.31), more specifically, S. cerevisiae (P < 0.05, LFC = −1.31, Fig. 3B; Supplementary Table S4). Additionally, melanoma patients had a lower abundance of Debaryomyces hansenii (P < 0.05, LFC = −1.53) and Sugiyamaella lignohabitans (<0.05, LFC = −0.07, Fig. 3B; Supplementary Table S4). On the other hand, melanoma samples were enriched in Candida (P < 0.05, LFC = 2.04), especially Candida dubliniensis and C. albicans (P < 0.05, LFC = 0.52 and LFC = 0.74, respectively, Fig. 3B; Supplementary Table S4). Analyses at the species level also showed that melanoma samples contained a significantly higher abundance of N. crassa (P < 0.05, LFC = 0.72, Fig. 3B; Supplementary Table S4) and Thermothielavioides terrestris (P < 0.05, LFC = 0.05, Fig. 3B; Supplementary Table S4).

To further assess changes in gut mycobiota of melanoma patients during treatment with anti–PD-1 therapy, we compared the samples taken at the three time points (before treatment, in the third month of treatment, and after 12 months) and healthy samples (Fig. 4A). None of the general α-diversity characteristics (richness, Shannon, and Pilou diversity) differed within selected time points of the treatment and between them and healthy samples. At the β-diversity level, only anti–PD-1 samples taken before treatment differed from healthy samples (pairwise adonis, q < 0.05; Fig. 4A).

The majority of statistically significant changes at the fungal species level were observed between the healthy samples and melanoma samples taken before anti–PD-1 therapy (Supplementary Table S5). Samples from patients before the treatment had a lower abundance of Ustilago maydis and D. hansenii, although the latter was confirmed by only one method, and a higher abundance of Aspergillus fumigatus, Scheffersomyces stipitis, C. albicans, C. dubliniensis, and N. crassa (confirmed by only one method; Fig. 4B; Supplementary Table S5). The abundance of C. albicans gradually increased during the observation timeline, reaching the highest log₂ of fold change after 12 months of the therapy (LCF before treatment = 1.90, LFC after 12 months = 3.19; P < 0.05, Supplementary Table S5, Fig. 4B; Supplementary Fig. S7). The abundance of C. dubliniensis was also increased when comparing healthy samples and samples taken in the third month of the therapy; however, in this case, the difference in abundance decreased over the timeline of observation, taking the healthy samples as a reference (Supplementary Table S5; Fig. 4B; Supplementary Fig. S7). In the case of D. hansenii, we observed that its abundance dropped even more between starting anti–PD-1 therapy and the third month of anti–PD-1 treatment compared with healthy samples (Supplementary Table S5; Fig. 4B; Supplementary Fig. S7). For A. fumigatus, we also observed a decrease in abundance in samples taken in the third month of therapy when compared with samples before anti–PD-1 administration; however, it should be noted that this difference was reported as significant by only one utilized statistical method (Supplementary Table S5, Fig. 4B; Supplementary Fig. S7).

The results also revealed that during anti–PD-1 therapy, the level of Malassezia restricta successively increased compared with melanoma samples before therapy and healthy samples (P < 0.05, Supplementary Table S5; Fig. 4B; Supplementary Fig. S7). This trend could be observed even though initially, the abundance of M. restricta was slightly lower in melanoma patients before anti–PD-1 therapy compared with healthy samples (P > 0.05). Additionally, we saw an increase in Zygosaccharomyces rouxii during anti–PD-1 therapy when comparing samples taken after 12 months of anti–PD-1 therapy and before anti–PD-1 injection (P < 0.05, one method) and samples taken after 12 months of anti–PD-1 therapy and healthy samples (P < 0.05; Supplementary Table S5, Fig. 4B,). It is worth mentioning that although changes in S. cerevisiae were not indicated as significant by any of the utilized methods, we observed that during anti–PD-1 therapy, its abundance gradually dropped (Fig. 4B; Supplementary Fig. S7).

When analyzing differences during anti–PD-1 therapy between patients with normal BMI and overweight patients, we observed that before starting the anti–PD-1 therapy, overweight subjects had a higher abundance of C. glabrata and a lower abundance of Talaromyces rugulosus than patients with normal BMI (Supplementary Table S6; Supplementary Fig. S8A). C. glabrata level was also elevated in overweight patients in the third month of the anti–PD-1 therapy, whereas C. albicans abundance was lower than in normal subjects (Supplementary Table S6; Supplementary Fig. S8B). There were no statistically significant differences when analyzing changes during anti–PD-1 therapy in patients with normal BMI (Supplementary Table S7). In contrast, in overweight patients, we observed that the level of Yarrowia lipolytica was significantly higher after 12 months of anti–PD-1 therapy when compared with its level before the therapy and in the third month of the ICI treatment (Supplementary Table S8; Supplementary Fig. S9).

Comparison of patients with CB versus NB from anti–PD-1 therapy showed that patients from the CB group had lower gut fungal richness and Shannon diversity than patients from the NB group (t test, P_adj < 0.05, Fig. 5A). Additionally, Tetrapisispora blattae was less abundant in patients who benefited from anti–PD-1 treatment (LFC = −0.93, P < 0.05; Supplementary Table S9; Fig. 5C). When R versus NR were compared, we observed that similar to CB versus NB, R had a lower fungal richness and Shannon diversity than NR (t test, P_adj < 0.05, Fig. 5B). Analysis of individual taxa differentiating R from NR showed a higher abundance of Saccharomyces paradoxus in R (LFC = 1.34, P < 0.05; Supplementary Table S9; Fig. 5D). Figures for other potentially important fungal species, although not significantly differentiating, are presented in Supplementary Figs. S10 and S11.

The results from CPH model showed that there is a significant effect of both M. restricta and C. albicans abundance levels on melanoma progression probability (including effect of LDH in the model, P < 0.05, Fig. 6A). Patients with a high level of M. restricta were 2.2 more likely to experience melanoma progression than patients with a low level of M. restricta, and the increase in the likelihood of progression for patients with a high level of C. albicans was two times the likelihood for patients with a low level of C. albicans. To further assess the importance of the specific fungal species on the survival of the melanoma patients treated with the anti–PD-1 therapy, we performed survival analysis for the selected fungal species. Subjects with a higher abundance of C. albicans and M. restricta showed a trend of shorter survival without melanoma progression, although the differences were not formally significant (P = 0.09 and P = 0.07, Fig. 6B and C). The median progression time was shorter for patients with an increased level of M. restricta (∼5 months) compared with patients with a lower abundance of M. restricta (∼13 months). A greater difference in median progression time also was observed in the case of patients with higher and lower abundance of C. albicans (∼3 vs. 14.5 months, respectively).

Here, we present a cohort study with the main objective of identifying differences between the gut mycobiota of melanoma patients and a healthy cohort, as well as characterizing changes in gut mycobiota during anti–PD-1 therapy and its association with the therapeutic efficacy of ICI therapy in a cohort of metastatic melanoma patients. As most metagenomics studies published thus far have concentrated only on the bacterial part of the microbiota, in this work, we are broadening this perspective and taking a closer look at another essential part of the human gut, fungi.

The influence of the gut mycobiota may be an important factor in the context of cancer progression and the efficiency of immunomodulation therapies. Consistent with Vitali and colleagues (4), we observed that the fungal community structure showed significant differences compared with healthy subjects. The mycobiota of patients with metastatic melanoma were mainly enriched in fungi from the Candida genus (C. dubliniensis and C. albicans), as well as N. crassa, and depleted in S. saccharomyces and D. hansenii. Our analysis showed that although there was no statistically significant difference at the level of richness between healthy subjects and melanoma patients, the maximum number of species identified in the melanoma sample was much higher than in the healthy cohort (11 vs. 22 fungal species). Vitali and colleagues observed a significantly higher richness of fungi in melanoma patients (4). Moreover, we observed many changes in fungal abundance during treatment with anti–PD-1, and identified fungi associated with a response to anti–PD-1 therapy. S. paradoxus was associated with a positive response to anti–PD-1 therapy, whereas T. blattae was associated with a lack of CB during treatment. Moreover, M. restricta and C. albicans were associated with a higher risk of melanoma progression and, therefore, a worse outcome of the anti–PD-1 treatment.

Fungi from the Candia genus, especially C. albicans, which we observed in higher amounts in the melanoma samples, are widely known as opportunistic pathogens. It is known that C. albicans infections are common in cancer patients, whose immune system is often impaired (37). Infections with C. albicans are connected to inflammatory cytokines such as TNFα and IL18, increased levels of which correlate with metastasis of many cancer types. In fact, it has been shown that C. albicans significantly increased the ability of murine B16 melanoma (B16M) cells to colonize the liver (37). In that study, treatment with an antifungal agent prevented the enhancement of hepatic metastasis in fungal-challenged mice, further indicating that C. albicans may be responsible for metastasis in melanoma patients. Moreover, candidalysin, the cytolytic toxin from C. albicans, activates innate epithelial immune responses via epidermal growth factor receptor (EGFR; ref. 38) indirectly through its effect on matrix metalloproteinases (MMP) and EGFR ligands, both independently linked to several cancers (39, 40). C. albicans can also activate epithelial MAPK and extracellular signal-regulated kinase (ERK) signaling pathways, which are related to cancer (41). Apart from the higher abundance of C. albicans in melanoma subjects in general, our study showed that C. albicans abundance increases during anti–PD-1 treatment, but the exact reason remains for further studies. Moreover, in our study, a higher abundance of C. albicans was associated with a higher risk of melanoma progression. The potential usefulness of C. albicans as a biomarker of patients’ prognosis is additionally supported by the fact that the Candida genus, although in saliva, has been found to be significantly associated with poor overall survival of oral squamous cell carcinoma patients in the Sudanese cohort (42).

Elevated expression and activation of a MAPK in response to cell wall synthesis inhibitors are frequently seen in filamentous fungi (43), such as N. crassa, which were significantly more abundant in melanoma samples. Noticeably, MAPK pathway overexcitement is related to cancer (41). Moreover, in N. crassa, regulation of tyrosinase transcription, an enzyme that catalyzes the production of the secondary metabolite l-DOPA melanin, is very similar to that known for mammalian melanoma cells (44). In humans, tyrosinase is a key enzyme for melanin biosynthesis and is associated with melanin hyperpigmentation and melanoma (45). It has also been shown for Neurospora that tyrosinase expression can be induced by nitrogen starvation and during sexual development (46). Taking into account the similarity of tyrosinase transcription pathway regulation between N. crassa and human cells, its induction may provide a link between N. crassa and melanoma progression. Moreover, the positive correlation between increased N. crassa abundance and elevated serum LDH level—LDH is considered a prognostic biomarker of melanoma progression (47)—may also give some hints to the putative mechanism of melanoma development.

M. restricta, which we found to be associated with a higher risk of melanoma progression, was previously reported also to be associated with colorectal cancer (48, 49) and promotion of pancreatic ductal adenocarcinoma (PDA; ref. 50). In the PDA study, fungal elimination with amphotericin B was tumor-protective, whereas repopulation with Malassezia (and not the other fungi) promoted oncogenesis. The contribution of Malassezia-synthesized aryl hydrocarbon receptor ligands to basal cell carcinoma through UV radiation-induced carcinogenesis has also been suggested (51), further supporting a role for M. restricta in skin cancers, and indirectly validating M. restricta as a potential biomarker of worse prognosis in melanoma patients. Aspergillus synthesizes aflatoxin B1 (AFB1), a well-known hepatocellular carcinogen (10), and patulin which also has the potential for causing DNA damage and plausibly facilitating carcinogenesis (52). Of potential relevance in this regard, four species of Aspergillus were enriched in colorectal cancer (49), and we identified A. fumigatus in higher amounts in melanoma patient samples.

In addition to the observed changes in gut mycobiota, our results are consistent with the view that serum LDH level can be an early marker for outcome in patients treated with anti–PD-1 therapy in metastatic melanoma (47). Moreover, in our study, increased serum LDH level was associated with higher gut fungal diversity and increased abundance of S. lignohabitans and N. crassa.

In our studies, we also observed changes in the abundances of species generally regarded as beneficiary, among which there are also probiotic strains (53). We noticed that S. cerevisiae, a common component of the human gut microbiota (32), was less abundant (P < 0.05) in melanoma patients compared with the control. This fungi is able to inhibit pathogenic bacteria and yeast and has been shown to have antioxidant properties (54). Anticancer properties of S. cerevisiae have been reported in many studies (55). Importantly, in vitro and in vivo murine studies have shown that oral administration of |$\beta $|-glucan extracted from S. cerevisiae can suppress the development of melanoma in a dose-dependent manner (56). A similar effect was observed for intravenous administration with |$\beta $|-glucan derived from mutated S. cerevisiae strain in metastatic colon carcinoma (57). Moreover, pretreatment of mice with the same type of |$\beta $|-glucan, two days before inoculation with the tumor cells, enhanced mice survival time. This was probably due to enhanced proinflammatory cytokine production and anticancer activity of peritoneal macrophages, and to increased natural killer (NK) cell cytotoxicity. A clinical trial of advanced breast cancer showed that oral administration of two 10-mg capsules of soluble |$\beta $|-glucan derived from S. cerevisiae enhanced the proliferation and activation of peripheral blood monocytes with no clinical adverse effects (58). This may indicate that S. cerevisiae–derived β-glucan may exert immunomodulatory effects that can help facilitate an answer to available anticancer therapies, a notion further supported by the reported increase in regulator of immune response CD95 expression on CD14⁺ monocytes.

Much less is known about D. hansenii, which we observed to be decreased in melanoma patients. Some of its strains have been studied for potential use as probiotics and may be beneficial (59). Probiotic D. hansenii CBS 8339 yeast enhanced immune responses in mice by upregulating gene expression of proinflammatory cytokines such as INF-γ, IL6, and IL1β (60). In addition, two novel exopolysaccharides produced by D. hansenii DH-1 significantly inhibited the proliferation of Hela, HepG2, and PC-9 cancer cells, thus exhibiting the potential to be utilized as natural antioxidants and candidates for cancer treatment (61). In our study, melanoma subjects not only exhibited a significant decrease in D. hansenii abundance compared with healthy subjects but a gradual drop in D. hansenii abundance was observed during anti–PD-1 therapy.

Our analysis has revealed interesting links between gut fungi and metastatic melanoma. Fungi generally regarded as beneficial were less abundant in melanoma samples, whereas fungi considered potentially harmful and presumably associated with cancer development were more abundant in melanoma. We also showed that specific fungal species may be associated with a cancer prognosis. Our results may allow for a better understanding of the complex interactions between the intestinal mycobiome and melanoma and may pave the way for new insights into the development and progression of this cancer. As a result, this could lead to the development of new effective yet safe therapeutic strategies (e.g., supplementation with probiotic Saccharomyces strains) and new biomarkers of cancer progression.

B. Pietrzak reports grants from the National Science Centre Poland during the conduct of the study. M. Schmidt reports personal fees from the National Centre for Research and Development (Poland) and grants and personal fees from National Science Centre (Poland) during the conduct of the study; in addition, M. Schmidt has a patent for method of isolating microbial nucleic acids from feces pending. Ł. Galus reports personal fees from BMS, MSD, Novartis, Pierre Fabre, and Roch during the conduct of the study. J. Mackiewicz reports travel grants, advisory board, and lectures for Bristol Myers Squibb and MSD. A. Philips reports grants from the Polish National Center of Research and Development (NCBiR), the European Regional Development Fund under the Smart Growth Operational Programme, Submeasure 4.1.2, Regional research agendas, and grants from Polish National Center of Research (NCN) during the conduct of the study. No disclosures were reported by the other authors.

N. Szóstak: Data curation, software, formal analysis, investigation, visualization, methodology, writing–original draft, writing–review and editing. L. Handschuh: Methodology, project administration. A. Samelak-Czajka: Methodology. K. Tomela: Data curation, methodology. B. Pietrzak: Data curation, methodology. M. Schmidt: Validation. Ł. Galus: Resources, data curation. J. Mackiewicz: Resources, data curation. A. Mackiewicz: Resources, validation. P. Kozlowski: Validation. A. Philips: Conceptualization, resources, supervision, funding acquisition, validation, investigation, project administration.

This work was carried out under contract no. POIR.04.01.02-00-0025/17-00 to A. Philips as “Polish Microbiome Map” cofinanced by the European Regional Development Fund under the Smart Growth Operational Programme, Submeasure 4.1.2, Regional research agendas. This work was also supported by a research grant from the Polish National Science Centre to A. Philips (2021/41/B/NZ9/03765). The authors would like to thank Ardigen SA, especially Błażej Szczerba, Łukasz Pruss, and Kaja Milanowska-Zabel for help with taxonomy profiling.