Bioinformatics‑Based Analysis Reveals Diagnostic Biomarkers and Immune Landscape in Atopic Dermatitis

Introduction

Atopic dermatitis, also referred to as atopic eczema, is a prevalent chronic inflammatory skin disorder that affects millions worldwide, with particularly high incidence in pediatric populations.¹ This multifactorial condition presents clinically with hallmark features including xerosis (dry skin), erythema, pruritus (intense itching), and recurrent eczematous lesions. The underlying pathogenesis involves a complex interplay of genetic susceptibility, environmental triggers, epidermal barrier dysfunction, immune system dysregulation, and altered microbial colonization.²

A pivotal factor in the pathogenesis of atopic dermatitis is compromised skin barrier function, which exacerbates cutaneous hypersensitivity to external irritants and allergens.^3,4 This impairment is primarily attributed to the downregulation of key genes responsible for the biosynthesis of intercellular lipids, such as long-chain fatty acids, as well as deficiencies in natural moisturizing factors (NMFs).^5,6

The immunopathogenesis of atopic dermatitis is characterized by a predominant Th2-type immune response, mediated primarily by IL-4, IL-13, and IL-31 cytokines.^7,8 These mediators drive IgE production, promote inflammatory cell recruitment, and perpetuate chronic inflammation. Furthermore, Th22 and Th17 immune pathways are upregulated in lesional skin, exacerbating epidermal barrier dysfunction and triggering acute disease exacerbations.⁹

The diagnosis of AD primarily relies on clinical examination and patient history, with treatment aiming to alleviate symptoms and restore the skin barrier. While previous therapeutic approaches were limited to topical corticosteroids and calcineurin inhibitors, current treatment strategies have undergone significant transformation. Biologics such as dupilumab (anti-IL-4Rα), which selectively inhibits Th2 cytokine signaling and correlates with reduced serum TARC/CCL17 levels,¹⁰ have revolutionized AD management, achieving a 50% improvement in disease severity (EASI-50) in refractory cases.¹¹ Janus kinase (JAK) inhibitors (eg, abrocitinib, upadacitinib) are oral small molecules that inhibit the JAK/STAT pathway by selectively blocking JAK1.¹² However, 30–40% of patients show suboptimal responses, underscoring the need for better predictive biomarkers and targeted therapies.¹³ FDA approvals of tralokinumab (anti-IL-13) and lebrikizumab (anti-IL-13Rα1) highlight the shift toward cytokine-specific biologics.¹⁴ Emerging multi-omics studies have revealed dysregulation of genes involved in epidermal differentiation (FLG, CLDN1) and neuroimmune signaling (IL31RA), though their diagnostic utility and mechanistic links with immune infiltration remain unexplored.¹⁵ Novel biomarkers like serum periostin, TSLP and CCL17/TARC may predict treatment response and disease progression.

Given the substantial impact of atopic dermatitis on patients’ quality of life, developing novel therapeutic strategies is imperative. We hypothesize that key genes driving immune dysregulation in AD can be identified through an integrated bioinformatics framework combining WGCNA, machine learning, and Mendelian randomization - an approach distinct from previous single-method analyses. This study seeks to elucidate critical genes governing immune dysregulation in AD by integrating bioinformatics with experimental validation, further investigating their roles in immune cell infiltration, with the ultimate goal of providing clinicians with novel diagnostic and therapeutic insights.

Methods

Data Download

