Bacterial DNA on the skin surface overrepresents the viable skin microbiome

Abstract

The skin microbiome provides vital contributions to human health. However, the spatial organization and viability of its bacterial components remain unclear. Here, we apply culturing, imaging, and molecular approaches to human and mouse skin samples, and find that the skin surface is colonized by fewer viable bacteria than predicted by bacterial DNA levels. Instead, viable skin-associated bacteria are predominantly located in hair follicles and other cutaneous invaginations. Furthermore, we show that the skin microbiome has a uniquely low fraction of viable bacteria compared to other human microbiome sites, indicating that most bacterial DNA on the skin surface is not associated with viable cells Additionally, a small number of bacterial families dominate each skin site and traditional sequencing methods overestimate both the richness and diversity of the skin microbiome. Finally, we performed an in vivo skin microbiome perturbation-recovery study using human volunteers. Bacterial 16S rRNA gene sequencing revealed that, while the skin microbiome is remarkably stable even in the wake of aggressive perturbation, repopulation of the skin surface is driven by the underlying viable population. Our findings help explain the dynamics of skin microbiome perturbation as bacterial DNA on the skin surface can be transiently perturbed but is replenished by a stable underlying viable population. These results address multiple outstanding questions in skin microbiome biology with significant implications for future efforts to study and manipulate it.

Keywords: Cutibacterium acnes; E. coli; Staphylococcus epidermidis; bacterial viability; human; infectious disease; microbiology; microbiome.

Figure 1.. Bacterial fluorescence in situ hybridization (FISH) staining of human tissue.

(A–C) Scale bar = 20 µm. The bottom-left corner of each diagram shows a schematic of the hair follicle in white and the anatomic location of each image frame in yellow. DAPI staining is shown in blue in all parts. EUB338 hybridization is shown in red for all images. (A) Human tissues stained with the pan-bacterial FISH probe EUB338 show little bacterial presence at the skin surface. (B) Human tissues stained with EUB338 show abundant bacterial signal that is concentrated in hair follicles, pilosebaceous units, and other cutaneous structures. (C) Human tissues stained with a C. acnes-specific FISH probe (in green) demonstrate the same overall spatial organization as those stained with EUB338. (D) Quantification enrichment scores showing the median and interquartile range. Significance was calculated using the Mann–Whitney test. *p≤0.05, **p≤0.01. N = 8 for human follicle, N = 6 for human follicle (C. acnes), N = 6 for human stratum corneum, N = 5 for human stratum corneum (C. acnes) where ‘N’ represents different follicles or stratum corneum sections. The human tissue samples shown here were obtained from adult facial tissues (cheek and forehead).

Figure 1—figure supplement 1.. Fluorescence in situ hybridization (FISH) on stationary-phase skin microbiome bacterial species.

All bacterial species were grown in appropriate conditions to stationary phase. The pan-bacterial FISH probe EUB338 was used. Hybridization is indicated in red.

Figure 2.. Propidium monoazide-droplet digital PCR (PMA-ddPCR) and viability scores for human skin and non-skin microbiomes.

(A) Schematic of the PMA-ddPCR workflow. (B) Sampling scheme showing each skin site that was sampled. Colors indicate site type (sebaceous in blue, moist in green, dry in red). (C) PMA-ddPCR on skin and non-skin microbiome sites shows that the viability score of the skin microbiome is significantly lower than other microbiome sites. ****p≤0.0001 for Student’s t-test on pooled skin and non-skin samples. Four volunteers contributed skin and non-skin microbiome samples. Additional samples were collected from some individuals and represent biological replicates. N = 8 for glabella, N = 6 for retroarticular crease, N = 5 for lower back, hair shaft, nares, and dorsal forearm, N = 3 for antecubital fossa, tongue, saliva, and plaque, N = 2 for popliteal fossa, and N = 1 for human feces. Each human skin sample site consists of samples from four different individuals. Some volunteers were sampled multiple times on different days (at least 2 wk apart). For glabella, one volunteer was sampled four times, one volunteer was sampled two times, and two volunteers were sampled one time. For retroauricular crease, two volunteers were sampled two times, and two volunteers were sampled one time. For lower back, one volunteer was sampled two times and three volunteers were sampled one time. For hair shaft, all samples came from one volunteer. For antecubital fossa, three volunteers were sampled one time. For popliteal fossa, two volunteers were sampled one time. For nares, one volunteer was sampled two times and three volunteers were sampled one time. For dorsal forearm, one volunteer was sampled two times and three volunteers were sampled one time. Tongue, saliva, and plaque all represent one sample from three different individuals. For raw ddPCR counts, see Figure 2—figure supplement 2A and B. (D) PMA-ddPCR on follicle contents and forehead swabs from five individuals. Mean viability score for follicle contents is 0.15 and for forehead is 0.013. All error bars indicate standard deviation.

