Exploration of Targets Potentially Linked to IL-17A Inhibitor Response in Psoriasis Using Machine Learning

Introduction

Psoriasis (Ps) is a chronic inflammatory disorder with multi-system involvement, characterized by cutaneous manifestations and significant associations with metabolic syndrome, psoriatic arthritis, and cardiovascular comorbidities^1,2 The pathogenesis of Ps is complex, intricately involving immune, genetic, and environmental factors.^3–5 Early research in the 1980s identified CD4⁺ T-helper cell subsets as key drivers, and subsequent studies have demonstrated that Th1 cells (secreting IFN-γ) and Th17 cells (secreting IL-17/IL-22) play dominant roles in coordinating keratinocyte hyperproliferation, dendritic cell activation, and neutrophil recruitment.^6–11 Based on this understanding of psoriasis pathogenesis, biologics specifically targeting core inflammatory cytokines, such as IL-17 and TNF-α, have been widely adopted in clinical practice, offering significant benefits to psoriasis patients.^12,13 While highly effective, the use of biologics in psoriasis treatment is not without risks. For instance, anti-IL-17 monoclonal antibodies may lead to disease progression during treatment by enabling psoriasis to bypass the blockade pathway. Alternatively, they may shift the inflammatory profile from a predominantly Th1-driven response to a Th2-driven response, resulting in eczema-like manifestations.^14,15 A major clinical challenge is the development of treatment resistance. Some patients show no clinical or molecular response to biologics upon initial treatment, termed primary non-response. Others may experience secondary loss of response, where initially effective treatment gradually becomes ineffective due to acquired resistance.^16,17 The transcriptomic contrasts before and after IL-17A inhibitor treatment may reveal: In primary non-response, baseline expression defects or inhibitory epigenetic marks in key genes of the IL-17A signaling pathway; In secondary loss of response, persistent activation of IL-17A inhibitor-induced negative feedback pathways or aberrant upregulation of genes associated with effector cell exhaustion.

Therefore, conducting in-depth research on the mechanism of action of IL-17 inhibitors in psoriasis holds significant research value and clinical importance. This study employed multiple machine learning diagnostic model analyses on transcriptome sequencing data obtained from skin lesions of psoriasis patients before and after treatment with IL-17A inhibitors. The primary objective was to investigate gene expression changes, immune cell infiltration characteristics, and pathway enrichment following the blockade of the IL-17-mediated pathway. This study aims to provide potential solutions and a theoretical foundation for addressing critical clinical challenges such as treatment non-responsiveness, drug resistance, and recurrence encountered in the current clinical use of IL-17 inhibitors.

Materials and Methods

Data Acquisition and Preprocessing

The specific workflow of this study is illustrated in Figure 1. This study utilized microarray data from the Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo), a public repository of gene expression datasets. This study was exempt from ethical review, per Article 32(1) and (2) of China’s Ethical Review Measures for Life Science and Medical Research Involving Human Subjects (promulgated on February 18, 2023). Specifically, the GSE226244 (n=33 psoriasis pre-treatment, n=28 psoriasis post-treatment, n=8 normal controls), GSE31652 (n=24 psoriasis pre-treatment, n=6 psoriasis post-treatment), and GSE201827 (n=95 psoriasis pre-treatment, n=113 psoriasis post-treatment) datasets were obtained from GEO, encompassing total RNA extracted from epithelial cells of psoriasis patients before and after IL-17A inhibitor treatment. The time points in the training set primarily correspond to 12 weeks post-treatment, while those in the validation set mainly focus on 52 weeks post-treatment. Datasets GSE226244 and GSE31652 served as the training set (n=91), and GSE201827 served as the testing set (n=208). Using a t-test to compare model performance, with an effect size (d) set at 0.5, a post hoc power analysis was conducted in GPower to confirm the adequacy of the sample size. Detailed information on these datasets is provided in Table 1. To ensure data consistency for analysis, the datasets were preprocessed using the SVA package (version 4.2.1) in R software. This preprocessing included gene ID conversion, clinical information extraction, batch effect correction and integration, and homogenizationon the psoriasis datasets GSE226244 and GSE31652. Homogenization transformed data from different batches into a common basis, followed by principal component analysis (PCA) to confirm the effectiveness of batch effect removal.

