Multimodal analysis identifies microbiome changes linked to stem cell transplantation-associated diseases

Artacho, Alejandro; González-Torres, Cintya; Gómez-Cebrián, Nuria; Moles-Poveda, Paula; Pons, Javier; Jiménez, Nuria; Casanova, María Jinglei; Montoro, Juan; Balaguer, Aitana; Villalba, Marta; Chorão, Pedro; Puchades-Carrasco, Leonor; Sanz, Jaime; Ubeda, Carles

doi:10.1186/s40168-024-01948-0

Research
Open access
Published: 07 November 2024

Multimodal analysis identifies microbiome changes linked to stem cell transplantation-associated diseases

Alejandro Artacho¹^na1,
Cintya González-Torres¹^na1,
Nuria Gómez-Cebrián²,
Paula Moles-Poveda³,
Javier Pons¹,
Nuria Jiménez¹,
María Jinglei Casanova³,
Juan Montoro³,
Aitana Balaguer³,
Marta Villalba³,
Pedro Chorão³,
Leonor Puchades-Carrasco²,
Jaime Sanz^3,4,5 &
…
Carles Ubeda^1,6

Microbiome volume 12, Article number: 229 (2024) Cite this article

2081 Accesses
8 Altmetric
Metrics details

Abstract

Background

Allogeneic hematopoietic stem cell transplantation (allo-HSCT) is one of the most efficient therapeutic options available to cure many hematological malignancies. However, severe complications derived from this procedure, including graft-versus-host disease (GVHD) and infections, can limit its success and negatively impact survival. Previous studies have shown that alterations in the microbiome are associated with the development of allo-HSCT-derived complications. However, most studies relied on single techniques that can only analyze a unique aspect of the microbiome, which hinders our ability to understand how microbiome alterations drive allo-HSCT-associated diseases.

Results

Here, we have applied multiple “omic” techniques (16S rRNA and shotgun sequencing, targeted and un-targeted metabolomics) in combination with machine learning approaches to define the most significant microbiome changes following allo-HSCT at multiple modalities (bacterial taxa, encoded functions, and derived metabolites). In addition, multivariate approaches were applied to study interactions among the various microbiome modalities (the interactome). Our results show that the microbiome of transplanted patients exhibits substantial changes in all studied modalities. These include depletion of beneficial microbes, mainly from the Clostridiales order, loss of their bacterial encoded functions required for the synthesis of key metabolites, and a reduction in metabolic end products such as short chain fatty acids (SCFAs). These changes were followed by an expansion of bacteria that frequently cause infections after allo-HSCT, including several Staphylococcus species, which benefit from the reduction of bacteriostatic SCFAs. Additionally, we found specific alterations in all microbiome modalities that distinguished those patients who subsequently developed GVHD, including depletion of anti-inflammatory commensals, protective reactive oxygen detoxifying enzymes, and immunoregulatory metabolites such as acetate or malonate. Moreover, extensive shifts in the homeostatic relationship between bacteria and their metabolic products (e.g., Faecalibacterium and butyrate) were detected mainly in patients who later developed GVHD.

Conclusions

We have identified specific microbiome changes at different modalities (microbial taxa, their encoded genes, and synthetized metabolites) and at the interface between them (the interactome) that precede the development of complications associated with allo-HSCT. These identified microbial features provide novel targets for the design of microbiome-based strategies to prevent diseases associated with stem cell transplantation.

Video Abstract

Introduction

Allogeneic hematopoietic stem cell transplantation (allo-HSCT) is a frequently used therapeutic option for hematologic and non-hematologic disorders [1]. However, this procedure comes with risks for several complications, including graft-versus-host disease (GVHD) and infections [2,3,4], that limit its success. Infections can be facilitated during allo-HSCT by conditioning regimens, such as chemotherapy or radiotherapy, that induce neutropenia and mucosal injury, which favors bacterial translocation from the intestinal tract⁴. On the other hand, the inflammatory GVHD syndrome develops when transplanted donor T cells react against recipient host tissues, including the gut epithelium [5]. The initiation of these two complications is associated with changes occurring during allo-HSCT in the bacterial communities inhabiting the intestinal tract—the gut microbiome [2]. Indeed, patients undergoing allo-HSCT frequently receive antibiotics that alter their microbiomes. These alterations include the expansion of opportunistic pathogens such as Enterobacteriaceae or Enterococcus strains [3, 6,7,8], which promotes bloodstream infections [3, 6, 9]. In addition, gut domination by Enterococcus facilitates GVHD development [8]. In parallel with the expansion of potentially detrimental bacteria, antibiotics eliminate commensal microbes that could play a protective role against GVHD development and infections. Using 16S rRNA taxonomical identification approaches, it has been demonstrated that the use of broad-spectrum antibiotics decreases the levels of Clostridiales in patients undergoing allo-HSCT [10], a microbial change that is associated in patients with the subsequent development of GVHD and the gut colonization by opportunistic pathogens [7, 11], and that promotes GVHD and pathogen colonization in mice [7, 11, 12]. Accordingly, lower levels of Blautia, a Clostridiales genus, have also been associated with GVHD development and lethality in multiple studies [11, 13,14,15,16,17,18]. Besides 16S rRNA sequencing, a few additional studies have started to apply other sequencing techniques (shotgun sequencing), that allow the identification of bacterial encoded functions that could define potential mechanisms by which commensal microbes impact disease outcomes. This approach has enabled the identification of changes in hundreds of bacterial-encoded genes before GVHD development [13, 19], although their potential role in disease progression is largely unknown. In addition, metabolomic techniques have identified changes in metabolites produced by the microbiome that could be playing a role in the progression of diseases associated with allo-HSCT. A targeted approach focusing on immunomodulatory metabolites identified reductions in butyrate, propionate, isovaleric acid, indole-3-carboxaldehyde and desaminotyrosine to be associated with increased risk of GVHD and transplant-related mortality [19], while an untargeted metabolomic approach detected a link between changes in tryptophan-derived metabolites and GVHD development [20]. In summary, multiple studies have used different omic approaches to define changes in bacterial taxa, bacterial-encoded genes, or derived metabolites that are associated with allo-HSCT complications. However, changes in all these different microbiome modalities (i.e., bacteria, encoded genes, and metabolites) and their relations within the same transplanted patient are largely unknown and have only recently started to be characterized [19]. This multimodal approach is required to fully understand how shifts in specific microbes during allo-HSCT affect final metabolite outputs that could impact disease development. Information that is required to design effective microbiome-based approaches to prevent allo-HSCT-associated diseases. To reach this broad view of microbiome dynamics, we have analyzed fecal samples from patients undergoing allo-HSCT using multiple “omic” techniques (16S rRNA sequencing, shotgun sequencing, and targeted/untargeted metabolomics) and applied machine learning and multivariant analytical approaches. Our results demonstrate that the fecal microbiome undergoes extensive shifts during allo-HSCT in the different studied modalities (i.e., bacteria, encoded genes, and metabolites). These changes include the decrease in SCFAs which facilitate the expansion of opportunistic pathogens that frequently cause bacteremia in patients (i.e., Staphylococcus aureus and S. epidermidis). Moreover, we identified in all studied modalities and at the interface between them (i.e., the bacterial-metabolite interactome) microbiome changes that distinguish those patients that later developed GVHD.

Results

Taxonomic changes associated with allo-HSCT

Fecal samples from patients with hematological malignancies undergoing allo-HSCT were collected to determine microbiome changes associated with the transplant procedure. Fecal samples were collected from each patient at two different time points: before allo-HSCT (average ± SD, 4.49 ± 4.24 days before) and after allo-HSCT, around the time of engraftment (14.01 days ± 1.15 days after). We were able to collect pre-transplant samples from every patient (N = 105) and post-transplant samples from most patients (N = 70, 66,7%). All samples were collected before GVHD development. Patient characteristics are described in Supplementary Table S1, while therapeutic agents received are indicated in Supplementary Table S2. Most patients were male (N = 59, 56.19%) with an average age of 45.76 ± 13.23 years at the date of transplant. The transplant source was either peripheral blood (N = 79, 75.24%), umbilical cord (N = 25, 23.81%), or bone marrow (N = 1, 0.95%). The most frequent hematological malignancies were either acute myeloblastic leukemia (N = 47, 45.71%), acute lymphoblastic leukemia (N = 19, 18.1%), non-Hodgkin lymphoma (N = 11, 10.46%), and myeloproliferative syndrome (N = 10, 9.5%). All patients received ciprofloxacin as prophylaxis upon hospital admission, while intravenous piperacillin-tazobactam, meropenem, vancomycin, or amikacin were administered as antibiotic therapies in a subset of patients (47.6%, 45.7%, 32.4%, and 38.1% of patients, respectively, Supplementary Table S2). Other less frequently used antibiotics, such as colistin or linezolid, were administered in less than 9% of the patients (Supplementary Table S2). We first characterized the composition of the microbiota through 16S rRNA gene sequencing, DADA2 analysis for the identification of amplicon sequence variants (ASVs), and naïve Bayesian-based taxonomic assignment. 83,571 ± 34,913 high-quality sequences (see “Methods” section) were obtained per sample which led to the identification of 155.7 ± 114 ASVs per sample. The identified ASV represented for most samples the maximum number of ASVs that could be detected based on rarefaction curves, which reach a plateau in most cases, and based on comparison with the Chao Index—the potential maximum number of ASVs—which was almost identical to the real number of detected ASVs (Supplementary Fig. S1). Thus, the obtained sequencing depth was sufficient to accurately determine the fecal microbiota composition in our patient’s cohort. Subsequently, we studied the microbiota changes induced after allo-HSCT. We first analyzed overall changes in the microbiota composition at the ASV level. As expected, considering previous studies [10], fecal samples collected after allo-HSCT were separated from samples collected before allo-HSCT in a principal coordinate analysis (PcoA) (Fig. 1A), and their microbiota composition significantly differed (permutational multivariate analysis of variance test (PERMANOVA), p = 0.002). Besides beta-diversity, a significant decrease in alpha-diversity (i.e., ASV richness and Shannon diversity index) was detected after allo-HSCT (Fig. 1B). Next, we applied the Boruta machine learning algorithm to identify bacterial taxa that could differentiate samples pre- and post-transplant. The Boruta algorithm identified 2 phyla, 8 classes, 7 orders, 10 families, 25 genera, and 30 ASVs whose abundance could discriminate between pre- and post-transplant samples (Fig. 1C, Supplementary Table S3). The changes detected could classify with high precision the time-point in which the sample was collected: Area under the receiver operating characteristic curve–AUC = 0.78–0.88, depending on the taxonomic level utilized (Fig. 1D). A result that suggests that significant and reproducible changes occur during allo-HSCT in patients. Most of the changes detected involved the depletion of taxa from the Clostridiales order. Indeed, 14 out of the 16 most significant changes detected at the genus level belonged to the order Clostridiales (Fig. 1E, log2FC >|1|, ANCOM-BC and Benjamini-Hochberg False Discovery Rate (FDR) q < 0.05). These include multiple commensal bacteria (e.g., Blautia, Faecalibacterium, Oscillibacter) that have been shown to confer benefits to the host due to their capacity to produce SCFAs, mainly butyrate or acetate [21]. Consistently with the changes detected at the genus level, most of the ASVs depleted belonged to the Clostridiales order and this order was depleted in fecal samples collected after allo-HSCT (Supplementary Table S3). While most changes detected consisted of a reduction of commensal bacteria upon allo-HSCT, a few bacterial changes were detected in the opposite direction (Supplementary Table S3 and Fig. 1E), including a highly significant increase in the abundance of Staphylococcus genus after transplant (Fig. 1E, ANCOM-BC, p = 1.2⁻¹¹). It has been proposed that gut colonization by certain Staphylococci species, such as S. aureus, can promote the subsequent development of bacteremia [22]. Therefore, we decided to confirm the obtained result using a more precise technology (shotgun sequencing) that could facilitate the determination of the species level. A total of 92 samples from 46 patients, for which we had enough material available, and sequencing depth (2.26 Gb ± 2.69 Gb high-quality sequences), were included in this analysis (Supplementary Table S1). An expansion of three staphylococci species (S. aureus, S. epidermidis, and S. simulans) was detected after allo-HSCT (Supplementary Fig. S2). Next, we investigated if the microbiota changes detected were associated with the antibiotic therapy received (Supplementary Fig. S3). All patients received ciprofloxacin, as prophylaxis, at hospital admission, before the collection of the first fecal sample. Thus, the effect of this particular antibiotic on the microbiota changes could not be evaluated. However, a subset of patients received other antibiotics as therapeutic agents after the first fecal sample was collected, including, as described above, amikacin, vancomycin, meropenem, or piperacillin-tazobactam (Supplementary Fig. S3A). Therefore, we analyzed if the administration of these antibiotics could explain changes in the microbiota composition detected after allo-HSCT. A higher reduction of microbiota richness and diversity was detected in those patients who received therapeutic antibiotics (Supplementary Fig. S3B–C). Albeit this association was only significant for those antibiotic regimens that contain amikacin in addition to other antibiotics (Mann–Whitney test, p < 0.02, Supplementary Fig. S3B–C, Supplementary Table S4). In contrast, no significant associations with the change in microbiota richness or diversity were found for additional variables analyzed, besides antibiotics (i.e., parenteral nutrition, mucositis, myelosuppression, gender, transplant source, or type of transplant, Supplementary Table S4). Next, considering that the length of the antibiotic therapy varies among patients (Supplementary Table S2), we evaluated if a higher number of days of antibiotic administration could be associated with a higher reduction in microbiota richness or diversity. Neither the length of amikacin, vancomycin, or piperacillin-tazobactam administration was found to be associated with the reduction of microbiota richness or diversity (one-sided Spearman correlation test, p > 0.3, Supplementary Table S4). In contrast, a higher number of days of meropenem administration was associated with a higher reduction in microbiota richness (one-sided Spearman correlation test, r = − 0.406, p = 0.053, Supplementary Fig. S3D). Subsequently, we analyzed if the length of antibiotic administration could also be associated with changes detected in specific taxa. To diminish potential false positive associations, we focused our analysis on those genera that developed a higher significant change (log2FC >|1|, q < 0.05), as depicted in Fig. 1E. Consistent with the association, detected for meropenem administration and microbiota richness, the number of days that meropenem was administered was also associated with a higher reduction of the Clostridiales taxa Blautia and Unclassified Lachnospiraceae and a higher expansion of Staphylococcus. In addition, we found that the length of piperacillin-tazobactam administration was associated with the higher depletion of Clostridium XIVb and the higher expansion of Staphylococcus (Supplementary Fig. S3D). In summary, major changes were detected in the fecal microbiota composition of patients undergoing allo-HSCT, including a reduction in microbial diversity and richness, depletion of potentially beneficial microbes, mainly Clostridiales, and the expansion of potential opportunistic pathogenic bacteria such as S. aureus. Moreover, some of these microbiota changes could be associated with the administration of specific antibiotics.

