Novel Genetic Associations for Skin Aging Phenotypes and Validation of Previously Reported Skin GWAS Results

Featured Application

We confirmed six genetic skin indexes, which will be applied to develop a novel Korean-specific Skin DTC genetic test and custom cosmetics.

Abstract

Facial skin characteristics are complex traits determined by genetic and environmental factors. Because genetic factors continuously influenced facial skin characteristics, identifying associations between genetic variants [single-nucleotide polymorphisms (SNPs)] and facial skin characteristics may clarify genetic contributions. We previously reported a genome-wide association study (GWAS) for five skin phenotypes (wrinkles, pigmentation, moisture content, oil content, and sensitivity) conducted in 1079 subjects. In this study, face measurements and genomic data were generated for 261 samples, and significant SNPs described in previous papers were verified. We conducted a GWAS to identify additional genetic markers using the combined population of the previous study and current study samples. We identified 6 novel significant loci and 21 suggestive loci in the combined study with p-values < 5.0 × 10⁻⁸ (wrinkles: 4 SNPs; moisture content: 148 SNPs; pigmentation: 6 SNPs; sensitivity: 18 SNPs). Identifying SNPs using molecular genetic functional analysis is considered necessary for studying the mechanisms through which these genes affect the skin. We confirmed that of 23 previously identified SNPs, none were replicated. SNPs that could not be verified in a combined study may have been accidentally identified in an existing GWAS, or the samples added to this study may not have been a sufficient sample number to confirm those SNPs. The results of this study require validation in other independent population groups or larger samples. Although this study requires further research, it has the potential to contribute to the development of cosmetic-related genetic research in the future.

1. Introduction

To investigate the variability in skin characteristics among individuals, genetic studies of skin aging are currently being conducted. Several GWASs are in progress to assess skin phenotypes, which have identified gene variants associated with pigmentation in OCA2, MC1R, ASIP, SLC45A2, TYR, IRF4, and BNC2 that are linked to hyperpigmentation spots and pigmentation characteristics such as skin color [1,2,3,4,5,6,7,8]. The first variants recognized as affecting pigmentation in animals were MC1R variants, which are also associated with human pigment phenotypes [9], including freckles [6,10] and sun sensitivity [11,12].

In addition, this phenotype contributes to the risk of melanoma, Basal cell carcinoma and Squamous cell carcinoma, and the genes contribute to skin disease risks, in addition to skin phenotypes [13,14,15,16,17,18,19]. IRF4, MC1, and SLC45A2 were identified as significant loci associated with sensitive skin [12], and these genes are also associated with pigmentation [20], melanogenesis [21,22], and skin-related cancer [23]. Studies have shown that the EGFG, MMP16 and COL17A1 genes are associated with an increased risk of developing wrinkles, in addition to skin pore size [24]. An SNP located proximal to TGIF1 is related to sagging eyelids [25], and the STXBP5L locus is related to facial photoaging [26]. An SNP located upstream of the Patched-1 gene was identified as being related to stress-induced skin conditions and aging [27], and Myc interacts with skin hydration [28].

Skin undergoes decreases in moisture, elasticity, and luster due to aging and external factors, resulting in increased wrinkle depth and the hyperpigmentation of blemishes [29]. Previous studies have confirmed similar trends in age-specific skin phenotype measurements. In a previous study, age increased the measured value of wrinkles and slightly increased pigmentation. However, the measured values of moisture and oil contents decreased with age [30].

Environmental changes and personal life patterns are considered external factors for skin aging [31,32,33], and a genetic approach to skin aging is necessary to better understand individual skin characteristics. A recent study reported new genomic loci [34,35] related to pigmented facial spots in Koreans and identified SNPs associated with melanin, gloss, hydration, wrinkles, and elasticity [29]. However, few studies have investigated the genetic basis of skin phenotypes among Koreans.

This research team previously performed GWASs of skin characteristics among 1079 Korean subjects and identified 23 genes related to skin characteristics [30]. Five cosmetic phenotypes (wrinkles, moisture content, pigmentation, oil content, and sensitivity) were evaluated. This GWAS identified SNPs in 23 genes for the five skin phenotypes, including two loci associated with both wrinkles and pigmentation (FCRL5 and OCA2), four wrinkle-related genes (FCRL5, REEP3, ADSS, and SPTLC1), five moisture-related genes (TBS4, TRPM3, CEMIP2, CTSH, and TTC39C), six pigmentation-related genes (OCA2, NCLN, TNS1, CDC42BPA, HS3ST4, and UNCX), five oil-related genes (SYN2, CNDP1, GAS6, INSR, and TNFRSF19), and the sensitivity-related gene CREB5. This present study aimed to determine whether these of 23 previously reported genes could be validated as significant genetic contributors in the new sample of 261 participants.

