Deep analysis of skin molecular heterogeneities and their significance on the precise treatment of patients with psoriasis

Background: Psoriasis is a highly heterogeneous autoinflammatory disease. At present, heterogeneity in disease has not been adequately translated into concrete treatment options. Our aim was to develop and verify a new stratification scheme that identifies the heterogeneity of psoriasis by the integration of large-scale transcriptomic profiles, thereby identifying patient subtypes and providing personalized treatment options whenever possible.

Methods: We performed functional enrichment and network analysis of upregulated differentially expressed genes using microarray datasets of lesional and non-lesional skin samples from 250 psoriatic patients. Unsupervised clustering methods were used to identify the skin subtypes. Finally, an Xgboost classifier was utilized to predict the effects of methotrexate and commonly prescribed biologics on skin subtypes.

Results: Based on the 163 upregulated differentially expressed genes, psoriasis patients were categorized into three subtypes (subtypes A–C). Immune cells and proinflammatory-related pathways were markedly activated in subtype A, named immune activation. Contrastingly, subtype C, named stroma proliferation, was enriched in integrated stroma cells and tissue proliferation-related signaling pathways. Subtype B was modestly activated in all the signaling pathways. Notably, subtypes A and B presented good responses to methotrexate and interleukin-12/23 inhibitors (ustekinumab) but inadequate responses to tumor necrosis factor-α inhibitors and interleukin-17A receptor inhibitors. Contrastly, subtype C exhibited excellent responses to tumor necrosis factor-α inhibitors (etanercept) and interleukin-17A receptor inhibitors (brodalumab) but not methotrexate and interleukin-12/23 inhibitors.

Conclusions: Psoriasis patients can be assorted into three subtypes with different molecular and cellular characteristics based on the heterogeneity of the skin's immune cells and the stroma, determining the clinical responses of conventional therapies.

1 Introduction

Psoriasis is a common chronic auto-inflammatory skin disease characterized by obviously erythematous and scaly skin lesions accompanied by prominent skin and joint manifestations (1, 2). Chronic plaque or psoriasis vulgaris is the most common form of psoriasis, where lesions arise from hyperproliferation and disrupted differentiation of epidermal keratinocytes provoked by innate and adaptive immune system dysfunction (3, 4). The American Academy of Dermatology proposes biologics as a first-line treatment option for moderate to severe plaque psoriasis owing to their superior efficacy and acceptable safety profile compared to other treatment options. Specifically, the most commonly prescribed biologics are inhibitors to tumor necrosis factor-α (TNF-α; including etanercept and adalimumab) and cytokine-targeting treatments such as the p40 subunit of interleukin (IL)-12/23 (ustekinumab) and IL-17 (brodalumab) (5). The systemic anti-inflammatory response caused by most of these biologics in treating psoriasis decreases the severity of comorbidities. Nevertheless, a considerable number of patients with psoriasis still present an inadequate response to these therapies (6–8), which may result from the heterogeneous pathophysiological background of psoriasis.

Gene expression profiling of psoriasis skin tissues has been used to provide pathophysiological insights for explaining the variations in treatment outcomes. Ainali et al. have revealed two distinct molecular subgroups from a single chronic plaque psoriasis skin transcriptomic profile and noted that one of these subgroups enriched in the transforming growth factor-β and ErbB signaling pathways might be more amenable to biologics therapies (9). Wang et al. have delineated two distinct immune phenotypes and constructed a psoriatic microenvironment score from a meta-analysis of psoriasis skin transcriptomic profiles (10). There is a correlation between improved clinical outcomes and overexpression of genes responsible for keratinocyte differentiation in the nonlesional phenotype. Some characteristics are common among the various subtypes, but their conclusions are inconsistent. Consequently, gaining more profound and comprehensive mechanistic insights into divergent and shared features of the psoriasis skin subtypes is necessary for determining the pathobiological approaches for treatment-resistant patients with psoriasis.

In this study, we aimed to elucidate the cellular, molecular, and clinical features of three newly identified psoriasis skin subtypes using an unsupervised machine learning method performing a comprehensive meta-analysis of publicly available microarray datasets published thus far. Finally, we sought to apply a machine-learning model to identify the psoriasis skin subtypes for evaluating the therapeutic efficacy of biologics.

2 Methods

2.1 Systematic search, data selection, and preprocessing

