Skin aging is characterized by observable major changes in the composition and organization of the skin, including the appearance of wrinkles, tissue sagging, dryness, and pigmentary disorders. These changes are also accompanied by alterations in the underlying skeletal and adipose structures of the face, which further contribute to the overall appearance of aging (Coleman and Rajiv, 2006; Ng and Chew, 2022). Both environmental and genetic factors significantly influence facial aging, and their respective contributions have been studied (Ng and Chew, 2022; Krutmann et al., 2017). While these factors have been studied individually, less is known about how they interact and contribute to the aging process collectively (Krutmann et al., 2017).

The interaction between genetic variations and environmental factors (G × E) have a significant impact on complex human traits and diseases, including skin and facial aging (Krutmann et al., 2017; Herrera-Luis et al., 2024). The skin serves as an excellent model for studying G × E interactions, as it undergoes both intrinsic aging and extrinsic aging in areas exposed to environmental factors, such as the face (Krutmann et al., 2021). Genetic predispositions may influence how skin responds to environmental stressors such as UV radiation and pollution, thus investigating these interactions could provide critical insights into the underlying biological processes of skin aging and overall progression of the aging process (Krutmann et al., 2021).

Among the high number of skin aging phenotype definitions, perceived facial age, a subjective assessment based on an individual’s appearance, stands out as a relevant phenotype for studying G × E interactions due to its inclusion of several features related to skin aging and its association with both genetic and environmental factors (Ng and Chew, 2022; Roberts et al., 2020; Avila et al., 2023; Voegeli et al., 2021). While perceived facial aging is mainly explained by skin aging features (wrinkles, sagging and hyperpigmentation), the muscles, the fat distribution, and the ligaments contribute to the complex underlying anatomy of the face that modulates perception of facial aging (Gunn et al., 2009; Nkengne et al., 2008; Merinville et al., 2015; Swift et al., 2021; Cotofana et al., 2016). Consequently, Perceived age is multifactorial process, not limited to skin changes. Moreover, perceived facial aging is a robust biomarker of ageing that predicts survival in elderly individuals and health (Christensen et al., 2009; Rippon and Steptoe, 2015). It therefore provides a holistic phenotype for G × E studies on facial aging and overall health.

To further explore the role of G × E interactions in perceived facial aging, we assessed perceived age in 226 women at two time points, 12 years apart, allowing defining a perceived facial aging progression (PFAP) as phenotype. Then, we conducted a two-step study focusing on perceived age and sun exposure, allowing reducing the multiple testing burden by focusing on SNP likely to be involved in interactions (Gauderman et al., 2017). In the first step, a genome-wide association study (GWAS) was performed to identify SNP with potential marginal genetic effect on PFAP, and in the second step, we tested these SNP for interaction with sun exposure.

To explore the interaction between genetic variations and sun exposure in perceived facial aging, we conducted a two-step study, using an initial screening to prioritize SNP that are more likely to be involved in an interaction and testing for interaction in a second step (Gauderman et al., 2017). We assessed perceived age in 226 women at two time points, 12 years apart, allowing us to define a perceived facial aging progression (PFAP) phenotype used as facial and skin aging phenotypes.

In the first step, we performed a classical GWAS focused on PFAP on 6,897,523 SNP (MAF >5%; imputation r ² > 0.3). We identified four significant SNP (rs114462394, rs77074982, rs567952380, rs75470026) associated with PFAP (Table 1; Figure 1). The four SNP showed imputation r ² between 0.3 and 0.4, suggesting a high level of uncertainty in the imputed genotypes. However, the MAF of these SNP in our sample were similar to those in the 1000 Genomes European population (Table 1). The genomic inflation factor did not show any deviation (λ = 1.019).

The strongest signal of the GWAS was detected at rs114462394 and rs77074982 SNP, which were in high LD (LD r ² > 0.8) with three SNP (rs114203715, rs114017876 and rs114002886) and they were all located near or in the CGGBP1 and ZNF654 genes on chromosome 3. Among them the rs114017876 and rs77074982 SNP presented several lines of evidence for their involvement in the regulation of gene expression: eQTL effects on the CGGBP1 gene in arterial aorta tissue and the CHMP2B gene in visceral adipose tissue and a high score in RegulomeDB (Table 2), (Sobin et al., 2024). The rs114017876-A allele is associated with a higher PFAP and a lower expression of CGGBP1 while the rs77074982-G allele is associated with a higher PFAP and a higher expression of CHMP2B (Tables 1,2; Figures 2A,B). The second signal was carried by the rs567952380 SNP and was located on chromosome 9 in an intron of the gene PGM5-AS1. The rs567952380-C allele was associated with a higher PFAP in our study and with an increase of PGM5 expression in cultured fibroblasts (Tables 1,2; Figure 2C). The last signal of the GWAS was observed for the rs75470026 SNP which was located on chromosome 8 in an intron of the CSMD1 gene, with the C allele associated with a higher PFAP (Tables 1,2; Figure 2D).

