Lasting Gammaproteobacteria profile changes characterized hematological cancer patients who developed oral mucositis following conditioning therapy

ABSTRACT Background Oral mucositis (OM) is a common side effect of conditioning therapy implemented before hematopoietic stem cell transplantation (HSCT). The role of oral microbiome in OM is not fully elucidated. Objective To determine oral microbiome profile changes post-conditioning in HSCT patients who developed moderate OM, or mild to no OM. Design Patient groups were: Muc0-1 with OM-score = 0–1 (43 paired samples) and Muc2 with WHO OM-score = 2 (36 paired samples). Bacterial DNA was isolated from oral samples (saliva, swabs of buccal mucosa, tongue, and supragingival plaque) at pre-conditioning (T0), post-conditioning mucositis onset (TMuc), and one-year post-conditioning (TYear). 16S-rRNA gene next-generation sequencing was used to determine the relative abundance (RA) of >700 oral species. Alpha-diversity, beta-diversity and linear discriminant analyses (LDA) were performed Muc2 versus Muc0-1. Results Muc2 oral microbiome alpha- and beta-diversity differed between T0 and TMuc. Muc2 alpha-diversity and Muc0-1 beta-diversity did not differ between T0 and TYear. T0 to TMuc LDA scores were significant in Muc2 for Gammaproteobacteria. For Muc2 patients, the average RA decreased for Haemophilus parainfluenza, a species known as mucosal surfaces protector, but increased for Escherichia-Shigella genera. Conclusions Post-conditioning OM might contribute to long-term oral microbiome changes affecting Gammaproteobacteria, in HSCT patients.

Patients with hematological cancers undergo conditioning therapy prior to hematopoietic stem cell transplantation. A myelosuppressive or myeloablative conditioning regimen provides the anti-cancer effect, while the transplant re-establishes hematopoietic functioning [1,2]. Over 50,000 patients undergo conditioning therapy prior to hematopoietic stem cell transplantation each year [2].
Oral mucositis (OM) is one of the most frequently occurring side effects associated with conditioning therapy. Patients with OM experience damage to the oral mucosa ranging from redness and soreness to ulceration [3,4]. OM grading defined by the World Health Organization (WHO) utilizes a 0 to 4 scale, based on the presence of erythema and soreness (score 1 and above) and/or ulcers (score 2 and above) and the ability for patients to consume solid food (score 2) or liquids only (score 3), or inability to consume any food orally (score 4), due to the painful condition of inflamed oral mucosa [5]. Sonis has developed a five-phase model describing OM pathophysiology [4].
Patients consistently report OM as the most painful and debilitating side effect of cancer treatment [6][7][8].
The incremental cost of OM-associated hospitalization can be as high as 70,000 USD for patients who develop ulcerative mucositis post-conditioning [9]. The substantial impact of OM, combined with a lack of evidence-based treatment protocols, creates a knowledge gap impacting patient care [10,11].
Various cancers, antibiotic therapy, and myelosuppression by conditioning therapy have been associated with microbial dysbiosis [12]. Nevertheless, little research has explored the relationship between specific microbes and OM. In 2013, Ye et al. analyzed the microbial diversity and richness of pediatric patients with malignancies and found a non-statistically significant increase in the relative abundance of the genera previously associated with OM, such as Enterococcus, Escherichia, Porphyromonas and Pseudomonas [13]. In agreement with these results, a recent study in chemotherapy patients with solid tumors reported that mucositis severity was positively correlated with three salivary gram-negative bacilli but negatively correlated with 24 commensal species, including Streptococcus, Actinomyces, Gemella, Granulicatella, and Veillonella genera [14]. Additionally, in patients subjected to conditioning therapy, changes in microbiome diversity and similarity occur, which may increase susceptibility to a proinflammatory state [15,16]. We have also proposed a model integrating the role of the oral microbiome in cancer therapy-induced OM, within the five-step model established by Sonis [4,17].
Approaches aimed at microbial species-level characterization, host/microbe interaction, and functional genomics hold the promise of providing insight into the mechanisms of OM development and the discovery of possible biomarkers of risk assessment.
The purpose of this study was to conduct a species-level longitudinal analysis of the oral microbiome in hematological cancer patients who developed moderate OM compared to those who developed none or mild OM, as a result of conditioning therapy. Next-generation sequencing (NGS) of the V3-V4 hypervariable region of the 16S rRNA gene was used to identify over 700 oral bacterial species [18,19].

