Introduction

Aging is a complex biological process characterized by a gradual decline in physiological functions and is closely associated with the onset of various degenerative diseases. Studies have shown that health conditions vary significantly among elderly individuals, even among those of the same chronological age [1, 2]. This variation may be attributed to differences between chronological age and biological age. Individuals whose biological age exceeds their chronological age are more susceptible to aging-related diseases, including cardiovascular disorders, cancer, Alzheimer’s disease, progeria, and even premature mortality [3]. Therefore, chronological age alone is insufficient for accurately predicting disease risk and guiding clinical interventions. In contrast, biological age serves as a better predictor of health status and aging-related disease outcomes [4].

Aging clocks utilize omics-derived features within linear or nonlinear frameworks to infer chronological age and aging speed. To quantify aging, aging clocks often rely on AgeDiff—the deviation between biological and chronological age estimates [5]. For example, an ensemble classifier model developed by Jason G. Fleischer et al. predicted human age from dermal fibroblasts with a mean absolute error (MAE) of 7.8 years [6]. A transcriptomic aging clock based on human peripheral blood achieved a mean absolute deviation (MAD) of 4 years [7]. The performance of DNA methylation clocks revealed an MAD of 4.9 years in peripheral blood and 3.6 years across various tissues [8]. An immune aging clock based on 50 blood cytokines resulted in an MAD of 15.2 years [9]. Individuals with higher AgeDiff (fast agers) tend to show accelerated physiological aging, whereas those with lower AgeDiff (slow agers) exhibit delayed age-related decline.

The age-related changes measured by biological clocks, and the resulting variations in metrics like AgeDiff, are direct reflections of fundamental biological processes. Indeed, aging is characterized by several molecular and cellular hallmarks, including genome instability, telomere erosion, epigenetic drift, loss of proteostasis, defective autophagy, nutrient sensing imbalances, mitochondrial decline, cellular senescence, stem cell exhaustion, impaired intercellular communication, chronic inflammation, and microbial dysregulation, with cellular senescence considered a central mechanism underlying systemic aging [10].

The interplay of these aging hallmarks contributes to the functional decline observed across different organs. As the largest organ in the human body, the skin is affected by intrinsic aging processes and simultaneously acts as a frontline defense against environmental aggressors [11], particularly ultraviolet (UV) radiation [12]. UV radiation is classified based on wavelength into UVA (320–400 nm), UVB (280–320 nm), and UVC (100–280 nm) radiation. UVC is absorbed by the stratospheric ozone layer and does not reach the Earth’s surface [13]. Although UVB radiation penetrates the skin less than UVA, it causes more severe skin damage [14]. UVB irradiation contributes to the development of various dermatological abnormalities, including pigmentary disorders and inflammatory reactions [15]. Prolonged UVB exposure leads to photoaging over time, manifested by wrinkling [16], loss of barrier integrity and elasticity, and an increased risk of skin cancers such as keratinocyte carcinoma, basal cell carcinoma, squamous cell carcinoma, and melanoma [17, 18].

Extracellular vesicles (EVs) are nanosized, lipid bilayer-bound particles secreted by cells that carry various bioactive components, including lipids, proteins, and nucleic acids such as microRNAs (miRNAs) [19, 20]. EVs play a pivotal role in intercellular communication, regulating both physiological and pathological processes through local and long-distance signaling [21]. Embryonic stem cells (ESCs), known for their unlimited proliferative capacity, pluripotency, and intrinsic resistance to senescence [22], are considered promising candidates for treating aging and aging-related diseases [23]. While the broader anti-aging potential of hESC-EVs is an emerging area of research, their specific effects on photoaged skin and the underlying mechanisms remain to be fully elucidated. Therefore, this study investigated the effects of stem cell-derived EVs on photoaged dermal fibroblasts.

In this study, we developed a rapid transcriptome-based DNN system to assess cellular senescence levels. We also deployed the pretrained model onto a web-based platform to enhance accessibility and usability. Furthermore, we isolated EVs from embryonic stem cells and evaluated their anti-aging effects on UVB-induced photoaged fibroblasts. The results demonstrated that hESC-EVs exerted significant anti-senescence effects, as reflected by the DNN model. Overall, our SkinAGE system demonstrated robust performance and potential applicability in regenerative medicine.

Methods

Data sources