Figure 2—figure supplement 1.. Propidium monoazide-droplet digital PCR (PMA-ddPCR) and sampling controls.

(A) PMA-ddPCR validation controls on known ratios of heat-killed and exponentially growing E. coli cultures. PMA-ddPCR performed on a population of exponentially growing cells resulted in a viability score of 0.903,. PMA-ddPCR performed on a population of 100% heat-killed cells in exponential growth resulted in a viability score of 0. Populations consisting of 50% (by volume) heat-killed and 50% exponentially growing cells exhibited an average viability score of 0.506. PMA-ddPCR performed on a population of stationary-phase cells resulted in a viability score of 0.965. PMA-ddPCR performed on a population of 100% heat-killed cells in stationary phase resulted in a viability score of 0. Populations consisting of 50% (by volume) heat-killed and 50% stationary-phase cells exhibited an average viability score of 0.621. Mean and standard deviation are shown. (B) Standard curves showing correlation between PMA-ddPCR counts and colony-forming units (CFUs) for four skin microbiome bacterial species when grown to stationary phase. Open data points indicate samples without the use of PMA, and closed data points indicate that PMA was used. Horizontal lines connect paired samples. (C) Standard curve showing correlation between PMA-ddPCR counts and CFU for exponentially growing Staphylococcus epidermidis. 95% confidence interval is shown in the green shaded region. (D) Standard curve generated using S. epidermidis cultures shown in green. Shading represents 95% confidence interval. Open and closed circles represent skin microbiome samples that did (closed circles) or did not (open circles) receive PMA treatment. Dark or light gray shading represents 95% confidence interval for skin microbiome samples. Paired samples are connected to show the downward shift in DNA abundance with the inclusion of PMA. (E) Comparison of the predicted ddPCR counts to measured ddPCR counts based on CFU for samples that were treated with PMA (closed circles) and samples that were not treated with PMA (open circles) (mean for PMA-treated samples is 1.31, mean for untreated samples is 82.2). (F) Comparison of the predicted CFU to measured CFU based on ddPCR for samples that were treated with PMA (closed circles) and samples that were not treated with PMA (open circles) (mean for PMA treated samples is 1.28, mean for untreated samples is 58.5). (G) The effect of including 0.1% Triton X100 in swabbing buffer used for skin microbiome sampling. (H) Comparison of viability scores obtained with the two methods of DNA isolation used. (E-F) Error bars indicate standard deviation.

Figure 2—figure supplement 2.. Copies per 20 µL droplet digital PCR (ddPCR) reaction without (A) and with (B) the use of propidium monoazide (PMA).

Data in (A) and (B) were used to calculate the viability score shown in Figure 2. Error bars represent standard deviation.

Figure 3.. Relative abundance and change in richness and diversity of traditional sequencing compared to PMA-seq.

(A) Relative abundance of all sequenced bacterial taxa at the family level. Paired bars represent data from traditional sequencing (left) and PMA-seq (right). Samples are ordered by increasing richness in traditional sequencing. Labels below each pair of bars indicate each sample’s donor, replicate, and site (e.g., HV1.1 RAC indicates Healthy Volunteer 1, replicate sample 1, retroauricular crease). HS: hair shaft; RAC: retroauricular crease; VF: volar forearm; PF: popliteal fossa; TW: toe web; AF: antecubital fossa. *Relative abundance data for Staphylococcaceae was determined using forward-read sequencing information only. Samples with fewer sequencing reads than PBS controls are not displayed. All identified bacterial taxa with corresponding colors can be found in Figure 3—figure supplement 1. (B) The richness changes between traditional sequencing (R_trad) and PMA-seq (R_PMA) are demonstrated by plotting the change in richness (∆R) against R_trad. Colors represent different site types and shapes represent different sample sites. The shaded gray region represents the 95% confidence interval for the linear regression. (C) The Shannon diversity changes between traditional sequencing (H_trad) and PMA-seq (H_PMA) are demonstrated by plotting the change in diversity (∆H) against H_trad.

Figure 3—figure supplement 1.. Full list of identified taxa with corresponding colors.

(A) Full list of identified bacterial groups shown in Figure 3.

Figure 4.. Propidium monoazide (PMA) index and change in relative abundance between traditional and PMA-seq.

(A) The PMA index for each bacterial taxon that was present in at least four samples is shown here as an average between samples of the same sample site shown in Figure 3. Color indicates PMA index value. Saturation indicates confidence (sigma) in the PMA index value and was calculated using the standard deviation of PMA index across the samples that went into that pixel. Bacterial taxa are ordered by decreasing overall relative abundance. Each square represents the average of at least four samples taken from different individuals. PMA index is calculated by comparing the relative abundance of a given taxon as measured by PMA-seq (A_PMA) to the sum of the relative abundance for that taxon in both traditional sequencing (A_trad) and PMA-seq. (B) Relative abundance at each body site for the top three most abundant (overall) bacterial taxa as assessed by traditional sequencing and PMA-seq. Colors of bars correspond to colors in Figure 2F. Error bars represent standard deviation. *Relative abundance data for Staphylococcaceae was determined using forward-read sequencing information only.