Metabolomic changes associated with allo-HSCT

Having shown that multiple commensal anaerobic bacteria from the order Clostridiales, which are major SCFAs producers, are depleted during allo-HSCT, we decided to evaluate the impact of allo-HSCT on the levels of the 3 major SCFAs produced by commensal microbes (acetate, propionate, and butyrate) using targeted-metabolomics. A highly significant reduction was detected in all 3 analyzed metabolites (Fig. 2A). Next, untargeted nuclear magnetic resonance (NMR) metabolomics was applied to identify changes in the levels of other metabolites. Orthogonal partial least square discriminant analysis (OPLS-DA) showed a pronounced shift in the metabolome upon allo-HSCT (Fig. 2B, R²Y = 0.759, Q²Y = 0.496). Subsequently, we identified NMR peaks/buckets with a higher contribution to the detected metabolome change (i.e., buckets with variables importance in projection (VIP) > 1; p < 0.05, Mann–Whitney test, Supplementary Table S5). Twenty-four of the detected buckets could be confidently assigned to 9 metabolites. Note that in some cases more than one bucket was assigned to the same metabolite (Supplementary Table S5). Confirming the results obtained with targeted metabolomics, a depletion of the 3 major SCFAs were detected in samples collected after allo-HSCT (Fig. 2C). In addition, an even higher reduction of malonate, a three-carbon dicarboxylic acid, was detected in post-transplant samples (Fig. 2C). In contrast, we detected an increase of several amino acids (threonine, glycine, serine, glutamate, and isoleucine) and the monosaccharide glucose upon allo-HSCT. To understand which microbiota changes could be causing these metabolomic changes, we first identified, in the pre-transplant sample, the commensal bacteria that were most positively associated with SCFA levels and therefore potentially responsible for their production. For this purpose, a multivariate partial least squares canonical (PLS-C) approach in combination with univariate analysis was applied (see “Methods” section). Most of the genera that were highly positively associated with SCFA levels in pre-transplant samples corresponded to Clostridiales genera that were highly depleted after allo-HSCT, including Roseburia, Blautia, Faecalibacterium, Oscillibacter, Intestinimonas, Ruminococcus, Clostridium_IV, Anaerostipes, Pseudoflavonifractor, and Clostridium_XIVb (Fig. 2D). In addition, most of the genera with a highest negative association with several amino acids (i.e., serine, isoleucine, glutamate and glycine) in pre-transplant samples corresponded also to Clostridiales genera highly depleted after allo-HSCT. These included Faecalibacterium, Anaerotroncus, Oscillibacter, Intestinimonas, Ruminococcus, Clostridium_IV, Anaerostipes, and Pseudoflavonifractor (Fig. 2D). Thus, in this case, the reduction in the levels of certain Clostridiales bacteria upon allo-HSCT was associated with a higher availability of amino acids. In relation to the metabolomic observed changes, we and others have previously shown that certain metabolite changes associated with the disruption of the microbiome, including a decrease in SCFAs and an increase of specific amino acids [7, 23], can promote the colonization of opportunistic multidrug-resistant pathogens within the Enterobacteriaceae family (MDRE). Although we did not analyze the levels of MDRE in this study, an expansion of another potentially opportunistic pathogenic bacteria was detected in this cohort upon allo-HSCT (i.e., Staphylococcus spp.) (Fig. 1E, Supplementary Fig. S2). Therefore, we wondered if, like with MDRE, the metabolomic changes detected could drive the expansion of Staphylococcus in the gut of allo-HSCT patients. The levels of the three major SCFAs were significantly and negatively associated with the Staphylococcus abundance after allo-HSCT transplant (PLS-C correlation value (r_pls) < − 0.4; Spearman correlation test and FDR q < 0.01). Accordingly, post-transplant samples with low SCFAs levels (< 2 μg/g) contained significantly higher levels of Staphylococcus (Fig. 2E, Mann–Whitney, p < 0.05). This result suggests that SCFAs could have a potentially negative effect on Staphylococcus growth and their depletion after allo-HSCT may favor Staphylococcus expansion. In contrast, other metabolites that changed after allo-HSCT (isoleucine, glucose, threonine, glycine, serine, glutamate, malonate) were not associated with Staphylococcus levels (r_pls <|0.25|; Spearman and FDR q > 0.05). To verify the potential role of SCFAs in inhibiting Staphylocococcus growth, an S. epidermidis strain isolated from a fecal sample collected from a hospitalized acute leukemia patient was grown in tryptic soy broth (TSB) in the presence of concentrations of SCFAs detected pre- or post-allo-HSCT. S. epidermidis growth was significantly more impaired in media containing SCFAs concentrations that mimic pre-transplant samples (Fig. 2F), as compared to media containing SCFA concentrations that mimic post-transplant samples. A similar result was obtained with an S. aureus strain (Fig. 2F). Altogether these results suggest that depletion of major microbiome commensals, including several Clostridiales taxa, leads to major changes in the metabolome of patients receiving an allo-HSCT. These changes included the depletion of the most abundant SCFAs, which promote the growth of potential opportunistic pathogens within the Staphylococcus genus.

Changes in encoded bacterial gene orthologs during allo-HSCT

Considering that depletion of major commensal anaerobic bacteria was detected during allo-HSCT, which potentially led to changes in the metabolome, we decided to identify bacterial encoded functions depleted during the transplant that could explain the switch detected in the levels of metabolites. To do so, we analyzed the metagenomes obtained through shotgun sequencing. Bacterial encoded open reading frames (ORFs) were identified and a total of 10,632 functional orthologs (5127 ± 991.1 per sample) were assigned using the Kyoto Encylopedia of Genes and Genomes (KEGG) database to identify changes in metabolic pathways and bacterial functions associated with allo-HSCT. In concordance with the results at the taxonomic level, the richness and diversity of KEGG orthologs (KOs) were lower in post-transplant samples (Fig. 3A), although the level of change was less pronounced compared with the taxonomic data. Next, we identified KOs whose abundance significantly differs depending on the collection time-point. The abundance of 264 KOs was found to be significantly different between both groups of samples (Mann–Whitney and FDR, q < 0.05, Fig. 3B). In most cases, the change consisted of a decrease in abundance of a particular KO upon allo-HSCT (82.58% of detected changes), a result that is consistent with the decrease in the bacterial diversity detected after the transplant. This result indicates that a reduction of certain taxonomic groups is accompanied by a depletion of specific bacterial functions. Moreover, 34 KOs were confirmed by the Boruta algorithm to accurately discriminate between pre- and post-transplant samples (Fig. 3C, D, Supplementary Table S6, AUC = 0.88). Thus allo-HSCT was associated with robust and consistent changes in functional ortholog abundance that could accurately classify samples based on the collection time-point. To gain a better understanding of the type of bacterial functions that were lost after transplant, a hypergeometric test was applied in order to identify metabolic pathways that were enriched with KOs whose abundance significantly changed. Using this test, we could not find any pathway significantly enriched with KOs that increase after transplant. However, 5 pathways were found to be significantly enriched with KOs that decreased after transplant (Hypergeometric test and FDR, q < 0.05). The most significant one was “cell growth” in which 14 out of 61 KOs analyzed were depleted upon allo-HSCT (Fig. 3E, F). The assigned function for all these genes was related to the production of spores, and consistent with being genes of the same metabolic pathway, their abundance significantly and positively correlated with each other (Fig. 3F). Clostridiales are major producers of spores in the human gut [24]. Thus, a reduction of specific taxa was accompanied by a depletion of a specific pathway mainly encoded by these particular taxa. Indeed, the 14 identified KOs were mainly encoded by Clostridiales (Supplementary Fig. S4). The second most significantly depleted pathway was porphyrin metabolism. Within this pathway, 11 out of 71 analyzed genes were significantly depleted (Hypergeometric test & FDR, q < 0.05; Fig. 3G). All of them were related to the synthesis of vitamin B12. Moreover, the abundance of other genes within this pathway was also significantly reduced albeit the test did not pass the adjusted p value (q value) threshold (Mann–Whitney test, p < 0.05, q > 0.05; Fig. 3G). These results suggest that allo-HSCT could have a detrimental effect on the ability of the microbiota to synthetize vitamin B12. A result that should be corroborated using other metabolomic techniques since none of the statistically significant NMR buckets could be assigned to vitamin B12. Other depleted pathways were energy metabolism, methane metabolism, and peptidases in which 8, 10, and 18 KOs were found to be significantly depleted after transplant, respectively. Finally, since the 16S rRNA analysis identified a reduction of potential SCFA producers upon allo-HSCT, which was associated with a decrease in SCFA levels, we decided to study the effect of allo-HSCT on key genes involved in the production of SCFAs. To this end, we focus our analysis on the last enzymes involved in the production of the 3 major SCFAs (acetate, propionate, and butyrate). This analysis identified a significant reduction in the abundance of functional orthologs encoding the enzymes of the last step for butyrate and acetate production (K00132, K01034, K01035, K00128, K01067, Fig. 3H, I). This result indicates that a likely explanation for the reduction of the levels of SCFAs detected after allo-HSCT was the depletion of commensal bacteria encoding enzymes required for the synthesis of these SCFAs. To verify this hypothesis, we analyzed the bacteria encoding the detected KOs and found that Clostridiales taxa were the most abundant bacteria encoding the KOs involved in butyrate production (i.e., K00132, K01034, and K01035) (Supplementary Fig. S5). Additionally, Clostridiales were also the taxa with a higher contribution to the levels of K00128, which converts acetaldehyde into acetate, and the second most relevant contributor to K01067, which converts acetyl-CoA to acetate (Supplementary Fig. S5). Altogether these results indicate that, during allo-HSCT, commensal bacteria, mostly from the order Clostridiales, are depleted, which leads to a decrease in the abundance of genes that are required to produce key metabolites, including SCFAs. In parallel, certain opportunistic pathogens, such as Staphylococcus, take advantage of these changes and expand in the gastrointestinal tract, which may facilitate their dissemination to the bloodstream.

Taxonomic changes associated with the future development of GVHD

Our previous results defined changes in the microbiome at multiple levels (i.e., taxonomy, genes, and metabolites) occurring after allo-HSCT. Some of the detected changes could be involved in the development of GVHD in this type of patient, including the reduction of Clostridiales and SCFAs levels, as previously reported [10,11,12, 25]. Forty-seven out of the 105 analyzed patients (44.76%) developed moderate/severe acute GVHD (grades II–IV, Supplementary Table S1). Thus, we decided to study if the detected changes upon allo-HSCT were exclusive or at least more pronounced in those patients who were going to develop GVHD. Since the post-transplant samples were always collected before GVHD development, the identified changes could potentially be involved in the initiation of the disease or could be used as biomarkers to predict patients at risk of developing GVHD. We first focused on changes in alpha diversity using the 16S rRNA sequencing data. A similar decrease in richness and Shannon diversity was detected in both groups of patients (Fig. 4A). Thus, at least in our patient’s cohort, microbial diversity dynamics was not a key determinant for GVHD development. Next, we evaluated if changes in specific taxa could distinguish both groups of patients. The Boruta algorithm identified changes in several bacteria, at different taxonomic levels, that could differentiate patients according to their future GVHD status (Fig. 4B, Supplementary Table S7). Moreover, most Boruta-detected changes were also found to be significantly different by the univariate non-parametric Mann–Whitney test (p < 0.05, Fig. 4D, Supplementary Table S7). At the order level, only one feature (Clostridiales) was found to differ between GVHD and non-GVHD patients. A decrease in Clostridiales was detected only in those patients who were going to develop GVHD (median log2FC = 0.01 in non-GVHD and − 3.36 in GVHD patients). Consistent with this change, at the family and genus levels, we found Clostridiales taxa whose reduction was more pronounced in those patients that were going to develop GVHD. These taxa include the family Ruminococcaceae, genera Clostridium_IV and Flavonifractor. In contrast, the opposite change was detected for the genus Parabacteroides (increased after allo-HSCT in patients that were going to develop GVHD and reduction in patients that were not going to develop the disease) (Fig. 4D). At the ASV level, the Boruta algorithm did not identify any ASV whose change was different depending on the future GVHD status. Next, we evaluated to what extent these changes could classify patients based on their future GVHD status. GVHD patients form a differentiated cluster in a PCoA analysis performed with the taxonomic features identified by the Boruta algorithm (Fig. 4C, p = 0.006, PERMANOVA test). In contrast, this patient’s clustering was lost in a PCoA implemented with all taxonomic features (Supplementary Fig. S6, p = 0.248, PERMANOVA test). Consistent with the PCoA results, a random forest model combining the 16 s taxonomic features identified by Boruta could classify patients according to their future GVHD status (AUC = 0.8). This result indicates that although both groups of patients follow similar overall changes in their microbiome, specific taxonomic changes distinguish them. Subsequently, we analyzed the shotgun sequencing data, which potentially allows taxonomic classification up to the species level, in order to identify additional microbiota changes that could be associated with GVHD development. The Boruta algorithm identified 5 taxonomic features at the species level with the capacity to discriminate patients according to their future GVHD status (AUC = 0.78; Supplementary Table S7). Among these features, 4 of them were confirmed to be significantly different among the two groups of patients also by a univariate analysis test (Fig. 4E, Mann–Whitney, p < 0.05). Consistent with the results obtained at higher taxonomic levels, 3 identified taxonomic features belong to the order clostridiales (Anaeromassilibacillus sp., Butyricicoccus sp, and unclassified Lachnoanaerobaculum) and underwent a higher depletion after allo-HSCT in those patients that will end up developing GVHD. In contrast, an expansion of the species Enterococcus faecium was detected in those patients who subsequently developed GVHD. Having shown that taxonomic dynamics of specific taxa could discriminate between both groups of patients, we tested if the information from a single time-point (either pre-transplant or post-transplant) could be sufficient to predict GVHD development. Random forest models built with Boruta selected features using 16S or shotgun sequencing data from either pre- or post-transplant samples predicted with lower accuracy future GVHD status than microbiota dynamics information (Supplementary Table S8, Supplementary Fig. S7). Thus, microbiota dynamics was the most accurate method to classify patients according to their future GVHD status and both 16S and shotgun sequencing data identified a higher depletion of multiple Clostridiales taxa upon allo-HSCT in those patients that will end up developing GVHD.

