Xinru Hu,1,* Shuang Du,1,* Meng Chen,1 Hao Yang,1 Jia He,1 Lei Zhang,1 Bowen Tan,1 Tao Wu,2 Xi Duan1

1Department of Dermatovenereology, Affiliated Hospital of North Sichuan Medical College, Nanchong, Sichuan, People’s Republic of China; 2Department of Urology, Affiliated Hospital of North Sichuan Medical College, Nanchong, Sichuan, People’s Republic of China

Correspondence: Xi Duan, Department of Dermatovenereology, Affiliated Hospital of North Sichuan Medical College, Nanchong, Sichuan, People’s Republic of China, Email [email protected]

Background: The aging of skin is a diversified biological phenomenon, influenced by a combination of genetic and environmental factors. However, the specific mechanism of skin photoaging is not yet completely elucidated.
Methods: Gene expression profiles for photoaging patients were obtained from the Gene Expression Omnibus (GEO) collection. We conducted single-cell and intercellular communication investigations to identify potential gene sets. Predictive models were created using LASSO regression. The relationships between genes and immune cells were investigated using single sample gene set enrichment analysis (ssGSEA) and gene set variance analysis (GSVA). The molecular processes of important genes were studied using gene enrichment analysis. A miRNA network was created to look for target miRNAs connected with important genes, and transcriptional regulation analysis was used to identify related transcription factors. Finally, merging gene co-expression networks with drug prediction shows molecular pathways of photoaging and potential treatment targets. Furthermore, we validated the role of key genes, immune cell infiltration, and the Adenosine 5‘-monophosphate (AMP)-activated protein kinase (AMPK) pathway in photoaging, which were identified through bioinformatics analysis, using in vivo reverse transcription quantitative PCR (RT-qPCR), immunofluorescence labeling, and Western blotting.
Results: This study discovered three key genes, including Atp2b1, Plekho2, and Tspan13, which perform crucial functions in the photoaging process. Immune cell infiltration analysis showed increased M1 macrophages and CD4 memory T cells in the photoaging group. Further signaling pathway analysis indicated that these key genes are enriched in multiple immune and metabolic pathways. The significant roles of Atp2b1, Plekho2, Tspan13, M1 macrophages infiltration, CD4 memory T cells infiltration and the AMPK pathway in photoaging was validated in vivo.
Conclusion: This research revealed the underlying molecular mechanisms of photoaging, indicating that key genes such as Atp2b1 and Tspan13 play crucial roles in the regulation of immune cell infiltration and metabolic pathways. These findings provide a new theory for the treatment of photoaging and provide prospective targets for the advancement of relevant drugs.

Keywords: skin photoaging, Atp2b1, Tspan13, Plekho2, immune cell infiltration, AMPK pathway

Introduction

The skin is a crucial barrier protecting the body from external threats. Skin aging leads to a decline in structure and function and can be sorted into endogenous and exogenous aging.1 Endogenous aging, or chronological aging, occurs naturally over time. Exogenous aging, primarily photoaging, as a result of prolonged exposure to ultraviolet (UV) radiation, resulting in skin roughness, laxity, sagging, wrinkles, capillary dilation, and pigmentation. Severe cases can lead to skin tumors like melanoma and squamous cell carcinoma.2

The mechanisms associated with photoaging have not been fully clarified. It is discovered that photoaged skin is linked to cellular and extracellular alterations. These changes comprise increased epidermal thickness, disorganized collagen fibers, accumulation of dystrophic elastin fibers, instability of the cellular genome, reduced viability, as well as morphological modifications in keratinocytes and human dermal fibroblasts.2 The molecular process preponderantly contains DNA damage, oxidation stress, inflammatory reaction, alterations in collagen architecture, and cell apoptosis.3 Nucleus kappa-B (NF-κB) signaling pathway, mitogen-activated protein kinase (MAPK) signaling pathway, nuclear factor erythroid 2-related factor 2 (Nrf2) / antioxidant-response element (ARE) signaling pathway, and transforming growth factor-β (TGF-β) signaling pathway are the main pathways associated with skin photoaging. These pathways are interconnected and related to each other.2 The pathogenesis of photoaging remains incompletely understood, making it challenging to recognize patients at high risk for photoaging and related tumors and to develop precise treatments. Dermatologists currently rely on complementary strategies to manage existing photoaging. These strategies are often aggressive, have significant side effects, and are expensive. They include topical treatments with retinoids, laser therapy, fillers, botulinum toxin, and chemical peels.4 Therefore, understanding the genes and mechanisms involved in skin photoaging is crucial for the precise prevention and treatment of photoaging and related tumors, holding significant clinical importance. Bioinformatics analysis is an emerging research method that emerged with the launch of the Human Genome Project and has strong advantages in exploring core genes and biological pathways. In recent years, breakthroughs in histological technologies, such as high-throughput sequencing and microarray technology, have allowed us to explore and identify disease features in depth. In the field of dermatology, it is now becoming an effective method to explore disease pathogenesis.5

In recent years, studies have indicated that certain Chinese herbal medicine compounds or monomers can alleviate photoaging by modulating autophagy. We identified Vicenin-2 (6-C-glucose-8-C-xyloconitin), a flavonoid derived from natural plants and Chinese herbs, as having significant anti-UVB toxicity and radical scavenging effects in human dermal fibroblast experiments. Additionally, Vicenin-2 counteract skin photoaging in mice through antioxidant and anti-inflammatory effects.6