The study design is outlined in Figure 1. We retrieved the microarray gene expression datasets of psoriasis skin datasets from the Gene Expression Omnibus database. Nine microarray datasets (GSE30999 (11), GSE13355 (12), GSE14905 (13), GSE41662 (14), GSE53552 (15), GSE34248 (14), GSE67853 (16), GSE47751 (17), and GSE50790 (18) of eligible psoriasis skin tissues were collected for this study (Supplementary Table 1). The independent cohorts for treatment include 66, 85, 73 and 15 psoriasis lesions that were treated with etanercept (14, 19), ustekinumab (19, 20), brodalumab (20), and methotrexate (21), respectively, to which therapeutic responses were measured using the Psoriasis Area and Severity Index (PASI) score endpoint after the initiation of therapy. Here, patients were considered responders when there was an improvement of at least 75% on the PASI score from baseline and were considered non-responders otherwise.

Figure 1

Figure 1 Overview of data processing steps. From the public databases, nine microarray datasets containing 250 patients with psoriasis were selected. According to the established methodologies, DEGs were filtered, and enrichment analysis, PPI network analysis, and supervised clustering were performed. Finally, the Xgboost classifier was constructed to predict the responses of stratified subtypes to commonly used biologics. DEGs, differentially expressed genes; GO, Gene Ontology; GSEA, Gene-set enrichment analysis; KEGG, Kyoto Encyclopedia of Genes and Genomes; MCODE, Molecular Complex Detection; PPI, protein–protein interaction.

Using the affy R package, microarray raw data from Affymetrix was processed by employing the robust multi-array average algorithm for background correction, quantile normalization, and probe-set summarization. Normalized matrix files are downloaded directly from Illumina into raw microarray data.

Residual technical batch effects from different data sets were corrected using the ‘ComBat’ function in the sva R package to reduce the systematic, dataset-specific bias (22). Quality assurance and distribution bias were assessed using the principal component analysis (PCA) of the same datasets before and after normalization and batch correction (Supplementary Image 1).

2.2 Filtering of differentially expressed genes and functional comments

The limma R package performed gene identification by differential expression analysis of microarray data from lesional and non-lesional psoriasis skin tissues using a linear model and modified t-test (23). To control for the proportion of false positive findings, we adjusted P-values using the Benjamini-Hohberg correction (24), and adjusted P-values <0.05 and absolute fold change values ≥1.5 were considered statistically significant. Then, the functional enrichment analysis of the upregulated DEGs lists in psoriatic lesions skin was performed via the Metascape online website. Adjusted p <0.05 was deemed to be significant for enriched functional pathways (25).

2.3 Construction of protein-protein interaction network and identification of the important modules

To visualize the interconnectedness of DEGs in the psoriasis skin samples, we constructed a protein-protein interaction (PPI) network based on the STRING and BioGrid databases and Cytoscape software (26, 27). The network is described by proteins (i.e., nodes) and their relationships (physical or functional interactions) (i.e., edges). Four significant modules were identified during the analysis using the Molecular Complex Detection (MCODE) algorithm (28). The biological functions of important module genes were identified by enrichment analysis. Based on the P value, the three highest-scoring terms were retained to describe the function of the corresponding module.

2.4 Gene-set enrichment analysis

Gene-set enrichment analysis (GSEA) ranked the list of upregulated DEGs in psoriasis-lesioned skin using a predefined gene set database to identify potential biologically critical pathways related to psoriasis (29, 30). Information on gene sets for signaling pathways or biological processes was obtained from the Kyoto Encyclopedia of Genes and Genomes (KEGG) and Reactome databases. Terms with a false discovery rate below 0.25 were considered significant.

2.5 Unsupervised clustering for psoriasis skin gene expression profiles

To classify patients with psoriasis into subtypes based on the skin transcriptomic signatures, we performed hierarchical agglomerative clustering using the ConsensuClusterPlus R package (31). The procedure was repeated for 1000 iterations to ensure the stability of the subtype stratification. The clustering approach adopted a Partitioning Around Medoids algorithm with the Euclidean distance and Ward-D2 linkage, and the optimal and stable number of clusters was selected by the cumulative distribution function (CDF). We then used the PCA to confirm the results of unsupervised clustering. DEGs were differentially expressed in the stroma- and immune-enriched subtypes compared to those differentially expressed in the other subtypes. An adjusted P-value <0.05 and absolute fold-change ≥1.5 were used to identify enriched pathways in the Gene Ontology Biological Process (GO-BP) and Reactome databases.

2.6 Inference of cell types and pathways activated in psoriasis skin subtypes

To explore the specific biological pathways associated with the subtypes, we used single-sample GSEA to condense information from gene expression profiles into pathway or marker gene sets (32). Here, the enrichment score is defined as the absolute enrichment of the gene set in each sample of a given dataset. Psoriasis skin-related pathways were collated from the published literature and GSEA results, with gene sets from KEGG and response group databases. We used the “xCell” algorithm to calculate the enrichment scores for immune and stroma cell-type signatures for elucidating the cellular composition of the three subtypes (33). Enrichment scores representing given cell types and pathway activities were compared among the three subgroups using the Wilcoxon and Kruskal-Wallis test.

2.7 Construction of multi-classification model for psoriasis skin subtype prediction

The Xgboost [extreme gradient boosting] (version 1.5.0.2) is a machine learning method that assembles weak prediction models to build more robust prediction models. We made a prediction model based on 163 gene features using the Xgboost-tree method with softmax objective function. The predictive power of the system was then measured using the average area under the receiver operating characteristic curve (AUC). For training the classifier, 250 psoriasis skin samples were divided into training (n=176) and testing (n=74) sets in 70% and 30% proportions, respectively, which used the caret R package. We controlled overfitting in the model using 10-fold cross-validation and used the fitted model to select subtypes for the new samples. Finally, we applied the 163-gene classifier for categorizing the samples into subtypes.

2.8 Statistical analysis

All statistical analyses were performed using the R software (version 4.0.3). We used the Kruskal–Wallis test to compare differences among two or more samples and the Wilcoxon test to compare two samples. The chi-square or Fisher’s test was applied to analyze the relationship between the psoriasis skin subtypes and clinical measures. Statistical significance was set at P<0.05 (two-tailed test).

3 Results

3.1 DEGs, PPI network, and enriched signaling pathways

On comparing the gene expression profiles of lesional and non-lesional skin samples from 250 patients with psoriasis, 163 upregulated DEGs were obtained (Figures 2A, B). Gene ontology [GO] annotation demonstrated that these upregulated DEGs in lesional skin were appreciably enriched in biological processes(BP) such as keratinization; inflammatory response; response to viruses, bacteria, and fungi; and type I interferon (IFN) production (Figure 2C). The KEGG and Reactome analysis indicated that the upregulated DEGs primarily enriched the IL-17 signaling, p53 signaling, neutrophil degranulation, and IFN-α/β signaling pathways (Figures 2D, E).

Figure 2

Figure 2 Upregulated DEGs analysis between lesional and nonlesional skin samples of patients with psoriasis. (A, B) Volcano plot and heatmap of all the DEGs. (C) Functional enrichment analysis of upregulated DEGs. The top 20 most significant biological processes in the GO-BP database. (D, E) Pathway enrichment analysis of upregulated DEGs. Top 5 most considerable signaling pathways in the KEGG and Reactome databases. (F) Protein–protein network of upregulated DEGs. The nodes and edges of the network represent genes and the functional or physical relationships between them, respectively. Four modules were found to be significant using the MCODE algorithm. GO-BP, Gene Ontology Biological Process.

In total, 188 interactions among 82 DEGs were identified in the PPI networks constructed using these DEGs, and 81 genes were isolated without any direct relation to each other. The MCODE analysis identified 29 hub genes and clustered them into four highly correlated protein clusters (Figure 2F). Enrichment analysis was performed separately for each module, and a functional descriptor was chosen for each module based on the three terms that scored best based on the P-value. In line with our expectations, combining the result of GSEA, IFN-related (IFN-α/β/γ signaling and response to viruses) genes, neutrophil-related (neutrophil degranulation and response to bacteria and fungi) genes, toll-like receptor (TLR), NOD-like receptor (NLR), RIG-I-like receptor (RLR) signaling pathways, and keratinization were significantly enriched in lesional skin, adequately explaining the molecular mechanism of psoriasis skin.

3.2 Identification of skin gene expression-driven subtypes

Using the “ConsensusClusterPlus” R package of 1,000 iterations, we evaluated the number of clusters from k = 2 to k = 6. The CDF value and the delta area were used to measure the robustness of clustering results. The results show that when K=3, the consensus matrix (Figure 3A), consensus CDF plot (Figure 3B), and delta area plot (Figure 3C) all showed stable results, and the clustering consistency scores of each subgroup exceeded 0.8 (Figure 3D). Psoriasis skin subtype segregation patterns were revealed using the PCA (Figure 3E). The heatmap plot showed upregulated DEG in the three subtypes, showing the gene-level variability of the three subtypes (Figure 3F).

Figure 3

Figure 3 Identification of psoriasis skin subtypes based on upregulated DEGs. (A) A recorded consensus matrix at k=3 for the psoriasis skin compendium. The values of the consistency matrix from white to dark blue are from 0 to 1. (B) Consensus clustering for the CDF for k=2–6. (C) Relative change in area under CDF curve for k=2–6. (D) The cluster consistency score for each subgroup of k=2–6, and the red horizontal line indicates a cluster consistency score of 0.8. (E) PCA for the DEG expression profiles shows the stability and reliability of the classification. (F) A heatmap of 250 patients with psoriasis with a red gradient illustrating the gene expression distribution for three psoriasis skin subtypes. In each column, patients are grouped based on cluster assignment. Red represents overexpression, while blue represents under-expression. CDF, cumulative distribution function; PCA, Principal components analysis.

We compared each cluster with the others (such as subtype A vs subtype B and subtype C) to determine the specific upregulation of DEG signatures in each psoriasis skin subtype and to evaluate the molecular signatures and biological processes of the corresponding subtype. Here, 144 DEGs were significantly upregulated in subtype A, only 1 in subtype B, and 305 in subtype C using the same filtering threshold (Figure 4A). Thereafter, the Metascape was used to enrich the most obviously dysregulated biological processes and signaling pathways of each subtype from the Gene Ontology Biological Process (GO-BP), KEGG, and Reactome databases. To be specific, there was a significant activation of canonical immune pathways in subtype A, such as neutrophil-related pathways (including neutrophil degranulation and response to bacteria or fungi) and IFN-related pathways (including response to IFN-α/β/γ and response to viruses) (Figures 4B–D). There was an enrichment of stroma proliferating pathways in subtype C, including cornified envelope formation, peroxisome proliferator-activated receptor (PPAR) signaling pathway, and Wnt signaling pathway (Figures 4E–G). Additionally, the mixed subtype (subtype B) shared features with both the immune-activating and stroma-proliferating subtypes, whereas upregulated DEGs of these subtypes exhibited no overlap. As it was a blend of the two subtypes, very few genes were unique to the mixed subtype.

Figure 4

Figure 4 Gene expression patterns of psoriasis skin subtypes. (A, B) Molecular pattern distribution of subtypes A and C concerning different biological processes. The top 20 most significant biological processes in the GO-BP database. (C–F) Molecular pattern distribution of subtypes (A, C) concerning different pathways. Top 5 most considerable signaling pathways in the KEGG and Reactome databases. (G) A Venn diagram shows upregulated DEGs in subtypes (A, C).

3.3 Molecular and cellular characterization of the three skin subtypes

There were three clustered subtypes: A (n=121), B (n=94), and C (n=35). To validate whether immune and stromal characteristics differ among these subtypes, we compared enrichment scores of important psoriasis-associated pathways and cell subsets. Twelve psoriasis-associated pathways or processes were curated from the literature and KEGG and Reactome databases. There were significant differences among these subtypes in enrichment scores for these signaling pathways associated with psoriasis. In subtype A, NLR, TLR, RLR, p53, IFN-α/β, IFN-γ, IL-17, IL-23, T-cell receptor, and B-cell receptor signaling were strongly activated, indicating immune activation. Subtype C was characterized by stroma proliferation and enrichment in the Wnt and PPAR signaling pathways, whereas all signaling pathways are moderately activated in subtype B (Figure 5A).

Figure 5

Figure 5 Pathway and cell subset-driven characterization of psoriasis skin subtypes. (A, B) Box plots for enrichment scores of pathways and xCell-inferred cell subsets for each psoriasis skin subtype. Differences across the three subtypes were analyzed using the Wilcoxon and Kruskal-Wallis test. ns, not significant; *P<0.05; **P<0.01; ***P<0.001; **** P<0.0001.

The xCell software was used to identify cell types leading to gene expression discrepancy among subtypes, and a machine-learning framework was built to estimate cell type enrichment. All three subgroups showed differential activation of immune and stromal cells, consistent with these expression patterns. Immune cells [including neutrophils, basophils, macrophages, plasmacytoid dendritic cells (pDCs), and monocytes] and stroma cells (including epithelial cells, keratinocytes, and sebocytes) were more activated in subtype A than in subtypes B and C, whereas mast cell activation declined sharply. Moreover, fibroblasts, adipocytes, and endothelial cells were significantly infiltrated in subtype C (Figure 5B). An analysis of functional pathways in subtypes of psoriasis confirmed these results.

3.4 Prediction of the treatment responses of the subgroups

An Xgboost machine learning algorithm developed a 163-gene classifier to verify the psoriasis subtyping scheme. After training the classifier with 176 psoriasis samples, we applied it to the testing set containing 74 psoriasis samples to confirm its practicality. We observed that the AUC of the training set is 0.932, which proves that the model effectively classifies it. Furthermore, the classifier demonstrated robustness by achieving an average classification performance of 0.905 in the testing set. As a consequence, the classifier was highly effective and applicable in assessing psoriasis skin subtypes.

To determine whether the therapeutic response to etanercept, brodalumab, methotrexate (MTX) and ustekinumab were associated with subtypes of psoriasis, we classified patients at baseline using a fitted Xgboost classifier in four datasets. In the etanercept and brodalumab treated groups, the proportion of good responders was significantly higher in subtype C [9/10 (90.0%) and 10/10 (100.0%)] than in subtypes A [25/39 (64.1%) and 26/33 (78.8%)] and B [11/17 (64.7%) and 25/30 (83.3%); Figures 6A, B]. In contrast, an opposite trend was observed in the MTX-treated group; the proportion of good responders was higher in subtypes A [3/8 (37.5%)] and B [1/3 (33.3%)] than in subtype C [0/4 (0.0%); Figure 6C]. An excellent response to ustekinumab was more frequently observed in subtypes A [32/43 (74.4%)] and B [26/31 (83.9%)] than in subtype C [6/11 (54.5%; Figure 6D]. In spite of this, the differences were not statistically significant, likely due to the limited sample size. In our analysis, subtypes A and B exhibited good responses to MTX and IL-12/23 inhibitors (i.e., ustekinumab) and inadequate responses to TNF-α inhibitors (etanercept) and IL-17A receptor (IL-17RA) inhibitors (i.e., brodalumab) compared with subtype C. Collectively, psoriasis molecular subtyping could have an impact on the benefits of drug treatment. Future clinical studies should integrate this information.

Figure 6

Figure 6 Treatment response of each psoriasis skin subtype. Distribution of psoriasis skin subtypes based on the treatment response of patients treated with (A) etanercept, (B) brodalumab, (C) methotrexate and (D) ustekinumab. Response: responded to the biologics; non-response: did not respond to the biologics.

4 Discussion

Studying the molecular features of skin stratification in psoriasis has improved our understanding of its biological complexity and clinical heterogeneity and promoted research on psoriasis. Although previous studies have emphasized that stromal-cell-rich subtypes may respond well to existing biological therapies (9, 10), adequate information on approaching treatment-resistant patients with psoriasis remains lacking. Here, we classified the psoriasis skin tissues based on their molecular signatures using unsupervised cluster analysis (31). We characterized the different features of the three clustered subgroups in cellular components and biological processes and explained the results from therapeutic perspectives. In detail, subtype A (named immune activation subtype) had a transcriptomic signature of immune cells and proinflammatory activation-related pathways. In contrast, subtype C (designated as the stroma proliferation subtype) exhibited more enriched transcripts in stroma cells (such as fibroblasts and endothelial cells) and tissue proliferative-related pathways, while the features of subtype B (named as mixed subtype) fell somewhere in between. Notably, the therapeutic responses were different among these three subgroups.

Currently, 20-30% of patients do not respond to FDA-approved treatments for psoriasis (34). Therefore, correct disease stratification and active exploration of the response of Psoriasis patients to different treatments are the first steps in the precision treatment of Psoriasis. In the past few years, Wang et al. analyzed genome-wide mRNA expression in psoriatic skin biopsies, discovered different immune cell infiltration patterns that distinguish psoriatic lesions from healthy skin, and successfully classified psoriasis into two different immune phenotypes. Among them, the nonlesional phenotype was characterized by overexpression of genes involved in keratinocyte differentiation and referred to as associated with better treatment response of biologics (10). This conclusion is highly consistent with our findings but has not provided an adequate explanation and address on how to approach it. Our results indicate that subtype C was strongly enriched with proliferative tissue pathways, including PPAR and Wnt signaling, and cellular components, such as fibroblasts and endothelial cells, but not for immune cells and proinflammatory signaling pathways. This finding indicates that immune cell activation was not necessarily involved in skin destruction and may provide a reasonable explanation for the conflict between the low level of inflammatory markers and continuous disease progression on the part of patients with psoriasis. Patients with subtype C presented good responses to TNF-α inhibitors (such as etanercept) and IL-17RA inhibitors (brodalumab). Previous studies have shown that a small amount of TNF-α and IL-17 can act on fibroblasts and endothelial cells to stimulate and produce numerous proinflammatory and proliferative cytokines, such as IL-6, IL-8, and CXCL-1, which cause excessive proliferation of psoriatic skin epidermis (35). Notably, Zaba et al. found that the efficacy of TNF-α inhibitors was associated with IL-17RA inhibition, as rapid downregulation of IL-17A pathway-related gene expression was observed in patients who responded to etanercept but not in the non-responders (36). These results suggest that the current therapeutic strategies targeting TNF-α and IL-17RA could almost wholly inhibit the skin pathology in the good responders. Interestingly, brodalumab’s effects on fibroblasts and endothelial cells have been demonstrated in systemic sclerosis (37). Brodalumab reduces fibroblast proliferation and collagen production by maintaining the regulatory T (Treg) cells/T helper (Th) 17 balance (38), resulting in reduced dermal thickness and improvement in modified Rodnan skin score. In addition, Takemichi Fukasawa et al. found that after six months of guselkumab (IL-23 inhibitor) being used to treat patients with psoriasis vulgaris complicated by systemic sclerosis, both Th2 and Th17 cells showed a decline, and the severity of the disease was also significantly improved (39). These findings further confirm that the IL-23/IL-17 axis is an essential pathway for targeted therapy of inflammatory skin diseases.

Activation of the inflammatory profile of subtypes A may account for the positive results of current psoriasis treatment targeting the upstream immune-inflammatory response. Patients with subtypes A were relatively sensitive to disease-modifying anti-rheumatic drugs (DMARDs) (i.e., MTX) and IL-12/23 inhibitors (i.e., ustekinumab). The anti-microbial peptide LL-37 is primarily secreted by keratinocytes and immune cells in the early stages of psoriasis development through direct activation of pDCs and myeloid dendritic cells to secrete IL-12 and IL-23 (40–42). Moreover, IL-12 and IL-23 stimulated TNF-α and IL-17 secretion in Th1 and Th17 cells (43, 44), resulting in a strong skin inflammatory reaction. Ustekinumab was designed to prevent the proliferation of the series of immune cells and the secretion of various proinflammatory cytokines by blocking IL-12/23 (45), an upstream target of the inflammatory signaling pathway enriched in subtypes A and B. The effects of MTX therapy on molecular signatures were largely restricted to proinflammatory pathways and immune cell subsets. In addition, MTX has recently been discovered to inhibit the JAK/STAT pathway, while many of its side effects are likely to be related to the folate pathway (46). The IL-23 receptor relies on a heterodimer of Janus kinase 2 (JAK2) and tyrosine kinase 2 (TYK2) for signal transduction, thereby highlighting the role of JAKs in the therapeutical potential of JAK inhibitors (47, 48). Moreover, Ishizaki et al. observed that TYK2-deficient mice injected with IL-23 showed significantly reduced ear erythema and epidermal hyperplasia compared to wild-type mice (49). A lack of TYK2 also impaired the infiltration of various immune cells into the skin and the production of IL-17 and IL-22. The JAK/STAT and IL-12/23 pathways have similar contributions to the progression of psoriasis. Thus, it ought to be prudent to research whether the loss of efficacy or delayed treatment resistance can be attributed to the conversion of skin molecular patterns.

Nevertheless, the weak response to TNF-α and IL-17RA inhibitors in patients with subtypes A compared with that for subtype C seemed to contradict their enrichment of numerous inflammatory pathways and underlying pathophysiology. This might be attributed to a paradoxical reaction, in which the disease worsens during treatment with targeted biological drugs. Two main hypotheses have been proposed for this contradictory response (50). One theory is the involvement of TNF-α and IFN-α cross-regulation. Palucka et al. demonstrated that blocking endogenous TNF-α results in increased IFN-γ production by pDC and subsequent T-cell activation resulting in increased TNF-α production (51). The skin lesions in patients with psoriasis induced by TNF-α inhibitors are characterized by IFN-α overexpression compared with those in patients with psoriasis vulgaris. Another hypothesis posits that following TNF-α inhibition, Th17 cells are enhanced and regulatory T cells are downregulated, leading to increased production of the Th17 cytokine IL-22 (52). Both pathways produce cytokines that act on keratinocytes and generate a proinflammatory cycle, leading to suboptimal treatment of psoriasis. This suggests that although targeted therapy for specific cytokines has achieved good efficacy, some people may experience negative feedback after long-term targeted therapy, leading to the loss of treatment effectiveness. Moreover, inhibiting single downstream cytokines (TNF-α and IL-17) to prevent inflammatory response and reverse the disease outcome entirely is challenging. Thus, active intervention in the dysregulation of upstream targets (IL-12/23) is warranted.

Our study has some limitations. First, this study was conducted in different clinical environments from different public datasets; more meta-data is required, but this may prove challenging. Second, owing to the lack of complete annotation for each psoriasis sample, we could not address the association of each psoriasis subtype with clinical factors, such as autoantibody levels. Finally, these differences across responders and non-responders of subtypes were not statistically significant owing to the limited sample size of the clinical response in patients with psoriasis.

5 Conclusions

In summary, we deconvoluted psoriasis skin tissues into pathologically discrete subsets by combining skin transcriptomic data with machine learning algorithms. We described their different molecular and cellular characteristics regarding treatment response. Our results provide critical insights into distinct and shared mechanistic features of psoriasis to determine the pathobiological approaches for benefiting drug-resistant patients.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding author.

Ethics statement

This study used publicly available data from the Gene Expression Omnibus (GEO) database. All subjects in the original studies have provided informed consent, so further ethical approval is not required.

Author contributions

SZ: Writing – original draft. MC: Writing – original draft. LZhe: Writing – original draft. CanW: Writing – review & editing. RZ: Writing – review & editing. SS: Writing – review & editing. JH: Writing – review & editing. LZha: Writing – review & editing. CaiW: Writing – review & editing. XL: Writing – review & editing.

Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article. This work was supported by the National Natural Science Foundation of China (No. 82001740).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2024.1326502/full#supplementary-material

Supplementary Image 1 | Principal component analysis (PCA) plots with or without the elimination of batch effects. (A) PCA before batch effect adjustment for the training microarray datasets. Samples from the different datasets cluster together. (B) PCA after batch effect adjustment for the training microarray datasets. Samples from different datasets overlap.

Abbreviations

DEGs, differentially expressed genes; MTX, methotrexate; PASI, Psoriasis Area and Severity Index; PCA, Principal components analysis; PPI, protein–protein interaction; MCODE, Molecular Complex Detection; GSEA, Gene-set enrichment analysis; KEGG, Kyoto Encyclopedia of Genes and Genomes; CDF, cumulative distribution function; IFN, type I interferon; IL, interleukin.

References

1. Kim WB, Jerome D, Yeung J. Diagnosis and management of psoriasis. Can Fam Physician. (2017) 63:278–85.

PubMed Abstract | Google Scholar

2. Rendon A, Schäkel K. Psoriasis Pathogenesis and Treatment. Int J Mol Sci. (2019) 20:1475. doi: 10.3390/ijms20061475

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Griffiths CE, Barker JN. Pathogenesis and clinical features of psoriasis. Lancet. (2007) 370:263–71. doi: 10.1016/S0140-6736(07)61128-3

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Boehncke WH, Schön MP. Psoriasis. Lancet. (2015) 386:983–94. doi: 10.1016/S0140-6736(14)61909-7

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Armstrong AW, Read C. Pathophysiology, clinical presentation, and treatment of psoriasis: A review. Jama. (2020) 323:1945–60. doi: 10.1001/jama.2020.4006

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Tong Y, Peranteau AJ, Nawas Z, Tyring SK. A review of brodalumab, an IL-17 receptor antagonist, for moderate-to-severe plaque psoriasis. Skin Ther Lett. (2017) 22:1–6.

Google Scholar

7. Armstrong AW, Puig L, Joshi A, Skup M, Williams D, Li J, et al. Comparison of biologics and oral treatments for plaque psoriasis: A meta-analysis. JAMA Dermatol. (2020) 156:258–69. doi: 10.1001/jamadermatol.2019.4029

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Bai F, Li GG, Liu Q, Niu X, Li R, Ma H. Short-term efficacy and safety of IL-17, IL-12/23, and IL-23 inhibitors brodalumab, secukinumab, ixekizumab, ustekinumab, guselkumab, tildrakizumab, and risankizumab for the treatment of moderate to severe plaque psoriasis: A systematic review and network meta-analysis of randomized controlled trials. J Immunol Res. (2019) 2019:2546161. doi: 10.1155/2019/2546161

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Ainali C, Valeyev N, Perera G, Williams A, Gudjonsson JE, Ouzounis CA, et al. Transcriptome classification reveals molecular subtypes in psoriasis. BMC Genomics. (2012) 13:472. doi: 10.1186/1471-2164-13-472

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Wang G, Miao Y, Kim N, Sweren E, Kang S, Hu Z, et al. Association of the psoriatic microenvironment with treatment response. JAMA Dermatol. (2020) 156:1057–65. doi: 10.1001/jamadermatol.2020.2118

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Suárez-Fariñas M, Li K, Fuentes-Duculan J, Hayden K, Brodmerkel C, Krueger JG. Expanding the psoriasis disease profile: interrogation of the skin and serum of patients with moderate-to-severe psoriasis. J Invest Dermatol. (2012) 132:2552–64. doi: 10.1038/jid.2012.184

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Nair RP, Duffin KC, Helms C, Ding J, Stuart PE, Goldgar D, et al. Genome-wide scan reveals association of psoriasis with IL-23 and NF-kappaB pathways. Nat Genet. (2009) 41:199–204. doi: 10.1038/ng.311

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Yao Y, Richman L, Morehouse C, de los Reyes M, Higgs BW, Boutrin A, et al. Type I interferon: potential therapeutic target for psoriasis? PloS One. (2008) 3:e2737. doi: 10.1371/journal.pone.0002737

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Bigler J, Rand HA, Kerkof K, Timour M, Russell CB. Cross-study homogeneity of psoriasis gene expression in skin across a large expression range. PloS One. (2013) 8:e52242. doi: 10.1371/journal.pone.0052242

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Russell CB, Rand H, Bigler J, Kerkof K, Timour M, Bautista E, et al. Gene expression profiles normalized in psoriatic skin by treatment with brodalumab, a human anti-IL-17 receptor monoclonal antibody. J Immunol. (2014) 192:3828–36. doi: 10.4049/jimmunol.1301737

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Kim J, Oh CH, Jeon J, Baek Y, Ahn J, Kim DJ, et al. Molecular Phenotyping Small (Asian) versus Large (Western) Plaque Psoriasis Shows Common Activation of IL-17 Pathway Genes but Different Regulatory Gene Sets. J Invest Dermatol. (2016) 136:161–72. doi: 10.1038/JID.2015.378

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Johnston A, Guzman AM, Swindell WR, Wang F, Kang S, Gudjonsson JE. Early tissue responses in psoriasis to the antitumour necrosis factor-α biologic etanercept suggest reduced interleukin-17 receptor expression and signalling. Br J Dermatol. (2014) 171:97–107. doi: 10.1111/bjd.12937

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Swindell WR, Xing X, Stuart PE, Chen CS, Aphale A, Nair RP, et al. Heterogeneity of inflammatory and cytokine networks in chronic plaque psoriasis. PloS One. (2012) 7:e34594. doi: 10.1371/journal.pone.0034594

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Brodmerkel C, Li K, Garcet S, Hayden K, Chiricozzi A, Novitskaya I, et al. Modulation of inflammatory gene transcripts in psoriasis vulgaris: Differences between ustekinumab and etanercept. J Allergy Clin Immunol. (2019) 143:1965–9. doi: 10.1016/j.jaci.2019.01.017

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Tomalin LE, Russell CB, Garcet S, Ewald DA, Klekotka P, Nirula A, et al. Short-term transcriptional response to IL-17 receptor-A antagonism in the treatment of psoriasis. J Allergy Clin Immunol. (2020) 145:922–32. doi: 10.1016/j.jaci.2019.10.041

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Ducreux J, Houssiau FA, Vandepapelière P, Jorgensen C, Lazaro E, Spertini F, et al. Interferon α kinoid induces neutralizing anti-interferon α antibodies that decrease the expression of interferon-induced and B cell activation associated transcripts: analysis of extended follow-up data from the interferon α kinoid phase I/II study. Rheumatol (Oxford). (2016) 55:1901–5. doi: 10.1093/rheumatology/kew262

CrossRef Full Text | Google Scholar

22. Leek JT, Johnson WE, Parker HS, Jaffe AE, Storey JD. The sva package for removing batch effects and other unwanted variation in high-throughput experiments. Bioinformatics. (2012) 28:882–3. doi: 10.1093/bioinformatics/bts034

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. (2015) 43:e47. doi: 10.1093/nar/gkv007

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Hochberg Y, Benjamini Y. More powerful procedures for multiple significance testing. Stat Med. (1990) 9:811–8. doi: 10.1002/sim.4780090710

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Zhou Y, Zhou B, Pache L, Chang M, Khodabakhshi AH, Tanaseichuk O, et al. Metascape provides a biologist-oriented resource for the analysis of systems-level datasets. Nat Commun. (2019) 10:1523. doi: 10.1038/s41467-019-09234-6

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Szklarczyk D, Franceschini A, Wyder S, Forslund K, Heller D, Huerta-Cepas J, et al. STRING v10: protein-protein interaction networks, integrated over the tree of life. Nucleic Acids Res. (2015) 43:D447–52. doi: 10.1093/nar/gku1003

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Oughtred R, Rust J, Chang C, Breitkreutz BJ, Stark C, Willems A, et al. The BioGRID database: A comprehensive biomedical resource of curated protein, genetic, and chemical interactions. Protein Sci. (2021) 30:187–200. doi: 10.1002/pro.3978

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Bader GD, Hogue CW. An automated method for finding molecular complexes in large protein interaction networks. BMC Bioinf. (2003) 4:2. doi: 10.1186/1471-2105-4-2

CrossRef Full Text | Google Scholar

29. Canzler S, Hackermüller J. multiGSEA: a GSEA-based pathway enrichment analysis for multi-omics data. BMC Bioinf. (2020) 21:561. doi: 10.1186/s12859-020-03910-x

CrossRef Full Text | Google Scholar

30. Subramanian A, Kuehn H, Gould J, Tamayo P, Mesirov JP. GSEA-P: a desktop application for Gene Set Enrichment Analysis. Bioinformatics. (2007) 23:3251–3. doi: 10.1093/bioinformatics/btm369

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Wilkerson MD, Hayes DN. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics. (2010) 26:1572–3. doi: 10.1093/bioinformatics/btq170

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinf. (2013) 14:7. doi: 10.1186/1471-2105-14-7

CrossRef Full Text | Google Scholar

33. Aran D, Hu Z, Butte AJ. xCell: digitally portraying the tissue cellular heterogeneity landscape. Genome Biol. (2017) 18:220. doi: 10.1186/s13059-017-1349-1

PubMed Abstract | CrossRef Full Text | Google Scholar

34. van Schouwenburg PA, van de Stadt LA, de Jong RN, van Buren EE, Kruithof S, de Groot E, et al. Adalimumab elicits a restricted anti-idiotypic antibody response in autoimmune patients resulting in functional neutralisation. Ann Rheum Dis. (2013) 72:104–9. doi: 10.1136/annrheumdis-2012-201445

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Furue K, Ito T, Tsuji G, Kadono T, Furue M. Psoriasis and the TNF/IL23/IL17 axis. G Ital Dermatol Venereol. (2019) 154:418–24. doi: 10.23736/S0392-0488.18.06202-8

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Zaba LC, Suárez-Fariñas M, Fuentes-Duculan J, Nograles KE, Guttman-Yassky E, Cardinale I, et al. Effective treatment of psoriasis with etanercept is linked to suppression of IL-17 signaling, not immediate response TNF genes. J Allergy Clin Immunol. (2009) 124:1022–10.e1-395. doi: 10.1016/j.jaci.2009.08.046

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Fukasawa T, Yoshizaki A, Ebata S, Fukayama M, Kuzumi A, Norimatsu Y, et al. Interleukin-17 pathway inhibition with brodalumab in early systemic sclerosis: Analysis of a single-arm, open-label, phase 1 trial. J Am Acad Dermatol. (2023) 89:366–9. doi: 10.1016/j.jaad.2023.02.061

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Chizzolini C, Dufour AM, Brembilla NC. Is there a role for IL-17 in the pathogenesis of systemic sclerosis? Immunol Lett. (2018) 195:61–7. doi: 10.1016/j.imlet.2017.09.007

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Fukasawa T, Yoshizaki-Ogawa A, Yoshizaki A, Sato S. Impact of guselkumab on three cases of SSc accompanying psoriasis. Rheumatol (Oxford). (2023) 63:e6–e8. doi: 10.1093/rheumatology/kead287

CrossRef Full Text | Google Scholar

40. Zasloff M. Antimicrobial peptides of multicellular organisms. Nature. (2002) 415:389–95. doi: 10.1038/415389a

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Lande R, Gregorio J, Facchinetti V, Chatterjee B, Wang YH, Homey B, et al. Plasmacytoid dendritic cells sense self-DNA coupled with antimicrobial peptide. Nature. (2007) 449:564–9. doi: 10.1038/nature06116

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Ganguly D, Chamilos G, Lande R, Gregorio J, Meller S, Facchinetti V, et al. Self-RNA-antimicrobial peptide complexes activate human dendritic cells through TLR7 and TLR8. J Exp Med. (2009) 206:1983–94. doi: 10.1084/jem.20090480

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Nestle FO, Conrad C, Tun-Kyi A, Homey B, Gombert M, Boyman O, et al. Plasmacytoid predendritic cells initiate psoriasis through interferon-alpha production. J Exp Med. (2005) 202:135–43. doi: 10.1084/jem.20050500

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Santini SM, Lapenta C, Donati S, Spadaro F, Belardelli F, Ferrantini M. Interferon-α-conditioned human monocytes combine a Th1-orienting attitude with the induction of autologous Th17 responses: role of IL-23 and IL-12. PloS One. (2011) 6:e17364. doi: 10.1371/journal.pone.0017364

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Yiu ZZ, Warren RB. Ustekinumab for the treatment of psoriasis: an evidence update. Semin Cutan Med Surg. (2018) 37:143–7. doi: 10.12788/j.sder.2018.040

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Alqarni AM, Zeidler MP. How does methotrexate work? Biochem Soc Trans. (2020) 48:559–67. doi: 10.1042/BST20190803

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Subhadarshani S, Yusuf N, Elmets CA. IL-23 and the tumor microenvironment. Adv Exp Med Biol. (2021) 1290:89–98. doi: 10.1007/978-3-030-55617-4_6

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Kvist-Hansen A, Hansen PR, Skov L. Systemic treatment of psoriasis with JAK inhibitors: A review. Dermatol Ther (Heidelb). (2020) 10:29–42. doi: 10.1007/s13555-019-00347-w

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Ishizaki M, Muromoto R, Akimoto T, Sekine Y, Kon S, Diwan M, et al. Tyk2 is a therapeutic target for psoriasis-like skin inflammation. Int Immunol. (2014) 26:257–67. doi: 10.1093/intimm/dxt062

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Miyagawa F. Pathogenesis of paradoxical reactions associated with targeted biologic agents for inflammatory skin diseases. Biomedicines. (2022) 10:1485. doi: 10.3390/biomedicines10071485

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Palucka AK, Blanck JP, Bennett L, Pascual V, Banchereau J. Cross-regulation of TNF and IFN-alpha in autoimmune diseases. Proc Natl Acad Sci U S A. (2005) 102:3372–7. doi: 10.1073/pnas.0408506102

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Ma HL, Napierata L, Stedman N, Benoit S, Collins M, Nickerson-Nutter C, et al. Tumor necrosis factor alpha blockade exacerbates murine psoriasis-like disease by enhancing Th17 function and decreasing expansion of Treg cells. Arthritis Rheumatol. (2010) 62:430–40. doi: 10.1002/art.27203

CrossRef Full Text | Google Scholar