The tumor microbiome as a predictor of outcomes in patients with metastatic melanoma treated with immune checkpoint inhibitors

Emerging evidence supports the important role of the tumor microbiome in oncogenesis, cancer immune phenotype, cancer progression, and treatment outcomes in many malignancies. In this study, we investigated the metastatic melanoma tumor microbiome and potential roles in association with clinical outcomes, such as survival, in patients with metastatic disease treated with immune checkpoint inhibitors (ICIs). Baseline tumor samples were collected from 71 patients with metastatic melanoma before treatment with ICIs. Bulk RNA-seq was conducted on the formalin-fixed paraffin-embedded (FFPE) tumor samples. Durable clinical benefit (primary clinical endpoint) following ICIs was defined as overall survival ≥24 months and no change to the primary drug regimen (responders). We processed RNA-seq reads to carefully identify exogenous sequences using the {exotic} tool. The 71 patients with metastatic melanoma ranged in age from 24 to 83 years, 59% were male, and 55% survived >24 months following the initiation of ICI treatment. Exogenous taxa were identified in the tumor RNA-seq, including bacteria, fungi, and viruses. We found differences in gene expression and microbe abundances in immunotherapy responsive versus non-responsive tumors. Responders showed significant enrichment of several microbes including Fusobacterium nucleatum, and non-responders showed enrichment of fungi, as well as several bacteria. These microbes correlated with immune-related gene expression signatures. Finally, we found that models for predicting prolonged survival with immunotherapy using both microbe abundances and gene expression outperformed models using either dataset alone. Our findings warrant further investigation and potentially support therapeutic strategies to modify the tumor microbiome in order to improve treatment outcomes with ICIs.


121
Advances in immunotherapy, including immune checkpoint inhibitors (ICIs), have transformed 122 the standard of care for many types of cancer, including melanoma. While ICIs have improved 123 outcomes for melanoma patients, many patients suffer from primary or secondary tumor 124 resistance. For example, at 6.5 years, the overall survival rates with ipilimumab plus nivolumab, 125 nivolumab, and ipilimumab were 49%, 42%, and 23%, respectively, as reported in the pivotal 126 CheckMate 067 trial (1). Furthermore, mechanisms of resistance to immunotherapy remain 127 poorly understood, and many treatments are associated with immune-mediated toxicities. 128 Therefore, there is an urgent need to develop and improve biomarkers predictive of benefit from 129 ICI therapy. 130 131 Numerous biomarkers that predict the response of melanoma to ICIs are under investigation, 132 including those based on clinical characteristics, genomics, transcriptomics, and epigenomics. 133 For genomics data, these predictive biomarkers include tumor mutational burden (TMB) (2), 134 neoantigen load (3) Recently, high-throughput transcriptome, genome, or amplicon-based sequencing data 148 demonstrated an abundance and variety of microbes' nucleic acids inside tumors (8). In some 149 cases, hundreds of negative controls and paraffin-only blocks were sequenced to ensure a 150 thorough understanding of the background signal and reagent contamination. Further, the 151 presence of microbes has been validated using fluorescence in situ hybridization (FISH) and 152 immunohistochemistry (IHC) (19). The microbes showed cancer specificity (9,12,13), and blood-153 based measurements could predict early-stage disease. These findings suggest that microbes 154 observed in high-throughput sequencing data may also correlate with treatment outcomes. 155 Recent efforts to use these microbes as biomarkers showed that while generally less predictive 156 of prognosis than gene expression, when combined with gene expression they increase the microbes correlates with these signatures, suggesting an interaction with the immune system, 165 and how including tumor microbes in these models improves their predictive accuracy. 166 167 Figure 1. Overview of the research strategy 168 169 RNA-seq data from tumor specimens are processed to microbe abundances and human gene expression. Each is 170 associated with IO response, and then integrative analyses combine them into a model to predict outcomes. was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made

Signatures and pathways analyses 236
Gene signature scores were calculated using the IOSig and tmesig R packages. In brief, for 237 each published gene signature, we collected and harmonized gene names using the NCBI 238 Entrez gene number. To quantify the published gene expression score, we first transformed the 239 gene expressions across samples within a cohort into a Z-score. Next, we averaged the 240 standardized Z-score across the number of genes in the signature as previously described 241 (15,23,24). This score is used to compare responders and non-responders of immunotherapies 242 within individual cohorts based on the AUROC as previously described (23). We performed 243 clustering of gene signatures based on the correlation of AUROC across multiple cohorts. 244 Within a cohort of patients, we stratified the patients into "high" or "low" groups based on the 245 mean of the Z-score. A Mann-Whitney U test was performed in comparing the two groups to 246 determine the difference, and the false discovery rate (FDR) of <0.05 was deemed to be 247 significant. Prediction of response to treatment outcomes 262 To assess the predictive ability of the RNA-seq and microbe data for tumor response to ICIs, 263 random forest classifiers were created using the randomForest R package. Models were based 264 on 5 sets of input data: (1) microbe data, (2) 31-gene signature Z-score, (3) immune-activated 265 gene signature Z-score, (4) microbe and 31-gene signature Z-score combined, and (5) microbe 266 and immune-activated gene signature Z-score combined. Models were constructed with 500 267 trees and fivefold cross-validation. Additionally, 5 seeds were used for each model resulting in 268 25 trained models based on each set of input features. The AUROC curve was used to assess 269 the overall performance of the trained models. This metric assesses the model classification 270 accuracy, where 1 is a perfect classifier and 0.5 is a random classifier.

