Transcriptomic profiling of chlorogenic acid and taurine treatment in human skin cells provides insights into cellular senescence mechanisms

Three primary human skin cell lines were used. Normal Human Epidermal Keratinocytes (Cat. no. PCS-200–010; ATCC, Manassas, VA, United States) were cultured in Keratinocyte Growth Medium (Lonza, KGM™ Gold, Cat. no. 00192060; BS, Switzerland) supplemented with 10% fetal bovine serum (Gibco, Waltham, MA, United States), 100 U/mL penicillin, and 100 μg/mL streptomycin (Gibco). Neonatal light-pigmented Human Epidermal Melanocytes (HEMn-LP, Cat. no. C0025C; Thermo Fisher Scientific, Waltham, MA, United States) were cultured in Cascade Biologics™ Medium 254 supplemented with a Human Melanocyte Growth Supplement (Gibco, S0025). Human Dermal Fibroblasts were cultured in Dulbecco’s modified Eagle’s medium (DMEM, Sol Bio Pharm, Gyeonggi-do, Korea) supplemented with 10% fetal bovine serum (Gibco), 100 U/mL penicillin, and 100 μg/mL streptomycin (Gibco). Cells were maintained in a humidified incubator at 37 °C with 5% CO₂. For the assay, Normal Human Epidermal Keratinocytes and HEMn-LP were seeded at 2 × 10⁵ cells/well, and Human Dermal Fibroblasts at 1 × 10⁵ cells/well in 6-well plates. The cells were maintained for 10 h in a humidified incubator at 37 °C with 5% CO₂. Subsequently, CGA (Arshine Pharmaceutical Co., Ltd., Changsha, China) and taurine (Qianjiang Yongan Pharmaceutical, Co., Ltd., Qianjiang, China) were administered individually or in combination at appropriate concentrations (CGA, 10 μg/ml; taurine, 1,000 μg/ml; combined treatment, CGA 10 μg/ml + taurine 1,000 μg/ml), followed by incubation for 24 h under the same conditions. The selected concentrations were based on our previous findings (Lee et al., 2024a), which demonstrated that these doses exhibited the most pronounced cumulative effects on skin aging under co-treatment. The treatment time in this study was set to 24 h to ensure comparability with prior literature and to detect the integrated transcriptional effects of each compound (Moghadam et al., 2017; Subramanian et al., 2017; Alves et al., 2019; Lee et al., 2024a). After 24 h, the culture medium was removed and 1 mL of RNAlater (Cat. no. AM7020; Thermo Fisher Scientific) was added to each well to preserve RNA integrity. The plates were then immediately stored at −80 °C in a deep freezer for high-quality RNA preparation.

Total RNA was extracted using the TRIzol reagent (Thermo Fisher Scientific), QIAzol® Lysis Reagent (Qiagen, Germany), and RNeasy® Mini Kit (Qiagen), according to the manufacturer’s instructions. The total RNA concentration was measured using the Quant-iT^TM RiboGreen RNA Assay (Thermo Fisher Scientific). Total RNA integrity was assessed using a TapeStation RNA ScreenTape (Agilent Technologies, CA, United States). Samples with RNA integrity number >7.0 were used for RNA library construction. A library was independently prepared using 0.5 μg of total RNA for each sample by Illumina TruSeq Stranded Total RNA Library Prep Gold Kit (Illumina, San Diego, CA, United States) following the instructions in the Illumina TruSeq Stranded Total RNA Reference Guide. The libraries were quantified using the KAPA Library Quantification Kit for Illumina Sequencing platforms according to the qPCR Quantification Protocol Guide (KAPA BIOSYSTEMS, MA, United States) and TapeStation D1000 ScreenTape (Agilent Technologies). Total RNA sequencing (RNA-seq) was conducted by Macrogen (Seoul, Korea) using the NovaSeq X platform with 2 × 100 bp paired-end read chemistry (Illumina) (GEO ID: GSE302932↗).

