Transcriptomic Analysis Reveals Mitochondrial-Driven Subgroups in Atopic Dermatitis Lesions

Summary

Atopic dermatitis (AD) is a chronic inflammatory skin disease with significant molecular heterogeneity. This study identifies two transcriptomically distinct AD subgroups, driven by differential mitochondrial gene expression and immune infiltration, and reveals four hub genes as potential biomarkers for patient stratification.

Abstract

Atopic dermatitis (AD) is a common and chronic inflammatory skin disease with global prevalence. Its clinical heterogeneity and complex molecular mechanisms pose significant challenges for the development of effective therapies. To explore AD molecular heterogeneity using lesional skin transcriptomic data, characterize their biological and immune profiles, and identify key genes underlying the differentiation. Differentially expressed genes were identified using DESeq2, followed by pathway and co-expression analyses via GSEA and WGCNA, respectively. Mitochondrial-related genes were extracted by intersecting DEGs and WGCNA modules with the MitoCarta3.0 database, and their functional relevance was assessed through GO and KEGG enrichment. Hub genes were identified by protein-protein interaction network analysis, which were subsequently used to construct a classification model. Transcriptional regulators were predicted using hTFtarget, while immune cell infiltration was quantified using CIBERSORT. Two molecular subgroups were identified. Cluster 1 was enriched in cell signaling and adhesion pathways, whereas Cluster 2 exhibited upregulation of oxidative phosphorylation and proteasome-related processes. A total of 85 mitochondrial-associated genes, primarily involved in energy metabolism, were differentially expressed between clusters. PPI network analysis identified four hub genes (BAD, BOLA1, CHCHD5, and ISOC2), which were significantly upregulated in Cluster 1. A hub gene-based classifier demonstrated strong discriminatory power (area under the curve > 0.7). Predicted key transcriptional regulators included ATF3, BRD2, BRD4, and CEBPA. Immune profiling revealed higher regulatory T cell infiltration in Cluster 1 and increased follicular helper T cells in Cluster 2. This study reveals two molecularly and immunologically distinct AD subtypes, characterized by differential mitochondrial function and immune microenvironment signatures.

Introduction

Atopic dermatitis (AD) is a common and chronic inflammatory skin disease affecting up to 20% of children and 10% of adults¹. It is characterized by intense itch and recurrent eczematous lesions². Clinically, AD arises from a dynamic interplay of polygenic susceptibility (e.g., FLG loss-of-function mutations), immune dysregulation, and environmental exposures such as low humidity and microbial dysbiosis³. While AD shares some pathophysiological features with psoriasis, their clinical presentations are mutually exclusive, with different genetic effects in shared pathways and distinct immune alterations⁴.

AD is a heterogeneous disease, characterized by diverse transcriptome profiles across different patient groups and disease symptoms⁵. Integrated analyses of skin tissues and peripheral blood mononuclear cells have shown that clinical features such as erythema and papulation are associated with distinct immunological signatures, reflecting the interplay between local skin and systemic immune responses⁶. Large-scale transcriptomic studies further highlight the role of IL-13 pathways in AD pathogenesis, while AD exhibits greater molecular heterogeneity than psoriasis, with variations in gene expression patterns linked to disease severity, age of onset, and genetic background^5,7. These differences underscore the complexity of AD pathogenesis and the need for personalized approaches in research and therapy.

Mitochondrial proteins play an important role in AD pathogenesis through dysregulated oxidative stress and metabolic pathways. Studies reveal increased mitochondrial complex I and II activity in nonlesional AD keratinocytes, driving excessive oxidation of long-chain fatty acids and elevated ROS production, which increases epidermal barrier dysfunction^8,9. Concurrently, proteomic analyses identify reduced NRF2-antioxidant pathway proteins and mitochondrial components in AD epidermis, reducing oxidative stress resolution¹⁰. Mitochondrial DNA damage further contributes to inflammatory responses, while interventions like using mitochondrial-targeted antioxidants demonstrate efficacy in restoring epidermal homeostasis by mitigating ROS^11,12. These findings underscore mitochondrial proteins as both drivers of AD pathology and potential therapeutic targets.

