Pre-Transplant Prediction of Acute Graft-versus-Host Disease Using the Gut Microbiome

Gut microbiota is thought to influence host responses to allogeneic hematopoietic stem cell transplantation (aHSCT). Recent evidence points to this post-transplant for acute graft-versus-host disease (aGvHD). We asked whether any such association might be found pre-transplant and conducted a metagenome-wide association study (MWAS) to explore. Microbial abundance profiles were estimated using ensembles of Kaiju, Kraken2, and DeepMicrobes calls followed by dimensionality reduction. The area under the curve (AUC) was used to evaluate classification of the samples (aGvHD vs. none) using an elastic net to test the relevance of metagenomic data. Clinical data included the underlying disease (leukemia vs. other hematological malignancies), recipient age, and sex. Among 172 aHSCT patients of whom 42 developed aGVHD post transplantation, a total of 181 pre-transplant tool samples were analyzed. The top performing model predicting risk of aGVHD included a reduced species profile (AUC = 0.672). Beta diversity (37% in Jaccard’s Nestedness by mean fold change, p < 0.05) was lower in those developing aGvHD. Ten bacterial species including Prevotella and Eggerthella genera were consistently found to associate with aGvHD in indicator species analysis, as well as relief and impurity-based algorithms. The findings support the hypothesis on potential associations between gut microbiota and aGvHD based on a data-driven approach to MWAS. This highlights the need and relevance of routine stool collection for the discovery of novel biomarkers.


Introduction
Acute graft-versus-host disease (aGvHD) is a frequent serious complication in patients undergoing allogeneic hematopoietic stem cell transplantation (aHSCT) [1]. In aGvHD, immune cells from the donor attack healthy tissue of the patient, causing an increased risk of morbidity and mortality. aHSCT offers a potential cure for patients suffering from hematological diseases (most commonly leukemia), where other treatment options have been exhausted or are not suitable to tackle the patient's aggressive disease. Finding ways to predict aGvHD prior to aHSCT in recipients is an important and still open research question [2].
With the advent of high throughput sequencing the potential role of the gut microbiome in human disease can be better understood. It is suggested that the risk of developing aGvHD could, in part, be associated with the recipients' genetic characteristics and/or cohabitating microorganisms within the gut microbiome [3]. The gut microbiome offers an intriguing source of information having recently been suggested to improve aGvHD prediction models [4,5]. Previous studies, mainly using post-transplant samples, have found that abundancies of specific gut bacteria may provide useful information in predicting aGvHD [6]. Likewise, several gut bacteria have been associated with an increased risk of aGvHD and mortality [2]. Based on previous findings suggesting associations between the gut microbiome and aGvHD [5,7,8], we hypothesized that potential pre-aHSCT biomarkers of aGvHD could be discovered through a metagenome-wide association study (MWAS).
To test this hypothesis, we focused on species abundances with respect to prospective aGvHD occurrence. We used three taxonomic classifiers: DeepMicrobes [9], Kraken2 [10], and Kaiju [11] based on a custom database of the human gut bacteria to estimate species abundances using high throughput shotgun sequencing data from stool samples of aHSCT recipients. There was a pre-processing step to choose a configuration to aggregate metagenomic data derived from the taxonomic classifiers for downstream analyses. We used a broad set of statistical and information-theoretic methods to assess potential associations between gut species abundances and aGvHD.

Data Collection and Categorization
Stool samples were collected by the patient or nursing staff using the OMNIgene.GUT (DNA Stabilized-frozen Inc., Ottawa, Canada) stabilization tube according to the manufacturer's instructions. All samples were frozen a maximum of eight days after sampling and stored at −80 • C until shipment for sequencing. Paired-end sequencing (150 bp) was done using shotgun metagenomic sequencing-Illumina HiSeq (Illumina, San Diego, CA, USA). The sequencing reads underwent pre-processing and quality control steps. Using Trimmomatic [12] and the Illumina HiSeq program, certain adapters were eliminated. Additionally, Trimmomatic was used for quality trimming, that is, reads were cut if the average quality of 4 nucleotides was below 30 and nucleotides were trimmed from the start or end of the read if they had a quality score of less than 30. It also included removal of reads with <50 base pairs and removal of all reads mapping to the human genome (GRCh38) using Bowtie2 [13]. Samples with <1 million reads remaining after quality control steps were removed from further analysis. Samples were categorized into two groups depending on whether the patient developed aGvHD (grades 2 and above) or not, post transplant (i.e., an aGvHD group and a non-aGvHD (NaGvHD) group, respectively). A graphical summary of the cohort and methodology of the study is provided in Figure 1.