Patient recruitment
Patients diagnosed with hematological cancers scheduled for conditioning therapy were recruited at Carolinas Medical Center-Atrium Health, Charlotte, North Carolina and enrolled in the prospective cohort study: 'Multicenter Study on the Burden of Illness of Oral Side Effects from Conditioning Therapy Before Stem Cell Transplantation: Ora-stem Study [10]. The study was approved by the Institutional Review Board (IRB) and patients provided informed written consent. Patients with non-cancer-related hematological disorder (e.g., immunodeficiency) were excluded. Demographics of the hematological cancer patients (n = 22) used in this study are shown in Table 1, which excludes loss-to-follow up (e.g., patient deceased) or the absence of NGS processing due to low sample quality. The patient cohort was stratified for the time points "baseline" [T 0 ] prior to conditioning, time of 'oral mucositis occurrence' [T Muc ] , one-year post-transplant [T Year ]), the number of paired samples obtained from four different oral sites, and the patient's OM status: Muc2 group (moderate OM) and Muc0-1 group (no OM or mild OM) ( Table 1).
De-identified clinical characteristics of the patient cohort (n = 22) included: (i) worst WHO OM scores, (ii) hematological stem cell transplant type (autologous (4 Muc0-1 and 2 Muc2 patients)) or allogeneic (9 Muc0-1 and 7 Muc2 patients), (iii) hematological cancer diagnoses, paired oral samples (T 0 to T Muc or T 0 to T Year ), and paired for all three time points, (iv) whether or not patients had received antibiotic prophylaxis treatment (levofloxacin) within 2 weeks prior to oral sample collection at T 0 pre-conditioning (antibiotics treatment before T 0 beyond two-weeks possible, but no antibiotics treatment between T 0 and T Muc ), and (v) chemotherapy with or without total body irradiation (TBI) (Supplementary Table 1).

Sample collection and processing
Oral samples were collected at least one following dental manipulations (e.g., eating, oral hygiene). Stimulated saliva (S) samples, and swab samples from buccal mucosa (B), superficial supragingival plaque (P), and tongue (T) were collected. Samples were obtained at baseline (T 0 : one to 8 weeks pre-conditioning), days of likely occurrence of OM (T Muc : day of transplant (day 0) postconditioning, day 7 and day 14 post-transplant) and 12 months post-transplant ±30 days (T Year ). Saliva samples were collected at T 0 and T Year only.
Stimulated saliva was collected while chewing unflavored and unsweetened gum base (The Wrigley Company, Mars, Inc., Chicago, IL). The samples (1-2 mL) were centrifuged (2,600 x g; 4°C; 15 min) for pelleting. Other oral samples (B, P, T) were suspended in nuclease-free PBS solution containing 0.04%  sodium azide and rotated (2 hat room temperature) to release bacteria. The suspensions were centrifuged (16,000 x g) [20]. All pellets were stored at −80°C. Bacterial genomic DNA was extracted from oral samples using the modified QIAamp DNA Mini Kit procedure (QIAGEN, Valencia, CA) per manufacturer's instructions. Human Oral Microbe Identification using Next Generation Sequencing (HOMINGS) was used to identify bacteria at the species and genus levels, and relative abundances were determined as previously described [18,19,21]. Briefly, the amplified 16S rRNA gene (V3-V4 region) was sequenced using a modified MiSeq NGS method (Illumina, Inc., San Diego, CA) [21]. Oral taxa identification and abundance were determined using the ProbeSeq program, in which sequence reads were first matched against ProbeSeq species probes, in a BLAST-type electronic 'e-hybridization' [18][19][20]. The number of sequence-reads matched to one probe of 737 probes total (620 ProbeSeq species probes and 117 genus probes) was counted. The relative abundance of ProbeSeq matched genus probes and species probes was determined for each patient sample.

