Introduction
Atopic dermatitis (AD) is a common, chronic, relapsing inflammatory skin disease that severely impairs patients’ quality of life, with global prevalence rising annually.1,2 Clinically, AD is marked by intense pruritus, eczematous lesions, barrier dysfunction, and chronic inflammation.3 Immunologically, a hallmark of AD is aberrant activation and infiltration of CD4⁺ T cells. These cells secrete various cytokines and interact with other cell types to orchestrate the initiation, maintenance, and amplification of inflammation, creating a highly dynamic and complex immune network.4,5
Immune-mediated mechanisms have been recognized as central to the pathogenesis of various skin diseases.6 These include psoriasis,7 atopic dermatitis,8 vitiligo,9 and cutaneous lupus erythematosus,10 all characterized by aberrant immune responses against skin tissues, resulting in chronic inflammation and barrier dysfunction.11 T cells have been shown to play pivotal roles in the initiation and maintenance of these inflammatory responses,12 yet traditional approaches often treat immune states as static, thereby missing the temporal regulation of disease-associated genes.
Mendelian randomization (MR) offers an efficient and economical strategy that has proven useful for identifying and prioritizing potential therapeutic targets.13 An essential tool for the functional annotation of genetic variation in non-coding regions is Expression Quantitative Trait Loci (eQTL) analysis, which establishes a causal connection between genetic variation and levels of gene expression.14,15 While static eQTL affect generic cellular processes, dymatic eQTL are more often involved in hematopoiesis and immune response.16,17 Combining dynamic eQTL data with MR may enable a broader evaluation of how gene expression in activated T-cell subsets influences AD risk. CD4⁺ T cells are central regulators of adaptive immunity and play a pivotal role in the pathogenesis of AD. Aberrant activation and differentiation of CD4⁺ T-cell subsets, particularly Th2 cells, promote the production of IL-4, IL-13, and IL-31, which contribute to skin barrier dysfunction, chronic inflammation, and pruritus.14,18,19 Therefore, understanding the molecular mechanisms governing CD4⁺ T-cell activation is essential for elucidating AD pathogenesis.
Nathan et al revealed through single-cell eQTL analysis that dynamic changes in T cell functional states critically influence the regulation of disease-associated genetic loci, demonstrating that conventional bulk analyses may overlook cell state-specific gene regulatory mechanisms.20
We first applied two-sample Mendelian randomization (MR) to infer causal relationships between dynamic eQTLs and AD. Next, we employed colocalization analysis, single-cell transcriptomic trajectories, protein-protein interaction (PPI) networks to functionally interpret the regulatory impact of prioritized genes on skin barrier, immune signaling, and translational pathways. Finally, to bridge the gap between genetic causal inference and functional validation, we utilized a single-cell virtual knockout algorithm.21
The novelty of our study lies in: (i) focusing on dynamic eQTLs across multiple activation timepoints in CD4⁺ T cells; (ii) being the first to integrate these dynamic eQTLs with AD GWAS using two-sample Mendelian randomization and colocalization to generate genetically-informed causal prioritization; and (iii) integrating single-cell simulation knockout with pseudotime, cell–cell communication, and drug/network mapping to provide activation- and cell-context dependent mechanistic hypotheses. We have not only integrated genetic causal inference with single-cell functional annotations, but also bridged a critical knowledge gap by linking dynamic eQTLs to AD pathogenesis. Our aim is to select those genes that are highly expressed and carry high risk under specific activation conditions, with the goal of uncovering potential therapeutic targets that may be more personalized and temporally specific.
Method
Acquisition of Dynamic eQTL for CD4 ⁺ T Cell Genes
We obtained dynamic eQTL data from a single-cell RNA-seq GWAS of 119 European-ancestry donors under resting and PHA-activated conditions, identifying 482,971 cis-eQTL (±1 Mb).14 To select conditionally independent cis-eQTL, we performed LD clumping in PLINK 1.9 (1000 Genomes EUR reference; r2 < 0.001; 10 Mb window), retaining the most significant SNP per gene. We then computed F-statistics (F = β2/SE2) and selected SNPs with F > 10. Finally, we applied the Steiger directionality test to remove potential reverse-causation instruments.
Acquisition and Processing of Non Dynamic eQTL Data
To ensure consistency, we applied identical filtering criteria to static eQTL data. We extracted cis-eQTL from resting-state OneK1K data across 14 immune subsets, whole-blood eQTLGen (>30,000 samples), and immune-cell-specific OneK1K datasets. All instruments underwent the same LD clumping (r2 < 0.001), and F > 10 selection.
Selection of Outcome Traits for AD
The outcome phenotype was derived from FinnGen R12 AD GWAS summary statistics, comprising 31,245 cases and 432,874 controls.
MR Analysis of Dynamic eQTL and Non Dynamic eQTL Dataset
For each dynamic and static eQTL SNP, causal effects on AD were estimated using the single-SNP Wald ratio;22,23 for multi-SNP instruments, inverse-variance weighted (IVW) meta-analysis was used.24 P-values were Benjamini–Hochberg FDR-adjusted, with FDR < 0.05 deemed significant.25 Analyses were conducted in R 4.4.1 with the TwoSampleMR package. And all the results underwent reverse Steiger testing.
Colocalization Analysis
We conducted colocalization analysis to evaluate whether genetic variants underlying gene expression quantitative trait loci (eQTLs) and disease risk share the same causal signal. Specifically, we applied the R package coloc to perform Bayesian colocalization tests between eQTL summary statistics and GWAS summary statistics. For each genomic region, we extracted and harmonized summary-level association data, including effect sizes, standard errors, minor allele frequencies, and sample sizes, to ensure consistency across datasets.
The coloc.abf function was then used to compute posterior probabilities for five hypotheses: no association with either trait (H0), association with only trait 1 (H1), association with only trait 2 (H2), association with both traits but driven by different causal variants (H3), and association with both traits driven by the same causal variant (H4). The prior probabilities adopted in the analysis were p1 = 1 × 10−4, p2 = 1 × 10−4, and p12 = 1×10−5. P1 represents the prior probability that a randomly selected genetic variant is associated with gene expression, p2 represents the prior probability that the variant is associated with AD, and p12 represents the prior probability that the variant is associated with both gene expression and AD.
Evidence of colocalization was defined as a posterior probability for H4 (PP.H4) greater than 0.9, a threshold that has been widely adopted in previous studies.26
Identification of Immune Cell Specific Expression of Genes Related to AD
We used TPM/CPM data from the OneK1K consortium for five major immune-cell types (CD4⁺ T, CD8⁺ T, B, NK, monocytes) and their subsets to compare candidate gene expression.15 Overall differences were assessed by Kruskal–Wallis tests, with Dunn’s post hoc tests for multiple comparisons, identifying genes enriched in specific cell populations.
Analysis of Single-Cell Sequencing Data in AD Group
Raw data were obtained from the GEO dataset GSE153760 and imported into Seurat for quality control: cells with nFeature_RNA < 4500, nCount_RNA < 25,000, and percent.mt < 30% were retained; data were normalized using the LogNormalize method and highly variable genes were identified.27 Cell clustering was performed via principal component analysis (PCA) followed by Louvain clustering at a resolution of 0.5, and results were visualized using UMAP.28 Pseudotime trajectories were reconstructed with Monocle 3, and CellChat was employed to predict intercellular ligand–receptor interaction networks.
Identification of Drug Targets and Construction and Analysis of Protein-Protein Interaction (PPI) Networks
We retrieved drug target information from DrugBank (accessed 20 April 2025) using keyword searches.29 We then input these drug targets alongside candidate genes (e. g., IL4, IL13, AHI1, IRF5) into STRING (Homo sapiens) to construct a PPI network, revealing key protein interactions.
Single-Cell Virtual Knockout Analysis
Single-cell virtual knockout analysis was performed using the R package scTenifoldKnk.21 Single-cell data were extracted from the Seurat object and converted into a raw count matrix for downstream analysis. Subsequently, virtual knockout of genes such as AHI1 (specified as gKO = “AHI1”) were conducted using the scTenifoldKnk function. This method constructs and compares gene regulatory networks based on a manifold alignment framework to assess the impact of transcriptional network perturbation.
In this analysis, 10 subnetworks were generated (nc_nNet = 10), each randomly sampling 500 cells (nc_nCells = 500), and three principal components (nc_nComp = 3) were used during network dimensionality reduction to capture the primary sources of variation. Following this, differential regulation scores (fold change, FC) were computed and ranked for each gene, and the top 20 most significantly affected genes were visualized using ggplot2. The important point is that these scores measure the disruption of the regulatory network (as the role of highly sensitive nodes or critical regulatory hubs), rather than the absolute or relative changes in transcript abundance.
Genes with an adjusted p-value < 0.05 were considered significantly dysregulated following virtual knockout. The final results were exported as CSV files and plots for subsequent analysis and interpretation.
Software Environment and Repeatability
All statistical analyses were conducted in R 4.4.1, primarily using the TwoSampleMR, COLOC, Seurat, Monocle, and CellChat packages.
Result
Selection of Instrumental Variables
Figure 1 provides a summary of the study design. We constructed MR instruments from dynamic eQTL data. The eQTLs were obtained from single-cell transcriptomic analysis of CD4+ T cells at five different time points, transitioning from a resting state to the fully activated state.
The time points include the resting state (0 hours), low-activity state following activation (LA), 16 hours post-activation (16 hours, prior to cell division), 40 hours post-activation (40 hours, after the first cell division), and 5 days post-activation (5 days, during the acquisition of effector functions). From single-cell transcriptomic profiles of CD4⁺ T-cell activation across 46 cell-type/time-point combinations (eg, CD4_Memory_5d), we initially identified 482,971 cis-eQTL involving 13,379 genes. A three-stage filtering pipeline—LD clumping (r2 < 0.001, 10 Mb window), F-statistic > 10, and Steiger directionality testing—yielded 9644 independent SNPs representing 1538 genes.Figure 1Workflow of dynamic CD4⁺ T-cell eQTL analysis and causal inference for atopic dermatitis.
This diagram illustrates the overall analytical pipeline of the study.A flowchart of analyses linking dynamic eQTLs to atopic dermatitis and downstream single cell steps.
To assess the unique contribution of dynamic eQTL, we cross-compared our significant genes with static eQTL data from OneK1K (n = 982) and eQTLGen (n = 1418–31,684). Applying the same filtering steps produced 9723 and 65,904 static eQTL instruments, respectively, highlighting regulatory variants specific to activation conditions.
The Causal Relationship Between Dynamic eQTL of CD4+T Cell Activation and AD
Two-sample MR with FDR correction identified 298 positive associations, encompassing 59 dynamic eQTL across 47 activation conditions in CD4⁺ T cells. FADS1 increased AD risk in naïve CD4⁺ T cells at 16 h (OR = 1.042; 95% CI: 1.025–1.059; FDR < 0.001), 40 h (OR = 1.053; 95% CI: 1.031–1.076; FDR < 0.001), and 5 d (OR = 1.059; 95% CI: 1.028–1.091; FDR = 0.011). Its paralog FADS2 showed strong associations in memory T cells at 40 h (OR = 1.040; 95% CI: 1.019–1.063; FDR = 0.012) and in naïve T cells at 5 d (OR = 1.069; 95% CI: 1.032–1.108; FDR = 0.012).
IRF5 was significantly associated after 40 h activation in memory (OR = 1.050; 95% CI: 1.027–1.072; FDR = 0.002) and TEM (OR = 1.062; 95% CI: 1.034–1.091; FDR = 0.002) cells. Conversely, IL1RN exhibited a protective effect in TCM cells at 5 d (OR = 0.927; 95% CI: 0.896–0.960; FDR = 0.003). Notably, GNGT2 showed the strongest association in long-term activated naïve T cells at 5 d (OR = 1.190; 95% CI: 1.143–1.239; FDR < 0.001), and RPS26 was broadly significant across TEMRA, TCM, and TN subsets (eg, TEMRA 5 d: OR = 1.041; 95% CI: 1.017–1.065; FDR = 0.029).
HDAC3 in resting TN cells reduced AD risk (OR = 0.948; 95% CI: 0.917–0.979; FDR = 0.037), while CTSS was protective in memory cells at 5 d (OR = 0.888; 95% CI: 0.837–0.941; FDR = 0.006) (Supplementary Figures 1–5 and Supplementary Tables 1, 2). We identified multiple genes with consistently directional effects in response to cellular activation states. Figure 2A displays the top genes ranked by absolute effect size.
Steiger tests confirmed no reverse causality (Supplementary Table 3). In total, after applying FDR correction (q < 0.05), 298 genes demonstrated a significant causal association with AD risk. These findings highlight potential therapeutic targets within dynamic CD4⁺ T-cell regulatory networks.Figure 2Causal inference and cross platform validation of CD4 ⁺ T cell eQTL in atopic dermatitis.
(A) This lollipop plot illustrates the top target genes with the largest absolute effect sizes (Beta). We identified multiple genes with consistently directional effects in response to cellular activation states: COX6B1, GNGT2, PPP5C, UBE2Z, ARL6IP4, STMN3, PPP2R3C, NDUFS3, SLC4A7, and KDELR2 genes were associated with an increased risk of AD, whereas MPHOSPH9, CTSS, TNFRSF14, NDUFAF1, and ME2 genes demonstrated a protective effect against AD. (B) The UpSet plot compared AD related genes identified using dynamic eQTL and two static eQTL datasets (OneK1K and eQTLGen).
The dynamic eQTL uniquely contributes a large number of genes related to AD. (C) The volcano plot displays the Mendelian randomization results of plasma eQTL. Each point represents a gene, with the x-axis representing effect size (Beta) and the y-axis representing significance(-log10 (FDR)).
Genes that significantly upregulate or downregulate the risk of AD are marked in red and blue, respectively. (D) The volcano plot displays the Mendelian randomization results of non dynamic single-cell eQTLs in OneK1K. Each point represents a gene, with the x-axis representing effect size (Beta) and the y-axis representing significance (-log10 (FDR)).
Genes that significantly upregulate or downregulate the risk of AD are marked in red and blue, respectively.Four plots comparing gene associations: a lollipop chart, an UpSet plot and two volcano plots.A) Lollipop chart titled, Top Genes by Absolute Effect Size. X axis label, Effect Size left parenthesis Beta right parenthesis, unit not shown, ticks at negative 0.2, 0.0, 0.2, 0.4. Y axis label, Gene Name, unit not shown.
Gene labels: COX6B1, GNGT2, PPP5C, UBE2Z, ARL6IP4, STMN3, PPP2R3C, NDUFS3, SLC4A7, KDELR2, MPHOSPH9, CTSS, TNFRSF14, NDUFAF1, ME2. Legend text includes, Effect Direction, Negative Effect, Positive Effect and minus log base 10 left parenthesis FDR right parenthesis with dot sizes 2.5, 5.0, 7.5, 10.0, 12.5, unit not shown. B) UpSet plot.
Upper bar chart y axis label, Intersection Size, unit not shown, with labeled bars 149, 127, 25, 21, 19, 9, 1. Left set labels: Dynamic, OneK1K, eQTLGen. Dot matrix shows intersections with vertical connectors.
C) Volcano scatter plot. X axis label, Beta, unit not shown, range about negative 1.0 to 1.0. Y axis label, minus Log base 10 left parenthesis FDR right parenthesis, unit not shown, range 0 to about 18.
The Causal Relationship Between Non Dynamic eQTL and AD
MR using static eQTL data identified nine genes with causal links to AD shared between dynamic and static CD4⁺ T-cell instruments (Figure 2B). In plasma eQTL data, GNGT2 (OR = 1.17; 95% CI: 1.11–1.23; FDR < 0.001) and NMB (OR = 1.175; 95% CI: 1.10–1.234; FDR < 0.001) showed the strongest positive associations. MPHOSPH9 (OR = 0.844; 95% CI: 0.785–0.908; FDR = 0.00156) and KDELR2 (OR = 1.054; 95% CI: 1.030–1.078; FDR = 0.001) displayed significant negative and positive effects, respectively (Figure 2C and Supplementary Figures 6–8, Supplementary Table 4).
Steiger tests excluded reverse causation (Supplementary Table 5). Cell-type–specific analyses revealed that MPHOSPH9 increased risk in NK cells (OR = 1.137; 95% CI: 1.083–1.194; FDR < 0.001), CD8et (OR = 1.177; 95% CI: 1.099–1.260; FDR < 0.001), and CD8nc (OR = 1.128; 95% CI: 1.067–1.191; FDR = 0.001). AHI1 was uniformly risk-increasing across immune subsets (eg, CD8s100b: OR = 1.164; 95% CI: 1.096–1.235; FDR < 0.001), suggesting broad activation functions.
FADS1 exhibited cell-type duality—protective in CD4et (OR = 0.901; 95% CI: 0.863–0.939; FDR < 0.001) and NK cells (OR = 0.863; 95% CI: 0.813–0.916; FDR < 0.001) but risk-increasing in CD8nc (OR = 1.109; FDR = 0.027). KDELR2 was protective globally (FDR = 6.42×10−7). Remarkably, GNGT2 showed the strongest protective effect in NK cells (OR = 0.785; 95% CI: 0.742–0.830; FDR < 0.001).
PPP5C increased risk in CD8et (OR = 1.183; 95% CI: 1.106–1.265; FDR < 0.001) but was protective in NK cells (OR = 0.922; 95% CI: 0.891–0.954; FDR < 0.001), demonstrating bidirectional regulation. TMEM258 (CD8 effector: OR = 0.903; 95% CI: 0.858–0.948; FDR = 0.003; CD4nc: OR = 0.956; 95% CI: 0.929–0.982; FDR = 0.034), CDK5R1 (CD8nc: OR = 1.121; 95% CI: 1.048–1.198; FDR = 0.025), and NMB (CD4nc: OR = 1.151; 95% CI: 1.078–1.228; FDR = 0.002) further highlighted cell-type heterogeneity in AD risk (Figure 2D and Supplementary Figures 9–13, Supplementary Tables 6, 7).
Colocalization Results of Eqtl and AD in CD4+T Cells
Our colocalization analysis (PPH4 > 0.9) revealed 11 genes with shared genetic signals between dynamic CD4⁺ T-cell eQTL and AD. FADS1 showed high colocalization in memory cells at 16 h (PPH4 = 0.989), 40 h (0.988), and 5 d (0.988), as well as in naïve cells at 16 h (0.989) and 40 h (0.989). AHI1 colocalized in unstimulated memory cells (PPH4 = 0.980) and in naïve cells at 16 h (0.982) and 5 d (0.959).
KDELR2 peaked in TEM cells at 40 h (PPH4 = 0.999) and was also significant in memory cells at 5 d (0.961) and naïve cells at 16 h (0.976). IRF5 colocalized in memory cells at 40 h (PPH4 = 0.955) and TEM at 40 h (0.975). MPHOSPH9 displayed colocalization in unstimulated memory (0.926), naïve (0.914), and TN at 0 h (0.919).
DCTN5 was stably colocalized in memory cells at 16 h, 40 h, and 5 d (all 0.947). CTSS was significant at 5 d in memory cells (0.924). ARL6IP4 colocalized at 5 d in memory (0.989) and TEM (0.989).
JAZF1 showed colocalization in naïve cells at 16 h (0.975), unstimulated (0.969), and TN at 0 h (0.936). PPP5C colocalized at 40 h in naïve cells (0.989) and at 16 h in TCM (0.964). PTGER4 was specifically significant at 16 h in TN cells (0.986) (Supplementary Table 8).
Dynamic CD4+T Cell eQTL Data Reveals Causal Effects of Gene Expression on AD During Activation
By integrating dynamic CD4⁺ T-cell eQTL with colocalization results, we systematically assessed the causal effects of gene expression in resting versus activated states on AD. Among genes with strong colocalization evidence, 4 (36.4%) showed significant causal signals at rest, whereas all 11 (100%) exhibited activation-dependent causal effects (some genes showing signals in both states) (Figure 3A). Stratification by activation time revealed distinct causal genes at 0 h and LA, underscoring the importance of immune-cell activation context in dissecting AD genetics; genes at 16 h, 40 h, and 5 d largely overlapped with other time points (Figure 3B).
Figure 3C demonstrated differential gene expression among immune cell types, with some genes exhibiting higher expression in CD4+ T cells, CD8+ T cells, and monocytes. The FUMA-based tissue-specific heatmap (Figure 3D) further showed that these key genes are actively expressed in brain tissues, peripheral blood, and immune-related organs, with ARL6IP4, PPP5C, and KDELR2 showing elevated expression in various tissues, suggesting that these genes may play important regulatory roles in immune responses and tissue-specific functions.Figure 3Causal effects of dynamic and static eQTL on AD and expression differences of colocalized genes across cells and tissues. (A) Venn plot comparing the intersection of static and dynamic eQTLs with causal effects on AD; (B) Network plots showing common and unique relationships across five time points of dynamic eQTLs with causal effects on AD; (C) The Onek1k database displays the expression of genes with strong colocalization evidence in different cell types, with a color gradient from blue (low expression) to red (high expression); (D) FUMA database shows the total expression of genes with evidence of colocalization in different tissues.Four-part infographic with a Venn diagram, a timepoint network and two gene expression heatmaps.A) Venn diagram titled, Resting vs Activation Gene Intersection.
Circles labeled, Resting State and Activation status. Values: Resting State only 0, overlap 4, Activation status only 7. B) Network plot with nodes connected by curved lines.
Labels near the network read, 0h, 16h, 40h, 5d and LA. C) Clustered heatmap titled, Marker Gene Expression Across Cell Types. Row labels: PPP5C, MPP5H9, DCTN5, AHI1, ARL6IP4, PTGER4, JAZF1, KDELR2, FADS1, IRF5, CTSS.
Column labels: Monocyte, CD4 T cell, CD8 T cell, NK cell, B cell. Scale labeled from minus 1.5 to 1.5. D) Heatmap with row labels: AHI1, ARL6IP4, CTSS, DCTN5, FADS1, IRF5, JAZF1, KDELR2, MPP5H9, PPP5C, PTGER4.
Tissue columns include Adipose Subcutaneous through Whole Blood, including Brain Spinal cord cervical c 1. Scale labeled 0, 2.84, 5.67.
The Expression of Identified Genes in Single-Cell Sequencing
This study analyzed the expression patterns of 11 genes across a variety of human cell types, including monocytes, CD8⁺ T cells, CD4⁺ T cells, epithelial cells, keratinocytes, induced pluripotent stem (iPS) cells, mesenchymal stem cells (MSCs), chondrocytes, melanocytes, cortical neurons, astrocytes, endothelial cells, adipocytes, dendritic cells (DCs), and unclassified (NA) cell populations (Figures 4 and 5).Figure 4Expression of 11 colocalized genes in different cell types from AD single-cell sequencing data.Multiple violin plots showing expression level by identity for genes MPHOSPH9 through PTGER4.A) Violin plot, MPHOSPH9. X-axis label, Identity. Categories: Monocytes, CD8 plus T cells, Epithelial cells, Keratinocytes, iPS, MSC, Melanocytes, CD4 plus T cells, Fibroblasts, Endothelial cells, Adipocytes, DC, NA.
Y-axis label, Expression Level, range 0 to 3. Highest point cloud density in Monocytes, CD8 plus T cells, Epithelial cells, Keratinocytes. B) AHI1 differs by y-axis range 0 to 5.
C) FADS1 differs by y-axis range 0 to 3 and strongest density in Monocytes. D) DCTN5 differs by y-axis range 0 to 3 with dense early identities and scattered later. E) IRF5 differs by y-axis range 0 to 3 with densest Monocytes.
F) KDELR2 differs by y-axis range 0 to 3 with broader spread across identities. G) CTSS differs by y-axis range 0 to 3 with strongest Monocytes. H) ARL6IP4 differs by y-axis range 0 to 3 with distributions across most identities including NA.
I) JAZF1 differs by y-axis range 0 to 3 with dense Monocytes. J) PPP5C differs by y-axis range 0 to 3 with dense early.Figure 5UMAP plots showing the expression of 11 colocalized genes in AD single-cell transcriptomic data. Each panel represents one gene, with color intensity indicating expression level (gray, low; red, high).Multiple scatter plots showing UMAP gene expression patterns for 11 genes including MPHOSPH9 and AHI1.A) Scatter plot titled, MPHOSPH9.
X axis label umap underscore 1, range negative 15 to 10 with ticks negative 15, negative 10, negative 5, 0, 5, 10. Y axis label umap underscore 2, range negative 15 to 15 with ticks negative 15, negative 10, negative 5, 0, 5, 10, 15. Vertical scale labeled 0, 1, 2, 3.
B) Scatter plot titled, AHI1. Same axes as A). Vertical scale labeled 0, 1, 2, 3, 4.
C) Scatter plot titled, FADS1. Same axes as A). Vertical scale labeled 0, 1, 2, 3.
D) Scatter plot titled, DCTN5. Same axes as A). Vertical scale labeled 0, 1, 2, 3.
E) Scatter plot titled, IRF5. Same axes as A). Vertical scale labeled 0, 1, 2, 3.
F) Scatter plot titled, KDELR2. Same axes as A). Vertical scale labeled 0, 1, 2, 3.
G) Scatter plot titled, CTSS. Same axes as A). Vertical scale labeled 0, 1, 2, 3, 4.
H) Scatter plot titled, ARL6IP4. Same axes as A). Vertical scale labeled 0, 1, 2, 3.
I) Scatter plot titled, JAZF1. Same axes as A).
The results revealed that MPHOSPH9, AHI1, DCTN5, KDELR2, ARL6IP4, JAZF1, PPP5C, and PTGER4 exhibited relatively high expression levels in CD4⁺ T cells, suggesting that these genes may play potential roles in the function or differentiation of CD4⁺ T cells. Notably, MPHOSPH9 and ARL6IP4 showed particularly prominent expression in CD4⁺ T cells compared to other cell types.
Although expression levels varied among cell types, all genes except for IRF5 and FADS1 were clearly expressed in CD4⁺ T cells. Specifically, IRF5 expression was mainly restricted to monocytes, with negligible expression in CD4⁺ T cells, indicating strong immune cell specificity. FADS1 displayed generally low expression and was not prominently expressed in CD4⁺ T cells, implying limited functional involvement in this cell type.
In addition, JAZF1 and PPP5C were broadly expressed across various cell types, but their expression levels in CD4⁺ T cells remained relatively high, suggesting they may be involved in common regulatory processes across multiple cell types, including T cells.
Cellular Communication Network and Pseudotemporal Ordering Analysis of Identified Genes in CD4+T Cells
To explore CD4⁺ T-cell function and fate in AD, we performed comprehensive analyses of cell–cell communication, pseudotime differentiation trajectories, and dynamic gene expression using single-cell transcriptomes. In the ligand–receptor network, CD4⁺ T cells exhibited active intercellular communication in both interaction number and strength, particularly with monocytes, fibroblasts, and dendritic cells, suggesting a bridging role in the immune microenvironment (Figure 6A). Pseudotime analysis (Figure 6B) placed CD4⁺ T cells at intermediate states between naïve and terminally differentiated endpoints.
Cells were arranged along a continuous developmental path with evident branching structures (Figure 6B). The left panel illustrates the progression of pseudotime from early to late states, while the right panels show the classification of cells into distinct transcriptional states and annotated cell types. Notably, hematopoietic and monocytic populations were positioned near the root of the trajectory, suggesting that they represent potential progenitor-like states, whereas macrophages and endothelial cells occupied later positions in pseudotime, consistent with differentiated functional phenotypes.
We further examined the dynamic expression of representative genes along the trajectory (Figure 6C). Genes such as ARL6IP4 and CTSS displayed increased expression along specific branches, indicating their involvement in lineage-specific differentiation programs. These results highlight distinct transcriptional programs driving the bifurcating cell fates within the trajectory.Figure 6Cell–cell communication networks and pseudotime analysis from single-cell sequencing.
(A) The left panel shows the number of interactions between each cell type and the right panel shows the strength of the interaction. The color of the node represents different cell types, and the thickness of the line represents the number or strength of interactions; (B) Pseudotime trajectory analysis. The left panel shows the pseudo-temporal distribution of cells, and the color gradient represents the pseudotime progression.
The right figure is labeled according to cell type, showing the distribution of each cell type on the proposed time series trajectory; (C) Single-cell expression patterns of candidate causal genes. The distribution of single-cell expression of candidate genes along the quasi-temporal trajectory is shown. The color represents the -log10(value + 0.1) transformation value of the gene expression level, and the color bar is graded from purple (low expression) to yellow (high expression).Different data visualizations such as two network graphs and multiple pseudotime scatter plots with gene maps.A) Two circular network graphs.
Text above the first reads, Number of interactions. Text above the second reads, Interaction weight strength. Both graphs place labeled cell types around a circle: Endothelial cells, DC, CD8 plus T cells, CD4 plus T cells, Adipocytes, Unclassified, Monocytes, Melanocytes, Keratinocytes, HSC, Fibroblasts, Epithelial cells.
Curved connecting lines link many pairs of cell types. Line thickness varies to encode interaction amount in the first graph and interaction strength in the second graph. B) Three scatter plots of a branching trajectory.
All plots use axes labeled, Component 1 and Component 2, with tick labels not legible. One plot includes a horizontal color scale labeled, Pseudotime, with tick labels not legible. A second plot labels points by numbered states, with the state numbers not legible.
Identification of Drug Targets and Construction and Analysis of Protein-Protein Interaction (PPI) Network
Using STRING to construct a PPI network, we focused on targets of AD therapeutic drugs—dupilumab, upadacitinib, abrocitinib, and tacrolimus—and assessed their interactions with candidate genes. IRF5 interacts with IL13, IL4, JAK1, and TYK2; CTSS and PTGER4 both interact with IL4; and AHI1 interacts with JAK2. This indicates that IRF5, CTSS, and PTGER4 are connected to the IL-4/IL-13 and JAK-STAT signaling axis, which is a key pathway targeted by approved biologics such as dupilumab and JAK inhibitors.
In addition, the interaction between AHI1 and JAK2 suggests a potential correlation with JAK targeted therapy. No known interactions were identified for the remaining candidates (Figure 7). Several of our prioritized genes may exert their effects through pathways already therapeutically exploited in AD, providing a computational rationale for further exploring these candidates as potential drug targets.Figure 7Protein–protein interaction (PPI) analysis of drug targets and colocalized genes.Protein interaction network linking IL4 and IL13 with JAK kinases, TYK2, IRF5, CTSS and PTGER4.Network diagram of labeled protein nodes connected by many colored interaction lines.
A dense central cluster links IL4 and IL13 with JAK1, JAK2, JAK3, TYK2 and IL4R, with multiple lines between most pairs in this group. IRF5 sits to the right and is connected by lines to IL13, IL4, JAK1 and TYK2. CTSS and PTGER4 are positioned below and each has a single line connecting to IL4.
AHI1 is positioned to the left and has a single line connecting to JAK2. Several other labeled nodes are shown with no connecting lines: JAZF1, MPHOSPH9, PPP3CA, DCTN5, ARL6IP4, FADS1, PPP5C and KDELR2. FKBP1A is connected by multiple lines to PPP3CA.
Simulated Gene Knockout in AD
To investigate the molecular regulatory mechanisms of core genes in AD, a single-cell virtual knockout algorithm was employed to analyze target gene knockouts in AD single-cell transcriptomic data, revealing significant expression perturbations in the top 20 differentially regulated genes downstream of AHI1, ARL6IP4, CTSS, DCTN5, JAZF1, KDELR2, PPP5C, and PTGER4: specifically, the knockout of AHI1 (Figure 8A) and DCTN5 (Figure 8D) predominantly affected skin barrier-related genes DSP (FC = 35.80 and 88.31, respectively) and DSC3 (FC = 34.04 and 85.24, respectively); the knockouts of ARL6IP4 (Figure 8B), CTSS (Figure 8C), JAZF1 (Figure 8E), and PTGER4 (Figure 8H) collectively triggered substantial fluctuations in immune regulatory genes CYBA (max FC = 57.70), TYROBP (max FC = 56.16), and HLA family members such as HLA-DRA, with CTSS and JAZF1 knockouts also significantly impacting GYPC (FC = 99.40 and 84.78, respectively); KDELR2 (Figure 8F) knockout primarily led to expression changes in HLA-DPA1 (FC = 28.46) and VIM (FC = 28.31); while the impact of PPP5C (Figure 8G) knockout was most pronounced in translation-related genes, including RPL31 (FC = 437.78) and RPS17 (FC = 272.64), indicating that these pivotal genes exert central regulatory roles in AD pathogenesis by fine-tuning biological processes such as skin barrier integrity, antigen presentation, and protein translation (Supplementary Table 9). For FADS1, IRF5, and MPHOSPH9, virtual knockout could not be performed due to limitations in data availability and expression levels. Taken together, these virtual knockout results suggest that the eight core genes may regulate AD pathogenesis through three distinct functional modules—skin barrier, immune signaling, and protein translation.Figure 8Downstream regulatory effects of core AD genes via single-cell simulation knockout.
(A–H) Fold changes of TOP 20 differentially regulated genes following simulated knockout of eight candidate genes. The results highlight three functional modules: skin barrier maintenance (DSP, DSC3) regulated by AHI1 and DCTN5; immune activation and antigen presentation (CYBA, TYROBP, HLA family) targeted by ARL6IP4, CTSS, JAZF1, and PTGER4; and protein translation (RPL31, RPS17) perturbed by PPP5C.A set of eight horizontal bar charts of top 20 differentially regulated genes by fold change.A) Horizontal bar chart titled AHI1, Top 20 Differentially Regulated Genes. X axis label, fold change, unit not shown, ticks at 0, 10, 20, 30, 40.
Y axis lists genes including DSP, DSC3, CTSS, VIM, TYROBP, CYBA, GYPC, HLA DRA, HLA DPA1, HLA DQA1. Bars descend from about 40 to about 20. B) Same layout; title ARL6IP4; fold change ticks 0 to 40; top genes include CTSS and multiple HLA genes; bars about 40 down to about 25.
C) Same layout; title CTSS; fold change ticks 0 to 100; top gene GYPC near 100; remaining bars mostly below about 60. D) Same layout; title DCTN5; fold change ticks 0 to 100; top gene DSP near 100; next DSC3 near 90; remaining bars mostly below about 60. E) Same layout; title JAZF1; fold change ticks 0 to 100; top gene GYPC near 85; remaining bars mostly below about 60.
F) Same layout; title KDELR2; fold change ticks 0 to 30; top genes include HLA DPA1 and VIM near 30.
Discussion
This study represents the first comprehensive integration of dynamic CD4⁺ T-cell eQTLs and MR to explore the activation-dependent genetic regulation of AD. Unlike traditional eQTL studies that focus on static or resting immune states, our analysis emphasizes the temporal heterogeneity of gene regulation during T-cell activation, revealing regulatory variants that would likely be missed in non-dynamic contexts.
By leveraging single-cell eQTL data across five activation time points (0 h, LA, 16 h, 40 h, and 5 d), we identified 298 significant causal associations, including 59 dynamic eQTLs across 47 activation conditions. This dynamic perspective unveiled stage-specific gene–disease relationships.20,30 FADS1 exhibited a significant association with AD only upon activation, particularly at 16 h, 40 h, and 5 d, suggesting a role in the effector phase rather than in T-cell differentiation. Similarly, IRF5 and IL1RN showed clear activation-dependent causal effects, consistent with their known roles in inflammation and cytokine regulation.31–33
Importantly, the overlap between dynamic and static eQTL signals highlights both the complementarity and unique contributions of each approach.34,35 KDELR2 showed consistent associations across both datasets and high colocalization probabilities (PPH4 > 0.9), suggesting it may serve as a potential core mediator of AD risk in both steady-state and activated immune cells. Conversely, GNGT2 exhibited opposing effects in plasma and NK cells, underscoring the necessity of cell-type- specific and state-specific analyses to decode complex immune regulation.36–38
The single-cell transcriptomic data further supported the functional relevance of identified candidate genes. MPHOSPH9, AHI1, DCTN5, KDELR2, ARL6IP4, JAZF1, PPP5C, and PTGER4 were highly expressed in CD4⁺ T cells, indicating their potential roles in T-cell function or differentiation. JAZF1 suppresses inflammation by inhibiting NF-κB activation via downregulating TAK1, reducing pro-inflammatory cytokines (eg, IL-1β, IL-6, TNF-α).39,40
To further elucidate the functional consequences of these prioritized genes, we performed single-cell virtual knockout analysis, which revealed distinct regulatory modules. DSP and DSC3 genes are primarily involved in regulating skin barrier function and mediation, among other biological processes.41–43 Notably, targeted knockout of AHI1 and DCTN5 led to a significant alteration in the expression profiles of the desmosome-associated genes DSP and DSC3. Meanwhile, knockout of ARL6IP4, CTSS, JAZF1, and PTGER4 significantly altered the expression of CYBA, TYROBP, and HLA family members (such as HLA-DRA) which are clustered at the hub of immune regulation,44–47 with CTSS and JAZF1 specifically modulating GYPC.
KDELR2 knockout primarily affected HLA-DPA1 and VIM, while the pronounced effect of PPP5C on ribosomal genes RPL31 and RPS17 suggests a critical role in protein translation.48 These findings provide computational evidence that these genes act as central nodes regulating barrier function, antigen presentation, and cellular metabolism in AD, but these predictions require confirmation through orthogonal experimental approaches such as CRISPR-mediated knockout or pharmacological inhibition in relevant cell models.
Genes such as PPP5C, DCTN5, and CTSS demonstrated dynamic expression patterns along pseudotime trajectories, indicating distinct roles at different stages of CD4⁺ T-cell activation—from initiation to effector differentiation. Overexpression of CTSS in transgenic mice led to spontaneous development of chronic atopic dermatitis–like skin lesions. Elevated CTSS drove PAR-2 expression in dendritic cells, enhanced CD4⁺ T-cell differentiation, and skewed toward Th1 cytokines, implicating CTSS in itch induction and chronic AD pathogenesis.49
Previous studies have shown that CD4⁺ T cells play a pivotal role in inflammatory responses, both by promoting immune activation and by maintaining immune homeostasis through regulatory mechanisms.50 Their various subsets—such as Th1, Th2, Th17, and Treg cells—have important roles in a range of inflammatory diseases. Th1 cells primarily mediate cellular immunity by producing interferon-γ (IFN-γ), contributing to the development of autoimmune diseases; Th2 cells, by secreting cytokines such as IL-4, IL-5, and IL-13, play a dominant role in allergic diseases. Th17 cells, characterized by the secretion of IL-17, are involved in the pathogenesis of chronic inflammation and autoimmune diseases; Treg cells, on the other hand, maintain immune tolerance by suppressing excessive immune responses.51,52 Dysregulation of the balance among these subsets is considered a key pathological mechanism in many inflammatory diseases.53 From the perspective of cell–cell communication, CD4⁺ T cells engage extensively with monocytes, dendritic cells, and fibroblasts within the immune microenvironment, indicating that they not only initiate inflammation but also coordinate its magnitude and duration via intercellular signaling networks.
These findings underscore that in AD, CD4⁺ T cells exert both cell-intrinsic functions and cell-extrinsic effects on the broader inflammatory microenvironment.
We then constructed a PPI network using the STRING database and found that several candidate causal genes show predicted functional interactions with the known targets of existing AD therapies (eg, dupilumab, upadacitinib, abrocitinib). IRF5, CTSS, and PTGER4 are predicted to be involved in the IL-4/IL-13 and JAK-STAT pathways, which are the major target of current AD therapies.54 The interaction of AHI1 with JAK2 further suggests a potential influence on response to JAK inhibitors. These interactions provide a computational basis for prioritizing these genes as candidate therapeutic targets for future investigation, while proteins lacking known interactions may serve as starting points for novel drug development efforts.
Functional validation through experimental approaches is required to confirm these predictions.
Several causal genes identified in this study (eg, IL1RN, IRF5) have been confirmed to play critical roles in other autoimmune diseases, such as rheumatoid arthritis, psoriasis,55–59 suggesting that they may represent a shared immune core pathway. FADS1, a lipid metabolism gene, is expressed in sebaceous glands of AD skin. It may influence local lipid composition and inflammation, suggesting a potential role in skin barrier dysfunction and immune regulation in AD.60 The dynamic eQTL approach demonstrated high sensitivity and specificity in this study, indicating its utility as an important tool for target identification in future precision medicine.
By selecting highly expressed, high-risk genes in specific activation states, it may be possible to develop more personalized and temporally targeted therapeutic strategies.
Although this study yielded several important findings, several limitations should be acknowledged. First, the CD4⁺ T-cell eQTL data used in this study were derived from a European population and have not yet been validated in other ethnic groups. Second, our findings are inherently dependent on the quality and characteristics of the public datasets used, although we strengthened causal inference through MR and single-cell virtual knockout, colocalization, and computational analyses (such as pseudotime ordering and PPI mapping), we acknowledge that experimental validation is lacking.
Specifically, no CRISPR KO/KI, small-molecule perturbations, or in vivo AD models were employed to confirm the functional relevance of the prioritized genes. These experimental approaches will be critical in future work to validate whether modulation of the identified targets truly influences CD4⁺ T-cell function and AD pathogenesis. Third, since our MR analysis only includes a single SNP, we are unable to directly test for the pleiotropy of the multivariate pattern using methods such as MR-Egger intercept, which results in our MR analysis potentially being affected by residual pleiotropy.
Fourth, the virtual knockout approach itself has inherent limitations: as a computational simulation, it models network perturbation by altering regulatory connectivity rather than directly knocking down gene expression, and it does not fully recapitulate the biological effects of true gene ablation. Fifth, although dynamic eQTLs provide temporal regulatory information, they cannot fully capture nonlinear and feedback regulatory effects present in complex immune responses; integrating time-series single-cell multi-omics will be an important direction for future exploration. Sixth, we explicitly state that future studies with larger and more diverse dynamic eQTL cohorts are warranted to confirm and extend our findings.
Conclusion
In conclusion, this study identified 298 activation-dependent causal genes for AD, with 11 colocalized core genes confirmed through integrative analyses. Through approaches such as single-cell transcriptomic, pseudotime, cell–cell communication, and PPI analyses, we characterized the dynamic regulatory programs, differentiation trajectories, and interaction networks of CD4⁺ T-cell–related genes involved in AD pathogenesis. Notably, single-cell virtual knockout analysis provided crucial functional validation, revealing that core genes such as AHI1, CTSS, and PPP5C orchestrate AD pathogenesis by precisely regulating skin barrier integrity (eg, DSP, DSC3), immune signaling (eg, CYBA, HLA-DRA), and protein translation (eg, RPL31, RPS17).
These findings provide a comprehensive framework linking genetic regulation to immune activation and offer promising targets for precision and temporally informed immunotherapy in AD.