Disturbances in microbial skin recolonization and cutaneous immune response following allogeneic stem cell transfer

Sampling

Adults (>18years; n=50) previously diagnosed with a hematologic malignancy and suitable allogeneic HSCT were recruited before conditioning at the Medical University of Vienna after ethical approval (ECS 1087/2016) and participated voluntarily and with fully informed written consent. Skin scrapes from different regions (palm, back and face) and stool were obtained before conditioning (n=180), the day of transplantation (Tx, n=171), at day 14 (n=174), day 100 (n=130) and day 365 (n=109). Additional samples were collected from healthy, age-matched volunteers (n=15) and empty tubes were used as negative controls (n=71). 6 mm skin biopsies from a subset of patients (n=9) were obtained at the same 5 time points, while additional biopsies were collected after aGVHD onset (n=12). aGVHD lesions and unaffected skin were also swabbed for isolation of Staphylococcus spp. (n=5). All samples from GVHD patients were collected before start of any treatment (Suppl. Table 13). GVHD was diagnosed according to the Glucksberg criteria for acute (13) and the NIH consensus criteria for chronic GVHD (14). All patients were administered prophylactic antibiotics (Suppl. Fig 3G) and other antimicrobial agents upon hospitalization (Acyclovir), during their neutropenic period until cell engraftment (Paromomycin, Metronidazol and Nystatin) as well as upon development of neutropenic fever (Piperacillin/Tazobactam).

16S rRNA gene amplicon sequencing

DNA was extracted using the QIAamp DNA Microbiome Kit with additional lysis steps using Mutanolysin (0.01mg/ml), Lysozyme (1mg/ml) and Lysostaphin (0.1mg/ml). Primer pairs 341F/785R (15) were used in a two-step PCR workflow (16). Amplicon sequence variants (ASVs) were inferred using the DADA2 R package (17)^,(18). After trimming, sequences were classified with SINA version 1.6.1 (19) and SILVA database SSU Ref NR 99 release 138 (20). Eukaryota, mitochondria, chloroplast, and known PCR contaminations (ASV_qdi_vob) were removed. Alpha diversity was calculated using the R ampvis2 v.2.6.5 (21) and Bray Curtis using the R package 2.5–6 (22). For differential abundance analysis only, skin/stool samples containing between 900/ 2000 and 27000/60000 reads were used from patients with a minimum of three skin and two stool samples. Additionally, samples with a GVHD after d100 were excluded. For pairwise differential abundance testing DESeq2 was run with default parameters except fitType=‘local’.

Immunofluorescence

Acetone-fixed cryosections from skin were hybridized with a 16S rRNA-probe (EUB338: 5′-GCTGCCTCCCGTAGGAGT-3′ (23)) as well as Staphylococcus-specific probe (STAPHY2: (5’-TCCTCCATATCTCTGCGC-3’ (24)). Reverse complements were used as negative controls. Hybridization was performed in 0.9M NaCl/20mM Tris-HCl/0.01% SDS for 2 h at 46°C (EUB338) and 49°C (STAPHY2). After washing at 48°C (EUB338) and 46°C (STAPHY2) for 15 min in 0.215 M NaCl/20 mM Tris-HCl/0.01% SDS, samples were stained with CD45 and DAPI. Images were acquired using the TissueFAXS software (V6.01.126) and bacteria were detected with StrataQuest (V6.01.167; TissueGnostics GmbH). Epidermis and dermis (<500 µm) were analyzed separately, while blood vessels and hair follicles were excluded.

Co-localization score

Co-localization of bacteria and CD45⁺ cells was calculated with an adapted version of the spatial arrangement analysis-function from DAIME (25). As our input data contained the centroid position of cells, we adjusted the formula by calculating the total number of hits at range r as the cumulative hits at or below r and r-1 (H(r) = H_cumm(r)-H_cumm(r-1)). Hits or misses are replaced by the product of possible dipoles at range r (the circumference of a circle of radius r: 2* π*r) and the total area (At) (T(r)= 2* π*r*At). In addition, population densities (D1 and D2) are calculated using the number instead of the area of one population (N1 or N2) through the total image area (At): D1 = N1/At, D2 = N2/At).

Fluorescence-activated cell sorting

Cutaneous single-cell suspensions were stained for surface markers after collagenase (0.5U/ml) and DNAse I (10mg/ml) digestion (6). All antibodies are indicated in Suppl. Tab. 17. For Smart-seq2 bulk-RNA sequencing 100 cells of both T cells (CD45⁺CD3⁺HLA-DR^-) and mononuclear phagocytes (CD45⁺CD3^-HLA-DR^+high), were sorted into lysis buffer containing RNase inhibitor (1.25 U/ml). All stainings were adapted to isotype-matched controls and gated accordingly.