Bioinformatics analysis
Beta-diversity analysis The overall analytical strategy is presented in Figure 1. Groups of patients for time period analyses were Set-All, Subset-Common, Subset-noAntibiotics, and Subset-TBI (Supplementary Table 1). The time period analyses T 0 to T Muc and T 0 to T Year were performed for Muc0-1 and Muc2. For Muc0-1 and Muc2, there were, respectively, 43 and 36 total paired sample counts for T 0 to T Muc and T 0 to T Year combined, i.e., overall 114 individual oral samples forming pairs in effect were collected and sequenced. Cross-sectional analyses were performed to compare Muc0-1 and Muc2 at T 0 , T Muc and T Year . For both OM groups, if the same highest WHO score was present during more than one time point at T Muc , relative abundance data were averaged. Relative abundance data were squared-root transformed and converted to Bray-Curtis similarity matrices followed by PERMANOVA analyses using a mixed-model with unrestricted permutation of raw data, 9,999 permutations, and type III partial sum of squares, in the PRIMER v7 program (PRIMER-E Ltd., Ivybridge, UK), as previously implemented by our  1. Global oral microbiome beta-diversity changes were determined for the time periods T 0 to T MUC and T 0 to T Year , based on species and genera relative abundances (RA) in oral samples of hematological cancer patients (n = 22) who developed oral mucositis with World Health Organization OM scores 0, 1, or 2 post-conditioning therapy. 2. Post-conditioning oral microbiome relative abundance fold increase or decrease (RA-FcD = (Final RA minus Initial RA)/Initial RA, with addition of a pseudo-count to raw abundance data) within Muc0-1 group (no OM or OM score 1) and Muc2 group (OM score 2) was determined for the time periods T 0 to T Muc and T 0 to T Year . 3. Linear discriminant analysis effect size (LEfSe) analysis was performed for Muc0-1 and Muc2 groups using 'Group' for the LEfSe input option 'Class' and 'Patient' as the LEfSe input option 'Subject' in the LEfSe online program. 4. RA-FcD data were used for post-LEfSe sub-analyses for most discriminative taxa to compare changes in the microbial communities of patients who developed OM score 0-1 to those in patients who developed OM score 2, following conditioning for the time periods T 0 to T Muc and T 0 to T Year . Time points were: 'T 0ʹ (time at which cancer is present, pre-conditioning); group [20]. Fixed factors in the PERMANOVA longitudinal design were 'Time' (T 0 to T MUC and T 0 to T Year ) and 'Site' (up to four levels: S, B, P and T). In this design, the random factor 'Treatment' was coded as a single factor with four possible outcomes (i.e., four levels: Antibiotics [yes/no] and TBI [yes/no]) to control for degrees of freedom and was nested into 'Site'. The random factor 'Subject' was nested in 'Treatment' and 'Site' in the longitudinal analysis but was replaced by 'Group' (Muc0-1, Muc2) in the cross-sectional analysis. Monte-Carlo corrected p-values (α = 0.05) were determined and PCoA comparisons of relative abundance sample data of T 0 vs. T MUC and T 0 vs. T Year for Set-All groups Muc0-1 and Muc2 were carried out in PRIMER v7 to visualize significant comparisons (PERMANOVA) (Supplementary Figure 1a and 1b).
Linear discriminant analysis (LDA) effect size (LEfSe) Taxonomy levels were added manually to ProbeSeq derived datasets for Set-All groups Muc0-1 and Muc2. The tabular text file was formatted to perform LEfSe using online tool Galaxy v1.0 [22]. LEfSe data input consisted of 'Group' (Muc0-1and Muc2) for the LEfSe input option 'Class' and 'Patient' for the LEfSe input option 'Subject' in the LEfSe online program Galaxy v1.0 [22]. Using the 'one-against-all' strategy for multi-class analysis [23], the factorial Kruskal-Wallis test and pairwise Wilcoxon signedrank test were set at a Monte-Carlo significance level α = 0.05 to calculate LDA scores. The log LDA score was set at a threshold >0 and used to generate a cladogram representing the hierarchy of all significant biomarkers and a histogram of the top biomarkers, plotted at the genus and species levels [22].
Post-LEfSe proteobacteria sub-analysis Genus and species probes belonging to the Proteobacteria phylum (n = 91 probes) were divided into Gammaproteobacteria (n = 30) and remaining (n = 61) probes. Longitudinal PERMANOVA analyses T 0 to T Muc and T 0 to T Year were performed for both Muc0-1and Muc2 groups using the species and genera probes for 'Gammaproteobacteria' and 'all Proteobacteria'. PCoA were generated in PRIMER v7 .
Further, the relative abundance fold change difference (RA-FcD), based on the formula [('final RA' minus 'initial RA') divided by 'initial RA'], for each oral site, was determined from the relative abundance data obtained after adding a pseudo-count of +1 to each raw data count, as previously described [20,24]. The average RA-FcD sum per each sample site for all the patients in the Muc0-1 or Muc2 group, based on the total abundance of the 30 Gammaproteobacteria probes (not total of 737 probes) per sample, was calculated for the T 0 to T Muc time period. The sample site representation for Muc0-1 and Muc2 groups was similar, i.e., 3/2 B, 6/4 P, and 12/8 T samples, respectively. A Fisher's exact test (α = 0.05) was then performed for each time point in R v3.4.3 [R 25] to determine if the average RA-FcD sum per oral site across time points were significantly different between Muc0-1 and Muc2. Specificity and sensitivity and Receiver Operating Characteristic (ROC) curves were determined [26].

