Integrative Transcriptomic and Single-Cell Analysis Reveals the Role of Autophagy-Related Genes in Psoriasis Pathogenesis

Introduction

Psoriasis is a chronic inflammatory skin disease characterized by abnormal proliferation and differentiation of keratinocytes (KCs) and over-infiltration of immunocytes in the epidermis and dermis. Clinically, psoriasis presents in various forms, primarily classified as psoriasis vulgaris, psoriatic arthritis, psoriasis erythrodermic and generalized pustular psoriasis, with psoriasis vulgaris being the most prevalent.^1,2 Globally, the estimated incidence of psoriasis is 2–3%, and it often coexists with other systemic inflammatory conditions such as depression, arthritis, and cardiovascular disease.^3,4 Notably, psoriatic patients, particularly those with moderate to severe disease, exhibit higher mortality rates and a reduced life expectancy of approximately 5 years compared to healthy controls.^5,6 Despite the relative ease of diagnosis, psoriasis currently lacks an effective treatment, necessitating long-term management and significantly impacting the quality of life.⁷ While research into the development and progression of psoriasis has advanced, the precise molecular mechanisms remain elusive. Therefore, the development of new therapeutic targets and the identification of potential biomarkers are pivotal for enhancing the assessment and treatment of psoriasis.

Autophagy, distinct from the ubiquitin-proteasome system, is the only known conserved protein-degradation pathway. It encompasses macroautophagy, chaperone-mediated autophagy, and microautophagy.^8,9 Autophagy and related markers are implicated in a wide spectrum of human diseases, including those affecting the lungs, liver, heart, as well as neurodegenerative disorders, myopathies, cancer, aging, and metabolic conditions.¹⁰ Notably, autophagy has been identified as a negative regulator of pathways crucial to psoriasis pathogenesis, such as TLR2/6-mediated NF-κB activation, SQSTM1 expression, and inflammatory cytokine secretion in Keratinocytes.¹¹ The autophagy-associated protein SQSTM1 is overexpressed in the epidermis of psoriatic skin,¹¹ and ATG5-dependent autophagy is considered a pivotal factor in psoriasis development and a potential target for treatment.¹² Additionally, genetic variations in the autophagy gene ATG16L1 have been linked to psoriasis susceptibility,¹³ and certain psoriasis treatments are known to induce autophagy.¹⁴ Increasing evidence suggests that autophagy and related proteins are involved in autoimmune diseases, including psoriasis, where autophagy plays critical roles in immune function such as pathogen clearance, secretory pathways, lymphocyte development, and pro-inflammatory signaling.¹⁵ Therefore, modulating autophagy may represent a promising therapeutic strategy for psoriasis, warranting further investigation.

Bioinformatics has been extensively applied in the study of psoriasis, leveraging databases, such as the Gene Expression Omnibus (GEO), to uncover key regulatory genes. In recent years, numerous studies have utilized bioinformatic approaches to mine these datasets to gain insights into the genetic underpinnings of psoriasis.¹⁶

By combining single-cell and mRNA data, our study provides a more comprehensive view of the complex biological processes of psoriasis, and sheds light on potential new targets for therapeutic interventions. To strengthen the reliability of our findings, we performed reverse transcription-quantitative polymerase chain reaction.

(RT-qPCR) validation in clinical cohorts and animal models translates transcriptomic screening results into reproducible evidence and preliminary exploration of their in vivo function, providing insights for the development of improved diagnostic and therapeutic strategies for this challenging disease.

Materials and Methods

Data Download and Data Preprocessing

