Cutaneous dysbiosis characterizes the post–allogeneic hematopoietic stem cell transplantation period

Skin graft-versus-host disease (GVHD) affects up to ∼70% to 80% of allogeneic hematopoietic stem cell transplantation (HSCT) patients,¹ and yet, the underlying pathogenesis is poorly understood. Abnormalities in the gut microbiome have been associated with the development of GVHD² and mortality,³ but the impact of the skin microbiome on HSCT outcomes is unknown. Chronically immunosuppressed patients have a disrupted skin microbiome,⁴ and skin dysbiosis has been linked to inflammatory cutaneous disorders,⁵ facilitating Th17 polarization,⁶ and skin disease progression.⁵ Skin commensals are critical for maintaining local immune homeostasis, keeping pathogenic organisms in check,⁷ and are necessary for the development of the skin’s regulatory T-cell compartment.⁸ Based on this, we hypothesized that the skin microbiome is altered among allogeneic HSCT patients and might be an important biomarker for outcomes. We also hypothesized that skin dysbiosis drives the innate, inflammatory response characteristic of acute GVHD (aGVHD) of the skin, and keratinocytes are active mediators of this inflammation. Therefore, we aimed to describe the diversity and composition of the cutaneous microbiome among allogeneic HSCT patients across the 90-day posttransplant recovery period. We obtained data on the development of aGVHD as well as mortality.

After obtaining approval from the Duke Institutional Review Board and written informed consent, we obtained skin swabs from the forearm (arm, n = 61) and retroauricular areas (ear, n = 63) of 12 adult patients (supplemental Table 1) and from 4 primary caregivers who consented to participate in the study to identify any potential microbial transmission patterns or risk reservoirs between caregivers and patients. Skin swabs from HSCT patients were obtained before conditioning chemotherapy/radiation and at 1 week, 1 month, and at 3 months after allogeneic HSCT, corresponding to landmarks in the immune recovery process after HSCT. The sites were selected to represent oily and dry skin and body regions commonly affected by aGVHD. We had a maximum of 1 sample from each body site and time point, with some missing samples. Shotgun metagenomic sequencing of samples was performed as previously described.⁹ Patient characteristics and outcomes, development of acute skin GVHD, and 1-year mortality were collected (supplemental Table 1). Other pertinent clinical characteristics such as the development of systemic infections and antibiotic exposure were also collected (supplemental Table 2). We had previously used archival, formalin-fixed aGVHD skin biopsy tissue from an independent set of patients to confirm the expression of key antimicrobial cytokines S100A8 and IL36 by immunohistochemistry and immunofluorescence in keratinocytes of patients with aGVHD (data not shown). Because keratinocytes are more likely to be in contact with the cutaneous microbiome and are capable of mounting an inflammatory response, we performed single-cell RNA sequencing (scRNA-seq) on the epidermal compartment from a patient who developed skin aGVHD and compared it with that of a healthy volunteer. Anonymized truncal skin from a healthy adult was obtained through an institutional review board–approved tissue bank protocol and was used as healthy control for scRNA-seq and immunohistochemistry.

We removed adapter sequences from raw metagenomics reads using Trim Galore¹⁰ and discarded human reads by mapping to hg19 genome using Bowtie2.¹¹ Low-quality reads were discarded using Sickle.¹² For taxonomic annotation, we mapped cleaned metagenomic reads to PanDB¹³ and reassigned ambiguous reads using Pathoscope 2.¹⁴ We compared the microbial composition of the skin microbiome in our cohort with a previously published data set from healthy volunteers.⁹ We compared the microbiome characteristics between patients in our cohort who developed acute cutaneous GVHD and those who did not, as well as between patients who survived 1 year after HSCT and those who did not.

To identify plasmid-borne antimicrobial resistant genes (ARGs), we began by concatenating all samples and assembled sequence reads into contigs using MEGAHIT 1.0.5.¹⁵ We discarded contigs <1 kb and predicted plasmids’ contigs using PlasFlow.¹⁶ Next, we used Prodigal¹⁷ to predict genes in plasmids and subsequently identified ARGs using DeepARG.¹⁸ We mapped each sample to the pool of ARGs using Bowtie2 and reassigned ambiguous reads using Pathoscope2.¹⁴ 

