Skin Microbiome Surveys Are Strongly Influenced by Experimental Design

Abstract

Culture-independent studies to characterize skin microbiota are increasingly common, due in part to affordable and accessible sequencing and analysis platforms. Compared to culture-based techniques, DNA sequencing of the bacterial 16S ribosomal RNA (rRNA) gene or whole metagenome shotgun (WMS) sequencing provides more precise microbial community characterizations. Most widely used protocols were developed to characterize microbiota of other habitats (i.e., gastrointestinal) and have not been systematically compared for their utility in skin microbiome surveys. Here we establish a resource for the cutaneous research community to guide experimental design in characterizing skin microbiota. We compare two widely sequenced regions of the 16S rRNA gene to WMS sequencing for recapitulating skin microbiome community composition, diversity, and genetic functional enrichment. We show that WMS sequencing most accurately recapitulates microbial communities, but sequencing of hypervariable regions 1-3 of the 16S rRNA gene provides highly similar results. Sequencing of hypervariable region 4 poorly captures skin commensal microbiota, especially Propionibacterium. WMS sequencing, which is resource and cost intensive, provides evidence of a community's functional potential; however, metagenome predictions based on 16S rRNA sequence tags closely approximate WMS genetic functional profiles. This study highlights the importance of experimental design for downstream results in skin microbiome surveys.

Fig 1. Study design for analyzing cutaneous bacterial communities

(A) Seven skin sites were sampled from a healthy human cohort of nine volunteers (B). DNA was isolated from cutaneous swabs and sequenced for downstream bioinformatics analyses. (C) Schematic illustrating the primers used for the two different targeted hyper-variable regions on the 16S rRNA gene.

Fig 2. Taxonomic profiles of cutaneous bacterial communities vary

(A) Heatmap depicting relative abundances of the 20 bacterial species in the MCC. Rows display bacterial species and columns denote sequencing technique used for analysis. “Actual” refers to the expected abundance based on the community composition. Dendrogram (x-axis) clusters each sequencing type by similar taxonomic profiles. (B) Pie charts depicting the mean relative abundances of the top 15 taxa in the cutaneous samples. The innermost circle represents the V4 samples, the middle circle the V1-V3 samples, and the outermost WMS samples

Fig 3. Species level classification of Staphylococcus sequences

(A) Barplot of the total number of OTUs identified in cutaneous samples by V1-V3 and V4 sequencing, highlighting the number of OTUs named at the species level (B) Mean relative abundance of Staphylococcus species identified by 16S tag sequencing of skin samples at greater than 1% (C-E) Relative abundance of Staphylococcus sequences able to be classified at the species level. (C) Staphylococcus species in the WMS dataset were classified using MetaPhlAn. (D,E) V1-V3 and V4 species level classifications were determined by pplacer. Pie charts depict the percentage of sequences classified as Staphylococcus at the genus level that were further classified at the species level.

Fig 4. 16S predictions differ from whole metagenome functional enrichment

(A) Heatmap depicting the log2 fold change of statistically significantly different KEGG categories (FDR corrected paired Wilcoxon test, p < 0.05) between 16S and WMS functional profiles. Purple depicts enrichment in WMS samples and green depicts enrichment in the 16S samples. Each column represents a different sample and each row a KEGG pathway. The colors above the columns indicate sequencing of the V1-V3 (gray) or V4 (black) region of the 16S rRNA gene. (B) Box plots depicting mean relative abundances of significantly different KEGG pathways (FDR corrected paired Wilcoxon test, p < 0.05) between 16S and WMS functional profiles.

Fig 5. Cutaneous taxonomic and functional diversity trends depend on sequencing method

Shannon diversity of (A) taxonomic and (B) functional profiles for each sequencing technique is presented by site microenvironment, with asterisks (*) indicating significance of p < 0.05 using the Kruskal-Wallis and multiple comparison post hoc test. Boxplots were calculated using the ggplot2 R package. (C) Procrustes analysis, revealing congruence between NMDS ordinations of the V1-V3 (target) and V4 (rotated) Bray Curtis dissimilarity matrices. Circles indicate V4 samples and diamonds indicate V1-V3 samples, with matched samples connected by a line. Shorter lines reflect greater clustering similarity.