Diversity of bacterial communities on the facial skin of different age-group Thai males

Sample collection

Study subjects were comprised of 40 healthy Thai male volunteers with no health disease and skin disorders, except for the presence of acne vulgaris, where stated. Subjects were classified as 10 normal face young adults (19–24 y; teenage.hea), 10 normal face middle-aged (32–38 y; middle.hea), 10 normal face elderly (51–57 y; elderly.hea) and 10 mild grade acne young adults (19–24 y; teenage.acn). Skin types were graded with classification support from a dermatologist (Table 1). Volunteers were recruited at Chulalongkorn University and Phramongkutklao Hospital, in Bangkok, during June–September 2013. The sample selection criteria and collection process followed the established procedures (Staudinger, Pipal & Redl, 2011; Ling et al., 2013). All samples were swabbed on the same sample areas of the face by the same investigator, and in the same swabbing manner to prevent possible variation across samples. Participants had used neither topical antibiotics in the past 7 d nor oral antibiotics in the past 3 d. Participants were refrained from washing their faces for approximately 8 h prior to sample collection. To minimize sample cross-contamination, subjects belonging to different skin groups were instructed to come for sample collection on a different week. For each individual, an area of 1 × 1 inch was swabbed on the forehead, as well as on both the right and left cheek. With the pooling of the two cheek swabs, this led to two clinical samples per volunteer: forehead and cheeks, and that two areas of the cheek (left and right) was taken into account in the estimates of DNA density per swab area. With 10 individuals per group and four groups in total, the study thereby contained 80 samples. Each skin swab sample was immediately stored in sterile saline-Tween solution (0.15 M NaCl and 0.1% (v/v) Tween 20) at −20 °C, and processed within 14 d. The study was approved by the Institute Review Board Royal Thai Army Medical Department Office, Phramongkutklao College of Medicine, Bangkok, Thailand (Q003h/56).

Total bacterial DNA extraction and DNA quality examination

Total bacterial DNA was extracted from each sampling using GF-1 Bacterial DNA Extraction Kit (Vivantis, Selangor Darul Ehsan, Malaysia), according to the manufacturer’s protocol. The quality and quantity of the extracted DNA were determined by NanoDrop spectrophotometry and agarose gel electrophoresis. This DNA was used for B-RISA and 16S rRNA MiSeq.

B-RISA

The methods of B-RISA, including primers 16S-1392F and 23S-125R (Table S1) and thermocycling parameters (95 °C for 4 min, and 35 cycles of 94° 1 min, 55 °C 1 min and 72 °C 10 s, followed by 72 °C 2 min), followed established protocols (Fisher & Triplett, 1999; Hewson & Fuhrman, 2004; Danovaro et al., 2006). B-RISA targets an intergenic region between the 16S and 23S rRNA genes, where different bacteria specie (or group) is denoted the different size length based on electrophoresis technique. Thus, the community of bacteria are analyzed by the different size lengths composition (Gobet, Boetius & Ramette, 2014; Jami, Shterzer & Mizrahi, 2014; Dehingia et al., 2015). Each 20 µl B-RISA reaction comprised 1× EmeraldAmp^® GT PCR Master Mix (TaKaRa, Shiga, Japan), 0.3 µM of each primer, and 50 ng of the total bacterial DNA. Independent replicate reactions were performed per samples. The samples were analyzed for the different size lengths compositions on 1.75% agarose gel electrophoresis.

16S rRNA gene library preparation and MiSeq sequencing

