Changes in the Human Skin Microbiome Over One Year’s Time

Human skin comprises a large number of distinguishable ecological niches. To describe fully the human skin microbiome, it will be necessary to identify the bacteria in each niche and to distinguish the commensal bacteria from the temporary residents. To contribute to the description of the human skin microbiome and employing a gene-based technology, we have identified the bacteria in two niches: the front and back of the base of the neck and over the course of one year. There were 50 volunteers and a total of 232 neck skin swabs. Roche 454 Tag pyrosequencing was employed to sequence a short hypervariable sequence region (V6) of the 16S ribosomal RNA gene. To identify the bacteria corresponding to the front and back of the neck for each volunteer, the “Classifier” software in the “Pyrosequencing” section of the Ribosomal Database Project was employed. The bacteria on virtually all 232 neck skin swabs were classified into bacterial Class. The skin microbiome of these two niches was composed principally of a mixture of five Classes of bacteria: Actinobacteria, Alphaproteobacteria, Bacilli, Betaproteobacteria and Gammaproteobacteria. The fraction of each Class could change over time. We could not distinguish the skin microbiome from the front of the base of the neck from the back of the base of the neck. At these two positions, we could not distinguish the male from the female skin microbiome. The principal variable was the time point. We concluded that the skin microbiome at the front and back of the base of the neck was composed principally of a mixture of five Classes of bacteria. The proportion of each Class could change over time.


INTRODUCTION
Human skin is a large organ with a multitude of distinguishable niches in contact with the external environment. In order to delineate the entire human skin microbiome, it will be necessary to identify the microbes in every niche, as a function of time and with people living in different geographical areas ). An excellent start has been made (Grice and Segre, 2011;Kong, 2011;HMPC, 2012).
Palo Alto, California, U.S.A., has a mild climate. However, that climate does change with the seasons. The season is part of the interactive environment of the skin bacteria and might have a significant impact on the quantitative composition of the skin microbiome. To investigate this possibility, we have identified the bacteria on two human skin niches (the front and back of the base of the neck) for fifty volunteers and over the course of one year.

Science Publications
AJM medical condition was ineligible to volunteer for this study. In total, there were 50 volunteers ( Table 1). Forty volunteers were enrolled at the first time point (February, 2010; Table 1A). An additional 10 volunteers were enrolled at the second time point (July, 2010; Table 1B). Although this was their first swab, we gave their swabs the designation "2" (for the time point) to keep the focus on the time point. No new volunteers were accepted at the third time point (February, 2011). The volunteers were not asked to change any part of their routine: e.g., neither their standard hygiene procedures nor the clothing that they chose to wear. This point may be important because, for example, Staudinger et al. (2011) demonstrated that the use of facial powder significantly increased the diversity of the human forehead skin microbiome. We have previously published the results from the first time point, albeit in a very different form (Hyman et al., 2012). We included those results here for the purpose of comparisons.

Amplifying and Purifying the V6 DNAs
Each volunteer took her/his own swabs from the base of the front of the neck (overlaying the suprasternal notch) and the base of the back of the neck (overlaying the posterior cervical vertebrae: the nape, below the hair line) (Hyman et al., 2012). Total DNA was isolated from each front and back neck swab employing a Qiagen DNeasy Blood and Tissue Kit. The total DNA was dialyzed and concentrated by the use of Amicon Ultra Centrifugal Filters (Millipore Corp). For PCR amplification of the V6 region of the 16S ribosomal RNA gene (rDNA), the forward and reverse primers were from Dethlefsen et al. (2008). Five identical reactions were run in parallel for each template. The PCR conditions were from Hyman et al. (2005). Following amplification, the five identical reactions for each template were pooled. The V6 DNA was purified by gel electrophoresis. We have employed Qiagen DNeasy Blood and Tissue Kits, Amicon Ultra Centrifugal Filters, PCR and purification of the amplicons by gel electrophoresis in previously reported experiments without problems of environmental contamination (Hyman et al., 2012;2005).