The nf-core pipeline (v.3.17.0) (Ewels et al., 2020) was used for the alignment, quantification, and quality control of the raw data. RNA-seq reads were aligned to the GRCh38 reference genome obtained from the Broad Institute (Consortium, 2020) using STAR (v.2.7.11b) (Dobin et al., 2013) after filtering out alternate loci (ALT), human leukocyte antigen (HLA), and decoy sequence (Decoy) contigs. Isoform expressions of known Ensembl transcripts were quantified using Salmon (v.1.10.3) (Patro et al., 2017) and GENCODE release 47 (Mudge et al., 2025). QC and generating read counts were performed using the nf-core/rnaseq pipeline with RSeQC, Preseq, Qualimap, dupRadar, DESeq2, Kranken2, and MultiQC (Wang et al., 2012; Daley and Smith, 2014; Love et al., 2014; Ewels et al., 2016; Okonechnikov et al., 2016; Sayols et al., 2016; Wood et al., 2019).

To obtain reliable results, we applied strict QC criteria to the RNA-seq data instead of relying on independent filtering during differential expression tests using DESeq2 (v.1.44.0) (Love et al., 2014). Of the 78,816 generated genes, 60,801 remained after excluding spike-in controls, duplicates, artificial regions, unconfirmed genes, and pseudogenes. For each cell type, genes with read counts less than ten in at least one sample were excluded, resulting in 14,357, 14,445, and 14,528 genes retained in the epidermal keratinocytes, melanocytes, and fibroblasts, respectively. For each treatment group (treated and control) within each cell type, differential expression analysis was performed by comparing each individual sample against the remaining samples within the same group. The Bayesian shrinkage estimator for log₂ fold change (log₂FC), derived from the approximate posterior estimation of generalized linear model coefficients of each DEG, was used as the log₂FC value for all analyses in this study. Genes that passed the Benjamini–Hochberg multiple testing correction (adjusted P-value <0.05) and showed differential expression exceeding the suggestive threshold (|log₂FC| >0.585) within the same group were defined as within-group DEGs. For each cell type, genes identified as within-group DEGs in at least one treatment group were excluded to reduce heterogeneity within the same condition. A total of 14,111, 13,986, and 14,436 genes from epidermal keratinocytes, melanocytes, and fibroblasts, respectively, were included in the subsequent analyses. Principal component analysis (PCA) of the variance-stabilized gene expression data was performed using DESeq2 (v.1.44.0).

To identify genes responsive to CGA, taurine, and their combined treatment (CGA + Tau), samples treated with each compound were compared to untreated controls using DESeq2 (v.1.44.0) for each skin cell type, without applying independent filtering. Genes that passed the Benjamini–Hochberg multiple testing correction (adjusted P-value <0.05) and showed differential expression (|log₂FC| >1) were defined as DEGs. Log₂FC values were estimated using the apeglm shrinkage estimation (Zhu et al., 2019).

To test for synergistic effects, we defined two binary indicator variables representing CGA and taurine exposure: C=1 for CGA or CGA + Tau treatment (0 otherwise), and T=1 for taurine or CGA + Tau treatment (0 otherwise). For each gene i and sample j within each cell type, a negative binomial generalized linear model was fitted using DESeq2 (v.1.44.0), logqij=βi0+βiCCj+βiTTj+βiCTCjTj, where qij denotes the normalized mean expression. The interaction coefficient βiCT corresponds to the deviation of the combined treatment from additivity: βiCTlog⁡2=log2⁡FCiCGA+Tau−log2⁡FCiCGA+log2⁡FCiTau, where βiCTlog⁡2 denotes log₂FC of the interaction term. For each cell type, we tested the null hypothesis H0:βiCT=0 using a two-sided Wald test. Genes were considered to exhibit potential synergistic effects of CGA + Tau if the interaction term reached nominal significance (P_interaction < 0.05) and satisfied a synergy criterion (log2⁡FCiCGA+Tau>log2⁡FCiCGA+log2⁡FCiTau, with concordant effect directions). Given the limited statistical power of interaction tests (McClelland and Judd, 1993; Leon and Heo, 2009), genes reaching a less stringent multiple-testing correction threshold (Benjamini–Hochberg adjusted P_interaction <0.1) were considered to exhibit significant synergistic effects.

