Individual and household attributes influence the dynamics of the personal skin microbiota and its association network
© The Author(s). 2018
Received: 11 April 2017
Accepted: 19 January 2018
Published: 2 February 2018
Numerous studies have thus far characterized the temporal dynamics of the skin microbiota of healthy individuals. However, there is no information regarding the dynamics of different microbial association network properties. Also, there is little understanding of how living conditions, specifically cohabitation and household occupancy, may be associated with the nature and extent (or degree) of cutaneous microbiota change within individuals over time. In this study, the dynamics of the skin microbiota, and its association networks, on the skin of urban residents over four seasons were characterized.
Similar to western cohorts, the individuals of this cohort show different extents of variations in relative abundance of common skin colonizers, concomitant with individual- and household-associated changes in differential abundances of bacterial taxa. Interestingly, the individualized nature of the skin microbiota extends to various aspects of microbial association networks, including co-occurring and excluding taxa, as well as overall network structural properties. Household occupancy is correlated with the extent of variations in relative abundance of Propionibacterium, Acinetobacter, and Bacillus over multiple skin sites. In addition, household occupancy is also associated with the extent of temporal changes in microbial diversity and composition within a resident’s skin.
This is the first study investigating the potential roles household occupancy has on the extent of change in one’s cutaneous microbiota and its association network structures. In particular, we show that relationships between the skin microbiota of a resident, his/her cohabitants, and those of non-cohabitants over time are highly personal and are possibly governed by living conditions and nature of interactions between cohabitants within households over 1 year. This study calls for increased awareness to personal and lifestyle factors that may govern relationships between the skin microbiota of one individual and those of cohabitants, and changes in the microbial association network structures within a person over time. The current study will act as a baseline for future assessments in comparing against temporal dynamics of microbiota from individuals with different skin conditions and for identifying residential factors that are beneficial in promoting the dynamics of the skin microbiota associated with health.
A community of microbial life forms, constituting of bacteria, fungi, viruses, and parasites make up the commensal skin microbiota responsible for modulating the host immunity responses and preventing colonization and invasion by pathogens . Various medical and allergic conditions are associated with alterations in one’s overall cutaneous microbial community and dysbiosis between specific skin colonizers [2–4], highlighting the importance of the microbiota in cutaneous health. Thanks to advances in sequencing technology, the scientific community has gained a wealth of information regarding the skin microbiota, including the biogeography of the skin microbiota, roles of host properties, activities, and the environment on cutaneous microbial communities, co-abundance and co-exclusion correlations between taxa, and changes of microbiota associated with health and disease [4–15].
Recent investigations have also reported on the temporal dynamics of the cutaneous microbiota [2, 10, 16–21], highlighting the highly personal nature of the skin microbiota across healthy individuals over weeks and months. Changes in the microbiota may also be associated with cutaneous or immunological conditions, or as responses to clinical interventions [7, 19]. Altogether, these studies reinforce the importance of understanding potential relationships between the environment, host characteristics, and microbiota variations over time, and whether microbiota changes over time is related to temporal variations in disease phenotypes and prognosis .
As with other ecosystems, microorganisms found on skin engage in co-abundance and exclusion relationships . While correlation does not equate causation, microbial association networks can be analyzed to infer potential ecological relationships between community taxa, as well as relationships between taxa and the environment. Recent works have performed correlation network analyses to characterize potential relationships between specific taxa within the skin microbiota [5, 6, 22]. However, understanding association networks at the structural level provides additional information on how microbial associations as a whole respond to different human or temporal variables, potentially revealing how this global association framework changes, and the factors that may drive such changes over time.
Previous studies have illustrated the effects of cohabitation on the skin microbiota within individuals. Specifically, our previous work  recapitulates the discoveries of other studies [19, 23] that the cutaneous bacterial assemblage of cohabitants are more similar compared to non-cohabitants. While these works complement earlier reports describing the effects of lifestyles and living conditions on occupant skin microbiota [13, 15], no study has extensively and directly evaluated the potential roles of cohabitation and household occupancy in shaping the nature and extent of changes in the skin microbiota and its association networks within individuals over time. Given that cohabiting individuals and their microbiota may facilitate the transmission of microorganisms , a greater insight into the potential roles of household properties such as occupancy on an individual’s skin microbiota and its dynamics may be of clinical significance.
Therefore, following our previous skin microbiota studies involving a single time point [5, 6], we now provide an account of the skin microbiota (while the term microbiota commonly encompasses all microbial domains, the term is hereafter referred to as the total bacterial and archaeal community) and its dynamics, of 24 individuals within Hong Kong (HK) households over the period of 1 year. Different aspects of the nature and extent of microbiota dynamics, with a focus on the roles of cohabitation and household occupancy, are described. We show that while our observations are consistent with the individuality of the dynamics of cutaneous microbiota shown in previous works [10, 20], we extend this personalized property to the dynamics of microbial association networks. At the same time, we suggest that the extent of changes in different characteristics of the skin microbiota may be associated with the number of cohabitants an individual resides with. By focusing on the changes of multiple aspects of the cutaneous microbiota at individual and household levels, this study sets the stage for comparative analyses pertaining to microbiota dynamics between healthy individuals and those with various skin conditions, thereby potentially identifying personal, household, and related factors associated with healthy microbiota dynamics.
Results and discussion
In this study, we first characterize and discuss the nature and degree of variability of various facets of the skin microbiota within healthy individuals. This is followed by a focus of how effects of cohabitation and residential occupancy may play roles in shaping the observed differences in the degree of temporal variability in one’s skin microbiota over time.
Taxonomic and sub-genus overview of the skin microbiota over time
Top genera (i.e., those with average relative abundance of ≥ 1% across the entire dataset) are detected across all individuals and skin sites and comprise under 60% of the skin microbial community, reflecting the overall taxonomic diversity of the cutaneous microbiota (Additional file 1: Figure S1 and Additional file 2: Table S1) . The top genera include Enhydrobacter, which has been detected in previous reports and is considered to be enriched on skin surfaces of Chinese individuals [5, 26–28], suggesting that this genus is a prevalent and stable skin colonizer in subjects of this cohort. Information is currently lacking on the ecological and physiological bases for the apparent enrichment of Enhydrobacter in the Chinese. However, an association between relative abundance of Enhydrobacter on skin and facial sebum level and hydration has been reported in a recent Asian study . Further examination of genetic and/or lifestyle factors that may affect skin physiology (such as living conditions, environmental exposure, food intake, and use of cosmetic and sanitary products [11, 25, 29–31]) may help elucidate the physiology of this colonizer on skins of Chinese and other host populations.
We sought to examine the extent of change in relative abundance of a given genus within an individual over the four seasons. Here, the extent of change in relative abundance for a particular genus is expressed as its coefficient of variation (CV). Regardless of skin site, each individual presents a distinctive signature of CV for each of the top genera (Additional file 2: Table S1 and Additional file 3: Figure S2). Previously, it was revealed that the relative abundance and the variations (as expressed by the standard deviation) of microbial taxa within the skin microbiota of a western cohort exhibit a second-order relationship, where low- and high-abundance microbial taxa vary less than those that are moderately abundant . In contrast, our HK cohort shows strong linear correlative relationships across individuals (Additional file 4: Figure S3). As no genus in our dataset represents an average of over 45% of any individual’s microbiota, our correlative relationship between relative abundance and standard deviation over time is expected to be different from the previous study, where the second-order relationship observed appears to be strongly driven by taxa detected in over 50% of a subject .
Microbial diversity variations within individuals over time
Shannon diversity, taking into account both community richness and evenness, was calculated to determine the nature and range of temporal variability in microbial diversity across individuals. When all samples were included for analysis, significant differences in Shannon diversity are observed between individuals, seasons, and skin sites (FDR-adjusted p ≤ 0.02 for all, Additional file 7: Table S3). Within each season, individuals and households are the major factors explaining significant differences in diversity. Linear mixed model also suggests that each of these major factors is potential drivers for differences in Shannon diversity (p < 0.001 for the three factors of season, individuals, and households. For each test, the other two factors, as well as sequencing batches, were considered as random variables). Significant differences in forehead and forearm diversities are observed between households, while forearms also exhibit significant differences in diversity between individuals (Additional file 7: Table S3).
To understand the personal rate of change in microbial diversity over the course of 1 year, we calculated the CV of Shannon diversity at a given site within an individual. Considering all samples, the Shannon diversity CVs are different between individuals, households, and occupancy (Kruskal-Wallis (KW) FDR-adjusted p < 0.005 for all). In contrast, CVs are not significantly different between samples of different skin sites, age group, and gender (FDR-adjusted p > 0.05, KW for skin site and age group, and Mann-Whitney (MW) for gender). Communities on left forearm sites varied in CVs between individuals and households (KW FDR-adjusted p < 0.04 for both, Additional file 8: Figure S5a-b). Shannon diversity CVs in forearms are also significantly higher in adults (average within-individual CV = 0.145) compared to the elderly (CV = 0.0871, KW post hoc pairwise p < 0.05, Additional file 8: Figure S5c).
Microbial compositional variations within individuals over time
We analyzed variations in microbial community compositions between subjects, households, and within individuals over time using weighted UniFrac distances, taking into account both the presence and abundance of OTUs. Skin microbiota significantly cluster by season, individual, household, and age groups overall, as well as within separate skin sites, as shown by both analysis of similarities (ANOSIM) and permutational multivariate analysis of variance (PERMANOVA) analyses (Additional file 9: Table S4). Microbiota do not appear to cluster strongly by skin sites within each season in our cohort (and only marginally in summer).
Spearman’s correlation based on pairwise weighted UniFrac distance within individual over time
As with previous investigations [10, 17], the microbiota of the same site within some individuals between seasons closer together are more similar, as suggested by the significantly negative time-decay relationships (Additional file 9: Table S4). As samples for this study were collected within a single year, it is currently not known whether the decay relationships observed will persist for longer time frames as in other ecosystems , or whether an annual or seasonal cyclical pattern will emerge over longer periods of time, as seen in the gut microbiota of different mammals  and potentially in humans . A greater understanding of time-decay relationships of cutaneous microbiota has forensic implications, where the identifiability of a person’s microbiota on a surface may decrease over time . Given that time-decay relationships are more likely to be a personal trait rather than one shared by members of a cohort or household level (based on our results here and those reported previously ), additional works will be required to understand why some subjects present time-decay relationships and some do not, and how inter-individual time-decay variability affects the feasibility of using microbiota data in forensic applications.
Overall, our observations are congruent with previous western studies looking at the temporal dynamics of the skin microbiota [10, 17, 20]: (1) individuals experience ranges in the extents of temporal changes in relative abundance of genera and overall microbial diversity and community composition, (2) time-decay relationships in microbiota are present but not a ubiquitous property on skin, and (3) within the same individual, multiple skin sites present significant correlations in the extent of microbiota variations at any two time points, and (4) there exists significant correlations between the extent of temporal variations between different skin sites. Therefore, despite known differences in microbial community compositions across continental populations [5, 11, 13–15, 30], some temporal properties of cutaneous microbiota may be more universal. Such knowledge regarding the behaviors of microbial communities within and between individuals over time may provide fundamental basis for subsequent and more focused microbiota analyses in special populations around the globe, such as those adopting different lifestyles and are living in drastically different environments .
Dynamics of microbial networks are individualized in associating taxa and overall network structure
While the dynamics of microbial diversity and community composition within individuals have been examined above and in previous works [7, 10, 17, 20], little information is currently available on the temporal properties of microbial correlative association networks. As correlative association may represent ecologically significant relationships between microbial constituents, tracking correlations over time in each individual will allow the assessment on the factors governing network structure property at a personal level. Here, insights into SParse InversE Covariance Estimation for Ecological Association Inference (SPIEC-EASI)  were used to (1) identify co-abundance and co-exclusion relationships between specific taxa within individuals and (2) determine whether within-individual network structure dynamics can potentially be explained by individual, household, or temporal factors.
While strong positive associations are detected within OTUs of Enhydrobacter across individuals of different households (QB3z and SW3z), this genus is also involved in positive associations with OTUs of other genera (OTU_4 of Enhydrobacter with OTUs of Chryseobacterium and Acinetobacter within households of STW and TMb, respectively Additional file 14: Figure S9). Consistent with our previous work , Enhydrobacter is negatively associated with OTUs of Streptococcus (HFC3z, Additional file 14: Figure S9). OTUs of Enhydrobacter and Corynebacterium can be involved in both positive and negative associations (within individual TMb3x), demonstrating that species and strains of Enhydrobacter and other genera deserve further analyses in not only its physiology on skin but also their potential interactions within the skin microbiota.
With the exception of individual MOS3x, most of the associations between OTUs of Propionibacterium and Staphylococcus are positive correlations. Propionibacterium acnes and Staphylococcus epidermidis (the dominant Staphylococcal species on skin) can be mutually inhibitory, but genetic variations between strains influence the extents to which one species affect the physiology and survival of the other . Previous species-level analysis of Staphylococcus in our cohort reveals that as detected in some healthy individuals in the western world , most OTUs of this genus are in fact Staphylococcus aureus . Intriguingly, S. aureus has been shown to respond to molecules secreted by P. acnes at particular conditions . While correlation analyses act as gateways for identifying interactive relationships between microbial members, it provides no information on the environmental and strain physiological states necessary for such interactions to occur. Coupled with the recent identification of separate evolutionary patterns between continental strains of S. aureus , the results here again beckon the need for shotgun metagenomics sequencing to examine microbial associations at higher resolution, not only to understand species- and strain-level patterns in microbial associative relationships but also the metabolic and physiological bases for such relationships .
Roles of cohabitation on microbiota changes within individuals over time
Roles of occupancy on degree of microbiota changes within individuals over time
Members within a household do not appear to share similar signatures of CV in relative abundance of specific genera, depending on the genus and skin site. However, significant correlations can be observed between CVs of relative abundance and household occupancy. Specifically, individuals within households of higher occupancy are associated with greater CVs of Propionibacterium and Acinetobacter over the period of the year (Additional file 17: Figure S11a-c). Conversely, a significantly negative correlation, as shown for Bacillus on palms, suggests that individuals with fewer cohabitants tend to have greater rates of changes in relative abundance of this genus (Additional file 17: Figure S11d). In addition to changes in relative abundance of specific genera, individuals residing in households with higher occupancy also present higher CVs of Shannon diversity on left forearm but not for the right forearm and other sites (Additional file 17: Figure S11e). The significance on the left forearm is intriguing and may perhaps be related to the combined effects of handedness (since all but one individual in this study are right-handed) and environmental exposure. However, this is speculative, and additional investigations in controlled settings will be required to assess the roles of these effects, and with a focus on how cohabitation, handedness, and environmental exposure affect the temporal dynamics of skin microbiota. Thus, it is currently not known whether the different correlation trends in different skin sites are ecologically significant. Metagenomic shotgun sequencing may prove beneficial in explaining any potential ecological and physiological rationale for the correlative observations here.
To examine whether the number of cohabitants within a household is also associated with variations in community composition of the skin microbiota of each occupant and site, Kendall τ correlation was performed between the mean UniFrac distances for a particular site within an individual between any two seasons, and the number of occupants present in his/her household. Forehead and left forearm sites from individuals living in residences of higher occupancy show more temporally variable microbial communities over the course of the year (Additional file 17: Figure S11f).
Overall, the number of residents within a household appears to be associated with changes in the dynamics of the personal skin microbiota over 1 year at multiple levels: (1) the extent of variation in the relative abundance of specific colonizers (CVs of Propionibacterium, Acinetobacter, and Bacillus relative abundance across multiple skin sites), (2) the variation in community diversity (as determined by CVs of Shannon diversity within each individual on forearm sites), and (3) the variation in community structure (through average weighted UniFrac distances between forehead and left forearm sites of the same individuals over seasons). Our households are limited to a maximum number of six occupants, and it is currently not known how residences with more individuals, or other built environment settings (such as workplaces and other public spaces of higher occupancy) would affect one’s skin microbiota and its dynamics.
However, it is unlikely that occupancy directly shapes microbiota changes within occupants, as individuals living in different households of the same occupancy still have different microbiota and variable CVs. Instead, we propose that the number of residents within a household indirectly affects an occupant’s microbiota dynamics by how an occupant interacts with the cohabitants and their microbiota. A household with more people may facilitate increased contact, which may facilitate direct transfer of personal skin taxa and potentially pathogenic microorganisms [12, 24]. Alternatively, occupants may release personal microbiota into the residential air and onto surfaces, and microbial members may subsequently be picked up by its cohabitants through passive exposure or active surface touching [21, 32]. We suggest that over single or multiple time points, the relationships between one’s own microbiota and that of its cohabitants and others (that is, how similar or different their microbiota are, Fig. 4 and Additional file 16: Table S6) depend on the biogeographical (across skin sites within a time point) and temporal (within a site over multiple time points) variability of a person’s microbiota (which is in itself now known to be highly personal [10, 20]) and the nature of interactions between the person and his/her cohabitants. Individual factors, in addition to household occupancy and cohabitation, may as a result influence the temporal stability of one’s baseline microbiota over time, which is consistent with our pairwise UniFrac (Additional file 16: Table S6) and SourceTracker (Fig. 5) observations. Therefore, occupancy and cohabitation represent part of a collection of household properties, along with living conditions and lifestyles  that may influence the nature and the extent of change in one’s skin microbiota over time. In addition, given that population density in HK is among the highest in the world, it is currently not known how household properties that take into account the crowdedness of households (such as occupant density or privacy indexes as defined previously for microbiota of residences ) will influence skin microbiota of residents. We believe that separating these household properties, perhaps via controlled chamber and inoculation experiments [7, 17], will be crucial in specifying the relative proportion of contribution of each household factor in shaping an individual’s microbiota dynamics.
While our findings corroborate previous western studies in identifying the personalized properties of skin microbiota dynamics, we have highlighted in this study possible associations between occupancy to multiple aspects of microbiota dynamics within individuals, including for the first time, the changes of microbial association network structures within every individual over 1 year. Our findings underscore the need for close examination on individual physiological and lifestyle properties in explaining changes in cutaneous microbiota but also an in-depth assessment of the person’s living environment, including number of cohabitants and their interactions over time. Also, while a daunting task, a more individualistic approach may therefore be beneficial to determine factors governing changes in microbial association network structures within an individual, thereby gaining understanding on properties associated with network dynamics associated with adverse skin conditions.
Finally, there is little information regarding whether individuals suffering from cutaneous ailments present different temporal patterns on their skin microbiota from that of healthy individuals, in ways similar to how gut microbiomes of individuals suffering from various gut conditions have distinctive microbiota temporal traits . Data presented here therefore can act as community-level baseline in understanding how the microbiota dynamics of healthy individuals compare with those with various skin conditions. Future longitudinal and comparative microbiota works will potentially provide perspectives on whether disease onset, progression, and recovery are related to changes in microbiota-host relationships over time and identify host and other factors that may be associated with these changes.
Cohort characteristics and sampling
A total of 480 skin samples from this study are part of a larger seasonal analysis of skin, air, and surface microbiomes of Hong Kong residences , and a continuation of previous single-season skin microbiota works [5, 6]. Ethics approval for subject sampling was granted by the City University of Hong Kong Ethics Committee (reference number 3-2-201,312 (H000334)). After being informed about the nature of the study, as well as their roles and responsibilities as subjects, written informed consent was given by all individuals.
From 24 healthy individuals (coded one of 3u to 3z) of 11 households (ADMa, FH, HFC, HHa, MOS, QB, STW, SW, TK, TMb, WKS), samples were collected in January (winter), April/May (spring), August (summer), and November/December (autumn) of 2014 (Additional file 18: Table S7). All samples were collected between 19:00 and 20:00, as individuals were most likely to be at home during those times and to minimize the disruption caused to the individuals due to sampling. For each individual, five skin samples (forehead, left/right forearm, left/right palm) were collected at his/her residence by samplers wearing gloves sterilized with 70% ethanol to minimize contamination or skin taxa transfer between sampler and the subject. All subjects included this study are ethnically Chinese and are long-term residents of Hong Kong. Subjects of this study had not taken antibiotics and antifungals at least 3 months prior to each sampling episode. The individuals in this study were living in households throughout rural and urban areas of Hong Kong to cover a broad local geographical scope. Individuals and households were selected to cover a range of age and lifestyle choices such as smoking, pet ownership, and allergies. All households involved in this study did not use pesticide or have purchased new furniture during the course of the sampling periods. While antibiotic ingestion may or may not have an effect on skin microbiota dynamics [19, 20], we note that all subjects had not taken antibiotics 1 month prior and after each sampling period. Subjects were instructed to not wash their hands or shower immediately prior to sampling. Samples were collected as previously described . Briefly, skin samples were collected using sterile cotton swabs moistened with 100 μL swab solution (0.15 M NaCl and 0.1% Tween 20) , and each sample was collected using a back-and-forth swabbing motion for 15 s. Samples were stored in − 80 C within 1 h of sampling and were stored until genomic DNA (gDNA) extraction.
Genomic DNA extraction, 16S rRNA gene amplification and sequencing
Following sampling, gDNA was extracted using the PowerSoil DNA Isolation Kit (MO BIO Laboratories, Inc., Carlsbad, CA, USA), with modifications as described previously . Sterile cotton swabs not exposed to skin surfaces were included in the extraction process as negative controls to account for possible contamination. Purified gDNA samples, including negative controls, were sent to Health Genetech Corporation (Taoyuan City, Taiwan) for PCR, library preparation, and sequencing. PCR, library preparation, and Illumina MiSeq sequencing were prepared as described [5, 51]. Briefly, the 515f/806r primer pair was used to target and amplify the V4 hypervariable region of the 16S rRNA gene . PCR amplification was performed in a 20-μL reaction volume containing 10 μl 2× Phusion HF master mix (New England BioLabs, Ipswich, MA, USA), 0.5 μM each forward and reverse primer, consisting of customized barcodes present on both primers for multiplex sequencing, and 50 to 150 ng DNA template. The PCR conditions consisted of an initial of 98 °C for 30 s, followed by 30 cycles of 98 °C for 10 s, 54 °C for 30 s, and 72 °C for 30 s, as well as a final extension of 72 °C for 5 min. Positive amplification was verified by agarose gel electrophoresis. Amplicons in triplicates were pooled and purified using AMPure XP beads (Agencourt, Brea, CA, USA) and quantified using a Qubit double-stranded DNA HS assay kit on a Qubit fluorometer (Invitrogen, Carlsbad, CA, USA), all according to the respective manufacturers’ instructions. For library preparation, Illumina adapters were attached to amplicons using the Illumina TruSeq DNA sample preparation kit v2. Purified libraries were applied for cluster generation and sequencing on the MiSeq platform using paired-end 150-bp reads.
A total of 12,212,096 forward sequences in .fastq format underwent quality filtering using the “fastq_filter” command in USEARCH , based on a maximum expected error of 1 error/read, reads trimmed to a length of 139 bp and shorter reads removed. Eighty percent of reads (9,792,397 reads) were of acceptable quality. Filtered reads were then de-multiplexed and binned into OTUs using UPARSE with 97% sequence identity threshold . OTUs were provided with taxonomic information using the “assign_taxonomy.py” command in QIIME  based on the Greengenes reference database (May 2013 version, 99,322 sequences). Chimeras were detected in UCHIME2 (high-confidence mode) against the Greengenes database to minimize risk of identifying false-positive chimera OTUs . OTU lineages present in an average of > 5% of reads in negative controls (each control averages 20,000 reads and 600 OTUs) were considered contaminants and were removed from all samples. OTUs classified as chloroplasts and mitochondria were also removed. Following read filtering and OTU removal, a total of 8,093,026 reads are retained for community analyses described below. Following quality control and OTU taxonomic classification, archaeal OTUs represent ~ 0.3% of reads detected.
The coefficient of variation (CV) of relative abundance for top genera for each individual over the four seasons is calculated by the standard deviation of the genera within that individual, divided by the mean of the relative abundance of that genera in the individual. A low CV represents a relatively stable relative abundance for that given genus on a given individual, and high CV represents less stable relative abundance. Total Shannon (alpha-) diversity for each sample was estimated using the “breakaway” package (version 4) in R available from CRAN . Shannon diversity values presented in Additional file 7: Table S3 are average values following 999 reiterations. Similar to previous skin microbiota works [10, 20], the CV for Shannon diversity within each individual, representing the change of within-sample microbial diversity of an individual over the four seasons, is calculated by the standard deviation of the Shannon diversity values for that individual, divided by the mean of the Shannon diversity of the same individual. Weighted UniFrac distances were computed using the QIIME script “beta_diversity.py.” Analyses of similarities (ANOSIM) were computed on based on β-diversity data using “vegan” package in R. The PERMANOVA pseudo-F statistic and significance were calculated using QIIME’s “compare_categories.py” script (with 999 permutations) based on the weighted UniFrac distance matrix generated from QIIME. Differential abundance of OTUs between sample groups were performed with DeSeq2 using the QIIME script “differential_abundance.py,” performed separately for each individual and each season pair. OTUs that showed significant DeSeq2 results in more than half of the individuals for a given season pair were included for negative binomial and zero-inflated negative binomial mixed model analyses (“glmmADMB” package in R) to confirm significant differences in abundance for these OTUs, using season as a fixed effect and count data as response variable. Variations between individuals, skin sites, and sequencing batches were considered as random effects. OTUs that appear in less than 20% of samples, and OTUs with less than 100 reads were not included for the analysis. OTUs of top genera that were detected at an average of ≥ 1% relative abundance and were present across all samples were subjected to oligotyping (version 2.1 from http://merenlab.org/software/oligotyping/) . To remove oligotypes with low read counts, a minimum substantive abundance (M) of 10 was adopted. A linear mixed effects model using the R package “lme4” was used to determine whether Shannon diversity differences can possibly be driven by seasonal, individual, and household factors, with individual, site, household, and/or sequencing batch differences as random effects when not considered as a predictor effect. The analysis of variance (ANOVA) on the hypothesized and null models was performed to determine statistical significance in Shannon diversity differences for each predictor factor.
Microbial source-tracking analysis
To examine the roles of cohabitants as potential microbial sources of an occupant’s skin in each residence and time point, Bayesian SourceTracker  analysis was used to estimate the proportional contribution of the cohabitants’ skin microbiome to the individual microbiome at the given season. OTUs present in less than 10% of the samples were excluded for the analysis. SourceTracker predictions were generated for each individual across four seasons, with the skin samples of each individual in one season as the sinks, and the skin samples of the remaining individuals in the same season as potential sources. The results were presented as the source proportions from cohabitants, non-cohabitants, and unknown sources based on whether the sink and source communities can be sourced from the same household.
Microbial association network analysis
SParse InversE Covariance Estimation for Ecological Association Inference (SPIEC-EASI) was used to assess potential ecological associations between microbial taxa. SPIEC-EASI has been recently applied to other microbiota investigations [58, 59] and is currently one of the preferred methods for association network analysis, as it is robust to issues of compositional bias, conditional independence, and dimensionality, all properties commonly encountered in microbiota data . Two sets of association network analyses were performed, where (1) microbiota data for all seasons within an individual was combined (to assess for overall microbial associations within an individual, set of 24 networks), and (2) where microbiota data for each season within each individual was analyzed separately (to assess for dynamics of association networks within individuals, set of 96 networks). For each group, abundance data over different samples underwent centered-log ratio transformation. OTUs with less than 150 and 10 reads were not included in the first and second analysis, respectively. Based on transformed data, the Stability Approach to Regularization Selection (StARS) method , suitable for high-dimensional data, was used to infer the network structure, employing the node-based neighborhood selection procedure, with a minimum lambda ratio of 0.01 and reiteration of 50 times as recommended . Networks as shown in Fig. 2 were constructed with Cytoscape (version 3.4.0) .
Using SPIEC-EASI, general network structural attributes, including degree distribution (distribution of the number of significant associations (i.e., edges) per each OTU (i.e., node), and natural connectivity (a proxy for network stability and robustness to targeted node removal, by assessing the extent of alternative paths present between any two nodes) for each network were also determined. Natural connectivity for network within each individual and season was assessed based on node removal of either decreasing node betweenness centrality, which is a measure of the paths between a node to all other nodes, or decrease node degree, which is a measure of the number of edges a particular node contains. Natural connectivity is presented as the proportion of the total network size, and a higher natural connectivity represents a more stable and robust association network upon node removal. Graphlet correlation distance approach as described by Yaveroğlu and colleagues  was also performed on the set of 96 networks to compare the roles of season, household, and individuals in shaping association network structure. Briefly, each network was broken down into subgraphs (graphlets), each consisting of up to four nodes, and performed the graphlet frequency distribution of each node across networks. A graphlet correlation matrix was generated for each of the 96 networks based on pairwise Spearman correlations of the 11 non-redundant orbits of all nodes. This matrix is then used to calculate the graphlet correlation distances, and relationships between networks which were then visualized into multidimensional scaling plots constructed using R. Each network is represented by a position on the plot, and networks that are closer to each other in the MDS plot space are considered to be more similar in network structure.
The nonparametric Mann-Whitney (MW) and Kruskal-Wallis (KW) tests were employed to determine significance when comparing between two or more comparison groups, respectively. p values were adjusted for multiple comparisons using the false-discovery rate (FDR) algorithm in R. Potential batch effect between sequencing runs were considered as a possible confounding factor when performing multiple comparisons throughout the study (Additional file 19: Table S8) and have been treated as an extra sample variable along with other sample variables for FDR-multiple comparison using “p.adjust” in R. Where indicated in the main text, post-hoc KW pairwise comparison tests for significance between individual groups were performed using the kruskalmc function in R package pgirmess (http:// cran.r-project.org/web/packages/pgirmess/index.html), following significant KW observations (adjusted p ≤ 0.05). Spearman’s and Kendall’s τ ranked correlations were computed using the “cor.test” function in R.
We thank William Chan for sample processing, staffs of the Health GeneTech Corporation for their expertise in sequencing, and members of the sampling team (Ellen Li, Fred Kong, Wai Shan Chow, Catherine Chung, Ka Yan Ng, Yuet Ying Wong, and Flora Yeh) for their assistance. We are grateful for the participation of the cohort in this study.
This research was supported by the Research Grants Council of Hong Kong through Project 11211815.
Availability of data and materials
Raw reads in fastq format and metadata file for all samples (including negative controls) included in this study have been deposited in the NCBI SRA archive as BioProject PRJNA390040. In-house scripts for data analysis have been uploaded onto FigShare (https://figshare.com/articles/Skin_16S_temporal_microbiome/3410134).
MHYL collected samples, analyzed the data, developed and wrote the manuscript. XT analyzed the data and wrote the manuscript. DW wrote the in-house scripts and provided bioinformatics and technical support. HHLC processed samples for sequencing. PKHL guided and assisted in study design and analysis, as well as provided support for writing the manuscript. All authors read and approved the final manuscript in its current form.
Ethics approval and consent to participate
Ethics approval for subject sampling was granted by the City University of Hong Kong Ethics Committee (reference number 3-2-201,312 (H000334)). After being informed about the nature of the study, as well as their roles and responsibilities as subjects, written informed consent was given by all individuals.
Consent for publication
Consent was given by subjects to release personal data for publication as needed.
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
- Belkaid Y, Tamoutounour S. The influence of skin microorganisms on cutaneous immunity. Nat Rev Immunol. 2016;16:353–66.View ArticlePubMedGoogle Scholar
- Oh J, Freeman AF, Park M, Sokolic R, Candotti F, Holland SM, et al. The altered landscape of the human skin microbiome in patients with primary immunodeficiencies. Genome Res. 2013;23:2103–14.View ArticlePubMedPubMed CentralGoogle Scholar
- Hanski I, von Hertzen L, Fyhrquist N, Koskinen K, Torppa K, Laatikainen T, et al. Environmental biodiversity, human microbiota, and allergy are interrelated. Proc Natl Acad Sci U S A. 2012;109:8334–9.View ArticlePubMedPubMed CentralGoogle Scholar
- Xu Z, Wang Z, Yuan C, Liu X, Yang F, Wang T, et al. Dandruff is associated with the conjoined interactions between host and microorganisms. Sci Rep. 2016;6:24877.View ArticlePubMedPubMed CentralGoogle Scholar
- Leung MHY, Wilkins D, Lee PKH. Insights into the pan-microbiome: skin microbial communities of Chinese individuals differ from other racial groups. Sci Rep. 2015;5:11845.View ArticlePubMedPubMed CentralGoogle Scholar
- Leung MHY, Chan KCK, Lee PKH. Skin fungal community and its correlation with bacterial community of urban Chinese individuals. Microbiome. 2016;4:46.View ArticlePubMedPubMed CentralGoogle Scholar
- van Rensburg JJ, Lin H, Gao X, Toh E, Fortney KR, Ellinger S, et al. The human skin microbiome associates with the outcome of and is influenced by bacterial infection. MBio. 2015;6:e01315.PubMedPubMed CentralGoogle Scholar
- Findley K, Oh J, Yang J, Conlan S, Deming C, Meyer JA, et al. Topographic diversity of fungal and bacterial communities in human skin. Nature. 2013;498:367.View ArticlePubMedPubMed CentralGoogle Scholar
- Oh J, Byrd AL, Deming C, Conlan S, NISC Comparative Sequencing Program, Kong HH, et al. Biogeography and individuality shape function in the human skin metagenome. Nature. 2014;514:59–64.View ArticlePubMedPubMed CentralGoogle Scholar
- Oh J, Byrd AL, Park M, Kong HH, Segre JA. Temporal stability of the human skin microbiome. Cell. 2016;165:854–66.View ArticlePubMedPubMed CentralGoogle Scholar
- Perez GIP, Gao Z, Jourdain R, Ramirez J, Gany F, Clavaud C, et al. Body site is a more determinant factor than human population diversity in the healthy skin microbiome. PLoS One. 2016;11:e0151990.View ArticleGoogle Scholar
- Meadow JF, Bateman AC, Herkert KM, O’Connor TK, Green JL. Significant changes in the skin microbiome mediated by the sport of roller derby. PeerJ. 2013;1:e53.View ArticlePubMedPubMed CentralGoogle Scholar
- Clemente JC, Pehrsson EC, Blaser MJ, Sandhu K, Gao Z, Wang B, et al. The microbiome of uncontacted Amerindians. Sci Adv. 2015;1:e1500183.View ArticlePubMedPubMed CentralGoogle Scholar
- Blaser MJ, Dominguez-Bello MG, Contreras M, Magris M, Hidalgo G, Estrada I, et al. Distinct cutaneous bacterial assemblages in a sampling of South American Amerindians and US residents. ISME J. 2013;7:85–95.View ArticlePubMedGoogle Scholar
- Hospodsky D, Pickering AJ, Julian TR, Miller D, Gorthala S, Boehm AB, et al. Hand bacterial communities vary across two different human populations. Microbiology. 2014;160(Part 6):1144–52.View ArticlePubMedGoogle Scholar
- Grice EA, Kong HH, Conlan S, Deming CB, Davis J, Young AC, et al. Topographical and temporal diversity of the human skin microbiome. Science. 2009;324:1190–2.View ArticlePubMedPubMed CentralGoogle Scholar
- Costello EK, Lauber CL, Hamady M, Fierer N, Gordon JI, Knight R. Bacterial community variation in human body habitats across space and time. Science. 2009;326:1694–7.View ArticlePubMedPubMed CentralGoogle Scholar
- Kong HH, Oh J, Deming C, Conlan S, Grice EA, Beatson MA, et al. Temporal shifts in the skin microbiome associated with disease flares and treatment in children with atopic dermatitis. Genome Res. 2012;22:850–9.View ArticlePubMedPubMed CentralGoogle Scholar
- Abeles SR, Jones MB, Santiago-Rodriguez TM, Ly M, Klitgord N, Yooseph S, et al. Microbial diversity in individuals and their household contacts following typical antibiotic courses. Microbiome. 2016;4:39.View ArticlePubMedPubMed CentralGoogle Scholar
- Flores GE, Caporaso JG, Henley JB, Rideout JR, Domogala D, Chase J, et al. Temporal variability is a personalized feature of the human microbiome. Genome Biol. 2014;15:531.View ArticlePubMedPubMed CentralGoogle Scholar
- Wilkins D, Leung MHY, Lee PKH. Microbiota fingerprints lose individually identifying features over time. Microbiome. 2017;5:1.View ArticlePubMedPubMed CentralGoogle Scholar
- Faust K, Sathirapongsasuti JF, Izard J, Segata N, Gevers D, Raes J, et al. Microbial co-occurrence relationships in the human microbiome. PLoS Comput Biol. 2012;8:e1002606.View ArticlePubMedPubMed CentralGoogle Scholar
- Song SJ, Lauber C, Costello EK, Lozupone CA, Humphrey G, Berg-Lyons D, et al. Cohabiting family members share microbiota with one another and with their dogs. elife. 2013;2:e00458.PubMedPubMed CentralGoogle Scholar
- Ly M, Jones MB, Abeles SR, Santiago-Rodriguez TM, Gao J, Chan IC, et al. Transmission of viruses via our microbiomes. Microbiome. 2016;4:64.View ArticlePubMedPubMed CentralGoogle Scholar
- Fierer N, Hamady M, Lauber CL, Knight R. The influence of sex, handedness, and washing on the diversity of hand surface bacteria. Proc Natl Acad Sci U S A. 2008;105:17994–9.View ArticlePubMedPubMed CentralGoogle Scholar
- Ling Z, Liu X, Luo Y, Yuan L, Nelson KE, Wang Y, et al. Pyrosequencing analysis of the human microbiota of healthy Chinese undergraduates. BMC Genomics. 2013;14:390.View ArticlePubMedPubMed CentralGoogle Scholar
- Zeeuwen PL, Boekhorst J, van den Bogaard EH, de Koning HD, van de Kerkhof PM, Saulnier DM, et al. Microbiome dynamics of human epidermis following skin barrier disruption. Genome Biol. 2012;13:R101.View ArticlePubMedPubMed CentralGoogle Scholar
- Mukherjee S, Mitra R, Maitra A, Gupta S, Kumaran S, Chakrabortty A, et al. Sebum and hydration levels in specific regions of human face significantly predict the nature and diversity of facial skin microbiome. Sci Rep. 2016;6:36062.View ArticlePubMedPubMed CentralGoogle Scholar
- Bouslimani A, Porto C, Rath CM, Wang M, Guo Y, Gonzalez A, et al. Molecular cartography of the human skin surface in 3D. Proc Natl Acad Sci U S A. 2015;112:E2120–9.View ArticlePubMedPubMed CentralGoogle Scholar
- Ying S, Zeng DN, Chi L, Tan Y, Galzote C, Cardona C, et al. The influence of age and gender on skin-associated microbial communities in urban and rural human populations. PLoS One. 2015;10:e0141842.View ArticlePubMedPubMed CentralGoogle Scholar
- Staudinger T, Pipal A, Redl B. Molecular analysis of the prevalent microbiota of human male and female forehead skin compared to forearm skin and the influence of make-up. J Appl Microbiol. 2011;110:1381–9.View ArticlePubMedGoogle Scholar
- Wilkins D, Leung MHY, Lee PKH. Indoor air bacterial communities in Hong Kong households assemble independently of occupant skin microbiomes. Environ Microbiol. 2016;18:1754.View ArticlePubMedGoogle Scholar
- Lax S, Smith DP, Hampton-Marcell J, Owens SM, Handley KM, Scott NM, et al. Longitudinal analysis of microbial interaction between humans and the indoor environment. Science. 2014;345:1048–52.View ArticlePubMedPubMed CentralGoogle Scholar
- Fitz-Gibbon S, Tomida S, Chiu B-H, Nguyen L, Du C, Liu M, et al. Propionibacterium acnes strain populations in the human skin microbiome associated with acne. J Invest Dermatol. 2013;133:2152–60.View ArticlePubMedPubMed CentralGoogle Scholar
- N’Diaye AR, Leclerc C, Kentache T, Hardouin J, Poc CD, Konto-Ghiorghi Y, et al. Skin-bacteria communication: involvement of the neurohormone calcitonin gene related peptide (CGRP) in the regulation of Staphylococcus epidermidis virulence. Sci Rep. 2016;6:35379.View ArticlePubMedPubMed CentralGoogle Scholar
- Shade A, Gregory Caporaso J, Handelsman J, Knight R, Fierer N. A meta-analysis of changes in bacterial and archaeal communities with time. ISME J. 2013;7:1493.View ArticlePubMedPubMed CentralGoogle Scholar
- Maurice CF, Knowles S CL, Ladau J, Pollard KS, Fenton A, Pedersen AB, et al. Marked seasonal variation in the wild mouse gut microbiota. ISME J. 2015;9:2423.View ArticlePubMedPubMed CentralGoogle Scholar
- Davenport ER, Mizrahi-Man O, Michelini K, Barreiro LB, Ober C, Gilad Y. Seasonal variation in human gut microbiome composition. PLoS One. 2014;9:e90731.View ArticlePubMedPubMed CentralGoogle Scholar
- Kurtz ZD, Müller CL, Miraldi ER, Littman DR, Blaser MJ, Bonneau RA. Sparse and compositionally robust inference of microbial ecological networks. PLoS Comput Biol. 2015;11:e1004226.View ArticlePubMedPubMed CentralGoogle Scholar
- Faust K, Lima-Mendez G, Lerat J-S, Sathirapongsasuti JF, Knight R, Huttenhower C, et al. Cross-biome comparison of microbial association networks. Syst Microbiol. 2015;6:1200.Google Scholar
- Shade A, Jones SE, Caporaso JG, Handelsman J, Knight R, Fierer N, et al. Conditionally rare taxa disproportionately contribute to temporal changes in microbial diversity. mBio. 2014;5:e01371–14.View ArticlePubMedPubMed CentralGoogle Scholar
- Layeghifard M, Hwang DM, Guttman DS. Disentangling interactions in the microbiome: a network perspective. Trends Microbiol. 2017;25:217–28.View ArticlePubMedGoogle Scholar
- Christensen GJM, Scholz CFP, Enghild J, Rohde H, Kilian M, Thürmer A, et al. Antagonism between Staphylococcus epidermidis and Propionibacterium acnes and its genomic basis. BMC Genomics. 2016;17:152.View ArticlePubMedPubMed CentralGoogle Scholar
- Wollenberg MS, Claesen J, Escapa IF, Aldridge KL, Fischbach MA, Lemon KP. Propionibacterium-produced coproporphyrin III induces Staphylococcus aureus aggregation and biofilm formation. MBio. 2014;5:e01286–14.View ArticlePubMedPubMed CentralGoogle Scholar
- Ward MJ, Goncheva M, Richardson E, McAdam PR, Raftis E, Kearns A, et al. Identification of source and sink populations for the emergence and global spread of the East-Asia clone of community-associated MRSA. Genome Biol. 2016;17:160.View ArticlePubMedPubMed CentralGoogle Scholar
- Ma B, Wang H, Dsouza M, Lou J, He Y, Dai Z, et al. Geographic patterns of co-occurrence network topological features for soil microbiota at continental scale in eastern China. ISME J. 2016;10:1891–901.View ArticlePubMedPubMed CentralGoogle Scholar
- Knights D, Kuczynski J, Charlson ES, Zaneveld J, Mozer MC, Collman RG, et al. Bayesian community-wide culture-independent microbial source tracking. Nat Methods. 2011;8:761–3.View ArticlePubMedPubMed CentralGoogle Scholar
- Leung MHY, Lee PKH. The roles of the outdoors and occupants in contributing to a potential pan-microbiome of the built environment: a review. Microbiome. 2016;4:21.View ArticlePubMedPubMed CentralGoogle Scholar
- Ruiz-Calderon JF, Cavallin H, Song SJ, Novoselac A, Pericchi LR, Hernandez JN, et al. Walls talk: microbial biogeography of homes spanning urbanization. Sci Adv. 2016;2:e1501061.View ArticlePubMedPubMed CentralGoogle Scholar
- Halfvarson J, Brislawn CJ, Lamendella R, Vázquez-Baeza Y, Walters WA, Bramer LM, et al. Dynamics of the human gut microbiome in inflammatory bowel disease. Nat Microbiol. 2017;2:17004.View ArticlePubMedPubMed CentralGoogle Scholar
- Leung MHY, Wilkins D, Li EKT, Kong FKF, Lee PKH. Indoor-air microbiome in an urban subway network: diversity and dynamics. Appl Environ Microbiol. 2014;80:6760–70.View ArticlePubMedPubMed CentralGoogle Scholar
- Edgar RC. Search and clustering orders of magnitude faster than BLAST. Bioinformatics. 2010;26:2460–1.View ArticlePubMedGoogle Scholar
- Edgar RC. UPARSE: highly accurate OTU sequences from microbial amplicon reads. Nat Methods. 2013;10:996–8.View ArticlePubMedGoogle Scholar
- Caporaso JG, Kuczynski J, Stombaugh J, Bittinger K, Bushman FD, Costello EK, et al. QIIME allows analysis of high-throughput community sequencing data. Nat Methods. 2010;7:335–6.View ArticlePubMedPubMed CentralGoogle Scholar
- Edgar R. UCHIME2: improved chimera prediction for amplicon sequencing. bioRxiv. 2016;2016:74252.Google Scholar
- Willis A, Bunge J. Estimating diversity via frequency ratios. Biometrics. 2014;71:1042–9.View ArticleGoogle Scholar
- Eren AM, Maignien L, Sul WJ, Murphy LG, Grim SL, Morrison HG, et al. Oligotyping: differentiating between closely related microbial taxa using 16S rRNA gene data. Methods Ecol Evol. 2013;4:1111–9.View ArticlePubMed CentralGoogle Scholar
- Mahana D, Trent CM, Kurtz ZD, Bokulich NA, Battaglia T, Chung J, et al. Antibiotic perturbation of the murine gut microbiome enhances the adiposity, insulin resistance, and liver disease associated with high-fat diet. Genome Med. 2016;8:48.View ArticlePubMedPubMed CentralGoogle Scholar
- Morris A, Paulson JN, Talukder H, Tipton L, Kling H, Cui L, et al. Longitudinal analysis of the lung microbiota of cynomolgous macaques during long-term SHIV infection. Microbiome. 2016;4:38.View ArticlePubMedPubMed CentralGoogle Scholar
- Liu H, Roeder K, Wasserman L. Stability approach to regularization selection (StARS) for high dimensional graphical models. Adv Neural Inf Process Syst. 2010;24:1432–40.Google Scholar
- Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13:2498–504.View ArticlePubMedPubMed CentralGoogle Scholar
- Yaveroğlu ÖN, Malod-Dognin N, Davis D, Levnajic Z, Janjic V, Karapandza R, et al. Revealing the hidden language of complex networks. Sci Rep. 2014;4:4547.View ArticlePubMedPubMed CentralGoogle Scholar