Data Analysis
The data analyses for MWAS were performed in R [14] (v4.1.2). The analyses included exploring potential associations between clinical variables including prospective aGvHD occurrence and (normalized) species abundances. Statistical significance was considered throughout this work at a level of 5% unless otherwise noted. Due to the limited sample size and our purpose of exploring a large number of candidate species in some of the analyses, the p values before adjustments for multiple testing are also reported.  Figure 1. An overview of the study depicting the data distribution for stool samples with transplantation day as reference (day 0). It also depicts the focus of the study and the methodology to benchmark taxonomic profiling configurations to describe associations between gut microbes before allogeneic hematopoietic stem cell transplantation (aHSCT) and occurrence of acute graft versus host disease (aGvHD). The blue block represents the taxonomic classifiers and their combined taxonomic profiles extracting species counts. Species counts are normalized, and two sets of species are analyzed: (1) unfiltered or complete species set and (2) filtered or reduced species set. The species set with the highest classification performance is chosen for metagenome-wide association study (the block in light gold color).

Clinical Variables
We used univariable and multivariable logistic regression to find statistically significant associations of the clinical variables, i.e., conditioning regimen (myeloablative or non-myeloablative), cyclophosphamide (used or not), disease groups (acute leukemia or other hematological malignancies), donor-recipient relationship (related or unrelated), donor sex, Karnofsky score at day 0 (<90 or ≥90), radiation groups (none or irradiated), recipient age (<60 or ≥60), recipient sex, transplant source (peripheral blood or bone marrow) with aGvHD as a binary outcome variable (NaGvHD = 0, aGvHD = 1). Myeloablative is a harsher form of conditioning including high-dose cytotoxic treatment, resulting in patients becoming neutropenic in the first few weeks post-transplant, whereas non-myeloablative is less harsh, resulting in higher levels of neutrocytes post-transplant [15]. For the analysis, Karnofsky score at day 0 was imputed for two patients with the median value of it across patients with the same outcome. The dichotomization of Karnofsky score and recipient age was applied to balance their distributions by aGvHD for potential deviation from normality.

Taxonomic Classification
Three taxonomic classifiers, i.e., DeepMicrobes (v1.0.0), Kraken2 (v2.1.2), and Kaiju (v1.7.3) were used to assign the species from a relevant custom reference database [9] to the reads involved in each metagenomic sample. The database contains 2505 bacterial species from a bacterial repertoire of human gut microbiota. These taxonomic classifiers use different methods to assign the species to reads. The main difference is that Deep-Microbes is based on deep learning and does not require phylogenetic inference for the classification as a key advantage [9]. We extracted three species abundance profiles given the three taxonomic classifiers that were used for the downstream analyses. The taxonomic classifiers provided species counts in each sample. The species counts from the three taxonomic classifiers were aggregated by arithmetic mean (AM) leading to four extra sets: AM(DeepMicrobes, Kraken2), AM(DeepMicrobes, Kaiju), AM(Kraken2, Kaiju), and AM(DeepMicrobes, Kraken2, Kaiju). The species counts were normalized to their genome sizes and total number of reads to compute (relative) species abundances to allow their comparisons across samples (see Supplementary Material). The methodological variance unravels potential effects of taxonomic classification on distinguishing samples labeled as aGvHD from NaGvHD. As a naming protocol, the name of the most recent taxa to each species is indicated followed by "|" and the taxa identifier from the reference database. The initials refer to the taxa ranks from which the name comes from (if any): s_: species, g_: genus, o_ order, f_:family, p_:phylum, c_:class.