To provide additional evidence for the identified DEGs, we utilized the Connectivity Map (CMap) database, which offers transcriptomic profiles of human cell lines treated with various perturbations (Subramanian et al., 2017). Using cmapR v.1.16.0 (Enache et al., 2019), we extracted level 5 L1000 signatures, which consist of moderated Z-score vectors (ModZ) as differential gene expression vectors, for each CGA treatment at doses of 10, 3.33, 1.11, 0.37, 0.125, and 0.04 µM in the melanocyte-derived human skin cancer cell line A375. Genes were considered to have CMap supporting evidence if they satisfied either of the following two criteria with consistent directionality in their expression changes: (1) |log₂FC| >1 and |ModZ| >1.67; or (2) |log₂FC| >0.585 and |ModZ| >2. Among the identified DEGs, those with significant adjusted P-value but modest fold change (0.585 < |log₂FC| ≤1) were also included if supported by CMap evidence. These genes were included into downstream analyses.

Functional enrichment analyses of the identified DEGs with canonical pathways (KEGG, REACTOME, and Wiki pathways) and Gene Ontology (GO) terms (molecular function [GO:MF], cellular component [GO:CC], and biological process [GO:BP]) were performed using gprofiler2 (v.0.2.3) (Kolberg et al., 2020). Pathways and GO terms with gene sets that passed the false discovery rate corrected P-value threshold (adjusted P-value <0.05) were considered significantly enriched. Among the identified pathways and GO terms, those related to antioxidative, anti-inflammatory, and anti-senescence effects were manually categorized based on functional descriptions and grouped into broader categories according to biological relevance and shared terminology: cellular senescence and oxygen response, cell cycle regulation, extracellular matrix organization, and immune and oxidative stress regulation.

We defined DEGs as AR-DEGs if they were supported by any of the following evidence of relevance to aging: (1) prior annotation in aging-related databases such as Aging Atlas (Aging Atlas, 2021), Aging Map (Mao et al., 2023), or GenAge (de Magalhaes et al., 2024); (2) proximity (within 500 kb) to genetic variants reaching genome-wide significance (P-value <5 × 10⁻⁸) in previous genome-wide association studies (GWASs) of perceived age (Roberts et al., 2020; Ingold et al., 2024); (3) inclusion in aging-related pathways or GO terms among those significantly enriched in the functional enrichment analysis; or (4) inferred to interact with drugs with antioxidant or anti-inflammatory activity, as annotated in DGIdb v.5.0.4 (Cannon et al., 2024). The regulons within the AR-DEGs were inferred using transcription factor–target interactions with the highest confidence level A from DoRothEA (Garcia-Alonso et al., 2019).