Pyrosequencing (Roche 454 Life Sciences)
For the first time point only (February, 2010), the front and back amplicons for each volunteer shared the same bar code (Hyman et al., 2012). Therefore, the combined front V6s were pyrosequenced separately from the back V6s. For the second time point, the front and back V6s had different barcodes , so they were combined before pyrosequencing, as was also the case for the V6s for the third time point.
The sequence reads were sorted by barcode and, thereby, assigned to the front or back of a specific volunteer. Then, the pyrosequencing reads were stripped of extraneous sequences. A data set was created that consisted of each unique sequence obtained for that sample and the number of times that sequence was represented in the sample. To identify the bacteria corresponding to the front and back of the neck for each volunteer, the "Classifier" software (Wang et al., 2007) in the "Pyrosequencing" section of the Ribosomal Database Project (RDP) was employed (Cole et al., 2009). Only reads that could be classified were considered further. There were reads too short to be classified. For example, for 3F, these short reads averaged 1.2% (n = 35) of the total. The DECIPHER software was employed to identify chimeras (Wright et al., 2012). In virtually all cases, the software identified the bacteria by Class. A very few reads in some sets were identified not by Class but by genus: e.g., TM7. In these cases, the reads were subsumed into Class. In a minority of cases, the software also identified the bacteria by Order. Actinobacteria were divided into Subclasses rather than Order. These were subsumed into Order.

RESULTS AND DISCUSSION
From the 50 volunteers, we achieved data for a total of 232 neck swabs, divided between the front (115) and the back (117) of the base of the neck ( Table 1). We have front and back swabs from three time points for 29 volunteers; two time points for 11 volunteers; and only one time point for 10 volunteers. The three time points were designated as 1, 2 and 3, respectively. Bacteria were identified by 454 Tag pyrosequencing of a short hypervariable sequence region (V6) of the 16S ribosomal RNA gene .