Results
Abundance data, species detection, and alpha-diversity Sequencing reads matched to 737 total probes (620 species and 117 genus probes) for all paired samples from our patient cohort at time points T 0 , T Muc and T Year (overall total of 114 samples forming pairs collected and sequenced) are summarized in Supplementary Table 2. Unmatched reads were excluded from relative abundance determinations. For all samples with sequencing data available, 397 of 620 species probes and 85 of 117 genus probes had at least one read matched to a genus or species probe. Derived relative abundances were used for analysis of patient subsets Set-All, Subset-Common, Subset-noAntibiotics and Subset-TBI (Figure 1). In the T 0 vs.T Muc comparison for Muc0-1, there was no alphadiversity difference for any subsets (Simpson and Shannon) (Set-All result is shown in Supplementary  Table 3). However, in the T 0 vs. T Year comparison, there was a significant alpha-diversity difference for Set-All (Shannon). In the T 0 vs. T Muc comparison for Muc2, there was a significant difference in alphadiversity for Set-All (Simpson and Shannon), while there was none T 0 vs. T Year .
We also found that the ranges of taxa detected per patient for Muc0-1 and Muc2 groups combined across T 0 to T Muc and T 0 -T Year , were 43.42 to 87.86 for genera and 129.5 to 281.14 (Supplementary Table 3).

Beta-diversity analyses
The baseline T 0 cross-sectional PERMANOVA analyses of Muc0-1 vs. Muc2 resulted in no significant betadiversity difference for any patient subsets, suggesting a similar oral microbiome 'starting point' (data not shown). There were beta-diversity differences for both time periods T 0 to T MUC and T 0 to T Year for Set-All and Subset-noAntibiotics, and Subset-TBI for T 0 -T MUC only, when combining data from Muc0-1 and Muc2 (p < 0.05) (Supplementary Table 4).
In the Set-All T 0 toT Muc comparison, differences in beta-diversity changes were noted when Muc0-1 and Muc2 groups were analyzed individually. There were no significant beta-diversity changes for Muc0-1 subsets. Although stratification reduces sample size, the Muc2 T 0 to T Muc results suggest that betadiversity changes in relation to moderate OM (i.e., p < 0.05 for Set-All, Subset-noAntibiotics, Subset-TBI, marginal p-value for Subset-Common) were possibly related to the absence of antibiotics treatment. Indeed, 7 of 13 Muc0-1 patients received TBI and were treated with antibiotics within two weeks prior to sampling at T 0 pre-transplant, compared to a single patient treated with antibiotics among 9 Muc2 patients including 6 treated with TBI for conditioning but not antibiotics within two weeks prior to T 0 (Supplementary Table 1). In the T 0 -T Year comparison, for similar numbers of patients and samples, there was no beta-diversity change for Set-All in Muc0-1 group, as opposed to Set-All in Muc2 (Table 2). Overall, although Muc0-1 and Muc2 groups had similar patient and sample counts representation, most of the comparisons for the Muc2 group (Set-All and Subsets) were significant or marginally significant, whereas the Muc0-1 comparisons generated only one significant p-value (Set-all T 0 vs.T Muc ) ( Table 2).

