Specific gut microbiome members are associated with distinct immune markers in pediatric allogeneic hematopoietic stem cell transplantation

Anna Cäcilia Ingham; Katrine Kielsen; Malene Skovsted Cilieborg; Ole Lund; Susan Holmes; Frank M. Aarestrup; Klaus Gottlob Müller; Sünje Johanna Pamp

doi:10.1186/s40168-019-0745-z

Research
Open Access
Published: 13 September 2019

Specific gut microbiome members are associated with distinct immune markers in pediatric allogeneic hematopoietic stem cell transplantation

Microbiome volume 7, Article number: 131 (2019) Cite this article

4260 Accesses
6 Citations
21 Altmetric
Metrics details

Abstract

Background

Increasing evidence reveals the importance of the microbiome in health and disease and inseparable host-microbial dependencies. Host-microbe interactions are highly relevant in patients receiving allogeneic hematopoietic stem cell transplantation (HSCT), i.e., a replacement of the cellular components of the patients’ immune system with that of a foreign donor. HSCT is employed as curative immunotherapy for a number of non-malignant and malignant hematologic conditions, including cancers such as acute lymphoblastic leukemia. The procedure can be accompanied by severe side effects such as infections, acute graft-versus-host disease (aGvHD), and death. Here, we performed a longitudinal analysis of immunological markers, immune reconstitution and gut microbiota composition in relation to clinical outcomes in children undergoing HSCT. Such an analysis could reveal biomarkers, e.g., at the time point prior to HSCT, that in the future could be used to predict which patients are of high risk in relation to side effects and clinical outcomes and guide treatment strategies accordingly.

Results

In two multivariate analyses (sparse partial least squares regression and canonical correspondence analysis), we identified three consistent clusters: (1) high concentrations of the antimicrobial peptide human beta-defensin 2 (hBD2) prior to the transplantation in patients with high abundances of Lactobacillaceae, who later developed moderate or severe aGvHD and exhibited high mortality. (2) Rapid reconstitution of NK and B cells in patients with high abundances of obligate anaerobes such as Ruminococcaceae, who developed no or mild aGvHD and exhibited low mortality. (3) High inflammation, indicated by high levels of C-reactive protein, in patients with high abundances of facultative anaerobic bacteria such as Enterobacteriaceae. Furthermore, we observed that antibiotic treatment influenced the bacterial community state.

Conclusions

We identify multivariate associations between specific microbial taxa, host immune markers, immune cell reconstitution, and clinical outcomes in relation to HSCT. Our findings encourage further investigations into establishing longitudinal surveillance of the intestinal microbiome and relevant immune markers, such as hBD2, in HSCT patients. Profiling of the microbiome may prove useful as a prognostic tool that could help identify patients at risk of poor immune reconstitution and adverse outcomes, such as aGvHD and death, upon HSCT, providing actionable information in guiding precision medicine.

Background

The microbiome has gained increasing attention as a crucial contributor in the course of various diseases, and as target of treatment [1,2,3]. A number of complications in allogeneic hematopoietic stem cell transplantation (HSCT) have recently been associated with the gut microbiota [4, 5]. HSCT is a curative treatment for various hematologic diseases including malignancies, such as acute lymphoblastic leukemia (ALL), as well as non-malignant diseases, such as metabolic disorders and immune deficiency syndromes [6]. One goal of HSCT is achieving a beneficial graft-versus-leukemia (GvL) effect where donor-derived T lymphocytes and natural killer cells target leukemic cells in the recipient [7]. Prior to allogeneic HSCT, the patients undergo a preparative conditioning regimen involving combinations of chemotherapeutic agents and total body irradiation (TBI) [8] to eradicate leukemic cells and induce immunosuppression. Immunosuppressive treatment (both prior to and post-HSCT) in the stem cell recipient prevents a graft-versus-host reaction caused by cytotoxic donor T lymphocytes that attack healthy cells in the recipient [8]. To limit infectious diseases due to immunosuppression, the patients are administered broad-spectrum antibacterial and antifungal compounds. Subsequently, the patients receive a stem cell graft originating from the bone marrow, peripheral blood or umbilical cord blood of a human leukocyte antigen (HLA)-matched sibling donor or an unrelated donor (i.e., allogeneic HSCT) [9, 10].

In patients undergoing HSCT, it has been previously observed that there are associations between the microbiome and clinical outcomes such as acute graft-versus-host disease (aGvHD) [5, 11, 12] and survival [4, 13]. GvHD after HSCT has been related to an expansion of the order Lactobacillales, especially Enterococcus spp. [11] and Lactobacillus spp. [14] and a loss of Clostridiales [14]. However, so far, few studies have monitored the microbiome longitudinally [11, 14, 15]. This indicates the need for more detailed investigations that take temporal monitoring of both the host immune system and the microbiome into account.

Several markers of the host immune system, including inflammatory markers (such as C-reactive protein (CRP) and interleukin 6 (IL-6)) and markers of intestinal toxicity (such as plasma citrulline), have been studied in HSCT [16, 17]. Potential novel markers, such as antimicrobial peptides (AMPs), have also been proposed to be involved in outcomes after HSCT, for example in immunomodulation and regulation of microbial homeostasis [18]. AMPs, especially defensins such as human beta-defensin 2 and 3 (hBD2 and hBD3), have previously been found to play a role in some inflammatory diseases [19, 20]. However, to our knowledge, AMPs have not yet been employed as markers in the context of HSCT. An interesting research question remains: How are known and novel markers within the host immune system associated with each other or with changes in the microbiome? A better understanding about these associations would provide a more holistic insight into the underlying mechanisms affecting clinical outcomes after HSCT.

Another crucial factor impacting on complications and clinical outcomes after HSCT is immune reconstitution. Immune reconstitution involves the essential cellular components of the adaptive immune system, namely T and B cells, as well as key cellular components of the innate immune defense, namely natural killer (NK) cells, monocytes, and neutrophils [21]. A microbial influence on immune cell differentiation has been observed previously, e.g., commensal Clostridiales were found to regulate Treg cell differentiation in the colon [22]. Therefore, an influence of the intestinal microbiome on immune reconstitution following HSCT is likely, but had not been investigated prior to this study.

Here, we monitor both host factors and the intestinal microbiome longitudinally. We included markers of inflammation (CRP and IL-6) and intestinal toxicity (plasma citrulline) as well as the antimicrobial peptides (hBD2 and 3), which in the following sections are collectively referred to as immune markers. To assess the prognostic potential of immune markers associated with gut microbial dynamics for immune reconstitution and clinical outcomes after HSCT (aGvHD and survival) in a holistic way, we implemented multivariate multi-table (also referred to as multi-way) approaches. This facilitated the integration of a variety of different factors that could influence the patients’ convalescence. We reveal distinct clusters of bacteria associated with sets of immune markers and clinical outcomes. Patients with rapid NK and B cell reconstitution that had no or mild aGvHD and low mortality exhibited high abundances of members belonging to the family of Ruminocaccaceae. In contrast, patients with moderate to severe aGvHD and high mortality showed high plasma concentrations of the antimicrobial peptide hBD2 already prior to HSCT. In these patients, we observed increased Lactobacillaceae abundances.

Results

To assess associations between immune markers and gut microbial dynamics in the context of immune reconstitution and clinical outcomes after HSCT, we monitored 37 children over time undergoing allogeneic HSCT (Fig. 1a, Additional file 1: Table S1). We gained insight into the patients’ immune reconstitution by determining T, B, NK, monocyte, and neutrophil cell counts in peripheral blood (Fig. 1a). We measured C-reactive protein (CRP) and plasma interleukin 6 (IL-6) as markers of inflammation, and human beta-defensin 2 and 3 (hBD2 and 3) as markers for potential innate immune activation and systemic infection (Fig. 1a, Additional file 1: Table S1). Plasma citrulline levels were measured as a marker of intestinal toxicity, being produced selectively by functioning enterocytes. To relate immune marker levels to members of the intestinal microbiota in patients undergoing HSCT, we characterized the longitudinal dynamics of the human intestinal microbiome in a subset of 30 patients by utilizing 16S rRNA gene profiling (Fig 1a).

Patient cohort and outcomes

At the time of HSCT, the 37 patients were on average 8.2 years old (age range 1.1–18.0 years). Twenty-five patients (68%) were diagnosed with at least one bacterial infection at median 75 days post-HSCT (range day − 19 to + 668). Twenty-six patients (70%) had no or mild aGvHD (grade 0–I) (Additional file 1: Table S1). Eleven patients (30%) developed moderate to severe aGvHD (grade II, III, or IV) at median 18 days (range day + 9 to + 45) after transplantation (Additional file 1: Table S1). Seven patients (19%) died during the follow-up period at median 266 days post-HSCT (range day + 9 to + 784) (five relapse-related and two treatment-related deaths) (Additional file 1: Table S1). In total, six patients (16%) relapsed, four of which underwent a re-transplantation. All patients received antibiotics pre- and post-transplantation (Additional file 1: Table S1). Prophylactic trimethoprim-sulfamethoxazole was administered to all patients from day − 7 until transplantation. During the period of neutropenia or latest from day − 1, patients received prophylactic intravenous ceftazidime. In case of infections indicated by fever or microbial culture, ceftazidime was substituted by intravenous meropenem, vancomycin, or other antibiotics, according to culture-based results.

Temporal dynamics of immune markers and the intestinal microbiota in HSCT patients

Prior to assessing the interplay between clinical variables (i.e., immune markers, immune reconstitution, clinical outcomes) and the intestinal microbiota, we characterized these components separately. In order to provide an overview of changes in immune markers and immune cell counts after HSCT in our cohort (n = 37), we assessed their temporal patterns (Fig. 1b, and Additional file 2: Figure S1). Of note, these supplemental univariate analyses mainly served the purpose of visually aiding our subsequent multivariate analyses approaches (Additional file 3: Figure S2). We characterized hBD2 for the first time in the context of HSCT by assessing plasma hBD2 concentrations over time from pre-HSCT to month + 3 post-HSCT in patients compared to healthy controls. The hBD2 concentration differed significantly between time points (P < 0.001, Kendall’s W = 0.6). It increased from pre-HSCT to the day of HSCT (P < 0.001), then slightly decreased in week + 1 (P = 0.038) before increasing again in week + 3 (P = 0.006) (Fig. 1b). HBD2 decreased again in month + 2 (P = 0.014) (Fig. 1b). CRP levels differed significantly between time points (P < 0.001, Kendall’s W = 0.33). They were high pre-HSCT and until week + 2, then decreasing significantly in week + 3 (P < 0.001) with the lowest levels in weeks + 4 to + 6 (P < 0.001) (Additional file 2: Figure S1A). Median plasma citrulline levels were significantly different between time points (P < 0.001, Kendall’s W = 0.32). They decreased from pre-HSCT to week + 1 (P < 0.001) and increased again in week + 3 (P < 0.001) (Additional file 2: Figure S1A). B cell counts (Kendall’s W = 0.5) as well as CD4⁺ T cell counts (Kendall’s W = 0.46) increased steadily from month + 1 to month + 6 (P < 0.001) (Additional file 2: Figure S1B).

To gain insight into intestinal microbial dynamics before, at the time of, and after HSCT, we obtained a total of 97 fecal samples from a subcohort of 30 patients. Using 16S rRNA gene sequence analysis, we identified 239 operational taxonomic units (OTUs) (see “Methods” section). Microbial alpha-diversity was lower at all time points post-HSCT compared to pre-HSCT (Fig. 1c). The median inverse Simpson index decreased from 3.27 (range 1.02–7.4) before HSCT to 2.89 (range 1.04–10.77) on the day of transplantation and further to 2.03 (range 1.0–16.51) post-HSCT (median of week + 1 to + 5) (Fig. 1c). When assessing the relative abundances of the eight most abundant taxonomic families over time, we observed a reduction in Lachnospiraceae and Ruminococcaceae abundances immediately after HSCT from 20.1 and 9.6%, respectively before HSCT to 6.1% and 1.7% on average in weeks + 1 to + 5 (Fig. 1d). In contrast, Enterococcaceae abundances rose from 17.4% before transplantation to 30.2% on average in weeks + 1 to + 5 (Fig. 1d).