Benchmarking the Taxonomy Profiling Configurations
The samples were classified to aGvHD and NaGvHD using an elastic net in a 2 × 5fold cross validation scheme in which the area under the curve (AUC) for receiver operating characteristic (AUROC) and precision-recall curves are reported to assess the predictability of aGvHD using each of the seven taxonomy profiling configurations. The classification was performed using SIAMCAT package (v1.14.0). We applied PERFect: PERmutation Filtering test [16], an unsupervised dimension reduction method, to find a reduced (filtered) species set in comparison with the complete (unfiltered) species set (2505 species from the database) to check whether the filtering increases signal to noise ratios quantified by AUC. This determines the taxonomy profiling configuration with the highest AUC to present the results for the downstream analyses.

Diversity Measures
We computed seven alpha and 10 beta diversity indices derived from species abundances as proxies to explore the patients' gut ecosystems. Alpha diversity indices specify how diverse the species communities are within each sample, whereas beta diversity indices quantify the community diversities between samples of the aGvHD and NaGvHD groups. Alpha diversity indices were Berger-Parker dominance, Simpson's dominance, Heip's evenness, inverse Simpson, richness, Shannon diversity, and Strong's dominance as implemented in abdiv package (v0.2.0). Beta diversity indices were Canberra distance, Chebyshev distance, Clark's coefficient of divergence, correlation distance, Euclidean distance, Geodesic metric, Hellinger distance, Horn-Morisita index, Jaccard nestedness, and Manhattan distance. We also computed F 0 [17] as a k-mer based measure of alpha diversity in addition to the mentioned alpha diversity indices. The differences in the indices were tested for statistical significance using Wilcoxon rank-sum test followed by BH adjustment [18] of the p values. To assess bias possibility in group difference test of significance for beta diversity indices due to different sample sizes, the test was conducted also for 1000 random subsets of the group with larger sample size against the group with smaller sample size (group-equalized testing).

Dissimilarity Measures
Species abundances were used to compare microbial dissimilarities in association with the occurrence of aGvHD. The distances between all pairs of samples were computed using vegdist function for "robust Aitchison" and "Bray-Curtis" measures. Analysis of similarities (ANOSIM) test using vegan package (v2.6.2) was conducted with 10,000 permutations for statistical significance of the dissimilarity measures. Permutational multivariate analysis of variance (PERMANOVA) using adonis2 function from vegan package (v2.6.2) was also performed with 10,000 permutations to assess the between-group dissimilarity of the relative abundances of samples based on the distance measures. The contributions of disease group, recipient age and sex to bacterial community composition were determined with PERMANOVA. The betadisper function was used to test the assumption of homogeneity of multivariate dispersion. To partly compensate for the sensitivity to unbalanced designs [19], the PERMANOVA and ANOSIM tests with robust Aitchison distance were also conducted for an adjusted balance setting based on Majority Weighted Minority Oversampling TEchnique (MWMOTE) algorithm [20] implemented by R package imbalance (v1.0.2.1), in which, oversampling is applied to the group with smaller sample size to generate synthetic samples resembling the distributions of the original samples.
In addition, indicator species analysis was performed using indicspecies package (v1.7.12) on species abundances to identify the species available more frequently in one group (aGvHD) compared to NaGvHD [21]. Using this analysis, group-equalized pointbiserial correlation coefficients (r.g) and indicator value (IndVal.g) were used to identify the species that could be used to distinguish groups [22,23]. IndVal.g quantifies associations between aGvHD and species abundances by indicating both fidelity (sensitivity or coverage within a group) and specificity (being exclusive to a group). The specificity and fidelity are two quantities acting as probabilities ranging from 0 to 1 in which the higher the value represents higher probability of a species being exclusively found in a specific group and the species being found in all samples of the group. Furthermore, a permutation-based test (10,000 permutations) implemented by signassoc function was applied to identify the species with higher abundance in the aGvHD group compared to NaGvHD. The adjusted p values by BH or Sidak's method for multiple testing are also reported.