Taxonomic changes associated with the grade of GVHD and the organs affected

The previous analysis detected taxonomic changes that were associated with GVHD development independent of the organ that was affected by GVHD. Considering that microbiota changes were analyzed in fecal samples, we reasoned that the detected changes would be associated with patients developing GVHD in the lower gastrointestinal tract (LGI-GVHD). Supporting this hypothesis, the previously detected microbiota changes, associated with overall GVDH development, were also identified when the microbiota dynamics of patients that will develop LGI-GVHD were compared to those patients that did not develop any sign of the disease in any organ (i.e., higher depletion of Clostridiales bacteria, including Ruminococcaceae, Clostridium IV, Flavonifractor, Anaeromassilibacillus sp., Butyricicoccus sp, unclassified Lachnoanaerobaculum and higher expansion of the genus Parabacteroides and species Enterococcus faecium, Supplementary Fig. 8A, B). To identify additional microbiota changes that could differentiate patients developing LGI-GVHD, we applied the machine learning Boruta algorithm. Using this approach, we identified a significant, albeit modest, higher expansion of the family Porphyromonadaceae (and its corresponding order Bacteroidales) in those patients that subsequently developed LGI-GVHD (Supplementary Fig. S8C). In addition, a significant depletion of an unclassified species from the Clostridiales genus Eisenbergiella was detected in those patients who subsequently developed LGI-GVHD when the shotgun sequencing data was analyzed (Mann–Whitney test, p = 0.007, Supplementary Fig. S8D). It is important to mention, for the interpretation of these results, that out of the 22 patients that developed LGI-GVHD included in the analysis, 19 did also develop GVHD in other organs (14 in the skin, 2 in the liver, and 3 in both the liver and skin, Supplementary Table S1). Thus, microbiota changes associated exclusively with LGI-GVHD development could not be studied in our patient’s cohort due to the low number of patients developing GVHD exclusively in the gut (N = 3). In contrast, we identified a larger number of patients (N = 13) that exclusively developed signs of the disease in the skin. Subsequently, we studied if microbiome changes identified in the gut could be associated with the development of GVHD in a distant localization such as the skin. Gut microbiota changes previously identified to be associated with overall GVHD were in general non-significant and less pronounced in patients that developed GVHD exclusively in the skin (mean difference between groups of 1.59 ± 0.79, see “Methods” section and Supplementary Fig. S9) as compared to changes detected in patients developing GVHD in the gut (mean difference between groups of 4.25 ± 2.5, Supplementary Fig. S8). To identify potential microbiome changes that could be associated with developing GVHD exclusively in the skin, we applied a machine learning approach as described above. This analysis identified a single bacterial family (Veillonellaceae), from the order Selenomonadales, whose dynamics significantly differ according to the future GVHD status in the skin (Mann–Whitney, p = 0.011). Specifically, a depletion of this bacterial family was detected in those patients who did not develop any signs of the disease as compared to patients who developed GVHD exclusively in the skin (Supplementary Fig. S9B). Accordingly, the same changes were detected at a higher taxonomical classification level of this family—the order Selenomonadales—(Mann–Whitney, p = 0.01, Supplementary Fig. S9B). These changes were specific to skin GVHD since no significant differences were detected in these taxa in those patients who developed LGI-GVHD (one-sided Mann–Whitney, p > 0.15).

Next, we studied if specific microbiota changes could be associated with the severity of GVHD symptoms. To this end, we performed a correlation analysis between the changes in microbiota taxa found to be associated with GVHD (Fig. 4) and the GVHD grade developed (I, II, III, and IV). This analysis revealed that Parabacteroides were positively associated with GVHD grade (one-sided Spearman correlation analysis, p = 0.028, r = 0.331). Those patients with a higher expansion of Parabacteroides after allo-HSCT subsequently developed more severe GVHD symptoms (Supplementary Fig. S10). In contrast, a negative correlation was found with Butyricicoccus sp. Patients with a higher depletion of this taxa developed a more severe GVHD (Supplementary Fig. S10., Spearman correlation analysis, p = 0.052, r = − 0.356). Subsequently, PLS regression in combination with a univariate correlation analysis (see methods), was applied in order to identify additional microbiota changes that could be associated with GVHD grade. This analysis identified a single genus (Escherichia) and a single ASV from this genus to be negatively associated with GVHD grade (Supplementary Fig. S10). In addition, several ASVs from the genus Bacteroides were found to be positively associated with GVHD grade (Supplementary Fig. 10; Spearman correlation and FDR, q < 0.05). Considering the different effects of different Bacteroides species on GVHD development [26], we used the shotgun sequencing data to investigate potential Bacteroides species associated with GVHD severity. We did identify several Bacteroides species in our patient’s cohort with high prevalence and abundance (B. fragilis, B. caccae, B. thetaiomicron, or B. stercoris, median > 0.1%). However, the change in their abundance of these and other less abundant Bacteroides species was not significantly positively associated with GVHD severity (one-sided Spearman, p > 0.05). In contrast, a significant negative association was identified between the change in E. coli abundance and GVHD severity (one-sided Spearman correlation, r = − 0.35, p = 0.046).

Metabolomic changes associated with the future development of GVHD

Having shown that taxonomic changes in the microbiome distinguish patients according to their future GVHD status, we decided to study if these taxonomic differences could lead to distinct metabolome dynamics that could impact GVHD development. Considering the relevance of SCFAs on immune system development and inflammation, we first performed a targeted metabolomics analysis of SCFA levels and analyzed changes between pre- and post-transplant samples. No differences in propionate dynamics were detected depending on the future GVHD status (Supplementary Fig. S11, p = 0.534). A higher decrease in butyrate levels after allo-HSCT, although not significant (p = 0.217), was detected in those patients who were going to develop GVHD (Supplementary Fig. S11). On the other hand, while acetate levels remain stable in most non-GvHD patients (log2FC median = − 0.4), a significant decrease in acetate levels was detected in patients who were going to develop the disease (log2FC median = − 2.36) (Fig. 5A, Man-Whitney test, p = 0.049). Consistently, a significant depletion of acetate (log2FC < − 1) was associated with a higher risk of developing GVHD (Fig. 5A, p = 0.002, OR = 9.5). Next, we made use of untargeted metabolomic NMR data to investigate other metabolomic changes that could differentiate both groups of patients. We first applied an effect projection by means of orthogonal partial least squares (OPLS-EP) analysis, that takes advantage of a matrix containing the effect of the intervention (e.g., pre- to post-transplant data) on each patient to characterize differences related to treatment effect. Although this analysis revealed a significant difference in the metabolic profile before and after transplant (CV-ANOVA, p-value = 6 × 10⁻⁵), no overall differences were observed based on future GVHD status (Supplementary Fig. S12). Also, OPLS-DA NMR models were not able to differentiate patients according to their GVHD status when pre- or post-transplant samples, instead of dynamics, were used to generate the models (Q²Y < 0). This result is consistent with the taxonomic analysis in which only subtle differences were detected between patients according to their future GVHD status (Fig. 4) but changes in the overall microbiome structure were not detected (Supplementary Fig. S6). For this reason, we decided to apply an alternative approach to identify differences in specific metabolites that could distinguish both groups of patients (machine learning). The Boruta algorithm identified 14 NMR buckets out of 792, whose change after transplant was able to discriminate both groups of patients (Supplementary Table S9). Nine out of 14 detected changes were confirmed by a univariate non-parametric test (p < 0.05, Man-Whitney, Fig. 5C). One of these 9 significant buckets was assigned to acetate, confirming both through a targeted and an untargeted approach that the change in this metabolite differentiates both groups of patients. The bucket with the highest fold change difference was malonate (Fig. 5C), a metabolite that may have anti-inflammatory properties (see discussion section) and that was found in our previous analysis to be highly depleted after allo-HSCT (Fig. 2C). The abundance of this metabolite significantly decreased in GVHD but not in non-GVHD patients (log2FC median = − 3.341 in GVDH vs 0.28 in non-GVHD). Malonate was positively associated with the levels of Clostridium_IV pre-transplant (Fig. 2), a Clostridiales tax on that was significantly more depleted in those patients that were going to develop GVHD (Fig. 4). Thus, the higher reduction of this taxa in GVHD patients could have contributed to the greater decrease in malonate detected in these patients. Finally, an additional bucket that could differentiate patients based on future GVHD development was assigned to a single metabolite (Histidine). This amino acid followed the opposite pattern (increase after allo-HSCT in GVHD patients and decrease in non-GVHD patients). As with taxonomic changes, the detected metabolomic changes followed a similar trend in those patients that developed GVHD in the lower gastrointestinal (Supplementary Fig. S13). No additional metabolomic changes were identified when the Boruta algorithm was applied to discriminate LI-GVHD patients, an analysis that was not pursued with patients exclusively developing GVHD in the skin due to the limited number of these types of patients with metabolomic data (N = 7).

Changes in encoded bacterial gene orthologs associated with the future development of GVHD

Having identified changes, both at the taxonomic and metabolite level, that discriminate between patients according to their GVHD status, we decided to study bacterial encoded functions (KOs) whose dynamics differ depending on these two groups of patients. Like taxonomic variables, a similar change in KO diversity and richness was detected after allo-HSCT independently on the future GVHD status (Fig. 6A). Nevertheless, we applied the machine learning approach in order to determine if the abundance of any specific KO differentially changed after allo-HSCT in patients according to the future development of GVHD. Only 2 KOs (K05919 and K18918) were identified by the Boruta algorithm to discriminate between GVHD and non-GVHD patients and confirmed by a univariate test to differ between the two groups of patients (Fig. 6B, Mann–Whitney, p < 0.05). K18918 encodes for a transcriptional regulator of a toxin/antitoxin system. Genes encoding this enzyme increased upon allo-HSCT in most patients that later developed GVHD in our cohort, while decreased in non-GVHD patients. This KO was mainly encoded by bacteria from the family Enterobacteriaceae (Supplementary Fig. S14). The other discriminating KO (K05919) is a superoxide reductase, an enzyme that is expressed by anaerobic bacteria to detoxify reactive oxygen species. Genes encoding this enzyme were significantly more reduced after allo-HSCT in patients who were going to develop GVHD (Fig. 6B). Moreover, a higher reduction of this KO was significantly associated with a higher grade of GVHD (one-sided Spearman correlation test, p = 0.046, r = − 0.352). This KO is mainly encoded by Clostridiales bacteria in our patient’s cohort, including bacteria from the family Ruminococcaceae (Supplementary Fig. S14), a bacterial family that was depleted to a higher extent after allo-HSCT in those patients that were going to develop GVHD (Fig. 4). Next, we verified if similar changes also occurred in those patients that developed GVHD in the lower gastrointestinal tract. A similar change in the abundance of the 2 identified KOs was also detected in those patients who developed LGI-GVHD (Supplementary Fig. S15). In addition, by applying the Boruta algorithm, we identified an additional KO that could also distinguish specifically those patients developing LGI-GVHD (K03366). In this particular case, an expansion of K03366, encoding a butanediol dehydrogenase, was detected in those patients who were going to develop GVHD in the lower gastrointestinal tract (Mann–Whitney test, p = 0.001, Supplementary Fig. S15B).

In the analysis described above (Fig. 3B, C), in which we studied changes at the gene ortholog level after allo-HSCT, independently on GVHD status, the Boruta algorithm identified 34 KOs that could differentiate pre- and post-transplant samples based on their differential abundance. In addition, a univariate test coupled with a false discovery rate identified hundreds of KOs whose abundance significantly changed after transplant. However, as described above, only the dynamics of 2 KOs was found to differ between GVHD and non-GVHD patients according to the Boruta algorithm, while a univariate test coupled with false discovery rate could not detect additional significant differences (i.e., Mann–Whitney and FDR, q > 0.05). Altogether, these results suggest that, although microbiome-encoded bacterial functions do change after transplant, most of these changes occur independently of the future GVHD status. To confirm this hypothesis and increase our understanding of the bacteria functional potential in each type of patient, we plotted the log2FC for every KO whose abundance significantly changed after transplant, either increasing (Fig. 6C) or decreasing (Fig. 6D). The dynamics of the abundance for each single KO analyzed followed the same trend in both groups of patients, independently on future GVHD development. However, the level of change detected was significantly more pronounced for those patients who were going to develop GVHD (p < 0.0001; Student’s t-test). This occurs both in the case of KOs whose abundance increased (average log2FC = 1.187 in non-GVHD vs 2.354 in the GVHD group, Fig. 6C) or those whose abundance decreased (average log2FC = − 1.34 in non-GVHD vs − 1.83 in the GVHD group, Fig. 6D). This last result indicates that a higher dysbiosis at the bacterial functional level precedes the development of GVHD.

Finally, we studied if, within those patients who developed any sign of GVHD, there could be specific bacterial functions that could be associated with the grade of GVHD (I-IV). To this end, we used an untargeted approach as described above (i.e., PLS regression in combination with a univariate analysis). This analysis identified 31 KOs whose change was significantly associated with GVHD grade (PLS-R r value >|0.4|, Spearman correlation, and FDR < 0.05, Supplementary Fig. S16). The majority of the detected associations (N = 29, 93.5%) were negative (i.e., a higher depletion of specific KOs was associated with a higher GVHD grade). Thus, in summary, allo-HSCT led mainly to a depletion of bacterial functions encoded by the microbiome (Fig. 3B), and a higher depletion of bacterial functions was associated with GVHD development (Fig. 6C) and its severity (Supplementary Fig. 16).