2. Materials & Methods

2.1. Participants

The sample population investigated in this study comprised 261 Korean women recruited between July 2021 and January 2022 at P&K Skin Research Center (Seoul, Korea). All recruited participants were women without skin-related diseases, and their average age was 45.10 years (

Table 1. Baseline characteristic and skin measurement information of the subjects.

			Phenotype Measurements				Target Phenotypes
			Reported Study	Present Study		Combined Study	Coded Phenotype
							Reported Study	Present Study	Combined Study
Characteristics			Average ± SD or N (%)
Population (n)			1079	261		1340
Gender (female, %)			100%
Age (year, mean ± SD)			40.81 ± 10.90	45.10 ± 11.92		41.65 ± 11.23
Phenotype (equipment)	Measure item	Code			p-value		Average ± SD or N (%)
Wrinkle (Primos CR)	Eye: Average roughness (Ra)	W101	20.81 ± 5.24	21.74 ± 5.06	<0.001	20.64 ± 5.22	7.75 ± 2.58	8.95 ± 2.23	7.74 ± 2.51
	Eye: Maximum depth (Rmax)	W102	198.47 ± 71.97	221.68 ± 60.84	<0.001	198.34 ± 69.75
	Glabela: Average roughness (Ra)	W103	25.18 ± 7.06	20.25 ± 5.96	<0.001	23.75 ± 6.92
	Glabela: Maximum depth (Rmax)	W104	193.75 ± 69.61	163.27 ± 50.32	<0.001	184.7 ± 62.81
Moisture (Corneometer)	Glabela: Moisture content (A.U.)	A101	64.84 ± 9.97	57.57 ± 11.75	<0.001	63.59 ± 10.55	3.98 ± 1.38	5.45 ± 0.72	4.24 ± 1.39
	Cheek: Moisture content (A.U.)	A102	70.13 ± 9.48	60.13 ± 11.16	<0.001	68.18 ± 10.34
Pigmentation (Mexameter) (CM-700d)	Pigmentation site: Melanin (M.I.)	M101	166.15 ± 35.83	135.38 ± 32.26	<0.001	160.56 ± 36.89	8 ± 2.35	6.39 ± 1.85	7.77 ± 2.31
	Nonpigmentation site: Melanin (M.I.)	M102	119.25 ± 28.35	98.49 ± 24.85	<0.001	116.51 ± 28.94
	Pigmentation site: Brightness (L*)	R201	59.95 ± 2.62	60.80 ± 2.91	<0.001	60.25 ± 2.68
	Nonpigmentation site: Brightness (L*)	R202	63.22 ± 2.33	64.27 ± 2.45	<0.001	63.49 ± 2.37
Oil (Sebumeter)	Glabela: Oil content (μg/cm2)	L101	85.19 ± 61.23	45.56 ± 30.69	<0.001	78.91 ± 59.68	4.07 ± 1.48	2.95 ± 1.3	3.94 ± 1.47
	Cheek: Oil content (μg/cm2)	L102	48.57 ± 42.49	38.46 ± 24.90	<0.001	47.55 ± 40.68
Sensitivity (10% lactic acid)	Sensitive/Nonsensitive	S101	534 (49.5%)/ 545 (50.5%)	206 (78.9%)/ 55 (21.1%)	<0.001	740 (55.22%) /600 (44.78%)

). All participants provided written informed consent, and this study was approved by the Institutional Review Board of the Theragen Bio Institute (IRB No.: 20200052001-210616-GP-006-001).

2.2. Measurement of Skin Phenotypes