Figure 4—figure supplement 1.. Viability score for three skin sites using lysostaphin and Staphylococcus-specific PCR primers.

(A) The three skin sites with the most abundant bacterial DNA are shown. Half of each sample was treated with lysostaphin prior to DNA isolation to assess how the viability score would change. Droplet digital PCR (ddPCR) was performed on each sample using both 16S primers (white bars) and Staphylococcus-specific (black bars) primers. The dashed line indicates the average viability score of non-skin microbiome sites (0.66). N = 3 for all. Error bars represent standard devation.

Figure 5.. Bacterial fluorescence in situ hybridization (FISH) staining of mouse tissue (A–D) and comparison of mouse viability scores and human viability scores.

(A). Tissues from a K14-H2B-GFP mouse stained with EUB338 show abundant bacterial signal in hair follicles but not on the skin surface. (B) Tissues from SKH1-Hrhr Elite nude mice also show bacterial presence concentrated to cutaneous structures and not at the skin surface. (C) E. coli applied to C57BL/6 mouse tissue was stained with either EUB338 (in red) or its complementary strand control probe NONEUB338 (in yellow). (D) Quantification enrichment scores showing the median and interquartile range. Significance was calculated using the Mann–Whitney test. *p≤0.05, **p≤0.01, N = 6 for hairy mouse follicle, nude mouse follicle, and hairy mouse stratum corneum, N = 5 nude mouse stratum corneum. (E) The PMA-ddPCR-based viability scores for mouse skin microbiomes are much lower than for mouse fecal microbiomes (0.66 and 0.98, respectively). These viability scores for mouse sites are very similar to those for humans (0.066 and 0.045 for skin microbiomes, 0.98 and 0.66 for fecal microbiomes). Error bars represent standard deviation.

Figure 6.. Skin microbiome perturbation and recovery.

(A) Bacterial relative abundance in each individual over the 48 hr following perturbation. 0 hr represents baseline, pre-perturbation community. Whether a sample was treated with propidium monoazide (PMA) is indicated by (-) and (+). (B) Quantification of bacterial DNA recovery over the 48 hr following perturbation. DNA was quantified using droplet digital PCR (ddPCR). (C) Bray–Curtis dissimilarity of each individual over the 48 hr following perturbation. Red data points are comparing PMA-treated samples to the PMA-treated baseline sample. Blue data points are comparing PMA-untreated samples to the PMA-treated baseline sample. Dashed vertical line indicates the point of perturbation.

Figure 6—figure supplement 1.. List of bacteria identified in perturbation recovery.

The bacterial groups listed here correspond to the entire sequencing dataset shown in Figure 6. Bacteria are listed in order of relative abundance.

Figure 6—figure supplement 2.. Skin microbiome perturbation and recovery.

(A) Bray–Curtis dissimilarity of each individual over the 48 hr following perturbation. Red data points are comparing propidium monoazide (PMA)-treated samples to the PMA-untreated baseline sample. Blue data points are comparing PMA-untreated samples to the PMA-untreated baseline sample. Dashed vertical line indicates the point of perturbation. (B) Amplicon sequencing variants (ASV)-level relative abundance converted into binary dataset in which each ASV present at greater than 1% appears in yellow and each absent ASV appears in black. Columns are ASVs and rows are PMA-treated or PMA-untreated perturbation recovery samples. ASVs indicated in red follow a specific pattern of repopulation in which the ASV is present in the live cells (+PMA) at T = 0, disappears from the live cells after perturbation (T = 3), and then recovers in the live population at either T = 24 or T = 48.

Figure 6—figure supplement 3.. Contamination removal performed on 600 nt sequencing data.

The Shannon diversity changes between traditional sequencing (H_trad) and PMA-seq (H_PMA) are demonstrated by plotting the change in diversity (∆H) against H_trad (A–C). The richness changes between traditional sequencing (R_trad) and PMA-seq (R_PMA) are demonstrated by plotting the change in richness (∆R) against R_trad (D–F). Analysis with no decontamination removal (A, D), decontamination removal using the Decontam program in R with a 0.1 threshold (B, E), and decontamination removal using the Decontam program in R with a 0.2 threshold (C, F) are all shown. Symbols in red represent samples associated with the perturbation recovery data (Figure 6). Symbols in blue represent PBS DNA isolation controls. Symbols in black represent follicle contents obtained from single follicles. The shaded red region represents the 95% confidence interval for the linear regression performed using only the perturbation recovery samples (red symbols).