Information-Theoretic Measures
Mutual information as implemented in infotheo package (v1.2.0.1) was used to quantify the association between species abundances and the occurrence of aGvHD conditioned by statistically significant confounders determined by PERMANOVA. This quantity allows measuring an estimated amount of information species abundances may contain with respect to the outcome given that the significant confounders are known.

Relief and Impurity-Based Algorithms for Metagenomic Feature Ranking
We used 13 feature ranking (selection) algorithms to identify most relevant features (species) with respect to our outcome (aGvHD vs. NaGvHD). The algorithms are based on extensions to Relief [24], impurity [25] or their combinations all implemented in CORElearn package (v1.56.0). The Relief-based algorithms differ on parameters such as number of neighbors (K) and strategies to calculate feature importance [26], all of which, lead to the following algorithms: equally assigned weights to K nearest neighbors (ReliefFequalK), exponentially decreasing weights with increasing ranks (ReliefFexpRank), K selection per feature (ReliefFbestK), myopic ReliefF (MyopicReliefF), merit-based ReliefF (ReliefFmerit), and weighting by inverse distance to selected instance (ReliefFdistance) in addition to the original Relief algorithm (Relief). The impurity-based algorithms are gain ratio (GainRatio), minimum description length (MDL), DKM index named after its authors (Dietterich, Kearns, and Mansour), information gain with uniform priors (UniformInf), Hellinger distance with equal weighting (EqualHellinger), and Euclidean distance (DistEuclid). The median rank of each species across algorithms was used to find the most relevant species in association with aGvHD.

Cohort Characteristics: Clinical and Metagenomic Data
Our cohort consisted of 181 pre-transplant stool samples from 172 patients receiving aHSCT between January 2016, and August 2020 at Rigshospitalet, Copenhagen University Hospital, Denmark ( Figure 1). Among the 181 samples, 42 were from patients who later developed aGvHD (grade II-IV according to the modified Glucksberg-Seattle criteria [27,28]) with a median (Q1-Q3) time from aHSCT to aGvHD of 33 (27-55) days within 100 days of post-aHSCT, where Q1 is the first quartile and Q3 is the third quartile. Of the 181 samples, 172 (95%) were collected within 30 days pre-aHSCT, with the median (Q1-Q3) sampling time being 20 (14)(15)(16)(17)(18)(19)(20)(21)(22) days prior to aHSCT. For an overview of patient characteristics, see Supplementary Table S1. There were nine patients (three who later developed aGvHD) who had two samples each, median 0 (range 0-61) days in between. Here, the sample closest to aHSCT was selected. For the patients with two samples on the same day, the samples were merged by the average of their normalized species abundances.

Clinical Variables in Association with aGvHD
Myeloablative conditioning was marginally associated with increased risk of aGvHD (p = 0.054) according to the multivariable logistic regression ( Table 1). None of the other clinical features were found to have a statistically significant association with aGvHD. The statistical significance did not appear for the conditioning in the univariable logistic regression (Supplementary Table S2).

Benchmarking the Taxonomy Profiling Configurations
The abundance of classified reads (mean% ± SD) across all three taxonomic classifiers in species level from lowest to highest were DeepMicrobes (21 ± 11), Kaiju (57 ± 7), and Kraken2 (74 ± 13). The classification performances of the best performing elastic net models based on the reduced species set and the complete species set across profiling configurations are outlined in Figure 2. The best performing models for both complete and reduced species sets were based on the combination of the three taxonomic classifiers, i.e., AM(DeepMicrobes,Kraken2,Kaiju). The highest classification performance was found for the reduced species set with AUROC = 0.672 consisting of 1443 species (57.6% of total number of species from the database). We achieved the highest performance, AUROC = 0.652 for AM(DeepMicrobes,Kraken2,Kaiju), across the seven taxonomy profiling configurations when using the complete set of species. The reduced species set had greater AUC of the precision recall curves compared with the complete set. Eubacterium biforme species was found to have the strongest association with aGvHD compared to other species (see Supplementary Figure S1). On the other hand, a species from Eggerthella genus (UMGS854) exhibited the greatest median weight towards NaGvHD in the elastic net model. The results in the following sections are based on the reduced set from AM(DeepMicrobes,Kraken2,Kaiju) as it was found to have the highest AUROC.