Universal prokaryote primers 338F and 803R for 16S rRNA gene V3–V4 with the appended eight-nucleotide barcode on the 5′ tail of the oligo sequences were listed in Table S1 (Meyer, Stenzel & Hofreiter, 2008; Somboonna et al., 2014; Somboonna et al., 2017). The PCR recipe and thermocycling parameters followed established protocols: 25-µl PCR reaction comprised 1× EmeraldAmp^® GT PCR Master Mix (TaKaRa, Shiga, Japan), 0.3 µM of each primer, and 50–100 ng of the total bacteria DNA; and the PCR conditions were 95 °C for 4 min, and 30 cycles of 94°C for 45 s, 50 °C for 55 s and 72 °C for 1 min 30 s, followed by 72 °C for 10 min (Meyer, Stenzel & Hofreiter, 2008; Somboonna et al., 2017). To prevent PCR stochastic bias, the template quantity and quality was adequate, and a minimum of three independent PCR reactions were performed per sample (Kennedy et al., 2014). The negative controls included a separate mock swabbing sample and no template control in the PCR, both of which showed no amplification. The 16S rDNA amplicons (∼466 bp) were purified using PureLink^® Quick Gel Extraction Kit (Invitrogen, New York, USA), quantified by Qubit (Life Technologies, Carlsbad, CA, USA), and determined a purity and size of the purified amplicon by Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Due to the heterogeneity of some samples of the groups as was demonstrated by B-RISA banding (Table 2), DNA from same-group subjects that shared similar B-RISA profiles were pooled together in equal amplicon amount generating a total of 11 subgroups for Illumina MiSeq sequencing. MiSeq sequencing adaptors were added to each of the 11 subgroups using Nextara Index Kit (Illumina, San Diego, CA, USA), and purified using Agencourt AMPure Technology (Beckman Coulter, Brea, CA, USA). The products were quantified and size checked by Qubit (Life Technologies) and Agilent 2100 Bioanalyzer (Agilent Technologies), respectively. The libraries were also quantified by quantitative PCR using Kapa Biosystems Library Quantification Kit (Kapa Biosystems, Woburn, MA, USA). Paired-end sequencing, 2 × 250, was performed using the Illumina MiSeq platform following the manufacturer’s protocols at Chulalongkorn Medical Research Center (Bangkok, Thailand).

Availability of DNA sequencing reads

All nucleic acid sequences in this study were deposited at the NCBI Sequence Read Archive (SRA) database (accession number SRR3656927). The data pertaining to each stage of the following bioinformatic analyses were available at https://doi.org/10.6084/m9.figshare.5288338 (for pre-processing data, and scripts in the phylotype pipeline) and https://doi.org/10.6084/m9.figshare.5288332 (the results for OTU classification data).

Bacterial composition and metabolic potential analyses

The raw sequences were categorized into groups based on the 5′ barcode sequences. The sequences were processed following mothur’s MiSeq SOP (Schloss et al., 2009; Kanehisa et al., 2016). The pre-processing steps included removal of (i) short read lengths of ≤100 nucleotides (excluding the primer and adaptor sequences), (ii) long homopolymers of ≥eight nucleotides, (iii) ambiguous nucleotides and (iv) chimera. The clean sequences were classified to the operational taxonomic unit (OTU) using the Ribosomal Database Project (RDP) Classifier (Wang et al., 2007) and Greengenes database (McDonald et al., 2012). A minimum bootstrap confidence score of 80% was used as a cutoff for taxonomic assignment. Genus and specie of OTU (GLOTU and SLOTU) followed the phylotype-based methods (Schloss & Handelsman, 2004). Alpha diversity by Shannon and Chao bacterial community richness, and the Good’s coverage index to estimate the data coverage of a community, were computed using mothur (Claesson et al., 2009; Schloss et al., 2009). Data normalization was performed to normalize the varying sequencing depth among 11 subgroups, using the mothur software (Schloss et al., 2009). For the groups with two B-RISA subgroups (Table 2: e.g., teenage.acn.foreheads1 and teenage.acn.foreheads2), the Student’s t-test (two-sided) and Boxplot variance analyses for these pair subgroups were computed. The data were merged if p > 0.05 (Table 3: i.e., teenage.acn.foreheads). The relative abundance of bacterial genera was visualized as Heatmap using the R statistics package. Additional mothur analyses included Venn diagram, and the community dissimilarity matrices by Morisita-Horn dissimilarity index (Schloss et al., 2009). The average-linkage dendrogram, used to cluster the community relationship among subgroups, was computed from the Morisita-Horn dissimilarity indices using the R statistics package. Non-metric multidimensional scaling (NMDS) along with analysis of molecular variant (AMOVA) statistics were computed using mothur, and visualized by R. Differentially abundant (representative) GLOTU and clinical features were detected using function Metastats in mothur with default parameters (White, Naharajan & Pop, 2009). The microbial metabolic potentials were predicted using PICRUSt software (Langille et al., 2013), which estimates the community metabolic potentials from the 16S rDNA sequencing data using the KEGG (Kyoto Encyclopedia of Genes and Genomes) database (Kanehisa et al., 2004; Kanehisa et al., 2016). Comparison of metabolic potentials between a pair of groups was computed in PICRUSt (Langille et al., 2013) based on a two-sided Welch’s t-test with adjusted Storey-Tibshirani’s (Storey & Tibshirani, 2003) and Benjamini–Hochberg’s (Benjamini & Hochberg, 1995) false discovery rates. Significant difference was accepted when the q-value was <0.05.

