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

Background Skin microbiome varies from person to person due to a combination of various factors, including age, biogeography, sex, cosmetics and genetics. Many skin disorders appear to be related to the resident microflora, yet databases of facial skin microbiome of many biogeographies, including Thai, are limited. Methods Metagenomics derived B-RISA and 16S rRNA gene sequencing was utilized to identify the culture-independent bacterial diversity on Thai male faces (cheek and forehead areas). Skin samples were categorized (grouped) into (i) normal (teenage.hea) and (ii) acne-prone (teenage.acn) young adults, and normal (iii) middle-aged (middle.hea) and (iv) elderly (elderly.hea) adults. Results The 16S rRNA gene sequencing was successful as the sequencing depth had an estimated >98% genus coverage of the true community. The major diversity was found between the young and elderly adults in both cheek and forehead areas, followed by that between normal and acne young adults. Detection of representative characteristics indicated that bacteria from the order Rhizobiales, genera Sphingomonas and Pseudoalteromonas, distinguished the elderly.hea microbiota, along the clinical features of wrinkles and pores. Prediction of the metabolic potential revealed reduced metabolic pathways involved in replication and repair, nucleotide metabolism and genetic translation in the elderly.hea compared with that in the teenage.hea. For young adults, some unique compositions such as abundance of Propionibacterium acnes and Staphylococcus epidermidis, with a minor diversity between normal and acne skins, were detected. The metabolic potentials of the acne vs. normal young adults showed that teenage.acn was low in many cellular processes (e.g., cell motility and environmental adaptation), but high in carbohydrate metabolism, which could support acne growth. Moreover, comparison with the age-matched males from the US (Boulder, Colorado) to gain insight into the diversity across national biogeography, revealed differences in the distribution pattern of species, although common bacteria were present in both biogeographical samples. Furthermore, B-RISA served as a crosscheck result to the 16S rRNA gene sequencing (i.e., differences between teenage and elderly microbiota). Conclusions This study revealed and compared the microbial diversity on different aged Thai male faces, and included analyses for representing the bacterial flora, the clinical skin characteristics, and comparison with the US age-matched. The results represent the first skin microbiota of Thai males, and helps the design of a large-scale skin microbiome study of Thais. The findings of the diversity among ages, skin type and national biogeography supported the importance of these traits in the skin microbiome and in developing a safe and sustainable treatment for acne and aging skin diseases.

skin microbiome study of Thais. The findings of the diversity among ages, skin type and national biogeography supported the importance of these traits in the skin microbiome and in developing a safe and sustainable treatment for acne and aging skin diseases.

