Taxonomic profiling of skin microbiome and correlation with clinical skin parameters in healthy Koreans

The interest in skin microbiome differences by ethnicity, age, and gender is increasing. Compared to other ethnic groups, studies on the skin microbiome of Koreans remains insufficient; we investigated facial skin microbiome characteristics according to gender and age among Koreans. Fifty-one healthy participants were recruited, the facial skin characteristics of each donor were investigated, their skin bacterial DNA was isolated and metagenomic analysis was performed. The donors were divided into two groups for age and sex each to analyze their skin microbiomes. Moreover, we investigated the correlation between the skin microbiome and clinical characteristics. The alpha diversity of the skin microbiome was significantly higher in the elderly, and beta diversity was significantly different according to age. The comparative skin microbials showed that the genus Lawsonella was more abundant in the younger age group, and Enhydrobacter was predominant in the older age group. Staphylococcus and Corynebacterium were more abundant in males, while Lactobacillus was more abundant in females. Lawsonella had a negative correlation with skin moisture and brown spots. Staphylococcus and Corynebacterium both had negative correlations with the number of UV spots and positive correlations with transepidermal water loss (TEWL). Furthermore, Staphylococcus aureus had a negative correlation with skin moisture parameters.

www.nature.com/scientificreports/ be limited to a specific gender. Paolo et al. suggested that understanding the physiological changes in the skin according to gender can help derive cosmetic improvement methods to prevent skin aging 14 .
In this study, we focused on the facial skin microbiome of Koreans to investigate the differences by age and gender. Taxonomic differences in skin microbiota were compared and analyzed between the younger and older age groups and between male and female groups. In addition, comparing and analyzing the skin microbiome and clinical information of healthy individuals was intended to lay a foundation for product development utilizing the skin microbiome.

