Skin microbiome-biophysical association: a first integrative approach to classifying Korean skin types and aging groups

Abstract

Introduction: The field of human microbiome research is rapidly expanding beyond the gut and into the facial skin care industry. However, there is still no established criterion to define the objective relationship between the microbiome and clinical trials for developing personalized skin solutions that consider individual diversity.

Objectives: In this study, we conducted an integrated analysis of skin measurements, clinical Baumann skin type indicator (BSTI) surveys, and the skin microbiome of 950 Korean subjects to examine the ideal skin microbiome-biophysical associations.

Methods: By utilizing four skin biophysical parameters, we identified four distinct Korean Skin Cutotypes (KSCs) and categorized the subjects into three aging groups: the Young (under 34 years old), the Aging I group (35-50), and the Old group (over 51). To unravel the intricate connection between the skin's microbiome and KSC types, we conducted DivCom clustering analysis.

Results: This endeavor successfully classified 726 out of 740 female skin microbiomes into three subclusters: DC1-sub1, DC1-sub2, and DC2 with 15 core genera. To further amplify our findings, we harnessed the potent capabilities of the CatBoost boosting algorithm and achieved a reliable framework for predicting skin types based on microbial composition with an impressive average accuracy of 0.96 AUC value. Our study conclusively demonstrated that these 15 core genera could serve as objective indicators, differentiating the microbial composition among the aging groups.

Conclusion: In conclusion, this study sheds light on the complex relationship between the skin microbiome and biophysical properties, and the findings provide a promising approach to advance the field of skincare, cosmetics, and broader microbial research.

Keywords: Korean skin cutotypes; core genera; dermatologic conditions; microbiome; skin microbiome-biophysical association.

Figure 1

Screening Korean Skin Cutotypes (KSC): Research strategies and BSTI analysis unveiled. (a) This diagram represents the sequence of steps in determining KSC type for a total of 756 female subjects. Based on 39 skin clinical measurement data parameters, four skin criteria were finally selected and four KSC types were determined. (b) Distribution of female subjects based on sebum production, skin sensitivity, pigmentation, and wrinkle severity categories. The x- and y-axis represent the calculated score in each category and the accumulated number of subjects, respectively. The red dashed line is the median. (c) Determined proportions of 16 Bauman types in 752 female subjects. Detailed BSTI results for each subject are presented in Supplementary Table 1 .

Figure 2

Assessment of key metrics selection and predictive performance using biophysical parameters. (a) Correlation analysis through the Pearson coefficient method between each of the 39 biophysics parameters. As shown in the scale bar on the right, blue and red colors represent positive and negative correlations, respectively, with higher correlations being darker. (b) Polygonal chart of age-related changes in six metrics representing skin condition (elasticity, skin tone, oil, pore, wrinkle, moisture). (c) Mean receiver operating characteristic (ROC) curves for age estimation using 39 biophysical parameters and figure (d) the representative six metrics. ROC curves were used to compare the predictive power of the six metrics. The Area Under the ROC Curve (AUC) scores for each age group are noted in the legend and represented by different colors. The green dashed line indicates the expected AUC for a random chance classifier. Plots depict the tradeoff between true-positive and false-positive. The closer the curve follows the left-hand border and the top border of the ROC space, the more accurate the test.

Figure 3

Determination and compartmentalization of the skin type for 705 Korean female. (a) Scatter plots of tone/elasticity (top) and oil/moisture (bottom) for the tertile of groups by each criterion. Each individual in the high, middle, and low values of the x- and y-axis criterion was plotted in green, blue, and orange, respectively. (b) Principal component analysis based on tone/elasticity and oil/moisture criterion values. The integrated coordinates of the skin criteria show a characteristic distribution across the four types based on the skin type determination strategy. Individuals in the middle range were marked in gray. For example, high values for both tone/elasticity and oil/moisture were annotated as HH and low values as LL. The mountain diagram indicates the distribution of number of the enriched individuals. (c) Age-dependent clockwise change in skin type. The green, orange, blue, and red contours represent the types HH, HL, LH, and LL, respectively. The triangles represent the number of Korean female subjects in each age group, demonstrating that the skin type determination criteria in this study are efficiently representative of skin clinical variation.

Figure 4

Turning points in skin condition and changes in skin type in aging groups. (a) The line chart represents the distribution of high and low skin type by age for 4 skin type criteria. The green and orange lines denote the proportion of subjects with high (positive) and low (negative) values in the age group. The red dotted lines indicate inflection points where tone, elasticity, and oiliness change. (b) Determination and the ratio of three Aging groups in skin turning points, based on four skin criteria. (c) Polygon chart for four skin types and biophysical variations by the Aging group. The green, orange, blue, and red polygons denote the types HH, HL, LH, and LL, respectively. The grey one is for the mean of each the Aging group.