INTRODUCTION
Following advances in metagenomics and next generation sequencing (NGS), the Human Microbiome Project (HMP) was established and several aspects of human microbial symbionts have been deciphered (HMP, 2012;Cosseau et al., 2016). Different organs and parts of the human body harbor an array of microbiota according to specific niche characteristics (Zhou et al., 2013). The anatomy and physiology are key factors in determining the skin bacterial diversity, such as the forehead, forearm, palm, finger, back of knee and foot (Costello et al., 2009;Zhou et al., 2013). Even on a particular niche of the body (e.g., face), the skin microbiota are still complicated by a combination of both internal and external factors, including but not limited to, age, sex, biogeography (climate and environment), genetics, cosmetics, diet, immune response, hormones and daily lifestyles (Cho & Blaser, 2012;HMP, 2012).
Generally, skin bacteria are dominated by members of the Actinobacteria, Proteobacteria, Firmicutes and Bacteroidetes phyla. Within these phyla, genera in Propionibacterium, Corynebacterium, Staphylococcus, Micrococcus and Brevibacterium are common (Van Rensburg et al., 2015;Mukherjee et al., 2016). The human skin bacteria are initiated since the different modes of delivery. For instances, Lactobacillus and Prevotella are predominant on the infant skin born naturally through the vagina, whereas Acinetobacter and Staphylococcus are predominant on those born by caesarian section (Dominguez-Bello et al., 2010). Bacterial diversity continues to vary through time, due in part to the internal and external combinatory factors of each individual (Capone et al., 2011;SanMiguel & Grice, 2015). Within the same age range, the skin type, such as moist and sebaceous skin vs. dry skin, also plays a role in the bacterial diversity (Grice et al., 2009). Mukherjee et al. (2016) analyzed the correlation between the sebum and hydration levels on the cheek and forehead microbiota of 30 Indian female volunteers, and found the sebum and hydration levels provided predictive power with respect to the nature and microbiota of forehead and cheek skins. This is consistent with earlier studies, which reported that the microbial compositions play roles in the ''normal,'' ''predisposed'' and ''disease'' states of the skin health, such as in atopic dermatitis and chronic skin ulcers (Costello et al., 2009;Cho & Blaser, 2012;Van Rensburg et al., 2015). It was hypothesized that a loss of protective bacteria led to the outgrowth of pathogenic bacteria (Van Rensburg et al., 2015;Barnard et al., 2016). Kong et al. (2012) reported that children with a severe level of atopic dermatitis had an increased level of Staphylococcus aureus but a decreased overall microbial diversity on the skin. Recently, based on 38 US acne patients and 30 US age matched controls, Barnard et al. (2016) found that the acne and health of the human facial skin was not associated with an abundant Propionibacterium acnes but with a disrupted microbiota. Accordingly, they proposed using probiotic and P. acnes phage as sustainable and safe strategies for maintaining skin health.
Possible reasons for the disrupted microbiota to cause acne can be many, such as diet habits (i.e., glycemic load and water intake), hours of quality sleep (i.e., affects emotional state), and genetics. This relationship is through the ''gut-brain-skin theory'' (Bowe, Patel & Logan, 2014). A healthy eating pyramid (types and proportions of food for a daily healthy eating), such as Mediterranean food that is rich in vegetable fibers, fruit antioxidants, vegetable and fish proteins, and low glycemic load, could change the gut microbiota (Wu et al., 2011;Queipo-Ortuno et al., 2012;Lopez-Legarrea et al., 2014). Clark, Haas & Sivamani (2017) reviewed the links between eating plant fibers and low glycemic food with altered gut microbiota and acne improvement. Another example was that the oral consumption of lactic acid bacteria has been shown to alleviate acne and mental stress (i.e., anxiety and depression). Lactic acid bacteria improve bowel function and reduce constipation, which is one underlying cause of acne and depression, and further, the diets and gut microbiota affect inflammatory signals in the human body (Bowe, Patel & Logan, 2014;Lopez-Legarrea et al., 2014). Other influence factors to acne include daily intake amount and types of beverages (e.g., milk vs. drinking water), and hours of sleep. Adebamowo et al. (2005) reported an association between teenage acne and milk (whole milk and skim milk) consumption, and Rigopoulos et al. (2007) reported that adequate hours of sleep could, in part, reduce teenage acne.
The present study utilized culture-independent metagenomics combined with bacterialribosomal intergenic spacer analysis (B-RISA) and 16S rRNA gene sequencing to disclose the bacterial communities of normal skin type (healthy face skin as inspected by dermatologist) Thai males of different age periods. These were categorized as young adults (teenage.hea), middle-aged (middle.hea) and elderly (elderly.hea) and, for young adults, the presence of mild-level acne vulgaris was also included (teenage.acn). Acne vulgaris, also called teenage acne, is a common skin disease among teenage and young adults. Approximately 80% of teenagers and young adults have or experienced this skin problem, because this age period is more prone to over active sebaceous glands, which cause plugged pores that promote the overgrowth of Propionibacterium acnes (James, 2005).
As for the methods, B-RISA served a pre-classification and cross-validation of the bacterial diversity without sequencing. B-RISA, denatured gradient gel electrophoresis (DGGE), and terminal restriction fragment length polymorphism (TRFLP) are established methods that classify bacterial diversity structures using electrophoresis technique (Fisher & Triplett, 1999;Danovaro et al., 2006;Sikaroodi & Gillevet, 2012;Gobet, Boetius & Ramette, 2014;Jami, Shterzer & Mizrahi, 2014). B-RISA relies on the characteristic DNA fragment lengths between the 16S and 23S rRNA genes upon bacterial species. Different species have different B-RISA lengths. Researches have compared the power of B-RISA and NGS in classifying bacterial community patterns, and reported the agreement in the term of the discriminating power among samples. Nevertheless, B-RISA yielded a much underestimated number of the species richness and provided no information regarding sequences or names of individual species (Grillevet, Sikaroodi & Torzilli, 2009;Sikaroodi & Gillevet, 2012;Gobet, Boetius & Ramette, 2014;Jami, Shterzer & Mizrahi, 2014;Dehingia et al., 2015). This study adopted B-RISA as an inexpensive and rapid bacterial community fingerprint tool prior to NGS.
For the 16S rRNA gene sequencing to obtain the bacterial taxonomic compositions, the study utilized mothur's MiSeq standard operating procedures (SOP) (Schloss et al., 2009) and Metastats (White, Naharajan & Pop, 2009) to detect differentially abundant features, including representative bacterial species and clinical features (skin, diet and lifestyle), in the datasets. The study also compared the Thai male face skin microbiota with age-matched US biogeography (reference data from Costello et al., 2009), and the bacterial diversity and metabolic potentials between these two biogeographic regions were reported. Together, the skin microbiota variation according to age, skin features, lifestyles, and biogeography, and thus various aspects of the skin, will help to complete the skin microbiota knowledge and databases (Van Rensburg et al., 2015).

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). Notes. a Fitzpatrick skin types I to VI are defined as follows: 'I', means always burns, never tans (pale white skin); 'II', means always burns easily, tans minimally (white skin); 'III', means burns moderately, tans uniformly (light brown skin); 'IV', means burns minimally, always tans well (moderate brown skin); 'V', means rarely burns, tans profusely (dark brown skin); and 'VI', means never burns (deeply pigmented dark brown to black skin). b Degree of roughness, wrinkle, mottled hyperpigmentation, flushing, pore, acne, and oil zone range from 0 (none) to 1 (maximal degree); except elasticity that ranges from −1 (minimal degree) to 0 (normal). All data represent the average ± standard deviation (S.D.) from 10 subjects per group. Note healthy daily food pyramid intake and below questions the data were collected as 'yes' or 'no', hence the percentage was calculated and no S.D. is available.

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.

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 R 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 R 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 (  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 260 /A 280 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.
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).