LEfSe analysis of Muc2 and Muc0-1 groups
LEfSe [22] identified 54 differential features for Muc0-1 and Muc2 in Set-All consisting of 737 total probes (620 species and 117 genus probes). Proteobacteria (log LDA ≈ −0.01) was identified as the leading differential feature between Muc0-1 and Muc2 (Figure 2(a, b)). Veillonella made up the largest differential feature for Muc0-1 (log LDA ≈ 0.003). When comparing the list of genera determined by LEfSe to the list representing 10 species/genera with the largest significant RA-FcD changes for T 0 to T Muc , Abiotrophia, Capnocytophaga, Gemella, Haemophilus, Lactobacillus, Prevotella and Veillonella genera overlap (Figures 2(b) and 3).   Figure 2 (e, f, g, h). In addition, Fishers' exact test was performed using the average of the sum of RA-FcD per oral sample site across Muc0-1 and Muc2 patients for each of the 30 Gammaproteobacteria probes. There was an overall average RA-FcD difference between the Muc0-1 and Muc2 group in the T 0 vs.T MUC (p = 4.03e-3) and T 0 vs. T Year (p = 1.17e-10) comparisons. Notably, from T 0 to T MUC and T 0 to T Year time periods the relative abundance of Haemophilus parainfluenza decreased, on average per sample site, while it increased for Escherichia and Shigella genera in Muc2 compared to Muc0-1. ROC curves showed that the specificity for the 30 Gammaproteobacteria probes was 0.75 regarding T 0 -T Muc RA-FcD changes and 0.724 for the T 0 to T Year comparison ( Supplementary Figures 3a and 3b). Linear discriminant analysis Effect Size (LEfSe) was performed to determine distinct microbiome features in oral samples (stimulated saliva and swabs of buccal mucosa, superficial supragingival plaque, and tongue) of the Set-All patient cohort (n = 22) for the T 0 to T Muc time period. T 0 corresponds to the time at which cancer is present, but prior to conditioning therapy and T Muc corresponds to the time at which patients developed oral mucositis (OM) following conditioning therapy; or the time at which patients did not develop OM the day of transplant, and day 7 and day 14 post-conditioning. Muc0-1 group represents patients with no OM or OM score 1 and Muc2 group represents patients with OM score 2, post-conditioning. Graphical representations of LEfSe results are shown: (a). Cladogram representing discriminant features (n = 54) that are potential biomarkers for Muc0-1 (shown in green) and Muc2 (shown in red) with a hierarchy based on taxonomy with 'Group' as the LEfSe input variable 'Class' and 'Patient' as the LEfse input variable 'Subject'. (b). Horizontal histogram of discriminant features (n = 54) for Set-All potential biomarkers for Muc0-1 (shown in green) and Muc2 (shown in red). All potential biomarkers are based on 737 probes total (620 ProbeSeq species probes and 117 genus probes).