Figure 5

The skin-aging trends for four skin criteria dimensions. (a) Aging associations of tone and elasticity and figure (b) personal skin condition variations with oil and moisture. The binned graph shows the distribution of tone/elasticity or oil/moisture values and the number of samples per aging group. The green and red dashed lines indicate the mean and median values, respectively. This figure describes that among the key skin criteria, Tone/elasticity shows a sharp difference as getting old compared to oil/moisture. (c, d) Scatter plots representing relationships between tone/elasticity and between oil/moisture. In the three determined aging groups, tone/elasticity is converging to sharply lower values, while oil/moisture indicates an overall decrease in the oil criteria but similar values for the moisture criteria. Based on skin type with tone/elasticity and oil/moisture, the high group is colored green, and the low group is colored blue. The mountain diagram indicates the distribution of number of the enriched individuals.

Figure 6

Comprehensive microbiome comparative analysis according to the decided Aging groups. (a) Box plot showing microbial alpha-diversity comparisons between three different age groups (Young, Aging 1, and Old groups). Shannon and Simpson alpha-diversity indices were applied to this alpha-diversity estimation, and it is a measure for confirming the microbial richness and evenness in each group (‘p’ means the p-value). (b) Relative abundance bar plot showing the difference of relative bacterial frequency within each group for 10 different core-bacterial genera (Cutibacterium, Streptococcus, Staphylococcus, Rothia, Corynebacterium, Neisseria, Actinocymes, Haemophilus, Fusobacterium, and Veillonella). Color notation information for each genus is indicated in the footnote on the right side of the figure. (c) LEfSe (Linear discriminant analysis Effect Size) analysis result confirming the distinct bacterial taxonomy showing statistically significant differences in relative frequency between each group. The threshold on the logarithmic LDA (Linear Discriminant Analysis) score for discriminative features was set to 2.0 (indicating significant differential abundance).

Figure 7

Determining KSC microbiome clusters using diversified clustering approaches. (a) Determine the optimal cluster with four clustering methods. The x-axis represents the number of k-clusters determined at each index, and the y-axis represents the score at the index. (b) Principal coordinate analysis (PCoA) indicates sample similarity based on genus-level taxonomic composition. A total of 726 subjects significantly grouped by two clusters at first. Applying the DivCom algorithm uncovered one hidden sub-cluster. The three determined clusters are labeled red (DC1-sub1), orange (DC1-sub2), and blue (DC2), respectively. (c) Midpoint-rooted maximum-likelihood phylogenetic tree based on taxonomy abundance of each subject. The tree contains 247 DC-sub1, 268 DC-sub2, and 211 DC2, representing genera of bacteria that are significantly enriched in each cluster. Figure (d, e) show the association between the two main clusters and three subclusters, respectively, and the clinical biophysics used to determine the KSC. The asterisk (*) represents the p-value of the statistical test (* < 0.05, ** < 0.01, *** < 0.001, and **** < 0.0001).

Figure 8

Heatmap of the Heatmap of the pairwise comparison analysis for the core genera. (a–c) The 10 of 15 core genera, the heatmap plots of negative log10-transformed p-values from all possible pairwise comparisons using ANCOM-BC2 were generated. Relative P-value differences between the two groups are marked with + (red) and - (blue). They all have different composition patterns, and the scale was set by their significance. The figure (a) shows a significant composition difference between the Young/Aging I group and the Old group, and the figure (b) shows that the LH and LL types in the Aging I group have the same composition pattern as the Old group. The figure (c) shows the specificity of certain KSCs regardless of the aging group, for example, HL in the younger group and LH in the older group with a higher composition than the other groups. (d) Predicting decision performance of KSCs by aging group using 15 core genera. The scale represents the z-score transformed values of the AUC score obtained from the CatBoost algorithm. The red and green colors indicate high and low feature importance, respectively.

Figure 9

Relative abundances of the functional pathways predicted by PICRUSt2. The legend at the top provides the aging group of Korean women, their age, KSC based on Tone/elasticity, oilness/moisture, and their scores gradually. Scope and description are in the left legend. The heatmap colors with + (red) and - (blue) indicates the z-score normalized relative abundances. The young and aging groups showed contrasting functional enrichment results with the old group, and even between the young and aging groups, the young group has a distinctive enrichment pattern in KSC HH and HL types. The 59 major functional prediction of KEGG pathways are listed in detail in the Supplementary Table 8 .