The metabolome-bacterial interactome distinguishes patients according to their future GVHD status

Our previous analysis identified specific microbiome features that changed after allo-HSCT and defined some microbial features whose dynamics differ according to the future GVHD status. However, in those analyses, we did not take into account the interactions between commensal bacteria and how they can affect final microbiome outputs (i.e., metabolites). These interactions constantly occur in the gut microbiome, which behaves as a community [27], rather than isolated individuals. Thus, we decided to investigate if interactions among bacteria and metabolites could differ according to the collected time-point (i.e., pre-/post-transplant) or the GVHD status. To this end, interaction networks between commensal bacteria and metabolites were built for four groups of samples: (i) pre-transplant non-GVHD, (ii) post-transplant non-GVHD, (iii) pre-transplant GVHD and (iv) post-transplant GVHD. An initial view of the generated networks (Fig. 7A) suggests that a higher number of interactions between bacterial and metabolites are found in pre-transplant samples, which seem to form more complex and robust networks, less dependent on specific interactions that support the topology of the network. For example, post-transplant non-GVHD samples are divided into two subnetworks that could fall apart by breaking two concrete interactions between bacteria and metabolites. To validate these observations, we calculated network metrics used to define the properties of a network. Moreover, to validate differences in network metrics among groups, the metrics were calculated on 50 different generated networks using a subsampling approach (see methods). In each network we can find two major components: nodes that in this case are bacterial genera and metabolites, and edges: the connections between each node that represent the interactions between bacteria and metabolites (either positive or negative, as represented by the color of the edge). We first measure the “degree centrality” of the network which quantifies the average number of edges connected to a given node. As shown in Fig. 7B, the degree of centrality was significantly higher in networks built with pre-transplant samples as compared to post-transplant samples, independent of their GVHD status. This result was expected since specific nodes (bacteria) were reduced in post-transplant samples, which may have diminished the number of possible interactions with the remaining bacteria. Similarly, we found that the “Strength” of the interactions, in other words, the level of correlations between two variables, was higher in networks built with pre-transplant samples of both GVHD and non-GVHD groups of patients (Fig. 7B). Thus, not only the number of interactions was reduced in post-transplant samples but also these interactions were weaker. In concordance with this result and the possibility that some of the potential interactions among bacteria/metabolites may have been lost post-transplant, the edge density (the number of interactions detected out of the total possible interactions that could have been detected) was higher in pre-transplant samples (Fig. 7B). Altogether these results suggest that microbiome interactions of pre-transplant samples are more complex and robust than microbiomes post-allo-HSCT.

Despite global metrics of microbial-metabolite interactions followed a similar change in GVHD and non-GVHD patients, we noticed that specific interactions between bacteria and metabolites were distinct according to the future GVHD status. As an example, we focused on a particular association that has been experimentally validated: Faecalibacterium and butyrate. Faecalibacterium is a major butyrate producer [21] and accordingly was positively correlated with butyrate levels in pre-transplant samples (Fig. 7C, Spearman test, p < 0.05). This positive association was also found in the network of pre-transplant samples, which in this case is defined by a positive r_pls value of the edge that connects the two nodes (Faecalibacterium and butyrate, Fig. 7A). To validate this positive association, we calculated and plotted the r_pls value in each of the 50 subsampling generated networks. As shown in Fig. 7D, this value is close on average to + 0.5, a high positive r_pls value, in pre-transplant samples of both GVHD and non-GVHD patients. In contrast, this interaction was lost in post-transplant samples from non-GVHD patients (Spearman test, p = 0.13, Fig. 7C; r_pls value close to 0, Fig. 7D). Most strikingly, in post-transplant samples from GVHD patients, this interaction completely changed. Instead of positive interaction, Faecalibacterium was significantly but negatively correlated with butyrate levels (Spearman test, p = 0.003, Fig. 7C; r_pls value < − 0.5, Fig. 7D). This prompted us to investigate the reason for this change in post-transplant GVHD samples. One possibility was that another butyrate producer, a potential competitor for Faecalibacterium, could have occupied the niche of Faecalibacterium in samples from GVHD patients after transplant. If this were true, butyrate, now produced by the competitor, could be indirectly but negatively correlated with Faecalibacterium levels. Accordingly, we found a potential competitor for Faecalibacterium in the GVHD network: Clostridium_XIVb. This bacterial genus is highly and negatively correlated with Faecalibacterium and highly positively correlated with butyrate in GVHD post-transplant samples (Supplementary Fig. S17). Thus, a complete novel interaction between a bacterium and a metabolite emerged after allo-HSCT in patients who were going to develop GVHD. We next evaluated the rest of the possible interactions between bacteria and butyrate but also between bacteria and the other two major SCFAs (acetate and propionate). We focused on stronger interactions (i.e., those with an r_pls value >|0.5| in at least one of the four networks). We first evaluated through a principal component analysis the overall nature of these interactions. The first principal component, which explains 51.2% of the variability among network interactions (the interactome), separated the GVHD post-transplant interactome from the rest (Fig. 7E). Additionally, the second principal component, explaining 28.7% of the variability, separated the non-GVHD post-transplant interactome from the rest. On the other hand, the pre-transplant interactomes from GVHD and non-GVHD patients clustered closely in the principal component analysis (PCA) analysis. These results suggest that upon allo-HSCT, interactions between bacteria and metabolites are shifted, and this change is more pronounced in those patients who later developed GVHD. To identify the specific differences within the interactomes, the average of r_pls values obtained from 50 network subsampling generated networks was used to build a heatmap, where positive r_pls values are shown in red, while negative r_pls values are shown in blue (Fig. 7F). In addition, a dendrogram was generated to group in different clusters the type of interactions detected between bacteria and metabolites (Fig. 7F). Consistent with the PCA analysis, the larger cluster containing the higher number of interactions (cluster II) differentiates the post-transplant GVHD samples from the rest. In this cluster, multiple negative interactions were detected in the post-transplant GVHD network between SCFAs and commensal bacteria frequently found in healthy individuals, mostly from the Clostridiales order. In contrast, many of these interactions were positive in the other networks, consistent with the role of some of these bacteria in SCFA production (e.g., Butyricicoccus vs butyric acid). Thus, the natural relationship of multiple bacteria with microbiome-derived metabolites was altered in the network built with post-transplant GVHD samples. Other clusters, which contained fewer number of interactions followed a distinct pattern that could either differentiate pre-transplant from post-transplant networks (cluster I and IV) or GVHD from non-GVHD samples (cluster III). Altogether, these results show that allo-HSCT is associated not only with major changes in the microbiome at multiple levels (taxa, genetic functional potential, derived metabolites) but also at the interface between them. Some of these changes occur independently of GVHD development and may promote the expansion of certain opportunistic pathogens. However, more pronounced changes at all microbiome levels and mainly at the interface between them are associated with the subsequent initiation of GVHD.

Discussion

Methodological advances have allowed the study of the microbiome through different angles, including the identification of commensal bacterial species that conform to these microbial communities, the bacterial functions encoded in their genomes, and the metabolites derived from their biological activity. In this study, we have combined these methodologies to broadly characterize the microbiome changes during allo-HSCT and their potential implication on diseases associated with this curative procedure. We have shown that the microbiomes of patients undergoing allo-HSCT suffer remarkable changes at the taxonomic, genetic, and metabolomic levels. At the taxonomic level, the most significant changes were mainly the depletion of commensal bacteria from the Clostridiales order that accompanied an overall depletion of microbiota richness and diversity. These changes could be a consequence, as previously shown [28,29,30,31], of the administration of multiple antibiotics and chemotherapeutic agents. In this regard, we found that antibiotic therapies containing amikacin were associated with a higher reduction of microbiota richness and diversity. This result was unexpected considering that amikacin treatment has a narrower in vitro activity against commensals than piperacillin-tazobactam or meropenem [32]. Nevertheless, it is important to mention for the interpretation of these results, that amikacin was always given in our patient’s cohort in combination with other antibiotics (Supplementary Fig. S3). Thus, the effect of amikacin administration could have been different in a scenario where other antibiotics have not been administered. In fact, in vitro studies indicate that aminoglycosides, including amikacin and gentamycin, can synergize with other antibiotics to kill specific bacterial pathogens [33, 34], albeit the potential synergistic effect against commensal bacteria has hardly been studied [35]. In addition, in our hospital unit, amikacin is mainly included in antibiotic regimens for treating severe sepsis, gastroenteritis, or enterocolitis, conditions that are associated with changes in the gut microbiome [36, 37]. Thus, another possibility is that the effect detected for amikacin was due to the patient’s condition, rather than the action of the antibiotic itself. Therefore, additional, more controlled studies, should be performed to evaluate the in vivo effect of amikacin on the gut microbiota composition, either alone or in combination with other antibiotics. In addition, we identified that longer treatments with meropenem and piperacillin-tazobactam were associated with a higher decrease in richness and Clostridiales taxa. This result is in agreement with previous studies that have demonstrated that antibiotics with activity against anaerobes, such as meropenem or piperacillin-tazobactam, are associated with a depletion of Clostridiales taxa [38, 39]. Besides, longer treatments of meropenem and piperacillin-tazobactam were associated with the increase of bacteria from the genus Staphylococcus that was detected after allo-HSCT. Moreover, shotgun sequencing analysis confirmed the increase after allo-HSCT of two potentially pathogenic Staphylococcus species (S. epidermidis and S. aureus). Staphylococci are one of the most frequent pathogens causing bloodstream infections in patients receiving HSCT [40]. Although nares colonization, especially for S. aureus, has been considered the origin of infections [41], S. aureus can also colonize the gut, and multiple reports have demonstrated that intestinal carriage is a source of infection [42,43,44]. Gut expansion of opportunistic pathogens such as Enterococcus or Enterobacteriaceae promotes their dissemination to the bloodstream and subsequent infection [3, 6, 9]. Thus, the increase that we have detected in the abundance of Staphylococci species upon allo-HSCT may increase the risk of bacteremia in these patients, a hypothesis that should be confirmed in future studies.

In addition to 16S rRNA sequencing, we also defined the metabolomic changes associated with allo-HSCT. Using both a targeted and untargeted metabolomic approach, we detected a depletion of the three major SCFAs (butyrate, propionate, and acetate). This reduction was also detected in a previous work in which the levels of SCFAs were analyzed using a targeted-metabolomic approach in patients receiving HSCT [19, 25]. SCFAs are major products derived from the bacterial fermentation of complex carbohydrates present in the diet. Butyrate is mainly produced by Clostridiales [45], thus the reduction of these order could explain the decrease in butyrate levels observed. Indeed, correlation analysis identified Clostridales taxa to be significantly associated with butyrate levels in our patient’s cohort, and genes encoding enzymes required for the production of butyrate, whose abundance was depleted upon allo-HSCT, were mainly encoded by Clostridiales taxa in our patient’s cohort. On the other hand, acetate is mainly produced by species from the phylum Bacteroidetes, although some Clostridiales genus (Blautia) are also acetate producers [45, 46]. Notably, one of the major changes detected upon allo-HSCT was the depletion of Blautia, a taxon that was found to be positively associated with acetate levels in pre-transplant samples, and whose depletion could have contributed to a lower capacity of acetate production by the microbiome in post-transplant samples. In contrast, we found that the levels of SCFAs were negatively associated with Staphylococcus abundance in post-transplant samples. We and others have shown that SCFAs can directly diminish the growth of opportunistic pathogens from the Enterobacteriaceae family [7, 47]. Thus, SCFAs derived from commensal bacteria could also be preventing Staphylococcus intestinal expansion. Indeed, we were able to demonstrate that the concentrations of SCFAs detected in fecal samples inhibit the growth of S. aureus and S. epidermidis. This inhibition was significantly lower when Staphylococcus was grown in the presence of the SCFA concentrations found in post-transplant samples. Thus, the depletion of SCFAs upon allo-HSCT may contribute to the gut expansion of Staphylococcus after the transplant. In line with our results, microbiome production of propionate has been shown to be important for ameliorating S. aureus skin infections in mice [48].