To investigate the biological attributes of psoriasis, computational biology methodology was adopted. All pertinent datasets were retrieved from GEO (https://www.ncbi.nlm.nih.gov/geo/), which encompasses genome-wide transcriptional profiles originating from either psoriatic or healthy control skin tissues of Homo sapiens. This included the acquisition of 11 samples of psoriatic skin lesions and 10 normal skin controls from the single-cell RNA sequencing (scRNA-seq) dataset GSE220116. Furthermore, 92 psoriasis cases and 82 control samples were extracted from the RNA-Seq dataset GSE54456 and sequenced on the GPL9052 Illumina Genome Analyzer platform for Homo sapiens. The primary dimensional matrix data were annotated with gene identifiers and underwent quality control measures employing the “limma” R package in anticipation of downstream analysis. A compendium of 358 autophagy-linked genes was sourced from the Molecular Signatures Database (MsigDB) (http://www.gseamsigdb.org/gsea/msigdb/index.jsp),¹⁷ as detailed in Supplementary Table 1).

Download and Processing of scRNA-Seq Data

The scRNA-seq dataset GSE220116 was accessed and incorporated into the analysis pipeline through utilization of the “Seurat” R package, incorporating 11 psoriatic lesion and 10 normal skin control samples. Low-quality cells and genes were excluded based on rigorous criteria: cells expressing fewer than 200 genes, those with a mitochondrial gene ratio of over 20%, those expressing more than 5000 genes, erythrocyte content exceeding 10%, and those with UMI counts below 50,000 were retained. Subsequently, data normalization was executed using the “NormalizeData” function within “Seurat”. Following normalization, highly variable genes in individual cells were identified by balancing mean expression and dispersion. Principal Component Analysis (PCA) was implemented, and principal components (PCs) were fed into graph-based clustering. Sample batch effects were mitigated using a harmony algorithm. Unsupervised and unbiased clustering was achieved via a shared nearest neighbor (SNN)-based algorithm, “FindNeighbors”, which yielded 20 clusters across 20 PCs at a resolution of 0.8. Dimensionality reduction was further accomplished with t-Distributed Stochastic Neighbor Embedding (t-SNE) utilizing “RunTSNE”, visualizing cell clusters through t-SNE-1 and t-SNE-2 coordinates. Cell clusters were annotated using cell type-specific biomarkers, and their proportions were quantified and assessed.

Autophagy-Related Gene Scores

For assessing autophagy-related activity in each cell from the scRNA-seq dataset, the “AUCell” R package was employed.¹⁸ By calculating the area under the curve (AUC) for autophagy-related genes, a feature-specific score for each cell was derived, generating a ranked list of gene expression per cell, facilitating the identification of cell types with the most prominent autophagy signatures. Thresholds defining active gene sets were established with “AUCell_exploreThresholds”. The AUC scores were then projected onto t-SNE embeddings using the “ggplot2” R package, visually highlighting clusters of autophagy-active cells.

Pseudo-Time Analysis to Construct Cell Trajectories

Pseudo-time ordering of cells was carried out using the “Monocle2” R package,¹⁹ generating a pseudo-time graph that elucidated branching and linear differentiation processes. For cells demonstrating heightened autophagic activity, the raw counts were normalized prior to trajectory inference. Pseudotime trajectories were constructed focusing on genes with high discretization and elevated expression levels (discretization estimate ≥ 1, mean expression ≥ 0.1).²⁰ The “DDRTree” algorithm was applied using default parameters. Branching events were meticulously analyzed using Branching Expression Analysis Modeling (BEAM) in Monocle 2, aiding in the recognition of genes displaying significant branch-dependent expression patterns,¹⁹ represented visually as heat maps.

Cell Communication Analysis and Ligand Receptor Expression

Investigation of cell-cell communications (CCCs) entailed examining the expression of ligand-receptor pairs and associated signaling cascades.²¹ “CellChat” R package was utilized to establish CellChat objects based on UMI count matrices for both psoriasis and control groups (https://www.github.com/sqjin/CellChat).²² Employing the “CellChatDB.human” interaction database as a reference, and intercellular communication for various cell types in psoriatic skin samples was quantified using default parameters. “mergeCellChat” function combined CellChat objects across groups to compare interaction frequencies and intensities. “netVisual_diffInteraction” function visualized variations in interaction counts or strengths across groups and cell types. Lastly, “netVisual_bubble” and “netVisual_aggregate” functions were used to illustrate the signaling gene expression distribution between groups, revealing key signaling pathway activities in psoriatic cells and selecting specific pathways for further inspection.

Gene Set Variation Analysis

To explore functional discrepancies between psoriasis patients and controls, Gene Set Variation Analysis (GSVA), an unsupervised, non-parametric approach for gene set enrichment, was implemented using the “c2.cp.kegg.v7.5.1. symbols” gene set from MsigDB (http://software.broadinstitute.org/gsea/msigdb) with the “GSVA” R package.

Screening of Differentially Expressed Genes

Differential gene expression analysis was executed with the “limma” R package,²³ identifying differentially expressed genes (DEGs) between the control group (n=82) and the psoriasis group (n=92), adopting a threshold of |log2Fold Change| > 1.5 and p-value < 0.05 for further examination. Heatmaps were produced by “pheatmap” R package, organized by Euclidean distance and hierarchical clustering. To decipher the expression profiles of candidate genes, autophagy-related differential genes in KC compared to other cells from the psoriasis scRNA-seq dataset and DEGs between the psoriasis and control groups were identified, and their overlap was depicted in a Venn diagram.

Protein-Protein Interaction Network Construction

To elucidate protein interaction dynamics, we forecasted protein-protein interactions (PPIs) employing the STRING database’s (http://string-db.org) online utility, focusing on gene pairs with a confidence score exceeding 0.4, constructing PPI networks for our gene cohort of interest.²⁴ Vital genes, characterized by the highest intra-genetic connectivity, were subsequently harvested utilizing the CytoHubba plugin within Cytoscape software. This process facilitated the detection and visualization of core genes within the PPI landscape. Specifically, we isolated the uppermost 8 genes based on interaction scores using CytoHubba’s maximal clique centrality (MCC) methodology, which was then expanded to include 20 first-order interacting genes, forming a candidate set of 28 genes. Seven pivotal key genes were identified based on comprehensive consideration of network centrality, significant differential expression, and preferential selection of genes not previously associated with psoriasis. Finally, the violin plot (by rank sum test analysis) further illustrated the differences in the comparative expression patterns between patients with psoriasis and controls.

Functional Enrichment Analysis

Function-centric enrichment evaluations were conducted on differentially expressed autophagy genes that distinguished between the autophagy-active and autophagy-inactive cell populations. This was achieved by leveraging the “FindAllMarkers” function in the “Seurat” Seurat R package, adopting a filter of |log2Fold Change| > 0.25 and p-value < 0.05. To decipher the biological implications, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were carried out through the “clusterProfiler” R package.²⁵ GO enrichment, a staple for large-scale functional annotation,²⁶ encompassed biological process (BP), molecular function (MF), and cellular component (CC) dimensions. KEGG, a repository housing genomic, pathway, disease, and pharmacological insights,²⁷ was instrumental in uncovering the significantly modulated metabolic pathways among the studied genes. Additionally, Gene Set Enrichment Analysis (GSEA) was deployed to probe the biological process discrepancies between patient groups in psoriasis expression profiles, capitalizing on its ability to statistically differentiate gene sets across distinct biological states, thereby shedding light on expression dataset-driven biological process shifts.

GeneMANIA

The GeneMANIA platform (http://genemania.org) was instrumental in predicting functional linkages among genes resembling the central genes, including of protein-protein and protein-DNA interactions, along with co-expression, co-localization, and physiological responses,²⁸ which guided the assembly of a PPI network centered around our key genes.

Construction and Analysis of mRNA-TF Regulatory Network

In the development and scrutiny of the mRNA-Transcription Factor (TF) regulatory framework, The manually curated TRRUST database (http://www.grnpedia.org/trrust/) plays a pivotal role in the development and scrutiny of the mRNA–TF regulatory framework. This database, archiving transcriptional regulatory networks completed with TF-target gene mappings and TF regulatory relationships, allowed us to minimize false-positive occurrences. Our investigation led to the analysis and depiction of a hub gene-TF (mRNA-TFs TF– interaction network using Cytoscape.

Animals and Ethical Statement

Twelve BALB/c mice (6 males and 6 females, 6–8 weeks old) were obtained from Huachuang Sino Co., Ltd. (Taizhou, China). The mice were housed under specific pathogen-free (SPF) conditions in individually ventilated cages (IVCs) with a controlled 12-hour light/dark cycle, temperature of 22 ± 2 °C, and humidity of 50% ± 10%. They were provided with free access to a standard laboratory diet and autoclaved water. All animal experiments were conducted in accordance with the National Institutes of Health (NIH) Guide for the Care and Use of Laboratory Animals and the protocol was approved by the Laboratory Animal Ethics Committee of the Affiliated Hospital of Guangdong Medical University (AHGDMU-LAC-B-202403-0009).

Construction of Animal Models

Adaptive feeding was performed for one week prior to the formal experiments. Twelve mice were randomly divided into the control (n = 6) and experimental (n = 6) groups. A 2×3 cm hairless area on the back was produced by topical application of depilatory cream. Five% imiquimod (IMQ) and Vaseline (62.5 g) were applied to the hairless area in the experimental and control groups every day for 7 days.²⁹ The skin lesions were photographed and scored by PASI. All mice were sacrificed by cervical dislocation on day 8, and skin tissues on the hairless area were preserved in RNAlater solution at −80 °C for RT-qPCR.

Collection of Clinical Samples

Skin samples were collected from 15 patients (14 males and 1 female; 25–57 years) with pathologically diagnosed psoriasis vulgaris. Normal skin samples were obtained from six patients (three males and three females; 9–30 years) with pigmented nevi who underwent surgical excision. The clinical characteristics of the patients are summarized in Table 1. None of the subjects had autoimmune diseases or other chronic inflammatory diseases and used no immunosuppressants, glucocorticoids, or biological agents within 1 month. This study was approved by the Ethics Committee of our hospital (PJKT2025-196) and all participants provided written informed consent. The skin samples were stored at −80 °C until use.

Table 1 Clinical Features of Cases and Controls in the Study

Reverse Transcription Quantitative Polymerase Chain Reaction Analysis

Total RNA was extracted from the collected animal clinical lesion samples using TRIzol reagent (Vazyme). One microgram of total RNA was reverse-transcribed into cDNA using a White Shark Reverse Transcription Kit (White Shark Biosciences, Shanghai, China) in a 20 μL reaction mixture. Quantitative real-time PCR was performed using the ChamQ™SYBRR^®qPCR Premix (High ROX Premix, Vazyme). The β-actin gene was used as an internal standard and the 2−ΔCT method was used to quantify the data. All experiments were repeated three times. All the primers listed in Supplementary Table 2 were synthesized by Guangzhou Aiji Biotechnology Co., Ltd. (Guangzhou, China).

Analysis of Clinical Sample Sequencing Results

The clinical samples collected by RT-qPCR were sent to Shenzhen BGI Co. Ltd. (Shenzhen, China) for transcriptome sequencing. The expression of differentially expressed genes obtained in RT-qPCR experiments was verified by transcriptome sequencing. Differential gene expression analysis was executed with the “limma” R package,²³ identifying differentially expressed genes (DEGs) between the control group (n=6) and the psoriasis group (n=15), adopting a threshold of |log2Fold Change| > 1.5 and p-value < 0.05 for further examination. Heatmaps were produced by “pheatmap” R package, organized by Euclidean distance and hierarchical clustering.

Statistical Analysis

Statistical Analysis for Bioinformatics

Statistical evaluations were performed using the R software version 4.1.2. Spearman correlation test was used to gauge the association between the dual parameters. Comparisons between two groups were performed using the Wilcoxon test, whereas distinctions between three or more groups were assessed using the Kruskal–Wallis test. Statistical significance was determined at p of less than 0.05.

Statistical Analysis of mRNA Expression Data

All mRNA expression data are presented as the mean ± standard error of the mean. GraphPad Prism 7 software was used for data analysis and graphs were generated. The Mann–Whitney U-test was used for comparisons between two groups, and statistical significance was defined as p < 0.05.

The study workflow is illustrated in Figure 1

Figure 1 Study workflow diagram.

Results

Establishing a Cellular Map of Psoriasis Using Single Cell Transcriptomic

To elucidate the origin of the highly expressed transcripts, we examined cellular populations from both psoriasis and healthy control samples using the scRNA-seq dataset, GSE220116. After the initial data quality filtration, a collective of 30,846 cells from single-cell transcriptomes was retained. Our investigation included 21 samples with all cells consolidated into 20 distinct clusters (Figure 2A). By leveraging gene expression patterns across these clusters, we annotated nine unique cell types by associating them with cell-specific biomarkers (Supplementary Table 3), including fibroblasts, KC-stratum basale (SB), KC-stratum corneum (SC), KC-stratum spinosum (SP), mature_DC, melanocytes, NK-cell, semimature_DC, and T-cell. The relative abundance of each cell type in the control and psoriatic samples is shown in Figure 2B. Dot plots in Figure 2C and D illustrate the distribution of cells and cell type-specific genes in the psoriasis cohort, whereas Figure 2E and F show equivalent visualizations of the control samples. The proportions of T cells, KC-SC, and KC-SB varied considerably between normal and psoriatic skin.

Figure 2 Identification of cellular types from scRNA-seq. (A) t-SNE plot visualizing the dispersion of single-cell subsets. (B) Bar graph showing the proportion of diverse cell types in control and psoriasis contexts. (C) t-SNE plot showcasing annotation outcomes for cells in the psoriasis cohort. (D) Expression profiles of marker genes across cell types in psoriasis. (E) t-SNE plot presenting annotation results for control group cells. (F) Marker gene expressions across control group cell types.

Establishment of Cell Subpopulations with Significant Autophagic Activity

To identify the active cells involved in autophagy at single-cell resolution, we analyzed the expression patterns of autophagy-related genes. The application of an optimal threshold criterion revealed 1129 autophagy-linked active cells, characterized by an AUC value surpassing 0.16 (Figure 3A). Figure 3B highlights autophagy-active cells in the t-SNE plot. The autophagic activity of all cells is color-coded in Figure 3C, in which more vivid hues denote higher autophagic activity, indicating that KC-SC have the most pronounced autophagic activity.

Figure 3 Autophagy-associated active cells. (A) Autophagy-related genes’ AUC scores with a threshold of 0.16. (B) t-SNE plot featuring cells with elevated AUC values, colored to distinguish highly autophagic active cells (blue) from inactive ones (gray). (C) t-SNE chromatic map reflecting cell activity scores, with brighter colors signifying higher activity.

Identification of Dynamic Key Pathways During Keratinocytes State Transitions

The reclustering of the KCs defines nine states, and their distributions are illustrated (Figure 4A). The pseudo-time trajectory positioned cluster 0 at the onset and cluster 4 at the terminus. Notably, node 2 experienced a conspicuous shift in the cellular composition (Figure 4B and C). To understand the molecular transitions, we investigated the genes defining keratinocyte branch point 2. GO analysis revealed that pre-branch genes were primarily enriched in skin development and KC differentiation pathways. Branch 2 was enriched with genes involved in keratinization and KC differentiation pathways, whereas branch 1 featured genes related to cytoplasmic translation and regulation of endopeptidase activity (Figure 4D and Supplementary Table 4).

Figure 4 Transcriptional trajectory analysis of keratocytes. (A) Disposition of 9 keratocyte cluster isoforms. (B) Pseudo-time progression from darker to lighter blue. (C) Distribution of keratocyte cluster isoforms along the pseudo-time trajectory. (D) Heatmap presenting differential gene expression across branches, with significantly enriched GO pathways on the left.

GSVA pathway enrichment for each cluster, including the significantly downregulated renin angiotensin system in cluster 0 and NOD-like receptor signaling pathway in cluster 4, is presented in Figure 5A–F and Supplementary Figure 1.

Figure 5 GSVA enrichment outcomes for diverse clusters. GSVA bar graphs depicting divergent pathway enrichments in keratocyte clusters 0 through 5, comparing psoriasis to normal conditions. Blue and green denote upregulated and downregulated pathways, respectively. (A) GSVA bar plot for Cluster 0, the renin-angiotensin system is significantly downregulated in this cluster. (B) GSVA bar plot for Cluster 1. (C) GSVA bar plot for Cluster 2. (D) GSVA bar plot for Cluster 3 (E) GSVA bar plot for Cluster 4, highlighting the significant downregulation of the NOD-like receptor signaling pathway (green bar) alongside other enriched pathways. (F) GSVA bar plot for Cluster 5.

Intercellular Communication Explores Ligand and Receptor Interactions of Keratinocytes

Compared to the controls, the psoriasis group displayed an increase in the number and intensity of cell-type interactions (Figure 6A and B). The communication signaling patterns between the cell types intensified, as shown in Figure 6C and D. The detailed outgoing and incoming signaling patterns for the control and psoriasis groups are shown in Figure 6E and F.

Figure 6 Intercellular communication dynamics. (A) A bar graph depicts the counts of interactions among cell types in control and psoriasis cohorts. (B) A bar graph depicts the intensities of interactions among cell types in control and psoriasis cohorts. (C) Network diagram showing changes in the strength of interactions between cell types in the control and psoriasis groups. (D) Network diagram showing changes in the quantities of interactions between cell types in the control and psoriasis groups. (E) Heatmap illustrating outgoing signaling patterns of cell types in the control group. (F) Heatmap illustrating incoming signaling patterns of cell types in the control group.

In addition, we found that HLA-E signaling was significantly enhanced in psoriasis patients and that HLA-E signaling was enhanced with CD94/NKG2A and attenuated with KLRK1, implying that there may be an interaction between Keratinocytes and NK cells (Figure 7A–D). HLA-E is a non-classical Major Histocompatibility Complex (MHC) Class Ib molecule with limited polymorphism, suggesting that the specific MHC-I pathway is a potential biological target of Keratinocytes with NK cells.

Figure 7 Ligand-receptor dynamics in cellular communication. (A) Ligand-receptor pairs with heightened communication efficiency between other cell types and keratocytes in the psoriasis context. (B) Ligand-receptor pairs with reduced communication efficiency between other cell types and keratocytes in the psoriasis context. (C) Distribution of ligand-receptor expression within the MHC-I signaling pathway in the psoriasis group. (D) Distribution of ligand-receptor expression within the MHC-I signaling pathway in the control group.

Identification of DEGs and Enrichment Pathways Using the mRNA Transcriptome

Contrasting psoriasis and healthy control samples from the mRNA dataset GSE54456 led to the discovery of 209 DEGs, with statistically significant differences (p < 0.05, |Log2 fold change| > 1.5) between the groups. This subset consisted of 134 upregulated and 75 downregulated genes in diseased tissues. A volcano plot (Figure 8A) visually represented all DEGs, and a heatmap (Figure 8B) emphasized the expression profiles of the 11 prominently altered genes.

Figure 8 Enrichment analysis of psoriasis differential genes. (A) A volcano plot visualizes the DEGs distribution. (B) A heatmap emphasizes the top 11 ranked DEGs. GSEA results underscore the significant enrichment of multiple pathways, such as proteasome (C), cytosolic DNA sensing pathway (D), graft versus host disease (E), focal adhesion (F), hypertrophic cardiomyopathy HCM (G), dilated cardiomyopathy (H).

To unravel the underlying mechanisms, GSEA was performed using the MSigDB dataset, leading to the selection of significantly enriched signaling pathways based on Normalized Enrichment Scores (NES) (Supplementary Table 5). Pathways such as the proteasome (NES = 2.3044, adjusted P =0.0322, FDR = 0.0234), cytosolic DNA-sensing pathway (NES = 2.2024, adjusted P =0.0322, FDR = 0.0234), and graft-versus-host disease (NES = 2.1871, adjusted P =0.0322, FDR = 0.0234) (Figure 8C–E), focal adhesion (NES = −1.7635, adjusted P =0.0322, FDR = 0.0234), hypertrophic cardiomyopathy (HCM) (NES = −1.9137, adjusted P =0.0322, FDR = 0.0234), and dilated cardiomyopathy (NES = −2.013, adjusted P =0.0322, FDR = 0.0234) (Figure 8F–H), were significantly enriched in psoriasis.

Combining scRNA-Seq and mRNA to Search for Key Pathways

Differential analysis of KCs from other cells in the scRNA-seq dataset yielded 3317 autophagy activity DEGs. Intersection with the previously identified DEGs resulted in 72 overlapping genes (Figure 9A and Supplementary Table 6). Lollipop diagrams present the GO (Supplementary Table 7) and KEGG (Supplementary Table 8) enrichment analyses of these genes, with BP linked to keratinization and KC differentiation, CC associated with cornified envelopes, and MF featuring RAGE receptor binding activities, among others (Figure 9B). The KEGG pathway enrichment analysis is shown in Figure 9C, and the most significantly enriched pathway was the IL-17 signaling pathway.

Figure 9 Enrichment analysis of intersected genes. (A) Venn diagram of the intersection of autophagy activity differentially expressed genes and DEGs. (B) Lollipop plot showing the top three pathways with the highest significance of BP, CC and MF enriched for the intersecting gene GO. (C) KEGG enrichment analysis.

Identification of Key Hub Genes

PPI networks were constructed using STRING to explore interactions between multiple genes (Figure 10A). CytoHubba was used to identify central genes, and seven hub genes were generated using the MCC algorithm (SPRR2B, KRT6B, PKP3, KRT75, KRT6C, LGI2, and KRT84)(Figure 10B). To place these hubs in the context of their functional modules, we expanded their scope to 20 first-order interacting genes, forming a candidate set of 28 genes (Figure 10C). Interactions between the 12 co-expressed genes and 18 TFs were mapped (Figure 10D). Based on a combination of network centrality, significant differential expression, and preferential selection of genes not previously associated with psoriasis, we refined this list to include seven pivotal genes (SPRR2B, KRT6B, PKP3, KRT75, KRT6C, LGI2, and KRT84) for subsequent experimental validation. Differential expression patterns were verified using the rank-sum test (Figure 10E).

Figure 10 PPI analysis for screening hub genes. (A) PPI network diagram was constructed for the intersected genes. (B) The MCC algorithm obtained the network graph of the 7 genes with the strongest interactions. (C) Candidate Gene Set of Hubs and Interactions. (D) Co-expressed Gene and Transcription Factor Regulatory Network. (E) Box plot demonstrating the difference in the expression levels of hub genes in psoriasis and control samples, and statistical significance of the difference was tested by rank sum test. *p < 0.05,***p < 0.001.

Investigating the Signaling Pathway of Hub Genes

The Signaling pathways of hub genes were analyzed in relation to 50 hallmark pathways (Supplementary Figure 2A), revealing 30 upregulated and 20 downregulated pathways in the psoriasis group compared with the controls (Supplementary Figure 2B). Fifty significant pathways were differentially enriched in both the disease and control groups, and the top five ranked hub genes were significantly correlated with all of these pathways.

Verification of Key Hub Genes in Animal and Clinical Skin Lesion Samples

To verify the effect of animal modeling, PASI score curves for scale, erythema, and skin thickness were generated (Figure 11A–D), and the attached photograph is a representative mouse back photograph on Day 8 (Figure 11E). In the psoriasis animal model, RT-qPCR was used to verify the expression of seven key hub genes, resulting in three genes that were highly expressed and statistically significant relative to the control group: SPRR2B, KRT 6B, PKP3 (Figure 12, p<0.01). RT-qPCR and transcriptome sequencing were used to verify the high expression of KRT6B and PKP3 in clinical samples (KRT6B, p<0.01; PKP3, p <0.05; Figure 13). SPRR2B was not statistically significant according to the analysis. We speculated that the results were not statistically significant owing to the large heterogeneity of individual clinical samples and the small sample size (n=21); however, the gene expression trend was still higher than that of the control group. This trend was consistent with our expectations.

Figure 11 Evaluation of the Psoriasis Animal Model Establishment. (A) Total PASI score curve of the psoriasis animal model over time. (B) Scaling score curve over time. (C) Erythema score curve over time. (D) Skin thickness score curve over time. (E) Representative photographs of the mouse dorsal skin at different time points (Day 1, Day 3, Day 5, Day 7, and Day 8) post-modeling.

Figure 12 RT-qPCR validation of hub genes in the animal model. SPRR2B, KRT6B, and PKP3 were significantly upregulated in the psoriasis model compared to controls (**p < 0.01; ns, not significant, p>0.05) (n=12).

Figure 13 Validation of hub gene expression in clinical samples. RT-qPCR and transcriptome sequencing revealed that the expression levels of KRT6B, PKP3, and SPRR2B were higher than those in the control group. (PKP3,*p < 0.05; KRT6B, **p < 0.01; SPRR2B, ns, not significant, p>0.05) (n=21).

Bioinformatics Analysis of Clinical Samples

To investigate the expression patterns of SPRR2B, PKP3, and KRT6B in our clinical samples, we performed unsupervised hierarchical clustering based on the transcript levels obtained from RNA sequencing. The heatmap (Figure 14A) clearly shows that, with the exception of a few outliers due to individual variations, the two clusters largely corresponded to the control and psoriasis groups, indicating a strong association between the expression profiles of these three genes and disease status. Compared to the control group, all three genes were generally highly expressed in the psoriasis group.

Figure 14 Analysis of Hub Gene Expression Patterns in Clinical Samples. (A) Unsupervised hierarchical clustering analysis based on the expression levels of SPRR2B, PKP3, and KRT6B shows that their expression profiles distinguish the psoriasis group from the control group. (B) Volcano plot of differentially expressed genes, indicating that SPRR2B is significantly upregulated in the psoriasis group (log2FC = 6.08). (C) Box plots showing the expression distribution of SPRR2B, PKP3, and KRT6B in the two groups (KRT6B: *p < 0.05; SPRR2B: **p < 0.01; PKP3: ns, not significant, p>0.05. Wilcoxon test).

According to the volcano plot (Figure 14B), SPRR2B was significantly upregulated in the psoriasis group with marked statistical significance (log ₂FC = 6.08, p < 0.01). KRT6B was also upregulated in the psoriasis group, with an approximately two-fold increase in expression; however, this change was not statistically significant (log ₂FC = 0.97, p > 0.05). In contrast, PKP3 exhibited neither notable changes in expression nor statistical significance (log ₂FC = −0.11, p > 0.05). Given the limitations in clinical sample size and quality, the RNA-seq results were regarded as supplementary validation for preliminary RT-qPCR experiments, rather than conclusive evidence.

Finally, to further verify and visualize the differential expression of the target genes, we examined their expression distributions using box plots (Figure 14C). KRT6B expression was significantly higher in psoriasis samples than in controls (p < 0.05, Wilcoxon test), whereas PKP3 expression was slightly elevated but not statistically significant according to the Wilcoxon test, and SPRR2B was significantly upregulated in psoriasis (p < 0.01, Wilcoxon test).

In summary, the integrated analysis of the heatmap, volcano plot, and box plots supports the conclusion that KRT6B, PKP3, and SPRR2B were highly expressed in the collected psoriasis clinical samples, even when considering patient heterogeneity and the lack of statistical significance for some genes in certain analyses.

Discussion

As a chronic inflammatory skin disease, psoriasis is characterized by immune disorders, abnormal proliferation, and differentiation of epidermal KCs and is often associated with psoriatic arthritis.^2,30 Autophagy is the basic process of cell degradation and recovery. Therefore, we comprehensively explored the role of autophagy-related genes in psoriasis using a combination of scRNA and seq-RNA. Based on our previous bioinformatics analysis, we found that the signal balance of NK cells was disrupted, and seven candidate hub genes were screened. The expression of RT-qPCR in animal models and clinical samples was verified, and it was confirmed that KRT6B, PKP3, and SPRR2B are important genes that regulate autophagy in psoriasis microenvironment.

According to the study by Reemann et al, most of the hub genes investigated in the present study are highly expressed in various cell types within normal skin samples, including keratinocytes, fibroblasts, and melanocytes. Notably, PKP3 exhibits pronounced cell type specificity, with expression primarily enriched in keratinocytes.³¹ Consistent with this, our AUcell analysis revealed that cells with high autophagic activity were predominantly enriched within keratinocytes, indicating that keratinocytes serve as the predominant site of autophagy in psoriasis, and may drive disease progression. Specifically, in psoriatic keratinocytes, IL-17A activates mTORC1, thereby inhibiting autophagy and subsequently promoting the release of pro-inflammatory cytokines (IL-1β, IL-6, and TNF-α).^32,33 Furthermore, impaired autophagic flux disrupts nuclear degradation and lysosome-mediated protein clearance during KC differentiation, which compromises terminal differentiation and epidermal barrier function,^34,35 thereby exacerbating psoriatic skin pathology. Correspondingly, previous studies by Wang et al have proved that the incidence of KC autophagy is positively correlated with the severity of psoriasis in human patients and mouse models.³⁶ In previous experiments, the expression of keratin-related proteins decreased in mice lacking Atg7, suggesting that autophagy regulates the expression of KC proteins involved in barrier function.³⁷ At the same time, autophagy damage is related to inflammation and KC differentiation and destruction.³⁸ Combined with our research, autophagy is active in KCs and participates in KC differentiation and skin barrier function. Regulation of autophagy may be a promising strategy for the treatment of inflammatory skin diseases.

To gain a deeper insight into the state changes of KCs in the disease, we performed a pseudotime analysis of their development. The results revealed that the RAS and NOD-like receptor signaling pathways were downregulated. Psoriasis is a skin disease, and local RAS is over-activated due to the increased expression of angiotensin converting enzyme (ACE).³⁹ It is reported that the serum ACE level of patients with psoriasis is elevated, and the elevated ACE level in psoriasis tissue is more likely to have disease specificity.⁴⁰ However, the NOD2 receptor has been proved to be able to induce autophagy by recruiting autophagy-related proteins to eliminate pathogens,⁴¹ and the data of epidermal transcription disorder show that,⁴² the activation of inflammatory corpuscles in keratinocytes benefits from the recognition signal of NOD-like receptor signal transduction and its related functions. However, the pathogenesis of the RAS and NOD-like receptor signaling pathways in psoriasis requires further investigation.

Dysregulated intercellular communication is a key characteristic of psoriasis, in addition to intrinsic pathological alterations within cells. Signal transduction between HLA-E and CD94/NKG2A was enhanced in the communication between keratinocytes and NK cells, whereas that of KLRK1 was weakened. Many studies have shows that,⁴³ the inhibitory receptor CD94/NKG2A can inhibit antigen-mediated T cell effector functions such as cytokine release and cytotoxicity. As an activating receptor, KLRK1 directly induces cytotoxicity and cytokine production in NK cells.⁴⁴ This imbalance may exacerbate inflammation. The imbalance of inflammation (especially the release of cytokines such as IL-1β, IL-17A, and TNF-α) leads to the over-activation of the PI3K/AKT/mTOR pathway, comprising phosphoinositide 3-kinase or phosphatidylinositol-3 kinase (PI3K), protein kinase B (Akt), and mechanistic target of rapamycin (mTOR), which reduces the expression of the autophagy marker ATg5, inhibits KC differentiation and autophagy, and enhances psoriasis-like lesions.^45,46 In addition, the IL-36 cytokine family is also involved in this positive feedback loop. Studies have shown that IL-36 cytokines are significantly overexpressed in both lesional and non-lesional skin of patients with psoriasis, suggesting a prerequisite role in the development of the disease.⁴⁷ Moreover, the expression of IL-36 cytokines can be induced by TNF-α, IL-17A, and IL-22, and in turn, IL-36 cytokines themselves can promote the expression of pro-inflammatory cytokines such as TNF-α, IL-6, and IL-8 in keratinocytes, thereby forming a positive feedback loop.^48,49 Therefore, we speculate that autophagy defects in NK cells further destroy their signal balance, restrict IFN-γ production,⁵⁰ affect immune surveillance, and weaken the clearance ability of hyperproliferative KCs, thus maintaining the local inflammatory-proliferative cycle.

Notably, autophagy-related genes have been implicated in various immune-mediated inflammatory diseases. Douroudis et al were the first to report an association between polymorphisms in the ATG16L1 gene and psoriasis vulgaris, revealing that certain alleles, such as rs13005285, and their corresponding haplotypes exhibit significantly abnormal distribution in patients with psoriasis.¹³ Specifically, functional defects in the ATG16L1 protein may disrupt the autophagic process, thereby affecting signaling pathways that regulate cytokine production, leading to the toxic accumulation of damaged proteins and organelles, and ultimately resulting in cell death, tissue damage, and chronic inflammation.¹³ Furthermore, dysregulation of the autophagy-apoptosis crosstalk is not confined to skin pathologies. In a study on recurrent miscarriage (RM), Rull et al found significantly elevated mRNA expression levels of TRAIL and S100A8 in placental tissues from RM patients. As an apoptosis-inducing ligand, the upregulation of TRAIL may reflect an imbalance between apoptosis and autophagy at the maternal-fetal interface, where autophagy plays a critical role in maintaining placental cell homeostasis and resisting inflammatory injury.⁵¹ In addition to autophagy-related genes, genetic variations in cytokines themselves also influence susceptibility to psoriasis. The study by Kingo et al provides preliminary evidence that polymorphisms in the genes encoding IL-19, IL-20, and IL-24 may affect susceptibility to palmoplantar pustulosis (PPP) and plaque-type psoriasis, suggesting that genetic variants of IL-19, IL-20, and IL-24 contribute to the pathogenesis of psoriasis and its subtypes by modulating epidermal function and inflammatory responses.⁵²

Against this established pathological backdrop, in order to further verify the relationship between KRT6B, PKP3 and SPRR2B and psoriasis, we analyzed the clinical skin lesions by RT-qPCR and transcriptome sequencing, and found that they were all up-regulated in patients’ samples, which was consistent with the results in previous bioinformatics and animal models. These genes play critical roles in epidermal homeostasis. As a type II keratin, the up-regulation of KRT6B expression will lead to dyskeratosis and scale accumulation.⁵³ PKP3 is a scaffold protein containing armadillorepeat domain, which has been proved to promote the proliferation of non-transformed keratinocytes.⁵⁴ SPRR2B is a member of the family of small proline-rich proteins, and its cross-linking properties determine the biophysical properties of cornifiedenvelope(CE), such as rigidity and flexibility.^55,56

At present, there are no experiments proving that these three genes play a role in regulating autophagy in psoriasis. We speculate that they inhibit autophagy through the accumulation of reactive oxygen species (ROS), which leads to chronic diseases. Transcriptomics has proved that KRT6B, PKP3 and SPRR2B play an important role in the synergistic biological process of “keratinization” and “keratinization envelope formation”,⁵⁷ which is consistent with the enrichment analysis results in our study. Studies have shown that NAPDH oxidase can produce ROS in psoriasis fibroblasts,⁵⁸ and the accumulated ROS can promote the expression of KRT6 and trigger phosphorylation of PKP3.^59,60 Notably, the structure of SPRR2B protein contains proline site which can be recognized by mTORC1, and the activation of mTORC1 has been proved to significantly inhibit the autophagy process of macrophages.^61,62 It has been reported in the literature that “skin barrier destruction -ROS accumulation-autophagy block” is a positive feedback event that promotes the progress of psoriasis,⁶³ from which we speculate that KRT6B, PKP3 and SPRR2B may be involved, and this vicious circle positive feedback mechanism is closely related to the occurrence and development of psoriasis.

Building on these mechanistic findings, it is important to explore their clinical translational potential. Drugs that restore autophagy and target PKP3, SPRR2B, and KRT6B during clinical transformation may become new therapies for psoriasis. Based on the over-activation of PI3K/AKT/mTOR pathway involved in our cell communication analysis, mTOR inhibitors (such as everolimus and tacrolimus) have been reported to be effective in the treatment of intractable psoriasis.⁶⁴ In addition, studies have confirmed that fenofibrate can restore autophagy by targeting IL-17A and up-regulating LC3 expression.⁶⁵ However, due to the short-term up-regulation of KRT6B and SPRR2B in normal wound repair,^66,67 complete occlusion may affect epidermal regeneration, so future therapy development needs to rely on targeted delivery to achieve therapeutic specificity.

This study has some limitations. Although we have verified the upregulation of genes in both animal models and clinical samples, owing to the heterogeneity of patients and the difficulty of sample collection, whether these findings can be extended to a wider population requires further investigation. To this end, we aimed to build a larger prospective clinical cohort. Second, the current work mainly focused on verifying the expression of the target gene, and no functional experiments have been carried out. This follow-up study aimed to further clarify the specific regulatory mechanisms in the pathogenesis of psoriasis by constructing a genetically modified animal model.

Conclusion

By integrating multiple omics data and experimental verification, this study successfully identified KRT6B, PKP3, and SPRR2B as the key nodes of autophagy regulation in psoriasis. Further analysis showed that these genes may disrupt the signal balance of NK cells in the psoriasis microenvironment by mediating the vicious circle of “skin barrier destruction -ROS accumulation-autophagy block”, leading to immune monitoring dysfunction and joint disease progression. Despite the limitations of the sample size, KRT6B, PKP3, and SPRR2B, as reliable diagnostic biomarkers and potential therapeutic targets, revealed a new mechanism of autophagy in psoriasis and provided a new direction for the accurate diagnosis and treatment of psoriasis.