Table 1 Basic Information on the Datasets Used in This Study

Figure 1 Flowchart Illustrating the Research Methodology.

Identification of Differentially Expressed Genes (DEGs)

This study utilized the Limma package in R to identify DEGs between 57 control cases and 34 treatment cases with IL-17A inhibitors. Limma, a type of generalized linear model, enables the identification of DEGs between different experimental groups. To control for multiple comparisons, the false discovery rate (FDR) was employed to adjust p-values. DEGs were defined by a |log2 fold-change (FC)| greater than 0.585 and a stringent FDR threshold of less than 0.05.

Weighted Gene Co-Expression Network Analysis (WGCNA)

WGCNA, a systematic biological approach, constructs gene co-expression networks to identify modules of genes associated with specific traits. This method investigates co-expression patterns and functional relationships among genes by analyzing gene expression data. It identifies genes with similar expression patterns and clusters them into distinct modules based on their co-expression profiles.¹⁸ Genes within these modules often share similar biological functions. The process involves several steps: 1) data quality assessment and preprocessing, including handling missing values, clustering samples, data matching, and transforming the expression matrix into a correlation matrix; 2) selection of an appropriate soft threshold (β) using the pickSoftThreshold function to ensure scale-free network properties in the adjacency matrix (In this study, β=8); 3) conversion of the adjacency matrix into a topological overlap matrix (TOM) based on the Topological Overlap Measure, which reflects the degree of association between genes; 4) hierarchical clustering of genes based on the TOM-based dissimilarity measure (1-TOM) to group genes with similar expression patterns into modules; and 5) identification of gene co-expression modules using dynamic tree-cutting algorithms. Finally, the method calculates the correlation between module eigengenes, modules, and phenotypes, selects the module with the highest correlation, and computes module identity and gene significance to extract key genes within the selected module.

Functional Enrichment Analysis

This study employed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses to elucidate the biological functions, associated signaling pathways, and gene-gene interactions of co-expressed genes. GO encompasses three domains: Biological Process (BP), Molecular Function (MF), and Cellular Component (CC), providing a comprehensive understanding of gene function. KEGG pathway enrichment analysis was performed to identify the cellular pathways potentially dysregulated in the context of differential gene expression.

Machine Learning

Machine learning encompasses a suite of algorithms and research methodologies that enable machines to learn from data. This involves analyzing vast datasets to identify patterns and subsequently utilize these patterns to make predictions.¹⁹ The study employs a supervised learning approach for analysis. In this study, 113 distinct machine learning methods were implemented to screen core model genes in both the control and treatment groups of patients. The specific algorithms used in this study are detailed in Supplementary Material 1. In machine learning, feature engineering is a critical step to enhance model performance, with its core processes primarily including data preprocessing, feature construction, feature selection, feature transformation, and feature dimensionality reduction using PCA. To ensure data integrity, we implemented stringent measures to prevent data leakage. All preprocessing steps (such as scaling and imputation) were performed independently on the training and test sets. Where applicable, temporal or spatial partitioning was applied, with explicit validation of independence. For model evaluation, we utilized a nested CV framework (5-fold outer CV, 3-fold inner CV) for model training and evaluation. To mitigate class imbalance, we employ SMOTE (Synthetic Minority Oversampling Technique) for synthetic sample generation and Tomek Links to eliminate borderline cases. Hyperparameter optimization was conducted via grid search. All code and scripts used in this study are publicly available in Supplementary Material 2. Subsequently, machine learning models were constructed based on these analyses. Through rigorous evaluation and validation of the performance of each model, the optimal and robust machine learning model for this study was selected.

Analysis of Signaling Pathways Related to Hub Genes