Besides the role of SCFAs in the defense against pathogens, one of these acids (butyrate) has been shown to play a protective role in the development of GVHD in mice [12]. Consistent with this protective role, lower levels of butyrate were detected in pediatric patients who were going to develop GVHD ²⁵, and we have detected a higher depletion of butyrate after allo-HSCT, albeit not significant, in those patients that will develop GVHD. Additionally, we detected a significant reduction of acetate in those patients who developed GVHD, while acetate levels remained stable in most patients who did not develop the disease, suggesting that acetate could have also a protective role against GVHD. Indeed, acetate reduction was also detected in patients with severe GVHD in a previous study [17], albeit this analysis was done at GVHD onset, thus it was unclear if the change was a consequence of the disease or could participate in its initiation [17]. Acetate induces IgA production through stimulation of free fatty acid receptors [49]. IgA coats intestinal bacteria which diminish intestinal inflammation [50]. Thus, theoretically, acetate could confer protection against GVHD through IgA induction. However, although a recent study showed that IgA plasma cell numbers are decreasing in patients developing GVHD [51], the potential direct role of IgA on GVHD development is unclear. Alternatively, acetate can be used to produce butyrate by specific microbes [52]. Thus, a reduction of acetate could indirectly affect GVHD progression by diminishing the potential capacity of commensal bacteria to produce butyrate. A higher reduction of another acid (malonate) was detected in patients who were going to develop GVHD. Malonate is a dicarboxylic acid that is present in the food and can be synthetized and catabolized by commensal microbes. Indeed, both diet and administration of probiotics can influence the levels of malonate in the gut [53, 54]. Malonate is a well-known inhibitor of the succinate dehydrogenase, an enzyme that is crucial for the activation of macrophages and its inhibition induces an anti-inflammatory state in macrophages [55]. Donor monocyte-derived macrophages activate allogeneic T cells and promote GVHD [56]. Thus, although the role of malonate on GVHD has not been studied, it is possible that its anti-inflammatory effect on macrophages could protect against GVHD and lower levels would be required for the initiation of the disease. Future studies should evaluate this possibility and also investigate the mechanism by which the microbiome can impact gut malonate levels. Besides metabolites, we also detected specific microbes, at different taxonomic levels, that were associated with GVHD development. Consistent with the positive role of Enterococcus, and more specifically Enterococcus faecium in enhancing GVHD⁸, we detected a higher expansion of Enterococcus faecium in those patients that subsequently developed GVHD. In addition, we detected a higher reduction of bacteria from the order Clostridiales in patients who were going to develop the disease. This change was detected at the order, family, genus, species, and gene ortholog levels. Consistent with our result, a reduction of Clostridiales taxa has been shown to be associated with GVHD in patients in previous studies [14, 38, 57]. Moreover, administration of Clostridiales to mice reduces GVHD severity in part by increasing butyrate levels which increases intestinal epithelial junction cell integrity [11, 12]. In addition, Clostridiales induce Treg cell differentiation mainly in the low gastrointestinal tract [58], which could also help to prevent GVHD. In this regard, depletion of Clostridiales taxa was mostly associated with patients who developed signs of GVHD in the lower gastrointestinal tract in our patient’s cohort. In contrast, changes in these taxa (i.e., depletion of Ruminococcaceae, Flavonifractor, and Clostridium IV) were minor in patients that exclusively developed GVHD in the skin. This result suggests that while depletion of Clostridiales taxa precedes LGI-GVHD development, its reduction may not be so critical for GVHD progression in other locations such as the skin. In contrast, we identified a change in a bacterial gut family (Veillonelaceae) that differentiates those patients who developed signs of GVHD exclusively in the skin but not in the gut. Although the number of patients that could be included in this last analysis was quite limited, if confirmed in future validation studies, it will be important to understand how a change in commensal bacteria inhabiting the gastrointestinal tract could have an effect on GVHD occurring in a different location. One possibility is that immune populations and/or cytokines induced by the gut microbiome could have an effect on GVHD, preferentially in the skin. In this regard, different types of T cells may be relevant for GVHD development in different locations. Th17 cells, producers of IL17 and IL22, are key for cutaneous GVHD [59,60,61], in part due to preferential migration to this organ [62]. BATF-dependent IL-7RhiGM-CSF + T cells, producers of granulocyte–macrophage colony-stimulating factor, are most relevant for GVHD in the gut [63]. Of note, Th17 cells can be induced by gut commensals [64], and still induce inflammatory disorders in distal sites [65]. Thus, certain gut commensals, by inducing Th17 cells, could preferentially be causing GVHD in the skin, a hypothesis that could be tested using germ-free mice colonized with Th17 inducer commensals.

In addition to studying microbiome changes associated with GVHD in different locations, we have also identified commensal bacteria that are associated with GVHD severity. We found that within those patients who developed GVHD, an expansion of ASVs from the genus Bacteroides was associated with the subsequent progression of a more severe GVHD. This positive association was only found at the ASV but not at the genus level, suggesting that only certain Bacteroides are enhancing GVHD. Consistent with this result, a recent study demonstrated that different species of Bacteroides can have opposite effects on GVHD [26]. B. thetaiotaomicron degrades host mucin glycans which reduces the mucus layer, facilitating bacterial translocation and GVHD. In contrast, B. ovatus diminish GVHD by liberating the monosaccharide xylose, through polysaccharide degradation. This sugar is preferentially used by B. thetaiotaomicron as a carbon source, preventing the degradation of the mucus. Thus, depending on the glycan preferences, a Bacteroides could be beneficial or detrimental for GVHD.

At the genetic functional level, we identified a KEGG ortholog (KO5919) that could discriminate patients according to their future GVHD status. This KO is encoded by Clostridiales taxa, including the family Ruminococcaceae, and consistent with the dynamics of the taxa encoding this KO, a higher depletion of this KO after allo-HSCT was detected in patients that end up developing GVHD. Changes in this KO may just reflect the changes detected in bacteria found to be relevant for GVHD progression. However, another possibility is that the function of this KO plays a protective role in GVHD. This KO encodes for a superoxidase reductase, an enzyme that protects commensals by converting superoxide to hydrogen peroxide [66]. Reactive oxygen species, including superoxide, are produced by neutrophils recruited to the gastrointestinal tract after allo-HSCT and are critical for enhancing GVHD through tissue damage [67]. Therefore, a higher reduction of genes encoding the superoxide reductase in the gut could contribute to higher superoxide levels which could affect both the survival of potentially beneficial anaerobic commensal bacteria and enhance GVHD through tissue damage. In addition, we identified multiple KOs whose change in abundance was significantly associated with GVHD severity. The most significant KO (K10000) encodes for an arginine transport system in bacteria, whose depletion was associated with a higher grade of GVHD. Arginine is required for allogeneic T-cell responses and depletion of arginine results in significant GVHD reduction [68]. Therefore, considering the impact that the microbiome has on serum amino acid levels [69], a reduction in the capacity of the microbiome to uptake arginine may increase the availability of arginine and promote GVHD.

Finally, since we obtained multiple “omic” data from the same patients’ samples, we were able to evaluate the relationship between commensal bacteria and metabolic-derived products. Moreover, we also evaluated the interactions among bacteria, since the gut microbiome is a community of microbes that interact with each other, an important point that should be considered when evaluating the impact of the microbiome on disease progression and that has not been evaluated in the context of GVHD. Based on network metrics, the interactome of bacteria and metabolites was more complex and the strength of the interactions was higher in pre-transplant samples. A result that was expected considering that specific commensals were depleted upon allo-HSCT. Albeit higher complexity and stronger interactions are not necessarily linked to the stability of bacterial communities [70], we notice that the interactomes of post-transplant samples were more likely to be disrupted since subnetworks within the complete network were connected by very few interactions. Thus, disaggregation of these key interactions could lead to major changes in the topology of the network and the relationship between bacteria and metabolites. In this sense, extensive differences in the associations between bacteria and metabolic products were detected according to the time-point of collection and the future GVHD status. As an example, we found that Faecalibacterium, a key butyrate producer, was negatively correlated with butyrate levels in the gut of patients who later developed GVHD. Since it is unlikely that Faecalibacterium reduces the levels of butyrate in GVHD patients, we hypothesized that this negative association was indirect and was dependent on the presence of other bacteria in GVHD patients. In fact, our network analysis identified Clostridium_XIVb to be negatively correlated with Faecalibacterium and positively correlated with butyrate levels. Thus, it is possible that Clostridium_XIVb occupied the niche of Faecalibacterium and is now a key butyrate producer bacterium in the gut of GVHD patients. Although this type of analysis should be confirmed with experimental data, it suggests that a different ecosystem has emerged upon allo-HSCT, being this new ecosystem more dissimilar in those patients that will develop GVHD. Understanding the functionality and resilience of this new ecosystem will be important for applying microbiome-based therapies to prevent disease development. For example, it is likely, considering the metagenomic and metabolomic data, that post-transplant bacteria, occupying the niche and metabolic activities of pre-allo-HSCT communities, are less efficient in producing beneficial molecules such as SCFAs. If this is the case, it will be desirable to restore the original bacterial populations inhabiting the gut. However, the new ecosystem could be resistant to the introduction of commensal microbes. Thus, understanding the interactions between bacteria in the newly generated microbiome communities after allo-HSCT will be important to be able to restore a healthy microbiome, less capable of promoting GVHD.

We acknowledge several limitations in our study. First, our study lacks a validation cohort and the number of participants was limited. Therefore, microbiome changes detected in this study associated with allo-HSCT or GVHD development should be confirmed in subsequent studies. Nevertheless, some of the microbiome changes detected (e.g., Clostridiales and SCFAs depletion upon allo-HSCT) have also been detected in previous studies [10, 25], suggesting that our methods and analytical approaches were capable of detecting meaningful microbiome changes in the context of allo-HSCT. Second, associations identified between microbiome changes and allo-HSCT-associated diseases should be validated with experimental models. In this sense, we were able to validate at least one of the detected associations (i.e., the inhibitory effect of SCFAs against Staphylococcus spp). Despite these limitations, by combining multiple “omic” approaches, including 16S rRNA gene sequencing, shotgun sequencing, and NMR-based metabolomics, our study provides a novel and broad view of the microbiome changes in the context of allo-HSCT, which could help the design of novel microbiome-based approaches to prevent associated diseases, such as infections or GVHD.

Methods

Study design and sample collection

Fecal samples were collected from November 2013 to September 2017 from all patients undergoing allo-HSCT at the Hospital La Fe (Valencia, Spain) who agreed to participate in this study. One fecal sample was obtained from each patient before allo-HSCT, while another sample was collected approximately 2 weeks after allo-HSCT. After collection, specimens were initially stored at 4 ºC. Within 24 h, specimens were aliquoted into cryovials and stored at − 80 ºC. The study was approved by the institutional review board at the Hospital La Fe (2015/0369) and the Ethics Committee of CEIC Dirección General de Salud Pública y Centro Superior de Investigación en Salud Pública (20130515/08). Patient-related metadata was prospectively collected and recorded in a computerized database in Access®. The metadata collected included: (i) type of hematological malignancy (see Supplementary Table S1), (ii) the type of transplant received, (ii) transplant source (iii) conditioning intensity, either reduced or complete [71]; (iv) antibiotics received; (v) parenteral nutrition; (vi) development of mucositis; (vii) gender and (viiii) age. Acute GVHD was defined and graded according to standard criteria [72]. Patients were divided into two groups for subsequent analysis: patients with moderate to severe GVHD (score II–IV) and patients with mild (score I) or no signs of GVHD. Additional analysis (shown in Supplementary Figs. S8–9, 13, 15) was performed in which the microbiota of patients that develop GVHD either in the lower gastrointestinal tract or specifically in the skin (Supplementary Table S1) was compared with that of one of the patients that did not develop any sign of the disease in any organ.

Fecal DNA extraction

Bacterial DNA from fecal samples was extracted using the QIAamp® DNA Fast Stool Mini kit (QIAGEN) as previously described [73]. Extractions were performed according to manufacturer instructions with the introduction of a previous mechanic disruption step with bead-beating. Briefly, samples (~ 200 mg) were resuspended in the first buffer of the extraction kit. Subsequently, the samples were shaken on a Vortex-Genie 2 equipped with Vortex Adapters (Mobio) at maximum speed for 5 min in the presence of 500 µl of glass micro-beads (acid-washed glass beads 150-–212 µm, Sigma®). Then, we followed the protocol of QIAamp® DNA Fast Stool Mini kit.

16S rRNA sequencing and analysis

The V3-V4 regions of the 16S rRNA gene were amplified (Kapa HiFi HotStart Ready Mix), indexed with Nextera® XT Index Kit (96 indexes, 384 samples), and sequenced as described in the manual for “16S Metagenomic Sequencing Library Preparation” of the MiSeq platform (Illumina) using MiSeq Reagent Kit V3. 16S rRNA gene sequences were denoised and processed with DADA2 v1.8 [74], as we have previously described [75]. Briefly, prior to denoising, as recommended in the DADA2 pipeline tutorial, the 3′ ends of the forward and reverse pair-end sequences (35 and 75 nucleotides, respectively), which contained lower-quality nucleotides, were removed from each sequence. The first nucleotides (17 in the forward read and 21 in the reverse read) were also trimmed to remove the V3–V4 primers. Sequences containing undetermined nucleotides (N) were discarded. After denoising, the forward and reverse pair-end sequences were merged together, with a minimum overlapping region of 15 bases and a maximum of 1 mismatch. DADA2 and the command removeBimeraDenovo were subsequently applied to remove chimeras. Taxonomy was assigned for each amplicon sequence variant (ASV) using the Silva 138 database [76] as a reference and the naive Bayesian classifier implemented in DADA2. Shannon index was obtained at the ASV level using the package vegan v2.5–3 and the R 3.6.0 software.

Metagenomic shotgun sequencing

Metagenomic shotgun sequencing was performed using a NextSeq platform (Illumina, USA) according to the manufacturer’s recommendations for obtaining paired-end sequences of 150 base pairs. Obtained sequences were filtered by quality with Fastp v.0.20.1 [77] using default parameters to eliminate low-quality bases that were trimmed from the tail of each read, low complex reads, and reads shorter than 35 bases. Next, all the reads passing the previous filters were mapped onto the Homo sapiens genome (GCA_000001405.28 GRCh38.p13 assembly) by means of Bowtie2 (version 2.4.2) [78]. Finally, the reads which did not align against the Homo sapiens genome were analyzed with the SqueezeMeta pipeline with coassembly mode (version 1.3.1; database built on September 2020) [79]. Only samples containing > 10⁸ bp after the filtering process were included in the SqueezeMeta analysis. Tabular outputs were generated from the SqueezeMeta results using the sqm2tables.py SqueezeMeta script. The generated tables contained the abundance of each ORFs identified, the functional annotation of each ORF using the KEGG database [80], the taxonomic annotation of each ORF using the GenBank nr database, and taxonomic classifications at the species level. The TPM approach was used for normalization of metagenomic count data.

Metabolomic analysis

Fecal samples were resuspended in phosphate saline buffer (PBS). Afterward, bacterial cells were centrifuged at 13,200 rpm and the supernatant was analyzed through NMR-spectrometry as we have previously described [81]. ¹H-NMR spectra were acquired using a Bruker AVANCE-TM 600 MHz Spectrometer with a 5-mm BBI probe. ¹H -NMR experiments were acquired at 310 K (37 ºC) for every sample. A one-dimensional (1D) Carr − Purcell − Meiboom − Gill (CPMG) spin-echo pulse sequence was collected for each sample with a total of 256 scans and 65 K data points over a spectral width of 20 ppm and a relaxation delay of 4 s between free induction decays (FIDs) [82, 83]. A water presaturation pulse of 25 Hz was applied throughout the relaxation delays to improve solvent suppression [82, 83]. Then, fecal spectra were multiplied by a line-broadening factor of 1 Hz and Fourier-transformed. Finally, all spectra were automatically phased, baseline corrected, and referenced to the methyl group signal of alanine (δ 1.47 ppm) using TopSpin 3.5 (Bruker Biospin).