We used the same method as in previous reports for the skin measurements [30]. To give a more detailed description, skin traits were measured using various measuring devices. A Primos CR (Canfield Scientific, Parsippany, NJ, USA) was used to obtain wrinkle area measurements, including the roughness and wrinkle depth of the skin on the side of the eye and glabella (Figure 1). A CM-825 Corneometer^® (EnviroDerm Services Ltd., Hedworth, Grange Court, Westbury-on-Severn, Gloucestershire, UK) was used to obtain moisture content measurements for the skin of the glabella and cheek. A Mexameter^® MX 18 (Courage + Khazaka Electronic GmbH, Köln, Germany) was used to obtain melanin measurements, and a CM-700d (Konica Minolta Inc., Tokyo, Japan) instrument was used to measure brightness. Oil content was measured on the skin of the glabella and right cheek using an SM 815 Sebumeter^® (Courage + Khazaka Electronic GmbH., Köln, Germany). Finally, a sensitivity analysis was performed by evaluating the sensitivity of skin treated with a lactic acid sting test. The lactic acid sting test (LAST) is a classical method to identify sensitive skin. The subjective assessment for sensitivity relies on the results of the subject’s skin response to stimuli, such as the LAST. The detailed method is described as follows: Subjects were rested for 1 h in a waiting room under constant temperature and humidity conditions (indoor temperature 20–25, humidity 40–60%) after face washing for equivalence of skin condition. The lactic acid sting test was performed immediately after the hydration process in which the subject applies a steam towel to the face for 5 min in a comfortable lying position. The ten percent lactic acid was dropped around the nasolabial fold of subject, and then rubbed gently with a cotton swab. After 1, 5, and 10 min of application, the subject’s subjective sense of self-inflicted feeling was evaluated with a score of 0 to 3 points. (0 = None, 1 = Mild, 2 = Moderate, 3 = Severe). We did not use a device to obtain sensitivity measurements because skin sensitivity differs between individuals; instead, we recorded the skin’s reaction.

2.3. Target Skin Phenotype Grading Scale

Skin phenotypes for the dependent variable in GWAS were defined using the same method as in a previous report [30]. To describe the method in greater detail, before performing GWASs, each measured value was divided into tertile groups using codes because the measured values of the same phenotype may differ depending on the measuring device. The classification of each phenotype is shown in Table S1, and the integrated score for each phenotype is shown in Table 1. The scoring process is shown as a flow diagram in Figure S1. Figure S1 is a schematic diagram showing that the integrated score is the sum of all four skin characteristic scales. This approach aims to demonstrate the generalizability of results based on tertiles, even when using different measurement devices.

2.4. Identification of SNP Genotype Based on SNP Array

Samples were obtained by oral swab. DNA was extracted using the Exgene™ Tissue SV (GeneAll, Seoul, Korea). To obtain the genetic variant information for 820,000 SNPs in the whole human genome, we used a partially customized microarray chip called a TPMRA array. The TPMRA array was partially customized from a commercialized microarray chip (Asian PMRA chip). The customized chip was produced by Thermo Fisher (Thermo Fisher Scientific, Waltham, MA, USA) for the design and production. The Theragen PMRA array is an SNP array that combines DNA hybridization, fluorescence microscopy, and solid surface DNA capture. The SNP array comprises three major components—an allele-specific oligonucleotide probe, fragmented nucleic acid sequences labeled with fluorescent dyes, and a detection system for records—and interprets the hybridization signal.

We determined the genotypes of 820,000 SNPs for each participant sample using the TPMRA. The amplified DNA was randomly digested into 25- to 125-bp fragments, purified, and resuspended. Fluorescently labeled allele-specific oligonucleotide probes were used to bind the human genomic DNA in the Theragen PMRA. To prepare for detection, the probe-bound DNA was hybridized. After hybridization, unbound DNA fragments were washed away to minimize the noise of nonspecific ligation.

The experimental results were initially analyzed by Affymetrix Power Tools to determine the quality of the genotyping results based on the manufacturer’s instructions. Briefly, the ratio of the target probe signals and background noise signals (DQC: Dish Quality Control) were higher than 0.82. We also applied the initial SNP filtering of FLD (Fisher’s linear discriminant) > 3.6, HetSO (heterozygous strength offset) > −0.1, OTV (off-target variant) > −0.3, Finally, we selected 778,329 recommended SNPs for further analyses.

For each of the 778,329 SNPs, quality control procedures were performed before association testing. SNP sets were filtered according to the genotype call rates (≥0.95) and MAF (≥0.01). The Hardy–Weinberg equilibrium (HWE) was calculated for individual SNPs using an exact test. All SNPs reported for this study this article exhibited HWE p-values > 0.01. After filtering, 313,740 SNPs were analyzed on chromosomes 1 through 22.

2.5. Imputation of SNPs

While previous reports only used the experimental genotyping data [30]. In this study, to expend the genome coverage of our genome-wide SNP dataset, we performed imputation analysis using minimac4 [36] at the Michigan imputation server with the 1000 G Phase 3 v4 reference panel [37]. Phased GWAS genotypes were uploaded to obtain imputed genomes. After imputation, 5,512,463 SNPs were analyzed on chromosomes 1 through 22.

2.6. Genome-Wide Association Scan