Associations between AR-DEGs and skin aging-related traits, such as perceived age and skin color (CIE LAB values: L*, a*, and b*) were tested using FUSION (released on 2022–02–01) (Taylor et al., 2019), a tool that performs transcriptome-wide association study (TWAS) based on GWAS summary statistics by mapping genes to traits through expression quantitative trait loci (eQTLs). Summary statistics of GWAS for skin color (48,433 individuals of East Asian ancestry) reported by Kim B. et al. (2024) were obtained from the NHGRI-EBI GWAS Catalog (GCST90320257 for L*, GCST90320258 for a*, and GCST90320259 for b*). GWAS summary statistics for perceived age (European ancestry) were obtained from studies by Roberts et al. (2020) (423,992 individuals) and Ingold et al. (2024) (403,945 individuals), available from the University of Bristol data repository (https://data.bris.ac.uk/data/dataset/21crwsnj4xwjm2g4qi8chathha↗) and Zenodo (https://doi.org/10.5281/zenodo.10554253↗), respectively. Using cis-eQTLs for AR-DEGs (defined as variants within ±1 Mb of the transcription start site), associations between genes and traits were tested using precomputed gene expression weights from the Genotype-Tissue Expression project (GTEx) v8 (Consortium, 2020) for non-sun-exposed suprapubic and sun-exposed lower leg skin tissues. Genes that passed the Benjamini–Hochberg multiple testing correction (adjusted P-value <0.05) for the number of genes in each test group (tissue-trait pair) were considered statistically significant. Among these, gene–trait eQTL mappings were considered reliable when both the Z-score of the top cis-eQTL for the gene and the corresponding GWAS Z-score of that variant were greater than 3.

To quantify the mRNA expression levels, we conducted quantitative real-time PCR (RT-PCR) analysis. Cell culture conditions were the same as those described in Section 2.1. Total RNA was extracted using the AccuPrep® Universal RNA Extraction Kit (Bioneer, Daejeon, Republic of Korea) according to the manufacturer’s instructions. The purity of the extracted RNA (A260/A280) was assessed using the NanoDrop spectrophotometer. Complementary DNA (cDNA) was synthesized by reverse transcription using the AccuPower® RocketScript™ Cycle RT PreMix (Bioneer) on a PCR thermocycler (R&D Systems), according to the manufacturer’s protocol. Quantitative RT-PCR was performed using cDNA obtained from control cells and cells treated with CGA and taurine. The following TaqMan probes were used: GAPDH (Assay ID: 4333764F) as an internal control, TGFB2 (Hs00234244_m1), ETS1 (Hs00428293_m1), IL1A (Hs00899844_m1), and IL1B (Hs01555410_m1). The TaqMan™ Universal Master Mix II, with UNG (Applied Biosystems, Waltham, MA, United States) was used for amplification. PCR reactions were performed on the ABI 7500 Real-Time PCR system according to the manufacturer’s protocol. Data were analyzed using ABI software (version 2.3).

The expression levels of p16 and p21 in fibroblasts were evaluated by Western blotting using incubated supernatants and cell lysates. Cells were washed with ice cold PBS and lysed on ice in M-PER buffer (Thermo Fisher Scientific, MA, United States) supplemented with Complete™ protease inhibitor cocktail and phosphatase inhibitor (Roche, Indianapolis, IN, United States). 40 μg of protein was analyzed by Western blotting with appropriate antibodies. Protein bands were detected using a chemiluminescence detector (iB-right FL1500, Thermo Fisher Scientific, MA, United States). Western blotting results were quantified by measuring band intensity using ImageJ software (NIH, MD, United States) for three individual trials. The expression levels of p16 and p21 were normalized to β-actin.

A total of 36 bulk RNA-seq samples were generated from three human skin cell types (i.e., epidermal keratinocytes, melanocytes, and fibroblasts) each treated with one of three compounds (i.e., CGA, taurine, or CGA + Tau), or left untreated as a control. Each treatment condition included three biological replicates per cell type. PCA of overall gene expression profiles showed clear separation among the three skin cell types (). Supplementary Figure S1A

We independently performed differential gene expression analysis after treatment with each compound in each human skin cell type (;). In total, 14,111, 13,986, and 14,436 genes from epidermal keratinocytes, melanocytes, and fibroblasts, respectively, were tested to identify DEGs in response to treatment (). Genes with inconsistent expression among samples within the same cell type and treatment conditions (within-group DEGs; see Methods) were excluded prior to analysis, resulting in a more homogeneous transcriptome profile for downstream analysis (). Supplementary Figure S2 Supplementary Table S1 Supplementary Table S2 Supplementary Figure S1B, C

We identified 190 DEGs based on an adjusted P-value <0.05 and |log₂FC| >1 (Figure 2A; Supplementary Figure S3A; Supplementary Table S3). In addition, of the 857 genes, seven genes (E2F2, CDC42EP3, NRG1, CEBPD, ABHD4, TUBB6, and CXCL10) that met a suggestive threshold (adjusted P-value <0.05 and |log₂FC| >0.585) had supporting evidence from the CMap database (Subramanian et al., 2017) and were also considered DEGs responsive to CGA (Supplementary Figure S4; Supplementary Table S3). Among the 197 DEGs, 147, 41, and 10 were identified in epidermal keratinocytes, melanocytes, and fibroblasts, respectively, with all DEGs being specific to a single cell type, except for ANGPTL4. A total of 174, 86, and 16 DEGs were responsive to CGA + Tau, CGA, and taurine, respectively. Of note, five DEGs showed |log₂FC| greater than 2, including CCN2 and KRTAP2-3 (responsive to both CGA + Tau and CGA in epidermal keratinocytes), PDE3B and MARCHF4 (responsive to CGA + Tau in epidermal keratinocytes), and a long non-coding RNA gene, ENSG00000280800 (responsive to CGA in melanocytes). The majority of DEGs were responsive, either specifically to CGA + Tau (52.7%, 104 DEGs) or to both CGA + Tau and CGA (29.9%, 59 DEGs). Considering the CGA + Tau treatment alone, 88.3% (174 DEGs) of all identified DEGs were detected. In addition, IL1B in epidermal keratinocytes and AK4, ANGPTL4, BNIP3, PKD1, and TAC1 exhibited potential synergistic effects of CGA + Tau treatment (P_interaction <0.05 and satisfying a synergy criterion; see Methods) (Supplementary Table S4; Supplementary Figures S5 and S6). Among these DEGs, AK4 and ANGPTL4 passed multiple testing correction under a less stringent threshold (adjusted P_interaction = 0.051). The DEGs identified in this study showed directional concordance in expression changes across treatments with the three compounds within a given cell type, with CGA + Tau inducing greater fold changes than either compound alone (Supplementary Figure S7). However, when a given compound was used to treat different cell types, distinct sets of DEGs and their highly cell-type-specific fold changes were identified (Supplementary Figure S8).

We identified DEGs encoding transcription factors, such as ETS1, EGR1, and E2F2, and their high-confidence targets using DoRothEA (Garcia-Alonso et al., 2019) (Supplementary Table S5). Among these, ETS1 and EGR1 were identified as regulators of other DEGs at the highest confidence level (Figure 2B). EGR1 was also identified as a high-confidence target of ETS1, a DEG responsive in epidermal keratinocytes (CGA + Tau, log₂FC = −1.05, standard error [SE] = 0.173, adjusted P-value = 8.12 × 10⁻⁹), whereas EGR1 was regulated specifically in melanocytes (CGA + Tau, log₂FC = 1.73, SE = 0.407, adjusted P-value = 1.25 × 10⁻⁴; taurine, log₂FC = 1.52, SE = 0.417, adjusted P-value = 5.03 × 10⁻³). Other targets of ETS1, including BMP4, DUSP6, and THBS1, showed decreased expression in response to CGA and taurine treatments, which was consistent with the downregulation of ETS1 in epidermal keratinocytes. Among the targets of EGR1, ANGPTL4 showed increased expression in response to compound treatment in melanocytes (log₂FC = 1.34, SE = 0.256, adjusted P-value = 2.34 × 10⁻⁶), consistent with upregulation of EGR1, but was downregulated in epidermal keratinocytes (log₂FC = −1.15, SE = 0.174, adjusted P-value = 2.52 × 10⁻¹²). Other EGR1 targets, including THBS1 and FGF2, showed decreased expression in epidermal keratinocytes.

To identify the potential mechanisms by which DEGs contribute to antioxidative, anti-inflammatory, and anti-senescence effects, we performed functional enrichment analysis using canonical pathways (KEGG, REACTOME, and WikiPathways) and GO terms. Of the 197 identified DEGs, 126 protein-coding genes were included in the enrichment analysis, which identified 71 canonical pathways and 405 GO terms as significantly enriched gene sets (adjusted P-value <0.05) (Supplementary Table S6). Among the enriched pathways and terms of the DEGs, several functional categories relevant to skin aging were identified, such as oxygen response (GO:1901700, GO:1901701, and GO:0070482), cell cycle regulation (cellular senescence [KEGG:04218], MAPK cascade [GO:0000165], and TGF-beta signaling [KEGG:04350, WP:WP560, WP:WP366]), extracellular matrix organization (GO:0062023 and REAC:R-HSA-1474244), and immune and oxidative stress regulation (cytokine-cytokine receptor interaction [KEGG:04060], vitamin D receptor pathway [WP:WP2877], and NRF2 pathway [WP:WP2884]) (Figure 2C). Of note, 48 DEGs were involved in one or more of these skin aging-related functional categories (Supplementary Table S7). Among these, 23 were involved in two or more functional categories. The transcription factors ETS1 and EGR1, along with their targets (DUSP6, BMP4, THBS1, FGF2, and ANGPTL4), formed a regulon whose components were not part of a single pathway but were instead distributed across multiple pathways related to cellular longevity. TGFB2 was involved in all major anti-senescence-related functional categories and was intricately interconnected with other DEGs within each category (Figure 2D). TGFB2 might activate the cell cycle by influencing PDK1 (GO:0070482 and WP:WP366) and ETS1 (KEGG:04218 and WP:WP366), which are involved in cellular senescence and oxygen response pathways. Additionally, TGFB2 is a reported downstream target of THBS1 (Schultz-Cherry et al., 1994), a DEG involved in the extracellular matrix organization pathway (REAC:R-HSA-1474244 and WP:WP366). PDE3B and CCN2 were both markedly regulated by CGA + Tau treatment in epidermal keratinocytes. PDE3B was upregulated (log₂FC = 2.28, SE = 0.316, adjusted P-value = 3.57 × 10⁻¹⁴) and involved in the functional category of cellular response to oxygen-containing compound (GO:1901701). CCN2 was downregulated (log₂FC = −2.21, SE = 0.122, adjusted P-value = 1.92 × 10⁻⁷⁰) and involved in both the MAPK cascade (GO:0000165) and collagen-containing extracellular matrix (GO:0062023).

We prioritized 62 AR-DEGs, including 48 involved in skin aging-related functional categories and 14 additional genes, supported by evidence of aging relevance (Figure 3; Supplementary Figure S3B; Supplementary Table S8). For instance, among the DEGs, 16 were previously reported in aging databases such as the Aging Atlas (Aging Atlas, 2021), Aging Map (Mao et al., 2023), and GenAge (de Magalhaes et al., 2024). In addition, 13 were located within 500 kb of previously reported GWAS loci associated with perceived age (Roberts et al., 2020; Ingold et al., 2024). We inferred drugs interacting with antioxidant-responsive DEGs using the DGIdb (Cannon et al., 2024) and identified nine genes, CCN2, CYP24A1, FST, IL1A, IL1B, PDE3B, PDE4B, PMP22, and TGM2, that interact with 18 known antioxidative or anti-inflammatory agents, including apremilast, calcitriol, crisaborole, apremilast, isotretinoin, theophylline, retinoic acid (tretinoin), and vitamins A, D, and E (Supplementary Table S9). A heatmap showing the standardized expression levels of the 62 AR-DEGs across samples within each cell type is provided in Supplementary Figure S9. Among the 62 AR-DEGs, 56 were CGA + Tau-responsive and 52 were differentially expressed in epidermal keratinocytes. Of note, 37.1% (23 of 62 genes) of the AR-DEGs were identified specifically under the CGA + Tau treatment. In addition, five of the six DEGs exhibiting potential synergistic effects of CGA + Tau treatment were prioritized as AR-DEGs, including AK4 and ANGPTL, which passed multiple testing correction in melanocytes (Supplementary Figures S5 and S6). All these synergistic AR-DEGs were involved in oxygen response-related GO terms. Although AR-DEGs identified exclusively under single-treatment conditions were not classified as CGA + Tau-responsive AR-DEGs, they exhibited suggestive expression changes (|log₂FC| > 0.585) under co-treatment.

Associations between the predicted expression levels of AR-DEGs and skin aging-related traits (perceived age and skin color [CIE LAB values: L*, a*, and b*]) were tested using a targeted TWAS. A total of 52 AR-DEGs were tested using cis-eQTLs shared between GTEx and the GWAS summary statistics for the skin aging-related traits (Supplementary Table S10). In this targeted TWAS, four AR-DEGs were found to be associated with skin aging-related traits. BNC2 was associated with perceived age in suprapubic skin tissue (Z-score of TWAS [Z_TWAS] = −6.20, adjusted P-value = 2.91 × 10⁻⁸) (Supplementary Figure S10). CDC42EP3 (Z_TWAS = 3.45, adjusted P-value = 0.010), NT5E (Z_TWAS = 3.40, adjusted P-value = 0.01), and NPNT (Z_TWAS = −2.83, adjusted P-value = 0.042) were associated with perceived age in lower leg skin tissue, but were not considered to have reliable gene–trait associations through eQTLs, as they did not meet the criterion for both the top cis-eQTL Z-score and the corresponding GWAS Z-score (Z-score >3). ADM (suprapubic, Z_TWAS = −4.77, adjusted P-value = 4.70 × 10⁻⁵; lower leg, Z_TWAS = −5.54, adjusted P-value = 1.35 × 10⁻⁶) and MIR3936HG (suprapubic, Z_TWAS = −3.62, adjusted P-value = 4.98 × 10⁻³; lower leg, Z_TWAS = −3.01, adjusted P-value = 0.029) were associated with perceived age in both skin tissues (Supplementary Figures S11 and S12). FST, an AR-DEG downregulated by CGA + Tau treatment, was associated with decreased L* (brightness, Z_TWAS = −3.28, adjusted P-value = 0.043) and increased a* (redness, Z_TWAS = 3.70, adjusted P-value = 8.90 × 10⁻³) in lower leg skin tissue (Supplementary Figure S13).

The potential anti-senescence effects of CGA and taurine were additionally validated through in vitro assays, including quantitative RT-PCR analysis of representative AR-DEGs for cellular senescence (TGFB2, ETS1, IL1A, and IL1B) and Western blot analysis of p16 and p21 proteins. In epidermal keratinocytes, treatment with CGA or CGA + Tau significantly suppressed expression of TGFB2, a known pro-senescence factor (CGA, 0.52-fold, P-value = 5.96 × 10⁻³; CGA + Tau, 0.48-fold, P-value = 5.17 × 10⁻³) (Supplementary Figure S14A). This result was consistent with the transcriptomic profiling, which showed TGFB2 downregulation under the same conditions. In fibroblasts, changes were observed in other senescence markers, including cell cycle inhibitors and pro-inflammatory cytokines. Specifically, the mRNA level of IL1A was decreased following treatment with taurine or CGA + Tau (taurine, 0.68-fold, P-value = 0.046; CGA + Tau, 0.59-fold, P-value = 0.011) (Supplementary Figure S14B). In parallel, Western blot analysis demonstrated that CGA, taurine, or their combination reduced the protein levels of p16 and p21 in fibroblasts (Supplementary Figure S14C). The effects of CGA and Tau were most pronounced under combined treatment, suggesting a potential enhancement in attenuating cellular senescence at the molecular level. In contrast, ETS1 expression exhibited changes in the opposite direction at the mRNA level in epidermal keratinocytes (CGA + Tau, 5.29-fold, P-value = 0.011) and IL1B expression showed modest differences in fibroblasts (CGA + Tau, 0.63-fold, P-value = 0.092) following treatment (Supplementary Figures S14A and S14B).

The summary statistics of GWAS for perceived age in UK Biobank European participants can be downloaded at https://doi.org/10.5523/bris.21crwsnj4xwjm2g4qi8chathha↗ (Roberts et al., 2020) and https://zenodo.org/records/10554253↗ (Ingold et al., 2024). The summary statistics of GWAS for skin color in East Asian populations can be downloaded at the NHGRI-EBI GWAS Catalog (GCST90320257; https://www.ebi.ac.uk/gwas/studies/GCST90320257↗, GCST90320258; https://www.ebi.ac.uk/gwas/studies/GCST90320258↗, and GCST90320259; https://www.ebi.ac.uk/gwas/studies/GCST90320259↗) (Kim B. et al., 2024). The Connectivity Map (CMap) resources can be downloaded at https://clue.io↗. The accession numbers for the processed RNA-seq data from the skin samples used in this study are available at the National Center for Biotechnology Information/Gene Expression Omnibus under repository accession number Gene Expression Omnibus (GSE302932↗; https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE302932↗). United Kingdom Biobank data were obtained under application no.33002.