To investigate the interaction networks and signaling pathways involving the six identified genes, we employed the GENEMANIA database (http://genemania.org/). Next, we utilized gene set enrichment analysis (GSEA) to explore gene enrichment patterns across different risk groups. Following this, we performed an analysis of GO term and KEGG pathway enrichment. Gene Set Variation Analysis (GSVA), in contrast to GSEA, is a non-parametric and unsupervised method. Unlike GSEA, which requires samples to be pre-grouped, GSVA calculates enrichment scores for specific gene sets within each sample. This enables the comparison of gene set expression patterns across different treatment groups, revealing enriched signaling pathways.

Immune Infiltration Analysis

The CIBERSORT algorithm, a deconvolution method, infers the relative abundance of 22 distinct human immune cell types within a mixed sample based on gene expression profiles. This is achieved by deconvolving the observed gene expression data with a reference matrix containing gene expression signatures of known immune cell types.^20,21 The algorithm estimates the cellular composition by minimizing the difference between the observed and predicted gene expression profiles, effectively determining the proportion of each immune cell type contributing to the overall gene expression signal. The results are typically visualized using bar charts to depict the proportion of each immune cell type in each sample. Furthermore, the “corrplot” R package can be utilized to generate a heatmap illustrating the correlations among the 22 identified immune cell types. Finally, violin plots were generated to visualize and compare the distributions of immune cell proportions between the control and treatment groups.

Mendelian Randomization Study

In this study, the exposure data were sourced from the GWAS website (https://gwas.mrcieu.ac.uk/datasets/), and the outcome data were obtained from the FinnGen database (https://www.finngen.fi/en/access_results). For the selection of instrumental variables, we utilized cis-eQTLs (tissue-specific P < 5×10⁻⁸) and applied strict linkage disequilibrium removal (r²=0.001, kb = 10,000). The harmonize data function was employed to align effect alleles, and palindromic variants were excluded. All instrumental variables had F-statistics > 10 to mitigate weak instrument bias. Additionally, Bonferroni correction was applied for multiple testing, and MR-PRESSO was employed to identify and correct for outliers. A Mendelian randomization analysis (http://app.mrbase.org/) investigating the causal relationship between GFRA3, CLCNKB, and cytokine levels was performed, utilizing the widely accepted Inverse Variance Weighting (IVW) method for its reliability in calculating the weighted average effect size across all instrumental variables.²² Following this analysis, horizontal pleiotropy and heterogeneity tests were implemented to ensure the validity of the results.^23,24

Cell Culture and Establishment of the Psoriatic Cell Model

The HaCaT cells purchased from ATCC were cultured in DMEM containing 10% fetal bovine serum, penicillin (100 U/mL), and streptomycin (100 μg/mL) at 37 °C in a humidified atmosphere with 5% carbon dioxide. Subsequently, the HaCaT cells were stimulated with M5 (TNF-α, IL-17A, IL-22, IL-1α, and oncostatin M) at a concentration of 10 ng/mL for 24 hours, as previously described.²⁵ For IL-17 inhibition, the HaCaT cells were first treated with M5 for 24 hours, followed by treatment with Secukinumab (MCE, USA, cat. No. HY-P9927, 10 ng/mL) for 1 hour.

qRT-PCR

Total RNA was extracted from HaCaT cells using an RNA Isolation Kit (Aidlab Biotechnologies, China) and subsequently reverse transcribed into cDNA using a reverse transcription kit (Aidlab Biotechnologies, China). q-PCR was performed on the BioRad CFX96 system. Each reaction was carried out in triplicate. β-actin was used as an internal reference for analysis. The primer sequences used in this study are as follows:

CLCNKB: F 5’- GTGGGCATAGTGCGAAGGG-3’,

CLCNKB: R 5’- CAAAGAGGTTGTGTGCCTCAT-3’,

GFRA3: F 5’- TGCACCTCTAGCATAAGCACC-3’,

GFRA3: R 5’- TGCAGCCTATCAGAGAGCTGT-3’,

β-actin: F 5’-GAAGATCAAGATCATTGCTCCT-3’,

β-actin: R 5’- TACTCCTGCTTGCTGATCCA-3’.

Statistical Analysis

All data were processed and analyzed using R software (https://www.r-project.org/) and Perl (Practical Extraction and Report Language). Differences between variables were assessed using the Student’s t-test and the chi-square test. For data that did not follow a normal distribution, non-parametric tests were employed. For gene-level analysis, independent t-tests were performed to compare expression levels between groups, after confirming normality via Shapiro–Wilk tests. For immune-cell proportion analysis, chi-square (χ²) tests were used to assess differences in proportions, as the data were categorical. For qPCR data, non-parametric tests (Mann–Whitney U) were applied for comparisons where normality assumptions were violated. A significance level of p < 0.05 was used throughout the analysis.

Results

Data Preprocessing and Identification of DEGs

In this study, datasets GSE226244 and GSE31652 served as the training set, and GSE201827 served as the testing set. Using the SVA package (version 4.2.1) in R, we performed gene ID conversion, clinical information extraction, batch effect correction and integration, and homogenization of the psoriasis datasets GSE226244 and GSE31652. The resulting training dataset comprised 57 control samples and 34 samples from patients treated with an IL-17A inhibitor. Evaluation of the data confirmed the effectiveness and reliability of the preprocessing. Box plots demonstrated significant differences in sample distributions across datasets before batch effect removal, indicating the presence of batch effects. Following batch effect correction, the sample distributions became consistent across datasets (Figure 2A and B). PCA confirmed consistency in the principal components after processing, demonstrating the suitability of the data for further analysis (Figure 2C and D). Using the Limma package in R, DEGs were identified. A total of 2200 DEGs were identified between the two patient groups, comprising 1126 upregulated and 1074 downregulated genes. Heatmaps and volcano plots were generated to visualize these DEGs (Figure 2E and F).

Figure 2 Data preprocessing and Identification of DEGs (A) Boxplot before batch correction. (B) Boxplot after batch correction. (C) Principal component analysis of the two original Psoriasis datasets before batch effect correction. (D) Principal component analysis of the corrected Psoriasis dataset. (E) Heatmap of DEGs between the control group and treatment group. (F) Volcano plot of DEGs between the control group and treatment group.

Construction of Weighted Gene Co-Expression Networks and Identification of Modules

WGCNA was performed on the training dataset. Figure 3A demonstrates that, with a power value of 8, the scale-free fit index reached 0.883 with relatively low average connectivity, consistent with the characteristics of a scale-free network. Following WGCNA, nine distinct characteristic modules were identified (Figure 3B). Among these modules, the module eigengenes in the pink module exhibited the highest correlation (r=0.91, p=4e-35) between the control and treatment groups, with statistical significance (Figure 3C). The importance of feature genes in the nine modules was evaluated. The feature genes in the pink module exhibited the highest importance scores and demonstrated a strong correlation with the module traits (Figure 3D and E).

Figure 3 Construction of Weighted Gene Co-expression Networks and Identification of Modules (A) Scale Independence and Mean connectivity plots of psoriasis. (B) Gene dendrogram distribution in different color modules. (C) Traits of different color modules in psoriasis control group and treatment group. (D) Importance scoring diagram of feature genes in different color modules. (E) The correlation and importance distribution map of characteristic genes in the pink module. Correlation plot between module genes (MM, X-axis) in the pink module and genes in traits (GS, Y-axis).

Functional Enrichment Analysis of Co-Expressed Genes

Figure 2F shows 2200 DEGs between the psoriasis control and treatment groups. WGCNA identified 217 hub genes within the pink module. The intersection of these two gene lists revealed 95 commonly expressed genes (Figure 4A). GO and KEGG pathway analyses were performed on these 95 genes, and the results are presented in Figure 4. GO analysis indicated enrichment of these genes in biological processes, including epithelial tube morphogenesis, embryonic organ development, gonad development, development of primary sexual characteristics, negative regulation of neurogenesis, negative regulation of axonogenesis, and the activin receptor signaling pathway. Cellular component analysis showed enrichment in the glutamatergic synapse, synaptic membrane, postsynaptic specialization membrane, an integral component of the postsynaptic membrane, an intrinsic component of the postsynaptic membrane, bicellular tight junction, tight junction, presynaptic membrane, and apical junction complex. Molecular function analysis revealed enrichment in DNA-binding transcription repressor or activator activity, RNA polymerase II-specific basal transcription machinery binding, basal RNA polymerase II transcription machinery binding, UDP-glycosyltransferase activity, and hexosyltransferase activity (Figure 4B-D). KEGG pathway analysis revealed significant enrichment of the co-expressed genes in several signaling pathways, including circadian rhythm, FoxO signaling, Apelin signaling, long-term depression, Cushing syndrome, lysine degradation, the cGMP-PKG signaling pathway, and renin secretion (Figure 4E and F).

Figure 4 Functional enrichment analysis of co-expressed genes (A) Venn diagram of DEGs and module eigengenes of the pink module from WGCNA analysis. (B) The barplot chart in GO analysis, where BP stands for Biological Process, CC stands for Cellular Component, and MF stands for Molecular Function. (C) Bubble Chart in GO Analysis. (D) Circular Graph in GO Analysis. (E) The barplot chart for KEGG analysis. (F) The bubble chart for KEGG analysis.

Construction and Validation of Machine Learning Models

This study employed 113 distinct machine-learning methods to screen core genes from psoriasis patients in both control and treatment groups and to construct predictive models. Through evaluation and validation of model performance, the optimal model was identified. Figure 5A demonstrates that the model constructed using the Least Absolute Shrinkage and Selection Operator (LASSO)-combined glmBoost algorithm exhibited the highest diagnostic performance. Receiver operating characteristic (ROC) curve analysis revealed that this model achieved an area under the curve (AUC) of 0.920 (95% CI: 0.862, 0.968) in the training set and 0.858 (95% CI: 0.798, 0.911) in the test set, indicating a strong diagnostic value (Figure 5B and C). Performance metrics for the LASSO-combined glmBoost model were as follows: in the training set, accuracy was 92.22%, precision was 88.89%, recall was 98.25%, and the F1 score was 0.9514; in the validation set, accuracy was 83.65%, precision was 76.99%, recall was 91.58%, and the F1 score was 0.8366 (Figure 5D and E). This model identified six differentially expressed genes. Compared to the control group, CLCNKB, GFRA3, ZNF471, CNKSR2, and SSTR1 were significantly upregulated in the treatment group, while GJA3 was downregulated (Figure 5F and G). Gene correlation analysis revealed significant correlations among these hub genes (Figure 5H). ROC curve analysis for each of these six differentially expressed genes demonstrated high AUC values for CLCNKB, GFRA3, SSTR1, ZNF471, and CNKSR2, indicating their high accuracy in distinguishing between patients in the psoriasis control and treatment groups (Figure 5I). This suggests a potential involvement of these genes in IL-17-related signaling pathways.

Figure 5 Construction and Validation of Machine Learning Models (A) Exhibition of multiple machine learning diagnostic models. (B) ROC curve of the training group under the optimal machine learning model. (C) ROC curve of the test group under the optimal machine learning model. (D) Confusion matrix of the training group. (E) Confusion matrix of the test group. (F) Volcano plot of the differential model genes. (G) Expression levels of differential model genes in the psoriasis control group and treatment group. (H) Correlation analysis of differential model genes. (I) ROC curve of differential model genes. ***p<0.001.

Functional Validation of CLCNKB and GFRA3 via Enrichment Analysis and Cellular Assays

Figure 5 demonstrates that the expression of CLCNKB and GFRA3 is significantly upregulated in psoriasis patients following treatment with an IL-17 inhibitor, correlating with higher AUC values. Consequently, relevant signaling pathways associated with CLCNKB and GFRA3 were further investigated. Using the GeneMania database, the functions of CLCNKB and GFRA3 were predicted, and a gene interaction network was constructed, revealing interactions with multiple proteins that contribute to cellular functions (Figure 6A). GSVA indicated that CLCNKB was primarily involved in pathways such as chemokine signaling, natural killer cell-mediated cytotoxicity, antigen processing and presentation, leukocyte transendothelial migration, cytokine-cytokine receptor interaction, base excision repair, and homologous recombination. GFRA3 was primarily associated with pathways such as T cell receptor signaling and the tricarboxylic acid (TCA) cycle (Figure 6B and C). GSEA suggested that CLCNKB and GFRA3 could participate in processes such as maintaining taste receptor activity (Figure 6D and E). Cellular experiments were conducted to assess the expression levels of CLCNKB and GFRA3 following treatment with M5 alone and in combination with the IL-17 inhibitor HYP9927. Consistent with previous findings, the expression levels of CLCNKB and GFRA3 were upregulated after treatment with the combined M5 and HYP9927 regimen (Figure 6F and G).

Figure 6 Functional Validation of CLCNKB and GFRA3 via Enrichment Analysis and Cellular Assays (A) Protein-protein interaction network of CLCNKB and GFRA3. (B) GSVA chart illustrating the enrichment of gene sets associated with CLCNKB. (C) GSVA chart illustrating the enrichment of gene sets associated with GFRA3. (D) GSEA chart illustrating the enrichment of gene sets associated with CLCNKB. (E) GSEA chart illustrating the enrichment of gene sets associated with GFRA3. (F) Changes in mRNA levels of CLCNKB in HaCaT cells after treatment with different factors; (G) Changes in mRNA levels of GFRA3 in HaCaT cells after treatment with different factors. *p<0.05, **p<0.01.

The Immune Cell Infiltration Landscape of Psoriasis Treatment with IL17 Inhibitors

CIBERSORT, a commonly used algorithm for analyzing the distribution of immune cell types in lesional tissues, was employed to analyze immune cell infiltration characteristics in psoriasis patients in both control and treatment groups. Compared with the control group, psoriatic skin lesions treated with IL-17A inhibitors exhibited higher levels of resting mast cells, M2 macrophages, CD8 T cells, and plasma cells, and lower levels of activated dendritic cells, resting dendritic cells, follicular helper T cells, activated CD4 memory T cells, and naïve B cells (Figure 7A and B). Correlation analysis of 22 immune cell types revealed several significant correlations: follicular helper T cells were negatively correlated with regulatory T cells (Tregs) (r = −0.65); resting dendritic cells were negatively correlated with M2 macrophages (r = −0.51); activated dendritic cells were negatively correlated with resting mast cells (r=−0.73); CD8 T cells were positively correlated with Tregs (r = 0.53); memory B cells were positively correlated with naïve CD4 T cells (r=0.68), and naïve CD4 T cells were positively correlated with Tregs (r=0.77) (Figure 7C). After adjusting the CIBERSORT p-value to 0.01, the proportions of key cell types showed no significant changes, demonstrating that the core conclusions are not affected by the threshold selection. The relevant graphics can be found in Figure S1.

Figure 7 The immune cell infiltration landscape of psoriasis treated with IL17 inhibitors (A) The distribution ratio of 22 types of immune cells in the control group and treatment group patients. (B) Comparison of the distribution of 22 immune cell subtypes between the control group and the treatment group of psoriasis patients (green represents the control group, and red represents the treatment group). (C) Correlation matrix diagram of 22 immune cells. (D) Correlation between 6 differential genes and 22 immune cell subtypes (red indicates positive correlation, and green indicates negative correlation). * p < 0.05; **p < 0.01; ***p < 0.001.

Further correlation analysis between the CIBERSORT results and the expression data of the hub genes revealed associations between the expression levels of CLCNKB, CNKSR2, GFRA3, and GJA3 and immune cell infiltration patterns (Figure 7D). Specifically, CLCNKB expression was negatively correlated with gamma delta T cell infiltration; CNKSR2 expression was positively correlated with neutrophil infiltration; GFRA3 expression was negatively correlated with gamma delta T cell infiltration; GJA3 expression was positively correlated with M0 macrophage and follicular helper T cell infiltration; and GJA3 expression was negatively correlated with resting mast cell infiltration (Figure 8A–F).

Figure 8 The correlation between CLCNKB, CNKSR2 and immune cell infiltration (A) The correlation between CLCNKB expression level and gamma delta T cells. (B) The correlation between CNKSR2 expression level and Neutrophils. (C) The correlation between GFRA3 expression level and gamma delta T cells. (D) The correlation between GJA3 expression level and M0 Macrophages. (E) The correlation between GJA3 expression level and resting Mast cells. (F) The correlation between GJA3 expression level and follicular helper T cells. The Orange area in the figure is the visualization of the Confidence Interval. The blue area is the visualization of the Regression Line.

Mendelian Randomization Analysis of CLCNKB, GFRA3 and Cytokines

Causal relationships between CLCNKB, GFRA3, and various cytokines were further analyzed using the IVW method implemented in the TwoSampleMR package. As shown in Table 2, CLCNKB exhibited causal associations with IL-13, IL-10, IL-12, and TNFα, while GFRA3 showed causal associations with IL-13, IL-12, and TGFβ. To facilitate visualization of the Mendelian Randomization (MR) analysis results, scatter plots were generated (Figure 9), displaying regression lines for IVW, MR-Egger, and weighted median analyses. Individual causal effect analysis revealed that increasing effects of individual single nucleotide polymorphisms (SNPs) on GFRA3 corresponded to increased promoting effects on IL-12 and TGFβ and an increased inhibitory effect on IL-13 (Figure 9A–C). Similarly, increasing effects of individual SNPs on CLCNKB corresponded to increased inhibitory effects on IL-10, IL-12, IL-13, and TNFα (Figure 9D–G). The Cochran’s Q test indicated no significant heterogeneity among the instrumental variables (P > 0.05), suggesting homogeneity at the genetic variant level. The MR-Egger intercept test, used to assess the presence of horizontal pleiotropy, showed intercept values close to zero with P > 0.05, indicating no significant pleiotropy. MR-PRESSO was employed to identify and correct for outliers, and the global test detected no significant outliers. Finally, after applying Bonferroni correction to the cytokine analyses, the results remained unchanged. These analyses collectively further validate the robustness of our findings.

Table 2 Mendelian Randomization Analysis of GFRA3, CLCNK and Cytokines

Figure 9 Scatter plot of causal effects of GFRA3, CLCNKB and cytokines in Mendelian randomization analysis (A) Scatter plot of causal effects between GFRA3 and IL-12. (B) Scatter plot of causal effects between GFRA3 and TGF-β. (C) Scatter plot of causal effects between GFRA3 and IL-13. (D) Scatter plot of causal effects between CLCNKB and IL-10. (E) Scatter plot of causal effects between CLCNKB and IL-12. (F) Scatter plot of causal effects between CLCNKB and IL-13. (G) Scatter plot of causal effects between CLCNKB and TNFα.

Discussion

Studies have demonstrated that the IL-23/IL-17 inflammatory pathway plays a critical role in the development and progression of psoriasis.²⁶ Consequently, targeting the IL-17A signaling pathway is considered an effective therapeutic strategy for psoriasis. Numerous clinical studies have demonstrated the rapid and sustained efficacy of IL-17A inhibitors in improving psoriatic lesions, with favorable safety profiles.²⁷ Although significant progress has been made in elucidating the mechanism of action of IL-17A inhibitors in the treatment of psoriasis, the precise molecular mechanism remains unclear; more importantly, the mechanism of resistance to IL-17A inhibitors also needs further exploration.

In the present study, 113 machine learning algorithms were applied to select hub genes differentiating control and treatment groups, and a predictive model was constructed. Through rigorous evaluation and validation of multiple models, the optimal model was selected. Comparison of the treatment and control groups revealed upregulation of CLCNKB, GFRA3, ZNF471, CNKSR2, and SSTR1, and downregulation of GJA3 in the treatment group. Correlation analysis demonstrated significant relationships among these genes. Notably, CLCNKB, GFRA3, SSTR1, ZNF471, and CNKSR2 exhibited high AUC values, suggesting their potential as accurate biomarkers for distinguishing between psoriasis control and treatment groups. These findings suggest that these genes and their associated pathways may play crucial roles in the efficacy of IL-17A inhibitors in psoriasis treatment. In vitro experiments further demonstrated that combined treatment with M5 and the IL-17 inhibitor HYP9927 significantly increased the expression of CLCNKB and GFRA3, suggesting that IL-17A upregulate the expression of these genes. Functional enrichment analysis indicated that CLCNKB is primarily involved in pathways such as chemokine signaling, natural killer cell-mediated cytotoxicity, antigen processing and presentation, leukocyte transendothelial migration, cytokine-cytokine receptor interaction, base excision repair, and homologous recombination. GFRA3 is primarily associated with T cell receptor signaling and the tricarboxylic acid (TCA) cycle.

Psoriasis is a T-cell-mediated inflammatory disease, where abnormal infiltration and increased activity of specific T-cell subsets, such as γδ and βδ T cells, play a critical role in its pathogenesis.²⁸ The expression of CLCNKB, CNKSR2, GFRA3, and GJA3 was significantly associated with immune cell infiltration patterns. Specifically, CLCNKB expression exhibited a negative correlation with γδ T cell infiltration, while GFRA3 expression showed a negative correlation with βδ T cell infiltration. These findings suggest that CLCNKB and GFRA3 may play a role in regulating the function and infiltration of these T cell subsets. Further Mendelian randomization analysis revealed causal associations between CLCNKB, GFRA3, and several cytokines. CLCNKB was associated with IL-10, IL-12, IL-13, and TNF-α, while GFRA3 exhibited associations with IL-12, IL-13, and TGF-β. Cytokines such as IL-13, IL-10, IL-12, TNF-α, and TGF-β play significant roles in psoriasis²⁹ IL-10, a key anti-inflammatory cytokine primarily secreted by regulatory T cells (Tregs), B cells, and macrophages, suppresses Th1 and Th17 inflammatory responses and is significantly reduced in patients with psoriasis^29,30 IL-12 promotes Th1-mediated immune responses and is critically involved in the development of psoriasis.³¹ TNF-α, a canonical pro-inflammatory cytokine, activates the NF-κB pathway, promoting keratinocyte hyperproliferation and immune cell recruitment.³² TGF-β, while exhibiting anti-inflammatory properties under normal conditions, can promote Th17 differentiation under pathological conditions, thereby exacerbating psoriatic inflammation.³³ Therefore, IL-17A, as the central effector cytokine in psoriasis, may interact with multiple other cytokines, including IL-13, IL-10, IL-12, TNF-α, and TGF-β, through various immune pathways. These interactions collectively influence the development and severity of psoriasis, forming a complex network of effects. Our findings suggest that IL-17A inhibitors may indirectly modulate psoriasis progression by targeting CLCNKB and GFRA3 as key immunoregulatory factors.

This study reveals a novel mechanism of IL-17A inhibitors in psoriasis treatment, involving the modulation of CLCNKB and GFRA3 pathways beyond direct IL-17A pathway inhibition. Mutations or dysregulation of these genes may influence drug efficacy and contribute to immune escape, suggesting their potential as combination therapy targets. The key differences between this study and previous multi-cohort WGCNA/ML analyses lie in two aspects: at the data level, it specifically focuses on psoriasis patients before and after IL-17A inhibitor treatment in the GEO database, resulting in a significantly distinct sample composition compared to prior studies; at the methodological level, it innovatively introduced MR analysis to achieve multi-dimensional validation, while strictly satisfying the key assumptions of MR analysis to demonstrate robust causal associations. This study has limitations that should be acknowledged, including the relatively small sample size and the absence of clinical efficacy data correlated with observed genomic changes. Further investigation of the functions and regulatory mechanisms of the identified differentially expressed genes in cellular and animal models is warranted. While IL-17A inhibitors have significantly advanced psoriasis treatment, their mechanisms of action remain complex, involving multiple pathways and immune cell subsets. A comprehensive understanding of these changes will facilitate the discovery of novel therapeutic targets and strategies, ultimately guiding the precise application of IL-17A inhibitors to optimize outcomes for psoriasis patients.

Conclusion

This study provides novel insights into the mechanisms of IL-17 inhibitor response in psoriasis, revealing that upregulation of CLCNKB and GFRA3, along with cytokine dysregulation (eg, IL-13, IL-10, IL-12, TGF-β, TNF-α), may contribute to on-treatment biological changes consistent with resistance/relapse hypotheses. While these findings offer a promising approach for predicting clinical outcomes, their clinical utility requires further validation through patient-level response phenotyping and longitudinal studies in real-world settings to further facilitate the clinical translation of novel psoriasis treatment strategies targeting CLCNKB and GFRA3.