In the second step, only the SNP that passed the GWAS threshold of 1 × 10⁻³ (n = 7,469 SNP) were tested for an interaction with sun exposure regarding the perceived facial aging progression phenotype. In order to determine the Bonferroni threshold, we computed the LD among the 7,469 SNP to evaluate the number of independent tests (LD r ² > 0.8). We identified 1,705 independent SNP, resulting in a threshold of 2.93 × 10⁻⁵. One SNP reached the Bonferroni threshold for the G × E interaction (Table 3).

The only significant signal involved a group of three SNP and one indel in high LD (r ² ≃ 0.8, rs10923186, rs10923187, rs59788048 and rs111978958), with the rs10923186 SNP reaching the Bonferroni threshold (interaction p-value = 2.42 × 10⁻⁵, interaction β = 0.008, FDR = 0.04). The joint test p-value of the signal carried by rs10923186 was 9.23 × 10⁻⁸, nearly reaching the genome-wide significance threshold (5 × 10⁻⁸). To further characterize the interaction signal, we explored the genotype effects on the PFAP across sun exposure categories. These effects were estimated using an additive linear model including genotype, sun exposure, their interaction term, and adjusting for population stratification via PC1 and PC2 (Figure 3A) (R Core Team, 2024). We observed that the rs10923186-T allele favored a higher PFAP in the categories 1, 2 and 3 and, no effect or opposite effect in the categories 4, 5, 6 and 7 (Figure 3A). Overall, the rs10923186-T favored a higher PFAP in less exposed population (Figure 3B). The four variants (rs10923186, rs10923187, rs59788048 and rs111978958) were located on chromosome 1, near or within introns of the VANGL1 gene. In addition to a high score for involvement in gene expression regulation (RegulomeDB), this SNP cluster was involved in an eQTL for CASQ2 gene in artery aorta and the VANGL1 gene in tibial nerve and in a mQTL with the cg20810993 CpG site (Table 2), (Sobin et al., 2024). In summary, through the LD analysis, we observed the rs10923186-T allele was associated with higher PFAP exclusively in the low sun-exposed population, reduced methylation at cg20810993 and higher expression of CASQ2 and VANGL1 genes (Figure 3A; Tables 2,3). In addition to the first signal, we obtained a trend of G × E interaction for the rs10181773 SNP, located on chromosome 2 in an intron of the TRAPPC12 gene. The rs10181773-A allele was associated with a higher PFAP exclusively in the low sun-exposed population (FDR = 0.14, interaction p-value = 1.61 × 10⁻⁴, joint test p-value = 1.31 × 10⁻⁶; Figure 3C; Table 3). The rs10181773 SNP was involved in an eQTL with the TRAPPC12-AS1 gene in skin notably (Table 2), (Sobin et al., 2024).

To explore the role of G × E interactions in facial aging, we conducted a two-step study on 226 women, focusing on perceived facial aging progression with sun exposure as the environmental factor. In the first step, we performed a GWAS to identify marginal genetic effects on PFAP and prioritize SNP for G × E interaction analysis. This initial screening identified four genome-wide significant signals and prioritized 7,464 SNP. In the second step, the 7,464 SNP were tested for interactions with sun exposure in the PFAP phenotype, resulting in one significant SNP and a second one SNP showing a trend towards association.

The strongest signal from the GWAS on the PFAP phenotype was carried by six SNP in high LD (r ² > 0.8: rs114462394, rs77074982, rs114203715, rs114017876, rs1576144130 and rs114002886), all located on chromosome 3, near or within the CGGBP1 and ZNF654 genes. The CGGBP1 gene encodes a CGG repeat-binding protein involved in DNA damage/repair, telomere metabolism, and mRNA metabolism (Singh and Westermark, 2015). CGGBP1 has been associated with telomere protection, and given the known connection between telomere, cellular senescence, aging, and skin aging, it suggests that CGGBP1 could play a role in facial aging via telomere protection (Singh et al., 2014; Rossiello et al., 2022; Nakamura et al., 2002). Notably, the rs114017876-A allele was associated with a higher PFAP and a lower CGGBP1 expression, suggesting that higher CGGBP1 expression could favor telomere protection and slow down facial aging process. The ZNF654 gene is a protein coding gene predicted to be involved in transcriptional regulation and it is expressed broadly across tissues. It is has not yet been involved in skin aging, but it could influence multiple tissue systems relevant to PFAP, including muscles and connective tissues. The second GWAS signal was carried by the rs567952380 SNP, located on chromosome 9 within the PGM5-AS1 gene, and was associated with the expression of PGM5 in cultured fibroblasts. The PGM5-AS1 gene codes for an antisense RNA of PGM5 while the PGM5 itself gene codes for a component of adherens-type cell-cell and cell-matrix junctions (Belkin et al., 1994). Besides, PGM5 was among the top ten upregulated gene expression between progeroid versus wild-type mice and showed a downregulation of PGM5 protein level ex vivo following the use of marine-derived skin care ingredients (Quintana-Torres et al., 2023; Hameury et al., 2019). Our findings were consistent with these observations, as the rs567952380-C allele was linked to increased PGM5 expression and higher PFAP, suggesting that higher PGM5 expression may contribute the facial aging process. The last significant signal was carried by the rs75470026 SNP located on chromosome 8, within the CSMD1 gene. CSMD1 encodes a large membrane-bound protein and was associated with various biological processes and conditions, including schizophrenia, the complement system, cancer, metastasis and cell migration (Akyuz and Bell, 2022). Several studies have linked CSMD1 to aging and skin-related processes, including associations with telomere length, cognitive performance in the elderly, and hypertrophic scars (Son, Cui, and Xi, 2022; Stepanov et al., 2017; Gu et al., 2021).

