Diverse signatures of skin microbe stability at the species level

We then sought to identify if specific species or community types are intrinsically unstable in the skin niche, or if increased variability is a phenomenon where many species flux simultaneously in a community. For example, if a species is poorly adapted to a skin niche, relative abundances could vary significantly over time with growth competition. Or, in antagonistic interspecies relationships, cyclic changes in relative abundance could also occur. Community composition, describing the relative abundances of bacteria, fungi, and DNA viruses across samples confirmed that skin communities are largely stable over time (Figure 2A, Figure S2A). Of communities that tended to be less stable, we asked whether certain species were more likely to contribute to instability. Correlation analysis of relative abundances of skin species between timepoints in these sites showed that ‘high flux’ species, those that fluctuated more over time, differed by site and individual (Figure 2B). For example, in the dry hypothenar palm (Hp), we observed that Propionibacterium acnes was a high-flux organism for multiple individuals. In the moist popliteal crease, Polyomavirus HPyV7 was a high flux species for one individual, Polyomavirus HPyV6 for another, and Staphylococcus epidermidis for yet another. In the high-instability plantar heel, abundances of S. hominis, S. epidermidis, and Alphapapillomavirus changed significantly over time. While these species may flux in certain individuals and sites, most skin microbes appeared intrinsically stable. Thus, a complex interaction of host status and environmental perturbations likely influences these local changes in skin communities.

Open in a separate window

Figure 2

Community Stability is Defined by Different Species Characteristics, Depending on Site and Individual

For all panels, sites are color coded by site characteristic. (A) Relative abundances of the most common skin bacteria, fungi, and viruses shown for 3 representative individuals. Full set of taxonomic classifications is shown in Figure S3. T1, T2, and T3 indicate order in time series. (B) Correlating taxon abundances for ‘Long’ vs ‘Short’ time series for select sites from (A) that showed higher variability. Different shapes represent different individuals; different colors represent kingdom of species, and size of point indicates species standard deviation between samples. Select abundant taxa (five most highly variable) are numbered as indicated. Partial Spearman correlations adjusting for repeated measures with associated P-values are in the upper right. (C) Phage and host dynamics at select sites. For Propionibacterium acnes and Staphylococcus, relative abundances of phage versus bacteria are plotted for each sample. A linear model was fitted to the data with shading showing the 95% confidence interval for the goodness of fit. To quantify correlation or anticorrelation, partial Spearman correlation coefficients and associated adjusted P-values were calculated. (D) The relationship between taxon abundance and variance, measured as standard deviation, is shown for taxa > 0.25% relative abundance in the population. A generalized additive model, showing a second order relationship, was fitted to the data to generate correlation coefficients (upper left). All P-values 0.05.

See also Figure S2.

We also examined phage-host relationships as a potential interspecies source of variability in microbial communities. While mixed phage-bacterial communities were typically stable (Figure 2A, Figure S2A), we observed instances of large-scale shifts in phage-bacteria ratios that could reflect a perturbation event, e.g., HV12. We investigated the relative abundances of the most common phages in skin (Figure 2C, Figure S2B) to search for relationships that might provide clues about bacteria-phage interactions. Both P. acnes and Propionibacterium phage are abundant in sebaceous sites, and we observed a strong anti-correlation in sebaceous communities that contain both P. acnes and its phage (Figure 2C). While the compositional (versus absolute) nature of the data means that taxon abundances are not perfectly independent, this anti-correlation together with the observed phage-host dynamics over time suggests antagonism.

We also observed multiple populations of P. acnes that had a neutral or positive correlation with Propioinibacterium phage abundance (Figure S2B, points along x-axis), potentially reflecting variability of P. acnes strains in skin communities in which some strains may be interacting with phages and others may not (Liu et al., 2015). Similar findings were also evident in dry and moist sites, most notably in the feet in which a positive correlation was observed. On the feet, Propionibacterium phages may interact with different Propionibacterium species (e.g., humerusii), or they may exist in a lysogenic versus a primarily lytic state in sebaceous sites. Relative abundances of Staphylococcus and Streptococcus phages were largely positively correlated with Staphylococcus and Streptococcus, suggesting that these may be prophages.

Finally, we examined species abundance as a predictor of stability over time (Figure 2D). We observed a second order, power-law relationship between taxon mean relative abundance across time and variance, excluding transient species (taxa that disappear from one timepoint to the next)(Ma 2015). Both low and high abundance species had low variance, and the largest variance was observed in intermediate abundance species. This second order relationship was present whether site or individual was used to track species, supporting that this intrinsic property of abundance is not determined solely by species identity. This suggests a broader ecological context in which high abundance species can be considered ‘fixed’ in the community, and a limited assortment of rare (but not transient) species can also stably persist. Overall, dominant community characteristics such as diversity and composition are largely stable over time, where primary determinants of stability are site and individual-specificity, with different species contributing to variation in the community over time.