Recent transcriptome analyses of AD have significantly advanced the understanding of its genetic architecture and molecular heterogeneity. The first RNA sequencing profiling of AD had revealed the increased expression in the TREM-1 pathway and the IL-36 cytokine¹³. Weighted gene co-expression network analysis based on gene expression profile has further revealed distinct molecular modules and hub genes, such as HSPA4, LCE3E, and LCE3D, that orchestrate inflammatory responses and keratinization, highlighting potential therapeutic targets¹⁴. These studies underscore the value of multi-tissue transcriptomics and polygenic risk modeling in refining disease prediction and uncovering the complex underpinnings of AD. However, existing studies have predominantly focused on overall AD transcriptomics without specifically resolving the role of mitochondrial genes in defining molecular subgroups. The present study advances beyond prior transcriptomic analyses by integrating multiple GEO datasets to identify consensus-clustering-based AD subgroups and systematically intersecting differentially expressed genes, WGCNA modules, and a curated mitochondrial gene catalog to pinpoint hub mitochondrial genes that define subgroup identity and may serve as novel biomarkers.

With the hypothesis that lesional AD skin harbors molecularly distinct transcriptomic subgroups characterized by differential mitochondrial gene expression, which may underpin the clinical heterogeneity of AD. To test this, this study integrated RNA sequencing data from published studies and collected gene expression data of lesion skin samples from 266 AD patients. Molecular subgroups were identified based on consensus clustering of gene expression, and the gene expression was compared between the two molecular subgroups. The genes driving this difference show functional enrichment in cell signaling and oxidative phosphorylation. Furthermore, mitochondrial genes that differentiate the two molecular groups were examined, and the key hub genes were identified within a protein-protein interaction network. These findings highlight the genetic heterogeneity of AD disease and expand the knowledge of the roles of mitochondrial proteins in AD pathology.

Protocol

This study used publicly available gene expression datasets from the Gene Expression Omnibus (GEO) database. No patient-identifiable data were accessed, and no new patient samples were collected. Therefore, no institutional review board (IRB) approval or patient consent was required for this secondary analysis of publicly available data. The software and the databases used are listed in the Table of Materials.

1 Data and resources

Transcriptome data of atopic dermatitis (AD) patients were obtained from the GEO database, including four studies: GSE121212 (N = 55)⁵, GSE157194 (N = 57)¹⁵, GSE193309 (N = 111)¹⁶ and GSE277961 (N = 43)¹⁷. All datasets contained raw or pre-normalized count data from lesional skin biopsies. Raw count matrices were downloaded and integrated across studies. Cross-study batch effects were corrected using the ComBat-seq method (sva R package, v3.44.0) to harmonize expression profiles across the four GEO datasets prior to downstream analysis. All four datasets are based on RNA-seq performed on human skin biopsy specimens and were aligned to the human reference genome GRCh38 using study-specific pipelines. Gene-level expression quantification was performed using Ensembl gene annotation (v105).

2 Consensus clustering

Unsupervised consensus clustering was conducted using the ConsensusClusterPlus R package (v1.64.0)¹⁸ to stratify 266 lesional skin samples from AD patients based on transcriptomic profiles. Prior to clustering, raw count data from all studies were merged, and cross-study batch effects were corrected using the removeBatchEffect function from the limma package. Gene expression data were then variance-stabilizing transformed (VST) using DESeq2 and filtered to retain the top 5,000 most variable genes. Clustering was conducted using hierarchical clustering with Pearson correlation and average linkage across 1,000 iterations, subsampling 80% of samples per iteration. The best number of clusters (k ranges from 2 to 10) was determined by evaluating the consensus cumulative distribution function (CDF), delta area plots, and cluster-consensus scores. The resulting clusters were validated using consensus heatmaps and PCA and used in downstream biological and clinical analyses.

3 Differential gene expression analysis

The gene expression levels were compared between molecular subgroups using the DESeq2 R package (v1.46.0)¹⁹. Raw count data were input to estimate gene-wise dispersion and fit a negative binomial model. Differentially expressed genes (DEGs) were identified using the Wald test, and results were filtered using a significance threshold of adjusted p-value < 0.01 and absolute log2 fold change > 1.

4 Gene set enrichment analysis

Gene set enrichment analysis (GSEA) was performed using the clusterProfiler R package (v4.12.6)²⁰. All genes were ranked by their log2 fold change from the differential expression analysis, and the GSEA function was applied with parameters eps = 0, minGSSize = 10, and maxGSSize = 500, while other settings were kept at default. Enrichment was performed against the MSigDB Hallmark gene sets. For each molecular cluster (Cluster 1 and Cluster 2), the top three enriched pathways were selected based on nominal p-value and normalized enrichment score (NES). Results were visualized using the GseaVis R package (v0.1.0)²¹.

