This study utilized SL skin samples that we previously collected from seven patients at Kyung Hee University Hospital at Gangdong, where all procedures were approved by the local Institutional Review Board (IRB). Written informed consent was obtained from all participants prior to sample collection, and the study was conducted in full accordance with the ethical principles for research involving human subjects. Paired skin samples (5 mm diameter) were collected from the SL and adjacent normal skin of seven patients at Kyung Hee University Hospital at Gangdong. The mean age was 82.0 ± 10.2 years (range, 67–96 years). The locations were on the chin (71.4%), cheek (14.3%), and temple (14.3%). Generation of raw sequencing data and subsequent bioinformatic processing were conducted according to our previously published methods (Choi et al., 2025). Key processing steps included read alignment to the human reference genome (GRCh38) and normalization of transcript expression to transcripts per million (TPM) for all downstream analyses.

Differential expression analysis between SL and control samples (n = 7 per group) was performed on the raw count matrix using the DESeq2 package (v1.48.0) in R. Genes were considered significantly differentially expressed if they met a False Discovery Rate (FDR) adjusted p-value <0.05 and an absolute log2 fold change >0.01.

Functional annotation of differentially expressed genes was conducted through Gene Ontology (GO) term enrichment using the clusterProfiler R package (v4.16.0) (Wu et al., 2021). To examine pathway-level changes, Gene Set Enrichment Analysis (GSEA) was applied using pre-ranked gene lists ordered by log2 fold change values, implemented through the fGSEA R package (v1.24.0) (Korotkevich et al., 2016). Metabolic and signaling pathway analysis was performed using KEGG database annotations via the gage package (v2.58.0). Statistical significance for all enrichment tests was set at FDR-adjusted p-value <0.05. Protein interaction networks were constructed by querying the STRING database (v12.0) (Szklarczyk et al., 2023) with the differentially expressed gene set. Network topology analysis and visualization were subsequently performed using Cytoscape (v3.10.3) (Shannon et al., 2003). Functional annotation enrichment analysis was conducted for each PPI network using Gene Ontology biological process terms to elucidate pathway-specific molecular functions.

To identify stable gene clusters based on co-expression patterns, consensus clustering was performed on the subset of genes related to the predefined metabolic pathways using the ConsensusClusterPlus (v1.64.0) package in R. The log2 fold change values for these genes were used as the input data matrix. Prior to clustering, the data was centered by the median of each gene. The k-means algorithm with Euclidean distance was applied for clustering. The analysis was repeated 30 times, and the number of clusters (k) was evaluated across a range from 2 to 10 to ensure the stability of the clustering results. The consensus matrix, generated from this analysis, served as a measure of co-expression similarity between all gene pairs.

Genes associated with energy metabolism were selected based on the model’s gene-reaction annotation framework, followed by differential expression analysis using the DESeq2 R package (v1.48.0). Expression pattern-based clustering of differentially expressed genes was executed using the ConsensusClusterPlus R package. The K-means clustering algorithm implemented in ConsensusClusterPlus facilitated the identification of distinct gene expression clusters and enabled the characterization of potential molecular subtypes within the dataset.