Transience and eukaryotic dsDNA viral stability

Because transient species also likely contribute to community variations over time, we next examined the fraction of the skin community that was transient, i.e. species that were below the detectable range between timepoints. Because technical limitations, including sampling depth, sequencing depth, and amplification bias, can influence detection particularly at low relative abundances, we measured transience at three abundance ranges (<0.1%, ≥0.1% to ≤ 1%, >1%) and required a species to occur in that range in at least one timepoint (Figure 3A). Of very low abundance organisms (<0.1% relative abundance), a significantly higher fraction consisted of transient organisms from all kingdom classes except bacteria, yet this could be partially attributed to reduced robustness in species classification of fungi and viruses at this abundance threshold. For higher abundance ranges, the proportion of stable species significantly exceeded the number of transient organisms for bacteria, fungi, and phages, suggesting that at higher abundances, species from these kingdoms were more likely to be stable residents.

Open in a separate window

Figure 3

Assessment of Transient vs. Non-transient Species Across Time by Kingdom, and Eukaryotic Viral Stability in the Skin

For panels in (A and B), sites are color coded by site characteristic. (A) Fraction of species that are transient or not transient for each complete longitudinal set of sites per individual. Transience, defined as present in only one of three time points, is measured at three ranges of abundances. All comparisons were adjusted P-value < 0.05 between transient/non-transient except for Bacteria, <0.1%, which is marked not significant (N.S.). (B) Relative abundances of the most common predicted eukaryotic viruses are shown for 3 representative individuals. Full set of classifications is shown in Figure S3. T1, T2, T3 indicate order in time series. (C) Predicted eukaryotic viruses tend to be individual-specific, with significantly more sharedness, calculated with the Yue-Clayton theta coefficient, within an individual than between individuals and with little site-specificity. P-value, Wilcoxon rank-sum test shown in left corner. (D) Coreness of predicted eukaryotic viruses. Box plots and overlaid points show the proportion of samples within an individual in which a virus is observed, with a cutoff of >50% occurrence. For lineage, color code is based on panel (B) legend. Right panel shows number of individuals in which a eukaryotic virus is observed. (E) Relative abundance of eukaryotic viruses arranged by site with each individual’s value overlaid and colored by site characteristic. Thick black line indicates mean relative abundance of total (eukayotic and bacterial) DNA viruses at that site.

See also Figure S3.

Noteworthy, we observed that double stranded eukaryotic DNA viruses had a significantly greater fraction of species that were transient at the higher abundance classes (0.1–1% and >1%), in striking contrast to the stable residences of bacteria, fungi, and phages (Figure 3A). This is consistent with previous studies where viral communities were more unstable than the whole metagenome (Hannigan et al., 2015). Types of eukaryotic viruses identified in skin were relatively few (<two dozen, Figure S3A), with large fluctuations in number of species observed in some individuals. Taxonomic characterization of these viruses showed three major classes of communities (Figure 3B, Figure S3B). We observed a ‘uniform’ phenotype, in which the viral communities at all sites were predominantly a single viral type (Figure 3B, e.g., HV08). In an ‘intermediate’ phenotype, we observed a dominant species in a large number of sites but high diversity at other sites (Figure 3B, e.g., HV05). In a ‘variable’ phenotype, viral communities tended to differ markedly from site to site and over time (Figure 3B, e.g., HV12) at low abundances. Viruses’ lack of site specificity was reflected in comparisons using the theta index (Figure 3C, e.g., HV12), where eukaryotic viral communities were highly dissimilar between sites. Moreover, individuals tended to harbor a unique virome (Figure 3C), with intrapersonal similarity significantly exceeding interpersonal similarity.

Unlike bacterial and fungal species like P. acnes and Malassezia globosa, which are nearly universal both between and within individuals, we observed little evidence of a core virome in our cohort, although human polyomaviruses and papillomaviruses were identified in ~half of our subjects (Figure 3D, top). Generally, the most abundant viruses in our cohort were specific to single individuals, including retroviruses, which may be stably incorporated into the host human genome. By contrast, we observed the presence of a strong core phageome (Figure 3D, bottom), in which Propionibacterium, Staphylococcus, and Streptococcus phages were nearly universal between and within individuals. This underscores differences in host specificity in which phages are likely obligate partners with a community’s abundant bacterial species and thus have site specificity while eukaryotic viruses do not (Figure 3C).. As viral databases improve and expand, we anticipate that strain resolution of phages and eukaryotic viruses will further refine these interkingdom interactions.