Species Diversity in Association with aGvHD
Alpha and beta diversity indices for the aGvHD and NaGvHD groups are illustrated in Figure 3. No significant differences were found between groups in any of the alpha diversity indices, according to Wilcoxon rank-sum tests. The most distinguishing alpha diversity albeit not statistically significant was k-mer diversity (F 0 ) in which greater diversity in median value for NaGvHD group exists compared to aGvHD group (approx. 14% lower in aGvHD group than the NaGvHD group by fold change of the mean values). There were statistically significant differences between groups in the beta diversity indices as outlined in 10 beta diversity indices where beta diversity was greater for NaGvHD group as compared with the aGvHD group. The greatest difference in beta diversity indices was found by Jaccard's nestedness in which it was 37% lower in the aGvHD group by mean fold change than the NaGvHD group. By simulating group-equalized tests of significant difference between beta diversity measures of the groups, we found that the statistical significance based on adjusted p-values with greater diversity in the NaGvHD group existed in at least 39.4% of the random subsets for Canberra distance, and the most for 76.9% of the subsets for Chebyshev distance.
Several species were found to associate with aGvHD based on Wilcoxon's tests (Figure 4a and Supplementary Table S3). The indicator species analysis based on r.g, also yielded the identification of 152 species in total in which 41 species were statistically significant indicators of aGvHD and 111 species being NaGvHD indicators (Figure 4b and Supplementary Table S4). Our analysis based on mutual information (Figure 4c) yielded the identification of the species UMGS807 from Clostridiales order to have the highest conditional mutual information with aGvHD as conditioned by the disease group (acute leukemia vs. other hematological malignancies). The species with second highest conditional mutual information with aGvHD was UMGS66 from Mollicutes class, in which it had higher mutual information with aGvHD than acute leukemia. As per Spearman correlation (Figure 4d), the species Eubacterium biforme had the highest positive correlation with aGvHD, and relatively low correlation with the disease group. The correlation coefficient for top 10 negatively correlated species resides in the interval of (−0.23 −0.20). The species UMGS602 from Fusobacterium genus stands out in Figure 4d for being highly correlated to aGvHD and Leukemia. The species of UMGS431 from Eggerthella genus was commonly found in the shortlists of top 20 highly correlated and highly associated species, respectively based on Spearman and conditional mutual information, in which it indicates that the species is more abundant in patients with no prospective aGvHD (NaGvHD).   Using the permutation-based test in species indicator analysis, we found 157 species were significantly associated with aGvHD (Supplementary Table S5). Based on the indicator species analysis using IndVal.g, 65 species were found to be significant indicator of aGvHD (Supplementary Table S6). The test was done on a filtered set of species to include species with the minimum fidelity (sensitivity) and specificity of 0.2 and 0.5, respectively, leading to 457 remaining species.

Relief and Impurity-Based Algorithms for Feature Ranking
The result of running 13 algorithms based on relief and impurity provided ranking for each species. The most relevant species in association with aGvHD are depicted Figure 5. The most relevant species from this analysis is UMGS2051 from genus Prevotella.