Smart-seq2 bulk RNA sequencing

cDNA synthesis was performed using the Smart-seq2 protocol (26) and sequenced on an Illumina HiSeq 3000/4000 platform at the 50-bp (base pair) single-read configuration. Read trimming was conducted using ReadTools (v.1.0.0) (27) and mapped using STAR (v2.5.3a) (28). Reads mapped in multiple genomic locations were excluded by SAMtools (v1.4) (29). Read counts were generated using the featureCounts function from the Subread package (v.1.22.1) (30). Gene expression analysis was done with DESeq2 (v.3.22.3) (31). Normalized count (log2) by DESeq2 was used for t-SNE and heatmaps. Pathway enrichment analysis was done by using the KEGGA and GOANA methods, respectively, from the limma package (32).

Co-cultures

Keratinocytes isolated from skin removed in abdominoplasty surgeries using Dispase II (4 U/ml) and PBMCs obtained by Ficoll-Paque reagent were co-cultivated with heat-inactivated S. epidermidis from three aGVHD patients at MOI 0.1, 1 and 10 for 8 h. Staphylococci were grown on a Columbia CNA agar with 5% sheep blood at 37°C and 10% CO₂. Identification to species level was performed by MALDI-TOF MS. From keratinocytes, RNA was isolated using the RNeasy Mini Kit and RT-qPCR was performed: 10 min at 55°C, 1 min at 95°C, 40 cycles of 10 sec at 95°C and 1 min at 60°C. All primers are listed in Suppl. Tab. 18. Stimulated PBMCs were labeled using surface antibodies followed by intracellular cytokine staining (see Suppl. Tab. 17) (33).

Whole-genome sequencing

DNA was extracted using the DNeasy PowerSoil Pro Kit, followed by library preparation (NEBNext® Ultra™ II FS DNA Library Prep Kit) and isolate sequencing on the Illumina MiSeq platform (v3 chemistry, 600 cycles). Sequencing reads were quality trimmed using cutadapt v.2.10 (34) and assembled by SPAdes v.3.14.1 (35). Contigs shorter than 1 kbp were removed with seqtk v.1.3 (github.com/lh3/seqtk). Assembly quality was assessed with ChechM v.1.1.1 (36) and genomes were compared to each other using fastANI v.1.3 (37). Genomes were annotated using prokka v.1.14.5 (38) and protein orthologs were detected using proteinortho v.6.0.24 (39).

Statistical analysis

For normally-distributed data, one-way ANOVA and Bonferroni pairwise comparisons were used, while for non-normally distributed data, Kruskal-Wallis and Dunn’s post hoc tests were applied. Competing risk analysis for GVHD with the competing events death and relapse were calculated using the R-package cmprsk. For patients with GVHD, covariates were compared between mild and severe GVHD-patients by t-, Mann-Whitney- and Chi-squared tests. Mixed linear models were analyzed including all covariates with random factor patient ID in the full patient set. Changes of mean alpha diversity in skin/stool were analyzed by generalized estimating equations. P values smaller than 0.05 were considered significant.

HSCT disrupts the skin and gut microbiome