Because eukaryotic viruses inhabiting the human body presumably infect human cells, we inquired if different skin sites had different abundances of eukaryotic viruses. For example, the thick skin of the foot might harbor fewer detectable viruses than a moist site, whose thinner layers may be more susceptible to skin barrier breaks. While shotgun sequencing methodologies prevent estimation of true abundances, we can gain some insight from the proportion of the community comprised of viral sequences. A significant increase in detectable viral DNA could signify a viral bloom relative to the bacterial and fungal community. Note that we observed no consistent trend with amount of human DNA recovered (data not shown), which could have skewed viruses recovered with higher or lower numbers of active human cells obtained in collection.

We created a ranked distribution of eukaryotic viral relative abundances stratified by skin site (Figure 3E). Certain sites, including alar crease, toeweb, toenail, retroauricular crease, and back, had a large fraction of samples with no viral representation and lower relative abundance than other sites. Sites with modestly higher viral relative abundances were moist skin creases like the popliteal and antecubital fossae (as well as adjacent sites on the forearm and palm), cheek, and occiput. But we observed little uniformity in what viruses tended to be higher in relative abundance at these sites; in most individuals, higher viral relative abundances could be attributed to a diversity of polyomaviruses or papillomaviruses. Interestingly, the composition and presence of eukaryotic viruses in skin communities appeared less dependent on physiologic characteristics that strongly shape bacterial and fungal communities; these dynamics may depend more on host status and topography than nutrient availability.

Individual-specific signatures are shared across time and space

As overall stability is likely an individual trait, with the degree of stability and involved species varying between individuals, we next investigated whether individuals retain unique taxonomic signatures that would allow them to be distinguished from other individuals at different points in time. Recent research has suggested that individuals might possess distinguishing microbial fingerprints that are stable over time (Franzosa et al., 2015). Here, we used a supervised random forests algorithm (Liaw and Wiener, 2002) for each timepoint to identify key taxa that could differentiate individuals, and to identify if these taxa remained the same at each timepoint. Comparing classifier strength (mean decrease in accuracy) over the different timepoints, we observed that, surprisingly, a relatively small number of species retained their discriminatory power over time (Figure 4). Low abundance (<5% relative abundance) organisms were typically those with the highest resolving power between individuals, for example Merkel cell polyomavirus, Gardnerella vaginalis, and Acheta domestica densovirus (top half of Figure 4, right panel). These organisms also tended to be stable in relative abundance over time with low coefficients of variation (Figure 4, right), and were present at most sites within an individual across time (Figure 4, center). Species with minimal discriminatory power between individuals (bottom half of Figure 4) were characterized as lacking ability to discriminate individuals at one or more timepoints, having low prevalence across body sites, and having high variance in their relative abundances over time. The limited ability of species to discriminate individuals over time has been previously described by (Franzosa et al., 2015) who found strain-level codes were more robust to temporal variability.

Open in a separate window

Figure 4

Individual-specific Signature Taxa are Shared across Time, Sites, and are Typically Low Abundance and Highly Stable

Left, variable importance plot of supervised random forest algorithm output. Top 30 taxa contributing to model are shown on the y-axis ranked by importance, quantified by the mean decrease in accuracy; i.e., degree to which inclusion of this predictor in model reduces classification error. Black line indicates the mean decrease in accuracy across the three timepoints; colored points indicate the score for each timepoint. Center, for each individual, proportion of body sites for which each taxa is present, averaged across timepoints. Right, variation in mean relative abundance across sites and timepoints (color) with coefficient of variation.

Individuals have distinct microbial SNV signatures that are stable over time

Shotgun metagenomic data empowers unprecedented resolution of microbial communities at the strain and SNV level. To explore stability at this high resolution, we focused on the prevalent skin bacterial commensals Propionibacterium acnes and Staphylococcus epidermidis, which have dozens of available genome sequences (Conlan et al., 2012; Tomida et al., 2013). Previously, by mapping shotgun metagenomic reads to this reference set of phylogenetically diverse genomes, we identified that skin sites harbor complex subspecies variation for both P. acnes and S. epidermidis. Strain heterogeneity could be unique to an individual (P. acnes) or specific to skin site (S. epidermidis) (Oh et al., 2014).