For the gene and sun exposure interaction regarding PFAP, we obtained one significant signal, carried by the rs10923186 SNP, located on chromosome 1, near the VANGL1 gene. The VANGL1 gene encodes a protein of the tetraspanins family, which contributes to the structural organization of the plasma membrane (Singethan and Schneider-Schaulies, 2008). Interestingly, VANGL1 was identified as a key gene for cutaneous function and was expressed during keratinization according to the Human Protein Atlas (Uhlén et al., 2015; DiTommaso et al., 2014). The keratinization process seemed relevant for skin aging, since this is the process of keratinocyte cells differentiation essential to constitute the stratum corneum, which forms the skin’s barrier against the external environment (Shetty and Gokul S 2012). In summary, our analysis showed that the rs10923186-T allele was associated with higher PFAP in the less sun exposed population, a reduced methylation level of cg20810993 CpG site and higher expression of VANGL1. Regarding the known impact of UV exposure on DNA methylation in the skin and the association of rs10923186 with methylation level and VANGL1 expression level, the absence of association with PFAP in sun exposed population could be explained by UV-induced methylation changes that could lower VANGL1 expression, thereby diminishing its effect on PFAP (Vandiver et al., 2015). Both our result and the possible role of VANGL1 in skin, underlined the keratinization as a relevant candidate for further investigations to better understand the impact of UV exposure on facial aging. In addition to the first signal, a second SNP, rs10181773, showed a trend toward association exclusively in the low sun-exposed population. This SNP was located on chromosome 2 in the TRAPPC12 gene was involved in eQTL effect with TRAPPC12-AS1 gene notably in skin and skeletal muscle. TRAPPC12 encodes a protein that is a component of the TRAPP complex, which play a role in endoplasmic reticulum to Golgi apparatus trafficking (Scrivens et al., 2011; Milev et al., 2017). It also contributes to chromosome congression, kinetochore assembly and stability (Milev et al., 2015).

The absence of association with major genes identified in previous skin aging GWAS, may reflect the difference in phenotype definitions (Ng and Chew, 2022). Most GWAS have focused on individual features of facial skin aging, such as lentigines, wrinkles, or sagging. While some have explored composite phenotypes, to the best of our knowledge, no one has ever investigated a dynamic facial aging phenotype as PFAP. The mechanisms underlying the perceived facial aging progression may involve different and complementary processes compared to those driving specific features like wrinkles. Indeed, these results may reflect the multifactorial process not limited to skin changes and involving components of face structure such as muscles, fat and ligaments (Swift et al., 2021; Cotofana et al., 2016). Several limitations could also account for these discrepancies with prior studies: the relatively small sample size of 226 individuals, the potential bias introduced by self-reported sun exposure questionnaire and the imputation uncertainty of the SNP identified in the GWAS. Sun-exposure questionnaires, however, have shown good concordance between self-reports and objective measures in several studies (Hillhouse et al., 2012; Jennings et al., 2013; Køster et al., 2017).

This study revealed four relevant genetic associations with perceived facial aging progression and one significant G × E interaction. These findings support the known link between telomere shortening/protection and aging, and suggest a potential role for keratinization in the context of facial aging and sun exposure. Although further validation in larger cohorts is necessary, these insights into the molecular etiology of PFAP through genetic approaches may help identify targeted interventions that could address not only cutaneous signs but also underlying facial structures through personalized approaches, notably in the context of a particular environment (Krutmann et al., 2021). Future studies could include functional assays targeting the identified variants (organotypic skin cultures) and integrate multi-omics data to get deeper into the understanding of the mechanisms at stake. Finally, by using a dynamic, holistic measure like PFAP, future genetic investigations may uncover novel biomarkers and mechanisms relevant for facial aging and aging.