5 Weighted gene co-expression network analysis

WGCNA was performed on the gene expression profile of AD patients to identify gene co-expression modules associated with transcriptomic subtype identity using the R package WGCNA (v1.73)²². Genes were filtered to retain the top 75% with the highest variance across all samples. A signed co-expression network was constructed using a soft-thresholding power selected from 1 to 30 to approximate scale-free topology. Gene modules were identified by hierarchical clustering and dynamic tree cutting. Module-trait associations were studied by correlating co-expressed module eigengenes with the molecular subtype labels, and modules showing significant correlation (p < 0.05) with molecular subtype were selected for downstream analysis.

6 Mitochondrial proteins in AD molecular subgroups

To identify mitochondrial-associated genes within the DEGs and WGCNA modules, these gene sets were overlapped with the curated list of mitochondrial proteins from MitoCarta3.0²³. The intersecting genes were considered putative mitochondrial proteins relevant to the AD disease context.

7 Gene ontology and KEGG enrichment analysis

To identify the major cell functions and biological processes that differentiate molecular subgroups. GO and KEGG enrichment analyses were performed for mitochondrial-associated genes using clusterProfiler. GO enrichment was conducted separately for Biological Process (BP), Cellular Component (CC), and Molecular Function (MF) categories using the enrichGO function with OrgDb = “org.Hs.eg.db”, ont = "ALL", and default parameters. KEGG pathway enrichment was performed using the enrichKEGG function with the organism set to "hsa". Enriched terms with adjusted p-value < 0.05 were considered as significant.

8 Protein-protein interaction networks

