Identification of key circadian rhythm genes in skin aging based on bioinformatics and machine learning

Microarray data

Transcriptomic datasets related to skin aging were downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geoprofiles/). Three datasets were selected for analysis: GSE85358, GSE120783, and GSE39170. The GSE85358 dataset, based on the GPL17077 platform (Agilent-039494 SurePrint G3 Human GE v2 8x60K Microarray), included 24 skin tissue samples from young volunteers (Young group) and 24 skin tissue samples from elderly individuals (Old group), and was used as the training set in this study. The GSE120783 dataset, based on the GPL18573 platform (Illumina NextSeq 500), consisted of 17 normal skin tissue samples and 17 skin tissue samples treated with glucocorticoids, a commonly used anti-aging and anti-inflammatory drug in dermatology. This dataset was used as a validation set. In addition, the GSE39170 dataset, based on the GPL9115 platform (Illumina Genome Analyzer II), contained 5 skin tissue samples from young individuals and 5 skin tissue samples from old individuals, and was also used for validation. 260 genes related to the circadian rhythm (CRGs) were downloaded from the MSigDB database (Supplementary Table 1).

Differential expression analysis and WGCNA analysis

Differential expression analysis was performed using the R package “limma” on the GSE85358 dataset to identify genes that were differentially expressed between the old and young groups. Weighted gene co-expression network analysis (WGCNA) was performed using the R package “WGCNA” to identify genes closely related to skin aging. Specifically, we constructed an adjacency matrix based on the expression matrix of differentially expressed genes (DEGs) to describe the strength of the correlation between nodes. An appropriate soft-thresholding value (β) was calculated, and the adjacency matrix was transformed into a topological overlap matrix (TOM). Hierarchical clustering and dynamic tree cut algorithms were used to identify modules. Module eigengenes (MEs), which represent the first principal component of a module, were used to calculate the correlation between the modules and clinical data to identify modules significantly associated with aging (p-value < 0.05). Gene significance (GS) and module membership (MM) were calculated for each module. Hub genes that were significantly associated with clinical traits (Skin aging) were selected from the co-expressed gene modules using the threshold “GS > 0.5 and MM > 0.5 and p-value < 0.05”.

Selection of circadian rhythm genes associated with skin aging

Hub genes identified from the WGCNA network were compared with genes related to the circadian rhythm downloaded from the MSigDB database. The overlapping genes were selected as the circadian rhythm genes associated with skin aging (Skin aging related-CRGs).

Functional enrichment analysis

Gene set enrichment analysis (GSEA) was performed to analyze the enrichment of circadian rhythm-related signaling pathways and biological processes between the Old and Young groups in the GSE85358 dataset using GSEA version 4.0.1 software with KEGG gene sets and GO biological processes gene sets as reference genomes. Enrichment terms with a p-value < 0.05 were considered significant.

The hub genes and skin aging related-CRGs were subjected to GO functional annotation and KEGG functional enrichment analysis using the DAVID database. GO functional annotation involved biological processes (BPs), cellular components (CCs), and molecular functions (MFs). Enrichment terms with a p-value < 0.05 were considered significant.

Protein-protein interaction (PPI) analysis