We estimated subspecies/strain-level growth rate using strain-level metagenomic estimation of growth rate.¹⁹ We retrieved 187 strains of C acnes from NCBI (National Center for Biotechnology Information) GenBank and constructed the species database using a single-nucleotide polymorphism (SNP) assignment threshold of 0.6, without iterative clustering, resulting in 8 clusters (supplemental Table 2). Because strains of C acnes have low within-species genetic diversity, we detected clusters at a 0.5× coverage cutoff and reported growth rates for clusters with coverage ≥1× and nonzero SNP count ≥100 in a given sample, as described previously.¹⁹ We used a cluster detection threshold of 0.2 to identify subspecies clusters in samples.

We downloaded 760 453 viruses or viral fragments and their annotations from integrated microbial genomes/virus²⁰ to create our viral database. We mapped each sample to the database using Bowtie2,¹¹ reassigned ambiguous reads using Pathoscope2,¹⁴ and excluded viruses present in <3 samples. This resulted in hits to 17 541 viruses/viral fragments annotated as bacterial or eukaryotic viruses. We performed additional analysis for polyomavirus, because we had a previous analysis (not shown) that suggested enrichment for polyomavirus but not for others, such as human papilloma viruses, in the skin of HSCT patients. Polyomaviruses are relatively limited in types compared with other cutaneous viruses and yet have relatively high prevalence and can persist in the skin of healthy individuals.²¹ In addition, cutaneous polyomaviruses can cause significant cutaneous disease in immunosuppressed individuals. For strain-level polyomavirus analyses, we downloaded all sequences annotated as “human polyomavirus” (n = 20) from NCBI GenBank and generated a phylogenetic tree using FastTree,²² based on whole-genome alignment generated using mafft.²³ 

We conducted all statistical analyses using R software. To calculate P values between groups for repeated measures, we used unrestricted permutation implemented using a custom R script. We first computed a test statistic from an observation and generated 500 permutations by shuffling the dependent variable. The F statistic was used as the test statistic in the permutation for analysis of variance. We expressed the P value as the proportion of permutations resulting in a larger test statistic than the observed statistics. P values were adjusted under the Benjamini-Hochberg procedure.

After 2 brief washes in cold phosphate-buffered saline, the freshly obtained skin biopsies were incubated in 4-mL dispase (1 unit/mL; Thermo-Fisher Scientific, Waltham, MA) at 4°C overnight. Next day, the epidermis and dermis were separated from each other and chopped into small pieces with a pair of scissors. The epidermal pieces were digested with 0.25% trypsin for 15 minutes, and the dermal pieces were digested with collagenase B (1 mg/mL; Sigma, St. Louis, MO) for 1 hour, followed by the addition of 0.25% trypsin and incubation for another 15 to 30 minutes. During digestion, the tubes are placed at 37°C water bath with agitation at 240 revolutions per minute. The digestion mixtures were neutralized in 10% fetal bovine serum in Dulbecco’s modified Eagle medium and filtered through 70-μm and then 40-μm strainers. The cell suspension was centrifuged at 300g for 5 minutes to obtain cell pellets, washed 3 times, and then resuspended in 0.04% bovine serum albumin at 1000 cells per μL. After quality check, cells were immediately used for complementary DNA amplification and chromium library construction using scRNA-seq library kit v2, according to the manufacturers’ recommendations (10X Genomics) by the Molecular Genomics Core facility at the Duke Molecular Physiology Institute. Complementary DNA libraries were sequenced on the Illumina HiSeq 4000, with a read length of 150 bp by the Duke Center for Genomics and Computational Biology core facility.

Analysis of the scRNA-seq data followed the manufacturers’ protocol (10X Genomics). Briefly, demultiplexed raw base call files were generated by Illumina sequencers into FASTQ files, and alignment was performed to the h19 human transcriptome. Filtering, barcode, and unique molecular identifier counting were performed using Cell Ranger software version 3.1 and default parameters (10X Genomics). Secondary statistical analysis was performed using the R package Seurat, following the standardized pipeline.²⁵ Data were first normalized and scaled after basic filtering for minimum gene and cell observance frequency cutoffs. We then examined the data and performed further filtering (nFeature 7500; nCount 100 000; and percent_mt 0.1). The removal of further technical artifacts was performed using regression methods to reduce noise.

After the quality control procedures were complete, we further performed linear dimensional reduction, calculating principal components using the most variably expressed genes in our data set. Significant principal components for downstream analyses were performed similar to the methods previously reported by Macosko et al.²⁶ Cells were groups into an optimal number of clusters for de novo cell type discovery using Seurat’s FindNeighbors and FindClusters functions. Graph-based clustering approaches with visualization of cells were achieved through the use of manifold learning technique Uniform Manifold Approximation and Projection, which reduces the information captured into 2 dimensions. In total, 8861 and 3393 cells from normal skin and aGVHD skin, respectively, were used for downstream analyses.