Associations between immune markers and immune cell reconstitution in HSCT patients

In order to identify patient baseline parameters and clinical outcomes (e.g., aGvHD, relapse, overall survival), as well as immune markers and immune cell types that might be important determinants in HSCT in relation to the microbiome, we performed variable assessment by permutational multivariate analysis of variance (adonis). The variables that were found to be significant (P ≤ 0.05), i.e., those that explained most variation in the microbial community distance matrix, were selected for subsequent analyses (Additional file 3: Figure S2, Additional file 4: Table S2). Of note, the occurrence of relapse as an indication of transplantation outcome was assessed but not found to be significant in adonis, and was therefore not included in follow-up analyses. We then assessed associations of the selected immune markers and immune cells in the data set comprising 37 patients by determining Spearman’s rank correlations. The hBD2 concentrations pre-HSCT, on the day of HSCT and in weeks + 1 and + 2 post-HSCT were positively correlated with each other (ρ = 0.73–1, P < 0.001) (Fig. 2a). NK cell counts in month + 1 exhibited a positive correlation with total B cell counts (ρ = 0.64, P = 0.0046) and mature B cells counts (ρ = 0.62, P = 0.0114) in month + 2. When we related immune cell reconstitution to outcomes, we observed significantly higher NK cell counts and total B cell counts in month + 2 in patients with no or mild aGvHD (grade 0–I) compared with patients with moderate to severe aGvHD (grade II–IV) (Wilcoxon rank sum test, NK cells, P = 0.011; B cells, P < 0.001) (Fig. 2b).

High plasma hBD2 and monocytes prior to HSCT in patients with high Lactobacillaceae

To gain insight into how the selected immune markers and immune cell counts co-vary with gut microbial abundances in patients with distinct outcomes, we implemented two multivariate multi-table approaches for the subcohort (n = 30), namely sparse partial least squares (sPLS) regression and canonical correspondence analysis (CCpnA). The sPLS regression models OTU abundances as predictors and clinical variables as response variables and explains the latter in an asymmetric (i.e., unidirectional) way. In contrast, the CCpnA assesses relationships between parameters of the immune system and microbiota bidirectionally. In the following paragraphs, the results of these two analyses are reported for one observed cluster at a time, respectively.

First, we performed sPLS regression to reveal multivariate correlation structures between immune markers, immune cell counts, and OTU abundances, modeling the latter as explanatory variables. The sPLS regression and subsequent hierarchical clustering suggested that the data separated into three clusters (Fig. 3a). High monocyte counts and high plasma hBD2 concentration prior to HSCT, high patient age at the time of transplantation, and high abundances of Lactobacillaceae independent of time point contributed the most to the formation of cluster 1. Of note, monocyte counts and hBD2 were positively correlated with each other, in agreement with the correlation analysis above (Fig. 2a). The Lactobacillaceae were represented by microaerophilic Lactobacillus sp. OTUs (e.g., AF413523.1, GU451064.1, KF029502.1) (Fig. 3b, Additional file 5: Table S3). These OTUs exhibited high loading weights in sPLS dimension 2 (Fig. 3c), indicating that they contributed strongly to the separation of clusters in dimension 2 (Fig. 3a, b).

Second, we applied CCpnA to model the canonical relationships between OTU abundances and clinical variables through the construction of common “latent” variables. The CCpnA confirmed the separation of the data into three clusters as observed in the sPLS regression (Figs. 3a, 4, and Additional file 6: Figure S3A), including the clustering of OTUs. In addition, the CCpnA facilitated the inclusion of categorical variables, such as the patients’ baseline parameters (e.g., recipient sex, donor type) and clinical endpoints (aGvHD grade, overall survival). Because CCpnA is an unsupervised method and upheld the results of the sPLS regression, it provides confidence in the cluster findings.

Cluster 1 in the CCpnA seemed to include patients who developed moderate to severe aGvHD (grade II–IV) and who died. As suggested by both, sPLS and CCpnA, these patients exhibited high levels of plasma hBD2 before HSCT (and in week + 1 and + 2) and high monocyte counts before HSCT. OTUs within this CCpnA cluster predominantly included members of the family Lactobacillaceae and were most abundant in fecal samples of these patients (Fig. 4 and Additional file 7: Figure S4). OTUs that were assigned to cluster 1 by the sPLS-based hierarchical clustering were congruently associated with the same clinical variables in the CCpnA (Figs. 3b and 4).

Temporal patterns of the Lactobacillaceae-dominated community state type

Upon revealing an association between high plasma hBD2 concentrations and monocyte counts prior to HSCT, moderate to severe aGvHD, and high mortality with high abundances of Lactobacillaceae in multivariate analyses, we assessed in more detail the importance of longitudinal changes in these components. We implemented an additional approach to identify distinct bacterial community patterns by employing partitioning around medoid (PAM) clustering (see “Methods” section). In this analysis, we detected four community state types (CSTs). Dominating taxa, similar to those identified by the sPLS-based hierarchical clustering, were revealed in the CSTs, e.g., Lactobacillaceae members dominated CST 1 (Additional file 5: Table S3 and Additional file 8: Figure S5). We then used this information to examine temporal community state changes in individual patients, i.e., their transitions between CSTs over the course of time (Fig. 5). Based on clinical outcomes, the patients can be divided into four groups: (1) patients who developed no or mild aGvHD (aGvHD grade 0–I) and survived compared to (2) patients who died, and (3) patients who developed moderate to severe aGvHD (grade II–IV) and survived compared to (4) patients who died (Fig. 5). We observed that six out of eight patients (75%) with moderate to severe aGvHD (groups 3 and 4) harbored the Lactobacillaceae-dominated CST 1 at least once during the monitored period (1×: n = 2, 2×: n = 2, 3×: n = 2) (Fig. 5). In comparison, only five out of 22 patients (23%) with aGvHD grade 0–I (groups 1 and 2) carried CST 1 one or two times (Fig. 5). Interestingly, high abundances of Lactobacillaceae in patients with aGvHD grade II–IV (groups 3 and 4) occurred predominantly at late time points (week + 1 and later) (Fig. 5 and Additional file 7: Figure S4).

Temporal association of Lactobacillaceae with aGvHD and immune markers

To relate temporal changes in immune markers to those in the bacterial community composition in patients with aGvHD grade II–IV who died (group 4), we assessed their individual longitudinal profiles (Fig. 6). In all three patients (P24, P26, P30), a large expansion of Lactobacillaceae abundances after the onset of aGvHD was observed. The average relative abundance of Lactobacillaceae after aGvHD onset was 72.92% (range 0.22–97.04%), as compared to before aGvHD (average 9.88%, range 0.25–37.83%) (Fig. 6). Furthermore, bacterial alpha diversity was lower at the time point after aGvHD onset compared to the time point before (Fig. 6). All three patients were treated with antibiotics for different durations between these two time points and prior to the time point before aGvHD onset.

In agreement with the results of the multivariate analyses, two of the patients (P24 and P26) exhibited between 1.95- and 14.56-times higher plasma hBD2 concentrations for all measured time points compared with the average of the whole data set (average plasma hBD2 concentration: 10,983.9 pg/ml, range 0–177,400.28 pg/ml). The plasma hBD2 concentration was the highest at the time of HSCT in patient P26 and before HSCT in patient P24 (Fig. 6). Of note, patient P24 received corticosteroid-based GvHD prophylaxis prior to HSCT and both patients received anti-thymocyte globulin (ATG) as part of their conditioning regimen. These two patients had underlying malignant diseases and their pre-HSCT CRP levels, providing insight into underlying inflammation, were 1.6 times lower and 1.29 times higher, respectively, compared with the average of the whole data set on that time point (average plasma CRP concentration before HSCT: 19.21 mg/L, range 1–60.36 mg/L), respectively. Monocyte counts before HSCT, i.e., recipient-derived cells, were 6.25 times higher in patient P24 compared with the average of the whole data set (average monocyte count before HSCT: 0.45 × 10⁹/L, range 0.08–2.83 × 10⁹/L) and close to average in patient P26. In both patients, monocyte counts were higher before (i.e., recipient-derived cells) than after HSCT (i.e., donor-derived cells) (Fig. 6).

In contrast to the Lactobacillaceae expansion after aGvHD onset in group 4, we observed high Lactobacillaceae abundances already before aGvHD onset in three survivors (P16, P22, and P29) who developed aGvHD grade II–IV (group 3) (Additional file 9: Figure S6). The average relative abundance of Lactobacillaceae before aGvHD onset was 53.18% (range 0.35–98.86%), as compared to after aGvHD (average 18.25%, range 0–54.78%). However, these observations were limited by a small number of patients per outcome group (groups 3 and 4) (Fig. 5) and therefore cannot serve to provide statistical evidence.

High NK and B cells and no or mild aGvHD in patients with high obligate anaerobes such as Ruminococcaceae

The association between NK and total B cell counts after HSCT observed in the Spearman’s correlation analysis (Fig. 2) was supported by the multivariate analyses (Figs. 3 and 4). In the sPLS regression, cluster 2 included high NK cell counts in months + 1 and + 2, as well as high immature, mature, and total B cell counts in month + 2 (Fig. 3a). Furthermore, NK and B cell counts were positively correlated with OTUs found in cluster 2 (Fig. 3a, b), in particular with members of the family Ruminococcacea (Order: Clostridiales) (Fig. 3b, Additional file 5: Table S3). Strong positive correlations between mature and total B cell counts in month + 2 and NK cell counts in month + 1 were observed, for example, with Faecalibacterium sp. (DQ804549.1) (Fig. 3b, c, Additional file 5: Table S3). Additionally, we observed positive associations of these immune cells with Lachnospiraceae (Order: Clostridiales), among those two Blautia spp. (DQ800353.1 and DQ802363.1) (Fig. 3b, c, Additional file 5: Table S3).

In support of the sPLS regression, the CCpnA also revealed high NK and B cell counts in cluster 2 in months + 1 and + 2. Based on the CCpnA, we found that NK and B cell reconstitution was associated with certain disease outcomes. For example, patients who had no or mild aGvHD and survived were predominantly represented in cluster 2. Cluster 2 included OTUs mainly belonging to Ruminococcaceae and Lachnospiraceae, exhibiting their highest abundances in patient samples associated with this cluster. OTUs with the highest scores in dimension 2 predominantly belonged to the family of Ruminococcacea and matched the sPLS-based hierarchical cluster assignment (Fig. 4).

As all patients received antibiotic treatment prior to and post-transplantation, we examined the potential effect of antibiotics on the intestinal bacterial community composition. A trending influence of the vancomycin treatment was revealed by adonis analysis (P = 0.055). Using a CCpnA model that also included information about antibiotic treatments, we found that patients exhibited higher abundances of Ruminococcaceae and Lachnospiraceae and lower abundances of Enterobacteriaceae at time points without the vancomycin treatment compared to with vancomycin (Additional file 6: Figure S3B). The same pattern was observed for treatment with ciprofloxacin (Additional file 6: Figure S3B).

Community state typing also revealed a state type dominated by Ruminococcaceae and Lachnospiraceae family members (CST 2) (Additional file 8: Figure S5). A total of ten out of 22 patients (45%) with no or mild aGvHD were assigned to CST 2 at least at one time point (1×: n = 6, 2×: n = 3, 3×: n = 1) (Fig. 5). Half of the patients (four out of eight) with moderate or severe aGvHD were also assigned to CST 2 once, but only at the time point before transplantation (Fig. 5). That is, a Ruminococcaceae- and Lachnospiraceae-dominated community only persisted in patients with no or mild aGvHD.