The skin aging related-CRGs were uploaded to the STRING database (https://string-db.org/) for protein-protein interaction (PPI) analysis. The PPI network was visualized using Cytoscape software.

Screening for key CRGs by machine learning

Two machine learning algorithms, namely LASSO regression analysis and SVM-REF analysis, were employed in this study. The gene expression matrix of skin aging related CRGs was used for the machine learning analysis. Firstly, the key genes were selected based on the LASSO regression algorithm using the R package “rms”. Then, the SVM-REF analysis was performed using the R package “e1071”. The intersection of the genes obtained from both analyses were considered as key circadian rhythm genes that are associated with skin aging.

Immune infiltration analysis

The impact of the immune microenvironment on skin aging was investigated by analyzing the GSE85358 dataset using the CIBERSORTx database (https://cibersortx.stanford.edu/). The gene expression data from each sample were uploaded to the CIBERSORTx database to quantify the infiltration level of 22 immune cell types. We compared the infiltration levels of 22 types of immune cells in skin tissues between the young and old groups using Wilcoxon rank-sum test. Correlations were then examined between immune cells and key CRGs.

Expression of key CRGs

The Wilcoxon rank-sum test was utilized to analyze the expression levels of the identified key CRGs. First, we measured the expression levels of key CRGs in the young and old groups of dataset GSE85358. Subsequently, these profiles were validated in the independent datasets GSE120783 and GSE39170.

Receiver operating characteristic (ROC) analysis

To validate the accuracy of the screened key CRGs and assess the capability of key CRGs to distinguish between aging and young skin, we conducted receiver operating characteristic (ROC) analysis on the training set (GSE85358) and validation sets (GSE120783 and GSE39170), and determined the area under the curve (AUC) values using the “pROC” package.

Screening for genes associated with skin aging

A total of 5726 DEGs were screened from skin tissues of old and young groups in GSE85358, including 3169 upregulated and 2457 down-regulated in the old group (Figure 1A). The expression heat map of DEGs was shown in Figure 1B. The expression matrix of DEGs was used for WGCNA analysis. The WGCNA results were shown in Figures 2 and 3. The optimal soft threshold β was set to 16 (Figure 2A, 2B). Based on the DynamicTreeCut algorithm, with a minimum module size of 30, deep split of 3, and a maximum module distance of 0.25, 9 gene modules were generated, including black, blue, brown, cyan, darkturquoise, greenyellow, grey, pink and royalblue (Figure 2C). We analyzed the connectivity of module eigengenes (MEs), and the results indicated that the distance between the modules was greater than 0.25, indicating that the modules and genes in each module were independent (Figure 2D). The correlation coefficients between these modules and the clinical feature (Skin aging) were calculated. The results showed that all modules were significantly associated with skin aging (Figure 2E). Moreover, according to the criteria of |GS| >0.5 and |MM| >0.5 and p value < 0.05, a total of 1162 hub genes highly associated with skin aging were identified from 9 significant modules (Figure 3A–3I).

To initially explore the potential molecular mechanisms of skin aging, GO and KEGG functional enrichment analysis of hub genes was performed. The top 20 most significantly enriched terms were displayed in Figure 4. GO analysis suggested that these hub genes exerted multiple functions such as protein binding, RNA binding, ubiquitin-protein transferase activity, transcription cofactor activity, RNA polymerase II core promoter proximal region sequence-specific DNA binding (Figure 4A). The hub genes were mainly enriched in many CCs such as cytosol, nucleoplasm, cytoplasm, nucleus, mitochondrion (Figure 4B). And they were involved in several BPs such as rhythmic process, circadian regulation of gene expression, positive regulation of transcription from RNA polymerase II promoter, regulation of circadian rhythm, proteasome-mediated ubiquitin-dependent protein catabolic process (Figure 4C). KEGG enrichment analysis findings revealed that these hub genes were primarily enriched by signaling pathways including circadian rhythm, mRNA surveillance pathway, insulin signaling pathway, ubiquitin mediated proteolysis, RIG-I-like receptor signaling pathway (Figure 4D). Notably, hub genes were significantly enriched to biological functions and signaling pathways associated with circadian rhythms, suggesting that circadian rhythms may play an important role in the skin aging process.

Screening of skin aging-related CRGs

GSEA analysis was used to assess changes in the expression of biological processes and KEGG signaling pathways associated with circadian rhythms in aged skin tissues. The results showed that, compared with the young group of skin tissue, the biological processes and signaling pathways associated with circadian rhythms were significantly downregulated in the old group of skin tissue, according to the GSE85358 dataset (Figure 5A, 5B). This further suggested a strong correlation between circadian rhythms and skin aging. Therefore, we cross-linked the hub gene with CRGs to yield 39 common genes, which are considered skin aging-related CRGs (Figure 6A, Supplementary Table 2). These CRGs mainly participated in various MFs, such as ligand-dependent nuclear receptor binding, promoter-specific chromatin binding, transcription coactivator activity. They were enriched in many CCs including nucleus, nucleoplasm, chromatin, etc. And they were involved in various BPs, such as rhythmic process, regulation of circadian rhythm, circadian regulation of gene expression and circadian rhythm (Figure 6B). The KEGG functional enrichment analysis of skin aging-related CRGs showed that they mainly participated in signaling pathways such as glucagon signaling pathway, insulin resistance, circadian rhythm, thyroid hormone signaling pathway, adipocytokine signaling pathway and PI3K−Akt signaling pathway (Figure 6C, 6D). In addition, the PPI network of skin aging-related CRGs was constructed. As shown in Figure 6E, these genes had a good interaction with each other, and SIRT1, ARNTL, FBXW7, ATF4, SIN3A, and PPARA were the top six genes with the highest degree in the PPI network.

Screening for key CRGs by machine learning

Two machine learning algorithms were used in our study to further screen the key genes from 39 skin aging-related CRGs, which may become potential biomarkers or therapeutic targets for skin aging. Firstly, a LASSO regression model was designed based on samples from aged and young skin tissues. When λ= 0.13, the LASSO regression model was able to accurately evaluate both old and young skin tissues. Therefore, LASSO analysis yielded 6 candidate genes (Figure 7A). Moreover, SVM-REF analysis was performed and identified 3 candidate genes (Figure 7B). Finally, the results from the two algorithms were combined, yielding 3 key CRGs that are associated with skin aging, including SIRT1, ARNTL, ATF4 (Figure 7C).

Immune infiltration analysis

Firstly, we analyzed the infiltration of 22 immune cell types in skin tissues of the young and old groups. The bar graph and heat map showed the type and number of immune cell infiltrates in each sample (Figure 8A, 8B). Subsequently, we assessed the correlation between these immune cell populations (Figure 8C). A strong positive correlation between T cells regulatory (Tregs) and NK cells resting (r = 0.89) and a strong negative correlation between T cells CD8 and T cells CD4 memory resting (r = −0.82) were found.

In addition, we compared the infiltration levels of 22 types of immune cells in skin tissues between the young and old groups. Immune cells such as NK cells activated, Macrophages M1 and Mast cells resting had a significantly higher expression in skin tissue in the old group, whereas T cells CD4 memory activated and Macrophages M2 showed a significantly lower expression in skin tissue in the old group (Figure 9).

Furthermore, we further investigated the correlation of 3 key CRGs with 22 immune cell infiltration (Figure 10A–10C). SIRT1 expression showed a positive link to the infiltration status of T cells CD4 memory activated and negatively correlated with the infiltration status of Macrophages M1, Mast cells resting and NK cells activated. ARNTL expression showed a positive link to the infiltration status of Macrophages M2 and T cells CD4 memory activated and negatively correlated with Mast cells resting, Macrophages M1 and Dendritic cells resting. Moreover, ATF4 expression was positively linked to the infiltration status of T cells CD4 memory activated, B cells naive and Macrophages M2 and negative correlation with the infiltration status of Macrophages M1, Mast cells resting and NK cells activated.

Expression of key CRGs

We further investigated the expression profiles of SIRT1, ARNTL and ATF4 in skin aging. The results showed that in the training set GSE85358, the expression of SIRT1, ARNTL and ATF4 mRNAs were significantly downregulated in the skin tissue of the old group compared to the young group (Figure 11A). Similarly, in the validation set GSE39170, the expression of SIRT1, ARNTL and ATF4 mRNAs were down-regulated in the skin tissues of the old group compared to the young group. Glucocorticoids are the most commonly used drugs for anti-skin aging in clinical practice (Figure 11B). In the validation set GSE120783, glucocorticoid treatment was used to mimic anti-skin aging. In GSE120783, the expression of SIRT1, ARNTL and ATF4 mRNA was significantly up-regulated in the glucocorticoid treated group compared to the control group (Figure 11C).

Receiver operating characteristic (ROC) analysis

To evaluate the diagnostic specificity and sensitivity of key CRGs, we constructed ROC curves and calculated the AUC for each item. In the training set GSE85358, the AUC values of SIRT1, ARNTL, and ATF4 were 0.852, 0.849, and 0.880, respectively (Figure 12A). In the validation set GSE39170, the AUC values of SIRT1, ARNTL, and ATF4 were 0.880, 0.880, and 0.840, respectively (Figure 12B). Moreover, in the validation set GSE120783, the AUC values of SIRT1, ARNTL, and ATF4 were 0.886, 0.962, and 0.737, respectively (Figure 12C). These results indicate that the expression of SIRT1, ARNTL, and ATF4 genes could effectively distinguish aging skin from young skin.