Class
In sampling an ecological niche, cost and time considerations prevent exhaustive sequencing. Therefore, statistical tests have been derived that indicate the closeness to saturation. Chao1 is one such test (Colwell and Coddington, 1994). The results of Chao1 analyses of the data are shown in Fig. 1 and demonstrate that the data are not close to saturation (see Discussion).
The bacteria on all 232 neck skin swabs were classified into Class. The complete data for all 232 swabs are presented on the Stanford Genome Technology Center (SGTC) web site (Table  S1; http://med.stanford.edu/sgtc/research/skin_microbio me.html). Table 2 presents examples of the bacterial identifications.
To make the data visually comprehensible, only those bacteria supported by, at least, 1% of the reads are reported in the tables. The overall average per swab (n = 232) of the total percent of the bacteria supported by more than 1% of the reads was 98.2+/-1.0% ( Table 2). We concluded that this very high percentage demonstrated that little information was lost when the focus was on those bacteria supported by, at least, 1% of the reads.
Fifteen Classes of bacteria were identified ( Table 3). The presence of five bacterial Classes was supported by more than 1% of the reads on virtually all of the swabs, albeit in very different amounts: Actinobacteria, Alphaproteobacteria, Bacilli, Betaproteobacteria and Gammaproteobacteria (Table 3). Gammaproteobacteria were supported by the most reads for 85% of the swabs, followed by Actinobacteria (10%) and Bacilli (5%). On the second and third time points, Flavobacteria and Sphingobacteria appear on virtually all swabs, also in very different amounts. On the other hand, the presence of some bacteria was supported by more than 1% of the reads only peripatetically. Bacteroidia were found on only 11 swabs (4.7%) and only in small amounts, as were Deinococci (one swab, 0.4%), Deltaproteobacteria (3 swabs; 1.3%), Fusobacteria (5 swabs; 2.2%) and Verrucomicrobia (one swab; 0.4%). Cyanobacteria were not found on swab 3, while Chloroplasts were found only on swab 3. The presence of Proteobacteria were supported by reads in the amounts of Gammaproteobacteria >>> Betaproteobacteria > Alphaproteobacteria. The blank spaces in the data tables ( Table  2  and  Table  S1; http://med.stanford.edu/sgtc/research/skin_microbiom e.html) signify that that particular bacterium's presence was supported by <1% of the sequence reads, if present at all. As one example, there were three blank spaces for Bacilli in the complete data table (Table  S1; http://med.stanford.edu/sgtc/research/skin_microbiom e.html). However, there were sequence reads supporting the presence of Bacilli on all three of those swabs: 33-1F, 0.43% of the reads; 33-1B, 0.70%; and 34-1F, 0.75%. Analogously, there were reads supporting the presence of

AJM
Alphaproteobacteria and Betaproteobacteria on all swabs, albeit sometimes at a frequency of < 1%. As additional examples for what was not reported in Table 2, for 01-3F, the presence of Deinococci was supported by three reads (0.41%) and the presence of Deltaproteobacteria was supported by four reads (0.54%).
Analyses were undertaken for four variables: location (front versus back), gender (male versus female), hair length in back (long hair versus short hair) and time point. For an analysis of the variable of location, we calculated the average percent per swab (front and back separately) for each of the five major bacteria at each of the three time points. The results of these calculations are presented in Fig. 2. For all five bacterial Classes, there was no statistically significant difference between the paired front and back swabs. As examples, for Gammaproteobacteria, 1F was not statistically significantly different from 1B (p = 0.88); 2F was not statistically significantly different from 2B (p = 0.62); and 3F was not statistically significantly different from 3B (p = 0.38). Table 2. Examples of Class data. The bacterial Class composition as percent of four skin swabs as a function of time and location.
Only those bacterial Classes supported by >1% of the sequence reads appear in the
For an analysis of the variable of gender, we compared the skin microbiome at the base of the neck of the 27 males to that of the 12 females. None of the gender comparisons yielded statistically significant differences (data not shown). As one example, we compared the percents of Actinobacteria and Gammaproteobacteria for men and women on swab 1F. There were no statistically significant differences (p = 0.15 and 0.12, respectively). In addition, PCA did not distinguish females from males ( Fig. S1; http://med.stanford.edu/sgtc/research/skin_microbiom e.html).
Thirteen women had long hair in back and six women had short hair in back. To determine the effect of long hair in back on the skin microbiome at the base of the back of the neck, we compared the percents of Actinobacteria and Gammaproteobacteria for women with and without long hair in back for swab 2B. The differences were not statistically significant (p = 0.70 and 0.85, respectively). In addition, PCA did not distinguish women with and without long hair in back ( Fig.  S1; http://med.stanford.edu/sgtc/research/skin_microbiom e.html).
Lastly, we undertook an analysis of the time points, which did produce statistically significant results. For example, 1F was statistically significantly different from 2F (p<0.001) and 2F was statistically significantly different from 3F (p<0.00001) and the same for the back swabs. PCA was undertaken for the variable of time point. The results are shown in Fig. 3. The time points clustered and 34% of the variation in P1 versus P2 and P1 versus P3 could be ascribed to time point.
Bacterial diversity was examined by calculating the widely used Shannon Diversity Index (SDI) for each swab (Shannon, 1948). For example, when only one bacterium is present, the SDI = 0; there is no diversity. The SDI data are presented in Table 4, which also contains the average values and standard deviations for each group of swabs. There was no statistically significant difference between the average SDIs of 1F (1.37+/-0.36) and 1B (1.38+/-0.29; p = 0.89), 2F (1.95+/-0.37) and 2B (1.94+/-0.33; p = 0.89) and 3F (2.50+/-.22) and 3B (2.46+/-0.16; p = 0.39). In addition, we were unable to distinguish the SDIs of males and females. As examples, the SDI of the male (n = 27) front swab 1F was not statistically significant different from the female (n = 12) front swab 1F (p = 0.23). The SDI of the male (n = 24) back swab 2B was not statistically significant different from the female (n = 18) back swab 2B (p = 0.46). When we compared the SDIs of the swabs as a function of time point, an interesting pattern emerged: the average SDIs of the time points were statistically different from each other. As examples, the average SDI of the front swabs from the first time point was statistically significantly different from the average SDI of the front swabs of the second time point (p<0.00001). The average SDI of the front swabs from the second time point was statistically significantly different from the average SDI of the front swabs from the third time point (p<0.00001). The analogous statistically significant differences hold true for the SDIs of the back swabs.
To investigate the change in diversity as a function of time in more detail, we determined that we had data for 129 transitions: e.g., from 01-1F to 01-2F and from 01-2F to 01-3F (Table  S2; http://med.stanford.edu/sgtc/research/skin_microbiome. html). We considered a change in the percentage of reads supporting the presence of a bacterium of less than an absolute 10% as "no change", an arbitrary, but reasonable amount which we have employed previously (Hyman et al., 2012). For example, for the 01-1 F to 01-2F transition, the percentage of reads supporting Gammaproteobacteria went from 63.4 to 58.2%. For these numbers to register as a change, 01-2F would have to be greater than 73.4% or less than 53.4%. On the other hand, within the same transition, the percentage of reads supporting Bacilli went from 27.1 to 10.9%. These numbers were different by an absolute 16.2%. Therefore, there was a "change". With this definition, 29 of 129 transitions (23%) register as "no change". An example is the 02-2F to 02-3F transition (Table  S2; http://med.stanford.edu/sgtc/research/skin_microbiome. html), where the largest change is for Actinobacteria, which went from 17.2 to 26.2%, an absolute change of 9%. There was no series of 01-to-02-to-03 where there was no change. We considered a change of an absolute 20% or more as a large change: again, an arbitrary, but reasonable amount that we have employed previously (Hyman et al., 2012). There were 47 out of 129 transitions (36%) that were more than an absolute 20%.  There were three cases where there was more than an absolute 20% change in going from swab 1 to swab 2 to swab 3: e.g., 15-1F-to-15-2F-to-15-3F (Table S2; http://med.stanford.edu/sgtc/research/skin_microbiome.h tml). For this swab, the percentage of reads supporting the presence of Gammaproteobacteria went from 49.3to-17.4-to-28.0%. There was no transition of greater than an absolute 30%.

AJM
The 129 transitions were divided into four groups: 1F-to-2F, 2F-to-3F, 1B-to-2B and 2B-to-3B. The changes in the composition of the three most abundant bacterial Classes were examined individually for the four types of transitions ( Table 5). The percent of Actinobacteria and Bacilli seldom changed (up or down), while the percent of Gammaproteobacteria decreased by more than 50% in all four transitions ( Table 5).

Order
The bacteria on 67 neck skin swabs (29% of the total swabs) were classified further into Order: 27 (40%) front swabs and 40 (60%) back swabs (Table S3; http://med.stanford.edu/sgtc/research/skin_microbiome.ht ml). Presumably, our data and the RDP software were the determinants for which samples the bacteria could be identified by Order. As for the Class data, to make the Order data visually comprehensible, only those bacterial Orders supported by, at least, 1% of the reads are reported (Table  S3; http://med.stanford.edu/sgtc/research/skin_microbiome.ht ml). Thereby, 31 Orders of bacteria were identified, compared to 15 Classes of bacteria. (More Orders than Classes were expected because each Class is composed of more than one Order). Only one bacterial Order, Actinobacteridae, was supported by, at least, 1% of the reads on all of the 67 swabs ( Table 6).
The average percentage of each bacterial Order supported by more than 1% of the reads on each neck swab is given in Table 6. Actinobacteridae were supported by the most reads for 58% of the swabs, followed by Enterobacteriales (29%), followed by Bacillales (10%), followed by Burkholderiales and Flavobacteriales (both at 1.4%). (Unfortunately, there were only three swabs with Order data for the third time point: two Front swabs and one Back swab. Therefore, for our further analyses, we ignored the Order data for the third time point.).
The average percentage of sequence reads for each swab supporting a bacterial Order by more than 1% of the reads is given in Table 7. None of the values was statistically significantly different from the others. For example, the value for 1F (95.0+/-1.0%) was not statistically significantly different from the value for 2F (96.0+/-2.4%; p = 0.27). The overall average percentage (n = 67) of each bacterial Order supported by more than 1% of the reads was 95.7+/-1.8%. This percentage was statistically significantly different (p<0.00001) from the analogous average for the Class data. This difference was not surprising as the reads were spread over 31 Orders compared to 15 Classes.
The SDI was calculated for the Order data for each swab ( Table 8). The average values and standard deviations for each group of swabs were also determined and reported in Table 8. The average SDI for 1F was not statistically different from the average SDI for 1B (p = 0.10). The average SDI for 2F was not statistically different from the average SDI for 2B (p = 0.52). In contrast, the average SDI of 1F was statistically significantly different from the average SDI of 2F (p<0.00001) and the average SDI of 1B was statistically significantly different from the average SDI of 2B (p<0.00001).
There were only eight transitions within the Order data: two 1F to 2F (volunteers 05 and 35) and four 1B to 2B (volunteers 07, 16, 34, 35, 37 and 38). For six of the eight transitions, the percentage of Actinobacteridae increased substantially, while for all eight transitions the percentage of Enterobacteriales decreased substantially.

CONCLUSION
The skin microbiome at the base of the front and back of the neck was composed principally of the same five Classes of bacteria over the course of one year: Actinobacteria, Alphaproteobacteria, Bacilli, Betaproteobacteria and Gammaproteobacteria. We could not distinguish the front microbiome from the back microbiome, nor men from women. However, as demonstrated by SDI and PCA analyses, the skin microbiome at each of the three time points differed from the other two time points. That is, each time a neck swab was taken, the microbiome was unique, not as to the bacteria present but rather the relative percentages of the five principal Classes of bacteria. However, it must be noted that the mix of volunteers was different at each time point.
We chose to investigate the skin microbiome at the front and back of the base of the neck because no one else had published the skin microbiome at these two body sites, nor did we know of anyone else investigating those two sites. That makes comparisons of our results to published data difficult. Nevertheless, some comparisons may be made.
Our two studied skin sites may be classified as sebaceous sites (Kong, 2011). The closest previously studied sebaceous sites are the upper chest and the upper back. (The occiput, the back of the scalp, is not comparable to the nape, as our back swabs were taken below the hairline and the occiput is above the hairline). The microbiome of the upper chest is composed almost entirely of Actinobacteria (Kong, 2011). The microbiome of the upper back is composed principally of Actinobacteria with some Proteobacteria (Grice et al., 2009). These comparisons again emphasize that small differences in skin location may give rise to substantial differences in the composition of the skin microbiome.
There were two important technical limitations in our data collection. The first limitation was the universality of the PCR primers. The PCR primers employed to amplify the V6 hypervariable region of the rDNA may have some mismatches to the rDNAs of relevant bacteria (Soergel et al., 2012;Frank et al., 2008). Presumably, such mismatch resulted in less amplification, especially in the critical first round. The second limitation was the depth of sequencing. A minimum of 1% of the sequence reads supporting the presence of a bacterium was required to report the presence of that bacterium. The Chao1 analysis (Fig. 1) suggested that there were bacteria present at concentrations below this limit. In addition, scientists working within the Human Microbiome Project have reported some very idiosyncratic results from employing rDNA sequence to identify bacteria (JCHMPDGWG, 2012;Haas et al., 2011).