Persistence of Enterococcaceae-dominated community state type

A subcluster of cluster 2 comprised two facultative anaerobe Enterococcus spp. (GQ1330038.1 and AJ272200.1), exhibiting positive correlations with high NK and B cell counts (Fig. 3b) and contributing to the separation of the clusters (Fig. 3c). Enterococcaceae was the most abundant family in the overall study population (Fig. 1d), and community state typing revealed a state type predominantly characterized by Enterococcaceae (CST 4, Additional file 8: Figure S5). A total of nine out of 22 patients (41%) with no or mild aGvHD were assigned to CST 4 at least at one time point (1×: n = 2, 2x: n = 2, ≥ 3×: n = 5) and this Enterococcaceae-dominated CST often persisted in the patient over time (Fig. 5). A quarter of the patients (two out of eight, both survivors) with moderate or severe aGvHD were assigned to CST 4 at least at one time point on the day of or post-HSCT (Fig. 5).

High inflammation in patients with high facultative anaerobic bacteria

In multivariate analyses, OTUs belonging to facultative anaerobic Enterobacteriaceae and Staphylococcaceae were characteristic for sPLS cluster 3 (Fig. 3a, b, Additional file 5: Table S3). Cluster 3 was further comprised of high plasma citrulline concentrations pre-HSCT (before conditioning) and week + 1, high monocyte counts in week + 3, high CD4⁺ T cell counts in month + 2, high CRP levels, particularly in weeks + 1, + 5, and + 6 (Fig. 3a, b). The projection of the variables in the sPLS suggested a weak positive correlation between these clinical variables in dimension 1 (Fig. 3a). However, a weak negative association between CRP levels and monocytes (week + 3) as well as citrulline and CD4⁺ T cell counts is indicated in dimension 2 (Fig. 3a), which is in agreement with some results of the Spearman’s rank correlation tests (Fig. 2).

The clinical variables in cluster 3, in particular high CRP levels post-HSCT, exhibited positive correlations with the OTUs predominantly affiliated with Proteobacteria (e.g., Enterobacteriaceae), Bacteroidetes, and Staphylococcus spp. (Fig. 3b, Additional file 5: Table S3). The strongest positive correlations occurred with several facultative anaerobe bacteria, e.g., Enterobacter sp. LCR81 (FJ976590.1), Escherichia coli (FJ950694.1), and Staphylococcus sp. (JF109069.1). The CCpnA supported observations from the sPLS regression regarding cluster 3, but indicated that they were only represented by a few patient samples (Fig. 4 and Additional file 6: Figure S3A). These samples were characterized by high CRP levels, especially in week + 1, + 5, and + 6. High CRP levels at these time points were not associated with aGvHD grade, i.e., they were not higher in patients with either aGvHD grade 0–I, or grade II–IV. OTUs of cluster 3, i.e., members of Enterobacteriaceae and Staphylococcaceae, exhibited their highest abundances in samples of this cluster (Fig. 4 and Additional file 6: Figure S3A). Of note, patients represented in cluster 3 were younger compared to those in the other two clusters (Fig. 4).

Discussion

The human immune system and host-associated microorganisms are closely interlinked and play central roles in health and disease. The underlying components and mechanisms facilitating interactions between the immune system and microorganisms are however not completely understood. Understanding these associations is particularly relevant in patients that receive components of a “foreign” immune system, such as in patients undergoing allogeneic HSCT. Because both the patients’ immune system and microbiome appears severely affected, they potentially jointly impact on clinical outcomes. Here, we perform an integrated analysis of immune markers, immune reconstitution data, clinical outcomes, and microbiota, and we provide evidence for the association between specific microbial taxa, host immune markers, immune cells, and clinical outcomes.

We observed that a predominance of Clostridiales, represented by Ruminococcaceae and Lachnospiraceae, in the intestine did not persist after transplantation in patients who either developed aGvHD grade II–IV, died, or both. Clostridiales are common colonizers of the healthy distal gut [23], and their loss is associated with microbial community disruption, reduced diversity, and GvHD [14]. A low diversity before conditioning and at the time of engraftment is associated with increased mortality [4, 5], in line with our findings.

In addition, in our multivariate analyses, we observed high plasma hBD2 levels before HSCT and in weeks + 1 and + 2 post-HSCT in patients who developed aGvHD grade II–IV, died, or both. The reason for the high hBD2 levels before HSCT in patients with increased mortality and GvHD is unknown at present. One could speculate that it relates to a higher burden of inflammation before and during the first 2 weeks post-transplantation. This was also reflected in higher counts of (recipient-derived) monocytes pre-HSCT that, by secreting hBD2, might contribute to an inflammatory reaction and thereby a potentially higher risk of aGvHD after HSCT. However, hBD2 concentrations did not correlate with CRP levels in any of our multivariate analyses, and accordingly our data do not support a clear association between levels of systemic inflammation and high hBD2 levels. High hBD2 concentrations in weeks + 1 and + 2 could be an indicator for an innate immune reaction involving donor-derived cells, for example, against opportunistic pathogens that may have translocated to the bloodstream [24]. This might have been promoted by the decrease in Ruminococcaceae and Lachnospiraceae abundances in these patients after HSCT, indicating a microbial community disruption. Plasma hBD2 concentrations were characterized for the first time here in HSCT patients, and our findings emphasize the importance for further investigations to exploit the potential of hBD2 as a novel candidate marker for outcomes after HSCT. Importantly, our data suggested that differences in hBD2 secretion levels between patients were highly dependent on both, the microbial community composition and the time point relative to HSCT. We suggest that further investigations take into account variations of microbial community patterns and temporal changes when assessing hBD2 in the HSCT context. Moreover, these findings may be refined with a sampling time point homogeneity that is higher than what we provide in this study.

An interesting observation was that high pre-HSCT hBD2 levels and monocyte counts were associated with high abundance of Lactobacillaceae independent of time point. Probiotic Lactobacillus spp. have previously been shown to enhance hBD2-secretion in immune cells, thereby contributing to the innate immune defense [25, 26]. This mechanism could play a protective role during blood stream infections by opportunistic pathogens in HSCT patients. The increase of Lactobacillaceae abundances, particularly after the onset of aGvHD in patients who died, may be explained by a previously proposed compensatory mechanism to reduce aGvHD severity after onset in mice and humans [14]. For instance, high abundances of Lactobacillaceae could indicate homeostasis in the gut microbiome of children and thereby prevent inflammation caused by opportunistic pathogen expansion [27, 28]. Reduced intestinal inflammation might then benefit the outcome of aGvHD. This might however also lead to a less effective graft-versus-leukemia (GvL) effect as suggested by the finding that all patients, for which we observed moderate to severe aGvHD and subsequent increased Lactobacillaceae post-HSCT, overcame aGvHD, but died following a relapse. Only a few patients in our study represented this combination of outcomes (moderate to severe aGvHD with subsequent increase of Lactobacillaceae and death), therefore no significant conclusion can be drawn at this point and further studies are needed to address these observations in more detail.

We provide a discussion on survival following high Lactobacillaceae abundances prior to the onset of aGvHD in Additional file 10.

In contrast to patients with moderate to severe aGvHD, we observed a higher overall survival in patients with no or mild aGvHD in multivariate analyses. The latter patients also had increased numbers of NK and B cells. In agreement with our study, previous studies have shown slower NK cell reconstitution after HSCT in patients with moderate to severe aGvHD compared to those with no or mild aGvHD [29]. Moreover, our results are in agreement with the previously described association of low NK cell numbers after HSCT and reduced overall survival [30], in line with NK cells’ crucial role in the GvL effect [7]. A poor recovery of B cells has been found to pose an increased risk of late infections [31], which might be an explanation for the association of high B cell counts and lower survival. The lower B cell numbers we observed in patients with more severe aGvHD could partially be a consequence of aGvHD or the treatment with corticosteroids, which is known to reduce the number of B cell precursors [21, 29]. However, due to the design of this study, we are not able to differentiate between the roles of these two potential causes that may also act in synergy. However, based on our findings, the intestinal microbiota could also play a contributing role, because faster NK and B cell reconstitution as well as lower aGvHD severity, and higher overall survival were associated with high abundances of obligate anaerobes belonging to Ruminococcaceae and Lachnospiraceae. Indeed, decreased abundances of these bacterial families have previously been associated with an overall microbial disruption and decreased diversity [5, 15, 32]. In line with our observations, reduced microbial diversity was shown to contribute to lower survival [4] and a reduction of Clostridiales was observed in patients presenting aGvHD [14]. Furthermore, in melanoma patients, a low diversity and decreased Ruminococcaceae abundance were associated with a poor response to immunotherapy [33]. A potential explanation for this association could be that members of the order Clostridiales can downregulate inflammation and might thereby prevent aGvHD. Anti-inflammatory components produced by Clostridiales include, for example, urinary 3-indoxyl sulfate (3-IS) [34] and butyrate [35]. Acute GvHD might also reinforce microbial community disruption, as it is known to be accompanied by a reduction of Paneth cell numbers in the intestine. Paneth cells are secretors of α-defensins, important modulators of gut microbial homeostasis [36].

In contrast to the patient group with high abundances of obligate anaerobic bacteria (e.g., Ruminococcaceae and Lachnospiraceae), patients with high abundances of facultative anaerobic bacteria (e.g., Enterobacteriaceae, Staphylococcus spp., and Streptococcus spp.) showed slow NK and B cell reconstitution. An increase in facultative anaerobic bacteria has previously been observed in the gut microbiome of HSCT recipients before the start of conditioning compared with donors [5]. Here, we additionally observed high levels of C-reactive protein (CRP), indicating high inflammation in these patients. This would be in line with the notion that systemic inflammation during the phase of gastrointestinal toxicity after chemotherapy is partially caused by pathogen associated molecular patters (PAMPs) released by translocating bacteria, inducing the release of inflammatory cytokines through binding to Toll like receptors expressed by immune cells, epithelial cells, and other tissue cells. However, the design of the current study does not allow drawing conclusions regarding causality here. Another possible explanation for this association might involve the shift to microbial growth conditions that favor facultative anaerobic bacteria during intestinal inflammation, such as an increased availability of oxygen caused by inflammatory products [37]. Our multivariate analyses indicated that the patients in which we observed these associations were younger compared to the rest of the cohort. Therefore, one could speculate that an immature intestinal microbiome might exhibit a higher susceptibility to opportunistic growth of facultative anaerobic bacteria. Further investigations will have to elucidate this relation and its potential consequences for adjustments of monitoring and treatment by age.

We provide a discussion on our findings regarding associations of adverse outcomes with Enterococcus compared with previous studies in Additional file 10.