Discussion
This is the first one-year follow-up study investigating oral mucositis-associated microbiome profiles in multiple oral sites of hematological cancer patients undergoing conditioning therapy prior to hematopoietic stem cell transplant. The difference between patients who develop OM with a WHO score of 2 and those who do not (scores 0-1) is the presence of ulcer(s) in the oral mucosa. Therefore, the comparison Muc0-1 vs. Muc2 provided the opportunity for better understanding oral microbiome shifts associated with the 'threshold' of presence of ulcers in oral cavity. More severe OM (scores of 3 and 4) might involve additional sequential, possibly confounding host responses affecting the oral microbiome. Increased risk for bacteremia or infection may then be associated with the development of pseudomembranes on mucosal ulcerations with additional colonization of bacterial species capable of promoting further inflammation [27]. However, previous studies have positioned the oral microbiome as an OM exacerbator rather than an initiator [4,15]. Here we sought to determine which microbial sub-community may be disrupted, thereby leading to OM exacerbation.
Using V3-V4 16S rRNA gene next-generation sequencing along with the ProbeSeq species and genera identification program, we were able to detect at least 482 bacterial species from the 737 species and genus probes (65.4%) for our patient cohort, considering that approximately 700 predominant taxa are represented in oral cavity [28]. Changes in alphadiversity consistently occurred during T 0 -T MUC time period for the Muc2 group (Simpson and Shannon). Conversely, there were no changes in alpha-diversity for Muc0-1 for that time period.
The PERMANOVA results (Supplementary Table 4) showed that by analyzing Muc0-1 and Muc2 patients combined, the beta-diversity significantly changed for the T 0 -T MUC time period with Set-All (including the AB 'treatment' factor as defined earlier) and Subset-noAntibiotics yielding identical p-values. When analyzing Muc0-1 and Muc2 separately (Table 2), there was a significant beta-diversity change for both groups for the Set-All subset at T Muc . Significant beta-diversity changes for the T 0 -T MUC time period were more frequent among Muc2 subsets. For this time period, none of the Muc2 patients were exposed to antibiotics treatment close to T 0 , while 58% Muc0-1 patients were (Supplementary Table 1, Table 2). However, large wellcontrolled studies did not demonstrate that antibiotics prophylaxis reduces incidence or severity of OM [11]. Additionally, there was a beta-diversity change at T Year for Muc2, but not for Muc0-1, in Set-All. These results taken together suggest that the oral microbiome community in our patients' cohort who developed moderate OM did not return to a normal state, in contrast to Muc0-1 patients. Long-term microbial shifts have been associated with an OM history in a patient population with Fanconi Anemia who underwent conditioning therapy [29].
Proteobacteria showed the largest linear discriminant analysis score (log LDA = 0.011), potentially explaining why microbial communities did not recover after one-year. LEfSe also showed that many more features distinguish Muc2 compared to Muc0-1. Furthermore, Gammaproteobacteria were major contributors to the differences observed between Muc0-1 and Muc2, as shown by PERMANOVA and ROC analyses.
Following the post-LEfSe sub-analysis on the 30 Gammaproteobacteria probes, we observed the largest average of RA = FcD sum per oral site ratios for the T 0 vs.T MUC time period, between Muc0-1 and Muc2, corresponding to a decrease in Muc2, for H. parainfluenzae and the Acinetobacter genus, whereas an increase was observed for the Pseudomonas genus, Pseudomonas fluorescens and Escherichia and Shigella genus.    [4.8]. Although reduced, these differences subsisted relatively similar, T 0 -T Year . Previous studies in patients with sarcoidosis and oral lichen planus have provided evidence that H. parainfluenza is protective of mucosal surfaces, especially strains that are strong biofilm producers [30,31]. Although our data indicate a clear trend in this respect, abundance data alone cannot predict risk for OM, since strains of H. parainfluenza would need to be identified and functional assays to assess biofilm production capacity be performed.
Considering the low LDA effect sizes in our study, more research with a larger sample size, bacterial strain level identification, and functional analyses related to biofilm production may be required to determine composite microbiome signatures providing adequate diagnostic value to predict OM severity. Indeed, we were unable to differentiate Escherichia and Shigella species identified by a genus probe and therefore we have no precise knowledge about the potential for invasiveness affecting the oral mucosa.
In addition, while this study focused on the characterization of the oral microbiome profiles in hematopoietic stem cell transplant (HSCT) patients who develop moderate OM compared to patients with OM scores of 0-1, understanding sequential changes occurring in patients who develop severe OM (scores 3-4) would provide additional clinical significance to our study. In particular, one might show that certain microbiome changes for patients who transition from moderate to severe OM would be characteristic of patients who experience the most pain. Furthermore, the absence of salivary flow data at T Muc , precluded an analysis including such variable, which along with other factors (e.g., conditioning regimen; engraftment type) can influence the microbiome composition and activity. Another limitation of the study was the loss of samples' data, due to patient death or absence of multiple visits, therefore, a larger patient cohort and further collaborative efforts would likely provide more comprehensive insights into OM development in HSCT patients.
In conclusion, based on beta-diversity results, hematological cancer patients who develop moderate OM following conditioning therapy experience a lasting change in Proteobacteria subcommunities up to one-year post-conditioning. More research is needed to define susceptibility to OM, including investigation of host responses. Oral microbiome profiling and functional analysis could provide new means to prevent or mitigate OM in this compromised patient population.