Discussion
We investigated potential associations of human gut species abundances in aHSCT recipients before transplantation with subsequent occurrence of aGvHD through MWAS. We further presented candidate taxa and models that yielded the best predictive performance for aGvHD. This is the first study, to our knowledge, that utilizes metagenomics on stool samples from aHSCT patients to test the association of the abundances of gut bacteria prior to transplantation with prospective aGvHD occurrences. The identification of potential biomarkers of aGvHD prior to transplantation may confer significant advantage to clinical planning. We found potential prognostic biomarkers of aGvHD consistently standing out across multiple methods in MWAS.
Based on the predictive performances of the models, the reduced species set based on the ensemble method (configuration) to aggregate the taxonomic profiles seems to improve signal to noise ratio with regard to the prediction of aGvHD. This could partly be due to the different algorithms used to identify and extract species counts in which one method could outperform others depending for example on genomic characteristics of each species such as reference genome lengths. On the basis that the metagenomic samples are informative for the prediction of aGvHD, the results may suggest that the ensemble method would provide robust predictions. Notably, we used the data collected prior to transplantation whereas other studies [4,29,30] involved post transplantation data that is temporally and clinically closer to possible aGvHD incidences. The relatively small sample size and imbalanced classes were perhaps the main bottlenecks to capture stronger predictive signals, although, this does not demolish the relative strength in the power gained from gut flora measures. In a future study focused on machine learning, more complex models could be utilized to explore the possibility to improve current prediction performance and develop a deployable model.
We found conditioning being associated with aGvHD with a marginally significant effect from the multivariable logistic regression model. It may reflect that the conditioning effect on aGvHD depends on the presence of the other covariates. The association is in agreement with previous findings in which myeloablative conditioning has been shown to increase the risk of aGvHD [31,32]. The study included a limited set of variables where for example diet and ethnicity were lacking. Along with the clinical variables, metagenomic data could be obtained in the pre-and post-aHSCT periods to reveal the role of the treatments in associations with microbial compositions, and other clinical variables. For example, in post-aHSCT period, computational components derived from metagenomic data such as metabolic functions have also shown to indicate the effect of conditioning on gut microbiome for aHSCT patients [33].
Our statistical analyses revealed several species potentially associate with aGvHD. In particular, significant differences in multiple beta diversity indices gave support to the hypothesized association of the gut microbes with aGvHD. The highest difference in alpha diversity was achieved based on the k-mer measure which underlines the importance of including reference-less methods to reduce risks of biased inference due to limited list of taxa in reference databases as a common limitation for taxonomy-based methods. Another notable finding regarding species compositions was that the gut microbiome was strongly affected by the underlying disease, that is acute leukemia. This was revealed by PERMANOVA in which the disease group significantly associated with the difference in the species abundances and the effect size was seemingly higher than the association of recipient sex with species abundances. This underlined the importance of the effect of leukemia as the underlying disease to a higher extent and recipient sex to a lesser extent in the assessment of the risk of aGvHD prior to aHSCT based on metagenomics data. The lack of significant association between the bacterial compositions and aGvHD reflected in PERMANOVA and ANOSIM, as observed mainly for Bray-Curtis distance, may not necessarily rule out the presence of individual bacteria in significant association with aGvHD. It might be explained by the sensitivity of PERMANOVA and ANOSIM to unbalanced designs. The oversampling approach as an adjustment for the unbalanced design as well as the use of robust Aitchison distance with demonstrated stability to subsetting and aggregation [34] seemed effective to reveal significant associations. As such, the oversampling results suggested presence of interactions of aGvHD with recipient sex, recipient age, and disease group on independent PERMANOVA tests. The inference on the interactions is constraint to the oversampling with limitations of synthetic data to fully represent the original samples.
There are previous findings in line with the current study in which a subset of 10 bacteria were consistently determined to associate with aGvHD. As such, bacterial species from the genera Prevotella [5,35], Bacteroides [36], Parabacteroides [37] were found to be associated with aGvHD. The studies referred to for the bacteria are not necessarily inclusive taking into account different names available for the same species and may require a systematic review to cover the details. For example, Holdemanella biformis is an alternative name for the species Eubacterium biforme according to national library of medicine (NCBI) that was found in this study to associate with aGvHD and has been studied in relation to antibiotic shifts in aHSCT patients [38]. Some of the species such as Lascolabacillus massiliensis [39] found in this studies to associate with aGvHD are discovered in recent years that may explain why they are not found in the literature frequently. Bacteria from genus Eggerthella has been previously shown to significantly change in abundance from pre-aHSCT to post-aHSCT in a similar population [8], this may underline the significance of the species from this genus found in relation to aGvHD. In addition, higher abundance of bacteria from genus Fusobacterium with its unique pattern found from the correlation analysis may occur to characterize a specific immune profile of a subset of individuals with Leukemia as underlying disease to have higher risk of the occurrence of aGvHD. The notion of statistical significance should be taken cautiously in this context with regard to false discoveries especially for the species that were not found consistently to associate with aGvHD and that calls for a larger sample size and more diverse population for generalization purposes.