Here, we observed that P. acnes strains are remarkably stable over time across body sites (Figures 5A and S4A). We quantified this stability with the Yue-Clayton theta similarity index (Figure 5B), taking into account both strain presence/absence and relative abundance. Temporal stability, short- or long-term, surpassed the similarity between individuals, indicating that P. acnes stability likely derives from the maintenance of an individual’s strains over time and less from the acquisition of new strains from the environment or other individuals.

Open in a separate window

Figure 5

Individual Specific Strain and SNV Signatures are Stable Over Time

(A) Dendrogram of P. acnes strain genome similarity based on core SNVs. (B) P. acnes strain relative abundance plots of 3 representative individuals’ manubrium, colors as in (A). Full set of strain classifications is shown in Figure S6B. (C) Boxplots of Yue-Clayton (left) and Jaccard (right) theta indices indicate similarity between strains (all body sites) and SNVs (manubrium and back) of P. acnes in a time series (θ = 1 is identical). ***P < .001 value, Wilcoxon rank-sum test. (D) Rarefaction curves demonstrate core SNV accumulation with read subsampling for manubrium sites. For remaining panels, SNVs are reported for samples subsampled to 1 million reads. Colors correspond to individuals as shown in Figure S1B. (D) Number of SNVs that are mono, di, and triallelic. T1, T2, and T3 indicate order in time series. HV is healthy volunteer. See also Figure S4.

To validate this hypothesis and to assess within-subject retention of strains, we also examined shared SNVs in the P. acnes genomes across the longitudinal samples. To power this analysis, we focused on two sebaceous sites, the manubrium and back, which have high sequencing depth and P. acnes abundance (Figure S4C, S4D). For these sites, we used SNVs specific to the P. acnes core genome (2,248,676 bps region of the genome shared between all 78 sequenced strains) to ensure even representation of SNVs in every metagenomic sample. We identified 83,081 variant positions in the P. acnes core or ~24,000 SNVs per sample after filtering variants for an allele frequency > 1% and > 4 read depth. To test if P. acnes sequencing depth was sufficient to identify the majority of variants, we generated rarefaction curves of SNVs discovered over increasing read depths (Figure 5C). We found that 1 million reads, 40X coverage of the P. acnes core, was sufficient for variant discovery.

The major advantage of an SNV-based approach is that given sufficient sequencing depth, temporal stability and genetic diversity can be estimated in the absence of a large number of sequenced reference strains. By comparing the sharedness of SNVs between timepoints or individuals using the Jaccard index, which measures similarity based on presence/absence of features, we found the stability of P. acnes SNVs mirrors that of P. acnes strains over time (Figure 5B). We found that regardless of duration, an individual shares significantly more SNVs with themself over time than with other individuals (P-value<0.001).

Finally, polyallelism can reflect genetic heterogeneity, i.e., the presence of multiple strains, in a community. While the number of diallelic sites does not scale linearly with the number of strains in a sample, low levels are indicative of a mono-colonized population. We focused on diallelic states, as triallelism was extremely rare (<1.0% in the population). We calculated the cumulative distribution of alternate alleles derived from the shotgun metagenomic reads as a function of distance from a defined reference genome (Figure 5D, Figure S4E). Through these analyses, we identified a putatively monocolonized individual with strikingly low diallelism (1.0% of sites), in contrast to other individuals, e.g., HV02 and HV09, in which 97.4% and 25.6% of alternate sites have reads reliably mapping to both reference and an alternate allele indicating the presence of multiple strains. For future disease studies, fluctuations in pollyallelism, for example, a dramatic decrease in the number of diallelic positions could indicate emergence of a dominant pathogenic strain.

P. acnes pangenome maximized across a multi-phyletic community

Since skin sites stably maintain the same P. acnes strains over time, we wanted to explore the community’s full gene content to generate hypotheses about evolutionary forces shaping community drift, resilience, and stability on a functional level. A species’ total functional repertoire is the ‘pangenome’, composed of core (conserved) and non-core (absent in at least one strain) genes. The P. acnes pangenome is composed of 3,774 non-redundant gene clusters, of which 1,685 were core (Figure S5A–B). Each additional genome adds 3 novel genes to the total (Figure S5C) (Tomida et al., 2013), implying that the majority of P. acnes functional capacity is captured within these 78 reference genomes. We mapped our reads to this P. acnes pangenome database, requiring 1 million reads for adequate sequencing depth (Figure S5D). Between 82.6 (3,117) and 99.9% (3,771) of the known P. acnes pangenome (Figure 6A) was represented in healthy individuals. 70% of samples had > 95% (3,585) of the pangenome.

