A systematic review of skin ageing genes: gene pleiotropy and genes on the chromosomal band 16q24.3 may drive skin ageing

Skin ageing is the result of intrinsic genetic and extrinsic lifestyle factors. However, there is no consensus on skin ageing phenotypes and ways to quantify them. In this systematic review, we first carefully identified 56 skin ageing phenotypes from multiple literature sources and sought the best photo-numeric grading scales to evaluate them. Next, we conducted a systematic review on all 44 Genome-wide Association Studies (GWAS) on skin ageing published to date and identified genetic risk factors (2349 SNPs and 366 genes) associated with skin ageing. We identified 19 promising SNPs found to be significantly (p-Value < 1E−05) associated with skin ageing phenotypes in two or more independent studies. Here we show, using enrichment analyses strategies and gene expression data, that (1) pleiotropy is a recurring theme among skin ageing genes, (2) SNPs associated with skin ageing phenotypes are mostly located in a small handful of 44 pleiotropic and hub genes (mostly on the chromosome band 16q24.3) and 32 skin colour genes. Since numerous genes on the chromosome band 16q24.3 and skin colour genes show pleiotropy, we propose that (1) genes traditionally identified to contribute to skin colour have more than just skin pigmentation roles, and (2) further progress towards understand the development of skin pigmentation requires understanding the contributions of genes on the chromosomal band 16q24.3. We anticipate our systematic review to serve as a hub to locate primary literature sources pertaining to the genetics of skin ageing and to be a starting point for more sophisticated work examining pleiotropic genes, hub genes, and skin ageing phenotypes.