DISCUSSION
Metagenomic approaches, such as B-RISA and NGS, allow the in-depth evaluation and database generation of the microbial diversity, without the limitation of cultivability. The high Good's coverage indices indicated the sufficiency of the NGS (16S rRNA gene) sequencing reads that covered over 98% of the diversity at the genus level (Table 3) (Kennedy et al., 2014). The small amount of unclassified bacterial phylum, less than 1.6% in average ( Fig. 1), was consistent with previous studies. As unclassified phyla are often denoted as derived from environmental exposure, the data inferred no significant contamination from environmental exposure or sequencing error (Salter et al., 2014). Our studies, therefore, for the first time, revealed the skin bacteria of Thai males; the data by B-RISA and 16S rDNA profiles elucidated the community difference between skin ages and types. The B-RISA finding was consistent with our NGS data, and revealed a correlation between the cheek and forehead subgroups, whereas a significant variation among the skin ages and types (Table 2 and Fig. 2). The higher bacterial nucleic acid yield (Table 2), community alpha diversity (Table 3) and number of GLOTUs (Fig. 3) in the elderly might support a previous report on increased microbial frequency in dry skin areas (Grice et al., 2009). The clinical records confirmed that the elderly.hea had the relatively lowest hydration and sebaceous levels (Table 1: average oil zone of 0.00). Skin characteristic analysis using Metastats detected the oil zone and acne as major influences to the bacterial community on the young adult skin (Fig. 4C).
Detecting representative bacterial species among the skin groups using Metastats indicated Rhizobiales and Sphingomonas, and the wrinkles and pore skin features, within the elderly.hea (Figs. 4B and 4C). On the other hand, the teenage.acn, and the middle.hea.foreheads2 that was clustered with the teenage.acn, shared the signatures of P. acnes and S. epidermidis. Additionally, sequencing the 16S rRNA V3-V4 region allowed for hierarchical taxonomic classification to the species level of bacteria (Somboonna et al., 2014;Mukherjee et al., 2016). The teenage.acn also had prevalent S. aureus. These bacterial species are associated with acnes and other skin diseases (Krutmann, 2009;Roudsari et al., 2015). However, P. acnes normally reside in the follicles around sebaceous glands, for which the special sampling method of a deep cleansing pore strip is required to sample (Barnard et al., 2016). Thus, the prevalence of P. acnes found in this study merely represents the general abundance on the superficial cheek and forehead areas. Nevertheless, superficial swabs used to sample the skin in this study are widely accepted and therefore supported our current research, exploring the bacterial diversity, while also providing insight about our future goal of helping to resolve the healthy bacteria balance of facial skin (Ling et al., 2013;Barnard et al., 2016;Mukherjee et al., 2016).
The microbial metabolic potentials were different between the elderly and young adult skin, and between the healthy and acne young adult skin. Compared with the elderly, the KEGG functions involved in diseases were less prevalent but the cell replication and repair was more prevalent in the normal young adults (Fig. 6). Several unclassified KEGG pathways (level 1) and functional groups (level 2) reflect the limited KEGG database, and hence limit the interpretation of the metabolic potentials by this method.
Ying et al. (2015) studied age (young adult vs. adult vs. elderly) and geography (urban vs. rural) factors associated bacterial communities on seven skin sites and 71 healthy volunteers. While body sites showed variance in the bacterial community structures, the communities also significantly differed by age and geographic factors. The database even derived the bacterial markers and allowed an opportunity to use microbiota data to predict whether an individual is a rural or urban resident (Ying et al., 2015). Skin microbiota changes over lifetime, not only by body sites and geography (Chu et al., 2017), but also the inflammatory states which are directly and indirectly affected by ages, diets, lifestyle, clinical problems, and etc. (Zapata & Quagliarello, 2015). Following the ''gut-brain-skin theory,'' the diverse gut microbiota in differing age groups support the diversity of the skin microbiota, consistent with the US male cohort that were found different skin bacterial communities among ages (Costello et al., 2009), and our recent findings in Thai females that found age is a factor that brings the bacterial community relationship factor closer, i.e., normal and acne teenage skin vs. normal elderly skin (Somboonna et al., 2017). Further, we re-subsampled the data at the sequencing depth of 8,000 per subgroup to obtain two replicate datasets. Analyzing these data for relative community structure and composition relatedness (Figs. S3A and S3B) and representing GLOTUs (Figs. S4B and S4E) and clinical features (skin, diet and lifestyle) (Figs. S4C and S4F) exhibited consistent findings. Together, these supported the bacterial community variation by age and relevant skin features.
For healthy eating diet, we did not detect any significant correlation, which might correspond to general healthy eating diets among volunteers (Table 1: 80%-90%). The healthy eating pyramid defines the appropriate types and proportions of daily food consumption (e.g., food that is rich in vegetable fibers, fruit antioxidants, vegetable and fish proteins, and low glycemic load). The lack of variation in healthy eating pyramid factor across our study groups might have obfuscated the Spearman's correlation findings, as we did not find the healthy eating pyramid related with any microbiota structures ( Fig. 4C and Figs. S4C and S4F). To investigate the diet factor to the skin microbiota, more detail and variations of eating habits, such as glycemic load, should be obtained and analyzed (Wu et al., 2011;Lopez-Legarrea et al., 2014). To date, oral and topical probiotics have been proposed a healthy skin therapy (Isawa et al., 2008;Kimoto-Nira et al., 2014;Huang & Tang, 2015;Kober & Bowe, 2015); our ongoing studies include tracking the effect of oral and topical probiotics on skin microbiota. Different facial skin bacterial structures between Thai and previously reported US males (Costello et al., 2009) is in accord with previous research that suggested the skin bacteria vary across ethnicities and geographies; diet, genetic and environment factors are possible causes of this difference (HMP, 2012;Zhou et al., 2013). Here, extremophiles, such as Firmicutes, were more enriched on the Thai skin compared to the US, whereas Actinobacteria were more prevalent on the US skin (Fig. 5). This difference might potentially be involved in environmental adaptation. In support, common bacteria on dry skin were previously reported to include Proteobacteria, Firmicutes and Bacteroidetes (Grice & Segre, 2011;Staudinger, Pipal & Redl, 2011). Of additional note are the variation of V3-V4 vs. V2 sequences annotation, and that we sampled skin at ∼8 h after face wash, which more hours allowed for natural bacterial flora establishment (Costello et al., 2009).

CONCLUSIONS
As sebum and hydration levels represent two significant predictors to the nature and diversity of the skin, this study firstly surveyed the facial skin microbiota of Thai males that differ in age and acne skin groups. In addition, the facial skin microbiome of elderly individuals has rarely been reported. Using young adults as the reference point, this study demonstrated that the facial skin bacteria population structure continued to change from young to middle age to elderly adulthood, in association with specific bacterial flora and clinical skin features. The predicted KEGG metabolic profiles of the microbiota supported the functional diversity among the different skin types. Comparison between Thai and US biogeographies showed the bacterial community differences that could be related with specific skin characteristics. Together, our findings on bacterial community structure characterizing Thai male facial skin represents the first preliminary findings to be used in other comparative studies in order to gain insights into healthy facial skin microbiota, with possible correlation with biogeography, gender, and skin types.