Open in a separate window

Figure 6

P. acnes Pan-genome Reaches Functional Saturation with Distinct Strain Combinations

(A) Boxplots show number of genes present at > 40% coverage across sites of individuals at 3 time points. Black dashed line indicates the 3,774 genes in P. acnes pangenome. Blue dashed line indicates the 3,005 genes present in every individual at every time point in at least 50% of body sites. Sites with < 1 million P. acnes reads were excluded. (B) Relative abundance plots of P. acnes strains across all body sites, color-coded by site characteristic. Strain colors defined in Figure 5A. T1, T2, and T3 indicate order in time series. (C) Boxplots of Jaccard theta indices indicate stability of gene presence over time (θ = 1 is identical). **P <0 .01, ***P < 0.001, Wilcoxon rank-sum test. (D) Boxplots of Yue-Clayton theta indices indicate the stability of gene copy numbers over time (θ = 1 means identical). (E) Pie chart indicates distribution of P. acnes genes between those present in all individuals “core” and those absent from some individuals “noncore”. Majority of core and noncore genes are pathway “unknown” when compared to a KEGG database. (F) Distribution of core and noncore genes for prevalent KEGG pathways. Colors indicate broader KEGG class of each pathway. Pathways with * are functionally enriched in noncore based on Fisher exact test with FDR <0.05. See also Figure S5.

Interestingly, combining functional capacity with strain signatures revealed that similar pangenomic capacity can be achieved with distinct strain combinations (Figure 6A, 6B). This suggests that functional niche saturation can occur through multiple combinations of a limited number of strains, rather than requiring a full phylogenetic complement. Thus, while individuals have distinct P. acnes strain signatures, their functional capacities are remarkably similar, only differing 5% between individuals, and an individual is more likely to retain those unique genes over time (Figure 6C). This percent difference between individuals increases to ~40 when relative copy numbers of genes are compared (Figure 6D), further illustrating that while individuals’ microbial communities have the same set of genes their relative abundances vary.

In addition to gene retention, we also observed an increase in pangenome size in three of the individuals (HV01, HV02, and HV03) over time (Figure 6A). Thus, although strain stability is more typical, individuals can acquire new gene content over time. Because of convergent functionality between individuals, we redefined the core genome to represent a ‘functional’ core that is characteristic of healthy communities, towards which a probiotic approach might strive. We defined that 2,982 genes were present in at least 50% of sites in all individuals at all times, exceeding the 1,860 genes that are derived using genomes alone. This functional core genome increases to 3,186 genes when we exclude HV06 whose strain community is predominated by a single P. acnes clade (Figure S4A).

To evaluate functional enrichment within the 769 genes not identified as core, we assigned KEGG pathway annotations to clusters with BLAST. Unsurprisingly, these noncore genes were statistically enriched for pathway “None” (702 of 769 genes, Figure 6E), underscoring that more extensive gene annotations are needed to better understand functional variation. After removing unannotated clusters, we found accessory genes to be enriched in functions associated with ABC transporters and cysteine and methionine metabolism (Figure 6F). Ubiquitous across bacteria, ABC transporters facilitate communication between bacteria and the environment through the active transport of substances such as ions, sugars, lipids, proteins, and drugs across membranes, contributing to nutrient sensing and other processes.

Identification of SNVs in the P. acnes core

For high coverage back and manubrium samples (Figure S4C), metagenomic reads were mapped against the P. acnes core genome and the alignment file was processed with GATK’s IndelRealigner (McKenna et al., 2010), samtools (Li, 2011), and bcftools to identify possible variants with criteria previously described (Lieberman et al., 2014).

Pangenome analyses of dominant skin species

As illustrated in Figure S5A, P. acnes and S. epidermidis nucleotide-coding sequences were clustered into non-redundant orthologs with usearch (Edgar, 2010) and validated with single copy marker genes (Greenblum et al., 2015). A gene was considered present when 40% of its length was covered with reads (Zhu et al., 2015).

Amynah Pradhan, Sharon Osgood, Clay Deming, Cynthia Ng, Brian Schmidt, Pamela Thomas and Sean Conlan provided underlying efforts; Segre lab and Mark Udey engaged in helpful discussions. Work was supported by NHGRI and NCI Intramural Research Programs, and Chanel/CE.R.I.E.S. research award to J.A.S. Study utilized the high-performance computational capabilities of the NIH Biowulf Linux cluster. Sequencing was funded by grants from National Institutes of Health (4UH3AR057504-02). IGS Analysis Engine at University of Maryland School of Medicine provided structural and functional annotation of genomes.