Methodology
Search strategy for skin ageing phenotypes. Identifying skin ageing phenotypes. In this section, we identify and compile a list of changes that occur to the skin which are also widely accepted to be changes that occur due to ageing. We consulted skin ageing scales because they are frequently used to identify and quantify skin ageing phenotypes.
We searched the Web of Science for high quality validated skin ageing scales using the following search term: ((skin ageing OR skin aging) AND (scale) AND (validat*)).
We first included five scales that meet the criteria set by the Consensus-based Standards for the selection of health Measurement Instruments (COSMIN) because they are multidimensional assessment scales determined to have high methodological quality for all skin types. These are the Skin Ageing Score (SAS) 16 , Score of Intrinsic and Extrinsic skin Ageing (SCINEXA) 17 , Merz Aesthetics Scale (MAS) 18 , an unnamed scale by Allerhand et al. 19 and an unnamed scale by Flament et al. 20 . These five scales identified 30 skin ageing phenotypes.
Previously, Dobos et al. 9 have addressed the hundreds of available skin ageing scales that grade one specific phenotype each. Their team performed a systematic review on 100 scales from 130 studies and shortlisted several scales with high methodology, high performance, and focused on interrater reliability. We built upon their findings by studying each shortlisted scale independently and checking the way each phenotype is defined and comparing them against the definitions of the 30 phenotypes in COSMIN. This yielded several duplicates which indicate a degree of agreement across the scales. We including 'forehead wrinkles' and 'low brow positioning' from the upper face ageing scale by Flynn et al. 21 , 'glabellar wrinkles' from the scale by Honeck et al. 22 , 'upper lip fullness' , 'lower lip fullness' , and 'sagging of jawline' from the lower face ageing scale by Narins et al. 23 , 'melomental folds' from the scale by Carruthers et al. 24 , and 'neck sagging' from the scale by Sattler et al. 25 . Altogether, these five scales added eight more skin ageing phenotypes.
Altogether, these sources identified 56 skin ageing phenotypes (Table 2), which we use as our definition of skin ageing phenotypes in all subsequent discussions. The 56 skin ageing phenotypes are discussed in detail in "Results" "Skin ageing phenotypes" section. Briefly, the 56 skin ageing phenotypes can be broadly grouped into four key phenotype categories by their morphology-Category A Phenotype are skin wrinkling and saggingrelated phenotypes; Category B Phenotype are skin colour-related phenotypes; Category C Phenotype are skin cancer-related phenotypes; and Category D Phenotype are skin global impression phenotypes.
Search strategy for skin ageing genes. Identifying skin ageing genes containing SNPs associated with skin ageing phenotypes. This review was conducted in accordance with the Preferred Reporting Item for Systematic Review and Meta-Analyses (PRISMA) guidelines (Supplementary Table S1). A search specifically for candidate gene papers as opposed to GWAS did not return records for wrinkles, and only returned records on skin cancer [30][31][32] and pigmentation changes 7,33,34 . All records on wrinkles are from GWAS. Therefore, a primary literature search on GWAS was performed using the Embase, PubMed and Web of Science databases in July 2021. Search results were restricted to English journal articles published any time before and including 2021. The search term for all databases included 'skin aging' or 'skin ageing' in the title or abstract, and 'GWAS' and 'genotype' in all index fields. To be exhaustive for East Asian results, we performed another primary literature search. The search term for all databases included 'skin aging' or 'skin ageing' in the title or abstract, and 'Chinese' and 'SNP' in all index fields. Full search terms for their respective databases are summarised in Table 1. To remain current, searches were reperformed in November 2021. Eligible articles from the primary search were determined using pre-defined eligibility criteria. A secondary search was conducted by manually searching references cited by the eligible articles from the primary search. Results obtained in the secondary search were deduplicated and screened using the same eligibility criteria as that applied in the primary search. The manual search process was repeated for results from the secondary search to ensure a thorough record was obtained (Fig. 1). A final search using specific skin ageing phenotypes as keywords only returned records that were duplicates of previous results or irrelevant to this review's aims, ensuring that the literature search illustrated in Fig. 1 was sufficiently thorough.
Eligibility criteria. We included articles into this review in accordance with our definition of skin ageing and stated aims. The included articles examined association of SNPs with skin ageing phenotypes, both intrinsic and extrinsic, in any human subject assessed via non-invasive methods, by means of a non-experimental observational study. All full-text journal articles available in English were included. Studies excluded from the review met at least one of the following criteria: study pertained to non-human subjects (e.g., swine or dog experiments); study did not examine any SNP; phenotypes examined by the study were disease or disorders (e.g., Parkinson's Disease, psoriasis); study pertained to epidemiological risk factors (e.g., smoking, nutrition); study focused on cosmetic products or treatments; studies were reviews or meta-analyses. Phenotypes quantified using digital image analysis were included. Using this eligibility criteria, a total of 240 records were screened. Screening was performed in two stages. In the first stage, articles with titles and abstracts irrelevant to this review's aims were discarded. In the second stage, full-length journal articles were retrieved for all remaining papers and each record was manually assessed for suitability. Full-text articles excluded from the review were either review papers or meta-analyses papers.
Data extraction. The following variables were obtained from the respective articles and supplementary materials: author(s) and year, country, sample size and subject demographics, genetic background (SNP rsID, effect allele, allele model), and phenotype details (phenotype, anatomical location, and definition if any). In several studies, participants were not recruited and only photographs of their faces were investigated; the number of photographs included in statistical analysis, if reported, was extracted, and considered as the sample size. In www.nature.com/scientificreports/ several studies, the number of recruited participants and the number of participants with genetic data available to the investigators differed; the number of participants included in statistical analysis, if reported, was extracted, and considered as the sample size. In several studies, the same SNP was tested for associations with multiple phenotypes and only one result has a significant p-Value (p-Value < 1E−05); to identify between phenotypes that have been tested for this SNP and phenotypes that have never been tested for this SNP in any literature, all results are retained regardless of their p-Value. To account for varying terminology used by different studies for identical phenotypes, we compared the definition of the phenotype in words or photographs, if any, across different papers and collapsed them into the same phenotype. For quantitative data, the following were extracted: effect allele frequency, p-Values, estimates of effect sizes, beta coefficient, odds ratios, confidence intervals, standard errors, whichever was reported. Where possible, beta coefficients and standard errors were converted to odds ratios and confidence intervals respectively. Some SNPs were reported in SNP coordinates from an old human assembly; to obtain the most recent information, all SNP coordinates were lifted over from their original assembly to the hg38 human assembly using UCSC LiftOver, and both rsID and genome coordinates are reported for every SNP to maximise usability of information. Some discrepancy exists between our cleaned dataset and the genes reported in primary literature because new knowledge indicated that the reported SNP belongs to a different gene instead of what had been reported in the primary literature. To ensure that the reassignment of SNPs to genes is sufficiently thorough, all published SNPs were converted to genomic coordinates using UCSC Table  Browser, loaded onto a custom track and matched to the most updated gene positions using UCSC Genome Browser on the hg38 human assembly. Some papers report ssID instead of rsID; the genome coordinate of the ssID SNP is first tracked in the old human assembly, lifted over to the hg38 human assembly, then, the new genome coordinate is tracked in the hg38 human assembly to identify the rsID. Both the rsID and hg38 human assembly genome coordinate are reported instead of the ssID and old genome coordinates. Some SNPs are intergenic; the genome coordinate of the SNP is reported together with the closest landmark (e.g., protein coding gene, non-coding gene, splice variant, pseudogene) upstream and downstream of the SNP. Some SNPs are in uncharacterised genes; the Ensembl State ID (ENSG) is reported as no gene symbol has been assigned yet. Some gene symbols change over time; to maximise usability of information, both the current gene symbol and Ensembl Stable ID (ENSG) are reported for every gene.
Statistical analysis. SNP-phenotype associations reported in only one cohort were excluded from pooling; the results were still grouped at the phenotype level and reported. The random-effect meta-analysis was performed by first extracting the Odds Ratio (OR) and 95% Confidence Interval (CI) reported in all the studies of interest. For papers with no OR and CI reported, we calculate them using the reported statistical beta and Standard Error (SE) respectively. As most associations in the field were reported based on the additive model (i.e., each copy of the allele changes the OR by a certain amount), in this review, we estimated the effect size based on the additive model; all results from other models (e.g., dominant model, allele model) are not included in the figures as the results from different allele models cannot be directly compared. However, all results from other models are still reported for completeness of information. Meta-analyses were conducted on R with the RStudio interface, using the metafor package 35 . Meta-analysis for each risk factor (i.e., each SNP) was conducted when at least two independent studies or cohorts report an effect size estimate for an association between a skin ageing phenotype of interest and the said SNP. We used the DerSimonian and Laird procedure for random effects meta-analysis 36 to account for between-study heterogeneity, as what had previously been done 1 . Pooled odds ratios (pOR) were calculated for meta-analysis of ≥ 2 studies or cohorts. To avoid a possible overlap of study populations and to avoid including duplicates in the meta-analytic models, each study was only represented once for a given SNP unless the study reported effect sizes after stratifying their populations. We also avoided including the same entry twice (i.e., duplication) such as when two studies investigated the same phenotype in the same study cohort. In such cases where there could be overlapping study subjects, only the larger study was considered. To assess heterogeneity presented in the pooled risk estimate, we computed the inconsistency index (the I 2 index). An I 2 value ≥ 50% and a p-Value < 0.05 was considered statistically significant for heterogeneity. Begg's funnel plots and Egger's test were used to assess publication bias; Begg's funnel plots were only drawn for meta-analysis of ≥ 2 studies or cohorts and Egger's test was only performed for meta-analysis of ≥ 3 studies or cohorts to estimate asymmetry of data; an Egger's test p-Value < 0.05 was considered as that publication bias possibly exist.
Identifying pleiotropic genes. Pleiotropic genes are genes associated with two or more skin ageing phenotypes. We defined pleiotropy strictly such that a gene is pleiotropic if it is associated with morphologically distinct skin ageing phenotypes (e.g., forehead wrinkles and pigment spots) but not pleiotropic if the associated phenotypes are morphologically similar (e.g., forehead wrinkles and Crow's feet wrinkles). There are 3703 SNP-phenotype associations that pass GWAS-suggestive significant threshold or better (p-Value < 1E−05). We mapped them to 366 genes and intergenic regions on the hg38 human assembly. Genes containing SNPs associated with two or more morphologically distinct skin ageing phenotypes are treated as pleiotropic genes.  www.nature.com/scientificreports/ Gene enrichment analyses. As described in "Identifying skin ageing phenotypes" section, We grouped 56 skin ageing phenotypes into four key phenotypic categories (A, B, C, and D) based on their morphology. Skin ageing genes were stratified into the four categories according to gene-phenotype associations. Next, we subjected (1) all skin ageing genes, (2) Category A genes, (3) Category B genes, (4) Category C genes, and (5) Category D genes to three different gene enrichment strategies (general enrichment, functional enrichment, and network enrichment). This ensures that the gene enrichment analyses were performed sufficiently thoroughly.
General enrichment. General enrichment was performed using the functional profiling feature (g:GOSt) on g:Profiler. Lists of genes identified by their Ensembl Stable ID (ENSG) were queried for enrichment from multiple databases-Gene Ontology (GO), Biological Pathways (KEGG, Reactome, WikiPathways), regulatory motifs in DNA (TRANSFAC (TF), miRTarBase), protein databases (Human Protein Atlas, CORUM), and Human phenotype ontology (HP). Detailed results were retrieved and all enrichments between two or more genes were reported in Supplementary Table S9.
Functional enrichment. Functional enrichment was performed using Version 11.5 of the Search Tool for Retrieval of Interacting Genes/Proteins (STRING) database (https:// string-db. org/) (Figs. 3, 5, 6). Lists of genes identified by their Ensembl Stable ID (ENSG) were queried for enrichment in Homo sapiens with interaction confidence score > 0.4. Proteins that may not physically bind to each other but jointly contribute to a shared function were identified from curated databases and co-expression experiments. Other sources, namely textmined data (i.e., co-mentions in publications) and co-occurrence of the proteins in different genomes (if any) were also reported. The network is further clustered using the Markov Cluster (MCL) Algorithm-a fast and scalable unsupervised cluster algorithm-with a MCL inflation parameter of 3.
Network enrichment. Network enrichment identifies common nodes with high connectivity in the network (i.e., hub proteins). High connectivity among proteins could illuminate key pathways driving the development of key skin ageing phenotypes such as wrinkling and sagging skin.
Network enrichment was performed using Version 2.0.0 of Molecular Complex Detection (MCODE) available as a plug-in of Version 3.9.0 of Cytoscape (National Institute of General Medical Sciences, Bethesda, MD, USA) (Fig. 4). The network data was first exported from STRING to Cytoscape and MCODE plug-in of Cytoscape was employed to identify the functional modules in the system. Locally dense protein-protein interaction (PPI) networks are scored according to their graphical interconnectivity 37 . The parameters used to identify these networks are a degree cutoff of 2, a node density cutoff of 0.2, a K-Core of 2, a maximum depth of 100, 'Haircut' turned off, and 'Fluff ' turned off. A sensitivity analysis was conducted to determine the best node score cutoff among 0.1, 0.2, and 0.3. Results with the best Network Scores are shown and the corresponding node score cutoff is reported. Modules with a Network Score ≥ 4 and nodes ≥ 6 are further analysed for general enrichment. Hub proteins were queried for enrichment from multiple databases-Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway, Gene Ontology (GO), and Transfac (TF) through the g:Profiler gene set analysis toolkit (g:Profiler; https:// biit. cs. ut. ee/g: Profi ler/ gost). The full set of network and general enrichment results including the results at different node score cutoffs and results with Network Scores below 4 are reported in Supplementary Table S10 Fig. S6) 43 , and GSE17046: Analysis of the Infrared-A radiation induced transcriptome of human skin fibroblasts ( Supplementary Fig. S7) 44 . We further stratified the 'chronological age' and 'UVR exposure' risk factors to improve the accuracy of our associations. Specifically, 'chronological age' was stratified into 'aged sun-protected skin' , 'young adult sun-protected skin' , and 'infant sun-protected skin' , while 'UVR' was stratified into 'single exposure' , 'repeated exposure' , and 'repeated exposure followed by recovery' .
We extracted the raw data and made comparisons across the GEO datasets according to our aims. Comparisons were made between the expression levels in aged skin tissues from older participants and expression levels of the same gene in young skin tissues from younger participants for three datasets (GSE85358, GSE108674 and GSE181022). Comparisons were made between the expression levels in skin tissue exhibiting a particular skin ageing phenotype (e.g., melanoma skin or lentigo skin) and the adjacent normal skin of the same participant for two datasets (GSE44805 and GSE109778). Comparisons were made between cells exposed to a treatment (e.g., UV radiation or infrared radiation) and untreated cells for two datasets (GSE102676 and GSE17046).
The compared data were reported as a ratio of the average expression level in the cases (i.e., aged skin tissue or skin tissue exhibiting skin ageing phenotypes or cells exposed to a treatment) against the average expression www.nature.com/scientificreports/ level of the reference group (i.e., young skin tissue or adjacent normal skin or untreated cells). Ratios above 1 indicate higher expression levels in the cases group while ratios below 1 indicate lower expression levels in the cases group. Significant (p-Value < 0.05) expression ratios were retrieved for further analysis.

Results
SNPs associated with skin ageing phenotypes: literature search. Forty-four eligible publications were identified during the primary search, secondary search, and the subsequent screening process (Supplementary References and Fig. 1). As the shape of the Beggs' funnel plots are mostly symmetrical, and Egger's test p-Value are mostly > 0.05, there does not appear to be obvious publication bias across all phenotypes. The risk factor studied in this paper is genetics, specifically, we compare the association between the effect allele of a given SNP on a given skin ageing phenotype against the association between the non-effect allele of that SNP on the same phenotype. Pleiotropic SNPs exist (i.e., the effect allele of a SNP (e.g., rs258322) is significantly associated with multiple different skin ageing phenotypes (p-Value < 1E−05). There are 2349 effect SNPs making 3703 SNP-phenotype associations with skin ageing phenotypes at the GWAS-suggestive significance level (p-Value < 1E−05). Table S2). Sample sizes and subject characteristics varied between studies, with samples sizes ranging from 87 females from Japan and Germany only 45 , to 287,137 genetic samples from both genders from the 23andme database 46 . The studies varied in scope, by examining specific races (e.g. Ashkenazi Jews in New York 47  Skin ageing phenotype definition. Most publications investigated one phenotype (e.g., forehead wrinkles). Publications which investigated multiple phenotypes often use different methods to identify different phenotypes. For example, Liu et al. 57 detects three skin ageing phenotypes-wrinkles, facial pigmented spots, and perceived age. 'Wrinkles' are detected both by image analysis algorithms and manually assessed using photonumeric scales. 'Facial pigmented spots' are detected by image analysis algorithms. 'Perceived age' data is collected from manually assessing the facial appearance of participants. The detailed breakdown and definitions used are in Supplementary Table S3. Briefly, nineteen papers defined the phenotype by describing the appearance of the feature 15,33,48,49,[52][53][54][57][58][59][60][61][62][63][64][65][66][67][68] , three papers employed image analysis algorithms to identify the feature 56,57,69 , fourteen papers identified the feature using photonumeric scales 14,45,47,51,52,57,[69][70][71][72][73][74][75][76] , seven papers detected the feature through self-reporting via questionnaires 46,50,67,[76][77][78][79] , three papers detected the feature through equipment 55,80,81 , two papers studied genomic data collected for a different condition (e.g. the inability to tan) and attempted to use it for their phenotype (e.g., carcinoma) 54,82 , and the way that the feature was defined was not reported in five papers 30,33,61,83,84 .

Overview of study characteristics. Study populations originated from 16 countries (Supplementary
Overall, most methods of identifying skin ageing phenotypes required visual assessment of the appearance of the skin feature either by one or more dermatologists or the research subject. The presence/absence of the skin feature and its severity are reported through text-based questions (e.g., 'What type of skin cancer did you have? Please check all that apply. ' 50 ) or photonumeric grading scales. A few studies 56,57,69 performed digital image analyses to obtain more objective measures (e.g., wrinkle length, area). There is agreement with Wong and Chew 1 that even though SCINEXA was validated as a scale which calculates an overall skin aging score from the individual scores of 23 skin ageing signs, papers reported individual phenotype scores, instead of overall skin aging score. As skin ageing is almost always assessed as individual phenotypes, in this review, we examine individual skin ageing phenotypes as proxy measures of overall skin ageing.

Skin ageing phenotypes.
There is currently no complete list of skin ageing phenotypes; the most recent attempt is the Skin Ageing Atlas 26 a decade ago. Using the method detailed in "Identifying skin ageing phenotypes" section, we are the first to bring together the findings from the Skin Ageing Atlas and multiple other sources (e.g., the 44 GWAS publications in the skin ageing field) to create a comprehensive list of all the ways that skin ageing phenotypes have ever been collected and measured in the skin ageing field ( Table 2) across ethnicities worldwide.
Fifty-six skin ageing phenotypes were identified using methods detailed in "Identifying skin ageing phenotypes" section. Next, we group the 56 skin ageing phenotypes into four key phenotypes categories based on their morphology-Category A Phenotype are skin wrinkling and sagging-related phenotypes; Category B Phenotype are skin colour-related phenotypes; Category C Phenotype are skin cancer-related phenotypes; and Category D Phenotype are skin global impression phenotypes. We therefore treat skin ageing as the result of changes to the skin in these four key areas. While the appearance of all 56 phenotypes have been described in literature, SNP-phenotype associations have not been explored for all 56 phenotypes. Supplementary Table S4 contains a  detailed breakdown of the Table 2 phenotypes which also appeared in SNP literature.
We observed that while the field has an abundance of SNP-phenotype associations that reach GWAS-suggestive significance level (p-Value < 1E−05) or better in one publication, only 19  For the former three phenotypes, checks against the three databases (Web of Science, PubMed, Ensembl) also retrieved no GWAS papers on solar elastosis, perceived skin darkness, and an older perceived age other than the three papers cited in this illustration 57,67,72 . For the fourth phenotype, while our literature search retrieved GWAS papers on facial pigmented spots, solar lentigines, and seborrheic keratosis, none of the retrieved papers 14,15,45,51,58,60,61,[64][65][66]70,73 , reported a significant association with rs4785704 other than the sole paper cited in this illustration 69 .
This illustration shows that there is some uncertainty in establishing whether a given SNP is truly associated with skin ageing. The associations made in the discovery cohort are often either less significant in the replication cohort or absent in other studies investigating the same phenotype. We therefore hypothesise that skin ageing phenotypes are not so much the result of individual SNPs but rather, is the overall effect of multiple SNPs on the same gene or intergenic region. We therefore focus on studying the genes associated with each skin ageing phenotype (i.e., Category A, B, C, and D Phenotypes) and how the genes and proteins are interrelated, through enrichment studies and differential gene expression studies.
We quantify the association between the SNPs and the phenotype using Odds Ratios. The additive allele model is the preferred model in the field; out of the 2349 SNPs making 3703 SNP-phenotype associations with skin ageing phenotypes at the GWAS-suggestive significance level (p-Value < 1E−05), the additive allele model accounts for 2073 of these SNPs and 3135 of these SNP-phenotype associations. Hence, the Odds Ratios for all SNPs illustrated in the figures follow the additive allele model.
We report, in Fig. 2a-c, the 19 SNP-phenotype associations which are corroborated by other studies at the GWAS-suggestive significance level (p-Value < 1E−05) or better. Other than these 19 SNP-phenotype associations, there are another 1,060 SNP-phenotype associations from publications which performed internal validation by having a discovery cohort and validation cohort(s) (Supplementary Dataset S1) and another 2056 SNP-phenotype associations found at the GWAS-suggestive significance level (p-Value < 1E−05) in only a single publication (Supplementary Dataset S2). All SNP-phenotype associations, regardless of allele model, are reported in Supplementary Table S5. Pooled odds ratios (pOR) were calculated for random effects meta-analysis of ≥ 2 studies or cohorts where appropriate; Fig. 2d illustrates the forest plot for one of the 19 SNPs (rs258322), the rest are reported in Supplementary Dataset S3. Figure 2e illustrates the Begg's funnel plot for rs258322 and a skin ageing phenotype (solar elastosis on the hand dorsum); the rest are reported in Supplementary Dataset S4.
Skin ageing genes are pleiotropic. We define pleiotropic skin ageing genes as genes associated with multiple morphologically distinct skin ageing phenotypes (i.e., Phenotypes A, B, C, and D). To the best of our knowledge, we are the first to put forth the idea that pleiotropic skin ageing genes account for the bulk of the genetic associations with skin ageing.
We are also the first to bring together the findings from all 44 GWAS publications in the skin ageing field to create a comprehensive list of genes found to be associated with skin ageing phenotypes across ethnicities worldwide. Being the first to have a comprehensive picture of the skin ageing gene landscape, we seized the opportunity to find relationships and patterns among the skin ageing genes that could not have been discovered by researchers who lack this comprehensive list.
We report that most of the pleiotropic skin ageing genes are skin colour genes, genes on the chromosomal band 16q24.3, and their immediate intergenic neighbourhood ( Fig. 3 and Supplementary Table S6). Specifically, the intergenic region between the IRF4 skin colour gene and the gene directly downstream (i.e., ENSG00000286364) is highly pleiotropic. Similarly, the intergenic region between SPATA2L and VPS9D1 on the chromosomal band 16q24.3 is highly pleiotropic. Details on the extent of pleiotropy of skin ageing genes and SNPs can be found in the Supplementary Tables S7 and S8 respectively.
We earlier proposed that (1) genes traditionally identified to contribute to skin colour have more than just skin pigmentation roles, and (2) further progress towards understand the development of skin pigmentation requires understanding the contributions of genes on the chromosomal band 16q24.3. The first proposal is well-illustrated by how the MC1R gene (a skin colour gene) could have more than just skin pigmentation roles (Fig. 3). The second proposal is well-illustrated by the fact that the majority of the 23 genes on the chromosomal band 16q24.3 are associated with skin colour-related phenotypes (colour-coded in blue in Fig. 3) despite only one gene (MC1R) being traditionally identified as a skin colour gene.
Skin colour genes are genes first described for their role in human skin pigmentation. While our understanding of the genetics of skin pigmentation is relatively clear, the genetics for skin wrinkling and sagging is still largely unknown. Recent discoveries propose that skin colour genes are associated with more than just skin colour-related phenotypes 14,63 and may drive the development of these phenotypes through melanogenesisindependent pathways which are still poorly understood 63 . Our list of pleiotropic skin ageing genes includes many skin colour genes and can guide future work in characterising the melanogenesis-independent pathways in which skin ageing genes operate, which may ultimately explain their association with pigmentation-independent phenotypes such as facial wrinkles.  www.nature.com/scientificreports/ We report that many skin colour genes are pleiotropic for skin ageing phenotypes: SHC4 is pleiotropic for Category A Phenotypes (skin wrinkling and sagging phenotypes) and Category B Phenotypes (skin colour phenotypes). BNC2, HERC2, OCA2, RALY, and TYR are pleiotropic for Category B Phenotypes (skin colour phenotypes) and Category C Phenotypes (skin cancer phenotypes). IRF4 and SLC45A2 are pleiotropic for Category A Phenotypes (skin wrinkling and sagging phenotypes), Category B Phenotypes (skin colour phenotypes), and Category C Phenotypes (skin cancer phenotypes). MC1R is pleiotropic for Category A Phenotypes (skin wrinkling and sagging phenotype), Category B Phenotypes (skin colour phenotype), Category C Phenotypes (skin cancer phenotype), and Category D Phenotypes (skin global impression phenotypes).
Skin colour and Category B Phenotypes (skin colour phenotypes) are dependent on melanin synthesis by tyrosinase. One form of melanin, phaeomelanin, generates free radicals in response to ultraviolet radiation that in turn, damages the skin 85 . It is possible that some non-skin colour phenotypes (e.g., skin wrinkling and sagging) may result from the skin damage sustained through this process.

Enrichment analyses.
We scrutinised the dataset using three different enrichment analyses-general enrichment (Supplementary Table S9), functional enrichment, and network enrichment (Supplementary Table S10). All three enrichment strategies revealed several novel and coherent results which offer different perspectives to the same story.
Firstly, many skin ageing phenotypes are associated with a small handful of 44 genes. These 44 genes are (1) associated with morphologically diverse skin ageing phenotypes (i.e., each of them is a pleiotropic skin ageing gene) (Fig. 4 and Supplementary Table S6) or (2) code for hub proteins: highly interconnected proteins in the gene network whose frequent interactions could illuminate relevant pathways driving the development of key skin ageing phenotypes (Fig. 5), or (3) are pleiotropic genes coding for hub proteins (Supplementary Table S11).
Secondly, most of the pleiotropic and hub genes fall into one of two categories-skin colour genes (Fig. 6), or genes on the chromosomal band 16q24.3 (Fig. 7). This gives two immediate implications: First, genes traditionally identified to contribute to skin colour take on more than just skin pigmentation roles, and second, to truly www.nature.com/scientificreports/ understand the development of skin pigmentation, we must also account for the contributions by genes on the chromosomal band 16q24.3.
Realising the important role that skin colour genes and genes on the chromosomal band 16q24.3 contribute to understanding many diverse skin ageing phenotypes, we searched the literature to obtain a thorough list of 32 skin colour genes (Supplementary Table S12). This list contains traditionally well-studied skin colour genes (e.g., BNC2), and less well-studied skin colour genes (e.g., KIAA0930) 52 . However, we expect that there will be even more skin colour genes because we discovered that many genes on the chromosomal band 16q24.3 are associated with skin pigmentation phenotypes but are not yet classified under skin colour genes by the larger community. Hence, to be as thorough as possible, we paid close attention to our established list of 44 pleiotropic skin ageing genes and 32 skin colour genes.
Gene expression levels change with age. Defining case and control groups. Aged skin and young skin have different expression levels in the pleiotropic skin ageing genes, hub genes, and skin colour genes. We retrieve this expression data from the Gene Expression Omnibus (GEO) repository and represent the data as a ratio of the average expression level in the cases (e.g., aged skin tissue or skin tissue exhibiting skin ageing phenotypes or cells exposed to a treatment) against the average expression level of the reference group (i.e., young skin tissue or adjacent normal skin or untreated cells). Ratios above 1 indicate higher expression levels in the cases group while ratios below 1 indicate lower expression levels in the cases group. Expression ratios of pleiotropic skin ageing genes, hub genes, and skin colour genes are reported in Tables 3, 4, 5, S11, S12, and Supplementary Figs. S1-S7.
As described in "Gene expression omnibus (GEO) repository analysis" section, we extracted the raw data and made relevant comparisons across the GEO datasets. Firstly, among the four GEO datasets on skin ageing phenotypes, three of them (GSE85358, GSE108674, and GSE44805) have significantly different expression levels between aged skin and young skin and these results will be further discussed. Secondly, among the three GEO datasets used to study various risk factors of skin ageing, two of them (GSE181022 and GSE102676) have significantly different expression levels between skin exposed to the risk factor and skin unexposed to the risk factor and these results will be further discussed in "Chronological ageing is associated with the downregulation of BNC2, ENSG00000198211, MC1R, SHC4, SLC24A5, SPIRE2, and TYR in the aged skin", "Chronological ageing is associated with the upregulation of AGR3, DSTYK, and TPCN2, and the downregulation of BNC2, EDNRB, MC1R, SHC4, and TYR in the aged skin", and "UV radiation exposure is associated with the downregulation of DBNDD1 in the aged skin" sections. www.nature.com/scientificreports/ Table 3. Differentially expressed pleiotropic and hub genes associated with chronological age. a GSE181022 controls for Sun exposure by using only sun-protected skin. The skin tissue in GSE85358 is skin aged normally, which will typically be exposed concurrently to numerous factors (e.g., chronological age and sun exposure). b The expression ratio is the average expression level in the aged tissue or tissue with skin ageing phenotypes (i.e., cases) divided by (/) the average expression level in the young tissue or tissue without skin ageing phenotypes (i.e., reference). Ratios above 1 indicate higher expression levels in the cases group while ratios below 1 indicate lower expression levels in the cases group.  Table 3 and illustrated using SPIRE2 in Fig. 8. The full set of results are in Supplementary Table S11. BNC2 (GSE85358, Expression Ratio = 0.98, p-Value = 5.57E−03) and ENSG00000198211 (GSE85358, Expression Ratio = 0.96, p-Value = 8.10E−03) are downregulated in aged skin tissue. This finding is coherent with study GSE181022 which found that BNC2 (Expression Ratio = 0.90, p-Value = 3.58E−03) is downregulated in Sunprotected non-infant skin tissue as compared to infant skin tissue, and that ENSG00000198211 (Expression Ratio = 0.93, p-Value = 3.18E−02) is downregulated in Sun-protected young adult skin tissue as compared to aged skin tissue. This suggests that chronological ageing from infant to old age is associated with the downregulation of BNC2 observed in aged skin and that chronological ageing from young adult to old age is associated with the downregulation of ENSG00000198211 observed in aged skin.
Melanoma is a skin ageing phenotype. Study GSE44805 reports that SHC4 (Expression Ratio = 0.70, p-Value = 2.07E−03), SLC24A5 (Expression Ratio = 0.37, p-Value = 3.77E−03), TYR (Expression Ratio = 0.39, p-Value = 1.65E−08), and MC1R (Expression Ratio = 0.89, p-Value = 1.13E−02) are downregulated in melanoma skin relative to adjacent normal skin. This finding is coherent with study GSE181022 which found that SHC4 (Expression Ratio = 0.90, p-Value = 4.58E−08), SLC24A5 (Expression Ratio = 0.91, p-Value = 1.53E−02), and TYR (Expression Ratio = 0.74, p-Value = 6.63E−03) are also downregulated in Sun-protected non-infant skin tissue as compared to infant skin tissue and MC1R (Expression Ratio = 0.96, p-Value = 1.23E−02) is also downregulated in Sun-protected young adult skin tissue as compared to aged skin tissue. This suggests that chronological ageing from infant to old age is associated with the downregulation of SHC4, SLC24A5, and TYR observed in aged skin and that chronological ageing from young adult to old age is associated with the downregulation of MC1R observed in aged skin. SPIRE2 (Fig. 8) is both downregulated in aged skin tissue (GSE85358, Expression Ratio = 0.96, p-Value = 1.82E−02) and in melanoma skin relative to adjacent normal skin (GSE44805, Expression Ratio = 0.77, p-Value = 4.84E−02). This finding is coherent with study GSE181022 which found that SPIRE2 is also downregulated in Sun-protected non-infant skin tissue as compared to infant skin tissue (Expression Ratio = 0.98, p-Value = 5.69E−03). This suggests that chronological ageing from infant to old age is associated with the downregulation of SPIRE2 observed in aged skin.
Chronological ageing is associated with the upregulation of AGR3, DSTYK, and TPCN2, and the downregulation of BNC2, EDNRB, MC1R, SHC4, and TYR in the aged skin. In this section, we analyse gene expression level changes to our established list of 32 skin colour genes. Seven GEO Datasets relevant to our study aims were consulted. Genes with significantly different (p-Value < 0.05) expression levels between case and control groups in ≥ 2 GEO Datasets are reported in Table 4 and the full set of results are in Supplementary Table S12. The downregulation of BNC2, MC1R, SHC4, and TYR have been explored in the previous section as they are also pleiotropic genes and/or hub genes in addition to being skin colour genes. Melanoma is a skin ageing is upregulated in Sun-protected non-infant skin tissue as compared to infant skin tissue. This suggests that chronological ageing from infant to old age is associated with the upregulation of AGR3 observed in aged skin.
DSTYK (GSE85358, Expression Ratio = 1.03, p-Value = 2.86E−02) and TPCN2 (GSE85358, Expression Ratio = 1.02, p-Value = 1.89E−02) are upregulated in aged skin tissue. These findings are coherent with study GSE181022 which found that DSTYK (Expression Ratio = 1.21, p-Value = 3.75E−05) and TPCN2 (Expression Ratio = 1.05, p-Value = 1.24E−05) are upregulated in Sun-protected non-infant skin tissue as compared to infant skin tissue. Additionally, study GSE102676 found that DSTYK (Expression Ratio = 1.12, p-Value = 5.78E−03) and TPCN2 (Expression Ratio = 1.21, p-Value = 3.98E−02) are also upregulated in HaCaT keratinocyte cells exposed to a single dose of 12 J/cm 2 solar stimulated UVR. These suggest that both chronological ageing from infant to Table 4. Differentially expressed skin colour genes associated with chronological age. a GSE181022 controls for sun exposure by using only sun-protected skin. The skin tissue in GSE85358 is skin aged normally, which will typically be exposed concurrently to numerous factors (e.g., chronological age and sun exposure). b The expression ratio is the average expression level in the aged tissue or tissue with skin ageing phenotypes (i.e., cases) divided by (/) the average expression level in the young tissue or tissue without skin ageing phenotypes (i.e., reference). Ratios above 1 indicate higher expression levels in the cases group while ratios below 1 indicate lower expression levels in the cases group.  Table 5. Differentially expressed pleiotropic and hub genes associated with UV radiation exposure. a GSE181022 controls for sun exposure by using only Sun-protected skin. The skin tissue in GSE85358 is skin aged normally, which will typically be exposed concurrently to numerous factors (e.g., chronological age and sun exposure). b The expression ratio is the average expression level in the aged tissue or tissue with skin ageing phenotypes (i.e., cases) divided by (/) the average expression level in the young tissue or tissue without skin ageing phenotypes (i.e., reference). Ratios above 1 indicate higher expression levels in the cases group while ratios below 1 indicate lower expression levels in the cases group. c Significant values (p-Value < 0.05) are bolded. UV radiation exposure is associated with the downregulation of DBNDD1 in the aged skin. Chronological age is not the only risk factor for skin ageing. Other risk factors include UV radiation exposure (Table 5). DBNDD1 is both a pleiotropic gene and a hub gene. DBNDD1 is significantly associated with UV radiation exposure, not chronological ageing (Table 5). Melanoma is a skin ageing phenotype; DBNDD1 is downregulated in melanoma skin relative to adjacent normal skin (GSE44805, Expression Ratio = 0.84, p-Value = 2.10E−02). This finding is coherent with study GSE102676 which found that DBNDD1 is also downregulated in HaCaT keratinocyte cells exposed to 5 doses of 12 J/cm 2 solar stimulated UVR followed by a 1-week recovery period (Expression Ratio = 0.65, p-Value = 4.10E−02). There are no significant associations between Sun-protected infant, young adult, or aged skin tissue (p-Value for comparison between infant and non-infant = 7.82E−02, p-Value for comparison between young adult and old age = 4.64E−01). This suggests that even with a sufficient recovery period, short periods of frequent UVR exposure, not chronological age, is associated with the downregulation of DBNDD1 observed in aged skin.

Limitations
While there is a diversity in assessment methods of skin ageing phenotypes across all papers, the bulk of them still rely on visual assessment methods. On one hand, image analysis algorithms such as MATLAB used in some studies 57 would provide superior analysis capabilities because skin ageing features will be indiscriminately detected and will have less bias. Nonetheless, these image analysis algorithms have technical challenges of its own; already, one paper reported difficulties in separating one type of facial pigmented lesions (solar lentigines) from many other types (seborrheic keratosis, melanocytic nevi, freckles, melasma) using their algorithm 69 . Furthermore, obtaining the high-resolution digital photographs required for an accurate analysis is challenging and may require specialised equipment, making it difficult to be recommended by default.
On the other hand, visual assessment methods are cost-effective but are inherently subjective and can be influenced by cultural perceptions, for instance, Chinese women tend to underestimate the age of Chinese faces while overestimating the age of Caucasian faces 86 .
Often, the association between a given SNP and a given skin ageing phenotype could not be replicated. This could be caused by several reasons, some of which are listed below: (1) many skin ageing phenotypes are understudied, or (2) the heritability of the SNP differs in different populations in different environments, or (3) skin ageing is the cumulative effect of multiple SNPs of the same gene.

Conclusion
There is currently no consensus on the definition of skin ageing, no comprehensive list of skin ageing phenotypes, and no clear ways to identify them. Such a list is important for an overview of the morphologies in which skin ageing phenotypes manifest and to identify potential associations among different phenotypes. We are the first to create this comprehensive list of skin ageing phenotypes. Creating this list was challenging due to two reasons. Firstly, there are multiple discrepancies in the phenotypes that different researchers consider important to measure, and secondly, there are over 100 scales of different methodological qualities claiming to quantify the same Figure 8. Gene expression of the pleiotropic and hub gene SPIRE2 is downregulated in non-infant skin, aged skin, and melanoma skin, when compared against expression levels in infant skin, young skin, and nonmelanoma control skin respectively. Melanoma is a skin ageing phenotype. Expression ratio is the ratio of the average expression level in the cases (e.g., aged skin tissue) against the average expression level of the reference group (e.g., young skin tissue). Ratios above 1 indicate higher expression levels in the cases group while ratios below 1 indicate lower expression levels in the cases group. www.nature.com/scientificreports/ item (e.g., forehead wrinkles). To achieve reliable data, we only selected the best photonumeric grading scales (i.e., the scales with high methodology, high performance, and focused on interrater reliability).
To the best of our knowledge, nobody in the field has compiled a list of skin ageing genes. We brought together the findings from all 44 GWAS publications in the field to create a comprehensive list of genes found to be associated with skin ageing phenotypes. Being the first to have a detailed picture of the skin ageing gene landscape, we found relationships and patterns among the skin ageing genes that could not have been discovered by researchers who lack this comprehensive list.
Through three different enrichment analyses (general enrichment, functional enrichment, and network enrichment) and gene expression data from GEO Datasets, we discovered four new findings: (1) pleiotropy is a recurring theme among skin ageing genes, (2) SNPs associated with skin ageing phenotypes are mostly found in a small handful of 44 pleiotropic and hub genes and 32 skin colour genes, (3) various gene expression studies paint a coherent picture that most of these genes are down-or up-regulated with chronological age, and (4) most skin ageing genes are spatially confined within chromosomal band 16q24.3, or are skin colour genes elsewhere in the genome.
As a result, we propose two important implications. First, genes traditionally identified to contribute to skin colour take on more than just skin pigmentation roles, and second, to truly understand the development of skin pigmentation, we must also account for the contributions by genes on the chromosomal band 16q24.3.
Specifically, through our systematic review we identified 19 promising SNPs which are found to be significantly (p-Value < 1E−05) associated with skin ageing phenotypes in two or more independent studies. Multiple lines of evidence (e.g., SNP, gene expression levels, pleiotropy, enrichment) were employed to support and strengthen the association between a named gene (e.g., SPIRE2) and skin ageing.
Finally, this systematic review acts as a hub to locate primary literature sources pertaining to the genetics of skin ageing and offers new perspectives (i.e., skin ageing genes are pleiotropic and most of which are either spatially confined within chromosomal band 16q24.3 or are skin colour genes elsewhere on the genome) to understanding skin ageing as a condition with multiple morphologically distinct but nonetheless related phenotypes.
Future work will examine skin ageing phenotypes more closely and explore environmental risk factors in Han Chinese skin. Moreover, since many skin colour genes are pleiotropic, future work in characterising the melanogenesis-independent pathways in which skin ageing genes operate may ultimately explain their association with pigmentation-independent phenotypes such as facial wrinkles.

Data availability
Supplementary Information is available for this paper. Correspondence and requests for materials should be addressed to F.T.C. Reprints and permissions information is available at http:// www. nature. com/ repri nts. The forty-four publications selected for the review were identified through the primary search or the secondary search and remained relevant to the aims of this review after a subsequent screening process. The publications are reported in Supplementary References. All 3703 SNP-phenotype associations with skin ageing phenotypes at the GWAS-suggestive significance level (p-Value < 1E−05), 1,079 Begg's funnel plots, and Egger's test used to assess publication bias, and Odds Ratios are documented in Supplementary Table S5, Supplementary Datasets S1, and S2 respectively. A review protocol was not prepared prior to the writing of this manuscript but the review methodology follows what was described in "Methodology" section.