Antimicrobial treatment of patients can significantly affect the gut microbiota, especially in children [38]. Early use of antibiotics in general has been found to reduce Clostridiales in the intestine of HSCT patients [39]. Furthermore, a link of high total amounts of antibiotic in GvHD development has been observed in pediatric stem cell recipients [40]. However, to our knowledge, the effects of specific antibiotics on the intestinal microbiome especially in pediatric HSCT patients have not yet been elucidated in detail. Here, we identified a number of specific antibiotics, including vancomycin and ciprofloxacin, associated with simultaneous reduction of Clostridiales (in particular Ruminococcaceae and Lachnospiraceae), similar to what has previously been observed in adult patients [41]. Most interestingly, we found that treatment with vancomycin and ciprofloxacin was not only associated with reduced Clostridiales, but also with increased abundances of facultative anaerobic bacteria, e.g., Enterobacteriaceae (gamma-Proteobacteria). In this patient cohort, these antibiotics were generally given intravenously. Oral vancomycin-treatment in a cohort of rheumatoid arthritis patients was previously shown to be associated with an expansion of Proteobacteria [42]. The use of a prophylactic ciprofloxacin treatment to prevent chemotherapy- and transplantation-related bloodstream infections is an established method [43], and it was proposed that fluoroquinolones (the antibiotic class comprising ciprofloxacin) can prevent intestinal domination of Proteobacteria in HSCT patients [15]. However, ciprofloxacin treatment in healthy subjects has been associated with decreased microbial diversity and decreased Ruminococcaceae and Lachnospiraceae abundances [44], indicating microbial community disruption. Therefore, our findings further challenge the choice of antibiotics, such as vancomycin and ciprofloxacin, in patients undergoing chemotherapy and HSCT. Interestingly, in the present study, a number of less frequently used antibiotic agents, e.g., ceftazidime (a cephalosporin), showed positive associations with high Clostridiales abundances. In agreement, another cephalosporin (cefepime) has previously been attributed clostridial sparing effects [45]. However, ceftazidime was associated with reduced bacterial alpha diversity to a similar degree as vancomycin and ciprofloxacin in a previous study [41]. Elucidating the effects of antibiotic agents potentially contributing to maintaining gut microbial homeostasis in HSCT therefore require further investigation. Of note, we did not take the mode of application of the antibiotics into account here, which might limit the strength of our conclusion since, e.g., an effect of vancomycin on the microbiota might be greater when applied orally compared to intravenously. It remains to be determined whether certain antibiotics modulate the patients’ microbiome and how this might lead to either positive or adverse clinical outcomes.

Conclusions

Our findings support the increasing evidence of microbial involvement in the context of HSCT in cancer patients. We provide evidence for the association between specific microbial taxa and host immune markers. In particular, we examined the prognostic potential of immune markers and gut microbial community dynamics for immune reconstitution and outcomes after HSCT by revealing multivariate associations. We observed increased human beta-defensin 2 in patients with moderate to severe aGvHD and high mortality. In those patients, NK and B cell reconstitution was slow compared to patients with low mortality. These associations only applied when distinct gut microbial abundance patterns were observed, namely low abundances of Ruminococcaceae and high abundances of Lactobacillaceae. Therefore, hBD2, in connection with longitudinal microbial community pattern surveillance, could be further evaluated as a potential novel candidate marker to identify patients at risk of adverse outcomes (e.g., aGvHD) and slow immune cell reconstitution after HSCT, contributing to improved clinical outcomes. Of note, our cohort comprised a relatively small number of patients with different primary diseases and conditioning regimens, and our findings would therefore benefit from being assessed in larger, more homogenous patient groups. Importantly, microbial abundances also depended on antibiotic treatment. Our findings suggest that certain antimicrobial agents might contribute to a shift from obligate to facultative anaerobes. This highlights the need to assess the usage of specific antibiotics in more detail and to take antibiotic treatment into consideration when describing microbial communities in HSCT recipients.

Methods

Patient recruitment and sample collection

We recruited 37 children (age range 1.1–18.0 years) undergoing their first myeloablative allogeneic hematopoietic stem cell transplantation at Copenhagen University Hospital Rigshospitalet, Denmark, from June 2010 to September 2012. Patients’ clinical characteristics are listed in Additional file 1: Table S1. Further information can also be found in previous studies where the cohort has been examined in relation to other questions [17, 46,47,48]. All patients received pretreatment with a myeloablative conditioning regimen, starting on day − 7 (Additional file 1: Table S1). Four patients were re-transplanted at day + 157, + 518, + 712, and + 1360 after the first transplantation, respectively.

Sampling time points were defined according to the following intervals: pre-HSCT (collected between day − 33 and day − 3), at the time of HSCT, preferably before graft infusion (collected between day − 2 and day + 2) and weekly during the first 6 weeks after transplantation (week + 1: day + 3 to day + 10, week + 2: day + 11 to day + 17, week + 3: day + 18 to day + 24, week + 4: day + 25 to day + 31, week + 5: day + 32 to day + 38, week + 6: day + 39 to day + 45) (Fig. 1a). Broader intervals applied to follow-up time points: month + 1 (between days + 21 and + 44), month + 2 (between days + 45 and + 70), month + 3 (between days + 77 and + 105), month + 6 (between days + 161 and + 197), and 1-year post-transplantation (between days + 346 and + 375).

Infections and antibiotics

Bacterial infections from before transplantation until 1 year post-HSCT were taken into consideration for downstream analysis. For each time point, it was recorded whether any bacterial infection occurred within the respective specified interval or not (1/0). Antibiotic treatment from before HSCT (from day − 90) until month + 2 was taken into consideration. We included only those time points corresponding to the time points of microbiota profiling into downstream analyses.

Analysis of T, B, and NK cells in peripheral blood

T, B, and NK cell counts were determined in months + 1, + 2, + 3, and + 6 post-transplantation. Lymphocyte subsets in peripheral blood were quantified using Trucount Tubes (Becton Dickinson, Albertslund, Denmark) together with the following panel of conjugated monoclonal antibodies and analyzed on a FC500 flow cytometer (Beckman Coulter, Copenhagen, Denmark): CD3-PerCP, CD3-FITC, CD4-FITC, CD8-PE, CD45-PerCP, CD16/56-PE, CD20-FITC, and CD19-PE (Becton Dickinson). CD3⁺ T cells, CD3⁺CD4⁺ T cells, and CD3⁺CD8⁺ T cells were determined. NK cells were differentiated by CD3⁻CD45⁺CD16⁺CD56⁺ phenotype. The following B cell phenotypes were distinguished: total B cells (CD45⁺CD19⁺), mature B cells (CD45⁺CD19⁺CD20⁺), and immature B cells (CD45⁺CD19⁺CD20⁻). Data of these immune cell populations have been published previously in a different context [46].

Analysis of monocytes and neutrophils

Leukocyte numbers and subsets were monitored daily during hospitalization and subsequently every week in the outpatient clinic using flow cytometry (Sysmex XN) or, in case of very low leucocytes, counted by microscopy (CellaVision DM96 microscope). Mean monocyte and neutrophil counts were calculated for further analysis per time point according to the intervals specified above.

Quantification of inflammatory markers

EDTA-anticoagulated and heparinized blood was sampled and then centrifuged within 2 h after collection. The plasma was isolated and cryopreserved in 0.5-ml aliquots at − 80 °C. IL-6 levels on day + 7 were determined in EDTA-anticoagulated plasma using the Human Th1/Th2/Th17 Cytometric Bead Array kit (Becton, Dickinson and Co., Denmark) and a FACSCalibur flowcytometer (Becton, Dickinson and Co), according to the manufacturer’s instructions with a detection limit of 2.5 pg/ml. IL-6 data have been published previously in another context for a larger cohort than the patients included here [17, 47]. CRP levels were measured daily by Modular P Modular (normal range, 0 to 10 mg/L) at the Department of Clinical Biochemistry, Copenhagen University Hospital Rigshospitalet, Denmark. Mean CRP levels were calculated for further analysis per time point according to the intervals specified above. Mean CRP levels pre-HSCT include measurements from day − 7 to day − 3, i.e., the days after the start of conditioning (except for day − 7).

Quantification of citrulline

Plasma citrulline concentrations pre-HSCT (before the start of conditioning (day − 7)) and at days + 7 and + 21 were measured by reverse-phase high-performance liquid chromatography of their phenylisothiocyanate derivatives from heparinized plasma, as described previously [17, 47]. Citrulline levels have previously been described for patients of this cohort in a different context [17, 47].

Enzyme-linked immunosorbent assay of human beta defensins

Human beta defensin 2 and 3 (hBD2 and hBD3) concentrations in heparinized plasma samples of 37 patients at eight time points (pre-HSCT (before start of conditioning, except for four patients sampled at day − 6 or − 5), on the day of transplantation, weeks + 1 to + 4, month + 2 and + 3) and ten healthy controls (sampled once each) were quantified by two-step sandwich enzyme-linked immunosorbent assay (ELISA) following the manufacturer’s instructions (Peprotech Human BD-2 and BD-3, Standard ABTS ELISA Development Kit, cat.no. 900-K172 and 900-K210, respectively). The healthy control group comprised ten young adults (four males, six females) with an average age of 22 (range 19–24 years). Samples of three out of ten healthy individuals were additionally spiked to a peptide concentration of 1000 pg/ml. Two replicates were measured per sample and their mean was used for further analysis. Samples were measured undiluted as well as in 1:4 and 1:16 dilutions to also cover concentrations potentially exceeding the upper detection limit. Samples with very high concentrations were additionally measured in 1:32 and 1:128 dilutions. Detection limits were 16–2000 pg/ml for hBD2 and 31–4000 pg/ml for hBD3. Absorbance was measured on a VICTOR™ X3 Multilabel Plate Reader (Perkin Elmer, Inc., USA) at 405 nm. Wavelength correction at 540 nm was used to prevent optical interference caused by the material of the microtiter plate. Concentrations of hBD3 were mostly below the limit of detection, except for a few exceptionally high measurements (average 1450.69 pg/ml, median 0 pg/ml, range 0–279038.71 pg/ml).

DNA isolation from fecal samples and 16S rRNA gene sequencing

Fecal samples for analysis of the intestinal microbiome were collected from a subset of 30 patients at seven time points: pre-HSCT (five patients were sampled after the start of conditioning (between day − 6 and day − 4)), at the time of HSCT and once weekly during the first five weeks after transplantation. The intestinal microbiome was characterized at 1–2 time points in eight patients (27%), at 3–4 time points in 15 patients (50%), and at 5–6 time points in eight patients (27%) (Additional file 1: Table S1). In patients who underwent re-transplantation, no feces samples collected after the second transplantation were included in this study. In total, 97 fecal samples were obtained.