The linear regression analyses were performed for each skin phenotype with controlling the age as the covariate by using PLINK version 1.90 (https://www.cog-genomics.org/plink/ accessed on 6 September 2022) [38]. Also, the phenotype characteristics were analyzed by the SPSS (IBM SPSS Statistics Inc., New York, NY, USA) [39] and R Statistical Software [40]. We calculated the beta coefficient and standard error (SE) values for the association test. We used the inverse-variance meta-analysis method, assuming fixed effects, with Cochran’s Q test to assess differences in heterogeneity between previous and present populations [41].

3. Results

3.1. Study Population and the Results of Skin Measurements

A summary of facial phenotype measurement data is shown in Table 1. The collection conditions and collector (P&K agency) were the same as those described in a previous study. The sample recruitment method was the same as the previous study, and the sample was collected randomly. However, in the randomly recruited sample this time, the measurements of several skin phenotypes showed different characteristics than before. However, since both populations were randomly collected, they were combined and analyzed. We analyzed the skin phenotypes of 1340 total participants (100% women) with an average age of 41.65 years. The study population comprised 1079 participants from the previous study, with a mean age of 40.81 years, and 261 participants from the present study, with a mean age of 45.10 years. The differences between cohorts were tested statistically and are listed in Table 1. The five facial phenotypes were analyzed using each measurement device. Additionally, phenotypic measurement values were classified by scoring, and the measured scores were integrated. Table 1 shows the scores of the corresponding tertile groups for each of the four scales summed up and with one mean integrated score. The measured image difference between the different values for wrinkles are shown in Figure 2.

3.2. Genome-Wide Association Results

We identified 6 novel significant loci and 21 suggestive loci in the combined GWAS (

Table 2. Summary of significant new skin phenotype loci according to GWAS.

										Combined (n = 1340)
Trait	SNP	CHR:BP	REF	ALT	Feature	Real by Gene	MAF			Minor Allele	BETA ± SE a	p
							This Study	KOR	ASN
wrinkle	rs72721056	1:93620281	C	A	IntronVariant	BCAR3	0.066	7 × 10−4	0	A	0.75 ± 0.15	8.90 × 10−7
	rs34682406	7:6376481	G	C	IntronVariant	RAC1	0.054	0.038	0.03	C	0.88 ± 0.17	1.30 × 10−7
	rs10739924	9:92072959	A	C	IntronVariant	SPTLC1	0.419	0.433	0.07	C	0.37 ± 0.08	2.62 × 10−6
	12:29159521	12:29159521	G	A	Upstream	FAR2	0.421	1	0	GA	−0.39 ± 0.08	2.01 × 10−6
	rs61970162	13:108997793	G	A	Intron	MYO16	0.083	0.01	0.01	A	0.69 ± 0.14	8.44 × 10−7
	rs7159152	14:43993202	G	T	Intergenic	-	0.32	0.319	0.34	T	0.38 ± 0.08	3.70 × 10−6
Moisture	rs149963203	2:202264235	A	G	Upstream Variant	NOP58	0.189	0.064	0.17	G	−0.42 ± 0.07	6.14 × 10−10
	rs201058	6:6725784	A	C	Upstream Variant	LY86	0.435	0.185	1	A	−0.35 ± 0.05	2.38 × 10−11
	rs16912205	9:107096948	C	T	Downstream	ZNF462	0.191	0.108	0.16	T	0.46 ± 0.07	8.25 × 10−12
	rs56133064	12:124263243	A	G	Intergenict	RFLNA	0.07	0.031	0	G	0.62 ± 0.1	8.57 × 10−10
	rs11640236	16:85376910	C	A	Intron Variant	GSE1	0.175	0.089	0	A	0.65 ± 0.07	1.40 × 10−21
Pigmentation	rs4233226	1:227330340	C	A	Downstream	CDC42BPA	0.473	0.49	0.52	C	0.42 ± 0.09	1.00 × 10−6
	rs1830202	2:78054305	A	G	Downstream	LRRTM4	0.463	0.476	0.44	G	0.39 ± 0.09	9.08 × 10−6
	rs56133064	12:124263243	A	G	Intergenic	RFLNA	0.07	0.031	0	G	−0.89 ± 0.17	1.55 × 10−7
	13:32311074	13:32311074	-	-	-	-	0.413	-	-	GT	−0.46 ± 0.09	1.39 × 10−7
	rs11640236	16:85376910	C	A	Intron Variant	GSE1	0.175	0.089	0	A	−0.62 ± 0.11	5.66 × 10−8
Oil	rs1918689	2:84679729	G	A	IntronVariant	DNAH6	0.09	0.068	0.08	G	0.43 ± 0.09	2.89 × 10−6
	rs113608863	6:492601	A	G	IntronVariant	EXOC2	0.054	0.059	0.07	G	−0.52 ± 0.12	9.44 × 10−6
	rs62617088	12:52487848	C	A	MissenseVariant	KRT6A	0.054	0.056	0	A	−0.54 ± 0.12	4.16 × 10−6
	rs2639681	15:53030749	T	C	Downstream	ONECUT1	0.254	0.267	0.19	C	0.27 ± 0.06	7.97 × 10−6
	rs7230869	18:49630634	T	A	Downstream	LIPG	0.388	0.377	0.41	T	0.25 ± 0.06	6.19 × 10−6
	rs8108544	19:45916195	C	A	Downstream	NANOS2	0.177	0.188	0.19	A	−0.32 ± 0.07	4.53 × 10−6
Sensitivity	SNP	CHR:BP	REF	ALT	feature	Real by gene	This study	KOR	ASN	A1	OR (95%CI) b	p
	rs55999874	7:12575508	G	T	Intron Variant	SCIN	0.308	0.212	0.37	T	1.57 (1.34–1.85)	5.01 × 10−8
	rs28526775	8:82317642	A	C	Downstream	SNX16	0.127	0.125	0.01	C	0.53 (0.42–0.67)	1.49 × 10−7
	rs41308918	9:107097334	C	T	IntronVariant	ZNF462	0.19	0.108	0.16	T	1.67 (1.36–2.05)	9.96 × 10−7
	rs1106351	17:55247276	G	A	Downstream	MIR548Q	0.187	0.119	0.13	A	1.66 (1.35–2.04)	1.21 × 10−6
	rs78295829	22:31575365	C	T	Missense Variant	SFI1	0.092	-	0	T	2.96 (2.15–4.07)	2.44 × 10−11

Chr: chromosome; BP: base pair; Pos: position; KOR: Korean; ASN: Asian; SE: standard error. The reference MAF was obtained from the KRGDB (http://coda.nih.go.kr/coda/KRGDB/index.jsp (accessed on 6 September 2022)) and Ensembl DB (https://asia.ensembl.org (accessed on 6 September 2022)). The important GWAS loci were mapped to genes based on physical proximity. ^a Calculated β and SE by linear regression. ^b We analyzed OR using logistics regression. The SNPs with the top signals based on GWAS satisfying the cutoff value of p < 5.0 × 10⁻⁸ are underlined.

). In the genome-wide array, 176 significant SNPs were identified, with p-values < 5.0 × 10⁻⁸ (wrinkles: 4 SNPs; moisture content: 148 SNPs; pigmentation: 6 SNPs; sensitivity: 18 SNPs). Among these, we identified six novel, significantly clustered loci with p-values < 5.0 × 10⁻⁸. For the moisture trait, we found 5 significant loci, rs149963203 (β ± SE = −0.42 ± 0.07; p = 6.14 × 10⁻¹⁰), rs201058 (β ± SE = −0.35 ± 0.05; p = 2.38 × 10⁻¹¹), rs16912205 (β ± SE = 0.46 ± 0.07; p = 8.25 × 10⁻¹²), rs56133064 (β ± SE = 0.62 ± 0.1; p = 8.57 × 10⁻¹⁰) and rs11640236 (β ± SE = 0.65 ± 0.07; p = 1.40 × 10⁻²¹). We found 1 significant locus for the sensitivity trait, rs78295829 (odds ratio: 2.96; 95% confidence interval [CI]: 2.15–4.07). Two SNPs, rs11640236 (β ± SE = −0.62 ± 0.11; p = 5.66 × 10⁻⁸) and rs55999874 (odds ratio: 1.57; 95% CI: 1.34–1.85), showed p-values approaching <5.0 × 10⁻⁸.

The combined GWAS identified 575 suggestive SNPs with p-values between 5.0 × 10⁻⁸ and 1.0 × 10⁻⁵ (wrinkles: 56 SNPs; moisture content: 343 SNPs; pigmentation: 81 SNPs; oil content: 33 SNPs; sensitivity: 59 SNPs). We identified 21 novel, suggestive, clustered loci with p-values between 5.0 × 10⁻⁸ and1.0 × 10⁻⁵. Specifically, we found 6 suggestive loci for the wrinkle trait, 5 suggestive loci for the pigmentation trait, 6 suggestive loci for the oil trait and 4 suggestive loci for the sensitivity trait, as displayed in Manhattan plots (Figure 3). Quantile-quantile (Q-Q) plots were performed for statistical analysis for each phenotype (Figure S2). The results indicated a lambda value of 1.024 in the oil phenotype and 1.043 in the sensitivity phenotype with low genomic inflation.

3.3. Replication Study of the Previous Report

The GWAS of the meta-analysis is shown in

Table 3. Summary of significant skin phenotype loci according to previous GWAS.

									Reported Study (n = 1079)			This Study (n = 261)			Meta−Analysis
						MAF
Trait	SNP	CHR:BP	REF	ALT	Gene	This Study	KOR	ASN	Minor Allele	BETA ± SE a	p	Minor Allele	BETA ± SE a	p	p-Value	I	Q
wrinkle	rs117381658	1:157353684	C	T	FCRL5	0.07	0.04	0.01	T	0.82 ± 0.19	1.17 × 10−5	T	−0.16 ± 0.32	6.06 × 10−1	4.07 × 10−4	86.05	7.40 × 10−3
	rs1961184	10:63733371	G	T	REEP3	0.06	0.05	0.07	T	0.68 ± 0.19	3.09 × 10−4	T	0.29 ± 0.35	4.01 × 10−1	3.38 × 10−4	0	3.26 × 10−1
	rs1929013	1:244230708	G	C	ADSS	0.38	0.35	0.33	G	−0.4 ± 0.09	1.50 × 10−5	G	−0.07 ± 0.16	6.64 × 10−1	6.65 × 10−5	68.83	7.33 × 10−2
	rs7042102	9:92001508	C	T	SPTLC1	0.42	0.44	0.45	T	0.38 ± 0.09	1.94 × 10−5	T	−0.01 ± 0.16	9.44 × 10−1	1.97 × 10−4	78.09	3.26 × 10−2
Moisture	rs9873353	3:31233850	C	T	STT3B	0.02	0.02	0.02	T	−0.55 ± 0.12	3.45 × 10−6	T	0.55 ± 0.26	3.48 × 10−2	7.43 × 10−4	93.29	1.00 × 10−4
	rs34567709	17:61492168	T	G	TBX4	0.12	0.12	0.14	G	0.42 ± 0.09	1.85 × 10−6	G	−0.06 ± 0.09	5.42 × 10−1	2.02 × 10−3	92.78	2.00 × 10−4
	rs1362404	16:51973264	G	T	-	0.25	0.28	0.29	G	−0.32 ± 0.07	5.52 × 10−6	G	0.03 ± 0.08	7.37 × 10−1	1.86 × 10−3	91.15	8.00 × 10−4
	rs7853290	9:71638804	G	A	TRPM3	0.05	0.07	0.09	A	0.55 ± 0.12	2.15 × 10−6	A	0.08 ± 0.16	6.01 × 10−1	3.37 × 10−5	82.77	1.60 × 10−2
	rs143938096	15:79098451	C	A	CTSH	0.08	0.08	0.06	A	−0.54 ± 0.12	3.78 × 10−6	A	0.15 ± 0.11	1.96 × 10−1	2.02 × 10−2	94.41	0
	rs12955989	18:24106190	A	G	TTC39C	0.25	0.21	0.15	G	0.35 ± 0.07	2.67 × 10−6	G	0.1 ± 0.08	1.98 × 10−1	1.60 × 10−5	81.27	2.08 × 10−2
Pigmentation	rs74653330	15:27983407	C	T	OCA2	0.06	0.07	0.03	T	−1.04 ± 0.19	8.91 × 10−8	T	0.16 ± 0.34	6.40 × 10−1	7.86 × 10−6	89.23	2.30 × 10−3
	rs34466224	19:3219644	G	A	NCLN	0.22	0.22	0.22	A	0.65 ± 0.13	3.29 × 10−7	A	−0.03 ± 0.2	8.62 × 10−1	2.07 × 10−5	87.96	3.90 × 10−3
	rs11685354	2:217996408	C	A	TNS1	0.46	0.41	0.45	A	−0.46 ± 0.1	7.99 × 10−6	A	0.32 ± 0.16	4.05 × 10−2	8.56 × 10−3	94.28	0
	rs4653497	1:227355326	T	C	CDC42BPA	0.46	0.49	0.44	T	0.44 ± 0.1	1.46 × 10−5	T	0.11 ± 0.16	4.74 × 10−1	2.78 × 10−5	57.22	1.26 × 10−1
	rs59784607	16:25774628	G	T	HS3ST4	0.17	0.19	0.13	T	−0.54 ± 0.12	1.64 × 10−5	T	−0.19 ± 0.21	3.58 × 10−1	2.84 × 10−5	51.58	1.51 × 10−1
	rs76548385	7:1291682	C	T	UNCX	0.09	0.09	0.07	T	0.79 ± 0.18	1.07 × 10−5	T	0.42 ± 0.27	1.27 × 10−1	5.54 × 10−6	21.4	2.59 × 10−1
Oil	rs308971	3:12075120	G	A	SYN2	0.24	0.26	0.19	G	−0.34 ± 0.07	2.70 × 10−6	G	0.08 ± 0.13	5.31 × 10−1	1.67 × 10−4	88.22	3.60 × 10−3
	rs151209785	18:74549791	T	C	CNDP1	0.03	0.03	0.02	C	−0.57 ± 0.12	4.38 × 10−6	C	0.18 ± 0.32	5.77 × 10−1	4.03 × 10−5	79.03	2.90 × 10−2
	rs9577919	13:113861036	C	T	GAS6	0.06	0.05	0.07	T	0.61 ± 0.14	9.09 × 10−6	T	0.13 ± 0.22	5.60 × 10−1	4.18 × 10−5	71.01	6.33 × 10−2
	rs8107564	19:6964536	A	G	INSR	0.09	0.11	0.09	A	0.42 ± 0.1	2.41 × 10−5	A	−0.02 ± 0.18	9.26 × 10−1	2.33 × 10−4	77.6	3.46 × 10−2
	rs6490805	13:23510670	T	C	TNFRSF19	0.15	0.16	0.15	C	0.36 ± 0.08	7.76 × 10−6	C	0.13 ± 0.15	3.81 × 10−1	1.25 × 10−5	47.08	1.69 × 10−1
	SNP	CHR:BP	REF	ALT	gene	This study	KOR	ASN	A1	OR (95% CI) b	p	A1	OR (95% CI) b	p	p-value	I	Q
Sensitivity	rs7334780	13:106182099	T	C	-	0.3	0.33	0.33	T	0.61 (0.61–0.5)	8.02 × 10−7	T	1.02 (0.64–1.63)	9.38 × 10−1	6.02 × 10−6	74.22	4.89 × 10−2
	rs41308	7:28636081	C	T	CREB5	0.39	0.36	0.35	C	1.58 (1.58–1.3)	3.37 × 10−6	C	0.95 (0.61–1.48)	8.25 × 10−1	1.65 × 10−5	67.07	8.14 × 10−2

Chr: chromosome; BP: base pair; Pos: position; KOR: Korean; ASN: Asian; SE: standard error. The reference MAF was obtained from the KRGDB (http://coda.nih.go.kr/coda/KRGDB/index.jsp (accessed on 6 September 2022)) and Ensembl DB (https://asia.ensembl.org (accessed on 6 September 2022)); Q: Cochran’s Q, I: Heterozygosity. ^a Calculated β and SE by linear regression. ^b We analyzed OR using logistics regression. For the SNPs that were more relevant than previous studies based on GWAS, the cutoff value of p < 1.0 × 10⁻⁵ are underlined.

. Skin phenotypes for the dependent variable in GWAS were defined using the same method as in a previous report. We confirmed that 23 previously identified SNPs, none of them were replicated. In the meta-analysis, we did not identify any genome-wide significant SNPs. However, rs76548385 (UNCX, p = 5.54 × 10⁻⁶) of the pigmentation phenotype showed more significant p-values than in the previous study. It is considered an important marker showing a significant trend in this replication study.

4. Discussion

Few genetic studies have been reported that have examined Korean skin characteristics. The present study identified the genetic markers able to predict facial skin characteristics. As a result, significant new markers were identified in the combined study, and of the several markers identified in previous studies, two SNPs were replicated in this study. These results are expected to contribute to understanding the mechanisms underlying Korean skin characteristics.

As the number of samples increased in the combined study, new significant and suggestive SNPs were detected. We used well-established imputation methods to expand the genetic variation. As a result, we identified 21 new suggestive loci and 6 significantly clustered loci in the combined GWAS. The combined study was conducted because both the study design and sampling investigators were the same in both the previous study and the current study. Among the identified significant SNPs, one SNP functionally related to skin characteristics and the other five SNPs were identified as statistically significant, but no association with skin was reported.

The SNP that was functionally related to skin characteristics was NOP58 (rs149963203), which was encoded as a key component of small nucleolar ribonucleoproteins [42] and is known to be associated with skin pigmentation problems [43,44]. Although the relationship with other skin characteristics has not yet been studied, it is thought to be a candidate that can be sufficiently utilized to predict skin characteristics as a marker that affects pigmentation, one of the main keywords of the skin traits.

SNPs that have not yet been reported in the association to skin characteristics were LY86 (rs201058), ZNF462 (rs16912205), RFLNA (rs56133064), SCIN (rs55999874), and SFI1 (rs78295829), and the contents are as follows.

First, LY86 (rs201058) regulates the lipopolysaccharide (LPS)-mediated signaling pathway [45] and ZNF462 (rs16912205) is known to encode a member of the C2H2-type zinc finger family regulated by several arranged zinc fingers [46]. RFLNA (rs56133064) was predicted to activate filamin-binding activity, including the regulation of bone maturation [47]. However, the relevance of such findings to the skin is unknown. In addition to these SNPs, the roles of both SCIN (rs55999874), which encodes a Ca (2+)-dependent actin filament cleavage protein [48], and SFI1 (rs78295829), a protein essential for cell cycle progression and mitotic spindle assembly [49], in skin are still unclear.

And we looked at whether the newly identified SNPs were also identified in previous studies. As a result, among the newly discovered SNPs, it was confirmed that the p-value of rs113608863, a suggestive SNP, was more significant than the previous study. In the previous study, the cutoff p-value was 1.0 × 10⁻⁵. In the previous study, the p-value of rs113608863, a novel oil trait marker, was 5.52 × 10⁻⁴, which were more significantly confirmed in the combined study. Except for this, no other novel loci and suggestive loci were found in the previous study.

In addition, it was confirmed that newly identified SNPs have not yet been studied in many fields, as shown in Table S2. Although there are few GWAS studies on newly identified SNPs, this study shows a significant p-value, so it is necessary to focus on it as a major skin-related marker. As can be seen from the current research situation, research on the cosmetic phenotype is still in its infancy, and research on relevant key markers is still lacking.

In addition, it was confirmed whether the SNPs identified in the previous study were also validated in this study. As a result, two SNPs with p-values of less than 0.05 were replicated. As a result, two SNPs with p-values < 0.05 were replicated. Two replicated SNPs were rs9873353 (STT3B, β ± SE = 0.55 ± 0.26; p = 3.48 × 10⁻²) and rs11685354 (TNS1, β ± SE = 0.32 ± 0.16; p = 4.05 × 10⁻²). Although the trend was confirmed differently from the previous study, the problem of sample size cannot be excluded, and it may still be a meaningful result; however, further studies are needed in the future.

To provide a detailed description, the previously identified SNP (rs98733534) for the moisture trait was replicated in the current study. STT3B (rs98733534) has not been identified as a functional gene, and the function of the skin is still unclear. However, STT3B showed replication in both previous and current studies on moisture and has potential as an indicator of moisture properties. Also, TNS1 (rs11685354) was replicated in the pigmentation trait, encodes localized to the plasma membrane region where cells attach to the extracellular matrix [50], and was involved in the formation of fibrillar adhesions [51].

To evaluated the consistency or inconsistency of the results, a meta-analysis was performed. We did not identify any genome-wide significant SNPs. However, UNCX (rs76548385) of the pigmentation phenotype showed more significant p-values than in the previous study. It is considered an important marker showing a significant trend in this replication study. UNCX (rs76548385) encodes a homeobox transcription factor involved in somitogenesis. Some homeobox genes regulate melanocyte biology [52]; thus, we consider UNCX a potential candidate for markers associated with pigmentation.

As a limitation of this study, a new SNPs were discovered and several SNPs were replicated, but the same SNPs showed that the opposite effect direction may be due to the insufficient number of samples. However, since meta-analysis was not performed in this study, the results may change after meta-analysis in the future. Given that facial skin measurements are factors that can be affected by various environmental effects, it seems that a number of confounding factors were included in the results of the genetic correlation analysis. Although the phenotype was analyzed by dividing it into tertiles due to such confounding factors, the GWAS for cosmetological phenotypes is still in the initial research stage, and it is thought that the statistical analysis method can be continuously developed. We hope that this study will serve as a starting point for studying the combination of genes and cosmetological phenotypes in the field of personal customization.

This study discovered new significant SNPs in the combined study and confirmed that some SNPs were replicated in the current study and previous studies. However, since the current study was performed using a smaller sample size than the previous study, some data that are not important in this study may be significant in studies with larger sample sizes. Although this study requires additional research, it has the potential to contribute to the development of cosmetic-related genetic research in the future.

5. Conclusions

Skin characteristics are more sensitive to external factors than other clinical features. While this study did not analyze all genetic contributors to skin characteristics, the findings are considered significant because they were newly discovered and replicated in two independent investigations focused on skin phenotypes. These results will serve as meaningful indicators in future functional research.

The results of this study may be used to predict Korean skin characteristics and investigate customized cosmetics. We hope that the markers identified in the larger study group will be replicated in a future study.