To simulate cellular metabolism alterations, we employed the custom, constraint-based flux simulation previously developed by our group using the Recon1 genome-scale model (da Silveira et al., 2020; Guarnieri et al., 2023; Duarte et al., 2007). The simulation was governed by two primary constraints: first, corresponding enzyme expression levels were used to set the bounds for reaction fluxes, and second, main energy-associated pathways were set to maximal optimization (specifically, Citric Acid Cycle, Oxidative Phosphorylation, CoA Synthesis, CoA Catabolism, Glycolysis and Gluconeogenesis, NAD Metabolism, Fatty Acid Synthesis, Fatty Acid Oxidation, and Biomass and Maintenance Functions). The Van der Waerden (VdW) test was then used to quantify the resulting metabolic alterations between the two conditions. Metabolic flux simulations utilized reaction identifiers from the BiGG Models database (http://bigg.ucsd.edu/).

In the context of SL, our analysis revealed considerable upregulation of critical phospholipid synthesis reactions. At the transcriptional level, gene co-expression network analysis (Figure 1A) revealed significant regulatory interconnections among key genes involved in glycerophospholipid metabolism. As shown in Figure 1A, these genes (highlighted in yellow), such as DGKQ, PTDSS1, and PLA2G7, anchored to distinct co-expression clusters that demonstrated strong intra-pathway co-regulation and significant inter-pathway co-expression. Major glycerophospholipid gene clusters centered around DGKQ (closely co-expressed with other glycerophospholipid enzymes such as PLD2) and PTDSS1 (co-expressed with PAFAH2) showed robust connections to genes involved in sphingolipid metabolism (DEGS1) and methionine metabolism. Additionally, PLA2G7 exhibited co-expression with genes integral to cholesterol metabolism (ACAT2) and fatty acid elongation (ELOVL). This intricate network view underscores extensive crosstalk with other critical lipid pathways and supports the observed upregulation of phospholipid metabolic processes in the SL.

Supporting these transcriptional findings, protein-protein interaction (PPI) analysis (Figure 1B; Supplementary Figure S1) demonstrated functional clustering among PLA2 family enzymes (EC:3.1.1.4), including PLA2G10, PLA2G2D, PLA2G3, PLA2G2A, PLA2G4C, PLA2G4A, and PLA2G5, which facilitate arachidonic acid release from phosphatidylcholine. These protein-protein interaction patterns reinforce the coordinated regulation of phospholipid metabolic machinery.

At the metabolic level, flux simulations of the glycerophospholipid pathway (Figure 1C) demonstrated these transcriptional changes translate into functional metabolic alterations. We observed a pronounced increase in Phosphatidylserine decarboxylase reaction (BiGG ID: PSDm_hs) activity, which converts phosphatidylserine (PS) into phosphatidylethanolamine (PE). PE serves as an essential phospholipid that contributes to membrane integrity and stabilization of mitochondrial proteins. Furthermore, substantial upregulation of the PETOHMr_hs reaction was detected, catalyzed by phosphatidylethanolamine N-methyltransferase. This process involves the methylation of phosphatidylethanolamine to produce phosphatidylcholine (PC), utilizing S-adenosylmethionine as a methyl donor and yielding S-adenosylhomocysteine as a byproduct. Given that phosphatidylcholine is a predominant membrane component vital for lipid homeostasis and signaling, its increased synthesis in the SL suggests concerted modifications in membrane phospholipid composition that align with the observed transcriptional regulatory networks.

As shown in Figure 2C, our metabolic flux simulation of fatty acid metabolism in the SL revealed a substantial downregulation of the stearoyl-CoA desaturase reaction (BiGG ID: DESAT18_3), suggesting impaired fatty acid desaturation and elongation processes that may alter membrane rigidity. This reaction converts stearoyl-CoA to oleoyl-CoA. Stearoyl-CoA desaturase is critical for fatty acid desaturation, an essential process for maintaining lipid composition and membrane fluidity (Ntambi, 1999). Another reaction, dihydroceramide desaturase reaction (BiGG ID: DHCRD2), of the sphingolipid pathway showed a predominant reduction in SL, where the reduction in the activity of dihydroceramide desaturase led to a decline in ceramide production (Figure 3C). Dihydroceramide desaturase plays a pivotal role in the enzymatic conversion of dihydroceramide to ceramide, which is crucial for maintaining the equilibrium between sphingolipids and dihydrosphingolipids. In contrast, the upregulation identified through PPI analyses can be attributed to the increased transcription and translation of dihydroceramide desaturase, which serves as a compensatory mechanism for diminished DHCRD2 reaction. As shown in Figure 3B (see also Supplementary Figure S2), our PPI analysis revealed the upregulation of DEGS1 and DEGS2, both annotated as EC:1.14.19.17 and EC:1.14.18.5, respectively, in the KEGG database. These enzymes likely function in compensatory mechanisms to mitigate the reduced activity of dihydroceramide desaturase, thereby preserving cellular homeostasis. Such compensatory responses are frequently observed in patients with dermatological disorders (Coelho et al., 2015).

The complex interplay within the lipid metabolism was further illustrated by gene network analyses, as shown in Figure 2A (fatty acids) and Figure 3A (sphingolipids). Figure 2A shows that ELOVL4, crucial for the synthesis of very long-chain fatty acids (VLCFAs) necessary for ceramide production, is connected to PLA2G7, a gene related to glycerophospholipid metabolism, and ACAT2, a gene related to cholesterol metabolism. Figure 3A depicts a broader network, particularly within sphingolipid metabolism, which includes interactions involving SMPD3 (involved in ceramide generation), DEGS1 (dihydroceramide desaturase, critical for synthesizing ceramide from dihydroceramide), and UGCG (which utilizes ceramide), along with their connections to genes involved in other lipid pathways. These network analysis results collectively highlight the integrated regulation of pathways essential for maintaining sphingolipid homeostasis and membrane integrity in SL.

Having recognized that oxidative stress disrupts the balance of the cellular redox state, affecting homeostasis and cellular responses in the SL, our analysis of sphingolipid metabolism revealed a predominant downregulation of the DHCRD2 reaction in the SL (Figure 3C). This reaction relies on FADH2 as a cofactor, and its enzymatic function is influenced by the presence of O₂ and NAD(P)H. These observations highlight the bidirectional interactions between the cellular redox environment and dihydroceramide desaturase activity (Fabrias et al., 2012). Similarly, in tyrosine metabolism, SL samples displayed a substantial downregulation of the hydrogen peroxide synthesis reaction (BiGG ID: H2O2syn), where the reaction is facilitated by the enzyme hydrogen peroxide synthase, which is dependent on NADPH, and converts oxygen (O₂), NADPH, and protons (H⁺) into hydrogen peroxide (H₂O₂) and NADP⁺(Figure 4C). Hydrogen peroxide is essential for modulating the cellular redox state and functions as a signaling molecule that is critical for preserving homeostasis and cellular responses. Further insights into these metabolic shifts were provided by network analysis (Figures 3A, 4A). The sphingolipid network analysis (Figure 3A), consistent with the alterations in dihydroceramide desaturase activity, reaffirmed the complex regulatory interactions among the key enzymes governing ceramide synthesis and turnover. More pertinent to the changes in tyrosine metabolism, Figure 4A depicts interactions involving DUOX2 (relevant to H₂O₂ balance) and COMTD (linked to L-DOPA metabolism), and importantly, visualized tyrosine hydroxylase (TH), the pivotal enzyme for L-DOPA production from tyrosine. PPI analysis further supported this upregulation by revealing increased translation of the TH gene (EC:1.14.16.2), highlighting its enhanced role in tyrosine metabolism, as shown in Figure 4B and Supplementary Figure S3. Increased tyrosine hydroxylase activity promoted melanogenesis, resulting in increased the production of reactive oxygen species (ROS) (Iyengar and Misra, 1987). These observations suggest that tyrosine metabolism plays a significant role in the management of oxidative stress in SL.

As illustrated in Figure 5C, our metabolic simulation analysis exhibited a considerable reduction in the PI(3,4,5)P₃ phosphatase reaction (BiGG ID: PI345P3Pn) in SL, where the reaction is facilitated by phosphatidylinositol-3,4,5-trisphosphate 3-phosphatase, an enzyme encoded by the PTEN gene. PTEN is integral to the regulation of a variety of cellular processes, including cell metabolism, by modulating Akt activity via the PTEN/PI3K/Akt signaling pathway (Downward, 2004). Further insights at the network level were provided by gene coexpression analysis (Figure 5A). This analysis highlights the key enzymes within inositol phosphate metabolism (greyish-green nodes) and their interactions. Notably, INPP5E, an inositol polyphosphate-5-phosphatase involved in modulating PI3K/Akt signaling by regulating PIP levels, is part of a large cluster interconnected with glycerophospholipid (e.g., PLCB3) and sphingolipid metabolism genes. The network also prominently featured enzymes critical for calcium signaling. Phospholipase C isoforms, such as PLCB3 (a key generator of IP₃ and diacylglycerol (DAG) from PIP₂) and PLCH2, formed a module. In contrast, ITPKC (inositol-trisphosphate 3-kinase), which influences IP₃ levels and calcium signals, was observed in a separate cluster interacting with sphingolipid and methionine metabolism genes. These visualized networks underscore the complex interplay and co-regulation of proteins involved in both PI3K/Akt pathway modulation and calcium signaling through inositol phosphates.

In SL samples, membrane fluidity maintenance and fatty acid elongation were predominantly inhibited. In the context of fatty acid metabolism, pronounced reduction of the fatty-acyl-CoA elongation reaction (BiGG ID: FAEL183) was observed, indicating a decrease in the elongation activity of fatty acids (Figure 2C). This reaction, catalyzed by a fatty-acyl-CoA elongation enzyme, converts linoleoyl-CoA to dihomo-γ-linolenoyl-CoA using malonyl-CoA as a substrate.

Additionally, ACACA and ACACB (EC:6.4.1.2), which convert acetyl-CoA to malonyl-CoA, showed suppressed expression in SL (Supplementary Figure S4). Malonyl-CoA not only acts as a key substrate for fatty acid elongation but also serves as a regulatory molecule for fatty acid oxidation (Bowman and Wolfgang, 2019). Therefore, a decreased activity of these enzymes is likely to impair both processes.

PPI analysis, as shown in Figure 2B, supported this observation by revealing that ACACA, ACACB, and multiple long-chain acyl-CoA synthetases, such as ACSBG1 and ACSL5, form a closely connected module within the acyl-CoA biosynthetic process. These enzymes, many of which were downregulated, are central to fatty acid activation, further supporting the functional suppression of SL.

Similarly, disruption of bile acid metabolism was observed. Figure 6A shows that CYP27A1, AMACR, and HSD17B4, which are essential enzymes in bile acid biosynthesis, exhibited notable interactions with genes involved in cholesterol and sphingolipid metabolism. This network suggests that bile acid metabolism does not operate in isolation but rather interfaces with broader lipid regulatory pathways.

Additionally, Figure 6B and Supplementary Figure S5 show the coordinated suppression of CH25H (EC:1.14.99.38), CYP27A1 (EC:1.14.15.15), and AMACR (EC:5.1.99.4), which act sequentially during bile acid biosynthesis. This pattern implies multi-step disruption of the pathway and potential downstream effects on cholesterol turnover.

Collectively, these findings indicate that SL is characterized by broad dysregulation of fatty acid metabolism, encompassing elongation, activation, and oxidation, likely leading to altered lipid composition and signaling.

Figure 7C shows a substantial downregulation of key cholesterol biosynthesis reactions in SL. Specifically, the squalene synthase reaction (BiGG ID: SQLSr) and squalene epoxidase endoplasmic reticular NADP reaction (BiGG ID: SQLEr) were diminished, indicating a decreased conversion of farnesyl pyrophosphate to squalene, and subsequently, squalene to 2,3-oxidosqualene. Squalene epoxidase (SQLE) is responsible for the stereospecific conversion of squalene to 2,3(S)-oxidosqualene, which is the initial oxygenation step in cholesterol biosynthesis. Its downregulation may hinder cholesterol synthesis and disturb cholesterol homeostasis (Padyana et al., 2019). In addition, DHCR71r reaction, DHCR72r reaction, DHCR241r reaction, DHCR242r reaction, and DHCR243r reaction, which are catalyzed by 7-dehydrocholesterol reductase and 24-dehydrocholesterol reductase, were reduced. These reactions are pivotal in the terminal phases of cholesterol biosynthesis and facilitate the conversion of 7-dehydrocholesterol and desmosterol to cholesterol. As shown in Figure 7A, gene network analysis further underscored the coordinated nature of cholesterol metabolism. This analysis revealed distinct clusters of functionally associated genes central to cholesterol biosynthesis (nodes predominantly colored light gray), including key enzymes, such as TM7SF2 (DHCR7), DHCR24, MVK, PMVK, and CYP51A1. These clusters demonstrate their interconnectedness within the cholesterol synthesis pathway and their links to other lipid metabolic processes, such as sphingolipid and bile acid biosynthesis.

Translation-level analysis indicated a considerable reduction in the levels of critical proteins involved in cholesterol transport, as shown in Supplementary Figure S6. Complementing these findings, PPI analysis was performed to elucidate the interaction networks of these key cholesterol transport-related proteins (Figure 7B). This network highlighted LDLR, LRP1, and LRP2 receptors crucial for LDL uptake and their association with amyloid-beta and plasma lipoprotein particle clearance (nodes colored blue and/or green, reflecting their pathway involvement). The same PPI network underscored the centrality of APOA1, the main high-density lipoprotein (HDL), in the regulation of plasma lipoprotein levels and ‘Plasma lipoprotein particle clearance (red and/or green nodes). The observed downregulation of these interconnected proteins, LDLR and LRPs for LDL uptake, and APOA1 for HDL formation and reverse cholesterol transport, significantly impaired cholesterol homeostasis. For instance, diminished APOA1 impairs HDL formation, reduces cholesterol efflux from peripheral tissues to the liver, leads to lipid dysregulation, and diminishes the ability of HDL to modulate immune cell function, potentially worsening inflammatory and autoimmune conditions (Trakaki and Marsche, 2020).