DNA from fecal samples and a blank control were isolated with the use of the Maxwell 16 Instrument (Promega Corporation) following the manufacturer’s instructions for the low elution volume blood DNA system. Alterations to the protocol included additional lysozyme treatment and bead beating with stainless steel beads for 2 min/20 Hz in a tissue lyser (Qiagen). In each sample including the blank control, the V4–V5 region of the 16S ribosomal RNA gene were amplified in PCR using the following barcoded primers: 519F (5#-CAGCAGCCGCGGTAATAC-3#) and 926R (5#-CCGTCAATTCCTTTGAGTTT-3#). Amplicons were then analyzed for quantity and quality in an Agilent 2100 Bioanalyzer (Agilent Technologies) with the use of an Agilent RNA 1000 Nano Kit. For library preparation, 50 ng of DNA from each sample was pooled with multiplex identifiers for 2-region 454 sequencing on GS FLX Titanium PicoTiterPlates (70675) with the use of a GS FLX Titanium Sequencing Kit XLR70 (Roche Diagnostics). Library construction and 454 pyrosequencing were performed at the National High-Throughput DNA Sequencing Center, University of Copenhagen.

16S rRNA gene sequence pre-processing

Raw 454 sequence reads stored in standard flowgram format (SFF) were extracted, converted to and stored in FASTA format with associated quality files (containing sequence quality scores) using the sffinfo command of the bioinformatics software tool mothur [49]. Trimming according to the clipQualLeft and clipQualRight values provided by the sequence provider was disabled because cut-off values are opaque and not customizable.

Analysis was continued in the Quantitative Insights Into Microbial Ecology (QIIME; version 1.9.0) bioinformatics pipeline [50]. FASTA files were demultiplexed and quality filtered using the script split_libraries.py (Mapping files are available from figshare (https://doi.org/10.6084/m9.figshare.6508250)). As the samples were sequenced bidirectional, each FASTA file was demultiplexed in two steps. Firstly, based on a mapping file containing the 519F primer as the “LinkerPrimerSequence” and the 926R primer as the “ReversePrimer,” both in 5′ to 3′ orientation. Secondly, based on a mapping file containing the 926R primer as the “LinkerPrimerSequence” and the 519F primer as the “ReversePrimer,” again both in 5′ to 3′ orientation.

Reads between 200 and 1000 bp length and a minimum quality score of 25 were retained (default). Sequences with homopolymers longer than 200 bp were removed from the data set. Removal of reverse primer sequences (-z truncate_only option) was disabled during demultiplexing. Subsequently, the demultiplexed FASTA files that were not yet primer-truncated were then used to denoise flowgrams (.sff.txt files also generated by mothur’s sffinfo) with QIIME’s denoise_wrapper.py script. Reverse primer-truncation had not been done yet to ensure compatibility between FASTA and .sff.txt files. The denoised FASTA output files were then inflated, i.e., flowgram similarity between cluster centroids was translated to sequence similarity, to be used for OTU picking. Reverse primers and subsequent sequences in the demultiplexed and denoised FASTA files were then truncated using the truncate_reverse_primer.py script. In the following step, the orientation of the primer-truncated reads that started with the 926R primer as the “LinkerPrimerSequence” was adjusted by reverse complementation (with the script adjust_seq_orientation.py). All trimmed reads were then concatenated to a single file for further analysis.

Chimeras were identified using the script identify_chimeric_seqs.py and method usearch61, which performs both de novo (abundance-based) and reference-based detection (by comparing the dataset to the chimera-free reference database Ribosomal Database Project (RDP; training database version 15)). Only those sequences that were flagged as non-chimeras from both detection methods were retained (option –non_chimeras_retention = intersection). Operational taxonomic unit (OTU) clustering was performed, using the script pick_otus.py (based on the SILVA database, Silva_119_rep_set97). OTU tables in BIOM format were created with make_otu_table.py (and subsequently converted to JSON BIOM format to be compatible with analysis in R [51] with the package phyloseq [52]). The OTU table and the taxonomy table are available from figshare (https://doi.org/10.6084/m9.figshare.6508187).

Statistical analyses

Statistical analyses and creation of graphs were performed with the program R (Version 3.4.0, R Foundation for Statistical Computing, Vienna, Austria) [51]. All R scripts documenting our statistical analyses are available from figshare (https://doi.org/10.6084/m9.figshare.6508238). Sequencing data and all related experimental and clinical data (data sets available from figshare, https://doi.org/10.6084/m9.figshare.6508232) were integrated for analysis with the R package phyloseq and its dependencies [52] (Additional file 11). The resulting phyloseq objects are provided through figshare (https://doi.org/10.6084/m9.figshare.6508235). Plots were generated with the packages ggplot2 [53], plotly [54], and mixOmics [55]. Dose-response analysis of the ELISA data was performed with four-parameter log-logistic models in the R package drc [56].

Alpha diversity (measured by inverse Simpson index), levels of human beta-defensin 2 (hBD2) concentration, citrulline and C-reactive protein (CRP), as well as monocyte counts, NK cell counts, total B cell counts, and CD4⁺ T cell counts at different time points were compared using Friedman tests with Benjamini-Hochberg correction for multiple testing (Additional files 3: Figure S2, Additional files 11 and 12). In addition, Kendall’s coefficient of concordance (Kendall’s W) was calculated on ranked data for each marker (Additional file 13). Like the Friedman test, the test for Kendall’s W allowed the comparison of marker levels between time points. In addition, the coefficient of concordance informs about the level of agreement between patients. Therefore, Kendall’s W can be interpreted as a measure of effect size for the Friedman tests. A Kendall’s W < 0.1 was considered as indicating a small effect, 0.1–0.5 as a moderate effect, and > 0.5 as a strong effect. As an exception, hBD2 in healthy controls was compared with hBD2 in patients at individual time points by using Wilcoxon rank-sum tests, because hBD2 was only measured once in the healthy control individuals and can therefore not be analyzed in a Friedman test designed for repeated measurements. Monocyte counts in Additional file 2: Figure S1B are depicted at more time points than indicated in Fig. 1a because not all time points were included into further analyses. The day of HSCT, week + 1, and week + 2 were excluded as data were missing for ≥ 40% of the patients. Months + 2, + 3, and + 6, as well as 1 year were chosen as representative follow-up time points for further analyses as indicated in Fig. 1a.

A core set of OTUs was obtained by retaining 256 OTUs (out of 756 OTUs) with ≥ 5 reads in ≥ 2 samples using the function kOverA() from R package genefilter [57]. Subsequently, 17 OTUs that were more abundant in the blank control than in the majority of samples were removed as potential contaminants prior to downstream analyses. The resulting count data set of 239 OTUs was transformed for subsequent analyses using the function varianceStabilizingTransformation() in the package DESeq2 [58] (Additional file 11). The function implements a Gamma-Poisson mixture model [59] to account for both library size differences and biological variability.

Median imputations were performed for continuous clinical and immune marker data with less than 20% missing values. Variables with more than 20% missing values were excluded from the analysis (Additional files 13 and 14). Central tendencies of immune cell counts at single time points in relation to clinical outcomes (maximum aGvHD grade 0–I vs. grade II–IV) were assessed in univariate analyses by Wilcoxon rank-sum tests and displayed in boxplots (Additional file 14).

A model selection procedure was implemented to find the relevant variables to be included in subsequent multivariate analyses of how microbiome patterns are associated with clinical outcomes, baseline parameters, and immune parameters during the course of transplantation: a Manhattan distance matrix of the variance-stabilized bacterial community data was calculated using the distance() function in phyloseq [52]. Subsequently, permutational multivariate analysis of variance using distance matrices (adonis) for model selection was performed by applying the adonis2() function in the package vegan [60] (Additional file 3: Figure S2, Additional file 15). Permutation design was set up with respect to repeated measurements within the same patients and the intact chronological order of time points. Besides immune marker levels and immune cell counts (pre- and post-HSCT, i.e., recipient- and donor-derived cells, respectively) at the time points described above, we included clinical outcomes (i.e., overall survival, aGvHD (grade 0–I vs. II–IV), and relapse) after transplantations, antibiotic treatment during the course of transplantation, and clinical patient characteristics in the model to account for possible effects of recipient age at the time of transplantation, recipient sex, donor type (sibling vs. unrelated), malignant vs. benign diagnosis, graft type (stem cell source: bone marrow, umbilical cord blood, or peripheral blood), and application of irradiation therapy (yes/no). Variables that were found to be significant (P ≤ 0.05) in the adonis analysis were included in subsequent multivariate multi-table analyses, i.e., sparse partial least squares (sPLS) regression and canonical correspondence analysis (CCpnA) (Additional file 15). Choosing variables with significant effects in adonis for follow-up statistical testing has been performed previously [61]. Even though validation of the set of selected variables through a data-splitting approach might be preferable, this was not feasible due to our relatively small data set. To account for this and to avoid post-selective inference, we renounce calculation of p values from the two analyses that directly depend on the pre-selection (sPLS and CCpnA).

Correlations among the selected clinical variables were assessed in correlation matrices based on Spearman’s rank correlation tests (Additional file 3: Figure S2, Additional file 14). Matrices were calculated using the rcorr() function of the R package Hmisc [62] and displayed with the package corrplot [63]. P values were calculated with the rcorr.adjust() function with correction for multiple testing (method “Holm”). The Spearman’s rank correlation tests were performed on the set of variables selected from the adonis analysis. However, here we assess correlations among those variables, and not between variables and microbial abundances.

Sparse PLS regression was performed by applying the spls() function in the package mixOmics [55] (Additional file 3: Figure S2, Additional file 15). The sPLS regression allows the integration of the microbial community data matrix and the clinical variable matrix for multiple regression. It is robust enough to handle collinearity and noise in the data and is suitable to model multiple response variables [64]. The number of clinical variables to be kept in the model for each component (keepY) was set to 23, corresponding to the number of variables pre-selected with adonis. We ran the sPLS regression with a range of numbers (20–40) of OTUs to be kept for each component (keepX). As the results were robust to this choice, keepX was set to 30. The number of components to choose was estimated with the perf() function and set to ncomp = 2. The sPLS model was run in regression mode. Thereafter, hierarchical clustering was performed within the mixOmics cim() function based on the sPLS regression model with the clustering method “complete linkage” and the distance method “Pearson’s correlation.” Coefficients of pairwise correlations between OTU abundances and clinical variables were thereby obtained. Furthermore, loading plots were generated with the function plotLoadings() (method = “mean”) to visualize loading vectors of specific OTUs that contribute most to the separation of variables in components 1 and 2.

Canonical (i.e., bidirectional) correspondence analysis (CCpnA), a multivariate constrained ordination method, was performed by using the cca() function in the package vegan [60] (Additional file 3: Figure S2, Additional file 15). In this method, the microbial community data matrix is Chi-square transformed and weighted linear regression on pre-selected constraining variables is performed. The resulting fitted values are used for correspondence analysis by singular value decomposition. CCpnA is a constrained method in the sense that it does not aim at depicting all variation in the data, but only the variation directly explained by the constraints (i.e., the provided set of pre-selected variables). The resulting triplot is not displayed as a square representation, but rather corresponds to the percentage of variance explained by axis 1 and 2, respectively, as previously suggested [65]. In contrast to the sPLS analysis, the CCpnA was performed in canonical mode, i.e., modeling bidirectional relations between OTU abundances and clinical variables. OTUs with a correlation of > 0.2/< − 0.2 in the sPLS analysis were included in the CCpnA model.

As another approach to distinguish between microbial community states of the intestinal microbiome, we assigned samples to community state types (CSTs) by partitioning around medoid (PAM) clustering (function pam() in package cluster [66]) based on a Jensen-Shannon distance of the variance stabilized microbial count data (R code modified after [67]) (Additional file 3: Figure S2, Additional file 16). The number of clusters was determined by gap statistic evaluation and silhouette width quality validation. We further assessed patients’ transitions between CSTs over time. OTUs were assigned to CST-based clusters (Additional file 4: Table S2) based on in which CST they exhibited the highest average abundance over all samples (within each CST). Furthermore, we showed detailed longitudinal profiles of the microbial community on family-level, and selected immune markers for individual patients with aGvHD (Additional file 17).

Availability of data and materials

The 16S rRNA gene sequences are available through the European Nucleotide Archive (ENA) at the European Bioinformatics Institute (EBI) under accession number PRJEB25221. The datasets generated and/or analyzed during the current study as well as the R code used to analyze the data are available from the figshare repository at https://figshare.com/projects/Specific_gut_microbiome_members_are_associated_with_distinct_immune_markers_in_allogeneic_hematopoietic_stem_cell_transplantation/35201.

Abbreviations

3-IS:: 3-indoxyl sulfate
aGvHD:: Acute graft-versus-host disease
ALL:: Acute lymphoblastic leukemia
AML:: Acute myeloid leukemia
AMP:: Antimicrobial peptide
ATG:: Anti-thymocyte globulin
CCpnA:: Canonical correspondence analysis
CRP:: C-reactive protein
CST:: Community state type
ELISA:: Enzyme-linked immunosorbent assay
GvL effect:: Graft-versus-leukemia effect
hBD2/hBD3:: Human beta-defensin 2/3
HLA:: Human leukocyte antigen
HSCT:: Hematopoietic stem cell transplantation
IL-6:: Interleukin-6
NK cell:: Natural killer cell
OTU:: Operational taxonomic unit
P:: Patient
PAM clustering:: Partitioning around medoid clustering
QIIME:: Quantitative Insights into Microbial Ecology
SFF:: Standard flowgram format
sPLS regression:: Sparse partial least squares regression
TBI:: Total body irradiation

References

1.
Maukonen J, Kolho K-L, Paasela M, Honkanen J, Klemetti P, Vaarala O, et al. Altered fecal microbiota in paediatric inflammatory bowel disease. J Crohn’s Colitis. 2015;9:1088–95. https://doi.org/10.1093/ecco-jcc/jjv147.
Article Google Scholar
2.
Lynch SV, Pedersen O. The human intestinal microbiome in health and disease. N Engl J Med. 2016;375:2369–79. https://doi.org/10.1056/NEJMra1600266.
CAS Article PubMed Google Scholar
3.
Lemon KP, Armitage GC, Relman DA, Fischbach MA. Microbiota-targeted therapies: an ecological perspective. Sci Transl Med. 2012;4:137rv5. https://doi.org/10.1126/scitranslmed.3004183.
CAS Article PubMed PubMed Central Google Scholar
4.
Taur Y, Jenq RR, Perales M, Littmann ER, Morjaria S, Ling L, et al. The effects of intestinal tract bacterial diversity on mortality following allogeneic hematopoietic stem cell transplantation. Transplantation. 2014;124:1174–82. https://doi.org/10.1182/blood-2014-02-554725.The.
CAS Article Google Scholar
5.
Liu C, Frank DN, Horch M, Chau S, Ir D, Horch EA, et al. Associations between acute gastrointestinal GvHD and the baseline gut microbiota of allogeneic hematopoietic stem cell transplant recipients and donors. Bone Marrow Transplant Adv online Publ. 2017;doi August:1–8. doi:https://doi.org/10.1038/bmt.2017.200.
CAS Article Google Scholar
6.
Appelbaum FR. Haematopoietic cell transplantation as immunotherapy. Nature. 2001;411:385–9. https://doi.org/10.1038/35077251.
CAS Article PubMed Google Scholar
7.
Barrett AJ. Understanding and harnessing the graft-versus-leukaemia effect. Br J Haematol. 2008;142:877–88. https://doi.org/10.1111/j.1365-2141.2008.07260.x.
CAS Article PubMed Google Scholar
8.
Gyurkocza B, Sandmaier BM. Conditioning regimens for hematopoietic cell transplantation: one size does not fit all. Blood. 2014;124:344–53. https://doi.org/10.1182/blood-2014-02-514778.
CAS Article PubMed PubMed Central Google Scholar
9.
Holtick U, Albrecht M, Chemnitz JM, Theurich S, Skoetz N, Scheid C, et al. Bone marrow versus peripheral blood allogeneic haematopoietic stem cell transplantation for haematological malignancies in adults. Cochrane Database Syst Rev. 2014. https://doi.org/10.1002/14651858.CD010189.pub2.
10.
Tiercy J-M. How to select the best available related or unrelated donor of hematopoietic stem cells? Haematologica. 2016;101:680–7. https://doi.org/10.3324/haematol.2015.141119.
CAS Article PubMed PubMed Central Google Scholar
11.
Holler E, Butzhammer P, Schmid K, Hundsrucker C, Koestler J, Peter K, et al. Metagenomic analysis of the stool microbiome in patients receiving allogeneic stem cell transplantation: loss of diversity is associated with use of systemic antibiotics and more pronounced in gastrointestinal graft-versus-host disease. Biol Blood Marrow Transplant. 2014;20:640–5. https://doi.org/10.1016/j.bbmt.2014.01.030.
Article PubMed PubMed Central Google Scholar
12.
Shono Y, Docampo MD, Peled JU, Perobelli SM, Jenq RR. Intestinal microbiota-related effects on graft-versus-host disease. Int J Hematol. 2015;101:428–37. https://doi.org/10.1007/s12185-015-1781-5.
CAS Article PubMed PubMed Central Google Scholar
13.
Jenq RR, Taur Y, Devlin SM, Ponce DM, Goldberg JD, Ahr KF, et al. Intestinal Blautia is associated with reduced death from graft-versus-host disease. Biol Blood Marrow Transplant. 2015;21:1373–83. https://doi.org/10.1016/j.bbmt.2015.04.016.
Article PubMed PubMed Central Google Scholar
14.
Jenq RR, Ubeda C, Taur Y, Menezes CC, Khanin R, Dudakov JA, et al. Regulation of intestinal inflammation by microbiota following allogeneic bone marrow transplantation. J Exp Med. 2012;209:903–11. https://doi.org/10.1084/jem.20112408.
CAS Article PubMed PubMed Central Google Scholar
15.
Taur Y, Xavier JB, Lipuma L, Ubeda C, Goldberg J, Gobourne a, et al. Intestinal domination and the risk of bacteremia in patients undergoing allogeneic hematopoietic stem cell transplantation. Clin Infect Dis. 2012;55:905–14. https://doi.org/10.1093/cid/cis580.
CAS Article PubMed PubMed Central Google Scholar
16.
Gosselin KB, Feldman HA, Sonis AL, Bechard LJ, Kellogg MD, Gura K, et al. Serum citrulline as a biomarker of gastrointestinal function during hematopoietic cell transplantation in children. J Pediatr Gastroenterol Nutr. 2014;58:709–14. https://doi.org/10.1097/MPG.0000000000000335.
CAS Article PubMed PubMed Central Google Scholar
17.
Pontoppidan PL, Jordan K, Carlsen AL, Uhlving HH, Kielsen K, Christensen M, et al. Associations between gastrointestinal toxicity, micro RNA and cytokine production in patients undergoing myeloablative allogeneic stem cell transplantation. Int Immunopharmacol. 2015;25:180–8. https://doi.org/10.1016/j.intimp.2014.12.038.
CAS Article PubMed Google Scholar
18.
Sporrer D, Gessner A, Hehlgans T, Oefner PJ, Holler E. The microbiome and allogeneic stem cell transplantation. Curr Stem Cell Reports. 2015;1:53–9. https://doi.org/10.1007/s40778-014-0006-9.
CAS Article Google Scholar
19.
Liu S, He LR, Wang W, Wang GH, He ZY. Prognostic value of plasma human beta-defensin 2 level on short-term clinical outcomes in patients with community-acquired pneumonia: a preliminary study. Respir Care. 2013;58:655–61. https://doi.org/10.4187/respcare.01827.
Article PubMed Google Scholar
20.
Meisch JP, Nishimura M, Vogel RM, Sung HC, Bednarchik BA, Ghosh SK, et al. Human β-Defensin 3 peptide is increased and redistributed in Crohn’s ileitis. Inflamm Bowel Dis. 2013;19:942–53. https://doi.org/10.1097/MIB.0b013e318280b11a.
Article PubMed PubMed Central Google Scholar
21.
de Koning C, Plantinga M, Besseling P, Boelens JJ, Nierkens S. Immune reconstitution after allogeneic hematopoietic cell transplantation in children. Biol Blood Marrow Transplant. . https://doi.org/10.1016/j.bbmt.2015.08.028.
CAS Article Google Scholar
22.
Furusawa Y, Obata Y, Hase K. Commensal microbiota regulates T cell fate decision in the gut. Semin Immunopathol. 2015;37:17–25. https://doi.org/10.1007/s00281-014-0455-3.
Article PubMed Google Scholar
23.
Donaldson GP, Lee SM, Mazmanian SK. Gut biogeography of the bacterial microbiota. Nat Rev Microbiol. 2015;14:20–32. https://doi.org/10.1038/nrmicro3552.
CAS Article PubMed PubMed Central Google Scholar
24.
Taur Y, Pamer EG. Microbiome mediation of infections in the cancer setting. Genome Med. 2016;8:40. https://doi.org/10.1186/s13073-016-0306-z.
CAS Article PubMed PubMed Central Google Scholar
25.
Schlee M, Harder J, Köten B, Stange EF, Wehkamp J, Fellermann K. Probiotic lactobacilli and VSL#3 induce enterocyte β-defensin 2. Clin Exp Immunol. 2008;151:528–35. https://doi.org/10.1111/j.1365-2249.2007.03587.x.
CAS Article PubMed PubMed Central Google Scholar
26.
Paolillo R, Romano Carratelli C, Sorrentino S, Mazzola N, Rizzo A. Immunomodulatory effects of Lactobacillus plantarum on human colon cancer cells. Int Immunopharmacol. 2009;9:1265–71. https://doi.org/10.1016/j.intimp.2009.07.008.
CAS Article PubMed Google Scholar
27.
Chung H, Pamp SJ, Hill JA, Surana NK, Edelman SM, Troy EB, et al. Gut immune maturation depends on colonization with a host-specific microbiota. Cell. 2012;149:1578–93. https://doi.org/10.1016/j.cell.2012.04.037.
CAS Article PubMed PubMed Central Google Scholar
28.
Bäumler AJ, Sperandio V. Interactions between the microbiota and pathogenic bacteria in the gut. Nature. 2016;535:85–93. https://doi.org/10.1038/nature18849.
CAS Article PubMed PubMed Central Google Scholar
29.
Fujimaki K, Maruta A, Yoshida M, Kodama F, Matsuzaki M, Fujisawa S, et al. Immune reconstitution assessed during five years after allogeneic bone marrow transplantation. Bone Marrow Transplant. 2001 2712. 2001;27:1275. https://doi.org/10.1038/sj.bmt.1703056.
CAS Article Google Scholar
30.
Bartelink IH, Belitser SV, Knibbe CAJ, Danhof M, de Pagter AJ, Egberts TCG, et al. Immune reconstitution kinetics as an early predictor for mortality using various hematopoietic stem cell sources in children. Biol Blood Marrow Transplant. 2013;19:305–13. https://doi.org/10.1016/j.bbmt.2012.10.010.
Article PubMed Google Scholar
31.
Corre E, Carmagnat M, Busson M, de Latour RP, Robin M, Ribaud P, et al. Long-term immune deficiency after allogeneic stem cell transplantation: B-cell deficiency is associated with late infections. Haematologica. 2010;95:1025–9. https://doi.org/10.3324/haematol.2009.018853.
Article PubMed PubMed Central Google Scholar
32.
Montassier E, Gastinne T, Vangay P, Al-Ghalith GA, Bruley Des Varannes S, Massart S, et al. Chemotherapy-driven dysbiosis in the intestinal microbiome. Aliment Pharmacol Ther. 2015;42:515–28.
CAS Article Google Scholar
33.
Gopalakrishnan V, Spencer CN, Nezi L, Reuben A, Andrews MC, Karpinets TV, et al. Gut microbiome modulates response to anti-PD-1 immunotherapy in melanoma patients. Science. 2018;359:97–103. https://doi.org/10.1126/science.aan4236.
CAS Article PubMed Google Scholar
34.
Weber D, Oefner PJ, Hiergeist A, Koestler J, Gessner A, Weber M, et al. Low urinary indoxyl sulfate levels early after transplantation reflect a disrupted microbiome and are associated with poor outcome. Blood. 2015;126:1723–8. https://doi.org/10.1182/blood-2015-04-638858.
CAS Article PubMed Google Scholar
35.
Mathewson ND, Jenq R, Mathew AV, Koenigsknecht M, Hanash A, Toubai T, et al. Gut microbiome-derived metabolites modulate intestinal epithelial cell damage and mitigate graft-versus-host disease. Nat Immunol. 2016;17:505–13. https://doi.org/10.1038/ni.3400.
CAS Article PubMed PubMed Central Google Scholar
36.
Salzman NH, Hung K, Haribhai D, Chu H, Karlsson-Sjöberg J, Amir E, et al. Enteric defensins are essential regulators of intestinal microbial ecology. Nat Immunol. 2010;11:76–82. https://doi.org/10.1038/ni.1825.
CAS Article PubMed Google Scholar
37.
Rivera-Chávez F, Lopez CA, Bäumler AJ. Oxygen as a driver of gut dysbiosis. Free Radic Biol Med. 2017;105:93–101. https://doi.org/10.1016/J.FREERADBIOMED.2016.09.022.
Article PubMed Google Scholar
38.
Fouhy F, Guinane CM, Hussey S, Wall R, Ryan CA, Dempsey EM, et al. High-throughput sequencing reveals the incomplete, short-term recovery of infant gut microbiota following parenteral antibiotic treatment with ampicillin and gentamicin. Antimicrob Agents Chemother. 2012;56:5811–20. https://doi.org/10.1128/AAC.00789-12.
CAS Article PubMed PubMed Central Google Scholar
39.
Weber D, Jenq RR, Peled JU, Taur Y, Hiergeist A, Koestler J, et al. Microbiota disruption induced by early use of broad spectrum antibiotics is an independent risk factor of outcome after allogeneic stem cell transplantation. Elsevier Inc. 2017. https://doi.org/10.1016/j.bbmt.2017.02.006.
CAS Article Google Scholar
40.
Simms-Waldrip TR, Sunkersett G, Coughlin LA, Savani MR, Arana C, Kim J, 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. https://doi.org/10.1016/J.BBMT.2017.02.004.
CAS Article PubMed Google Scholar
41.
Weber D, Hiergeist A, Weber M, Dettmer K, Wolff D, Hahn J, et al. Detrimental effect of broad-spectrum antibiotics on intestinal microbiome diversity in patients after allogeneic stem cell transplantation: Lack of commensal sparing antibiotics. Clin Infect Dis. 2018. https://doi.org/10.1093/cid/ciy711.
Article Google Scholar
42.
Isaac S, Scher JU, Djukovic A, Jiménez N, Littman DR, Abramson SB, et al. Short- and long-term effects of oral vancomycin on the human intestinal microbiota. J Antimicrob Chemother. 2017;72:128–36. https://doi.org/10.1093/jac/dkw383.
CAS Article PubMed Google Scholar
43.
Yeh T-C, Liu H-C, Hou J-Y, Chen K-H, Huang T-H, Chang C-Y, et al. Severe infections in children with acute leukemia undergoing intensive chemotherapy can successfully be prevented by ciprofloxacin, voriconazole, or micafungin prophylaxis. Cancer. 2014;120:1255–62. https://doi.org/10.1002/cncr.28524.
CAS Article PubMed Google Scholar
44.
Dethlefsen L, Relman DA. Incomplete recovery and individualized responses of the human distal gut microbiota to repeated antibiotic perturbation. Proc Natl Acad Sci U S A. 2011;108 Suppl Supplement 1:4554–4561. doi:https://doi.org/10.1073/pnas.1000087107.
Article Google Scholar
45.
Shono Y, Docampo MD, Peled JU, Perobelli SM, Velardi E, Tsai JJ, et al. Increased GVHD-related mortality with broad-spectrum antibiotic use after allogeneic hematopoietic stem cell transplantation in human patients and mice. Sci Transl Med. 2016;8:339ra71. https://doi.org/10.1126/scitranslmed.aaf2311.
CAS Article PubMed PubMed Central Google Scholar
46.
Kielsen K, Jordan KK, Uhlving HH, Pontoppidan PL, Shamim Z, Ifversen M, et al. T cell reconstitution in allogeneic haematopoietic stem cell transplantation: Prognostic significance of plasma interleukin-7. Scand J Immunol. 2015;81:72–80. https://doi.org/10.1111/sji.12244.
CAS Article PubMed Google Scholar
47.
Jordan K, Pontoppidan P, Uhlving HH, Kielsen K, Burrin DG, Weischendorff S, et al. Gastrointestinal toxicity, systemic inflammation, and liver biochemistry in allogeneic hematopoietic stem cell transplantation. Biol Blood Marrow Transplant. 2017;23:1170–6. https://doi.org/10.1016/j.bbmt.2017.03.021.
CAS Article PubMed Google Scholar
48.
Kielsen K, Ryder LP, Lennox-Hvenekilde D, Gad M, Nielsen CH, Heilmann C, et al. Reconstitution of Th17, Tc17 and Treg cells after paediatric haematopoietic stem cell transplantation: Impact of interleukin-7. Immunobiology. 2018;223:220–6. https://doi.org/10.1016/j.imbio.2017.10.023.
CAS Article PubMed Google Scholar
49.
Schloss PD, Westcott SL, Ryabin T, Hall JR, Hartmann M, Hollister EB, et al. Introducing mothur: open-source, platform-independent, community-supported software for describing and comparing microbial communities. Appl Environ Microbiol. 2009;75:7537–41. https://doi.org/10.1128/AEM.01541-09.
CAS Article PubMed PubMed Central Google Scholar
50.
Caporaso JG, Kuczynski J, Stombaugh J, Bittinger K, Bushman FD, Costello EK, et al. QIIME allows analysis of high-throughput community sequencing data. Nat Methods. 2010;7:335–6. https://doi.org/10.1038/nmeth.f.303.
CAS Article PubMed PubMed Central Google Scholar
51.
R Core Team. R: A language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing; 2017. https://www.r-project.org/
Google Scholar
52.
McMurdie PJ, Holmes S. phyloseq: An R package for reproducible interactive analysis and graphics of microbiome census data. PLoS One. 2013;8:e61217. https://doi.org/10.1371/journal.pone.0061217.
CAS Article PubMed PubMed Central Google Scholar
53.
Wickham H. ggplot2: Elegant graphics for data analysis. New York: Springer Verlag; 2016.
Google Scholar
54.
Sievert C, Parmer C, Hocking T, Chamberlain S, Ram K, Corvellec, Marianne Despouy P. plotly: Create interactive web graphics via “plotly. js.” 2016. https://plot.ly/r, https://cpsievert.github.io/plotly_book/, https://github.com/ropensci/plotly#readme.
55.
Cao K Le, Rohart F, Gonzalez I, Dejean S, Gautier B, Bartolo F, et al. mixOmics: Omics data integration project. R package version 6.1.3. 2017. https://cran.r-project.org/package=mixOmics.
56.
Ritz C, Baty F, Streibig JC, Gerhard D. Dose-response analysis using R. PLoS One. 2015;10:e0146021. https://doi.org/10.1371/journal.pone.0146021.
CAS Article PubMed PubMed Central Google Scholar
57.
Gentleman R, Carey V, Huber W, Genefilter HF. Methods for filtering genes from microarray experiments. R package version. 2017;1:58.1.
Google Scholar
58.
Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:550. https://doi.org/10.1186/s13059-014-0550-8.
Article PubMed PubMed Central Google Scholar
59.
McMurdie PJ, Holmes S. Waste not, want not: why rarefying microbiome data is inadmissible. PLoS Comput Biol. 2014;10:e1003531. https://doi.org/10.1371/journal.pcbi.1003531.
CAS Article PubMed PubMed Central Google Scholar
60.
Oksanen J, Blanchet FG, Friendly M, Kindt R, Legendre P, McGlinn D, et al. Vegan: Community ecology package. R Package Version. 2.4-3. https://cran.R-project.org/package=vegan. 2017. doi:https://doi.org/10.1029/2006JF000545.
61.
Pérez-Losada M, Freishtat RJ, Kwak C, Authelet KJ, Hoptay CE, Crandall KA. Pediatric asthma comprises different phenotypic clusters with unique nasal microbiotas. Microbiome. 2018;6:1–13.
Article Google Scholar
62.
Harrell Jr FE, Dupont C. Hmisc: Harrell miscellaneous. R package version 4.0-3. 2017.
63.
Wei T, Simko V. corrplot: Visualization of a Correlation Matrix. R package version 0.77. 2017. https://github.com/taiyun/corrplot.
64.
Lee D, Lee W, Lee Y, Pawitan Y. Sparse partial least-squares regression and its applications to high-throughput data analysis. Chemom Intell Lab Syst. 2011;109:1–8. https://doi.org/10.1016/J.CHEMOLAB.2011.07.002.
CAS Article Google Scholar
65.
Callahan BJ, Sankaran K, Fukuyama JA, McMurdie PJ, Holmes SP. Bioconductor workflow for microbiome data analysis: from raw reads to community analyses. F1000Research. 2016;5:1492. https://doi.org/10.12688/f1000research.8986.1.
Article PubMed PubMed Central Google Scholar
66.
Maechler M, Rousseeuw P, Struyf A, Hubert M, Hornik K. cluster: Cluster analysis basics and extensions. R package version 2.0.6. 2017.
67.
DiGiulio DB, Callahan BJ, McMurdie PJ, Costello EK, Lyell DJ, Robaczewska A, et al. Temporal and spatial variation of the human microbiota during pregnancy. Proc Natl Acad Sci. 2015;112:11060–5. https://doi.org/10.1073/pnas.1502875112.
CAS Article PubMed Google Scholar

Download references

Acknowledgements

We thank the patients and their families for participating in the study, and the National High-Throughput DNA Sequencing Centre at the University of Copenhagen for library construction and sequencing. Sequence pre-processing described in this study was performed using the DeiC National Life Science Supercomputer at DTU. Furthermore, we are grateful for helpful discussions about the statistical analyses with Pratheepa Jeganathan and Kris Sankaran from the Department of Statistics at Stanford University. We thank Jeffrey Edward Skiby for reviewing the manuscript.

Funding

This work was supported by the European Union’s Framework program for Research and Innovation, Horizon2020 (643476), and by the National Food Institute, Technical University of Denmark.

Author information

Affiliations

Research Group for Genomic Epidemiology, Technical University of Denmark, Kongens Lyngby, Denmark
Anna Cäcilia Ingham, Frank M. Aarestrup & Sünje Johanna Pamp
Institute for Inflammation Research, Department of Rheumatology and Spine Disease, Copenhagen University Hospital, Rigshospitalet, Copenhagen, Denmark
Katrine Kielsen & Klaus Gottlob Müller
Department of Pediatrics and Adolescent Medicine, Copenhagen University Hospital Rigshospitalet, Copenhagen, Denmark
Katrine Kielsen & Klaus Gottlob Müller
Comparative Pediatrics and Nutrition, Department of Clinical Veterinary and Animal Science, University of Copenhagen, Frederiksberg, Denmark
Malene Skovsted Cilieborg
Department of Bio and Health Informatics, Technical University of Denmark, Kongens Lyngby, Denmark
Ole Lund
Department of Statistics, Stanford University, Stanford, USA
Susan Holmes
Department of Bacteria, Parasites and Fungi, Statens Serum Institut, Copenhagen, Denmark
Anna Cäcilia Ingham

Authors

Anna Cäcilia Ingham
View author publications
You can also search for this author in PubMed Google Scholar
Katrine Kielsen
View author publications
You can also search for this author in PubMed Google Scholar
Malene Skovsted Cilieborg
View author publications
You can also search for this author in PubMed Google Scholar
Ole Lund
View author publications
You can also search for this author in PubMed Google Scholar
Susan Holmes
View author publications
You can also search for this author in PubMed Google Scholar
Frank M. Aarestrup
View author publications
You can also search for this author in PubMed Google Scholar
Klaus Gottlob Müller
View author publications
You can also search for this author in PubMed Google Scholar
Sünje Johanna Pamp
View author publications
You can also search for this author in PubMed Google Scholar

Contributions

ACI, KK, KGM, and SJP designed the research; ACI, KK, MSC, and SJP performed the research; ACI, SH, and SJP contributed analytic tools; ACI and SJP analyzed the data; ACI and SJP wrote the manuscript; and KK, SH, FMA, OL, and KGM edited the manuscript. All authors have read and approved the manuscript as submitted.

Corresponding author

Correspondence to Sünje Johanna Pamp.

Ethics declarations

Ethics approval and consent to participate

Written informed consent was obtained from the patients and/or their legal guardians. The study protocol was approved by the local ethics committee (H-1-2010-009) and the Danish Data Protection Agency.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

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

Additional files

Additional file 1:

Table S1. Clinical patient characteristics. General patient characteristics, conditioning regimens, complications, and outcomes for the pediatric cohort (n = 37) and the subcohort (n = 30) for which the intestinal microbiome was characterized. Abbreviations: HLA, human leukocyte antigen; TBI, total body irradiation; CY, Cyclophosphamide; VP16, Etoposide; BU, Busulfan; MEL, Melphalan; GvHD, graft-versus-host disease. (PDF 379 kb)

Additional file 2:

Figure S1. Temporal patterns of immune markers and immune cells in HSCT patients. (A) C-reactive protein (CRP) and plasma citrulline levels in HSCT patients (n = 37) over time. CRP levels were significantly higher prior to HSCT and until week +2 compared to all following time points, e.g. at the day of HSCT (median: 16.93 mg/L, range: 1.22 - 85.28 mg/L) compared to week +3 (median: 3.92 mg/L, range: 1.22 - 55.89 mg/L) (P < 0.001). Plasma citrulline levels were significantly lower in week +1 compared to pre-HSCT (P < 0.001) and week +3 (P < 0.001). (B) Immune cell counts in HSCT patients over time. Monocyte counts are depicted at more time points than indicated in Figure 1a, because not all time points were included into further analyses (see Methods). NK cell counts were higher in months +2 to +6 compared to in month +1 (P < 0.001). B cell counts as well as CD4+ T cell counts increased steadily from month +1 to month +6 (P < 0.001). Y -axes in all plots, except for citrulline, were log10-transformed for better visualization. Zeros were replaced with 1 to avoid undefined values on the log-transformed axes. Asterisks indicate whether the component at each respective time point was significantly different from any of the other time points (showing the maximum significance level). * P < 0.05, ** P < 0.01 and *** P < 0.001. (PDF 419 kb)

Additional file 3:

Figure S2. Workflow of the statistical analysis approach. The diagram displays the major steps of the statistical analyses and their dependencies. Multivariate analyses (blue box) constitute the main approach, especially the multi-table analyses and clustering analyses (green box). To unravel the complexity of the multivariate analyses, these were supplemented with univariate analyses (upper grey box). (PDF 911 kb)

Additional file 4:

Table S2. Results of Permutational Multivariate Analysis of Variance Using Distance Matrices (adonis). Adonis was employed for model selection to identify relevant immune markers and immune cell types to be included in downstream analyses (See Methods for details). Significant variables (P<0.05) are marked in bold. Abbreviations: hBD2_sim, plasma human beta-defensin 2 levels at time points simultaneous to microbiome characterization; CRP_sim, C-reactive protein levels at time points simultaneous to microbiome characterization; Lymphocyte_count_sim, total lymphocyte counts at time points simultaneous to microbiome characterization; pIL6, plasma interleukin 6 concentration; Citr, plasma citrulline concentration; CD3+, CD3+ T cell counts; CD4+, CD3+CD4+ T cell counts; CD8+, CD3+CD8+ T cell counts; B, total B cell (CD45+CD19+) counts; mat_B, mature B cell (CD45+CD19+CD20+) counts. immat_B, immature B cell (CD45+CD19+CD20-) counts; NK, natural killer cell counts; mean_mono, mean monocyte counts at indicated time point; mean_neutro, mean neutrophil counts at indicated time point; Timepoints: pre, prior to transplantation; w0, on the day of transplantation; w1, w2, w3, w4, w5: one, two, three, four and five weeks after transplantation, respectively; m1, m2, m3, m4, m6: one, two, three, four and six months after transplantation, respectively; 1y, 1 year post-transplantation. (PDF 461 kb)

Additional file 5:

Table S3. Taxonomy and cluster affiliation of OTUs strongly associated with host-related variables based on sPLS analysis and community state typing (CST). List of the 57 OTUs correlated strongest with variables in the sPLS analysis (>0.2/<-0.2) . SPLS-based clusters were determined by applying the mixOmics cim() function to the sPLS regression model (hierarchical clustering method: complete linkage, distance method: Pearson’s correlation) (see Methods). Four community state types (CSTs) were defined by clustering of fecal samples with similar microbial community compositions by partitioning around medoid (PAM) clustering (see Methods). OTUs were then assigned to the CST-based clusters in which they exhibited the highest average abundance over all samples. The same taxonomic families dominated in sPLS- and CST-based clusters, respectively. Cluster 1 was dominated by Lactobacillaceae. Cluster 2 was characterized mainly by Ruminococcaceae and Lachnospiraceae. Cluster 3 harbored Proteobacteria (P), e. g. Enterobacteriaceae. CS-typing revealed one additional cluster (4), characterized by a high abundance of Enterococcaceae and Staphylococcaceae. OTU numbers refer to the SILVA database (silva_119_rep_set97). Phyla abbreviations: F, Firmicutes; B, Bacteroidetes; A, Actinobacteria; P, Proteobacteria; FU, Fusobacteria. (PDF 505 kb)

Additional file 6:

Figure S3. Canonical correspondence analysis (CCpnA) of immune markers and intestinal bacterial taxa in patients undergoing HSCT. Triplots showing dimension 1 and 2 of the CCpnA that includes continuous clinical variables (arrows), categorical variables (+), and OTUs (circles). Samples are depicted as triangles. OTUs with a correlation of >0.2/<-0.2 in the sPLS analysis were included in the CCpnA model. Only the variables and OTUs with a score >0.2/<-0.2 in at least one CCpnA dimension are shown. The OTUs in the CCpnA plots are colored according to the cluster they were affiliated with in the sPLS-based hierarchical clustering analysis, and the ellipses present an 80% confidence interval, assuming normal distribution. (A) Full size visualization corresponding to the CCpnA model shown in Figure 4. Plot dimensions correspond to the explained variances of each component. (B) CCpnA including antibiotic treatment at time points simultaneous to microbiome characterization. Antibiotics were added as categorical variables. Depiction of the antibiotic’s name (in red) indicates administration of the particular antibiotic, and the extension “_0” indicates no administration of the respective antibiotic. Abbreviations of variables are the same as in Figure 2. Further abbreviations: graft_BM: stem cell source bone marrow; graft_UC: stem cell source umbilical cord blood. (PDF 1356 kb)

Additional file 7:

Figure S4. Clustered image map (CIM) of OTU abundances by patient in the first two sPLS dimensions. Hierarchical clustering of OTU abundances (bottom) and patients’ fecal samples (right) (clustering method: complete linkage, distance method: Pearson’s correlation) was performed within the mixOmics cim() function based on the sPLS regression model. High abundance of an OTU in a sample is represented as positive correlation in the map (red) and low abundance as negative correlation (blue). The sampling time points of the fecal samples are displayed in the side bar on the left (blue gradient from pre-HSCT time point (light blue) to week +5 post-HSCT (dark blue)). The top side bar shows taxonomic information on family level. Sample names on the right indicate patient (P) and time point (pre: pre-HSCT, d0: day of HSCT, w: week). An “a” or “b” indicates that two samples were collected from the respective patient at the same time point, but on two different days. (PDF 491 kb)

Additional file 8:

Figure S5. Community state types and gut microbial patterns. Heat map of variance stabilized counts of the 50 most abundant OTUs of the intestinal microbiome over all samples, grouped into community state types (CSTs). Based on their OTU-composition, samples were assigned to community state types (CSTs) by partitioning around medoid (PAM) clustering using Bray-Curtis distance. The optimal number of clusters (k = 4) was estimated from the gap statistic and Silhouette width validation. Members of the Lactobacillaceae family dominated the abundance profiles within CST 1. CST 2 exhibited domination by Lachnospiraceae, Erysipelotrichaceae and Ruminococcaceae members. Enterobacteriaceae, Streptococcaceae and Staphylococcacea were characteristic for CST 3. CST 4 was characterized by a high abundance of Enterococcaceae. Average Silhouette width was s(i) = 0.16 (range: -0.02 – 0.36), with CST 1 and CST 4 being the best defined clusters (s(i) = 0.23 and 0.36, respectively). A Silhouette coefficient s(i) close to 1 indicates appropriate clustering of the respective samples. Sample names at the bottom indicate patient (P) and time point (pre: pre-HSCT, d0: day of HSCT, w: week). An “a” or “b” indicates that two samples were collected from the respective patient at the same time point, but on two different days. (PDF 516 kb)

Additional file 9:

Figure S6. Longitudinal profiles of microbial community composition and immune markers in patients with aGvHD who survived. In three representative patients with moderate to severe aGvHD who survived, high abundance of Lactobacillaceae was observed already before aGvHD onset. None of the depicted patients had a bacterial infection recorded during the monitored period. InvSimpson, inverse Simpson diversity index; hBD2, human beta-defensin 2. (PDF 358 kb)

Additional file 10:

Supplementary Discussion. Discussion concerning survival following high Lactobacillaceae abundances prior to the onset of aGvHD, and associations of adverse outcomes with Enterococcus compared with previous studies. (PDF 732 kb)

Additional file 11:

R data analysis report 1. Data preparation, filtering and transformation. (HTML 2337 kb)

Additional file 12:

R data analysis report 2. Bacterial alpha-diversity over time and rank abundance curve for the gut microbiome of HSCT patients. (HTML 1461 kb)

Additional file 13:

R data analysis report 3. Temporal patterns of immune markers and immune cells in HSCT patients. (HTML 1457 kb)

Additional file 14:

R data analysis report 4. Correlations between immune markers, immune cell counts, and outcomes in HSCT patients. (HTML 1238 kb)

Additional file 15:

R data analysis report 5. Variable selection and multivariate analyses of immune parameters and intestinal bacterial taxa in HSCT patients. (HTML 2214 kb)

Additional file 16:

R data analysis report 6. Clustering of samples into Community State Types (CSTs) based on Jenson-Shannon divergence. (HTML 2338 kb)

Additional file 17:

R data analysis report 7. Longitudinal profiles of microbial community composition and immune markers. (HTML 4520 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Cite this article

Ingham, A.C., Kielsen, K., Cilieborg, M.S. et al. Specific gut microbiome members are associated with distinct immune markers in pediatric allogeneic hematopoietic stem cell transplantation. Microbiome 7, 131 (2019). https://doi.org/10.1186/s40168-019-0745-z

Download citation

Received: 27 March 2019
Accepted: 29 August 2019
Published: 13 September 2019
DOI: https://doi.org/10.1186/s40168-019-0745-z

Keywords

Gut microbiota
Pediatric cancer
HSCT
16S rRNA gene profiling
Data integration
Immune reconstitution
Ruminococcaceae
Human beta-defensin 2
acute GvHD
B cells and NK cells

Specific gut microbiome members are associated with distinct immune markers in pediatric allogeneic hematopoietic stem cell transplantation

Abstract

Background

Results

Conclusions

Background

Results

Patient cohort and outcomes

Temporal dynamics of immune markers and the intestinal microbiota in HSCT patients

Associations between immune markers and immune cell reconstitution in HSCT patients

High plasma hBD2 and monocytes prior to HSCT in patients with high Lactobacillaceae

Temporal patterns of the Lactobacillaceae-dominated community state type

Temporal association of Lactobacillaceae with aGvHD and immune markers

High NK and B cells and no or mild aGvHD in patients with high obligate anaerobes such as Ruminococcaceae

Persistence of Enterococcaceae-dominated community state type

High inflammation in patients with high facultative anaerobic bacteria

Discussion

Conclusions

Methods

Patient recruitment and sample collection

Infections and antibiotics

Analysis of T, B, and NK cells in peripheral blood

Analysis of monocytes and neutrophils

Quantification of inflammatory markers

Quantification of citrulline

Enzyme-linked immunosorbent assay of human beta defensins

DNA isolation from fecal samples and 16S rRNA gene sequencing

16S rRNA gene sequence pre-processing

Statistical analyses

Availability of data and materials

Abbreviations

References

Acknowledgements

Funding

Author information

Affiliations

Contributions

Corresponding author

Ethics declarations

Ethics approval and consent to participate

Consent for publication

Competing interests

Additional information

Publisher’s Note

Additional files

Rights and permissions

About this article

Cite this article

Keywords

Microbiome

Contact us