Integrated genomic and functional analyses of human skin-associated Staphylococcus reveal extensive inter- and intra-species diversity

Abstract

Human skin is stably colonized by a distinct microbiota that functions together with epidermal cells to maintain a protective physical barrier. Staphylococcus, a prominent genus of the skin microbiota, participates in colonization resistance, tissue repair, and host immune regulation in strain-specific manners. To unlock the potential of engineering skin microbial communities, we aim to characterize the diversity of this genus within the context of the skin environment. We reanalyzed an extant 16S rRNA amplicon dataset obtained from distinct body sites of healthy volunteers, providing a detailed biogeographic depiction of staphylococcal species that colonize our skin. S. epidermidis, S. capitis, and S. hominis were the most abundant staphylococcal species present in all volunteers and were detected at all body sites. Pan-genome analysis of isolates from these three species revealed that the genus-core was dominated by central metabolism genes. Species-restricted-core genes encoded known host colonization functions. The majority (~68%) of genes were detected only in a fraction of isolate genomes, underscoring the immense strain-specific gene diversity. Conspecific genomes grouped into phylogenetic clades, exhibiting body site preference. Each clade was enriched for distinct gene sets that are potentially involved in site tropism. Finally, we conducted gene expression studies of select isolates showing variable growth phenotypes in skin-like medium. In vitro expression revealed extensive intra- and inter-species gene expression variation, substantially expanding the functional diversification within each species. Our study provides an important resource for future ecological and translational studies to examine the role of shared and strain-specific staphylococcal genes within the skin environment.

Keywords: 16S rRNA amplicon sequencing; growth curve; pan-genomics; skin microbiota; staphylococci.

Fig. 1.

Diversity of staphylococcal species present on healthy human skin. (A) Barplots display the relative abundance of staphylococcal species at 14 body sites. Each color represents a distinct species as shown in the legend. Each bar represents one Subject/Healthy Volunteer (HV). Empty bars represent missing data. (B) Beta diversity analysis of staphylococcal communities present at different body sites (displayed as colored dots) using principal-coordinate analysis (PCoA) plot based on Bray–Curtis dissimilarity (PERMANOVA; R² = 0.26804, P value < 0.0001). First two coordinates are shown accounting for 40.62% of the total variance. Individual staphylococcal ASVs driving the largest separation of body sites are shown as black arrows and are labeled by the corresponding species. Note that the body site legend is shared between figures B and C. (C) Relative abundance and prevalence of staphylococcal species detected in our dataset. Each dot represents proportion of a staphylococcal species relative to all other staphylococcal species in a sample. Dots are colored by body sites as shown in the shared legend. Percent prevalence threshold for a species to be considered part of the core staphylococcal community is shown in the upper gray bar (see text for details). Refer to SI Appendix, Fig. S1 for body site details.

Fig. 2.

Genus-level pan-genome summary to quantify gene sharing between species. (A) Histogram of core, soft-core, shell, and cloud genes of the genus-pangenome. Pie chart shows the number of genes in each category as defined in the legend; percent of genes in each category are displayed in brackets. (B) Venn diagram displaying the distribution of 7,744 genus pan-genome genes into seven subsets based on their presence in one, two or three species. Number of total genes in each subset are shown outside the bracket and the core genes in each sub-set (genes detected in ≥95% genomes present in the sub-set) are shown inside brackets. (C) Clustering of all 126 genomes using Euclidean distance and Ward hierarchical clustering based on the presence/absence patterns of 7,744 genes detected in the pan-genome. Sidebars display species, volunteer, and body site of isolation for each genome. (D) Pairwise Jaccard index between genomes displaying the proportion of shared gene content. Note: Distance of 1.0 depicts complete gene overlap, with lower numbers representing lesser degree of gene sharing. Hierarchical clustering of genomes and the displayed sidebars are the same as that shown in figure C. Note that the legend is shared between figures C and D. Refer to SI Appendix, Fig. S1 for body site details.

Fig. 3.

Relationship between phylogeny and pan-genome in conspecific staphylococcal genomes. (A–C) Tanglegrams depicting the relationship between the phylogenetic tree based on core genes sequence alignment (Left) and hierarchical clustering (binary distance; average linkage) of the pan-genome presence/absence matrix (Right). Clades A and B, as determined by rooting each species tree with an outgroup using a genus tree (SI Appendix, Fig. S9), are shown in purple and orange color, respectively. Colored lines are used to connect the same isolate on two trees of the tanglegram. Healthy volunteers and body sites of isolation for isolates are displayed as colored strips within the tanglegram. (D) Correlation between phylogenetic distance and pan-genome distance for all pair-wise genome comparisons within a species. Each dot represents one pair-wise comparison. Dot colors represent the clades to which the genomes being compared belong, with yellow color dots representing clade A versus clade B comparison. Refer to SI Appendix, Fig. S1 for body site details.

Fig. 4.

Kinetic growth analyses of select staphylococcal isolates in BHI-YE and artificial skin media, ES and ESL. (A–C) Growth curves of select S. epidermidis, S. capitis, and S. hominis isolates grown independently in three culture media, as measured by absorbance at 600 nm over 24 h. Each growth curve is an average of four independent biological replicates; error bars for each curve are shown. Isolates marked with an asterisk (*) were selected for RNA-Seq analyses. Purple and orange color of facets indicate the phylogenetic clade to which each isolate belongs. Growth curve color represents the growth medium. (D–F) The AUC was used to quantify the growth of each isolate shown in A–C. Boxplots depict combined AUC values for all isolates of a species in each growth medium, as depicted by the color of the boxplot. The center black line within each boxplot represents the median value, with edges showing the first and third quartiles. Individual AUC values are shown as dots colored by clade of the isolate. Note that the legend is shared between all figures.

Fig. 5.

RNA-Seq analysis of nine distinct staphylococcal isolates in BHI-YE and in artificial skin media, ES, and ESL (N = 3 biological replicates). (A) Principal-component variance stabilizing transformation-normalized reads from all RNA-seq samples with 1,647 genus-core genes. Colors represent different growth media, and shapes represents staphylococcal species. (B) A stacked bar plot representing the proportion of genus-core, species-restricted-core, and accessory within the top 10 down- and up-regulated genes in each isolate in the ES medium relative to BHI-YE. (C) A discrete heatmap of log₂-fold changes of all the genes that were differentially expressed (≥ or ≤ twofold change, adjusted P value < 0.05) in ES relative to BHI-YE. Columns are individual isolates and rows are DEGs. Genes up-regulated in the ES medium are shown in orange, down-regulated genes are in purple. Light gray cells represent genes that were encoded within an isolate genome, but not differentially expressed in that isolate. Dark gray cells denote genes that were absent in the isolate genome. The left sidebar depicts the category to which each gene belongs. Top and side dendrograms were generated by unsupervised clustering of the expression data (Euclidean distance; Ward). (D) KEGG pathway enrichment analysis of DEGs in ES medium relative to BHI-YE in each isolate. The X axis and Y axis show the isolate and pathway, respectively. Bubble size shows GeneRatio and color indicates adjusted P value.