After NMR acquisition, NMRProcFlow v.1.2.28 was used. NMRProcFlow is an open-source software for data processing prior to multivariate statistical analysis, including, among other tools, solvent signal suppression, internal calibration, phase, baseline, and misalignment corrections, bucketing, and normalization. Briefly, spectra were binned into 0.01 ppm wide rectangular buckets over the spectral region δ 9.10 − 0.60 ppm, and the residual water signal region (δ 5.11 − 4.47 ppm) was excluded from further analyses to avoid interferences arising from differences in water suppression. Signal intensities obtained for each spectrum were imported into SIMCA-P 14.0 software (Umetrics AB). Metabolites were assigned to peaks that could discriminate groups of samples (see below the statistical methods for their identification) using the Bruker NMR metabolic profiling database BBIOREFCODE 2.0.0 database (Bruker Biospin), in combination with publicly available databases and the Chenomx NMR Suite v.8.6 [84, 85].

For the quantification of SCFA, calibration curves were obtained using standard compounds, butyric acid: 0–6 mM, propionic acid: 0–5 mM, and acetic acid: 0–10 mM. MestreNova software (ver. 8.0) was used to apply a deconvolution algorithm for the integration of overlapping signals and to integrate the specific spectral region of each SCFA, butyric acid: 0.88 ppm, propionic acid: 1.05 ppm and acetic acid: 1.90 ppm. Finally, SCFA concentrations were calculated based on peak intensities, and normalized to the amount of fecal sample (μmol/g). For a few samples (Supplementary Table S10), there was an overlap of the butyrate and propionate signal with other NMR signals that impeded the correct identification of these metabolites. These samples were not included in the analysis shown in Fig. 2 and Supplementary Fig. S11.

In vitro culture experiments

The S. epidermidis strain CU867, isolated from the feces of a leukemia patient, or the S. aureus strain RN4220 [86], were grown on tryptic soy agar overnight. Grown bacteria were resuspended on PBS and 10⁴ colony-forming units (CFUs) were inoculated into fresh tryptic soy broth (TSB) containing acetate, butyrate, and propionate. Two different concentrations of SCFAs were used according to the average of the levels detected in patients before and after allo-HSCT. As a control, the Staphylococci strains were inoculated into TSB that did not contain SCFAs. The pH of the media was adjusted to 6,5 before bacterial inoculation (a pH within the range of 6–7 found in the human large intestine [87]). Two hundred microliters of TSB were used per replicate which were located in a well from a 96-well plate. All procedures were performed on an anaerobic chamber and all media was pre-reduced before bacterial inoculation. To maintain the anaerobic conditions, 30 µl of autoclaved mineral oil was added to the top of the 96-well plate, as we have previously described [73]. Subsequently, the absorbance was measured outside the anaerobic chamber in a spectrophotometer Tecan Infinite F200 at 570 nm of wavelength, using Magellan program v.6.6. The plate was incubated at 37 °C. The absorbance was measured every 30 min with the previous shaking. The area under the growth curve was calculated using Prism 9.5.0.

Statistical analysis

To compare the abundance of microbial taxa between pre- and post-allo-HSCT samples, the analysis of compositions of microbiomes with bias correction (ANCOM-BC) was applied [88]. Only prevalent taxa (present in at least 40% of samples within at least one of the groups under comparison) were included in this analysis. Information regarding the patient donor of the sample was included as a fixed effect via the “fix_formula” parameter of the “ancombc” function. To adjust for multiple hypothesis testing and obtaining q values, we used the FDR approach by Benjamini and Hochberg [89], implemented in the p.adjust function of the R stats package v.3.6.0. The non-parametric Mann–Whitney test and the FDR approach was applied for the identification of significant differences in the abundance of KOs among groups of samples. A hypergeometric test and FDR approach were applied in order to identify metabolic pathways that were enriched with KOs whose abundance significantly increased or decreased after transplant (q < 0.05). A pathway was considered significant if more than 5 genes within the pathway were found to increase or more than 5 genes were found to decrease after the transplant and the hypergeometric q value was lower than 0.05.

To compare differences in dynamics of microbiome features and metabolites between patients who did not develop GVHD and those who will end up developing GVHD, the log2FC between post- and pre-transplant samples was calculated by dividing the relative abundance of a particular feature in the post-transplant sample by the relative abundance of that taxa in the pre-transplant sample. Subsequently, the non-parametric Mann–Whitney test was applied to identify significant differential changes among groups of samples. Only prevalent features (present in at least 40% of each group under comparison) were included in this analysis. For shotgun-identified species, besides the “prevalence” threshold, due to the low abundance of specific bacteria that pass this threshold, a minimum abundant average threshold was also applied (> 0.001%).

To identify associations between antibiotic administration and microbiome changes, we classified each patient according to the antibiotics received. Patients in which the antibiotic therapy was initiated before the pre-transplant sample were excluded from this analysis (N = 6). Subsequently, using Dunn’s test, the log2FC of microbiota diversity or richness was compared between patients receiving an antibiotic therapy containing the antibiotic of interest vs those patients who only received the prophylactic antibiotic (ciprofloxacin). In addition, within those patients who received an antibiotic therapy besides the prophylaxis, we performed a Spearman correlation analysis between the number of days that a particular antibiotic was administered and the log2FC of the abundance of specific microbiome features described in the results section.

The Boruta algorithm was applied [90], as we have previously described [81], to identify microbiome features that could discriminate between groups of samples with a minimum of 10 samples. Subsequently, to calculate the power of classifying groups of samples using the Boruta confirmed features, a random forests model was fitted using these features and all the samples except one. This model was fitted with at least 3 features confirmed by the Boruta algorithm. Next, the likelihood of belonging to each group of samples was calculated in the remaining sample using the function predict.randomForest from the R package randomForest. This process was repeated for all the available samples. Subsequently, the area under the ROC curve was calculated with the prediction values obtained for each patient using pracma R package 2.2.2. Boruta selected features that could differentiate the abundance or the dynamics of abundance between groups of samples were further confirmed using a univariate test, as specified above.

To analyze global differences in the microbiota at the taxonomic level between pre- and post-transplant samples, the Bray–Curtis distance between pairs of samples was calculated using the vegan package. Bray–Curtis distance was calculated using the ASV count sequence table generated with DADA2 but normalized using ANCOM-BC [88]. Subsequently, PCoA analysis was performed using the Bray–Curtis distances between a pair of samples that were calculated using the package vegan v2.5–3 and the R 3.4.0 software. In order to analyze community-level differences in the microbiome among groups of samples, a non-parametric test, permutational multivariate analysis of variance (PERMANOVA), was applied using the Adonis function from the R vegan package.

Prior to multivariate metabolomic analyses, NMR data were Pareto scaled by dividing each variable by the square root of 1/SD, where SD represents the standard deviation value of each variable [91]. Orthogonal partial least squares (OPLS) discriminant analyses were conducted to determine if groups of samples could be separated based on the NMR data. OPLS is a supervised approach used to produce models of “best fit” by minimizing the possible contribution of intergroup variability. The default method of sevenfold internal cross-validation of the software was applied, from which Q²Y (predictive ability parameter, estimated by cross-validation) and R²Y (goodness of fit parameter) values were extracted. Those parameters were used for the evaluation of the quality of the OPLS models obtained. To identify the variations responsible for the separation between groups of samples, the VIP list of each OPLS model was carefully inspected and used to identify which variables were relevant for the discrimination between groups (VIP > 1). Mann–Whitney test was applied to confirm the differences (p < 0.05) in the abundance of identified variables among groups of samples.

To identify overall metabolome differences associated with future GVHD development, we performed multivariate modeling using effect projections by means of orthogonal partial least squares (OPLS-EP), an analytical approach developed for dependent or matched samples that allows the multivariate analysis of paired samples over time [92]. An effect matrix was calculated with Excel 2019 in which the pre-transplant values were subtracted from the matching post-transplant values. The effect matrix was modeled in relation to a response vector (Y = 1). Prior to modeling in Simca, the 792 variables from 48 patients were Pareto scaled (ParN) but not cantered, and cross-validation groups were set to 7 (default). The model estimates Yhat (Yhat = 1, model target) were used to visualize the magnitude of the effect for each matched pair. The significance of the OPLS-EP model was calculated using CV-ANOVA.

To identify associations between microbiome features (i.e., metabolites and bacterial genera), partial least squares canonical analysis (PLS-C) was applied using MixOmics 6.10.9 R package [93]. The associations of two features in PLS-C are defined as the inner product of their corresponding vectors in the principal component plane inferred by PLS-C. The inner product values range from + 1 as the maximum positive association to − 1 as the maximum negative association. These values are indicated in the text as r_pls values to differentiate them from the correlation values (r) derived from univariate analysis. Associations identified with PLS were further confirmed using a univariate test (i.e., Spearman Correlation test) coupled with Benjamini and Hochberg correction to adjust for multiple hypothesis testing.

To identify associations between microbiome features and GVHD grade (I–IV), partial least squares regression was applied. Associations were confirmed using a univariate test (i.e., Spearman correlation test) coupled with Benjamini and Hochberg correction to adjust for multiple hypothesis testing.

To build networks based on PLS-C-identified associations, the mixOmics 6.10.9 R package was used. Only associations with r_pls values >|0.5| were used to build the networks. Network metrics were calculated using the igraph software package for complex network research [94]. To validate differences among groups on network metrics, 50 independent PLS-C analyses were performed on random subsets of data containing 2/3 of the sample from each group in each subset. With the obtained data, 50 independent networks were built, and the metrics of each network were calculated and plotted. In addition, the average of the r_pls values was used to build a heatmap to visualize the different associations between bacterial genera and SCFAs shown in Fig. 7F and perform the principal component analysis shown in Fig. 7E. Both the heatmap, dendrogram, and principal component analysis were performed using ClustVist and default parameters [95].

To study if gut microbiota changes were less pronounced in patients who developed GVHD exclusively in the skin as compared to those who developed GVHD in the lower gastrointestinal tract, we calculated the difference between (i) the median log2FC (post-/pre-transplant) of the abundance of one particular bacterium in patients that were going to develop GVHD and (ii) the median log2FC of the abundance of that particular bacterium in patients that were not going to develop GVHD. This difference was calculated for all those taxa whose dynamics were found to be significantly different depending on the future GVHD status (Fig. 4D). Finally, the average and standard deviation of the results obtained for all taxa analyzed was calculated.

With the objective of identifying significant differences in the Shannon diversity index, ASV richness, in vitro growth of Staphylococci, the number of KOs increased/decreased post-transplant or network metrics, a parametric Student’s t-test or a non-parametric Mann–Whitney test was applied, as specified in each figure legend, depending if the populations under study followed or not a normal distribution which was analyzed applying the Shapiro–Wilk test. All statistical tests were two-sided unless we specified that the test applied was one-sided. One-sided tests were only applied when we had prior information about the expected direction of the microbiota change. For example, when the test was performed to confirm in a subset of the data (e.g., change associated with LGI-GVHD) a significant result that was obtained with the whole cohort (e.g., change associated with overall GVHD).

Data availability

All microbiome-generated data used for the analysis described in this manuscript, including information about bacterial taxa, counts at different levels (ASV, genus, family, order, class, phylum), KEGG orthologs, and bacterial species abundance derived from shotgun sequencing data and NMR spectra abundance are included in the Supplementary Table S11. Metadata used for the analysis performed is described in Supplementary Table S1. Raw sequencing data has been submitted to the SRA under the following accession number: PRJNA1088414.

References