Gene expression analysis and its association with response to ICIs 285
Gene expression profiles for the 71 patients with metastatic melanoma treated with ICIs were 286 obtained from ORIEN. We performed DEG analysis and identified five 5 genes (CLEC12A, 287 GBP1P1, CD96, CCL4, IDO1) that were over-expressed in the responders as compared to the 288 non-responders with log2FC >1 and adjusted p-value <0.1 (Figure 2A). Interestingly, these 5 289 genes were involved in immune modulation and have been previously identified in other studies 290 as predictive biomarkers associated with responders to ICIs. For example, CCL4 has been 291 previously identified as a biomarker in the 12-chemokine signature (13,14), as well as other 292 gene signatures predictive of neoadjuvant ipilimumab response (27). IDO1 has been identified 293 as a key marker in the IFN-γ signature (12) and gene signature predictive of response to ICIs in 294 lung cancer (28). CD96 is a marker that estimates CD8+ T cell infiltration (29,30). CD96 and 295 TIGIT along with the co-stimulatory receptor CD226 form a pathway that affects the immune 296 response in an analogous way to the CD28/CTLA-4 pathway (31). CLEC12A (32,33) and 297 GBP1P1 (

309
Next, we asked what gene sets and pathways were enriched or depleted in responders to ICIs. 310 We performed GSEA using the MSigDB Hallmark gene sets on the RNA-seq and found that 311 several immune-related gene sets were significantly enriched in responders (Figure 2B) non-responders as shown in Figure 2B. The GSEA results identified in this cohort are similar to 317 previously published studies (23). 318 319 We next hypothesized that tumor-infiltrating immune cells could associate with responses to 320 ICIs. To test this hypothesis, we performed cell-type deconvolution of the bulk RNA-seq using 321 CIBERSORT. From CIBERSORT results, we observed that responders had significantly (p-322 value <0.05) higher abundances of CD8+ T-cells, activated CD4+ memory T-cells, activated NK 323 cells, and M1 macrophages relative to non-responders, who were shown to have a significantly 324 higher amount of resting mast cells (Suppl. Figure 1). Similarly, when we performed GSEA 325 using TIMEx gene sets, we observed that 13 CD4+, CD8+, and NK-related cell types were 326 enriched in responders (FDR < 0.1), whereas the stromal cell type was enriched in non-327 responders ( Figure 2C). This suggests that the tumor microenvironment of responders had an 328 "immune-inflamed" phenotype, whereas non-responders had either "immune-excluded" or 329 "immune-desert" TME phenotypes.

The melanoma tumor microbiome and its association with response to ICIs 367
Exogenous taxa were identified in the tumor RNA-seq, including bacteria, fungi, and viruses. A 368 total of 54 phyla were observed, with Firmicutes being the most abundant phylum, followed by 369 Uroviricota (Figure 3A). Within the tumors responsive to immunotherapy, we found a significant 370 enrichment of several microbes, including Fusobacterium nucleatum, Porphyromonas 371 asaccharolytica, Nocardia mangyaensis, and Mollivirus sibericum. Comparatively, the cohort of 372 non-responsive tumors was found to have significant intratumoral enrichment of fungi and the 373 bacteria Delftia lacustris, Enterobacter hormaechei, Pseudomonas fluorescens, and Moraxella 374 osloensis ( Figure 3B). We observed no significant differences between alpha diversity metrics 375 of responders and non-responders (Welch 2 sample t-tests p-value >0.4) ( Figure 3C). We 376 found that the random forest classifiers based on microbe diversity measures with 5 rounds of 377 5-fold cross-validation performed poorly relative to our other microbe-based classifers ( Figure  378 3D).  390 We next asked whether microbe abundance in the tumor could be associated with tumor 391 intrinsic pathways or the composition of cell types in the tumor immune microenvironment. We 392 focused on the 15 microbes identified to be differentially expressed in relation to immunotherapy 393 response in melanoma. For the 15 microbes, 7 and 8 were associated with responders and non-394 responders of immunotherapy, respectively ( Figure 4A). To investigate the intrinsic pathways 395 that correlated with the microbes, we performed ssGSEA on melanoma patients using MSigDB 396 . CC-BY 4.0 International license available under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made

414
To further dissect the association of microbe abundance and the composition of cell types in the 415 context of immunotherapy responses in melanoma, we performed cell type deconvolution using 416 the bulk RNA-seq with TIMEx. We found that the two microbes, Fusobacterium nucleatum and 417 Porphyromonas asaccharolytica, were highly correlated with the enrichment of tumor-infiltrated 418 immune cell types, including CD8+ T cells, which are known predictors of immunotherapy 419 response ( Figure 4C). In contrast, the lack of tumor-infiltrated immune cell types was correlated 420 with microbes associated with non-responders. In particular, we observed that malignant and 421 stromal cell types were enriched in association with the 2 tumor microbes noted in non-422 responders, Theileria annulata and Moraxella osloensis (Figure 4C). The tumor immune cell 423 composition corroborated our previous findings (23). 424 425 Next, we asked whether the microbe abundance was associated with any gene signatures 426 predictive of immunotherapy responses. To investigate this question, we utilized 31 previously 427 published gene signatures that have been indicated to be associated with immunotherapy 428 responses (23). We correlated microbe abundance with these signatures, and found that gene 429 signatures associated with inflammation or immune activation were highly associated with 430 microbes abundant in responders ( Figure 4D). On the other hand, gene signatures associated 431 with immune-suppressive or intrinsic signaling were highly associated with microbes abundant 432 in non-responders ( Figure 4D). These results suggest that microbe abundance could provide a 433 different dimension in understanding the tumor immune microenvironment in predicting 434 immunotherapy responsiveness in melanoma. 435 436 We further hypothesized that combining microbe abundance features with gene expression 437 signatures could improve response prediction of melanoma to immunotherapy. To test this 438 hypothesis, we developed an ensemble learning random forest classifier using microbe 439 abundance and gene signatures identified to be associated with immunotherapy responses in 440 melanoma. We first developed the random forest classifier based on microbe abundance with 441 15 input features (microbe) and performed 5 rounds of 5-fold cross-validation on the melanoma 442 cohort (Figure 5). The average AUROC for the microbe classifier was 0.651. We also 443 constructed a random forest classifier based on 31 gene signatures (GeneSig_Z_score) or the 444 16 immune-activated gene signatures (Imm_Act_Z_score), and the AUROC values for GeneSig 445 or Imm_Act classifiers were 0.72 and 0.744, respectively ( Figure 5). Notably, when we 446 combined the microbe abundance and gene signatures to develop the random forest classifier, 447 the ensemble learning random forest classifiers for gene signatures plus microbe 448

Prediction of response using tumor gene expression and microbe abundance
(GS_Z_microbe) and immune-activated gene signatures plus microbe (Imm_Act_Z_microbe) 449 achieved 0.772 and 0.805, respectively (Figure 5). This suggests that microbe abundance 450 features provide a distinct layer of information in predicting response to immunotherapy and, 451 when combined with gene expression signatures, can improve the prediction of response to 452 immunotherapy in melanoma. 453 . CC-BY 4.0 International license available under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made although it is an established pathogen that has been linked to colorectal cancer (54). In our 500 study, it is associated with the same immune expression pathways as Fusobacterium 501 nucleatum, suggesting it acts through a similar mechanism. The diversity of mechanisms and 502 taxa suggests that additional mechanisms are likely. Furthermore, recent studies have identified 503 bacteria-derived human leucocyte antigen (HLA)-bound peptides in melanoma presented by 504 tumor cells could elicit immune reactivity. This intratumoral bacteria peptide repertoire could be 505 further explored to understand the mechanism by which bacteria modulate the immune system 506 and responses to therapy (55). The demonstration of the utility of high-throughput sequencing to 507 explore these correlations warrants a broader search. 508 509 Efforts have been made to identify predictors of response and resistance to ICIs. As previously 510 discussed, expression signatures have been established as predictors of ICI response in 511 metastatic melanoma (9,12,14,15,23,56). One such study assessing the model combining IFNγ 512 and TMB found that it was predicitve of response but not resistance (56). Another such study 513 developed a multi-omic-based classifier that successfully predicted response, but was also 514 unable to predict resistance (20). We showed significantly enriched taxa in both response 515 groups. We also showed that microbes alone are predictive of response/resistance to 516 immunotherapy and, when combined with gene expression, enhance the model's predictive 517 ability. Further studies are warranted to combine tumor microbiome abundance with other 518 clinical and "omics" (e.g., genomics and pathomics) for developing an accurate classifier for 519 predicting immunotherapy responses in melanoma. Our findings also warrant further research to 520 evaluate whether these correlations are causally associated with outcomes and their effect on 521 the tumor immune microenvironment and immune cell infiltration. 522 523 . CC-BY 4.0 International license available under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (which this version posted May 25, 2023. ; In conclusion, we found that the tumor microbiome in patients with metastatic melanoma was 524 significantly different in those that responded (>24 months survival) to treatment with ICIs from 525 those who didn't respond. Furthermore, the microbial communities had the ability to predict 526 response when incorporated into machine learning models. The tumor microbiome further 527 enhanced models to predict response when combined with gene expression data. Future 528 research has the potential to support therapeutic strategies to modify the tumor microbiome to 529 improve ICI treatment outcomes. 530

531
The authors acknowledge the support and resources of the Ohio Supercomputer Center 532 (PAS1695). We would like to thank Angela Dahlberg, Editor, Division of Medical Oncology at 533 The Ohio State University Comprehensive Cancer Center, for editing and proofreading the 534 manuscript. 535