Conclusions
In summary, our study provides an MWAS approach to enrich abundance data derived from taxonomic classification of the human gut microbiome to improve pre-transplant risk estimation of developing aGvHD. The results reveal the potential of the methods to process complex metagenomic data to extract relevant information regarding the associations of the gut microbiome and clinical factors in the study of aGvHD from pre-aHSCT perspective. The work provides prognostic measures of aGvHD prior to aHSCT highlighting the possibility to take advantage of this new information to improve treatment and prevention strategies. In this sense, microbial information provides an enticing opportunity to enhance individualized prognostication and stratification of aGvHD risk groups.
Supplementary Materials: The supporting information can be downloaded at: https://www.mdpi. com/article/10.3390/cells11244089/s1, Figure S1: abundances of the top-20 species contributing to the prediction of aGvHD as well as their associations with aGvHD and NaGvHD quantified in terms of generalized fold change, prevalence shift, and feature AUCs based on SIAMCAT package. Species names are formatted as beginning with their closest known taxonomy level assigned to them from s_:species, g_:genus, o_:order, f_:family, and c_:class or none when there is no assigned name for the species followed by the species identifier from the reference database.; Table S1: study population including 172 patients of two groups: allogeneic hematopoietic stem cell transplantation (HSCT) recipients diagnosed with acute-versus host disease (aGvHD) and HSCT recipients without aGvHD (NaGvHD); Table S2: Univariable logistic regression with aGvHD as dependent variable (NaGvHD:0, aGvHD:1) and clinical variables as predictors (features). Base (reference) and alternate levels for categorical features are mentioned. Number of samples were N = 172 (39 aGvHD and 133 NaGvHD); Table S3: list of 143 differentially abundant species according to the Wilcoxon's test. p values are not adjusted for multiple testing. Species are ordered by p values in ascending order; Table S4: list of 152 differentially abundant species according to indicator species analyses based on group-equalized point-biserial correlation coefficients (r.g). Species are ordered by p values in ascending order; Table S5: list of 157 differentially abundant species according to the permutation test against the null hypothesis that the abundance of each species is not higher in one group (aGvHD) than others (NaGvHD). p values are adjusted for multiple testing using Sidak's method. Mean fold change in percentage is also presented as an effect size. Species are ordered by p values in ascending order; Table S6: list of 65 differentially abundant species according to indicator species analyses based on group-equalized indicator value (indval.g). Species are ordered by p values in ascending order; Table S7: R packages used for all the analyses described in the article. The codes were compiled based on R version 4.1.2 in RStudio.  IrsiCaixa is supported by the RED de SIDA RD16/0025/0041 and co-financed by the ISCIII and the European Regional Development Fund (ERDF), "Investing in your future". This research was also supported by CIBER-Consorcio Centro de Investigación Biomédica en Red-(CB 2021), Instituto de Salud Carlos III, Ministerio de Ciencia e Innovación and Unión Europea-NextGenerationEU.

Institutional Review Board Statement:
The study was conducted in accordance with the Declaration of Helsinki, and approved by the Danish Ethics Committee of the Capital Region (Journal nr H-16047481) as well as the Danish Data Protection Agency (RH-2017-67, I-suite 05320).
Informed Consent Statement: Written informed consent was obtained from all subjects prior to participating in the study.

Data Availability Statement:
The datasets generated and analyzed during this study are derived from patients treated in Denmark. The datasets contain sensitive patient data governed by the General Data Protection Regulation (GDPR) and Danish law. Due to Danish legislation (Act No. 502 of 23 May 2018) and approvals granted by the Danish Data Protection Agency (RH-2017-67, I-suite 05320), it is not possible to upload raw data to a publicly available database. However, access to these data can be made available from the corresponding author on a reasonable request, provided a data transfer agreement is entered according to current regulations. All R packages used in this study are listed in Supplementary Table S7.