Chabannon C, Kuball J, Bondanza A, et al. Hematopoietic stem cell transplantation in its 60s: A platform for cellula vr therapies. Sci Transl Med. 2018;10(436):eaap9630.
Article PubMed Google Scholar
Chang CC, Hayase E, Jenq RR. The role of microbiota in allogeneic hematopoietic stem cell transplantation. Expert Opin Biol Ther. 2021;21:1121–31.
Article CAS PubMed PubMed Central Google Scholar
Taur Y, Xavier JB, Lipuma L, et al. Intestinal domination and the risk of bacteremia in patients undergoing allogeneic hematopoietic stem cell transplantation. Clin Infect Dis. 2012;55(7):905–14.
Article CAS PubMed PubMed Central Google Scholar
Gratwohl A, Brand R, Frassoni F, et al. Cause of death after allogeneic haematopoietic stem cell transplantation (HSCT) in early leukaemias: an EBMT analysis of lethal infectious complications and changes over calendar time. Bone Marrow Transplant. 2005;36(9):757–69.
Article CAS PubMed Google Scholar
Ghimire S, Weber D, Mavin E, et al. Pathophysiology of GvHD and Other HSCT-Related Major Complications. Front Immunol. 2017;8:79.
Article PubMed PubMed Central Google Scholar
Ubeda C, Taur Y, Jenq RR, et al. Vancomycin-resistant Enterococcus domination of intestinal microbiota is enabled by antibiotic treatment in mice and precedes bloodstream invasion in humans. J Clin Investig. 2010;120:4332–41.
Article CAS PubMed PubMed Central Google Scholar
Djukovic A, José GM, Canlet C, et al. Lactobacillus supports Clostridiales to restrict gut colonization by multidrug- resistant Enterobacteriaceae. Nature Communications. 2022;13:5617.
Article CAS PubMed PubMed Central Google Scholar
Stein-Thoeringer CK, Nichols KB, Lazrak A, et al. Lactose drives Enterococcus expansion to promote graft-versus-host disease. Science. 2019;366:1143–9.
Article CAS PubMed PubMed Central Google Scholar
Shimasaki T, Seekatz A, Bassis C, et al. Increased relative abundance of Klebsiella pneumoniae carbapenemase-producing Klebsiella pneumoniae within the gut microbiota is associated with risk of bloodstream infection in long-term acute care hospital patients. Clin Infect Dis. 2018;68:2053–9.
Article PubMed Central Google Scholar
Weber D, Jenq RR, Peled JU, et al. Microbiota disruption induced by early use of broad-spectrum antibiotics is an independent risk factor of outcome after allogeneic stem cell transplantation. Biol Blood Marrow Transplant. 2017;23(5):845–52.
Article CAS PubMed PubMed Central Google Scholar
Simms-Waldrip TR, Sunkersett G, Coughlin LA, et al. Antibiotic-induced depletion of anti-inflammatory clostridia is associated with the development of graft-versus-host disease in pediatric stem cell transplantation patients. Biol Blood Marrow Transplant. 2017;23:820–9.
Article CAS PubMed Google Scholar
Mathewson ND, Jenq R, Mathew AV, et al. Gut microbiome-derived metabolites modulate intestinal epithelial cell damage and mitigate graft-versus-host disease. Nat Immunol. 2016;17:505–13.
Article CAS PubMed PubMed Central Google Scholar
Ilett EE, Jørgensen M, Noguera-Julian M, et al. Associations of the gut microbiome and clinical factors with acute GVHD in allogeneic HSCT recipients. Blood Adv. 2020;4:5797–809.
Article CAS PubMed PubMed Central Google Scholar
Golob JL, Pergam SA, Srinivasan S, et al. Stool microbiota at neutrophil recovery is predictive for severe acute graft vs host disease after hematopoietic cell transplantation. Clin Infect Dis. 2017;65(12):1984–91.
Article CAS PubMed PubMed Central Google Scholar
Biagi E, Zama D, Rampelli S, et al. Early gut microbiota signature of aGvHD in children given allogeneic hematopoietic cell transplantation for hematological disorders. BMC Méd Genom. 2019;12(1):49.
Article Google Scholar
Han L, Jin H, Zhou L, et al. Intestinal microbiota at engraftment influence acute graft-versus-host disease via the Treg/Th17 balance in Allo-HSCT recipients. Front Immunol. 2018;9:669.
Article PubMed PubMed Central Google Scholar
Payen M, Nicolis I, Robin M, et al. Functional and phylogenetic alterations in gut microbiome are linked to graft-versus-host disease severity. Blood Adv. 2020;4:1824–32.
Article CAS PubMed PubMed Central Google Scholar
Jenq RR, Taur Y, Devlin SM, et al. Intestinal blautia is associated with reduced death from graft-versus-host disease. Biology of blood and marrow transplantation : journal of the American Society for Blood and Marrow Transplantation. 2015;21:1373–83.
Article PubMed Google Scholar
Orberg ET, Meedt E, Hiergeist A, et al. Bacteria and bacteriophage consortia are associated with protective intestinal metabolites in patients receiving stem cell transplantation. Nat Cancer. 2024;5(1):187–208.
Article Google Scholar
Michonneau D, Latis E, Curis E, et al. Metabolomics analysis of human acute graft-versus-host disease reveals changes in host and microbiota-derived metabolites. Nat Commun. 2019;10(1):5695.
Vital M, Karch A, Pieper DH. Colonic butyrate-producing communities in humans: an overview using omics data. Msystems. 2017;2(6):e00130-e217.
Article CAS PubMed PubMed Central Google Scholar
Senn L, Clerc O, Zanetti G, et al. The stealthy superbug: the role of asymptomatic enteric carriage in maintaining a long-term hospital outbreak of st228 methicillin-resistant Staphylococcus aureus. mBio. 2016;7(1):e02039-15.
Article CAS PubMed PubMed Central Google Scholar
Yip AYG, King OG, Omelchenko O, et al. Antibiotics promote intestinal growth of carbapenem-resistant Enterobacteriaceae by enriching nutrients and depleting microbial metabolites. Nat Commun. 2023;14(1):5094.
Article CAS PubMed PubMed Central Google Scholar
Paredes-Sabja D, Setlow P, Sarker MR. Germination of spores of Bacillales and Clostridiales species: mechanisms and proteins involved. Trends Microbiol. 2011;19(2):85–94.
Article CAS PubMed Google Scholar
Romick-Rosendale LE, Haslam DB, Lane A, et al. Antibiotic exposure and reduced short chain fatty acid production after hematopoietic stem cell transplant. Biology of blood and marrow transplantation : journal of the American Society for Blood and Marrow Transplantation. 2018;24:2418–24.
Article CAS PubMed Google Scholar
Hayase E, Hayase T, Mukherjee A, et al. Bacteroides ovatus alleviates dysbiotic microbiota-induced graft-versus-host disease. Cell Host Microbe. 2024;32(9):1621-1636.e6.
Article CAS PubMed Google Scholar
Culp EJ, Goodman AL. Cross-feeding in the gut microbiome: Ecology and mechanisms. Cell Host Microbe. 2023;31(4):485–99.
Article CAS PubMed PubMed Central Google Scholar
Nielsen KL, Olsen MH, Pallejá A, et al. Microbiome compositions and resistome levels after antibiotic treatment of critically ill patients: an observational cohort study. Microorganisms. 2021;9(12):2542.
Article CAS PubMed PubMed Central Google Scholar
Pettigrew MM, Gent JF, Kong Y, et al. Gastrointestinal microbiota disruption and risk of colonization with carbapenem-resistant Pseudomonas aeruginosa in intensive care unit patients. Clin Infect Dis. 2018;69(4):604–13.
Article PubMed Central Google Scholar
Gu S-L, Gong Y, Zhang J, et al. Effect of the short-term use of fluoroquinolone and β-lactam antibiotics on mouse gut microbiota. Infect Drug Resist. 2020;13:4547–58.
Article CAS PubMed PubMed Central Google Scholar
Montassier E, Gastinne T, Vangay P, et al. Chemotherapy-driven dysbiosis in the intestinal microbiome. Aliment Pharmacol Ther. 2015;42(5):515–28.
Article CAS PubMed Google Scholar
Maier L, Goemans CV, Wirbel J, et al. Unravelling the collateral damage of antibiotics on gut bacteria. Nature. 2021;599:120–4.
Article CAS PubMed PubMed Central Google Scholar
Nakamura A, Hosoda M, Kato T, et al. Combined effects of meropenem and aminoglycosides on Pseudomonas aeruginosa in vitro. J Antimicrob Chemother. 2000;46(6):901–4.
Article CAS PubMed Google Scholar
Bayatinejad G, Salehi M, Beigverdi R, et al. In Vitro antibiotic combinations of Colistin, Meropenem, Amikacin, and Amoxicillin/clavulanate against multidrug-resistant Klebsiella pneumonia isolated from patients with ventilator-associated pneumonia. BMC Microbiol. 2023;23(1):298.
Article CAS PubMed PubMed Central Google Scholar
Brook I. Synergistic combinations of antimicrobial agents against anaerobic bacteria. Pediatr Infect Dis J. 1987;6(3):332–5.
Article CAS PubMed Google Scholar
Zeng MY, Inohara N, Nuñez G. Mechanisms of inflammation-driven bacterial dysbiosis in the gut. Mucosal Immunol. 2017;10(1):18–26.
Article CAS PubMed Google Scholar
Miller WD, Keskey R, Alverdy JC. SePSIS AND THE MICROBIOME: A VICIOUS CYCLE J Infect Dis. 2020;223(Supplement_3):S264–9.
Article PubMed Central Google Scholar
Tanaka JS, Young RR, Heston SM, et al. Anaerobic antibiotics and the risk of graft-versus-host disease after allogeneic hematopoietic stem cell transplantation. Biology of blood and marrow transplantation : journal of the American Society for Blood and Marrow Transplantation. 2020;26:2053–60.
Article CAS PubMed Google Scholar
Messaoudene M, Ferreira S, Saint-Lu N, et al. The DAV132 colon-targeted adsorbent does not interfere with plasma concentrations of antibiotics but prevents antibiotic-related dysbiosis: a randomized phase I trial in healthy volunteers. Nat Commun. 2024;15(1):8083.
Article CAS PubMed PubMed Central Google Scholar
Balletto E, Mikulska M. Bacterial infections in hematopoietic stem cell transplant recipients. Mediterr J Hematol Infect Dis. 2014;7(1):e2015045–e2015045.
Google Scholar
von Eiff C, Becker K, Machka K, Stammer H, Peters G. Nasal carriage as a source of Staphylococcus aureus Bacteremia. N Engl J Med. 2001;344(1):11–6.
Article Google Scholar
Srinivasan A, Seifried SE, Zhu L, et al. Increasing prevalence of nasal and rectal colonization with methicillin-resistant Staphylococcus aureus in children with cancer. Pediatr Blood Cancer. 2010;55(7):1317–22.
Article PubMed PubMed Central Google Scholar
Bhalla A, Aron DC, Donskey CJ. Staphylococcus aureus intestinal colonization is associated with increased frequency of S aureus on skin of hospitalized patients. BMC Infect Dis. 2007;7(1):105–105.
Article PubMed PubMed Central Google Scholar
Squier C, Rihs JD, Risa KJ, et al. Staphylococcus aureus rectal carriage and its association with infections in patients in a surgical intensive care unit and a liver transplant unit. Infect Control Hosp Epidemiology. 2002;23(9):495–501.
Article Google Scholar
van der Hee B, Wells JM. Microbial regulation of host physiology by short-chain fatty acids. Trends Microbiol. 2021;29(8):700–12.
Article PubMed Google Scholar
Liu X, Mao B, Gu J, et al. Blautia—a new functional genus with potential probiotic properties? Gut Microbes. 2021;13(1):1875796.
Article PubMed PubMed Central Google Scholar
Sorbara MT, Dubin K, Littmann ER, et al. Inhibiting antibiotic-resistant Enterobacteriaceae by microbiota-mediated intracellular acidification. J Exp Med. 2019;216(1):84–98.
Article CAS PubMed PubMed Central Google Scholar
Han SH. Propionate ameliorates staphylococcus aureus skin infection by attenuating bacterial growth. Front Microbiol. 2019;10:1–13.
Google Scholar
Wu W, Sun M, Chen F, et al. Microbiota metabolite short-chain fatty acid acetate promotes intestinal IgA response to microbiota which is mediated by GPR43. Mucosal Immunol. 2017;10:946–56.
Article CAS PubMed Google Scholar
Gupta S, Basu S, Bal V, Rath S, George A. Gut IgA abundance in adult life is a major determinant of resistance to dextran sodium sulfate-colitis and can compensate for the effects of inadequate maternal IgA received by neonates. Immunology. 2019;158(1):19–34.
Article CAS PubMed PubMed Central Google Scholar
Scheidler L, Hippe K, Ghimire S, et al. Intestinal IgA positive plasma cells are highly sensitive indicators of alloreaction early after allogeneic transplantation and associate with both, graft-versus-host disease and relapse related mortality. Haematologica. 2023;108(11):2993–3000.
Article CAS PubMed PubMed Central Google Scholar
Khan MT, Duncan SH, Stams AJM, et al. The gut anaerobe Faecalibacterium prausnitzii uses an extracellular electron shuttle to grow at oxic–anoxic interphases. ISME J. 2012;6(8):1578–85.
Article CAS PubMed PubMed Central Google Scholar
Zheng H, Yde CC, Clausen MR, et al. Metabolomics investigation to shed light on cheese as a possible piece in the french paradox puzzle. J Agric Food Chem. 2015;63(10):2830–9.
Article CAS PubMed Google Scholar
Sun B, Ma T, Li Y, et al. Bifidobacterium lactis Probio-M8 adjuvant treatment confers added benefits to patients with coronary artery disease via target modulation of the gut-heart/-brain axes. mSystems. 2022;7(2):e00100-22.
Article PubMed PubMed Central Google Scholar
Mills EL, Kelly B, Logan A, et al. Succinate dehydrogenase supports metabolic repurposing of mitochondria to drive inflammatory macrophages. Cell. 2016;167(2):457-470.e13.
Article CAS PubMed PubMed Central Google Scholar
Jardine L, Cytlak U, Gunawan M, et al. Donor monocyte-derived macrophages promote human acute graft versus host disease. J Clin Investig. 2020;130(9):4574–86.
Article CAS PubMed PubMed Central Google Scholar
Legoff J, Resche-Rigon M, Bouquet J, et al. The eukaryotic gut virome in hematopoietic stem cell transplantation: new clues in enteric graft-versus-host disease. Nat Med. 2017;23:1080–5.
Article CAS PubMed Google Scholar
Atarashi K, Tanoue T, Shima T, et al. Induction of colonic regulatory T cells by indigenous clostridium species. Science. 2011;331:337–41.
Article CAS PubMed Google Scholar
Ito R, Katano I, Kawai K, et al. A novel xenogeneic graft-versus-host disease model for investigating the pathological role of human CD4+ or CD8+ T cells using immunodeficient NOG mice. Am J Transplant. 2017;17(5):1216–28.
Article CAS PubMed Google Scholar
Ito R, Katano I, Otsuka I, et al. Exacerbation of pathogenic Th17-cell-mediated cutaneous graft-versus-host-disease in human IL-1β and IL-23 transgenic humanized mice. Biochem Biophys Res Commun. 2019;516(2):480–5.
Article CAS PubMed Google Scholar
Cheng H, Tian J, Zeng L, et al. Halofugine prevents cutaneous graft versus host disease by suppression of Th17 differentiation. Hematology. 2012;17(5):261–7.
Article CAS PubMed Google Scholar
Yi T, Chen Y, Wang L, et al. Reciprocal differentiation and tissue-specific pathogenesis of Th1, Th2, and Th17 cells in graft-versus-host disease. Blood. 2009;114(14):3101–12.
Article CAS PubMed PubMed Central Google Scholar
Ullrich E, Abendroth B, Rothamer J, et al. BATF-dependent IL-7RhiGM-CSF+ T cells control intestinal graft-versus-host disease. J Clin Investig. 2018;128(3):916–30.
Article PubMed PubMed Central Google Scholar
Ivanov II, Atarashi K, Manel N, et al. Induction of intestinal Th17 cells by segmented filamentous bacteria. Cell. 2009;139:485–98.
Article CAS PubMed PubMed Central Google Scholar
Wu H-J, Ivanov II, Darce J, et al. Gut-residing segmented filamentous bacteria drive autoimmune arthritis via T helper 17 cells. Immunity. 2010;32:815–27.
Article CAS PubMed PubMed Central Google Scholar
Martins MC, Romão CV, Folgosa F, et al. How superoxide reductases and flavodiiron proteins combat oxidative stress in anaerobes. Free Radic Biol Med. 2019;140:36–60.
Article CAS PubMed Google Scholar
Schwab L, Goroncy L, Palaniyandi S, et al. Neutrophil granulocytes recruited upon translocation of intestinal bacteria enhance graft-versus-host disease via tissue damage. Nat Med. 2014;20(6):648–54.
Article CAS PubMed Google Scholar
Highfill SL, Rodriguez PC, Zhou Q, et al. Bone marrow myeloid-derived suppressor cells (MDSCs) inhibit graft-versus-host disease (GVHD) via an arginase-1–dependent mechanism that is up-regulated by interleukin-13. Blood. 2010;116(25):5738–47.
Article CAS PubMed PubMed Central Google Scholar
Li T-T, Chen X, Huo D, et al. Microbiota metabolism of intestinal amino acids impacts host nutrient homeostasis and physiology. Cell Host Microbe. 2024;32(5):661-675.e10.
Article CAS PubMed Google Scholar
Coyte KZ, Schluter J, Foster KR. The ecology of the microbiome: networks, competition, and stability. Science. 2015;350(6261):663–6.
Article CAS PubMed Google Scholar
Bacigalupo A, Ballen K, Rizzo D, et al. Defining the intensity of conditioning regimens: working definitions. Biol Blood Marrow Transplant. 2009;15(12):1628–33.
Article PubMed PubMed Central Google Scholar
Rowlings PA, Przepiorka D, Klein JP, et al. IBMTR Severity Index for grading acute graft-versus-host disease: retrospective comparison with Glucksberg grade. Br J Haematol. 1997;97(4):855–64.
Article CAS PubMed Google Scholar
Isaac S, Flor-duro A, Carruana G, et al. Microbiome-mediated fructose depletion restricts murine gut colonization by vancomycin-resistant Enterococcus. Nature Communications. 2022;13:7718.
Article CAS PubMed PubMed Central Google Scholar
Callahan BJ, McMurdie PJ, Rosen MJ, et al. DADA2: high-resolution sample inference from Illumina amplicon data. Nat Methods. 2016;13:581–3.
Article CAS PubMed PubMed Central Google Scholar
Yogev N, Bedke T, Kobayashi Y, et al. CD4+ T-cell-derived IL-10 promotes CNS inflammation in mice by sustaining effector T cell survival. Cell Rep. 2022;38: 110565.
Article CAS PubMed Google Scholar
Quast C, Pruesse E, Yilmaz P, et al. The SILVA ribosomal RNA gene database project: Improved data processing and web-based tools. Nucleic Acids Research. 2013;41:D590-6.
Article CAS PubMed Google Scholar
Chen S, Zhou Y, Chen Y, Gu J. fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics. 2018;34(17):i884–90.
Article PubMed PubMed Central Google Scholar
Langmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nature Methods. 2012;9(4):357–9.
Article CAS PubMed PubMed Central Google Scholar
Tamames J, Puente-Sánchez F. SqueezeMeta, A highly portable, fully automatic metagenomic analysis pipeline. Front Microbiol. 2019;9:3349.
Article PubMed PubMed Central Google Scholar
Ogata H, Goto S, Sato K, et al. KEGG: Kyoto Encyclopedia of Genes and Genomes. Nucleic Acids Res. 1999;27:29–34.
Article CAS PubMed PubMed Central Google Scholar
Artacho A, Isaac S, Nayak R, et al. The pre-treatment gut microbiome is associated with lack of response to methotrexate in new onset rheumatoid arthritis. Arthritis & Rheumatology. 2020;73(6):931–42.
Meiboom S, Gill D. Modified spin-echo method for measuring nuclear relaxation times. Rev Sci Instrum. 1958;29(8):688–91.
Article CAS Google Scholar
Foxall PJD, Spraul M, Farrant RD, et al. 750 MHz 1H-NMR spectroscopy of human blood plasma. J Pharm Biomed Anal. 1993;11(4–5):267–76.
Article CAS PubMed Google Scholar
Markley JL, Ulrich EL, Berman HM, et al. BioMagResBank (BMRB) as a partner in the Worldwide Protein Data Bank (wwPDB): new policies affecting biomolecular NMR depositions. J Biomol NMR. 2008;40(3):153–5.
Article CAS PubMed PubMed Central Google Scholar
Wishart DS, Feunang YD, Marcu A, et al. HMDB 40: the human metabolome database for 2018. Nucleic Acids Res. 2018;46(Database issue):D608–17.
Article CAS PubMed Google Scholar
de Azavedo JC, Foster TJ, Hartigan PJ, et al. Expression of the cloned toxic shock syndrome toxin 1 gene (tst) in vivo with a rabbit uterine model. Infect Immun. 1985;50(1):304–9.
Article PubMed PubMed Central Google Scholar
Yamamura R, Inoue KY, Nishino K, Yamasaki S. Intestinal and fecal pH in human health. Front Microbiomes. 2023;2:1192316.
Article Google Scholar
Lin H, Peddada SD. Analysis of compositions of microbiomes with bias correction. Nat Commun. 2020;11(1):3514.
Article CAS PubMed PubMed Central Google Scholar
Benjamini Y, Hockberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J Roy Stat Soc. 1995;57:289–300.
Article Google Scholar
Kursa MB, Rudnicki WR. Feature selection with the Boruta Package. J Stat Softw. 2010;36(11). https://doiorg.publicaciones.saludcastillayleon.es/10.18637/jss.v036.i11.
Worley B, Powers R. Multivariate analysis in metabolomics. Curr Metabolomics. 2013;1(1):92–107.
CAS PubMed PubMed Central Google Scholar
Jonsson P, Wuolikainen A, Thysell E, et al. Constrained randomization and multivariate effect projections improve information extraction and biomarker pattern discovery in metabolomics studies involving dependent samples. Metabolomics. 2015;11(6):1667–78.
Article CAS PubMed PubMed Central Google Scholar
Rohart F, Gautier B, Singh A, Cao K-AL. mixOmics: An R package for ‘omics feature selection and multiple data integration. PLoS Comput Biol. 2017;13(11):e1005752.
Csardi G, Nepusz T. The igraph software package for complex network research. InterJournal of Complex Systems. 2006;1695(5):1–9.
Metsalu T, Vilo J. ClustVis: a web tool for visualizing clustering of multivariate data using Principal Component Analysis and heatmap. Nucleic Acids Res. 2015;43(Web Server issue):W566–W570.