The National Center for Biotechnology Information Gene Expression Omnibus (NCBI-GEO) (https://www.ncbi.nlm.nih.gov/geo/) is a publicly accessible database established in 2000, hosting transcriptomic data from various global institutions.¹⁶ The GSE32924 mRNA expression profile was retrieved from the GEO database for analysis. Data normalization and log2 transformation were performed on this dataset, which contains gene expression data from 25 AD samples and 8 normal control samples, generated using GPL570 [transcript (gene) version].

Differential Expression Analysis

The normalized expression matrix from GSE32924 was used to identify DEGs. Probes were annotated according to the dataset’s annotation file, and the R package limma was used to identify DEGs. Genes with an adjusted P-value < 0.05 and an absolute fold-change (FDR > 1) were considered DEGs, following established thresholds for transcriptomic studies.^17,18 Heatmaps, volcano plots, and box plots were generated using the R packages heatmap and ggplot2.

Functional and Pathway Enrichment Analysis

Functional enrichment analysis of DEGs was conducted using the clusterProfiler package in R, with annotations from the GO and KEGG databases.¹⁹ GO analysis categorized enriched terms into biological processes (BP), cellular components (CC), and molecular functions (MF), offering insights into the biological roles of the DEGs. KEGG pathway analysis highlighted relevant signaling pathways.

Weighted Gene Co-Expression Network Analysis

WGCNA was conducted to explore gene co-expression networks.²⁰ The analysis environment was set up with necessary R packages, and the gene expression matrix and sample information were loaded. Data preprocessing involved log transformation, normalization, and exclusion of genes with low variability. Outlier samples were identified and removed, and a similarity network among genes was constructed. The dynamic tree-cutting algorithm (minModuleSize = 50, deepSplit = 2) was employed to initially partition the modules, followed by merging highly similar modules based on module eigengene clustering (with a merging threshold of MEDissThres = 0.25), with module-trait correlations visualized in heatmaps. Gene module membership (MM) and gene significance (GS) were calculated, and their relationships were illustrated in scatter plots. GS, MM data, and gene lists for each module were outputted Regression Method to Identify Key Genes.²¹

LASSO Regression Method to Identify Key Genes

Key genes were identified using the glmnet package for LASSO regression. After data preprocessing, the glmnet function performed LASSO regression fitting, with parameters set to “binomial” for binary classification and alpha to 1 for L1 regularization. Cross-validation was conducted using cv.glmnet, and plots for LASSO regression and cross-validation were saved as LASSO.pdf and cvfit.pdf. Genes with non-zero LASSO regression coefficients were considered potential key genes associated with AD.

SVM-RFE Method to Screen Feature Genes

SVM-RFE was conducted using the e1071, kernlab, and caret packages to screen feature genes. After data preprocessing, the rfe function implemented SVM-RFE with cross-validation (cv) and a radial basis kernel function (svmRadial). The feature selection results were visualized, with the minimum Root Mean Square Error (RMSE) value marked to indicate the optimal number of features.

Immune Cell Infiltration

The CIBERSORT method, using linear support vector regression, was employed to analyze immune cell infiltration in 22 human immune cell subtypes, comparing AD samples and healthy controls. Immune cells showing significant infiltration differences were identified, and correlations between immune cells and key gene expression were assessed using the Spearman correlation method.

RT-qPCR Validation of Key Genes

To validate bioinformatics findings, blood samples from 3 non-AD patients and 3 AD patients were collected for RT-qPCR. This study was approved by the Ethics Committee of Changzhou Children’s Hospital (approval number: 2023–002), and informed consent was obtained from all participants. Peripheral blood mononuclear cells (PBMCs) were isolated, and total RNA was extracted using an RNA Extraction Kit (Fastagene, Shanghai, China). RT-qPCR was performed using SYBR Green Real-Time PCR Master Mix (Toyobo, Beijing, China) following the manufacturer’s instructions. GAPDH served as the internal control, and the relative mRNA expression levels were calculated by the 2^−ΔΔCt method. Statistical significance was assessed using one-way ANOVA (P < 0.05).

ceRNA Network

miRanda was used for miRNA-mRNA target prediction, incorporating sequence matching and secondary structure to predict miRNA-mRNA interactions. Additionally, miRDB was used to provide experimentally supported miRNA-mRNA interactions using machine learning algorithms. TargetScan was used to predict miRNA binding sites on lncRNAs. The integration of miRanda, miRDB, and TargetScan results enabled the construction of a comprehensive mRNA-miRNA-lncRNA ceRNA regulatory network.

Mendelian Randomization

Two-sample Mendelian Randomization (MR) analysis was performed to investigate causal relationships between 731 immune cell traits and AD.²² MR uses genetic variation as a proxy for risk variation and employs validated instrumental variables (IV) to meet three key assumptions for causal inference: (1) exposure is directly related to genetic variation; (2) no genetic confounders exist between exposure and outcome; (3) genetic factors do not influence the outcome through unrelated channels to exposure. The IV tools (version v1.90) and single nucleotide polymorphisms (SNPs) revised with an LD r² threshold <0.1 within 500 kb were used, with significance threshold set at 5 × 10⁻⁸. To evaluate IV strength, the F-statistic was calculated. This analysis examined 731 immune phenotypes, including relative cell counts, morphological parameters, absolute cell counts, and median fluorescence intensity representing surface antigen levels from the GWAS catalog (accession numbers GCST0001391 to GCST0002121).

Statistical Analysis

All statistical analyses were conducted using R Statistical Software (version 4.3.3). The significance of fold changes in the microarray data was assessed using t-tests. Unless otherwise specified, results with a p-value < 0.05 were considered statistically significant.

Results

Differential Analysis of Atopic Dermatitis

To identify DEGs in AD, we analyzed mRNA expression profiles from AD and normal tissue samples in the GEO database (GSE32924). Comparison with normal tissues revealed 381 DEGs in AD samples, with 164 genes upregulated (log₂ FC > 1) and 217 genes downregulated (log₂ FC < −1). These DEGs are illustrated in the volcano plot and heatmap (Figure 1A and B), clearly showing the overall expression patterns.

Figure 1 Identification of DEGs in Atopic Dermatitis. (A) The volcano plot displays the expression of DEGs between atopic dermatitis and normal samples. (B) The heatmap illustrates the top 50 upregulated DEGs and top 50 downregulated DEGs. (C) The PCA plot demonstrates clear separation between atopic dermatitis and normal groups.

To evaluate data consistency within the sample groups, Principal Component Analysis (PCA) was performed. The PCA results demonstrated strong reproducibility of data in the GSE32924 dataset, showing a clear separation between samples from the AD group and the control group (Figure 1C), indicating robust grouping and reliable data distribution for further analysis.

Enrichment Analysis of Differentially Expressed Genes

Enrichment analysis was conducted to explore the functional roles of DEGs in atopic dermatitis. GO analysis covered three main categories: Biological Process, Cellular Component, and Molecular Function. Our findings revealed significant enrichment in multiple biological processes associated with skin and immune function, including epidermis development, skin development, skin barrier establishment, epidermal cell differentiation, and mitotic nuclear division. These processes highlight the genes’ involvement in skin structure and barrier integrity, crucial in AD pathology. Within the Cellular Component category, DEGs were significantly enriched in components such as the cornified envelope, lysosomal lumen, and endocytic vesicle lumen, which are integral to maintaining skin barrier function and immune response. The Molecular Function analysis revealed enrichment in gene sets associated with growth factor receptor binding, peptidase regulator activity, and copper ion binding, suggesting roles in cellular signaling and immune regulation in AD (Figure 2A).

Figure 2 Functional Enrichment Analysis of DEGs. (A) GO analysis of DEGs. The top 10 functional enrichments in BP, CC, and MF analyses. (B) KEGG analysis of DEGs. Abbreviations: NOD-like, Nucleotide-binding oligomerization domain -like; ABC transporters, ATP-binding cassette transporters.

KEGG pathway enrichment analysis further highlighted several pathways implicated in AD, including sphingolipid metabolism, cell cycle, NOD-like receptor signaling pathway, arachidonic acid metabolism, and Toll-like receptor signaling pathway (Figure 2B). Additional pathways, such as apoptosis, ABC transporters, and cytokine-cytokine receptor interactions, underscore the complex interplay of immune response, metabolic regulation, and cell cycle processes in AD onset and progression. Collectively, these enriched pathways provide insight into the underlying molecular mechanisms of immune response and barrier dysfunction in AD, offering potential targets for therapeutic intervention.

Construction of Weighted Gene Co-Expression Network

Using the WGCNA package in R, we constructed a scale-free co-expression network to explore gene modules associated with AD and healthy samples. The optimal soft-threshold power was determined to be 10, achieving a scale-free topology index of 0.85, which ensured good network connectivity and reliability (Figure 3A and B).

Figure 3 WGCNA Analysis and Identification of Candidate Hub Genes. (A) The soft-threshold power for WGCNA. (B) The average connectivity for WGCNA. (C) The clustering dendrogram for WGCNA. (D) The Module-trait heatmap displays the correlation of each gene module with AD and normal states. (E) The Venn diagram demonstrates the interaction between DEGs and genes in the MEbrown module, with 223 genes shared.

The clustering dendrogram reveals the hierarchical clustering of genes into 11 distinct modules (Figure 3C). Among these, the MEbrown module exhibited a strong positive correlation with AD (correlation coefficient = 0.77), while the MEgrey module displayed a strong negative correlation (correlation coefficient = −0.74), both with high statistical significance (p-value = 1e-07). These results, visualized in the heatmap (Figure 3D), highlight the correlations between gene modules and traits of normal and AD-affected tissues. The strong association of MEbrown and MEgrey with AD suggests that the genes within these modules may play essential roles in AD pathogenesis.

The MEbrown module, containing 749 genes, was identified as the primary module associated with AD. Further analysis of the overlap between DEGs and genes in MEbrown revealed 223 intersecting genes, as shown in Figure 3E. These intersecting genes are potential candidates for further investigation into their specific functions and mechanisms in AD, advancing our understanding of disease etiology and identifying potential therapeutic targets.

Selection of Feature Genes Through LASSO and SVM-RFE Algorithms

Two machine learning algorithms, LASSO regression and SVM-RFE, were employed to identify key feature genes associated with atopic dermatitis from the pool of candidate genes. In the LASSO analysis, 10 feature genes were selected, as shown in Figure 4A and B. The SVM-RFE analysis identified 28 feature genes with a relative importance score above 0.25 (Figure 4C). The intersection of genes selected by both algorithms yielded 5 core feature genes: ATP6V1A, CLDN23, ECSIT, LRFN5, and USP16 (Figure 4D). These five genes represent potential biomarkers or therapeutic targets due to their consistent identification as key features by both methods, indicating their relevance to atopic dermatitis pathogenesis.

Figure 4 Machine Learning Algorithms for Feature Genes. (A) The penalty plot of the LASSO model, with error bars representing standard error. (B) The LASSO plot shows the shrinkage of parameter coefficient sizes as the k penalty value increases. (C) The error rate confidence interval of the random forest model. (D) The interaction between the LASSO and random forest algorithms.

Diagnostic Performance of Feature Genes in Predicting Atopic Dermatitis

In this study, we identified significant differences in the expression of key genes by comparing the gene expression profiles of the control and AD groups. The expression levels and statistical analysis results of key genes are shown in Table 1. Specifically, the expression levels of LRFN5, CLDN23, ATP6V1A, and USP16 were significantly lower in the AD group compared to the control group, whereas ECSIT exhibited significantly higher expression in the AD group by statistical analyses with P-values < 0.05 (Figure 5A–E). To assess the diagnostic potential of these genes, Receiver Operating Characteristic (ROC) curves were generated. The Area Under the Curve (AUC) values, which range from 0.5 to 1.0, reflect the diagnostic accuracy of each gene, with higher AUC values indicating better diagnostic performance. The ROC analysis of the five core feature genes showed AUC values greater than 0.8, indicating high accuracy in distinguishing AD from control samples: LRFN5 (AUC=0.825), ECSIT (AUC=0.880), CLDN23 (AUC=0.810), ATP6V1A (AUC=0.925), and USP16 (AUC=0.895) (Figure 5F–J). These findings suggest that ATP6V1A, CLDN23, ECSIT, LRFN5, and USP16 could serve as valuable biomarkers for diagnosing atopic dermatitis and may also be promising therapeutic targets for future treatments.

Table 1 Expression Levels and Statistical Analysis Results of Key Genes in Control and AD Groups

Figure 5 Diagnostic performance of Key Genes. (A–E) The expression of the feature genes LRFN5, CLDN23, ATP6V1A, USP16 and ECSIT in atopic dermatitis and healthy cohorts. (F–J) ROC curves of Key Genes.

Immune Cell Infiltration in Atopic Dermatitis

Recent studies highlight the significant role of immune cell infiltration in the pathogenesis of atopic dermatitis.^23–25 We analyzed immune cell profiles between AD and control groups. Our findings revealed an increased proportion of activated dendritic cells and neutrophils in the AD group (Figure 6A), suggesting heightened immune activation and a potential role in disease progression. These shifts in immune cell populations point to the critical involvement of immune surveillance and inflammatory responses in AD.

Figure 6 Immune Infiltration Analysis. (A) The comparison of the number of different immune cell types between the control group (blue bars) and the AD group (yellow bars). (B) The heatmap showing the correlation between ATP6V1A, CLDN23, ECSIT, LRFN5 and USP16 genes and various immune cell types. (C) The scatter plot displaying the expression levels of the ECSIT gene in different immune cell types (* p < 0.05, **P<0.01). Abbreviations: NK cells, Natural Killer cells; Macrophages M0/M1/M2, Non-activated/Classically activated/Alternatively activated macrophages; aDCs, Activated dendritic cells; iDCs, Immature dendritic cells; pDCs, Plasmacytoid dendritic cells; APC, Antigen-presenting cell; CCR, Chemokine receptor; HLA, Human leukocyte antigen; MHC, Major histocompatibility complex; Tfh, Follicular helper T cells; Th1/Th2, Type 1/Type 2 helper T cells; TIL, Tumor-infiltrating lymphocytes; IFN, Interferon.

Further investigation into the correlation between the five key genes (ECSIT, LRFN5, USP16, ATP6V1A, and CLDN23) and various immune cell types revealed distinct patterns. ATP6V1A was positively correlated with CD4 memory-activated cells and neutrophils, while CLDN23 and ECSIT exhibited positive correlations with naive B cells and helper T cells. LRFN5 and USP16 showed unique relationships across multiple T cell and myeloid cell subtypes (Figure 6B).

Additionally, we observed the impact of ECSIT expression on immune cell subtype activities. In samples with elevated ECSIT expression, there was a significant increase in glycolytic cells, M1 macrophages, pro-inflammatory cells, monocytes, and neutrophils, with the most pronounced differences in M1 macrophages and pro-inflammatory cells (Figure 6C). These results suggest that higher ECSIT expression may exacerbate inflammation by promoting the activation of specific immune cell populations, which could contribute to the chronic inflammation seen in atopic dermatitis.

RT-PCR Validation of Key Genes

We performed RT-qPCR to validate the mRNA expression of core genes identified through bioinformatics analysis in blood samples from patients with and without atopic dermatitis. All experimental parameters, from primer design (Table 2) to reaction optimization (Table 3), have been systematically documented. The results demonstrated significant differences in the expression levels of LRFN5, ECSIT, and USP16 between the AD and control groups (Figure 7A–C), aligning with the findings from microarray data. However, no statistically significant differences were found for ATP6V1A and CLDN23 between the two groups (Figure 7D and E).

Table 2 Primer Sequences Used for RT-qPCR Validation of Candidate Genes in AD

Table 3 Optimization Procedures for RT-qPCR Gene Analysis

Figure 7 RT-PCR validation of the key gene relative mRNA expression between AD and normal controls. All experiments were performed in triplicate and results were presented as M ± SD (* p < 0.05, ***P<0.001). (A–C) ECSIT gene was significantly higher in AD group, while LRFN5 and USP16 genes were lower in AD group than in healthy group in the blood samples. (D-E) ATP6V1A and CLDN23 gene in AD patients had no significant difference compared with the healthy group in the blood samples. Abbreviation: ns, non significant.

Causal Relationship Between Atopic Dermatitis and Immune Cells

Atopic dermatitis is a complex chronic inflammatory skin condition influenced by genetic and environmental factors, with immune cells playing a pivotal role. Using Mendelian randomization analysis, we identified several immune cells and their markers significantly associated with an increased risk of AD. For example, elevated percentages of CD25hi %T cells, basophil %CD33dim HLA DR− CD66b−, naive DN (CD4−CD8−) AC, and CD40 expression on monocytes were all linked to an increased risk of AD. In contrast, a decreased risk of AD was associated with CD28− DN (CD4−CD8−) %DN (Figure 8). These findings suggest critical immune cell-based mechanisms underlying AD pathogenesis, providing potential targets for diagnosis and treatment of the disease.

Figure 8 Forest Plot Showing the Causal Relationship between Immune Cell Characteristics and Atopic Dermatitis. Abbreviations: nsnp, Number of Single Nucleotide Polymorphisms; pval, P-value; OR, Odds Ratio; CI, Confidence Interval; hi/dim, High/diminished expression; HLA DR, Human Leukocyte Antigen DR; DN, Double negative; AC, Absolute count; CCR2, C-C chemokine receptor type 2.

ceRNA Network

In the ceRNA network, non-coding RNAs such as miRNAs and lncRNAs interact to regulate gene expression through molecular competition. ECSIT emerged as a central player in the AD-related ceRNA network. Specifically, ECSIT interacts with miRNAs such as hsa-miR-1181 and hsa-miR-1205, which in turn connect to lncRNAs like RP11-469N6.1 and RP3-470B24.5 (Figure 9). This complex interaction suggests that ECSIT might regulate immune responses and inflammatory processes by modulating the activity of miRNAs and lncRNAs, which can control the expression of inflammation-related genes. Furthermore, lncRNAs may act as “molecular sponges”, absorbing miRNAs to prevent the suppression of other mRNA targets, thus refining the network of gene regulation in AD.

Figure 9 ceRNA Network of ECSIT Gene-Associated lncRNAs and miRNAs.

Discussion

Atopic dermatitis is a chronic inflammatory skin disease that severely impacts the quality of life for millions of individuals globally.^26,27 It typically manifests in childhood and may persist into adulthood, causing symptoms such as intense itching, skin erythema, and inflammation. The pathogenesis of AD is multifactorial, involving genetic predisposition, environmental triggers, immune dysregulation, and skin barrier dysfunction. While the exact cause remains unclear, there is significant evidence that immune responses, particularly involving T cells, and interactions with the microbiome and environmental factors, play crucial roles in AD development and progression.^21,22,27

To investigate the underlying mechanisms of AD, we identified DEGs associated with AD and performed enrichment analysis, which highlighted key biological processes and signaling pathways involved in the disease, such as skin barrier dysfunction, immune responses, and lipid metabolism. GO analysis underscored biological processes related to skin barrier establishment, epidermal cell differentiation, and cornified envelope organization—findings that align with AD’s hallmark compromised barrier function. Significant enrichment in cellular components like lysosomes and endocytic vesicles suggests an essential role of intracellular transport and degradation in AD pathophysiology, while molecular functions such as growth factor receptor binding and copper ion binding imply a potential regulatory role of cell signaling and metal ions in inflammation. KEGG pathway analysis further revealed significant alterations in pathways central to AD, including sphingolipid and arachidonic acid metabolism, which are integral to lipid metabolism and the generation of inflammatory mediators.²⁸ The observed changes in cell cycle and NOD-like receptor signaling pathways further underscore a disrupted balance of cell proliferation and immune response regulation in AD.²⁹ Additionally, significant enrichment of viral carcinogenesis and hepatitis C pathways might reflect the high sensitivity of AD patients to external infections and corresponding immune responses.³⁰

Using WGCNA and LASSO regression, to pinpoint critical genes related to AD, we identified five genes—ATP6V1A, CLDN23, ECSIT, LRFN5, and USP16—most strongly associated with AD. Notably, ECSIT was upregulated, while LRFN5 and USP16 were downregulated, a pattern confirmed through RT-qPCR in blood samples from AD patients. In contrast, ATP6V1A and CLDN23 showed no significant expression changes, suggesting they may be less central in AD pathogenesis. ROC curve analysis of these genes indicated high diagnostic value (AUC > 0.8), highlighting their potential as biomarkers.

Immune cell infiltration and activation are critical events in the pathology of AD, profoundly affecting the disease phenotype and clinical progression.³¹ Notably, the significant increase in activated dendritic cells and neutrophils in the AD group underscored the direct impact of AD on the activity of key cells in the immune system. This change not only reflected enhanced immune response but also likely influenced the regulation of the inflammatory response.³² Specifically, the activation state of dendritic cells suggested their central role in antigen presentation and initiating T-cell-mediated immune responses, while neutrophil increases align with the inflammatory nature of acute AD.³³ Furthermore, correlation analysis between ECSIT, LRFN5, Uarious immune cells suggested that these genes influence immune cell function in distinct ways. ECSIT, for instance, was significantly associated with M1 macrophages and pro-inflammatory cells, suggesting a pro-inflammatory role that may worsen AD pathology.^34,35 In contrast, antigen-presenting cell subgroups did not show significant expression differences, suggesting that more targeted strategies might be needed to regulate this process in AD treatment.³⁶ Additionally, the significant activation of glycolytic cells and M1 macrophages in relation to high ECSIT expression highlighted the potential role of metabolic pathways in regulating immune cell functions.³⁷ This metabolic reprogramming might be a way for inflammatory cells to adapt to the rapid response demands, providing a biological basis for disease chronicity and tissue damage.³⁸ LRFN5 and USP16 exhibited unique correlation patterns across multiple T cell and myeloid cell subtypes. LRFN5 is involved in nerve and immune cell signaling, which may be activated in chronic inflammatory states, potentially affecting immune cell responses to inflammation.³⁹ USP16 is involved in the activation and proliferation of T cells by regulating ubiquitination levels, especially deubiquitination of calcineurin A (CNA). In the case of an abnormal immune system, loss of function of USP16 may lead to limited proliferation of T cells, making them less responsive to inflammatory diseases.⁴⁰ The key roles of USP16 in immunomodulatory and inflammatory signaling pathways suggest that it may influence the pathophysiological mechanisms of AD, especially in regulating skin barrier function or inflammatory response.

The central role of ECSIT in AD was further supported by our construction of a ceRNA network, where ECSIT was found to interact with various miRNAs and long non-coding RNAs. These interactions suggest that ECSIT may influence gene expression regulation in immune responses, contributing to AD pathogenesis. Specifically, the involvement of ECSIT in the activation of M1 macrophages and pro-inflammatory cells suggests that it might exacerbate inflammation, while the regulatory interactions with miRNAs and lncRNAs highlight the complexity of gene expression regulation in AD.

Through Mendelian randomization, we further explored the causal relationships between immune cell markers and AD. An increase in CD28− DN (CD4−CD8−) %DN was associated with a reduced AD risk, indicating that this subgroup might play a protective role in regulating immune responses, maintaining immune tolerance, or skin barrier function.⁴¹ This finding offers new insights into the immune regulation of AD and suggests potential therapeutic targets.^42,43 These results provide important clues for the immunopathological mechanisms of AD and highlight the significance of intervention strategies targeting specific immune cells and molecular markers for AD prevention and treatment. Further research should aim to elucidate the specific mechanisms of these immune markers in AD and how they interact with other genetic and environmental factors to lead to disease onset and progression.

Notably, while our study focused on intrinsic genetic and immune drivers of AD, environmental factors (eg, microbial colonization, pollution) likely modulate the activity of key genes. For instance, mitochondrial dysfunction mediated by ECSIT may exacerbate inflammation under high oxidative stress conditions or pollutant exposure,⁴⁴ while psychosocial stress could impair LRFN5-dependent neuro-immune homeostasis.⁴⁵ Our study identifies USP16 as a critical regulator of T-cell activation, its interaction with microbiome-derived metabolites (eg, butyrate) remains unexplored. Recent evidence shows that butyrate enhances deubiquitinase activity, suggesting a potential mechanism by which the microbiome modulates USP16 function.⁴⁶ Future studies integrating exposome data and multi-omics profiling are needed to dissect these interactions. Despite this limitation, our findings provide a foundational framework for understanding AD pathogenesis, independent of environmental confounders.

In terms of treatment, immune-targeted therapies, such as biologics (eg, dupilumab) and JAK inhibitors, have revolutionized AD management by targeting key cytokine pathways involved in the disease.^47–50 However, as AD is a heterogeneous disease, individualized treatment approaches are crucial, with ongoing efforts to identify additional molecular targets and combination therapies to improve outcomes. Notably, our identification of ECSIT as a novel pro-inflammatory regulator in AD aligns with emerging evidence on mitochondrial dysfunction in chronic inflammation.⁵¹

While our study offers novel insights into the pathophysiology of atopic dermatitis and identifies potential diagnostic biomarkers, we recognize the necessity for future validation in larger cohorts. Expanding the sample size in future research would enhance the statistical robustness and generalizability of our results. To address this, we prioritized multicenter collaboration to expand our sample size and plan to incorporate additional independent datasets in subsequent studies to validate both the expression patterns and diagnostic efficacy of the key candidate genes. Furthermore, stratified studies among populations with different factors such as gender and age are crucial for assessing the consistency and applicability of our research findings across diverse demographic groups.

Additionally, while our findings highlight specific immune cell subsets and genes that may play pivotal roles in atopic dermatitis, further investigation is required to elucidate their precise mechanisms in disease progression. The potential involvement of immune cell infiltration in central genes also merits deeper exploration. Future planned experiments include functional analyses (eg, gene knockdown/overexpression in cell cultures) and in vivo experiments (AD Mouse models) to further clarify the roles and mechanisms of these genes. Addressing these limitations will strengthen the clinical relevance and translational potential of the biomarkers and therapeutic targets identified in this study.

Conclusion

In conclusion, this study identifies ECSIT, LRFN5, and USP16 as pivotal biomarkers linking immune dysregulation to AD pathogenesis. The upregulation of ECSIT and its association with M1 macrophages underscore its role as a potential therapeutic target for chronic inflammation. While our findings are mechanistically robust, future studies must validate these genes in diverse populations and integrate environmental exposome data to achieve personalized AD management.