Our analysis revealed a total of 24 clusters in the epidermis, which includes 11 keratinocyte clusters visualized using Uniform Manifold Approximation and Projection. Keratinocyte clusters (0, 1, 2, 3, 8, 9, 10, 11, 12, 15, and 18) were further divided into distinct subtypes, including basal, suprabasal, and terminal keratinocytes. Terminal keratinocytes were characterized by the expression of involucrin (IVL), loricrin (LOR), and filaggrin (FLG). Differentially expressed genes for a particular cell type were defined by a P value threshold <.05 and log fold change >0.25. Data were subjected to Ingenuity Pathway Analysis (Qiagen, Redwood City, CA) to identify the top differentially expressed genes.

The skin microbiota among our HSCT patients was enriched for predicted pathobionts, gram-negative bacteria, fungi, and viruses, from before transplant and throughout the 90-day post-HSCT period compared with the microbiome expected in healthy individuals⁹ (Figure 1A). The skin microbiome composition among our HSCT cohort was different from a previously published data set from healthy volunteers, making direct comparisons difficult (supplemental Figure 1).⁹ Bray-Curtis dissimilarity index (BCI) demonstrates separate clustering of microbial communities from HSCT patients and control healthy caregivers (supplemental Figure 2A). There was no separation of clusters when HSCT patients were stratified based on diagnosis or age (supplemental Figure 2B-C). Taxa of potential pathobionts such as beta papillomaviruses, actinomyces, and gram-negative bacteria are enriched in the skin of our HSCT patients before HSCT (supplemental Figure 2D). There was no association between skin microbiome alpha diversity (within-sample diversity; Shannon diversity index) and 1-year mortality or the development of aGVHD at any time point (Figure 1B). We calculated the BCI (0 = identical species composition and abundance; 1 = completely dissimilar) between the pretransplant microbiome and other time points to describe the stability of each individual HSCT patients’ microbiome across the posttransplant recovery period. The average BCI comparing the skin microbiota of patients using the forearm specimens across the 90-day post-HSCT recovery period with their pretransplant microbiota was significantly lower among patients who survived the first year after HSCT than those who died (P = .046; Figure 1C). This suggests that patients who survived 1 year after HSCT had a forearm skin microbiome composition and diversity that was more similar to their pretransplant skin microbiome than those who died, but it is not clear what the clinical significance of this is. There was no significant difference in the BCI between those who developed aGVHD and those who did not develop aGVHD.

Comparing between individuals, the BCI of the arm microbiome across the 90-day posttransplant period was lower among patients who survived than those who did not (P < 1.0 × 10^–20; Figure 1D). The pre-HSCT arm microbiome BCI among patients who survived was half the BCI among patients who died (0.47 vs 0.93; P = 1.1 × 10^–5; Figure 1D) and was also significantly lower than healthy controls and patients who survived (0.95; P = 4.0 × 10^–7), healthy controls alone (0.92; P = 5.0 × 10^–4), and healthy controls and patients who died (0.98; P = 9.0 × 10^–6; Figure 1E). We observed differential expressions of bacteria and viruses but not fungi (supplemental Figure 3), based on our outcome.

Analysis of the arm microbiome across our entire observation period demonstrated that the delta/epsilon and gamma classes of Proteobacteria, a major phylum of gram-negative bacteria, which includes Escherichia coli and Pseudomonas, were enriched in patients who survived compared with those who died (P = .022; Figure 2A). No differential expression was observed among those who developed aGVHD of the skin. We investigated a possible role of ARGs as an explanation for microbiome stability among 1-year survivors. We identified and classified patients according to the presence of plasmid-borne ARGs present at the pre-HSCT time point (Figure 2B). Patients clustered based on survival, primarily driven by fosmidomycin resistance genes and unclassified classes of ARGs.

Because commensals are important in the immunologic homeostasis of the skin, we estimated the growth rate of individual C acnes strains using strain-level metagenomic estimation of growth rate.¹⁹ Starting at 1 month after HSCT, the growth rate among patients who developed aGVHD was significantly lower than those who did not (P = .016; Figure 2C). By 90 days after HSCT, patients in the aGHVD group had lost C acnes strains, specifically strains belonging to cluster 1 (Figure 2D; supplemental Figure 4; supplemental Table 3).