For the Thai vs. US biogeography microbiota comparison, we compared the data with the culture-independent bacterial diversity databases of 6 normal age-matched male forehead-skin resided in Boulder, Colorado, USA (Costello et al., 2009): 5 males aged 30–35 y to compare with our middle.hea and 1 male aged 60 y to compare with our elderly.hea. The methods of sample collection and library preparation were similar: used 16S rRNA amplicon sequencing, except that Costello et al. (2009) did not require hours between bathe and sampling (ranged from 0.5 to several hours) vs. we sampled at ∼8 h (similar to Staudinger, Pipal & Redl, 2011; Ling et al., 2013, to allow establishment of natural flora on the volunteer faces), amplified the V2 region, which might overestimate species richness while V3–V4 region generally provides species richness comparable to the nearly full-length 16S rRNA gene (Youssef et al., 2009), and utilized 454 FLX sequencer. The Costello et al. (2009) sequencing data were downloaded and processed bioinformatic analyses the same procedures as our data for parallel data comparison.

Physical and descriptive skin-related features and bacterial DNA extraction

The physical and descriptive skin-related characteristics of individuals were recorded, and the average results for the teenage.hea, teenage.acn, middle.hea and elderly.hea groups were displayed in Table 1. Among teenage.hea, middle.hea and elderly.hea, an increased intensity of skin roughness, wrinkles, mottled hyperpigmentation and pores were recorded with age, whereas skin flushing and the oil zone decreased. The acne vulgaris of the teenage.acn group was of a mild level, covering 40 ± 26% of the entire face and none on the trunk (Table 1). To control for contamination, a separate mock swabbing of clean air inside the biosafety cabinet (Van Rensburg et al., 2015), was included and no PCR amplicon occurred. The diet and lifestyles that might affect acne were also recorded. Participants reported generally healthy eating pyramids (80–90%) and no makeup use (Table 1).

Following bacterial DNA extraction of all 80 clinical forehead and cheeks samples, the low yield of bacterial DNA (Table 2: 0.73–4.04 ng/µl) indicated the light number of bacteria on the facial skin in nature (Costello et al., 2009; Zhou et al., 2013). The amount of extracted bacterial DNA was mostly greater from the forehead swabs than the cheeks, and from the elderly than the teenage skins (Table 2), where teenage.hea foreheads yielded 133 ng/swab area, middle.hea foreheads 159 ng/swab area and elderly.hea foreheads 404 ng/swab area. Some samples had an unsatisfactory DNA quality and so could not be successfully amplified and analyzed by B-RISA and by 16S rRNA gene. These were 10 samples from middle.hea cheeks, due to an average A₂₆₀/A₂₈₀ ratio of 1.10, and one sample each from elderly.hea cheeks and forehead, leaving nine out of 10 samples for this group (Table 2).

Intragroup correlation and intergroup variation of bacterial communities by B-RISA

Of individual cheeks and foreheads from the same group showed a degree of relatedness in the B-RISA electrophoresis results, with uniform or just two banding patterns observed among the corresponding samples of cheek and forehead subgroup. This supported bacterial community similarity among intragroup samples, and led to 11 B-RISA classified subgroups (Table 2). The B-RISA results suggested the different bacterial community structures among the four groups (teenage.hea, teenage.acn, middle.hea and elderly.hea), and that the community structure of the teenage.acn foreheads were somewhat similar between the teenage.hea and middle.hea groups (Table 2).

Comparative taxonomic profiles by 16S rRNA gene sequencing

The 16S rRNA gene hypervariable V3–V4 sequencing was successful, as the sequencing depth of all 11 subgroups were sufficient to satisfy the >98% Good’s coverage index, and the taxonomic compositions could be annotated (Table 3). The V3–V4 region is accepted as most capable of bacterial species identification and hence of worldwide use (Youssef et al., 2009; Kennedy et al., 2014). Of groups that the B-RISA results suggested each consisted 2 subgroups (Table 2: teenage.hea foreheads, teenage.acn cheeks and foreheads, and middle.hea foreheads), the Student’s t-test (two-sided) and Boxplot variance of the 16S rDNA microbiota between the corresponding pairs were analyzed and the p statistics supported the community diversity in the teenage.hea foreheads1 and foreheads2, teenage.acn cheeks1 and cheeks2, and middle.hea foreheads1 and foreheads2 were significantly different (p ≤ 0.05) and so were retained as separate subgroups, whereas the teenage.acn foreheads1 and foreheads2 were not significantly different (p = 0.17) and so were merged to the single “teenage.acn.foreheads” subgroup. For the teenage.acn foreheads1 and foreheads2, the Boxplot analysis showed the overlapped variance, supported the merging of the “teenage.acn.foreheads” (Fig. S1). This led to 10 separate subgroups in the following analyses (e.g., Table 3 and Fig. 1).