Six transcriptomic datasets were obtained: four (GSE113957, GSE226189, GSE51518, GSE119009) from the Gene Expression Omnibus (GEO)(https://www.ncbi.nlm.nih.gov/geo/) and two (E-MTAB-7032 and E-MTAB-3037) from ArrayExpress(https://www.ebi.ac.uk/arrayexpress/)(Table 1) (Table 1) [24]. The first four datasets in Table 1, comprising a total of 300 samples, were combined to form the training set for our machine learning model. The E-MTAB-3037 dataset was used as the external validation set, while the GSE119009 dataset was used to demonstrate the application of the SkinAGE tool. To address interdataset variability, gene identifiers were first standardized to official gene symbols via the org.Hs.eg.db R package, followed by batch effect correction with the ComBat function from the sva R package. Skin aging-related transcriptomic data were obtained from the Human Ageing Genomic Resources (HAGR) [25].

Table 1 Summary of the transcriptomic datasets used
Full size table

Screening for Differentially Expressed Genes (DEGs)

To generate an initial list of high-confidence candidate genes, we utilized the integrated public dataset. We compared 59 samples from younger individuals (≤ 30 years old) with 128 samples from older individuals (≥ 60 years old). Differential expression analysis was performed using the limma R package (v3.40.6). To focus only on genes with the most substantial age-associated changes, we applied a highly stringent threshold: an adjusted p-value < 0.05 and an absolute fold change > 4 (equivalent to |log2FC| > 2). This stringent criterion was intentionally chosen to drastically narrow down the pool of potential features for model building.

In the second stage, we analyzed an in-house RNA-seq dataset generated for this study. For this dataset, we applied a standard threshold (adjusted p-value < 0.05 and absolute fold change > 1.5) to conduct a more comprehensive differential expression analysis. This allowed us to validate the trends observed in the public data and to explore a wider range of moderately changed genes within our own experimental context. Volcano plots were used to visualize the DEGs from both analyses.

Construction of WGCNA network and module detection

We performed WGCNA using the R package [26]. Genes were filtered by retaining the top 75% ranked by median absolute deviation (MAD > 0.01), and eight outlier samples were removed by hierarchical clustering. We used the pickSoftThreshold function to select β = 8 for a scale-free topology. This was converted to a Topological Overlap Matrix (TOM) to reduce noise. Modules (min 10 genes) were identified from the TOM’s dissimilarity (dissTOM) using a dynamic tree-cutting algorithm. Module eigengenes were calculated, similar modules were merged (height cutoff 0.25), and Pearson correlation (P < 0.05) identified the module most strongly associated with age.

Identification of hub genes via LASSO regression and SVM

To identify the most robust age-discriminative gene features, we implemented two machine learning selection methods: LASSO-penalized logistic regression [27] and Support Vector Machine-Recursive Feature Elimination (SVM-RFE) [28]. Both models were trained on the full set of genes remaining after low-expression filtering. The models performed a binary classification task, with the dependent variable being the class label distinguishing “young” (age ≤ 30) and “old” (age ≥ 60) samples. Optimal model parameters were tuned via cross-validation: a 10-fold cross-validation for the LASSO regularization parameter (λ) and a 5-fold cross-validation for the SVM-RFE feature subset evaluation. The final set of hub genes was determined by the intersection of features selected by both optimally-tuned models, as identified by a Venn diagram [29].

Development and validation of machine learning models

In this study, we developed an aging clock model using nine different machine learning algorithms to accurately predict biological age. The models included Bootstrap Aggregating (Bagging), Light Gradient Boosting Machine (LightGBM), Extreme Gradient Boosting (Xgboost), Elastic Net Regression (Enet), Logistic Regression (LR), DNN, Random Forest (RF), Gradient Boosting Decision Tree (GBDT) and SVM.

Following data harmonization and normalization of the entire dataset, we employed stratified random sampling to split the data into 80% for the training set and 20% for the testing set. To ensure robust model evaluation and generalizability, model performance was further validated using five-fold cross-validation on the training set.

Model performance comparison

To gauge the effectiveness of each model, we focused on three key metrics: the mean absolute error (MAE), the root mean square error (RMSE), and the coefficient of determination (R²). The MAE measures the average magnitude of the errors between the predicted and actual values. The RMSE reflects the standard deviation of prediction errors and is particularly sensitive to larger deviations. The coefficient of determination (R²) serves as a metric for the proportion of variance in the observed data that is explained by the model, where a greater R² is indicative of a more robust model fit [30].

All three metrics were calculated on both the training and testing datasets to comprehensively assess each model’s fitting accuracy and generalizability. Based on the performance across these evaluation criteria, we selected the DNN model as the final predictive tool, owing to its superior accuracy and robustness.

Model interpretability

Interpreting machine learning models remains a significant challenge. Shapley additive explanations (SHAP), a method grounded in game theory, addresses this issue by ranking input features according to their importance and providing explanations for model predictions, thereby mitigating the “black box” nature of ML models.

By quantifying the contribution of each feature to specific predictions and the overall model output, SHAP facilitates both local interpretability and global interpretability. In this study, SHAP values were employed to assist in feature selection [31]. Based on feature importance rankings, we employed a stepwise reduction approach to select the optimal number of features, evaluating counts ranging from 32 to 10. After evaluating the performance of the models at each stage, we determined that using all 32 features yielded the best predictive performance and thus retained them for the final model.

Construction of the aging clock

The final aging clock was constructed as a Deep Neural Network (DNN) using TensorFlow. The model’s input consisted exclusively of the 32 selected hub genes. The harmonized dataset was partitioned into training (80%) and testing (20%) sets using stratified sampling. The model architecture comprised four hidden layers (512, 256, 128, and 64 neurons) with ReLU activation and dropout (rate = 0.3). It was trained using the Adam optimizer (learning rate = 0.0001) to minimize mean squared error, with an early stopping patience of 2000 to prevent overfitting.

External validation

To validate the model externally, 12 evenly distributed samples (aged 30–74 years, median age 52.5) were obtained from the E-MTAB-3037 dataset (ArrayExpress < BioStudies < EMBL-EBI) and used to estimate biological age. The prediction error was then calculated by comparing the estimated age with the actual donor age. Additionally, transcriptomic data from 10 HGPS patient samples in the GSE113957 dataset were input into the model to evaluate their biological age and quantify the extent of premature aging. Furthermore, we employed the GSE119009 dataset to assess the utility of our web application.

A Streamlit-based web application for model deployment

To improve the applicability of the predictive model in practical and experimental contexts, the finalized model was deployed as a web application built with the Streamlit Python framework (URL: https://skinage-2025.streamlit.app/) [32]. Upon uploading a gene expression matrix CSV file, the application automatically performs feature extraction, and the application provides the original data values, the data after batch effect correction, and the normalized data. Subsequently, it predicts physiological age based on these normalized data. To enhance interpretability and credibility, we evaluated the predicted values as aging scores rather than absolute chronological age.

Cell culture

H1 human embryonic stem cells (hESCs) were purchased from the WiCell Research Institute (Madison, WI, USA). The cells were maintained on Matrigel-coated 10-cm culture dishes in PGM1 medium. For routine passaging, cells were seeded at a density of 1 × 10⁶ cells per dish and were cultured for approximately 5 days to reach 85% confluency.

The HFF-1 cells, derived from human foreskin fibroblasts (HFFs), was obtained from the Cell Bank of the Chinese Academy of Sciences (Shanghai, China). These cells were cultured in T-25 culture flasks. For routine passaging, cells were seeded at a density of approximately 2 × 10⁵ cells per flask and cultured in high-glucose (4.5 mg/mL) Dulbecco’s Modified Eagle Medium (DMEM, Thermo Fisher Scientific, Gibco) supplemented with 15% fetal bovine serum (FBS, Thermo Fisher Scientific, Gibco) and 1% penicillin-streptomycin (Thermo Fisher Scientific, Gibco).

Isolation of extracellular vesicles

The cell culture of hESCs was maintained until at least 60% confluence. The collected culture supernatant was then purified by sequential centrifugation (300 × g for 10 min, then 2000 × g for 20 min) to discard the cells and larger debris. Following centrifugation, the resulting supernatant was filtered through a 0.22 μm sterile filter. Subsequently, EVs were isolated from the filtered supernatant using the exoEasy Maxi Kit (20) (QIAGEN, Hilden, Germany) in accordance with the manufacturer’s protocol.

Transmission Electron Microscopy (TEM)

For morphological characterization by TEM, EV samples were first fixed with 4% paraformaldehyde (PFA). A small aliquot of the fixed EV suspension was deposited onto 300-mesh copper grids. The grids then underwent sequential procedures: washing with phosphate-buffered Saline (PBS), negative staining with uranyl acetate, and rinsing with distilled water. After air-drying at room temperature for 2 h, the grids were examined under a transmission electron microscope to observe the morphology of hESC-EVs, and representative images were captured and saved.

Nanoparticle Tracking Analysis (NTA)

The particle size distribution and concentration of hESC-EVs were determined using a qNano system (Izon Science, Christchurch, New Zealand). Prior to sample analysis, the system was calibrated by establishing a standard curve with CPC100 calibration particles. Diluted EV samples were then introduced into the nanopore, and measurements for each sample were performed in triplicate. Data analysis was subsequently conducted based on the generated standard curve. To ensure the accuracy and repeatability of the measurements, standard operating procedures, including Nanopore conditioning (activation), baseline calibration, and PBS washing, were strictly followed throughout the process.

Induction of cellular photoaging

HFFs were first rinsed with PBS. Following this, a thin PBS film was placed over the cells to prevent dehydration during the irradiation process. The cells were then irradiated with a UVB dose of 300 mJ/cm² using a UVB lamp (311 nm, Philips, Netherlands) positioned 10 cm above the culture surface. Following irradiation, the HFFs were incubated for 48 h in fresh culture medium supplemented with or without hESC-EVs at a resulting concentration of 20 µg/mL.

Quantitative Real-time Polymerase Chain Reaction (qRT-PCR)

Total RNA was isolated from cells using the RNAex PRO RNA Isolation Reagent (Accurate Biology, Changsha, China) according to the manufacturer’s instructions. The RNA quantity and purity were evaluated by measuring the absorbance at 260 and 280 nm with a spectrophotometer. Complementary DNA (cDNA) synthesis was performed using the Evo M-MLV RT Kit (Accurate Biology). Quantitative real-time PCR was carried out with SYBR Green qPCR Master Mix (Accurate Biology). The primer sequences are provided in Table S1. Relative gene expression was determined using the 2^−ΔΔCt method and normalized to that of GAPDH.

Western Blot (WB)

Proteins were extracted from cells using RIPA buffer (Solarbio, Beijing, China) supplemented with PMSF (Beyotime, Shanghai, China) at a 100:1 ratio. Protein concentrations were quantified using a BCA assay kit (New Cell & Molecular Biotech, Suzhou, China). An equal amount of protein per sample was separated by 10% SDS-PAGE and subsequently transferred onto a PVDF membrane. The membranes were blocked with 5% non-fat milk for 1 h and then incubated overnight at 4 °C with the following primary antibodies: anti-p21 (Proteintech, Wuhan, China, #10355-1-AP), anti-p53 (Proteintech, Wuhan, China, #10442-1-AP), and anti-GAPDH (Proteintech, Wuhan, China, #60004-1-Ig) as a loading control. The protein bands were visualized via enhanced chemiluminescence (ECL) detection (New Cell & Molecular Biotech).

Immunofluorescence staining

The hESC-EVs (50µL) mixed with DiD (0.3µL) were ice-incubated in dark for 20 min, terminated with 1mL serum-free medium. After 100,000 g centrifugation for 90 min, pellet was resuspended in PBS, recentrifuged, resuspended in basic medium. Cells treated at optimal concentration, incubated 12 h, then PBS-washed 3x, fixed with 4% PFA in dark for 20 min, permeabilized with 0.2% Triton X-100 for 5 min, washed again. Phalloidin (1:200) added for 20 min shaking incubation in dark. After final PBS washes, cells on coverslips mounted with DAPI-containing anti-fade medium, observed via confocal microscopy.

Senescence-associated β-galactosidase (SA-β-Gal) staining

The cells were washed with PBS and fixed in fixative solution for 15 min at room temperature, followed by three rinses with PBS. The samples were then incubated with freshly prepared SA-β-Gal staining solution (Solarbio, Beijing, China) at 37 °C for 48 h. Images were captured using an inverted microscope (Zeiss, Germany), and the proportion of SA-β-Gal positive cells was quantified relative to the total number of cells.

Scratch-wound assay for cell migration and proliferation

To assess cell migration and proliferation, a scratch-wound assay was performed. HFFs were seeded into 6-well plates at a density of 5 × 10⁵ cells per well and cultured at 37 °C until a confluent monolayer was formed. A linear “scratch” was then created in the center of the cell monolayer using a sterile 200 µL pipette tip. The plates were washed with PBS to remove any detached cells and debris. Images of the scratch area were captured at 0, 12, 24, and 48 h post-scratch using a microscope. The area of the cell-free gap was measured at each time point using ImageJ software. The migration rate was quantified as the percentage of scratch closure, calculated using the formula: [(Area at 0 h - Area at time X) / Area at 0 h] * 100%.

Cell proliferation assay

The proliferation of HFFs was assessed via a Cell Counting Kit-8 (CCK-8, New Cell & Molecular Biotech, Suzhou, China), in accordance with the manufacturer’s guidelines. HFFs were initially seeded into 3.5 cm cell culture dishes and subsequently subjected to UVB-induced photoaging, as detailed previously. At 24 h postirradiation, CCK-8 reagent was added to the culture medium of each dish, and the cells were incubated for 1 h at 37 °C. The resulting absorbance at 450 nm was then recorded using a SpectraMax iD3 microplate reader (Molecular Devices, San Jose, CA, USA). This CCK-8 measurement procedure was repeated 24 h later (48 h postirradiation). The cell proliferation rate was calculated based on the absorbance values, and the results were plotted.

Statistical analysis

Each experiment in this study was conducted independently a minimum of three times. The data are expressed as the means ± standard deviations (SDs). Graphical presentation and statistical analyses were performed using GraphPad Prism software (version 8.0.2, GraphPad Software, San Diego, CA, USA). Comparisons between groups were made using one-way or two-way analysis of variance (ANOVA), followed by an appropriate post-hoc test. A p-value less than 0.05 was considered statistically significant. Significance levels are indicated as follows: ns, not significant; *p < 0.05; **p < 0.01; ***p < 0.001; and ****p < 0.0001.

Results

Data integration and identification of DEGs

Four datasets were integrated to form a consolidated dataset comprising 300 samples (Fig. 1A-C). Differential expression analysis of this dataset identified 185 DEGs (Fig. 1D). Separately, a list of 1556 skin-related genes was retrieved from the HAGR database. The intersection between the identified DEGs and this list of skin-related genes revealed 8 overlapping genes (Fig. 1E and Table S2).

Fig. 1
Fig. 1
Full size image

Data integration, batch effect removal, and identification of DEGs. A Principal component analysis (PCA) plots before and after batch effect correction. B Uniform manifold approximation and projection (UMAP) plots before and after batch effect removal. C Box plots depicting the data distribution before and after batch effect adjustment. D Volcano plot displaying DEGs in the integrated dataset following batch effect correction. E Venn diagram illustrating DEG-HAGR database gene overlap

Construction of WGCNA and identification of candidate hub genes via ML

WGCNA was employed to construct a scale-free co-expression network for pinpointing modules most pertinent to aging. After removing 8 outlier samples, a soft-thresholding power of 8 was selected based on the criteria of scale independence and mean connectivity to achieve a nearly scale-free topology (Fig. 2A-D). The analysis revealed that the chartreuse module exhibited the highest correlation with aging (r = 0.33, p = 6 × 10− 6) (Fig. 2E). As a result, a total of 20 genes exhibiting a Module Membership (MM) exceeding 0.7 and a gene significance (GS) above 0.3 were chosen from this chartreuse module for further analysis (Fig. 2F and Table S2).

Fig. 2
Fig. 2
Full size image

WGCNA and machine learning for candidate gene screening. A Dendrogram of gene clustering. B Assessment of the scale-free fit index for various values of soft-thresholding power (β). C Assessment of mean connectivity at various soft-thresholding power values (β). D Distinct gene co-expression modules are shown by various colors below the gene dendrogram. E Correlation heatmap between module eigengenes and aging status, emphasizing the chartreuse module’s significant association with skin aging. F Selection of 20 hub genes from the chartreuse module. G Cross-validation error for lambda selection. H Results from the SVM-RFE algorithm, which likely show feature ranking or cross-validation performance. I Venn diagram illustrating the intersection of genes identified by the LASSO and SVM-RFE machine learning algorithms

We subsequently employed LASSO regression and SVM machine learning algorithms to identify potential candidate genes associated with skin aging. LASSO regression analysis identified 70 genes highly correlated with the condition (Fig. 2G and Fig. S2). Subsequently, the SVM algorithm identified 87 genes (Fig. 2H). A Venn diagram depicting the overlap between the two gene sets revealed 15 shared genes (Fig. 2I and Table S2).

Development of an aging clock model

The genes identified through differential expression analysis, WGCNA, and machine learning approaches were combined, resulting in a set of 43 characteristic genes for model construction. From the initial 43 candidate genes, we performed a final quality control step. We excluded 11 genes that exhibited near-zero variance, specifically those with over 50% duplicated expression values across the samples after batch correction. This resulted in a final set of 32 robust hub genes for subsequent modeling. Nine different machine learning algorithms were then employed to develop age prediction regression models (Fig. 3A-B and Table S3).

Fig. 3
Fig. 3
Full size image

Performance comparison of nine different machine learning models for age prediction. A Fitting curves and performance evaluation metrics (R2, RMSE, MAE) for the nine machine learning models. B Comparison of performance metrics among the nine machine learning models

Out of all models evaluated, the DNN model consistently emerged as the top performer. It achieved high accuracy, evidenced by an R² value of 0.75, a root mean square error (RMSE) of 10.91, and a mean absolute error (MAE) of 8.26. These metrics collectively underscore its superior predictive capability within this comparative study.

SHAP-based feature interpretation and optimization of feature gene number

SHAP analysis was employed to interpret the DNN model. The SHAP summary plot and dependence plots (Fig. 4A-B and Fig. S3A-B) illustrate the importance of each feature (gene) in age prediction and assess the positive or negative impact of each feature on the model output (Fig. S4A-H). Subsequently, to optimize the number of features, genes were iteratively removed based on their SHAP values (from lowest to highest importance), the DNN model was retrained, and its predictive performance was assessed iteratively. This process revealed that a model using 32 features yielded the best performance, followed by models with 30 and 25 features (Fig. 4C). To mitigate the influence of chance, models constructed with these three feature set sizes (32, 30 and 25 features) were repeatedly trained and compared. The 32-feature model consistently demonstrated the most stable and optimal performance, and was therefore selected and saved for subsequent analyses.

Fig. 4
Fig. 4
Full size image

SHAP-based interpretation of the aging clock model and optimization of the feature gene number. A Ranking of the 32 gene features by importance according to SHAP values. B SHAP summary plot illustrating the impact of each of the 32 features on the model output for age prediction. C Performance evaluation (R2, RMSE, MAE) of the DNN model with varying numbers of input features. D Box plots showing the distribution of age prediction errors across different age groups (1–30 years, 31–60 years, and ≥61 years). E Box plots showing the distribution of age prediction errors for male and female subjects

To further evaluate the prediction error across different demographic groups, subgroup analyses were conducted based on age and sex. Samples were stratified into three age groups: 1–30 years, 31–60 years, and ≥ 61 years, as well as into male and female sex groups. Box plots of prediction errors indicated that younger individuals (1–30 years) exhibited the smallest prediction errors. While some middle-aged individuals (31–60 years) showed larger errors, there was a tendency for the predicted age to be underestimated in older individuals (≥ 61 years) (Fig. 4D). Error distributions were similar between the male and female samples, with some individuals in both groups displaying larger prediction errors (Fig. 4E).

External validation of the aging clock model and web application deployment

To assess the generalizability of the developed aging clock model, external validation was performed using two independent datasets. First, transcriptomic data from 12 samples with a relatively uniform age distribution (age range: 30–74 years; median age: 52.5 years) from dataset E-MTAB-3037 were processed and inputted into the model to estimate biological age. This validation yielded an MAE of 11.08 years (Table S4), indicating the reliability of our aging clock model on external data.

Furthermore, transcriptomic data from 10 HGPS patients from dataset GSE113957 were also inputted into the model. The predicted ages for these HGPS samples were substantially greater than their chronological donor ages (Table S5). This finding further supports the model’s capability to reflect accelerated biological aging and underscores its reliability in predicting physiological age.

To enhance the practical utility of our predictive model, we deployed it as a web-based application (Fig. 5). Users can upload a CSV file containing transcriptomic count data, upon which the application automatically extracts feature values, removes batch effects, performs normalization, and promptly outputs the sample identifiers along with their predicted biological aging scores upon clicking “Predict.” To improve interpretability and biological relevance, we evaluated model outputs as aging scores, which quantitatively reflect the transcriptomic senescence status of each sample rather than absolute chronological age.

Fig. 5
Fig. 5
Full size image

Web application interface for the aging clock model, facilitating experimental use. The user interface allows for the upload of a transcriptomic count matrix (CSV format). Upon upload, the application automatically extracts the expression data for the 32 feature genes, performs necessary preprocessing steps such as batch effect correction and data normalization, after the user clicks the "Predict" button, outputs the predicted biological age for each sample, identified by its original column name from the input file

To demonstrate the model’s usability and reliability, we further tested it on the GSE119009 dataset, which includes both normal and UVB-induced photoaged skin samples. The average aging score of the normal group was 42.698, whereas that of the photoaged group increased to 53.353 (Fig. 5), highlighting the model’s effectiveness in distinguishing aging phenotypes and its ease of application in experimental settings.

hESC-EVs mitigated senescence in UV-induced HFFs

To translate the predictions generated by the SkinAGE model into functional validation, we next investigated whether hESC-EVs could ameliorate the UVB-induced aging phenotype at the cellular and molecular levels. This transition from in silico prediction to in vitro experimentation enables testing of whether interventions that improve predicted biological age also impact age-associated cellular functions. A comprehensive characterization of hESC-EVs was conducted using TEM, NTA, and Western blot to assess their structure, particle size, and surface markers. TEM imaging revealed that the hESC-EVs exhibited the classic cup-shaped morphology (Fig. 6A). According to NTA, the size distribution of the EVs is approximately 80–150 nm in diameter (Fig. 6B). Western blot analysis confirmed that the EVs were enriched with characteristic exosomal marker proteins, including CD63 and CD9, while being negative for the Golgi apparatus marker GM130 (Fig. 6C) [19].

UVB radiation drives skin aging, largely through its effects on dermal fibroblasts which regulate the extracellular matrix (ECM) [33]. To establish the photoaging model, HFF viability was assessed after exposure to various UVB doses, and proliferation significantly decreased above 100 mJ/cm² (Fig. 6D). We used HFFs to study the protective effects of hESC-EVs against UVB-induced photoaging. First, hESC-EVs labeled with DiD were efficiently taken up by HFFs, and this uptake was unaffected by UVB irradiation (Fig. 6E). A dose of 300 mJ/cm² was chosen for subsequent experiments. The optimal EV concentration for promoting the proliferation of UVB-irradiated HFFs was found to be 20 µg/mL, which was used for all further treatments (Fig. 6F). The therapeutic effects of hESC-EVs (20 µg/mL) on UVB-irradiated HFFs were then investigated. EV treatment effectively mitigated the UVB-induced accumulation of SA-β-gal-positive senescent cells (Fig. 6G-H). Cell migration, which is crucial for skin repair [34], was impaired by UVB but significantly enhanced by EV treatment (Fig. 6I-J). Concurrently, it also attenuated the elevated expression of the senescence markers p21 and p53 (Fig. 6K-L). At the gene level, EV treatment tended to attenuate the UVB-induced upregulation of CDKN2A (p16), CDKN1A (p21), and TP53 (p53) mRNA (Fig. 6M). Furthermore, EVs markedly ameliorated the UVB-induced increases in MMP3 and IL-6 levels (Fig. 6N-O) and mitigated the UVB-induced reduction in collagen secretion (Fig. 6P).

Fig. 6
Fig. 6
Full size image

hESC-EVs attenuate photoaging in HFFs. A A representative TEM micrograph depicting the morphology of hESC-EVs. Scale bar: 100 nm. B Size distribution profile of hESC-EVs as determined by NTA. C Exosomal marker protein (CD63, CD9) expression by Western blot and a cell-specific marker (GM130, Golgi apparatus) in hESC-EVs lysates and cell lysates. D The influence of varying UVB doses on HFF proliferation was measured by CCK-8 assay. E Representative fluorescence microscopy images showing the uptake of DiD-labeled hESC-EVs by UVB-irradiated HFFs. Scale bar: 50 µm. F Effects of different concentrations of hESC-EVs on the viability of UVB (300 mJ/cm²)-irradiated HFFs, as assessed by CCK-8 assay. G Representative images of SA-β-gal staining in HFFs under different treatment conditions (NC, UVB, and UVB+EV). Scale bar: 200 µm. H Percentage of SA-β-gal positive HFFs. I Representative images of the scratch assay in HFFs at the indicated time points post-scratch under different treatment conditions. Scale bar: 400 µm. J HFF migration rate in scratch assays. K Representative Western blot analysis showing the protein levels of p21 and p53 in HFFs under different treatment conditions. [GAPDH is typically shown for normalization of protein loading]. L Western blot quantification of p21 and p53 protein levels. M Relative mRNA expression levels of CDKN2A (p16), CDKN1A (p21), and TP53 (p53) in HFFs, quantified by qRT-PCR. N Relative mRNA expression levels of IL-6 secreted by HFFs, quantified by qRT-PCR. O Relative mRNA expression levels of MMP3 in HFFs, quantified by qRT-PCR. P Relative mRNA expression levels of COL3A1 in HFFs, quantified by qRT-PCR. *p<0.05; **p<0.01; ***p<0.001; ****p<0.0001. Data are presented as mean ± SD (n=3)

Collectively, these results demonstrate that hESC-EVs can effectively counteract UVB-induced cellular senescence, impaired proliferation and migration, and pro-inflammatory responses in HFFs.

Transcriptomic analysis of the hESC-EVs effects on photoaged HFFs and integration with the aging clock model

To further explore the mechanistic link between SkinAGE predictions and the biological effects of hESC-EVs, we examined the expression patterns of key model-derived aging genes in UVB-treated versus hESC-EV-treated cells: RNA sequencing (RNA-seq) analysis was performed on three groups: a normal control group (NC), a UVB-irradiated group (UVB), and a UVB-irradiated group treated with hESC-EVs (UVB + EV), to systematically investigate the transcriptomic alterations underlying UVB-induced skin photoaging and EV-mediated therapeutic effects. Differential expression analysis revealed that UVB treatment induced significant upregulation or downregulation of a large number of genes, including key genes such as IGFBP1, MMP1, MMP3, IL1B, and CDKN1A. EV intervention partially ameliorated this altered expression pattern (Fig. 7A and Fig. S1). Gene Ontology (GO) functional enrichment analysis of genes differentially expressed in both comparisons (UVB vs. NC and EV vs. UVB) indicated significant enrichment in biological processes such as cell migration, cell population proliferation, extracellular matrix organization, and signal transduction. Notably, many genes associated with these pathways were downregulated following UVB irradiation and subsequently upregulated after EV treatment, suggesting that EV administration significantly rectified the expression of these dysregulated pathways (Fig. 7B and Fig. S5A-B). Analysis of gene expression heatmaps and trend clustering revealed distinct patterns. Genes exhibiting an initial upregulation by UVB followed by downregulation with EV treatment (up-down cluster) were primarily enriched in pathways such as the estrogen signaling pathway, B-cell receptor signaling pathway, and extracellular matrix-related signaling pathways. Conversely, KEGG pathway analysis revealed that the “down-up” gene cluster, characterized by initial downregulation in response to UVB and subsequent upregulation after EV treatment, was primarily associated with pathways involved in cell cycle inhibition, the p53 signaling pathway, cellular senescence, and inflammatory responses (Fig. 7C and Fig. S6A-B). GSEA revealed significant activation of the p53 signaling pathway in the UVB-induced photoaging group compared to the NC group. In line with this, the expression of key inhibitory factors such as CDKN1A (p21) and TP53 (p53) was markedly upregulated, indicating that UVB exposure induces cellular senescence (Fig. 7D). In contrast, GSEA further demonstrated significant activation of the cell cycle pathway in the EV-treated group. Key regulatory factors within this pathway, including CDK2, CCNB1, CCNB2, and CCNA2, were significantly upregulated following EV intervention (Fig. 7E), suggesting that EVs may promote cell cycle progression, thereby ameliorating UVB-induced cellular dysfunction and reducing proliferative capacity.

Fig. 7
Fig. 7
Full size image

Transcriptomic analysis of the effects of hESC-EVs on UVB-irradiated HFFs and biological age prediction by the SkinAGE. A Volcano plots displaying DEGs in the comparisons of UVB vs. NC and UVB+EV vs. UVB groups. B Ridge plots illustrating Gene Ontology (GO) enrichment analysis results for DEGs identified in the two comparison groups (UVB vs. NC and UVB+EV vs. UVB). C Heatmap of clustered gene expression profiles across the NC, UVB, and UVB+EV groups, with trend analysis and associated Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment for distinct gene clusters. D GSEA plot showing enrichment of the p53 pathway in the UVB group compared with the NC group. E GSEA plot showing enrichment of the cell cycle pathway in the UVB+EV group compared with the UVB group. F Bar chart displaying the predicted biological ages of HFFs in the NC, UVB, and UVB+EV groups, as determined by the skin aging clock model

Subsequently, the transcriptomic data from the NC, UVB, and UVB + EV groups were inputted into our previously developed web-based skin aging clock model for aging score assessment. The predicted aging scores for the three replicates in the NC group were 43.718, 40.605, and 47.858, with an average of 44.1. For the UVB group, the predicted scores were 75.855, 67.119, and 61.365, averaging 68.1. In contrast, the UVB + EV group exhibited predicted scores of 59.245, 44.391, and 37.059, with an average of 46.898. Compared to the NC group, UVB irradiation increased the average aging score of the cells by approximately 24 units. Remarkably, hESC-EV treatment reduced the aging score by an average of 21.22 units relative to the UVB group (Fig. 7F).

These SkinAGE score predictions were highly consistent with the trends observed in molecular experiments, further substantiating the therapeutic potential of hESC-EVs in combating skin aging. Moreover, these findings validate the effectiveness and applicability of our novel transcriptome-based aging clock in evaluating intervention outcomes, offering a promising assessment tool and conceptual framework for regenerative medicine.

Discussion

The aging process is intimately associated with the development of a multitude of chronic diseases. As the global aging population is rapidly increasing, effectively managing age-related diseases presents a significant hurdle for society and the healthcare system [35]. Therefore, discovering reliable biomarkers for aging is a critical endeavor. Advances in high-throughput technologies have propelled the extensive use of multi-omics approaches—such as genomics, exome sequencing, transcriptomics, and epigenomics—in diagnostics, understanding disease mechanisms, guiding treatment strategies, and predicting prognoses [36]. The vast datasets produced by these omics methodologies offer substantial potential for accurate age prediction [37, 38].

This study presents a skin aging clock built upon an DNN, a cutting-edge deep learning approach, designed to predict skin sample aging status with high precision. This model was subsequently deployed as a web application to enhance its practical utility. DNNs emulate the architecture of biological neural networks, leveraging multi-layer non-linear mappings to discern and learn complex relationships between gene expression patterns and age [39]. Compared with traditional statistical models or shallow machine learning methods, DNNs possess superior feature learning capabilities, enabling them to extracting important features from high-dimensional data, thereby improving prediction accuracy and generalizability. Furthermore, to enhance model interpretability and applicability, we significantly optimized the number of input features. This resulted in a more lightweight model that maintains high accuracy, making it suitable for a broader range of research and potential clinical applications. The robustness of this assessment system was validated through rigorous testing, including independent evaluation of external datasets, which demonstrated strong predictive performance and generalizability.

The application of deep learning in skin aging research and precision medicine has garnered considerable attention in recent years [40]. Various studies have employed machine learning and deep neural networks to analyze imaging [40], genomic [41], and proteomic data related to skin aging, aiming to construct more accurate predictive models [42]. For instance, convolutional neural networks (CNNs) have been utilized to assess aging features from skin images [43]. However, many existing studies are often constrained by limited dataset sizes, suboptimal feature selection, and lower model interpretability [44]. Our research addresses some of these limitations by integrating multiple skin transcriptomic datasets and employing robust feature selection methods to optimize the DNN model. This approach has led to the development of a more stable and user-friendly skin aging assessment system, the efficacy of which was further corroborated through external dataset validation.

Our investigation into the therapeutic potential of stem cell-derived extracellular vesicles against UVB-induced photoaging yielded compelling evidence of their efficacy. The observed amelioration of key senescence-associated molecular markers, coupled with enhanced cell viability and a marked reduction in senescence-associated β-galactosidase accumulation, strongly suggests that these EVs can effectively counteract the cellular damage and senescence pathways triggered by UVB irradiation [45, 46]. These findings align with a growing body of literature highlighting the regenerative and anti-senescent properties of EVs derived from various stem cell sources in the context of skin aging and other age-related pathologies [47]. A particularly noteworthy aspect of our study was the innovative application of our pretrained DNN-based aging clock to further quantify and validate these anti-aging effects. The significant decrease in the predicted aging score for EV-treated HFFs, when compared to UVB-exposed controls, provides robust, systems-level confirmation of the EVs’ rejuvenating capacity. This outcome serves a dual purpose: it not only reinforces the therapeutic promise of the specific EVs investigated but also powerfully demonstrates the sensitivity and practical utility of our aging clock as a tool for assessing the impact of anti-aging interventions [48]. The ability of the clock to discern such treatment effects underscores its potential for future applications in screening novel therapeutic compounds and personalizing anti-aging strategies. This dual validation—confirming EV efficacy while simultaneously showcasing the clock’s evaluative power—represents a significant contribution to both the fields of EV-based therapeutics and aging biomarker development.

The performance of our model is comparable to that reported by Fleischer [6], yet our approach utilizes a more extensive range of datasets and a feature set that is two orders of magnitude smaller. This substantial reduction in feature dimensionality greatly reduces the computational burden, thereby enhancing the system’s usability. However, several limitations should be acknowledged. Batch effects arising from diverse data sources could potentially impact model stability and generalizability [49]. Although we applied within-pipeline normalization, inter-dataset biases may still remain. Future work could explore more advanced batch correction techniques, such as Harmony [50] or deep learning-based dimensionality reduction [51], to improve data comparability. Another consideration is the direct merging of raw count values during data integration. We acknowledge that applying batch correction and machine learning models directly to raw, un-transformed RNA-seq counts is a non-standard approach, as log-transformation (e.g., log2(TPM + 1)) is typically recommended to stabilize variance. However, in our specific predictive modeling context, we conducted empirical comparisons. We observed that the model trained on Min-Max scaled raw counts consistently yielded superior predictive performance compared to models trained on log-transformed TPM or FPKM values. This finding suggests that for our deep learning architecture and specific set of gene features, the raw count distribution, once scaled, may have retained more discriminative information relevant to age prediction. This suggests that raw counts, despite potential scale issues between datasets if not handled carefully, may better reflect nuanced gene expression changes critical for the model. While within-pipeline normalization was applied, some inter-dataset biases might remain. Intriguingly, while the 32 genes prioritized in this study do not have extensive literature directly associating them with aging, the model’s exceptional efficacy can be attributed to the powerful capacity of the DNN to detect subtle shifts in gene expression and their underlying relationships.

While our findings highlight the rejuvenating potential of hESC-EVs in photoaged skin, cellular senescence also functions as a critical tumor-suppressive mechanism by preventing the proliferation of damaged cells. Thus, any intervention that mitigates senescence must be evaluated with caution. Our qPCR results support this balanced perspective: EVs caused only a slight reduction in P16 and P21 expression in normal fibroblasts, whereas the suppression was markedly stronger in UVB-induced photoaged cells (Fig. S7A-B). This pattern suggests that hESC-EVs preferentially alleviate stress-induced senescence without substantially altering basal senescence programs. Nevertheless, further mechanistic and long-term studies are required to fully assess the safety and implications of modulating senescence pathways with hESC-EVs.

Despite these limitations, our study makes important contributions to both the assessment of skin aging and the preclinical evaluation of novel anti-aging therapies such as stem cell-derived extracellular vesicles. We developed a simplified yet high-performing transcriptomic aging clock that not only accurately predicts biological age but also serves as a practical tool for evaluating therapeutic interventions. This dual functionality highlights its potential for future translational research. Looking ahead, we plan to refine our data integration pipelines, adopt more advanced batch correction strategies, and broaden the model’s application across different biological systems. These improvements will further enhance the accuracy and real-world utility of our aging assessment framework, supporting more effective and individualized approaches to skin rejuvenation.

Conclusion

This study presents a robust and interpretable transcriptome-based aging clock powered by an artificial neural network, capable of accurately predicting the biological age of skin-derived fibroblasts. Beyond its predictive performance, the SkinAGE provides a scalable tool for quantitatively assessing anti-aging interventions such as hESC-EVs. By integrating this computational framework with transcriptomic data, we demonstrate the biological relevance of our predictions. Our results also underscore the translational value of this framework for regenerative medicine and personalized treatments of age-related skin disorders.