The effect of allo-HSCT on the cutaneous and gut microbiome as well as skin immunity was assessed in 50 patients up to one year after transplantation (Fig. 1A, Suppl. Tab. 1-4). Compared to healthy volunteers (n=15), we identified HSCT-specific clustering, particularly in the skin (Fig. 1B, Suppl. Fig. 1A-E). While both acute (skin n=10, gut n=6) and chronic GVHD (skin n=4, gut=0) cases were recorded during the one-year observation period, only acute GVHD (aGVHD) patients were considered for this analysis. Alpha diversity of the gut microbiome was disrupted early after HSCT (day of transplantation until d14), which was associated with aGVHD diagnosed by Glucksberg criteria (13) (mean [CI] d100 aGVHD free 3.65 [0.55, 6.74; aGVHD 9.96 [8.17, 11.74]; p=0.0338; Fig. 1C). However, Bray Curtis beta diversity showed a long-lasting shift of both the skin and intestinal microbiome compared to controls 14 days after transplantation (mean [CI] skin control 0.4 [0.28, 0.52] vs. skin HSCT 0.63 [0.57, 0.7]; p=0.0016 and gut control 0.54 [0.37, 0.66] vs. gut HSCT 0.9 [0.8, 1.0]; p=0.011; Suppl. Fig. 3A-B). While alpha diversity measures microbial diversity within a single sample, beta diversity estimates the similarity/ dissimilarity in-between samples.

We next investigated changes in relative abundance of bacteria after HSCT. In skin, commensals such as Corynebacterium, Cutibacterium and Dermacoccus were decreased (Fig. 1D, Suppl. Fig. 2A), while Pseudomonas were increased 14 days after transplant in HSCT recipients (mean [CI] baseline 0.7% [0.37, 1.01]; d14 2.12% [1.12; 3.15]; p=0.012; Suppl. Fig 3C). Further, Staphylococcus spp. showed first a lower abundance at d14 (mean [CI] control 50.66% [37.19, 64.13]; skin aGVHD 24.58 [14.02, 35.13]; p=0.034) and later an expansion at d100 and d365 in skin aGVHD patients (Fig. 1D-E, Suppl. Fig. 2A). Overabundance of staphylococci was not due to a particular ASV (amplicon sequencing variant, Fig. 1F). In stool, loss of butyrate-producing bacteria after conditioning was paralleled by an expansion of Enterococcus, Lactobacillus and Streptococcus (Fig. 1D, Suppl. Fig. 2, 3D-F).

Low alpha diversity and high staphylococci in skin aGVHD

We next investigated correlations between clinical outcomes and microbial signatures, however, only acute GVHD cases were included in the analysis, except for statistical models (Suppl. Tab. 5-12). Alpha diversity in skin was lower after cutaneous aGVHD onset (mean [CI] before onset 7.46 [5.56, 9.28]; after onset 3.84 [2.71, 4.97]; p=0.0242), while no differences were observed in the gut, although this comparison was highly underpowered given the low number of intestinal GVHD patients (n=6, Fig. 2A). Generalized estimated equations (GEE) confirmed these results, while also showing a significant decrease in diversity in stool samples (estimate [CI] skin -4.16 [-6.811, -1.287]; p=0.006 and gut -5.1901 [-7.17, -3.21]; p<0.000; Suppl. Tab. 5-6). As all patients were administered antibiotics prior and during the transplantation (Suppl. Fig 3G), we next investigated the effect of different antimicrobial agents on alpha diversity. GEE revealed a negative effect of Piperacillin/Tazobactam on diversity in skin GVHD (estimate [CI] -4.93 [-8.684, -1.182]; p=0.001; Suppl. Tab. 7), while no effect was observed for stool samples. In addition, competing risk analysis showed no cofounders for GVHD development, including variables such as age, gender, diagnosis prior to HSCT, conditioning regimen as well as number and type of antibiotics (Suppl. Tab. 8-11).

In skin, the degree of alpha diversity reduction before disease onset was predictive of disease severity (mean [CI] mild skin aGVHD 8.66 [6.33, 10.98], severe skin aGVHD 3.93 [2.44, 5.42], p=0.039; Fig. 2B). We again observed an increased relative abundance of staphylococci after aGVHD onset (mean [CI] before onset 39.5% [32.4, 46.59]; after onset 57.06% [43.69, 70.42]; p=0.02; Fig. 2C, Suppl. Tab. 12), which was not due to one particular ASV (Fig. 2D). While a significant increase in staphylococci abundance was observed after onset of mild skin aGVHD (mean [CI] before onset 39.2% [30.66, 47.8]; after onset 58.52% [44.17, 72.87]; p=0.0431), no difference was observed in severe cases (Fig. 2E).

Functional evaluation of S. epidermidis in skin aGVHD

To address our findings of staphylococcal overgrowth in skin aGVHD, we performed FISH targeting bacteria (EUB338) and only Staphylococcus spp. (Staphy2) in aGVHD lesions (n=12). Together with a high infiltration of CD45⁺ cells (Suppl. Fig. 4A), colonization with staphylococci was increased in the epidermis of lesional skin compared to unaffected controls (mean [CI] aGVHD free 12.04 [9.05, 15.03]; skin aGVHD 54.96 [41.47, 68.45]; p=0.012; Fig. 3A, Suppl. Fig. 4B).

Based on these observations, we next performed whole-genome sequencing of Staphylococcus spp. isolates from untreated aGVHD lesions and unaffected skin sites. Cultivation of skin swabs prior to sequencing confirmed increased growth of staphylococci in lesional skin (Fig. 3B). Interestingly, S. aureus only grew in one out of 6 lesions (Suppl. Tab. 14), while S. epidermidis was the only species found in all samples. We therefore performed sequencing of 6 lesional and 6 non-lesional S. epidermidis strains (Suppl. Tab. 15-16) and found an average nucleotide identity of 99% across all samples except one, with a 97% sequence identity (Suppl. Fig. 4C). Around 75% of genes were shared between all isolates (1983 core vs. 717 flexible genes), while most differences were found in-between patients (Fig. 3C). Only few genes were unique to lesional (36) and non-lesional skin sites (35). Co-cultivation experiments of the same S. epidermidis strains with primary human keratinocytes and peripheral blood mononuclear cells (CD14⁺ monocytes and CD3⁺ T cells) confirmed this observation, showing no differences in regulation of antimicrobial peptides and inflammatory cytokines between lesional-and non-lesional isolates (Fig. 3D-E, Suppl. Fig. 4D).

Skin recolonization in HSCT

To investigate bacterial colonization and the cutaneous immune response following HSCT, we combined rRNA-targeted FISH (EUB338) with antibodies for CD45 (Fig. 4A). Interestingly, not only CD45⁺ cells (Fig. 4B), but also bacteria were decreased after conditioning in the epidermis (mean [CI] EUB338 baseline 35.25 [27.71, 42.78] vs. d14 15.51 [11.03, 19.9]; p=0.0005; Fig. 4C). A similar reduction was observed for staphylococci by using a FISH-probe targeting only Staphylococcus spp. (mean [CI] Staphy2 baseline 22.56 [18.39, 26.73] vs. d14 5.37 [3.47, 7.27]; p=0.005; Fig. 4D). However, numbers of CD45⁺ cells, total bacteria and staphylococci were independent of skin aGVHD development (Suppl. Fig. 5A-C).

Further, spatial arrangement analysis adapted from the DAIME software(25) (Suppl. Fig. 6A-C) revealed closer proximity of bacteria and CD45⁺ cells 14 days after transplantation in the epidermis of patients that later developed aGVHD (mean [CI] GVHD free 0.09 [-0.02, 0.19]; aGVHD 1.5 [0.18, 2.83]; p=0.0211, Fig. 4E left), even though both bacteria and leukocytes were reduced in overall numbers at that time (Fig. 4B-C). A similar co-localization in the epidermis, but not in the dermis, was observed for Staphylococcus spp. prior to aGVHD onset (mean [CI] GVHD free 0.11 [0.02, 0.19]; aGVHD 3.01 [-0.41, 6.55]; p= 0.0341; Fig. 4E right, Suppl. Fig. 5D).

Activation of host immunity after conditioning

To assess the consequences of microbial disruption on skin immunity, we evaluated the transcriptome of isolated skin T cells (CD45⁺CD3⁺HLA-DR^-) and mononuclear phagocytes (MNPs, CD45⁺CD3^-HLA-DR^+high) at the 5 study timepoints by bulk-RNA seq (Fig. 5A). Upon unsupervised clustering, transcriptomes of MNPs and T cells were clearly distinguishable, while no time-specific clustering was observed (Suppl. Fig. 7A-B). However, the highest number of differentially regulated genes were observed from Tx until d14 (MNPs:728& 822; T cells: 448& 411, Fig. 5B), suggesting an early influence of conditioning on host immunity. In MNPs, we observed enrichment of cytosolic DNA-sensing and NOD-like receptor signaling pathways at d14, including CASP4, an essential regulator in NLRP3 inflammasome-dependent IL-1 release (40) and phago-lysosomal fusion genes PSME1, P4HB and PRCP (41, 42). Additionally, antigen receptor-signaling was upregulated, including genes PSMB3 and PSMD4 as well as TAPBP, which facilitates peptide loading (43) (Fig. 5C-D). Further, enrichment of type 1 interferon signaling (Fig. 5D) was observed, which promotes dendritic cell activation (44, 45).

In parallel, T cells showed enriched responses to LPS and positive regulation of NF-kappaB signaling. In addition, MYO1G, a gen involved in T cell interaction with antigen-presenting cells (46) was increasingly expressed at d14. At d100, T cell receptor signaling was enriched together with granzymes (Fig. 5C-D). Assessment of receptor-ligand pairs (47) between T cells and MNPs confirmed increased TCR co-stimulation by upregulation of CD28-CD80 and CTLA4-CD80 at d100 and d365 (Fig. 5G). Additionally, receptor-ligand pairs involved in chemotaxis and cell migration were increasingly expressed at d14, including CCR7-CCL21, CXCR3-CXCL9, CCL2-CCR2 and CCR1-CCL5 (Fig. 5G).