The alpha diversity at the genus, family and order hierarchical levels of OTUs (GLOTU, FLOTU and OLOTU, respectively) of the 10 subgroups, using Chao’s taxon richness and Shannon’s diversity evenness (combining taxon richness and evenness), and the estimated sequencing coverage index (Good’s coverage), were summarized in Table 3. The community was comprised of 109–167 GLOTU, 69–101 FLOTUs and 34–48 OLOTUs. The number of OTUs and sample diversity indices indicated the diversity was relatively high in the elderly. For example, elderly.hea.cheeks was ranked in the top two highest numbers of OTUs at all three taxonomic levels, and had the highest Chao’s genus richness (Table 3). Moreover, estimating the genus richness from the rarefaction curve of the number of GLOTUs (Y-axis) against the number of random subsampling sequences (X-axis) suggested that for the cheeks the highest genus richness was found in the elderly.hea.cheeks and teenage.acn.cheeks1, and in the foreheads was found in the elderly.hea.foreheads and teenage.hea.foreheads1 (Figs. S2A and S2B). Together, the diversity indices highlighted the microbial diversity in the elderly and teenage groups.

Intragroup correlation and intergroup variation of bacterial communities by 16S rRNA gene sequencing

The dominant phyla across the skin groups were Proteobacteria (37.9 to 86.2%) and Firmicutes (3.3% to 42.5) (Fig. 1A). These phyla have all been reported previously on the human face skin (Grice et al., 2009). The unclassified phyla averaged <1.6% (Fig. 1). Figure 2 showed the community distributions at the GLOTU level, in order of their community relatedness (clustering dendrogram on the left). For the young adults, abundant genera included Staphylococcus in phylum Firmicutes (elderly.hea 9.5%, middle.hea 23.4%, teenage.acn 19.6% and teenage.hea 17.6%), Propionibacterium in phylum Actinobacteria (elderly.hea 0.9%, middle.hea 7.0%, teenage.acn 23.3% and teenage.hea 1.8%) and Corynebacterium in phylum Actinobacteria (elderly.hea 2.4%, middle.hea 3.8%, teenage.acn 3.9% and teenage.hea 9.4%). For the elderly adults the abundant genera were an unclassified genus in the order Rhizobiales in class Alphaproteobacteria, phylum Proteobacteria (elderly.hea subgroups comprised 32.1% in average, middle.hea 31.3%, teenage.acn 17.5% and teenage.hea 15.3%) and Sphingomonas in phylum Proteobacteria (elderly.hea 6.4%, middle.hea 9.4%, teenage.acn 1.3% and teenage.hea 2.24%). Other prevalent bacteria included Pseudoalteromonas in phylum Proteobacteria (elderly.hea 13.9%, middle.hea 3.0%, teenage.acn 3.1% and teenage.hea 10.0%) and Stenotrophomonas in phylum Proteobacteria (elderly.hea 3.2%, middle.hea 8.3%, teenage.acn 3.7% and teenage.hea 2.7%) (Fig. 2).

Furthermore, the 16S rDNA V3–V4 sequences could phylotype Staphylococcus and Propionibacterium to the species level. Staphylococcus was found to mainly belong S. epidermidis (elderly.hea 3.13%, middle.hea 17.4%, teenage.acn 14.4% and teenage.hea 8.58%) and an unclassified OTU species (elderly.hea 5.54%, middle.hea 5.80%, teenage.acn 4.36% and teenage.hea 8.29%). The remaining species (<1%), in decreasing order of abundance, were S. saprophyticus, S. aureus, S. haemolyticus, S. lugdunensis, S. pettenkoferi, S. sciuri and S. succinus. Analysis at the species level found community differences between teenage.hea and teenage.acn, that the abundance of S. epidermidis was about two-fold greater in teenage.acn, and the third most common Staphylococcus species after the remaining unclassified in the teenage.acn was S. aureus, whereas that in the teenage.hea was S. saprophyticus. For Propionibacterium, the main species was P. acnes in teenage.acn (elderly.hea 0.76%, middle.hea 6.62%, teenage.acn 22.5% and teenage.hea 1.53%), while the rest (<1%) were an unclassified species and P. granulosum.