We investigated the relative abundance of viruses in the skin microbiome of our HSCT patients. Patients who did not survive the 1-year post-HSCT period had an enrichment of bacteriophages on the arm, specifically Proteobacteria phages but not eukaryotic viruses (P = .012; Figure 3A). We did not find a significant relationship between the presence of phages or eukaryotic viruses and aGVHD (Figure 3B). There was a positive correlation between phages and delta/epsilon Proteobacteria at week 1 but a negative correlation at 3 months in the arm skin microbiome among patients who died before 1 year after HSCT (Figure 3C).

Because polyomaviruses have been linked to cutaneous tumors and eruptions among immunosuppressed patients,²⁷ we investigated the relative abundance of polyomaviruses in the skin microbiome among our cohort. The relative abundance of polyomaviruses in both the arm and retroauricular skin microbiome was significantly higher among patients who died than those who survived the 1-year post-HSCT period (arm, P = .006; ear, P = 2 × 10^–16; Figure 4A). Several of polyomavirus strains were identified among patients who died before 1 year after HSCT (Figure 4B). As proof of concept to identify whether polyomavirus strains were shared between patient and caregiver, we used single-nucleotide polymorphism analysis to compare polyomavirus genomes. We identified shared polyomavirus strains between 1 caregiver and patient pair (Figure 4C; supplemental Figure 5). Overall, polyomavirus strain composition tended to be individual specific and temporally stable over different time points (Figure 4D).

Using scRNA-seq, we identified 24 clusters of cells from the epidermis obtained from lesional skin with new onset aGVHD and a healthy volunteer (Figure 5A). We distinguished terminally differentiated keratinocytes among total keratinocytes based on their expression of involucrin (IVL), loricrin (LOR), and filaggrin (FLG). The top 3 differentially expressed genes in total keratinocytes (Figure 5B) and terminal keratinocytes (Figure 5C) in aGVHD skin compared with cells obtained from the skin of a healthy volunteer were the antimicrobial peptides S100A8, S100A9, and S100A7 (supplemental Table 5). Elafin (peptidase inhibitor 3/PI3), a serine protease inhibitor with antimicrobial properties, which has been proposed as both a biomarker for cutaneous aGVHD²⁸ and a predictor indicator of poorer outcomes,²⁹ was the fourth most highly differentially expressed genes in both keratinocytes and terminal keratinocytes. interleukin-36γ, whose expression is strongly induced by microbial triggers and whose overexpression is linked to the immunopathology of psoriasis,³⁰ was strongly differentially expressed, as were keratins 6a and 16, which are induced in response to epidermal barrier stress³¹ (Figure 5B-C).

The skin microbiome has a profound effect on skin immunity, and many protective aspects of this immune modulation are compartmentalized from systemic immune effects.^32,33 Disruption of skin microbiota is intimately associated with inflammatory skin diseases such as eczema and psoriasis, and this process is exacerbated by epidermal barrier failure.³⁴ Several factors could contribute to skin dysbiosis among allogeneic HSCT patients. Disruption in the local skin environment, expected after radiation and chemotherapy, such as changes in sebocyte activity and moisture of the skin, as well as systemic immune deficiency, can alter the landscape of the cutaneous microbiome composition.^35,36 The ubiquitous presence of pathogenic microbes, which are known triggers of the innate immune response, coupled with an increased keratinocyte antimicrobial inflammatory response in patients with aGVHD, suggest a role for the host–microbiome interaction in either the initiation or propagation of cutaneous aGVHD. Skin commensals are important for the development of the local, regulatory immune compartment of the skin.⁸ On the contrary, some authors have observed that commensals can increase skin inflammation and impair wound healing.³⁷ C acnes is part of the normal flora of the skin. Because of the limitations of this study, we cannot draw any conclusions regarding either the cause or impact of the loss of C acnes in the post-HSCT period. Larger, prospective studies that can correlate host changes such as the disruption skin barrier integrity and environment, changes in local and systemic immune compartments, and antibiotic use are needed. It remains to be seen what the overall impact of the cutaneous microbiome is on the local immune homeostasis and reconstitution after allogeneic HSCT.

We observed that patients who survived were more likely to maintain their pre-HSCT microbiome composition and diversity in the 3-month posttransplant recovery period, although there may be some fluctuation. In our cohort, the skin microbiome composition among patients who survived were more similar to each other before transplant than those who did not survive or healthy/caregiver controls. In addition, the microbial community composition of those who survived was dissimilar from those of healthy/caregiver controls. The exact composition of this “survival signature” of the skin microbiome, which seems to be evident even in the pretransplant period, and the immunologic function or status that it reflects need to be explored further. Because systemic immunodeficiency can alter skin microbiome composition, the changes we observed may be a reflection of the patient’s overall immune status. Infection is a leading cause of non-relapse–related mortality among allogeneic HSCT patients, and the skin may also serve as a reservoir of pathogens.