Download references

Acknowledgements

We would like to thank all the participants who agreed to participate in this study. We also thank R. Jenq for the helpful comments on the manuscript.

Funding

C.U. was supported by the InfectERA-ERANET-Acciones complementarias grant [PCIN-2015–094] from the Spanish Ministerio de Economía y Competitividad and the 7th Research framework program from EU, grants from Conselleria d’Innovació, Universitats, Ciència i Societat Digital [CIPROM/2021/053], a grant from the Spanish MICINN [PID2020-120292RB-I00] and an intramural grant from FISABIO. JS was supported by Instituto de Salud Carlos III through the project “PI10/00280” (Co-funded by the European Regional Development Fund “A way to make Europe”). P.M.P was supported with a post-residency contract (2015/0369) from Instituto de Investigación Sanitaria del Hospital Universitari y Politècnic La Fe and with a grant to emerging groups (GV/2016/066) from Conselleria de Educación, Cultura y Deporte de la Generalitat Valenciana. L.P.C was supported by Instituto de Salud Carlos III through a Miguel Servet contract (CP22/00005), co-funded by the European Union.

Author information

Alejandro Artacho and Cintya González-Torres contributed equally to this manuscript.

Authors and Affiliations

Fundación Para El Fomento de La Investigación Sanitaria y Biomédica de La Comunitat Valenciana-FISABIO, Valencia, Spain
Alejandro Artacho, Cintya González-Torres, Javier Pons, Nuria Jiménez & Carles Ubeda
Drug Discovery Unit, Instituto de Investigación Sanitaria La Fe, Valencia, Spain
Nuria Gómez-Cebrián & Leonor Puchades-Carrasco
Hematology Department, Hospital Universitari I Politècnic La Fe, Valencia, Spain
Paula Moles-Poveda, María Jinglei Casanova, Juan Montoro, Aitana Balaguer, Marta Villalba, Pedro Chorão & Jaime Sanz
Departament de Medicina, Universitat de Valencia, Valencia, Spain
Jaime Sanz
CIBERONC, Instituto Carlos III, Madrid, Spain
Jaime Sanz
Centers of Biomedical Research Network (CIBER) in Epidemiology and Public Health, Madrid, Spain
Carles Ubeda

Authors

Alejandro Artacho
View author publications
You can also search for this author inPubMed Google Scholar
Cintya González-Torres
View author publications
You can also search for this author inPubMed Google Scholar
Nuria Gómez-Cebrián
View author publications
You can also search for this author inPubMed Google Scholar
Paula Moles-Poveda
View author publications
You can also search for this author inPubMed Google Scholar
Javier Pons
View author publications
You can also search for this author inPubMed Google Scholar
Nuria Jiménez
View author publications
You can also search for this author inPubMed Google Scholar
María Jinglei Casanova
View author publications
You can also search for this author inPubMed Google Scholar
Juan Montoro
View author publications
You can also search for this author inPubMed Google Scholar
Aitana Balaguer
View author publications
You can also search for this author inPubMed Google Scholar
Marta Villalba
View author publications
You can also search for this author inPubMed Google Scholar
Pedro Chorão
View author publications
You can also search for this author inPubMed Google Scholar
Leonor Puchades-Carrasco
View author publications
You can also search for this author inPubMed Google Scholar
Jaime Sanz
View author publications
You can also search for this author inPubMed Google Scholar
Carles Ubeda
View author publications
You can also search for this author inPubMed Google Scholar

Contributions

A.A., J.P., M.J.C., C.G.T and C.U. performed sequencing analysis. N.J. processed samples and obtained sequences. N.G.C and L.P.C performed metabolomic analysis. C.G.T performed in vitro experiments. P.M.P, J.M., A.B., M.V., P.C. and J.S. collected faecal samples and obtained and analysed metadata from patients. A.A., C.U and C.G.T performed statistical analysis. J.S. and C.U. designed research studies and wrote the manuscript.

Corresponding authors

Correspondence to Jaime Sanz or Carles Ubeda.

Ethics declarations

Ethics approval and consent to participate

The study was approved by the institutional review board at the Hospital La Fe (2015/0369) and the Ethics Committee of CEIC Dirección General de Salud Pública y Centro Superior de Investigación en Salud Pública (20130515/08). All patients gave their consent to participate in this study.

Consent for publication

Not applicable. No data on individual persons is reported.

Competing interests

The authors declare no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

40168_2024_1948_MOESM1_ESM.pdf

Supplementary Material 1: Figure S1. 16s rRNA gene sequencing depth was sufficient to identify most ASVs in all samples. Supplementary Figure S2. Multiple Staphylococcus species expand during allo-HSCT. Supplementary Figure S3. Antibiotic treatment is associated with microbiota changes detected after allo-HSCT. Supplementary Figure S4. Taxa encoding KOs within the “Cell growth” pathway whose abundance significantly decrease after allo-HSCT. Supplementary Figure S5. Taxa encoding KOs within the final step of the pathway of synthesis of SCFAs. Supplementary Figure S6. Overall microbiota changes after allo-HSCT do not differentiate patients according to their future GVHD status. Supplementary Figure S7. Microbiota dynamics classify to a better extent patients that will develop GVHD. Supplementary Figure S8. Specific microbiota taxonomic changes are detected in patients that will develop GVHD in the lower gastrointestinal tract. Supplementary Figure S9. Specific microbiota taxonomic changes are associated with exclusively developing GVHD in the skin. Supplementary Figure S10. Bacterial taxonomic changes correlate with GVHD severity. Supplementary Figure S11. Changes in butyrate and propionate after alloHSCT according to the future GVHD status. Supplementary Figure S12. Samples from GVHD and non-GVHD patients do not differ in their overall metabolome changes according to OPLS-EP analysis. Supplementary Figure S13. Specific metabolomic changes are detected in patients that will develop GVHD in the lower gastrointestinal tract. Supplementary Figure S14. Taxa encoding the K05919 and the K18918. Supplementary Figure S15. Changes in specific KOs are associated with GVHD development in the lower gastrointestinal tract. Supplementary Figure S16. Changes in specific KOs are associated with the severity of GVHD. Supplementary Figure S17. Clostridium XIVb is negatively associated with Faecalibacterium and positively associated with butyrate in post-transplant GVHD samples

Supplementary Material 2: Table S1

Supplementary Material 3: Table S2

Supplementary Material 4: Table S3

Supplementary Material 5: Table S4

Supplementary Material 6: Table S5

Supplementary Material 7: Table S6

Supplementary Material 8: Table S7

Supplementary Material 9: Table S8

Supplementary Material 10: Table S9

Supplementary Material 11: Table S10

Supplementary Material 12: Table S11

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.

Reprints and permissions

About this article

Cite this article

Artacho, A., González-Torres, C., Gómez-Cebrián, N. et al. Multimodal analysis identifies microbiome changes linked to stem cell transplantation-associated diseases. Microbiome 12, 229 (2024). https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s40168-024-01948-0

Download citation

Received: 14 March 2024
Accepted: 11 October 2024
Published: 07 November 2024
DOI: https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s40168-024-01948-0

Multimodal analysis identifies microbiome changes linked to stem cell transplantation-associated diseases

Abstract

Background

Results

Conclusions

Introduction

Results

Taxonomic changes associated with allo-HSCT

Metabolomic changes associated with allo-HSCT

Changes in encoded bacterial gene orthologs during allo-HSCT

Taxonomic changes associated with the future development of GVHD

Taxonomic changes associated with the grade of GVHD and the organs affected

Metabolomic changes associated with the future development of GVHD

Changes in encoded bacterial gene orthologs associated with the future development of GVHD

The metabolome-bacterial interactome distinguishes patients according to their future GVHD status

Discussion

Methods

Study design and sample collection

Fecal DNA extraction

16S rRNA sequencing and analysis

Metagenomic shotgun sequencing

Metabolomic analysis

In vitro culture experiments

Statistical analysis

Data availability

References

Acknowledgements

Funding

Author information

Authors and Affiliations

Contributions

Corresponding authors

Ethics declarations

Ethics approval and consent to participate

Consent for publication

Competing interests

Additional information

Publisher’s Note

Supplementary Information

Rights and permissions

About this article

Cite this article

Share this article

Microbiome

Contact us