In summary, at the phylum level, Proteobacterium was prevalent in the elderly and middle aged, and Actinobacteria and Firmicutes in the young adults (both healthy and acne skins). At the genus level, Rhizobiales (unclassified genus) and Sphingomonas were prevalent in the elderly and middle-aged, Pseudoalteromonas in the elderly and normal young adults and Stenotrophomonas in the middle-aged and acne young adults. Propionibacterium (P. acnes) were prevalent in the acne young adults. Staphylococcus (S. epidermidis, S. saprophyticus, S. aureus and unclassified species) and Corynebacterium were prevalent in the middle-aged, and normal and acne young adults (Fig. 2). The similarity in the microbial population structures across groups caused the phylogenetic analysis (Fig. 2) to cluster cheeks and foreheads of the same corresponding groups (intragroup correlation), and young adults (in particular teenage.hea) separate from the elderly (intergroup variation).

The Venn diagram displaying number of shared and unique GLOTUs across groups highlighted that teenage.acn shared a higher degree of GLOTU similarity with middle.hea than with teenage.hea (Fig. 3: teenage.acn & middle.hea 282 of 375 = 75.2%, teenage.hea & middle.hea 249 of 375 = 66.4%). This accounted for 65.7% (282 of 429) of the teenage.acn genus compositions. Furthermore, non-metric multidimensional scaling (NMDS) of the community correlation among subgroups against the two-dimensional distance located the teenage.acn communities nearby to the middle age (Fig. 4A).

The GLOTU compositions were tested by AMOVA to determine the statistical difference of their variation. No statistical significant intravariation between the forehead and cheeks of the same group was found, but significant intervariation was found between elderly.hea and teenage.hea (p = 0.019), and between teenage.hea and teenage.acn (p = 0.014) (Fig. 4A).

Metastats along AMOVA statistics was performed to identify any particular species of bacteria responsible for the direction of the community structure. The statistics suggested that the unclassified genus in the order Rhizobiales (p = 1e − 6) and Sphingomonas (8e−4) were distinguished to the elderly.hea community, and in an opposite direction to that driven by Propionibacterium (2.5e−2) and Staphylococcus (7e−6) for the young adult communities. Pseudoalteromonas (2.5e−2) drove a partial effect to the elderly.hea and teenage.hea, and a negative effect to the teenage.acn communities. Corynebacterium and Stenotrophomonas had no significant statistic impact on the skin community of middle-aged volunteers (9.8e−2 for Corynebacterium and 0.12 for Stenotrophomonas) (Fig. 4B).

The intervariation in the facial skin bacteria between age-matched Thai and US biogeographies were then depicted. A significantly different microbial diversity was observed at the phylum level, both for the middle-aged and elderly age groups (Fig. 5). While the major phyla were all present in both biogeographies, members of Proteobacteria were more prevalent in Thai males (Fig. 5),which is the same pattern that differentiated Thai elderly.hea and middle.hea from Thai teenage.hea (Fig. 4A). In contrast, the US males had relatively abundant members of Actinobacteria, and these were even more dominant in the elderly (Fig. 5). Note that comparison of the teenagers was not performed since the age-matched data were not found.

Diversity correlation by skin characteristics

Metastats was also performed to analyze the microbial community compositions with the clinical features in Table 1. The statistics demonstrated significant Spearman’s correlation (p < 0.05) of the microbiota patterns with the wrinkles, pores and mottled hyperpigmentation for the elderly, and flushing, oil zone and acne for the teenage. The directions were similar to those driven by Rhizobiales (Sphingomonas and Stenotrophomonas) for the elderly.hea and middle.hea.foreheads1, and by Propionibacterium and Staphylococcus for the teenage groups (Figs. 4B and 4C). Note that skin roughness did not show a significant correlation (Fig. 4C; p = 0.13). Adequate drinking water intake and hours of sleeps per day also exhibited significant correlation (Fig. 4C).

Comparison of metabolic potentials

As each bacteria species possess specific KEGG pathways, the functional potentials of the bacteria community could be predicted (Langille et al., 2013). Significant differences in the predicted KEGG functional pathways were found primarily between elderly.hea vs. teenage.hea, followed by between teenage.acn vs. teenage.hea (Fig. 6). Examples of different microbial functional pathways in aging were metabolism, genetic information processing, and human diseases (Figs. 6A and 6B). With respect to the acne, examples of different KEGG metabolic pathways were cell motility, environmental adaptation and carbohydrate metabolism. While several of the depicted metabolisms were lower for the teenage.acn than the teenage.hea, the carbohydrate metabolism was higher (Figs. 6C and 6D).