We noted an inverse relationship between the enrichment of delta/epsilon Proteobacteria in survivors on one hand and the increased presence of bacteriophages in the skin microbiome of those who did not survive on the other. We inferred that the bacteriophages may be in the early, nonvirulent lysogenic phase at week 1 and in the virulent lytic phase at 3 months, and the enrichment of delta/epsilon Proteobacteria in survivors may be attributed, partly, to the absence of lytic phages.³⁸ 

Our observation that skin microbiome instability and increased viral load may be associated with 1-year mortality after HSCT implies a relationship between systemic immune status and the skin microbiome homeostasis. In healthy individuals, the skin microbiome at most skin sites is relatively stable over time.⁹ This stability may be an inherent characteristic of a healthier community. This is despite the fact that the configuration of the pre-HSCT skin microbiome is not normally observed among healthy individuals but is common among the pre-HSCT patients. All this suggests that stability, rather than a recapitulation of a “healthy” microbiome, may be an important biomarker for survival. A shifting skin microbiome may be a sign of worsening skin barrier or deteriorating immune surveillance, but this remains to be proven. The impact of skin barrier integrity in the development of aGVHD deserves additional studies. The impact of antibiotic and antiseptic use and the post-HSCT recovery environment on dysbiosis also need further exploration.

Our study is limited by the small number of patients and controls and the availability of sampling from other anatomic sites. In addition, it is difficult to standardize antibiotic exposure and aGVHD treatments among HSCT patients. Nevertheless, larger prospective studies are warranted. The relationship between the skin microbiome and gut microbiome may also be worth examining. Non-antibiotic–based interventions, such as changing the patient’s post-HSCT care to a home rather than a hospital-based environment, could potentially affect the skin microbiome and patient outcomes.

The authors acknowledge the assistance of the Duke Molecular Physiology Institute Molecular Genomics Core for the generation of single-cell RNA sequencing data for the manuscript.

This study was supported by the Duke Medicine Scholars award (A.R.C.), National Cancer Institute (NCI) grant R21 CA223419 (A.R.C. and J.Z.), NCI grant R01CA203950 (A.D.S., L.M.B., and N.J.C.), National Institute on Aging (NIA) grant R21AG066388 (A.D.S., L.M.B., and N.J.C.), the American Society of Hematology Scholar award (A.D.S.), and Duke Pepper Center/NIA grant P30-AG028716-13 Mini#6 (A.D.S.). J.O. is supported by the National Institutes of Health (NIH) grants DP2 GM126893-01, 1U54NS105539, 1 U19 AI142733, and 1 R21 AR075174, National Science Foundation grant 1853071, the American Cancer Society, and the Leo Foundation. A.E. is supported by NIH grant K99GM135540, the Jackson Laboratory Scholars Program, and the Maximilian & Marion Hoffman Foundation Postdoctoral fellow.

The funders had no role in the conceptualization, design, data collection, analysis, decision to publish, or preparation of the manuscript.

Contribution: A.R.C. and J.O. conceptualized, designed, and supervised the article; N.J.C., R.P.H., and A.D.S. designed and supervised critical aspects of this study; A.R.C., A.D.S., N.J.C., and J.O. obtained funding; A.R.C., C.P., K.R.N., A.D.S., E.F., J.Z., Y.J.J., M.V.L., L.M.B., and A.T.B. collected the data; A.R.C., A.E., R.P.H., E.F., J.O., A.D.S., A.J.P., V.J., S.G.G., and J.Z. analyzed and interpreted the data; A.R.C., A.E., and A.J.P. wrote the original draft of the manuscript; and A.R.C., A.E., R.P.H., A.D.S., J.Z., N.J.C., K.M.S., S.G.G., and J.O. reviewed the manuscript and provided substantial contribution to the analysis and interpretation of the data and editing of the manuscript.

Conflict-of-interest disclosure: The authors declare no competing financial interests.

Correspondence: Adela R. Cardones, Division of Dermatology, Department of Internal Medicine, University of Kansas Medical Center, 3901 Rainbow Blvd, Delp 4094, Mailstop 2025, Kansas City, KS 66160; email: acardones@kumc.edu.