Methods
Sample collection. Fifty-one healthy volunteers were recruited for collecting skin samples to investigate their skin microbiome characteristics. Recruitment was conducted only among those who signed written consent for the collection of material of human origin. The selection criteria were applied after IRB approval to collect human materials at Kyunghee University's Skin Biotechnology Center (IRB No. KHUSBC 2018-MB). All experiments were performed in accordance with relevant guidelines and regulations. The participants were divided into two age groups: 25 subjects in the younger group ("Young, " 21-36 years, mean age: 26.4 ± 3.8) and 26 subjects of the older group ("Old, " 49-67 years, mean age: 58.1 ± 4.8). Subjects were recruited regardless of gender; 25 males and 26 females. Before the sample collection, all participants washed their faces with the same facial cleanser that is a common type of foam cleanser that includes sodium lauryl sulfate (SLS) and samples were collected between 30 minutes and 1 hour after washing without any skin conditioner. Although we were aware that the skin microbiome could be affected by facial washing, we wanted to identify the skin microbial flora that are not washed out and remain constant on the skin surface even after facial washing. Skin samples were collected in five replicates per participant from two sites on the cheek and forehead, with sterile cotton swabs. To examine the overall microbial community for each facial skin area and reduce the bias between samples, we collected skin samples from five different areas of the forehead from each donor and performed the same procedure for the cheek. The swab tips were cut and transferred to collection tubes that contained 2 mL of 0.45% monopotassium phosphate (KH 2 PO 4 ), 0.6% Disodium phosphate (Na 2 HPO 4 ), and 0.05% L-cysteine HCl•H 2 O, 0.05% Tween 80, immediately stored at 4℃ before DNA extraction.
Clinical skin parameter measurement. Participants' skin condition was investigated 30 min after facial washing for skin moisture, transepidermal water loss (TEWL), sebum level, skin texture (smoothing), periorbital wrinkle (average wrinkle depth, wrinkle volume, and wrinkle area), skin redness, moles, UV spots, brown spots, porphyrin, and skin tone 5 . Skin moisture was measured by Corneometer® (Courage + Khazaka electronic GmbH, Germany). TEWL was measured by Vapometer® (Delfin Technologies, Kuopio, Finland), and the skin sebum level was measured by Sebumeter® (Courage + Khazaka electronic GmbH, Germany). The skin texture and the periorbital wrinkles were measured by PRIMOS® (Canfield Scientific, USA). The skin redness was measured by a Spectrophotometer CM700d (Konica Minolta, Japan), and the measurements indicated a * value. All spot measurements and skin tones were measured by VISIA (Canfield Scientific, USA) skin analysis. All methods were carried out according to the equipment manufacturer's instructions.
Genomic DNA extraction. The total bacterial gDNA was extracted from 510 samples collected at two different sites (cheek and forehead), conducted by a partially modified method using the QIAamp DNA Mini Kit (Qiagen, Germany, cat. no. 51306). Initially, the skin samples' collection tubes were dispensed with 20 μL Protease K, 540 μL Phosphate-buffered Saline (Corning, USA), and 60 μL AL buffer (Qiagen, Germany). After incubation at 56℃ for 3 hours, sample tubes were vortexed for 15 seconds. Then 600 μL of the AL buffer was added, vortexed for 15 seconds, and heat-treated at 56℃ for 10 minutes using a heat block. Six hundred microliters of 100% ethanol was added to the heat-treated sample, vortexed for 15 seconds, and 600 μL of the mixed sample was placed in the tube containing the spin column. After centrifugation at 14,000 rpm for one minute, the buffer that came out under the filter was discarded and the remaining sample was put back into the spin column. The above process was carried out for 510 samples, collected using the same spin column process for each facial skin area from each donor to concentrate into 102 samples. The process was carried out following the manufacturer's instructions (QIAamp DNA Mini and Blood Mini Handbook). The concentration and purity of the extracted genomic DNA were measured using a nanodrop (Thermo Scientific™, USA) according to the manufacturer's instructions. PCR amplification and 16S rRNA sequencing. The amplification of the V3-V4 region of the bacterial 16S rRNA gene was performed using barcoded 341F (5′-TCG TCG GCA GCG TC-AGA TGT GTA TAA GAG ACA G-CCT ACG GGNGGC WGC AG-3′) and 805R primer (5′-GTC TCG TGG GCT CGG-AGA TGT GTA TAA GAG ACA G-GAC TAC HVGGG TAT CTA ATC C-3′) primers 15 . PCR conditions were performed as follows: initial denaturation at 95℃ for one min, 34 cycles of 95℃ for 30 s, 55℃ for 30 s, and 72℃ for 30 s, followed by a final extension at 72℃ for 5 min. The PCR-completed sample was cleaned using a PCR Purification Kit (Qiagen, Germany), and DNA concentration was confirmed using a nanodrop. The 16S rRNA gene library construction was carried out by Illumina's Demonstrated Protocol. After a quality check using Quanti-iT pICOgREEN dsDNA assay kit (Invitrogen), Illumina Miseq was performed under 500 + 7 cycles conditions using the MiSeq Reagent Kit (Illumina, San Diego, CA, USA) according to the manufacturer's instructions.
Skin microbiome analysis. The 16S rRNA gene sequence was identified in EzBioCloud using the Microbiome Taxonomic Profiling (MTP) pipeline provided by ChunLab, Inc. 16 . Quality control of raw data was performed according to ChunLab in-house process 17 www.nature.com/scientificreports/ tool, which searches and clusters algorithms that calculate sequence similarity against reads in the EzBioCloud database (https:// www. ezbio cloud. net) 16 . The bacterial OTUs were identified by UCLUST, clustering of the 16S rRNA sequences with a ≥97% identity threshold, for taxonomic profiling analysis. Skin microbiome analysis based on the classification and identification of 16S rRNA sequences was performed using the MTP platform of the EzBioCloud database. Based on the OTUs data obtained through taxonomic profiling, differences in alphaand beta-diversity between groups were analyzed using the R program (version 4.0.4, http:// www.R-proje ct. org/) 19 . The differences in relative abundance of skin microflora between groups were also confirmed using the same program. The R program was performed in RStudio (1.4.1), an integrated development environment.

Statistical analysis.
Statistical significance was demonstrated using the Wilcoxon signed-rank test to compare the two groups; the P values of the results are indicated in each chart. A comparative analysis of the clinical skin evaluation between the younger and older groups was conducted using the Student's t-test. Non-metric multidimensional scaling (NMDS) was applied to confirm the difference in distance between groups using the Multiple Response Permutation Procedure ("mrpp") function 20 , of Vegan package (v 2.5-7) and "indval", indicator species analysis 21 , function in "labdsv" package in R program. Permutational Multivariate Analysis of Variance was used to confirm the significance of beta diversity between groups. Linear discriminant analysis (LDA) effect size tool-(LEfSe) was utilized to identify genera with relative differential abundance between groups using the web-based application "Galaxy" version 1.0 (https:// hutte nhower. sph. harva rd. edu/ galaxy/) 22 . Pearson correlation analysis was conducted using correlation test in R program to verify the correlation coefficient and significance level within P < 0.05. Consent to participate/consent to publish. Consent was obtained from all subjects to participate prior to conducting the study and posting consent was obtained from all clinical research participants. The consent for publication of articles was obtained from all authors.

Results
Skin characteristics. A total of 51 Koreans were examined for skin dermatological properties, all without skin-related diseases, none having used antibiotics or antifungal drugs within the past 3 months, and none had performed medical skincare within the past 6 months. Before sample collection, a clinical survey was conducted of each subject's skin characteristics, which were compared by dividing participants into a younger and older group (Table 1). Additionally, we divided participants into males and females to determine whether the clinical characteristics of their skin differ by gender (Table 1). Significant differences were found between the younger and older group in the moisture content (cheek), periorbital wrinkles, and the number of spots (moles, UV spots, brown spots, and porphyrin). Moreover, the elderly group showed high values in previously mentioned clinical parameters with significant differences. The clinical skin characteristics were compared by gender; there were significant differences in many skin parameters. The skin moisture content and skin tone were significantly higher in the female group; the number of UV spots and brown spots were higher in women. Simultaneously, oil content, transdermal moisture loss, average skin roughness, skin redness, amount of moles, and porphyrin levels were significantly higher in the male group.
Alpha diversity and beta diversity. The bacterial 16S rRNA gene V3-V4 region in the skin collection sample was sequenced using the Illumina Miseq platform, and the total number of reads confirmed after quality check for the raw sequence was 5,514,434. The average number of sequencing reads for 102 samples was 55,144, and the average number of reads for each group was 56,089.98 (younger group) and 54,198.7 (older group), respectively. A total of 36,344 OTUs were identified using the UCLUST tool, clustering the 16S rRNA sequences with a ≥ 97% identity threshold for taxonomic profiling analysis. The total number of OTUs in the younger group was 16,963 and the total number in the older group was 19,381. Before analyzing the difference between the two age groups, all reads were normalized to reduce the bias due to the number of samples. We compared the species richness, alpha-diversity indices such as Abundance-based Coverage Estimator (ACE) and Chao1-estimated OTU number between the two groups classified by age ( Fig. 1A and B). There was no significant difference in OTU counts between the younger and older groups (Fig. 1C). The Wilcoxon rank-sum test evaluated the two species richness indices. The ACE index was significantly higher in the older group than the younger group; (P = 0.023). In addition, the older group had significantly higher Chao1-estimated OTUs (P = 0.048 www.nature.com/scientificreports/ diversity indices, Shannon's diversity, Phylogenetic diversity, and Simpson's diversity index showed no significant differences ( Fig. 1D-F). The analysis of beta diversity was executed at the species level of phylogenetic results to confirm the difference in distance between the two age groups. The UniFrac Principal Coordinates Analysis (PCoA) showed a significant difference between the two age groups ( Fig. 2A, P = 0.009). The non-metric multidimensional scaling (NMDS) plot based on the Bray-Curtis dissimilarity also showed a significant difference. However, the delta variance was similar between the two age groups; that is, the difference within groups was not significant (Chance corrected within-group agreement A = 0.0094, P = 0.015).
In addition, we investigated the difference in alpha and beta diversities according to gender, but we could not find any significant difference ( Supplementary Fig. 1).
Differences in relative abundance. First, we focused on differences in the relative abundance at the phylum level between the two age groups. Figure 3 indicates the distribution of the four dominant phyla: Actinobacteria, Bacteroidetes, Firmicutes, and Proteobacteria. The most predominant taxonomic group at the phylum level between the two age groups was Actinobacteria and the average abundances were 62.3% and 73.4% in the younger and older groups, respectively (P = 0.025). Supplementary Figs. 2 and 3 indicate the analysis results for all skin microbial compositions. Regarding the distribution of Proteobacteria, which is also known to dominate human skin bacteria, the average abundance was significantly higher in the older group than Actinobacteria. At the family level analysis of skin samples, Fig. 4 indicates the distribution of the six most abundant for younger and older groups. The compositions of Lawsonella and Morganellaceae were significantly different according to age groups (P < 0.05). www.nature.com/scientificreports/ In addition, differences in the distribution of skin microorganisms according to gender were confirmed at the phylum level, family level, and species level. There were no significant differences in the four major skin phyla ( Supplementary Fig. 4). Figure 5 indicates the relative abundance of the four dominant families between the male and female groups. Staphylococcaceae and Corynebacteriaceae showed a high distribution in the male group; Lactobacillaceae, which belongs to the lactic acid bacteria family, showed a higher relative abundance in the female group. The difference in distribution at the family level was similar to the genus level; Staphylococcus and Corynebacterium were significantly higher in the male group and Lactobacillus was significantly higher in the female group (Fig. 6).
The linear discriminant analysis (LDA) model identifies differently abundant taxa between groups and estimates each significantly different taxon's effect size 22,23 . We evaluated the LDA effect size (LEfSe) among groups to search for a statistically significant biomarker at the genus level. The LEfSe of all genera showed 24 bacterial taxa with significant differences (Fig. 7). LEfSe analysis revealed that the genus Lawsonella was the most abundant in the younger group, and Enhydrobacter was the most abundant in the older group. Additionally, we analyzed  www.nature.com/scientificreports/   www.nature.com/scientificreports/ LEfSe by including categories for each skin site (cheek and forehead) as a subclass. It was confirmed between the two groups that Lawsonella and Enhydrobacter showed significant differences (|LDA score|> 4.0) as the results of LEfSe analysis of genera classified by age (Fig. 8A). Figure 8B indicates the relative abundances of the dominant genus in the skin site classified by age. LDA was conducted according to gender; Table 2 indicates the LEfSe of all taxonomic biomarker (taxon) without < 0.1% abundance (|LDA score|> 2.5). Likewise, our previous analysis results of the distribution of skin microbiota by gender, Corynebacteriaceae, was widely distributed in males and the LDA effect size was the largest at the family level. At the genus level, the LDA effect size of Corynebacterium and Staphylococcus was also larger in the male group. Notably, the largest LDA effect size was the distribution of the genus Xanthomonas between the male and female groups. LDA analysis confirmed that Xanthomonas citri were present at a higher proportion in the female group and Corynebacterium tuberculostearicum strains were present in a higher proportion in the male group. Moreover, Enterococcus faecalis, which is a type of lactic acid bacteria, was widely distributed in women and the LDA effect size of Lactobacillus helveticus was 2.34, which had a higher distribution in the female group (Table 2 and Supplementary Fig. 5).

Correlation between skin symbiotic bacteria and clinical skin parameters. Based on the results
of LEfSe analysis, we investigated which skin clinical parameters correlated with the relative abundance of skin microflora. We analyzed Pearson correlation, which was classified by skin site (cheek and forehead) between the bacterial taxa at a species level and the subjects' clinical skin parameters. The cheek skin parameters, moisture content, and amount of brown spots were negatively correlated with the Lawsonella clevelandensis group belonging to the genus Lawsonella of subjects with statistical significance. (Fig. 9). The skin parameter-related correlation of Lawsonella clevelandensis in the younger group showed that the gradient between the cheek moisture content and abundance of Lawsonella clevelandensis was more negatively inclined. Moreover, regarding the gradient between the amount of brown spots and the abundance of Lawsonella clevelandensis, both groups showed similar patterns (Figs. 9B and 9D). The Enhydrobacter aerosaccus group belonging to the genus Enhydro- www.nature.com/scientificreports/ bacter was positively correlated with the cheek moisture content with statistical significance. (Fig. 10). However, there was no significant correlation with the forehead skin parameters of either of the two species. Figure 11 shows the skin parameter-related correlation with the relative abundance of Cloacibacterium haliotis. There was a high correlation coefficient between the abundance of Cloacibacterium haliotis and the periorbital wrinkle area (R = 0.48). We analyzed the correlation between gender and clinical skin parameters for the genus that showed significant differences. The genera Staphylococcus and Corynebacterium, which were highly dominant in male skin, were examined for their correlation with skin characteristics on subjects' cheeks and foreheads. Staphylococcus was negatively correlated with the number of UV spots in cheeks with statistical significance (R = − 0.42, P = 0.0023). These results suggest that the reduction of the UV spots could be predicted by Staphylococcus abundance (Fig. 12). Corynebacterium was positively correlated with the TEWL in cheeks with statistical significance (Fig. 13). Moreover, we focused on the correlation between skin characteristics and the relative abundance of species of Staphylococcus aureus, which is known to be closely related to atopic dermatitis 24 . It was confirmed that the higher the distribution of S. aureus, the lower the skin moisture content in the cheeks and the higher the transdermal moisture loss (Fig. 14).

Discussion
This study aimed to profile the distribution of facial skin microbiota according to age and gender in healthy Koreans and to confirm their association with skin characteristics. There was a significant difference in periorbital wrinkles and number of facial spots when comparing the clinical skin parameters between the younger and older age groups. These results are similar to the previous reports that reveal that chronological increase in age was related to the clinical appearance of facial wrinkles and facial hyperpigmentation 25,26 . However, the cheek moisture content was significantly higher among the elderly. This was a contrasting result with some previous reports that skin hydration decreases with age 27 . However, from other observations, there was no association between skin hydration and age 26,28 . There was no significant difference between the age groups in terms of forehead skin hydration in this study. The differences in skin characteristics according to gender were found to be more significant than the differences by age without periorbital wrinkles. In particular, the male group had higher sebum contents than the female group, which could be associated with the relative abundances of Staphylococcus and Corynebacterium that use facial sebum as nutrients 29,30 . The overall skin characteristics results confirm that the female skin had a more positive result from a cosmetic perspective but the number of UV spots was found to be higher in the female group.
In the results of skin microbiome analysis, the alpha-diversity index was significantly higher in the elderly, which was similar to previous reports that the older group showed a tendency toward a higher alpha diversity than the younger group for all skin microbiomes in Japanese women 5 . Nevertheless, there was no difference in alpha diversity by gender (data not shown). A study investigating the interactions between the host and the skin www.nature.com/scientificreports/ microbiome reported that the contribution of gender to skin microbial diversity likely arises as a downstream effect of male and female steroid production 31 . Since we had no gender restrictions, in contrast to the results of age-related beta diversity of Japanese women 5 and Chinese women 10 , the distance differences between the two age groups were significant, but it was difficult to find a largely distinguished difference. The relative abundance of skin microbiota between the younger and older groups confirmed that Actinobacteria and Proteobacteria, the two known major human skin phyla, are dominant in the younger age group and older age group, respectively. Indeed, the genus of Lawsonella was more predominant and Lawsonella clevelandensis was identified as a major species in the younger group. Recent studies showed that L. clevelandensis was among the most common species on the human skin, scalp, and nostrils 32-36 . Additionally, there was a significant negative Table 2. Skin characteristics of the subjects by gender. Values in the table were presented through linear discriminant analysis (LDA) effect size (LEfSe). FDR false discovery rate; Female, an average of relative abundance in the female group; Male, an average of relative abundance in the male group. www.nature.com/scientificreports/     www.nature.com/scientificreports/ correlation between L. clevelandensis and the cheek moisture contents and brown spots on the face, although the correlation was not high. Enhydrobacter aerosaccus, which has a high relative abundance in the elderly's function, has been discovered through a recently published skin microbiome study in China 37 . In that study, E. aerosaccus and M. osloensis were taxonomically considered to be the same species. Furthermore, the M. osloensis group is dominant in less-hydrated skin, which contrasts with the results of this study. These results can be explained by differences in the participant selection criteria, their living environment and ethnicity. We confirmed that the genus Staphylococcus was significantly associated with UV spot number on cheeks and Corynebacteria were significantly associated with TEWL. The genera Staphylococcus and Corynebacterium dominated in males 38 , known as major normal skin bacteria, and are reportedly related to the sebum or hydration levels of the facial skin 39 . Subsequently, the distribution of Staphylococcus aureus, which generally colonizes human skin and mucosa 40 , was predominant in the male group. We found that an increase in Staphylococcus aureus was significantly associated with an increase of TEWL and decreased skin moisture levels in subjects. With Staphylococcus aureus, TEWL was reportedly higher and hydration level was lower in atopic dermatitis (AD) patients 41 , and low hydration was associated with high S aureus growth in AD patients 24 . Characteristically, in this study, the distribution of Lactobacillus helveticus was higher in women, which is similar to previous studies 38 . Several recent studies have investigated the beneficial effects of Lactobacillus helveticus fermentation within human skin epidermal cells 42,43 . Moreover, there have been cases in which Lactobacillus helveticus was orally administered as a probiotic to atopic patients to improve symptoms 36 , but its role as a human skin symbiotic bacteria has not yet been confirmed.
In conclusion, this study attempted to identify the microbiome on human facial skin in Korea through 16S rRNA analysis using next-generation sequencing analysis technology. We proved the age-related distribution of facial skin microbiome and the difference in the distribution of facial skin microbiome according to gender. We found out how these differences in facial skin microbiota distribution correlate with the skin's clinical characteristics. Correlation with clinical indicators in the skin microbiome at the species level was derived, and this approach is necessary to determine the balance of skin microbiome that we need to control for skin health. Nowadays, as various anti-aging cosmetics are on the market and various skincare industries develop, there can be a difference between actual age and biological skin age. It will be better to infer skin-improving microbials in consideration of the skin characteristics classified by sex hormones or men and women's body characteristics 14 . Future studies must examine the facial skin microbiome and clinical skin parameter changes by considering the biological skin age and specifying cosmetics used in daily life.

Data availability
The subject's information on this study and the results of the survey on subjects' skin type and skin clinical characteristics are available as Supplementary Data 1 and 2, respectively. OTUs results and taxonomic composition www.nature.com/scientificreports/ results for each skin sample derived from 16S rRNA gene V3-V4 region sequencing are summarized in Supplementary Data. The dataset for the older and younger age groups are categorized for each skin site (forehead and cheek) and presented in Supplementary Data 3-6. The 16S rRNA gene sequencing data for all human-derived skin samples were registered as PRJNA723064 in the NCBI SRA (Sequence Read Archive) database. www.nature.com/scientificreports/