In this study, we used single-cell transcriptomics and bioinformatics methods to deeply analyze gene expression profiles associated with photoaging and validated them through in vivo experiments. First, GSE58915 and GSE173385 data were acquired from public databases, and the single-cell expression profiles were preprocessed and analyzed using the Seurat package. Next, the CellChat package was used to analyze the intercellular communication network and identify key cell subtypes. Characterized genes for photoaging were determined by Lasso regression and random forest algorithm. Immune cell infiltration was further analyzed by the CIBERSORT algorithm, and enrichment pathways were analyzed by GSEA and GSVA. Finally, the molecular mechanism of photoaging and possible therapeutic targets were revealed by combining gene co-expression network and drug prediction. We also established mouse photoaging models to validate the key genes, M1 macrophages infiltration, CD4 memory T cells infiltration and signaling pathways, providing important evidence to further clarify the pathogenesis of skin photoaging.

Methods

In brief, this research began with single-cell analysis. A disease prediction model was developed by using LASSO regression to identify important genes. Following that, immune infiltration and functional enrichment analysis were carried out, a miRNA network was built, major gene regulatory networks were analyzed, and drug-gene interactions were investigated (Figure 1).

Figure 1 Workflow of the study.

Data Acquisition

The GEO database (https://www.ncbi.nlm.nih.gov/geo/info/datasets.html) is a gene expression database created and maintained by the National Center for Biotechnology Information (NCBI), and its full name is GENE EXPRESSION OMNIBUS. The Series Matrix File data file for GSE58915 was downloaded from the NCBI GEO public database. The annotation file is named GPL1261. In total, expression profile data from 21 patients were used, including 6 in the normal group and 15 in the illness group. The single cell data file for GSE173385 was retrieved from the NCBI GEO public database. Additionally, sample data from two patients with complete single cell expression profiles were obtained for single cell analysis.

Single Cell Analysis

First, the expression profile is read through the Seurat software package and low-expression genes are filtered out (nFeature_RNA > 500 and percent.mt < 3*MAD & nFeature_RNA < 3*MAD & nCount_RNA < 3*MAD), where outliers are defined as the same as the median. The three median absolute deviations (MAD) of the package annotates the cluster and is tagged to some cells that are significant in the appearance of the disorder.

Cell Communication

With single-cell data, CellChat is a gadget that can quantitatively infer and analyze intercellular communication networks. To predict the primary signal inputs and outputs of cells as well as how these cells and signals harmonize functions, it uses network analysis and pattern recognition techniques. The standardized single-cell expression profile was utilized as the input data in this work, and the cell information came from the cell subtypes discovered by single-cell research. The activity and impact of each type of cell in the disease was then seen when we examined the cell-related interactions and used the interaction strengths (weights) and counts to determine the degree of closeness of the relationship.

Machine Learning

We employ Lasso logistic regression and random forest algorithms to pick characteristic diagnostic markers of disease. The Lasso algorithm utilizes the “glmnet” package. Random forest is an ensemble learning algorithm with decision trees as its core. It uses the sampling and replacement method to select multiple samples from the sample set as the training set. The sample set obtained through sampling is used to create a decision tree. On each generated node, non-repeating features are randomly selected, and these features are used to divide the sample set respectively. Find the best dividing features and determine the prediction results. This study uses the random forest algorithm to evaluate the importance of features based on %IncMSE, and selects the top 15 features for the following analysis.

Immune Infiltration Analysis

The CIBERSORT method is a widely used method to assess immune cell types in the microenvironment. This technique uses support vector regression to deconvolve expression matrices of immune cell subtypes. It contains 547 biomarkers that distinguish 25 mouse immune cell phenotypes, including T-cell, B-cell, plasma cell, and myeloid cell subsets. This study used the CIBERSORT method to analyze patient data, determine the relative proportions of 25 different types of immune invading cells, and conduct correlation analysis between gene expression and immune cell content.

GSVA Analysis (Gene Set Variation Analysis)

Gene set variation analysis (GSVA) is a non-parametric, unmonitored method for determining the enrichment of gene sets in transcriptomes. GSVA translates gene-level changes into pathway-level changes by exhaustively identifying gene sets of interest and establishing the biological function of the sample. In this study, we will download gene sets from the Molecular Characteristics Database (version 7.0) and apply the GSVA algorithm to fully score each gene set to analyze possible biological function changes in various samples.

GSEA Analysis

Patients were divided into two groups based on gene expression levels: high group and low group. GSEA was then used to investigate differences in signaling pathways between the two groups. The background gene set is a well-known gene set obtained from the MsigDB database for labeling subtype pathways. We continue to investigate differential expression of pathways across subtypes. Significantly enriched gene sets were identified using a concordance score (adjusted p value less than 0.05). Sort. GSEA analysis is often used in studies that combine disease classification with biological significance.

Construction of Key Gene Non-Coding RNA Network

MiRNAs (MicroRNAs) are short non-coding RNAs that have been found to influence gene expression by aiding mRNA degradation or inhibiting mRNA translation. Thus, we investigated if certain miRNAs in essential genes affect transcription or degradation of some hazardous genes. The TargetScanMouse 8.0 (https://www.targetscan.org/mmu_80/) database provided us with miRNAs associated with important genes, which we then visualized using Cytoscape software.

Transcriptional Regulatory Network

This study used the R package “RcisTarget” to predict transcription factors. All calculations performed by RcisTarget are motif-based. The Normalized Motif Enrichment Score (NES) is determined by the total number of motifs in the database. In addition to the motifs identified from the source data, this study inferred additional annotation files based on motif and gene sequence similarity. The first step in determining the overexpression of each motif in a gene set is to calculate the area under the curve (AUC) for each motif-motif pair. This was achieved by calculating a gene set recovery curve with respect to motif order. The NES of each motif is calculated based on the AUC distribution of all motifs in the gene set.

Differential Expression Analysis

The Limma package is an R package for differential expression analysis of expression profiles to identify genes with significantly differential expression between group comparisons. Use the R package “Limma” to analyze differences in the molecular mechanisms of disease data, identifying differentially expressed genes between control and disease samples. Gene differentiation screening conditions were P.Value < 0.05 and |logFC| > 0.585 and draw differentiating gene volcanoes Graphs and heat maps.

CMap Drug Prediction

Connectivity Map (CMap) is an interventional gene expression profile database developed by the Broad Institute; it is mainly used to reveal functional connections between small molecule compounds, genes and disease states. Contains microarray data of 1309 small molecule drugs before and after treatment of 5 human cell lines. There are different treatment conditions, including different drugs, different concentrations, different treatment times, etc. This study predicts targeted therapeutic drugs for diseases based on differentially expressed disease genes.

Animal Experiment

Female BALB/c hairless mice (Sibeifu Beijing) pondering 18–25g were used in this research. All animal experiments were sanctioned by the Animal Care and Use Committee of North Sichuan Medical College and follow the guidance of the Medical and Ethics Committee of North Sichuan Medical College (Sichuan) (NSMC2024-DW-034). All the mice were randomly sorted into the following 3 sets (5 mice per group): (1) control group, (2) model group (UVB+UVA irradiation), (3) Vicenin-2(6-C-glucose-8-C-xyloconitin) group (UVA+UVB+Vicenin-2), We use UVB and UVA to induce photoaging, using the method described before.7 The dorsal skin of BALB/c hairless mice was exposed to a UVB+UVA lamp six per week for 8 weeks to induce photoaging of the skin. The minimal erythema dose (MED) was set at 150 mJ/cm2. The UVA dose was ten times that of UVB. During this 8-week period, apply 0.1mg 0.05% Vicenin-2 ointment topically to the back of the mice 30 minutes after each irradiation. At the end of eight weeks, the photo area of each group was imaged using a digital camera. The cervical dislocation method was employed to sacrifice the mice. Skin tissues were quickly separated from the mice and put in a 4% paraformaldehyde solution or stored at −80 °C for further experiments.

RNA Extraction and Gene Expression Analysis

Total RNA extracted from the skin were evaluated by RT-qPCR. TRIzol was used to extract RNA from the skin. The gene expression level was normalized, and the relative quantification of gene expression was determined using the 2^ (− ΔΔCt) method.

Western Blot Analysis

Skin tissues were placed in lysis buffer and homogenized.The suspension was gathered and subjected to sonication for 1.5 minutes on ice. Left the lysate on ice for 10 minutes. Pre-cooled the centrifuge at 4°C, and centrifuge the lysate for 15 min. Use a BCA protein quantification kit to determine protein concentration. Proteins were separated using SDS-PAGE gels and then transferred onto nitrocellulose membranes. Subsequently, the membranes were incubated in TBST along with 5% skim milk powder at 25°C for 1 hour. Next, the membranes were placed in a solution with the indicated primary antibodies and incubated at 4°C. After that, the membranes were rinsed three times with TBST and then inoculated with a secondary antibody. Use the Bio-Rad Chemi DocXRS+ imaging system to visualize protein bands. Image J (NIH, Bethesda, MD) software was used to quantify band intensities. All assays were performed in five times.

Immunofluorescence Labeling

The protein expression of macrophage CD68 was identified by immunofluorescence single label. The protein expressions of CD68 and CD86 in M1 macrophages, CD68 and CD206 in M2 macrophages, CD4 and CD62L in CD4 memory T cells were detected by double immunofluorescence labeling. The tissue was mixed in 10% neutral buffered formalin. Then it was placed in paraffin and sectioned into slices with a thickness of 5 µm. Once deparaffinized and rehydrated through a series of alcohols, the skin sections were placed in sodium citrate solution and heated in microwave oven to 95 °C for 20 minutes. Subsequently, they were incubated with primary antibody CD68 at a 1:100 dilution (Proteintech, 28058-1-AP) at 4 °C overnight. The slides were then incubated in the same solution for one hour at room temperature on the following day. After removing the primary antibody solution through repeated washes with PBST, the sections were incubated with a Goat Anti-Rabbit IgG (H+L) Secondary Antibody, Cy3 Conjugate (BOSTER, BA1032) and Goat Anti-Rabbit IgG (H+L) Secondary Antibody, FITC Conjugated (BOSTER, BA1032) at a 1:100 dilution at room temperature for one hour. After repeated washes in PBST to eliminate residual antibodies, the nuclei were stained with DAPI at a 1:100 dilution (Biyuntian, C1002).

Double-immunofluorescence labeling was also carried out. Sections were processed as previously described. Subsequently, the sections were incubated with a mixture of two primary antibodies (CD68 and CD86, CD68 and CD206, CD4 and CD62L) at 4 °C overnight at the following dilutions: rabbit anti-CD68 (Proteintech, 28058-1-AP), mouse anti-CD86 (Santa, sc-58986), and mouse anti-CD206 (novus, NBP2-25208), mouse anti-CD4 (Proteintech, 67786-1-Ig), mouse anti-CD206L (Proteintech, 26477-1-AP). Then the sections were incubated with the mixtures of two secondary antibodies (Goat Anti-Mouse IgG (H+L) Secondary Antibody, Cy3 Conjugated (BOSTER, BA1031) and anti-rabbit Alexa Fluor 488 (Invitrogen, A11008) for 2 h at 37 °C at a 1:100 dilution. Nuclei were stained.

Statistical Analysis

In this analysis, R language (version 4.3.0) was used. The results were showing as the mean ± standard deviation (Mean ± SD) by using (ANOVA) one-way analysis of variance software tools. Values of P < 0.05 were indicative of significant differences.

Results After processing the expression profiles with the Seurat package, we identified a total of 12,328 cells, filtered violin plots, scatter plots (Supplementary Figure 2A and B), and the 10 genes with the highest standard deviations are shown (Supplementary Figure 1C). Data were sequentially subjected to normalization, normalization, PCA and concordance analysis (Supplementary Figure 2D–F). The positional relationship between each cluster was obtained through UMAP analysis, and there were a total of 15 cell clusters (Figure 2A). Each subtype was further annotated in this study, with all clusters annotated into 7 cell categories: DiffKera cells, BasalKera cells, macrophages, fibroblasts, NKT cells, endothelial cells, and melanocytes (Figure 2B). Bubble plot of classical markers for 7 cells (Figure 2C) and histogram of cell proportions corresponding to groups (Figure 2D).

Figure 2 Annotation of cells (A) UMAP showing cell clustering results (based on RNA expression, resolution 0.2). Each dot represents a cell, and the color distinguishes different clustering groups. The X-axis and Y-axis are UMAP dimensionality reduction coordinates, indicating the heterogeneity between cells (B) UMAP showing cell type annotation results. Based on the clustering results, the cells are further annotated into different cell types, such as Stem cells, Differentiated cells, Immune cells, etc. Different colors and labels indicate the distribution of different cell types. (C) Bubble maps show the expression patterns of specific genes in various cell types. The X-axis represents genes and the Y-axis represents cell types. The size of the bubble indicates the percentage of cells expressed that the gene was detected in the corresponding cell type, and the color indicates the average amount of gene expression (red for high expression, blue for low expression). (D) Histogram of cell type distribution under different sample conditions. Each column represents a sample, and the color inside the column distinguishes the composition ratio of cell types. The X-axis is the sample condition, the Y-axis is the cell composition ratio, and the numerical value shows the percentage of each cell type. The bar chart shows the P-value of significance, which is used to compare changes in cell type distribution under different conditions.

We used the software package cellchat to examine ligand-receptor relationships characterized in single-cell expression profiles. We found complex interactions between these cell subtypes (Figure 3A and B). Finally, we statistically found potential interactions between macrophages and other cells to be closer (Figure 3C). Therefore, we selected macrophage marker genes as a candidate gene set (Findall.markers.cellType_1.csv). Based on the candidate gene set obtained in the previous step, we used lasso regression and RF algorithms to screen out photoaging characteristic genes. The results showed that Lasso regression identified a total of 12 genes as characteristic genes of photoaging (Figure 4A and B). On the other hand, we evaluated the characteristic genes of photoaging through the RF algorithm (Importance of Variables.txt). The results showed that by screening the top 15 characteristic genes with the highest accuracy in the photoaging data set (Figure 4C), and crossing them with the characteristic genes screened by the Lasso regression algorithm, a total of 3 crossover genes were screened out (Figure 4C and D). Three genes will serve as key genes for our subsequent research: Atp2b1, Plekho2 and Tspan13. The immune microenvironment is predominantly constituted by immune cells, extracellular matrix, a variety of growth factors, inflammatory factors, and specific physicochemical characteristics. Significantly, it exerts a profound influence on disease diagnosis and the sensitivity of clinical treatment. We further explored the possible molecular mechanisms by which key genes affect the progression of photoaging by analyzing the relationship between key genes and immune infiltration in the photoaging data set. The proportion of immune cell content in each patient and the correlation between immune cells are shown in the figure (Figure 5A and B). In addition, the results showed that compared with patients in the control group, patients in the disease group had M1 Macrophage, T Cells CD4 Memory was significantly higher than that of patients in the control group (Figure 5C). We further investigated the connection between key genes and immune cells. The findings indicated that Atp2b1 had a significantly positive correlation with Gamma delta T cells and a significantly negative correlation with M1 Macrophage. (Figure 5D); Plekho2 was markedly positively correlated with T Cells CD4 Memory, and was significantly negatively correlated with Mast There was a markedly negative correlation between T cells and Th17 cells (Figure 5E); there was a markedly positive correlation between Tspan13 and Th17 cells, and a significant negative correlation with T Cells CD4 Memory (Figure 5F). We retrieved the correlation between these key genes and diverse immune factors from the TISIDB database, which encompassed immunosuppressive factors, immunostimulatory factors, chemokines, and receptors. These analyses implied that key genes are closely associated with the degree of immune cell infiltration and play a crucial role in the immune microenvironment (Figure 6A–E). The results are consistent with expectations.Next, we studied the specific signaling pathways enriched in three key genes to explore the possible molecular mechanisms by which the key genes affect the photoaging process. GSEA results showed that Atp2b1-enriched pathways include AMPK signaling pathway, Hippo signaling pathway and other pathways (Figure 7A); Plekho2-enriched pathways include peroxisomes, steroid biosynthesis and other pathways (Figure 7B). Tspan13-enriched pathways included ECM-receptor interactions, lipoic acid metabolism, and other pathways (Figure 7C). In addition, GSVA results show that highly expressed Atp2b1 is enriched in signaling pathways such as TNFA_SIGNALING_VIA_NFKB and UV_RESPONSE_DN (Figure 7A); highly expressed Plekho2 is enriched in signaling pathways such as BILE_ACID_METABOLISM and INTERFERON_ALPHA_RESPONSE (Figure 7B); highly expressed Tspan13 is enriched in HEDGEHOG_SIGNALING, TGF_BETA_SIGNALING and other signaling pathways (Figure 8C). This indicates that key genes might affect the progression of photoaging via these pathways. We used these three key genes as the gene set for this analysis and found that they are controlled by common mechanisms such as multiple transcription factors. Therefore, enrichment analysis of these transcription factors was performed using cumulative recovery curves. Motif-TF annotation and selection analysis results of important genes revealed the motifs with the highest normalized enrichment scores. (NES:5.8) is dbcorrdb__STAT1__ENCSR000DZM_1__m6. We then display all enriched motifs of key genes and corresponding transcription factors (Figure 8A and B). In addition, we performed reverse prediction on 3 key genes through the TargetScanMouse database and obtained 51 miRNAs, with a total of 57 pairs of mRNA-miRNA relationships. We then visualized them using Cytoscape. (Figure 9C). We obtained photoaging-related disease regulatory genes through the GeneCards database (https://www.genecards.org/). We analyzed the expression levels of the 20 genes with the highest Relevance score and found that the expression of Fbn1, Il1r1, Jun, Ptgs2, Rara, Tgfbr2, Trpv1, and Vcan was different between the two groups of patients (Figure 10A). Moreover, we performed correlation analysis on key genes and disease-regulated genes. The expression levels of key genes were significantly correlated with the expression levels of disease-regulated genes. Among them, Tspan13 and Ptgs2 were significantly positively correlated (r= 0.768), and Plekho2 and Ptgs2 were significantly negatively correlated (r =−0.717) (Figure 10B). We then used the limma package to calculate differential genes between disease controls and two groups of patients. The differential gene screening conditions are: P.Value<0.05 and |logFC| > 0.585. Finally, a total of 591 differential genes were screened out, including 130 up-regulated genes. 461 down-regulated genes. Then, the TOP50 genes that were differentially up-regulated and down-regulated in photoaging diseases were divided into two groups, and drug targets were predicted for the differential genes through the Connectivity Map database. The results showed that the expression profiles of drugs such as zebularine, temozolomide, valproic-acid, 9-methyl-5H-6-thia-4, 5-diaza-chrysene-6,6-dioxide were negatively correlated with the expression profiles. Plagued by disease, drugs may be able to alleviate or even reverse the disease state (Figure 11A–D). In addition, we also analyzed the expression of key genes in DiffKera cell, BasalKera cell, Macrophage cell, Fibroblast cell, NKT cell, Endothelial cell, and Melanocyte cell in single cells (Figure 12A and B). Next, we visualized the gene co-expression of photoaging-related disease regulatory genes (Eln, Fbn1, Jun) and 3 key genes in 7 types of cells (Supplementary Figure 2). Next, we used AUCell to quantify immune and metabolic pathways, and analyzed and displayed the correlation between key genes and immune and metabolic pathways (Figure 12C). To verify the above results, we conducted animal experiments. In the photoaging mice models, the irradiated skin showed thickening, deeper wrinkles, dryness, and scaling, indicating successful establishment of the photoaging model. After topical application of Vicenin-2, the skin became thinner, smoother, and more refined. (Figure 13A). The qPCR results showed that Atp2b1 and Tspan13 were significantly different in mRNA expression levels and were consistent with the bioinformatics analysis. However, the Plekho2 did not show any heterogeneity (Figure 13B). Among numerous signaling pathways that are involved in the regulation of autophagy, the AMPK/mTOR (mammalian target of rapamycin) pathway is regarded as important. Therefore, the expression level of the AMPK/mTOR pathway was determined by Western blotting. As the data showed, UV increased the expression of phosphorylated ULK1 but reduced that of AMPK and mTOR.Vicenin-2 enhanced the trends of AMPK/mTOR pathway in photoaging under UVR exposure (Figure 13C–F). The protein expression of macrophage CD68 was measured by immunofluorescence single labeling, The protein expressions of CD68 and CD86 in M1 macrophages, CD68 and CD206 in M2 macrophages, CD4 and CD62L in CD4 memory T cells were detected by double immunofluorescence labeling. The results found that CD4 memory T cells were significantly increased in the photoaging group. The treatment with Vicenin-2 reverses this phenomenon (Figure 13G–H). The expression of CD68 protein was comparable between the model and treatment groups, the expression of CD86 increased in the model group, the expression of CD206 increased in the Vicenin-2 treated group, and the ratio of CD86 / CD206 decreased in the Vicenin-2 treated group. It preliminarily showed that the pro-inflammatory M1 macrophages were more polarized in the skin after irradiation, and the M1 macrophages changed to anti-inflammatory M2 macrophages after Vicenin-2 treatment, and the total number of macrophages was basically unchanged (Figure 13I–N).

Figure 3 Cell communication (A) cell-to-cell interaction network analysis. Left: The Number of interactions between cell types. The nodes represent cell types, the thickness of the lines indicates the number of interactions, and the colors label different cell types. Right: Interaction weights/strength between cell types. The thickness of the line indicates the intensity of the interaction, and the depth of the color indicates the intensity. The analysis demonstrated the centrality and critical role of different cell types in cell communication. (B) Heat maps of cell-cell interactions. The x and y axes represent different cell types. The color of each square indicates the intensity of the interaction between the two cell types (the darker the color, the higher the intensity). This diagram shows the communication characteristics and possible biological connections between specific cell types. (C) Histogram of cell type distribution. The x axis represents the different cell types and the y axis represents the number or proportion of cells. Different colors indicate different cell states or conditions. The graph shows the abundance distribution of various cell types and their differences under a particular sample or condition.

Figure 4 Construction of LASSO and random forest models (A) Variable screening path graph of LASSO regression model. The X-axis represents the logarithmic values of the regularized parameters (Log Lambda), and the Y-axis represents the regression Coefficients (Coefficients). Each color curve represents the variation trend of the regression coefficient of different variables. As the regularization intensity increases (Lambda increases), some variables are gradually shrunk to zero. The dashed line represents the best Lambda value for the cross-validation selection. (B) Cross-validated Mean Squared Error (MSE) curve of LASSO regression model. The X-axis represents the logarithm value of the regularized parameter (Log Lambda), and the Y-axis represents the mean square error of cross-validation. The red dot represents the MSE corresponding to each Lambda value, and the error bar is the standard deviation. The dashed line represents the optimal Lambda value chosen, corresponding to a model with minimal bias or superior complexity. (C) Ranking charts of important characteristic variables. Left: Variable Importance based on random forest algorithm with importance scores on the X-axis and feature names on the Y-axis. Right: Absolute ranking of the variables screened based on LASSO algorithm and their regression coefficients, X-axis is regression coefficient, Y-axis is feature name. The two methods show the sorting results of key variables under different algorithms. (D) Wayne diagrams of important features screened by random forest and LASSO regression. The yellow area represents the features independently screened by the random forest model, the blue area represents the features independently screened by the LASSO model, and the overlapping part is the feature variables jointly screened by the two methods. The diagram shows the similarities and differences between the two approaches to important features.

Figure 5 Immune infiltration (A) Accumulation histogram of major gene expression levels in different groups. The figure shows the relative expression distribution of major genes in each group (Control, Disease), classified by genotype. (B) Heat maps of inter-gene correlation. The chart shows the correlation between major gene expression, showing the strength of the positive correlation (red) versus the negative correlation (blue) in color and correlation coefficient values. (C) Boxmaps of the distribution of gene expression levels in different groups. The X axis of gene name and Y axis of gene expression were compared between Control group and Disease group.(*p<0.05, **p<0.01, **p<0.01). (D–F) Expression significance of three key genes Atp2b1, Plekh02 and Tspan13 in different cell populations. The figure shows statistical significance (P-value) and gene expression trends in different cells, with the size and color of the dots indicating the level of significance and the amount of expression.

Figure 6 Relationship between key genes and immune factors (A–E) Bubble map of correlation between key genes and chemokine, Immunoinhibitor, Immunostimulator, MHC and receptor. Each subgraph shows Pearson correlation coefficient and significance level between different types of genes and key genes. Blue indicates a negative correlation and red indicates a positive correlation.

Figure 7 GSEA analysis of key genes (A–C) Combined results of gene set enrichment analysis (GSEA) and interaction network analysis for three key genes, Atp2b1, Plekh02, and Tspan13. The above section shows GSEA results based on correlated signaling pathways: different curves in each subgraph represent enrichment scores (ES) for correlated signaling pathways (eg Wnt, Hippo, Notch, etc.). The horizontal axis shows the ranking of the gene in the expression data, and the vertical axis shows the cumulative enrichment score. The degree of enrichment of the signaling pathway reveals the biological processes that the target gene may regulate. The lower part is divided into the network diagram of the target gene and other genes: the circular network diagram shows the relationship between the target gene and its significantly related genes and signaling pathways. The color represents the direction of gene expression (red is up-regulated, blue is down-regulated), and the width of the connecting lines reflects the strength of the correlation.

Figure 8 GSVA analysis of key genes (A–C) Gene Set Enrichment analysis (GSVA) of Atp2b1, Plekh02, and Tspan13 genes. Each subgraph shows the sequencing results of functional pathways (up-regulated and down-regulated) that are significantly associated with the target gene, with blue representing up-regulated pathways in which the gene is significantly involved, green representing down-regulated pathways, and gray representing non-significant pathways with low NES.

Figure 9 Key gene-related transcriptional regulation and miRNA network (A) Transcriptional regulatory network of key genes. Green represents key genes and red represents transcription factors. (B) All enriched motifs and corresponding transcription factors of key genes are displayed. (C) miRNA network of key genes, triangles represent mRNA and circles represent miRNA.

Figure 10 Correlation between key genes and disease progression genes (A) Box chart of differential expression analysis: showing the expression distribution of genes in the Control group and the Disease group: The horizontal axis is the gene name and the vertical axis is the expression quantity. The blue box diagram represents the Control group, and the pink box diagram represents the Disease group. “ns” means that the difference is not significant, and an asterisk indicates the level of significance (*p<0.05, ***p<0.001). The results showed that several genes were significantly upregulated or downregulated in the Disease group, suggesting a possible disease correlation. (B) Gene correlation: In the middle is the correlation bubble map between genes, the horizontal axis is the target gene, and the vertical axis is its potential regulatory factors; Bubble color represents Pearson correlation coefficient (red is positive correlation, blue is negative correlation), and bubble size represents significance level. The scatterplots of linear correlations of two genes (eg Timp1 and Tspan13) are shown on the left and right sides, respectively. The lines are fitted curves with correlation coefficients and significance P-values.

Figure 11 Potential therapeutic drugs for diseases (A) 9-methyl-5H-quinolino[8,7-c][1,2]benzothiazine 6.6-dioxide: (A) complex heterocyclic compound containing a quinoline and benzothiazine core structure with a dioxide modification. Its structural characteristics suggest possible pharmacological activity, especially in the field of anti-infection. (B) Temozolomide: a broad class of alkylating agent chemotherapy drugs. Its purine skeleton and carbonyl functional group enable it to cause methylation damage in DNA and play an anti-disease role. (C) Valproic Acid: a drug with a short chain fatty acid structure and a carboxylic acid group. It exerts its therapeutic effects by regulating the release of neurotransmitters and epigenetic mechanisms. (D) Zebularine: a nucleoside analogue with a pyrimidine and ribose structure. It is widely studied as an inhibitor of DNA methyltransferase and plays an important role in epigenetic research.

Figure 12 Expression profile of single cells (A) UMAP showed the expression distribution of Atp2b1, Plekh02 and Tspan13 genes in single-cell transcriptome data. The color depth of the dots represents the level of gene expression, and UMAP_1 and UMAP_2 represent the two principal components after dimensionality reduction. The expression location and intensity distribution of different genes in the cell population reveal their possible cell-specific functions. (B) The bubble map shows the expression of Atp2b1, Plekh02, and Tspan13 genes in major cell types (such as macrophages, endothelial cells, NK cells, etc.): genes on the horizontal axis and cell types on the vertical axis; The bubble size represents the proportion of genes expressed, and the color indicates the average expression level (dark blue is high expression). The results showed that these three genes had an expression advantage in specific cell types. (C) Bubble map of pathway enrichment analysis: It shows the significant enrichment results of Atp2b1, Plekh02 and Tspan13 genes in different functional pathways: the horizontal axis is the functional pathway, and the vertical axis is the gene; The bubble size represents the enrichment significance level (FDR value), and the color represents the enrichment effect direction (red is up-regulated, blue is down-regulated). The classification labels below further group pathways (eg, metabolism, immunity, signaling, etc).

Figure 13 Experimental validation. (A) The treatment effect of Vicenin-2 in photoaging mice model; (B) Relative mRNA expression levels of Atp2b1, Plekho2, Tspan13 in controls, UVR and Vicenin-2 mice (n = 5); (C) Photoaging model was induced UVA and UVB, pre-treated with Vicenin-2 and the expression levels of protein were determined by Western blot; (D) Relative expression levels of p-ULK1/ULK1 (n = 5); (E) Relative expression levels of p-AMPK/AMPK, (n = 5); (F) Relative expression levels of P-MTOR/MTOR (n = 5); (G and H) Immunofluorescence double labeling results of CD4+CD62L+ cells in different groups (n = 5, scale bar =20 μm); (I–N) Immunofluorescence single labeling of CD68 and immunofluorescence double labeling of CD68/CD86 and CD68/CD206 in different groups (n = 3, scale bar =20μm). Data represent mean ± SD, *P < 0.05, **P < 0.01, ***P < 0.001 ****P < 0.0001 ns: no significance.

Discussion

Investigating the genes associated with skin photoaging is of significant importance. For example, understanding the functions and interactions of genes abnormally expressed during the photoaging process can shed light on the underlying molecular mechanisms. Potential biomarkers related to photoaging may be identified, which could be utilized for early diagnosis and therapeutic monitoring. Uncovering genetic differences among individuals can aid in developing personalized prevention and treatment strategies. Moreover, inhibitors or activators targeting specific genes may become potential therapeutic options in the future. In the past few decades, researchers have found some genes associated with photoaging.7 For example, miRNA-22-5p can directly act on growth differentiation factor 11 (GDF11), which is closely related to skin aging. Downregulation of miRNA-22-5p in HDF-EVs can treat skin photoaging by targeting GDF11.8 Wenyi Ho et al discovered that LncSPRY4-IT1--targeting protein were primarily abundant in metabolic pathways, pathways in cancer, endocytosis, human papillomavirus infection and Phosphoinositide 3-kinase (PI3K)/protein kinase B (AKT) signaling pathway which were of great importance in skin photoaging process.9 In our study, using Single-Cell Sequencing Combined with Transcriptome Sequencing, we screened three signature genes Atp2b1, Plekho2 and Tspan13, which are strongly associated with photoaging.

Plekho2 belongs to the superfamily of proteins containing PH domains. Mice lacking Plekho2 show an obvious reduction in the macrophage population. Deficiency of Plekho2 in macrophages leads to increased apoptotic cell death and elevated activity of caspase-3.10 It may participate in the signaling pathways associated with cellular stress, inflammatory response, and indirectly affect the progression of photoaging. Tspan13 is a member of the Tetraspanin family and is involved in various biological processes such as differentiation, cell development, activation, growth factor, motility, proliferation, and metastasis.11 Atp2b1 is among the SERCA Ca2+-ATPases. Currently, its role in the skin is reported to be upregulated along with three other calcium signaling pathway genes in microcystic adnexal carcinoma.12 During photoaging, UV radiation may interfere with the normal transport and balance of Ca+ ions, and this imbalance state can be further exacerbated by changes in Atp2b1 gene expression and function.12 This may lead to a range of cellular stress responses, increased oxidative damage as well as abnormal cell signaling, which accelerating the skin aging process.13 The generation of ROS induced by UVB irradiation is also sensitive to Ca2+. Apoptosis induced by UVB irradiation is caused by increased ROS. Atp2b1 involves many signaling pathways, such as calcium ion signaling pathway, and its abnormal function may trigger a series of downstream signaling pathways related to cell proliferation, differentiation and apoptosis, which can directly affect the development of photoaging.14 These genes not only help us better understand the molecular mechanisms of photoaging but also may provide new biomarkers and targets for the early diagnosis and personalized treatment of photoaging. By further studying the functions and regulatory networks of these genes, we can develop more effective strategies for the prevention and treatment of photoaging. Specifically, Atp2b1, Tspan13, and Plekho2 genes may play a significant role in the pathogenesis of photoaging by regulating key biological processes such as intracellular calcium ion concentration, transmembrane signaling molecule transport, cell-to-cell interactions, and cell apoptosis. Although we found three genes associated with photoaging in our biogenic analysis, Atp2b1 and Tspan13 were confirmed to be associated with photoaging in vivo, which may require larger sample experiments to further clarify this result. Next, we investigated the specific signaling pathways related to the three key genes to uncover the underlying biological mechanisms of photoaging. The main pathways identified through single-gene GSEA of the three characteristic genes were the Hippo signaling pathway, ECM-receptor interaction, and lipid metabolism. The enrichment of these pathways is in line with the clinical fact that photoaged skin is more likely to develop tumors.15–17 This indicates that the three crucial genes may determine the process of cellular photoaging and tumorigenesis by regulating cell growth, intercellular adhesion, and substance metabolism. GSEA pathway enrichment analysis revealed that Plekho2 was enriched in the peroxisome and steroid biosynthesis signaling pathway, which has important roles in lipid metabolism, adipogenesis, and maintaining metabolism. Tspan13 is enriched in the ECM-receptor interaction and lipoic acid metabolism signaling pathway, which is involved in many significant biological processes such as antioxidant defense, differentiation, and apoptosis. Atp2b1 is enriched in the AMPK signaling pathway and Hippo signaling pathway, which is associated with cellular energy balance, metabolic regulation, apoptosis, cell proliferation, and stem cell self-renewal.18 Autophagy is an important cellular recycling system that mediates the turnover of dysfunctional organelles and aggregated or damaged proteins to maintain cellular homeostatic.19,20 Upstream pathways of autophagy include AMPK /mTOR signaling.21 Autophagy regulates DNA damage and inflammatory caused by UVB irradiation, thereby reducing UVB-induced photodamage.22 However, UVB irradiation can also regulate autophagy in skin cells. A previous study showed that autophagy is inhibited in skin keratinocytes irradiated by UVB, suggesting that autophagy induction may be a potential intervention target for anti-photodamage therapy.22 Indeed, autophagy activation is critical for reducing keratinocyte UVB-induced cell death.23,24 We choose and confirmed that AMPK /mTOR pathway is associated with skin photoaging in vivo.

Non-coding RNAs (ncRNAs) are functional RNA molecules that do not translate into proteins. Abnormal expression of ncRNAs is closely associated with the onset and progression of various diseases.25 It has been reported that repeated UVA and UVB irradiation can alter ncRNA expression profiles in human primary skin fibroblasts, keratinocytes, and primary melanocytes, impacting collagen metabolism, extracellular signaling pathways, apoptosis, and inflammatory responses.26,27 Kimball et al conducted transcriptome sequencing on exposed and non-exposed skin in Caucasian women of different ages, identifying mRNA changes related to oxidative stress, energy metabolism, and aging.28 Wang et al found that miR-4732-5p suppresses Tspan13 expression at both mRNA and protein levels, identifying Tspan13 as a direct target of miR-4732-5p.29 Additionally, research has shown that the downregulated genes in the high miR-4469 expression group are mainly involved in inflammatory cell activation, suggesting key targets in miR-4469 target genes, with Plekho2 being a specific target.30 Furthermore, miR-330-5p has been identified as a regulator of Atp2b, affecting cell viability, apoptosis, and cytokine production.31 These findings offer new insights into photoaging treatment.

The skin’s immune system, including mast cells, dendritic cells, lymphocytes and macrophages, is essential for defense beyond the physical barrier.32 These cells maintain skin homeostasis and respond to injury or infection. Macrophages, prevalent in the skin, are classified into M1 and M2 subtypes, influenced by stimuli and environment. The balance of M1 and M2 is crucial for tissue repair and regeneration.33,34

Zhang et al found that aging affects the balance of M1/M2, leading to a significant increase in the pro-inflammatory M1 phenotype in elderly skin, which further leads to chronic skin inflammation.35 It has been showed that M1 macrophages not only have the ability to reduce the expression of type I collagen, but also induce the fragmentation of type I collagen in dermal fibroblast.35 In addition, M1 can also induce the expression of matrix metalloproteinases to accelerate collagen degradation.36 Instead, M2 macrophages promote the expression of collagen and proteins in dermal fibroblasts, thereby inducing the formation of collagen fibers.37 In addition, the significant increase in the M1/M2 ratio is closely related to ultraviolet light-induced skin photoaging.38 Disruption of the balance between M1/M2 macrophages induced by ultraviolet radiation can lead to chronic inflammation, accelerate the aging of dermal fibroblasts, and damage the collagen homeostasis function, thereby promoting skin aging.36 Meanwhile, recent studies have shown that regulating M1/M2 polarization can promote the increase of collagen and elastin fibers.39 After the biogenic analysis found that M1 macrophages increased in patients with skin photoaging, we conducted in vivo to further verify this result and further clarify the M1/M2 relationship. It was found that the polarization of M1 macrophages in the skin of photoaged mice increased, M1 macrophages transformed into anti-inflammatory M2 macrophages after treatment, and the total number of macrophages remained basically unchanged. Therefore, reversing the trend of transformation from M1 to M2 may be the key to improving skin photoaging. In addition, We further explored the relationship between key genes and M1 Macrophage, and the result showed that Atp2b1 was significantly negatively correlated with M1 Macrophage. Atp2b1 may be a potential therapeutic target to control the activity and function of M1 macrophages by regulating calcium ion balance.

We found another type of cell- CD4 memory T cells, increased in skin photoaging. Actually, normal skin harbors recirculating CD4+memory T cells located around superficial dermal blood vessels.40 Studies have found that photoaged and acutely UV-damaged skin is infiltrated with more CD4+ memory T cells, while circulating CD4+ memory T cells are significantly reduced.41,42 This phenomenon of infiltration of memory T cells may be related to the activation of endothelial cells, which express E-selectin, which is a migration medium for lymphocytes to infiltrate skin tissue.43 In addition, after UV irradiation, fibroblasts secrete a large number of senescence-associated secretory phenotype(SASP) factors, including pro-inflammatory cytokines, chemokines, matrix metalloproteinases (MMPs) and growth factors.44 SASP may affect the localization and residence time of CD4 memory T cells in skin by inducing the expression of intercellular adhesion molecules.45 We found that Plekho2 was positively correlated with CD4 memory T cells, and Tspan13 was negatively correlated with CD4 memory T cells. The correlations suggest that they may play an important role in regulating immune memory and T cell function. In the future, investigating the roles of Plekho2 and Tspan13 in endothelial cell activation and fibroblast secretion of SASP, and their effects on T cell infiltration, may enhance our understanding of photoaging development. In drug targeting prediction studies, zebularine, temozolomide, valproic acid and 9-methyl-5H-