The 85 DEGs and AD-associated module overlapped mitochondrial genes were queried in the STRING database (https://string-db.org)²⁴ to retrieve known and predicted PPIs. Network topology analysis was performed using seven measures (Degree, Closeness, Betweenness, Eigenvector, PageRank, Hub, and Authorities) to rank the node importance of each protein. The top 30-ranked genes of each measure were selected, and the intersections across all seven approaches were visualized using an UpSet plot. The 4 intersecting genes from the measures were used to build a classification model based on their gene expression levels. Model accuracy, sensitivity, specificity, and area under the ROC curve (AUC) were calculated to evaluate classification performance.

9 Transcriptional regulation analysis

The hTFtarget database (http://bioinfo.life.hust.edu.cn/hTFtarget)²⁵ was queried to retrieve experimentally supported TF-target interactions for four intersecting hub genes. The resulting TF-gene regulatory network was constructed and visualized using the igraph²⁶ and ggraph²⁷ R packages.

10 Bulk RNA-seq immune cell infiltration analysis

The CIBERSORT²⁸ algorithm (https://cibersort.stanford.edu/) was used to estimate the relative proportions of 22 immune cell types in lesional skin transcriptome data. Normalized gene expression data were input into the CIBERSORT (v0.1.0) along with the LM22 signature matrix. The analysis was run with 1,000 permutations and quantile normalization disabled. Samples with CIBERSORT output p-values < 0.05 were considered for downstream analysis. The estimated immune cell fractions were compared between molecular subgroups using the Wilcoxon rank-sum test with FDR correction for multiple comparisons across the 22 immune cell types (p.adjust.method = “FDR”), and results were visualized with boxplots using the ggplot2²⁹ R package.

Results

Transcriptional subgroups in atopic dermatitis

RNA-seq data from 266 AD patient samples were analyzed to investigate transcriptional heterogeneity within the disease. After quality control and correction for batch effects across multiple studies, unsupervised consensus clustering revealed two distinct molecular subgroups (Figure 1A). Cluster stability and optimal cluster number were assessed using the cumulative distribution function (CDF) plot (Figure 1B), delta area plot (Figure 1C), and consensus matrix heatmap (Figure 1D). Together, these results support the existence of two robust transcriptional subtypes in AD, reflecting underlying genetic heterogeneity.

Differentially expressed genes between AD subgroups

The t-SNE plot based on the normalized gene expression matrix further validated the transcriptional subgroups identified by consensus clustering. The t-SNE plot revealed two well-separated clusters, each corresponding to one of the previously defined subgroups (Figure 2A), supporting the presence of distinct molecular profiles among AD patients. Differential expression between the two subgroups was then analyzed using DESeq2, with a threshold of adjusted p < 0.01 and |log₂ fold change| > 1. The resulting volcano plot (Figure 2B) showed differentially expressed genes (DEGs), indicating strong transcriptional divergence. The top 10 up-regulated genes are ABHD2, ADAR, ADCY3, ADCY9, ADD1, ADIPOR2, AFF1, AGFG1, AGRN and AHNAK in cluster 1 and C2orf68, CTTN, GPR108, HERPUD1, LRPAP1, MAP1LC3B2, NKIRAS2, NR1H2, PDE5D and PMPCA in cluster 2 (Figure 2C).

AD subgroup associated gene set

Gene Set Enrichment Analysis (GSEA) revealed distinct functional enrichment profiles between the two transcriptomic clusters (Figure 3A). Cluster 1 exhibited significant enrichment for pathways involved in cell signaling and adhesion, including focal adhesion (Figure 3B) and MAPK signaling pathway (Figure 3C), suggesting an active state characterized by enhanced cell-cell-extracellular matrix interactions and proliferation. In contrast, Cluster 2 showed strong enrichment for oxidative phosphorylation (Figure 3D) and proteasome function (Figure 3E), suggesting an active oxidative and proteolytic metabolic phenotype.

Co-expressed genes in AD molecular groups

To identify co-expression modules associated with transcriptomic subtypes, WGCNA was performed after data preprocessing. Outlier samples were first identified and removed based on hierarchical clustering of sample distances to ensure the robustness of downstream network construction (Figure 4A). A soft-thresholding power was then selected using the scale-free topology criterion, with a power of 6 chosen to achieve a scale-free R² > 0.85 (Figure 4B). Gene modules were identified through hierarchical clustering and dynamic tree cutting, followed by eigengene clustering to merge closely related modules (Figure 4C and Figure 4D). The resulting gene network was visualized using a heatmap of topological overlap, confirming the presence of distinct gene co-expression patterns (Figure 4E). Module-trait relationship analysis revealed a strong and significant correlation between the MEyellow module (N_gene = 743) and the molecular sub-group (Figure 4F).

Functional enrichment of the AD subgroup associated mitochondrial genes

GO and KEGG enrichment analysis was then performed on the intersected genes (N = 85) among DEGs, genes in the MEyellow module, and a list of mitochondrial proteins from MitoCarta3.0 (Figure 5A) to investigate the functional roles of mitochondrial-associated genes driving the transcriptional differences between clusters. KEGG pathway analysis identified enrichment in oxidative phosphorylation and metabolic pathways (Figure 5B). GO enrichment analysis revealed significant overrepresentation of terms associated with mitochondrial function, including proton motive force-driven mitochondrial ATP synthesis, respiratory chain complex, and NADH dehydrogenase activity (Figure 5C), suggesting that these genes are mainly associated with mitochondrial energy metabolism and energy regulation.

Hub mitochondrial genes in AD molecular differentiation

To identify the key genes in the 85 mitochondrial transcriptome-subgroup-associated genes, a PPI network was constructed (Figure 6A). Based on the ranking in each of the seven topological measures (See methods), the top 30 genes were selected, and their intersections were analyzed and visualized in an UpSet plot (Figure 6B). This analysis resulted in four hub genes consistently identified as central nodes based on all ranking criteria (BAD, BOLA1, CHCHD5, ISOC2). Their expression profiles were examined across the two transcriptomic clusters and found that all four hub genes were significantly upregulated in Cluster 1 compared to Cluster 2 (Figure 6C). Pairwise gene expression correlation analysis showed positive correlations among all four genes, indicating coordinated regulation, with CHCHD5 and ISOC2 exhibiting the strongest correlation (Figure 6D). Furthermore, a classification model we built using the expression of these four genes demonstrated robust discriminatory power between the two clusters, with the ROC curve showing an area under the curve (AUC) > 0.7 (Figure 6E). Additionally, transcription factor (TF) regulatory network analysis was performed to study the regulatory mechanisms governing the expression of the four identified hub genes. All known and predicted transcription factors that potentially regulate these hub genes were queried from hTFtarget, and the results were integrated and visualized as a transcriptional regulatory network (Figure 7). In the TF-hub gene regulatory network, BAD had the highest number of TFs, and ATF3, BRD2, BRD4, and CEBPA each interacted with all four hub genes, suggesting a shared regulatory mechanism.

Comparison of immune cell infiltration between AD subgroups

To investigate the immunological landscape associated with the transcriptomic subgroups, immune cell infiltration analysis using CIBERSORT was performed, which estimates the relative proportions of 22 immune cell types based on bulk transcriptome data (Figure 8). Among the immune subsets, regulatory T cells (Tregs) were found to be significantly more abundant in Cluster 1, suggesting an immunosuppressive microenvironment potentially associated with mitochondrial activity and signaling pathways upregulated in this group. In contrast, follicular helper T cells were significantly enriched in Cluster 2, indicating a potentially more active adaptive immune response in this subgroup.

DATA AVAILABILITY:

The transcriptomic data analyzed in this study are publicly available in the Gene Expression Omnibus (GEO) repository under accession numbers GSE121212, GSE157194, GSE193309, and GSE277961 (https://www.ncbi.nlm.nih.gov/geo/).

Figure 1: Consensus clustering of atopic dermatitis lesion samples based on transcriptomic profiles. (A) Heatmap and hierarchical clustering of the consensus matrix across atopic dermatitis (AD) samples. (B) Consensus cumulative distribution function (CDF) plot used to determine the optimal number of clusters (k = 2–10). (C) Delta area plot showing the relative change in the area under the CDF curve for each k. (D) Consensus cluster assignment for k = 2. Each column represents an individual sample, and colors indicate cluster membership (Cluster 1, red; Cluster 2, teal). Please click here to view a larger version of this figure.

Figure 2: Differential gene expression of atopic dermatitis molecular subtypes. (A) t-distributed stochastic neighbor embedding (t-SNE) plot of AD samples. Each point represents a sample projected into two dimensions and is coloured according to cluster assignment. (B) Volcano plot of differentially expressed genes (DEGs) between Cluster 1 and Cluster 2. Each point represents a gene, plotted by log2 fold change (x-axis) and −log10 adjusted p-value (y-axis). Red and blue points indicate significantly upregulated genes in Cluster 1 and Cluster 2, respectively, while grey points indicate non-significant genes. (C) Heatmap of the top 10 highly expressed genes in Cluster 1 and Cluster 2. Please click here to view a larger version of this figure.

Figure 3: Gene set enrichment analysis of atopic dermatitis molecular subtypes. (A) Two-sided bar plot showing GSEA results between Cluster 1 and Cluster 2, with the top significantly enriched pathways. Pathways enriched in Cluster 1 are shown on the right, and those enriched in Cluster 2 are shown on the left. (B–E) Representative enrichment plots for focal adhesion, MAPK signaling pathway, oxidative phosphorylation, and proteasome pathways. Please click here to view a larger version of this figure.

Figure 4: Weighted gene co-expression network analysis of atopic dermatitis samples. (A) Sample clustering dendrogram based on gene expression profiles. (B) Scale-free topology fit index and mean connectivity across soft-thresholding powers (1–30). (C) Clustering and heatmap of module eigengenes, with colors indicating pairwise correlations. (D) Hierarchical clustering dendrogram showing genes grouped into co-expressed modules. (E) Heatmap of the topological overlap matrix (TOM), representing co-expression similarity between gene pairs. (F) Heatmap of module–trait relationships, showing correlations between module eigengenes and clinical traits. Correlation coefficients are displayed within each cell, and color intensity indicates the strength and direction of the correlation (red, positive; blue, negative). Please click here to view a larger version of this figure.

Figure 5: Gene ontology and KEGG pathway enrichment of subtype-associated mitochondrial genes. (A) Venn diagram showing overlap among DEGs, subtype-associated module genes, and mitochondrial genes. (B) Bubble plot of the top 20 enriched KEGG pathways of intersecting genes. (C) Top 10 enriched Gene Ontology (GO) terms for Biological Process (BP), Cellular Component (CC), and Molecular Function (MF). All enrichment analyses were performed using an adjusted p-value < 0.05 (false discovery rate) as the significance threshold. Please click here to view a larger version of this figure.

Figure 6: Protein–protein interaction analysis and identification of hub genes. (A) Protein–protein interaction (PPI) network of 85 intersecting genes. Nodes represent proteins, and edges indicate predicted or experimentally validated interactions from the STRING database. (B) UpSet plot showing intersections among the top 30 ranked genes based on seven network centrality measures. (C) Boxplots showing expression levels of four hub genes in Cluster 1 and Cluster 2. (D) Pairwise correlation analysis among the four hub genes. (E) Receiver operating characteristic (ROC) curve showing classification performance, with sensitivity plotted against specificity. The area under the curve (AUC) indicates overall accuracy. Statistical significance in panel (C) was assessed using the Wilcoxon rank-sum test (*p < 0.05). Please click here to view a larger version of this figure.

Figure 7: Regulatory network of hub genes. Red circles represent hub genes and blue circles represent associated transcription factors (TFs). Edges indicate regulatory interactions. The size of each hub gene node reflects the number of interacting TFs. Please click here to view a larger version of this figure.

Figure 8: Comparison of immune cell infiltration between atopic dermatitis subgroups. Boxplots showing the estimated proportions of 22 immune cell types in each cluster (Cluster 1, red; Cluster 2, teal). Asterisks (*) indicate statistically significant differences between clusters, assessed using the Wilcoxon rank-sum test with false discovery rate correction for multiple comparisons. Please click here to view a larger version of this figure.

Discussion

Atopic dermatitis is a prevalent and genetically heterogeneous inflammatory skin disorder, posing significant challenges for personalized treatment strategies. This study provides insights into the molecular heterogeneity of AD by identifying two distinct transcriptional subgroups in 266 lesional skin samples of AD patients. These subgroups exhibit different biological function enrichment and immune cell infiltration, with Cluster 1 characterized by active cell signaling and adhesion pathways and an enrichment of regulatory T cells (Tregs), while Cluster 2 is defined by increased oxidative phosphorylation and proteasome function, alongside an increase in follicular helper T cells. Importantly, four mitochondrial-associated genes (BAD, BOLA1, CHCHD5, and ISOC2) were identified that are upregulated in Cluster 1. They are hub genes in the PPI network of 85 genes that define the two transcriptomic subtypes and can be regulated simultaneously by transcription factors, including ATF3, BRD2, BRD4, and CEBPA.

The two distinct clusters were identified based on the gene expression pattern, which suggests that AD exhibits significant gene expression heterogeneity. Previous studies have shown that this heterogeneity could be driven by genetic, epigenetic, and immune-mediated mechanisms that define distinct molecular endotypes^5,30,31. Between the two clusters, different immune cell infiltration was also observed, cluster 1 characterized by regulatory T cells and cluster 2 enriched with follicular helper T (Tfh) cells, implying underlying immunological heterogeneity and potentially different disease mechanisms in the AD cohort. A Treg-dominated cluster suggests an immune environment where regulatory mechanisms are prominent, possibly reflecting attempts to control inflammation or a compensatory response to chronic immune activation^32,33. However, in AD, Treg function can be impaired despite increased numbers, which may not effectively suppress inflammation. In contrast, a Tfh-enriched cluster points toward enhanced B cell help, increased germinal center activity, and likely elevated IgE production, which are hallmarks of allergic responses and more severe or extrinsic forms of AD^34,35. Tfh cells are known to support B cell differentiation and antibody class switching, and their expansion correlates with disease activity and allergic sensitisation³⁶. The presence of these distinct clusters may reflect different clinical phenotypes, disease severities, or responses to therapy, and highlights the importance of personalized approaches in AD research and treatment.

It is known that mitochondrial dysfunction plays a significant role in the pathogenesis of AD through mechanisms such as maintaining the epidermal barrier³⁷, producing reactive oxygen species³⁸, and regulating immune cell response³⁹. Four mitochondrial proteins were identified as central nodes of genes associated with transcriptomic differentiation, which can predict the molecular subtype with high accuracy (AUC>0.7). Although no previous study has shown the direct link of these genes to AD, these results suggest they are potential biomarkers for stratifying AD patients and assessing mitochondrial dysfunction in AD patients. BAD (BCL2-associated agonist of cell death) is a pro-apoptotic member of the BCL-2 family involved in mitochondrial apoptotic signaling; its upregulation in Cluster 1 may reflect enhanced mitochondrial apoptotic priming in this subtype. BOLA1 is a mitochondrial protein implicated in iron-sulfur cluster biogenesis and oxidative stress regulation. CHCHD5 (coiled-coil-helix-coiled-coil-helix domain containing 5) is a mitochondrial inner membrane protein associated with cristae organization and electron transport chain efficiency. ISOC2 (isochorismatase domain containing 2) has been linked to metabolic processes and mitochondrial function. The coordinated upregulation of these four genes in Cluster 1 suggests a state of heightened mitochondrial metabolic activity and apoptotic signaling in this AD subgroup.

This study has several limitations. While this study identified two molecular subgroups with distinct gene expression in mitochondrial genes and separated immune cell infiltration, we lack detailed clinical phenotype data that would allow us to correlate the stratification with disease severity and longitudinal treatment responses. To fully understand what the high expression of these four hub genes in Cluster 1 means, future studies should integrate comprehensive clinical metadata and ideally perform functional validation in keratinocyte or mouse models. In addition, this study relies on publicly available bulk RNA-seq data; although it is powerful in large-scale analysis, it lacks resolution in specific cell types. This limitation makes it challenging to determine whether the observed differential expression of mitochondrial genes originates from keratinocytes, infiltrating immune cells, or other skin-resident cell populations. With the wide use of single-cell and spatial transcriptomic approaches, future research would benefit a lot from the power of deconvoluting the expression signals and providing a more precise overview of the transcriptomic profile at the cellular level. Additionally, all data were obtained from punch biopsy specimens, which sample full-thickness skin. Future studies incorporating tape-strip RNA-seq may offer a complementary non-invasive approach to profiling the superficial epidermal transcriptome and could validate the identified subgroup signatures in a minimally invasive setting. Future integration of single-cell RNA-seq and spatial transcriptomic data would further advance the resolution of cell-type-specific mitochondrial signatures in AD skin.

In conclusion, this study examined the transcriptomic heterogeneity of AD in 266 samples, identified key mitochondrial genes and unique immune cell infiltration, differentiating the subgroups. These findings improve the understanding of AD molecular heterogeneity and highlight the potential for developing precision treatment approaches based on specific molecular profiles.

Materials

List of materials used in this article

Name	Company	Catalog Number	Comments
CIBERSORT	Stanford University	v0.1.0; immune cell deconvolution; https://cibersortx.stanford.edu
clusterProfiler	Bioconductor	v4.12.6; GSEA and GO/KEGG enrichment analysis; https://bioconductor.org/packages/clusterProfiler
ConsensusClusterPlus	Bioconductor	v1.64.0; unsupervised consensus clustering; https://bioconductor.org/packages/ConsensusClusterPlus
DESeq2	Bioconductor	v1.46.0; differential gene expression analysis; https://bioconductor.org/packages/DESeq2
Gene Expression Omnibus (GEO)	NCBI	Public transcriptomic data repository; datasets GSE121212, GSE157194, GSE193309, GSE277961; https://www.ncbi.nlm.nih.gov/geo
ggplot2	CRAN	Data visualisation; https://ggplot2.tidyverse.org
ggraph	CRAN	Graph and network visualisation; https://ggraph.data-imaginist.com
GseaVis	GitHub (junjunlab)	v0.1.0; GSEA visualisation; https://github.com/junjunlab/GseaVis
hTFtarget	Huazhong University of Science and Technology	Human transcription factor target database; http://bioinfo.life.hust.edu.cn/hTFtarget
igraph	CRAN	Network construction and visualisation; https://igraph.org
limma	Bioconductor	removeBatchEffect function; https://bioconductor.org/packages/limma
MitoCarta3.0	Broad Institute	Curated mitochondrial protein database; https://www.broadinstitute.org/mitocarta
MSigDB (Hallmark gene sets)	Broad Institute	Gene set database for GSEA; https://www.gsea-msigdb.org/gsea/msigdb
org.Hs.eg.db	Bioconductor	Human genome annotation database; https://bioconductor.org/packages/org.Hs.eg.db
R	R Core Team	Statistical computing environment; https://www.r-project.org
STRING	EMBL	v12.0; protein-protein interaction database; https://string-db.org
sva (ComBat-seq)	Bioconductor	v3.44.0; batch effect correction; https://bioconductor.org/packages/sva
WGCNA	CRAN	v1.73; weighted gene co-expression network analysis; https://cran.r-project.org/package=WGCNA