SPRR1B is Related to the Immune Microenvironment and Can Be Used as a Biomarker for the Diagnosis of Psoriasis

Introduction

Psoriasis is a common inflammatory skin disease that is classified as an autoimmune disorder. It is characterized by abnormal differentiation of epidermal keratinocytes (KCs) and excessive proliferation. The pathogenesis of psoriasis is complex and involves various factors such as genetic susceptibility, infection, immune abnormalities, and psychoneurology.¹ The prolonged course of the disease and recurrent flare-ups often cause significant physical and mental distress for most patients. The atypical nature of early lesions often leads to diagnoses of moderate-to-severe psoriasis, highlighting the need for the development of new diagnostic biomarkers that can facilitate timely diagnosis and treatment.²

Recent advancements in high-throughput gene microarray analysis have enabled the investigation of various diseases at multiple levels, including somatic mutations, copy number variations, transcriptome-level genomic expression, and epigenetic variation.^3,4 This technology has provided a platform for exploring new target genes for disease control. In the case of psoriasis, bioinformatics has identified specific genes that contribute to its progression. For instance, ANGPTL4 expression was found to be significantly elevated in murine models of psoriasis-like skin induced by imiquimod (IMQ). Injection of exogenous recombinant ANGPTL4 protein exacerbates epidermal hyperproliferation and psoriasis-like inflammation.⁵ OAS1, USP18, and IFIT3 also exhibited high expression levels in typical psoriatic tissues, and inhibition of IFIT3 hindered the proliferation of keratinocytes and the secretion of CXCL1, CCL20, IL-1β, and IL-6.⁶ However, the diagnostic potential of many psoriasis genes remains poorly understood. Further research is required to explore the diagnostic utility of these genes.

This study aims to identify new diagnostic genes for psoriasis through the analysis of psoriasis-related datasets and conducting biological experiments using raw signal analysis. We analyzed four GEO datasets to compare gene expression between psoriasis and healthy samples, and identify differentially expressed genes (DEGs). We then used two machine learning algorithms, WGCNA analysis, and regression analysis, to screen for novel psoriasis biomarkers with diagnostic potential. To validate the results of bioinformatics analysis, we performed RT-PCR and immunohistochemistry on clinical samples. Additionally, we utilized the immune infiltration analysis algorithm CIBERSORT to explore the relationship among immune cells, key genes, and psoriasis. The findings of our study provide insights into the underlying immunomodulatory mechanisms of psoriasis and identify novel target genes for psoriasis diagnosis.

Materials and Methods

Data Processing

We conducted a search for relevant gene expression datasets in the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) public database using “psoriasis” as a keyword. The following were the inclusion requirements: a) Homo sapiens as the top organism; b) inclusion of both healthy human skin samples and skin lesion samples from psoriasis patients; c) selection of datasets with larger sample sizes. We ultimately downloaded four microarray datasets: GSE13355, GSE109248, and GSE153007 as the training group, and GSE14905 as the validation group.

The GSE13355 dataset was obtained from the GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array platform, comprising 58 skin lesion samples from psoriasis patients (PP) and 64 normal healthy human skin samples (NN). The GSE109248 dataset, derived from the GPL10558 Illumina HumanHT-12 V4.0 expression beadchip platform, included 17 psoriasis skin lesion samples. The GSE153007 dataset, sourced from the GPL6480 Agilent-014850 Whole Human Genome Microarray 4x44K G4112F (Probe Name version) platform, featured 14 psoriasis lesions and 5 healthy control skin samples. The training group contained a total of 89 psoriasis cases, 89 psoriatic lesions, and 83 healthy human skin samples. The validation group, GSE14905, was obtained from the GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array platform and included 28 psoriatic lesion samples and 21 healthy control skin samples. Figure 1 shows a schematic illustration of the study workflow.

Figure 1 The study design flowchart illustrates the analytical pipeline for the investigation of gene expression data. The analysis includes several bioinformatics methods, such as GSEA and GSVA, to identify significant gene sets. DEGs are also detected using statistical methods, including LASSO and SVM-RFE. Furthermore, PPI networks are constructed to explore functional relationships among the identified genes. GO, KEGG, and DO databases are used to annotate the enriched gene sets and reveal their biological functions. WGCNA is employed to identify coexpressed gene modules and their associations with clinical traits. Using the CIBERSORT method, immune cell infiltration is analyzed. The accuracy of the analysis is evaluated using the ROC curve.

Data Pre-Processing and Differentially Expressed Gene (DEG) Detection

R software (version 4.3.1) (https://www.r-project.org/) was used to analyze the data. The original GSE13355 dataset’s CEL file was read using R software’s affy package and converted into a probe matrix file. Based on the annotation information for each corresponding platform, probe matrix files in the GSE13355 and GSE153007 datasets were converted into the respective gene names using Perl software (https://www.perl.org/). The average of all corresponding probes was taken for cases where a gene name corresponded to multiple probes. The validation group GSE14905 directly downloaded matrix files with gene names.

The training group’s expression data were merged using the R package sva, followed by the application of the Combat function to eliminate batch effects. The training group was subjected to differential analysis using the R package limma, which identified 57 differentially expressed genes, including 10 down-regulated and 47 up-regulated genes, with |logFC| > 2 and an adjusted P-value 0.05. The pheatmap R package was used to create heat maps.

Functional Enrichment

The biological functions and pathways enriched by the 57 DEGs were identified through Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses, and Disease Ontology (DO) enrichment analysis to determine the major diseases associated with these genes. GO enrichment analyses encompassed biological processes (BP), cellular components (CC), and molecular functions (MF). An adjusted P value below 0.05 was considered statistically significant.

PPI (Protein-Protein Interaction) Network Construction

The Search Tool for Retrieval of Interacting Genes (STRING) online library was used to create a protein-protein interaction (PPI) network. The selected differential genes were imported into STRING, with Homo sapiens chosen as the organism and a minimum interaction score of 0.4. The generated PPI network files were exported and uploaded to Cytoscape 3.9.1 software for enhanced visualization and analysis. Using the MCODE plugin for Cytoscape, tightly coupled modules in the PPI network were subsequently retrieved. Additionally, the cytoHubba plugin was used to locate important genes in the PPI network, and the MCC method in cytoHubba was used to select the top 10 genes in rank score.

Screening and Validation of Diagnostic Biomarkers

Employing two machine learning techniques, least absolute shrinkage and selection operator (LASSO) regression and support vector machine-recursive feature elimination (SVM-RFE), facilitated feature selection and the identification of crucial diagnostic genes. LASSO, a regression analysis algorithm using regularization for improved accuracy, was executed with R’s “glmne” package. SVM-RFE, a supervised machine learning method, employed R’s “e1071” package to identify optimal variables by eliminating feature vectors generated by SVM. Both techniques are commonly applied in classification and regression tasks. The genes resulting from the two classification model algorithms, the genes in module 1 calculated by MCODE, and the key genes determined by the MCC algorithm in CytoHubba were intersected and considered potential psoriasis biomarkers for further investigation. The intersected genes were then tested for differences in the validation group GSE14905, with P < 0.05 indicating significant differences. Using receiver operating characteristic (ROC) curves, the sensitivity and specificity of the diagnostic prediction model were evaluated. The results were quantified by determining the area under the ROC curve (AUC).

Weighted Gene Co-Expression Network Analysis

A weighted gene co-expression network was constructed using the R package WGCNA. The comprehensive gene expression matrix comprised 15,654 genes from both psoriasis lesion samples and healthy control samples in the training group. Data were not subjected to log2 transformation, as genes with minor fluctuations across all samples were excluded. To ensure the reliability of co-expression network results, the hclust function was employed to remove sample outliers, and the soft threshold was set to 8. The weighted adjacency matrix was subsequently converted into a topological overlap matrix (TOM) to measure each gene’s network connectivity. The correlation coefficients between genes were calculated, and a clustering dendrogram was constructed. Branches containing module eigengenes (ME) of over 50 genes were assembled into co-expression modules with designated color labels. After calculating the correlations between the modules, the modules that had the most similarities were combined. Finally, Clinical phenotypes (normal and psoriasis) and module eigengene correlations were visualized in a heatmap.

Identification and Correlation of Disease Immune Infiltrating Cells

The immune cell composition of 22 types in the psoriasis and control groups within the training set was determined using the CIBERSORT algorithm in R software. Histograms were used to display the immune cell composition across all samples. The correlations between immune cell subtypes and between important genes and anticipated immune cell levels were calculated using the R package corrplot, with p < 0.05 indicating statistical significance. To display the results, the R tools ggplot2 and ggpubr were used.

Gene Set Variance Analysis (GSVA) and Gene Set Enrichment Analysis (GSEA)

Using two different techniques, we discovered dysregulated marker pathways in psoriasis to further improve the enrichment analysis. First, the h.all.v2022.1.Hs.symbols.gmt and c2.all.v2022.1.Hs.symbols.gmt gene sets were retrieved from the Molecular Signature Database. Next, the differential expression of SPRR1B in each pathway was detected using the R packages GSVA and limma, with |t| > 2 as the cutoff value. The R package clusterProfiler was utilized to perform gene set enrichment analysis (GSEA) in this study. The enrichment study of 22 immune cells used the normalized enrichment score (NES) derived from GSEA as the immune infiltration value.

Acquisition of Skin Tissues

Forty human skin tissues were collected, including 20 psoriasis patient skin lesions and 20 healthy control skin tissues (obtained from discarded surgical skin samples). No significant differences in age and gender were observed among all participants, and no significant differences in Psoriasis Area and Severity Index (PASI) scores (ranging from 1.2 to 36.2) were found among psoriasis patients. The research content and process adhered to international and national ethical requirements for biomedical research. They were reviewed by the Scientific Ethics Review Committee of the Second Affiliated Hospital of Harbin Medical University and were in accordance with medical ethics.

Immunohistochemistry

A rabbit anti-SPRR1B (C-term) polyclonal antibody (cat. no. abs111836; Absin) was used. A humidity chamber was used to incubate the slides with the main antibody for a whole night at room temperature. A modified labeled avidin-biotin reagent, followed by a phosphate-buffered saline wash, was used to detect the bound antibody for 20 minutes. A 0.1% diaminobenzidine solution served as the chromogen and was applied for 5 minutes. Subsequently, slides underwent counterstaining with Mayer hematoxylin for a period of 5–10 minutes. Staining intensity and staining positive rate scores were used to reflect the immunohistochemistry (IHC) score. The five staining intensity grades, from 0 (negative), 0.5 (0.5+), 1 (1+), 2 (2+), and 3 (3+), were used to categorize the staining intensity scores. Low to high staining (0 to 100%) staining positivity values were available. In order to determine the overall IHC staining score, the staining intensity score and the staining positive rate score (0–300%) were multiplied. The total IHC score for the cytoplasm and membrane was divided into two categories: low expression and high expression.

Cell Culture and Inflammation Model Construction

Human immortalized keratinocytes (HaCaT cells) were cultured in high-glucose DMED containing 10% fetal bovine serum (FBS), 100 u/mL penicillin, and 50 u/mL streptomycin. The cells were cultured in a 37°C, 5% CO2 incubator until they reached the exponential growth phase. Subsequently, HaCaT cells were seeded in a 6-well plate at a cell density of 1×10^6 cells/well and incubated overnight. Once the cells were firmly attached, the old culture medium was discarded and 2 mL of sterile PBS was used for washing twice. Next, the HaCaT cell inflammation model was constructed by adding fresh medium containing 30ng/mL IFN-γ to each well, which was followed by stimulation for 24 hours. After discarding the drug-containing medium, complete medium was added to continue the culturing process.

Transfection

To prepare for transfection, HaCaT cells in the logarithmic growth phase were adjusted to a concentration of 1×105 cells/mL and seeded in a 6-well plate. This allowed the cell density to reach 30–50% by the time of transfection, which was performed 24 hours later. The transfection process followed the instructions provided in the lipofectamine 2000 transfection kit (Invitrogen Company, USA). The cells were transfected with si-NC and si-SPRR1B separately. The cells transfected with si-NC were cultured conventionally to construct the inflammation model for the NC group and INF-γ group. The cells transfected with si-SPRR1B were used to construct the inflammation model for the INF-γ + si-SPRR1B group. After 24 hours of continued culture, the cells were collected for subsequent experimental testing.

mRNA Isolation and Quantitative Real-Time PCR

The M5 Universal Plus RNA Mini Kit from Mei5 Biotech was used to extract total RNA from tissues, and the TAKARA reverse transcription kit was used to create cDNA. A Bio-Rad CFX96 Real-Time PCR Detection System was used to perform quantitative real-time PCR. Standardized relative mRNA expression levels were determined using GAPDH values. Supplementary Table 1 contains the primers used in this test.

Results

Identification of Differential Genes

The detailed experimental flowchart of this study is presented in Figure 1. Three datasets in the training group were pooled and subjected to differential gene expression analysis, comprising a total of 89 psoriatic skin lesion samples and 83 healthy human skin samples. Following the screening, 57 DEGs were found, with a ratio of 1:4.7 between down-regulated and up-regulated genes, consisting of 10 down-regulated genes and 47 up-regulated genes. All DEGs were represented by a volcano plot. (Figure 2A), revealing a greater number of up-regulated genes with more noteworthy adjusted P values and logFC absolute values compared to down-regulated genes. Additionally, a heatmap was used to display all differential genes (Figure 2B).

Figure 2 Analysis of genes with variable expression. (A) In the training group difference analysis volcano plot, the genes are represented by dots. Black dots indicate genes that exhibit no significant difference, while genes with red and green dots show up- and down-regulated expression, respectively. The respective names of the genes are used to identify those that exhibit statistically significant changes between the groups; (B) To see how differently expressed genes were expressed in the training group, a heat map was created. Up-regulated genes were represented in red, while down-regulated genes were shown in blue. The data were filtered using the criteria |logFC| > 2 and adjusted P-value < 0.05.

Enrichment Analysis

These 57 differentially expressed genes (DEGs) were analyzed using GO, KEGG pathway enrichment analyses, and DO enrichment analyses. The enriched GO terms primarily included leukocyte chemotaxis, cell chemotaxis, granulocyte chemotaxis, granulocyte migration, keratinocyte differentiation, and keratinocyte proliferation in biological processes (BP). In cellular components (CC), the enriched terms mainly comprised secretory granule lumen, cytoplasmic vesicle lumen, vesicle lumen, cornified envelope, and specific granule lumen. Molecular functions (MF) primarily included receptor ligand activity, signaling receptor activator activity, cytokine activity, chemokine activity, and chemokine receptor binding. The main diseases analyzed by DO enrichment encompassed female reproductive organ cancer, lung disease, skin disease, integumentary system disease, and hepatitis (Figure 2B). The enriched KEGG pathways featured the IL-17 signaling pathway, viral protein interaction with cytokine and cytokine receptors, chemokine signaling pathway, cytokine-cytokine receptor interaction, and nucleotide metabolism (Figure 3C).

Figure 3 GO, KEGG, and DO were used to detect the biological functions, enrichment pathways, and main disease types of 57 differential genes. (A) Improved version: Bubble plots based on differential genes were created to display the top 10 significant GO analysis results in BP, CC, and MF. The vertical axis displays the GO term names, while the horizontal axis represents the proportion of genes. The size of the circle reflects the number of enriched genes in each GO term, whereas the color denotes the significance of enrichment. A redder circle implies greater enrichment of differential genes in the GO term. The filtering condition for selecting GO terms was an adjusted p-value of less than 0.05; (B) Improved version: A bubble plot was constructed to illustrate the enrichment analysis of differential genes using the DO. The vertical axis denotes illness type, and the horizontal axis represents gene proportion. The size of the circle indicates the number of enriched genes associated with each disease type, while the color denotes the significance of enrichment. A redder circle signifies a higher degree of differential gene enrichment; (C) Improved version: A bubble plot was used to visualize the results of the KEGG enrichment study on differential genes. The names of the pathways are shown on the vertical axis, and the number of enriched genes in each pathway is shown on the horizontal axis. The color indicates the degree of enrichment, with a darker shade suggesting a greater degree of enrichment of differential genes in the pathway.

PPI Network Construction, MCODE Cluster Module Analysis, cytoHubba Identification of Intersecting Genes

The STRING database was employed to generate PPI networks for further exploration of the interactions between these DEGs. The resulting data were processed using Cytoscape to create a PPI network graph consisting of 34 nodes and 117 edges (Figure 4A). The MCODE plugin was utilized to identify key gene cluster modules (Figure 4B, C). With filtering conditions set to degree cutoff = 2, node score cutoff = 0.2, k-core = 5, and max depth = 100, two cluster module identification networks were obtained. Module 1 had the highest score of 5.667, consisting of 7 nodes and 17 edges, featuring genes primarily including CNFN, PI3, LCE3D, SPRR2G, SPRR3, SPRR1B, and S100A7A. Module 2 comprised 5 nodes and 10 edges, with a score of 5.000. Subsequently, the cytoHubba plugin was used to identify the hub genes. The top 10 key genes ranked by the cytoHubba’s MCC algorithm were S100A7, SPRR1B, SPRR3, LCE3D, SPRR2G, CXCL10, PI3, CCL20, CNFN, and ISG15, which were regarded the most critical genes in the PPI network that may play a substantial role in the pathogenesis of psoriasis.

Figure 4 PPI network construction, MCODE cluster module analysis, and the identification of the key gene using cytoHubba. (A) This image displays a PPI network graph based on differential genes, consisting of 34 nodes and 117 edges. Each node corresponds to a protein, and each edge corresponds to an interaction. The color of the node reflects the Log2FC value, with red indicating gene up-regulation, blue indicating gene down-regulation, and the deeper the color indicating a larger Log2FC fold. The degree value is reflected in the size of the node; the larger the node, the greater the degree value. The thickness of the border line symbolizes the size of the interaction score, with a thicker border line representing a higher score. Additionally, this graph identifies two gene cluster modules using the (B and C). MCODE plugin; (B) This image displays gene cluster 1, which has a score of 5.667, 7 nodes, and 17 edges; (C) This image displays gene cluster 2, which has a score of 5.000, 5 nodes, and 10 edges; (D) This image displays the key genes identified by the MCC algorithm in the cytoHubba plugin.

Screening of Key Characteristic Genes and Validation

We used two additional machine learning algorithms to look for signature genes in order to more correctly identify potential psoriasis biomarkers. The LASSO regression algorithm identified 13 signature genes, including SERPINB4, S100A7A, SPRR2C, LCE3D, AKR1B10, GJB2, OASL, RRM2, TYMP, SPRR1B, C5orf46, CCL18, and THRSP (Figure 5A). Similarly, the SVM-RFE algorithm identified 22 signature genes from the DEGs, such as S100A9, SPRR1B, OASL, SERPINB4, SERPINB3, S100A12, TYMP, SPRR2C, S100A7A, AKR1B10, S100A7, KYNU, PI3, GJB2, CXCL1, CNFN, KRT77, LCE3D, KRT16, SPRR3, C5orf46, and LTF (Figure 5B). We intersected the genes obtained from the two algorithms, the genes in module 1 calculated by the MCODE plugin, and the genes obtained by CytoHubba’s MCC algorithm utilizing R’s Venn package to find two intersecting genes, LCE3D and SPRR1B(Figure 5C), which can serve as marker genes for psoriasis diagnosis. Based on the results of |logFC| and adjusted P values of DEGs and PPI protein interactions, we selected the SPRR1B gene for further study. We first validated the expression of the SPRR1B gene using the validation data GSE14905(Figure 5D), and then performed ROC curve analysis, yielding a P-value <0.05(Figure 5E, F), which indicated a notable disparity in the diagnostic gene between the validation groups. Additionally, the validation dataset confirmed the strong discriminatory power of the key gene with an AUC of 1 in the training and validation groups, as well as the high accuracy of SPRR1B as a diagnostic gene. To further evaluate the value of SPRR1B as a biomarker, we conducted qRT-PCR (n = 4) to detect its expression in psoriasis patients’ and healthy controls’ skin lesions. SPRR1B levels were also considerably greater in the skin. (Figure 5G). Collectively, these findings imply that SPRR1B could be used as a diagnostic biomarker for psoriasis.

Figure 5 Screening for psoriasis signature genes and validation. (A) The image represents the results obtained from the LASSO logistic regression algorithm used for screening psoriasis signature genes. The horizontal axis reflects the log(λ) value, while the vertical axis represents the cross-validation error. The gene corresponding to the point with the smallest vertical coordinate is the characteristic gene that has diagnostic effect; (B) The image shows the results obtained from the SVM-RFE algorithm used for screening characteristic genes. The horizontal axis indicates the number of genes, and the vertical axis is the cross-validation error. The number of feature genes is the value corresponding to the smallest point in the vertical coordinate; (C) The VENN plot depicts the intersection of important genes identified by MCODE, cytoHubba, LASSO, and SVM-RFE. The plot was generated using the venn package; (D) The image represents the validation of the intersected genes in the validation group. The horizontal axis shows the sample type, with blue representing the control group and red representing the psoriasis group. The vertical axis depicts the expression of the distinctive gene SPRR1B. A P value of less than 0.05 indicates that the diagnostic gene is different in the validation group; (E and F) The ROC curves depict the diagnostic test’s false positive and true positive rates. The horizontal axis shows the false positive rate expressed as 1 minus specificity, while the vertical axis represents the true positive rate expressed as sensitivity; (G) The graphic depicts the results of real-time PCR investigations of the expression levels of SPRR1B in psoriasis sufferers’ and healthy controls’ skin lesions. The ***p<0.0001 indicates a highly significant difference between the two groups.

WGCNA Module Screens for Gene Modules Associated with Psoriasis

Utilizing WGCNA, we assigned genes exhibiting analogous expression patterns to identical modules through mean linkage clustering. The mean linkage and Pearson correlation techniques were used to cluster the training group sample dendrograms (Figure 6A), while a soft threshold of 8 was chosen by calculating the scale-free model fit and average connectivity. (Figure 6B, C). In total, eight modules were identified (Figure 6D), containing varying numbers of genes. The correlation between module trait genes and clinical characteristics was calculated, revealing through the module-trait correlation heatmap (Figure 6E) that the blue-green module exhibited a robust association with psoriasis. Notably, the intersection genes SPRR1B and LCE3D, identified in this study, were present within this module (Figure 6F).

Figure 6 scWGCNA constructs co-expression modules. (A) The cluster dendrogram of samples is accompanied by a trait heatmap; (B) The connection between scale-free fit metrics and each soft threshold is shown; (C) The relationship between mean connectivity and each soft threshold is shown; (D) The clustering dendrogram displays the branches of highly connected genes, resulting in eight gene co-expression modules; (E) The relationships between module feature genes and clinical phenotypes (normal and psoriasis) are depicted in a heatmap; (F) A scatter plot showcases the relationship between MM and GS in blue-green modules, where MM represents the correlation of module signature genes with individual genes, and GS refers to the correlation of gene expression levels with psoriasis.

Analysis of Immune Infiltration

Utilizing CIBERSORT, We looked at how immune cells differed in psoriatic lesions versus healthy human skin. The immune cell histogram depicts the proportions of 22 immune cell types in each sample (Figure 7A), which highlights the three most abundant cell populations: resting mast cells, resting dendritic cells, and M2 macrophages. In the correlation heatmap between immune cells (Figure 7B), the strongest positive correlation coefficient was found between activated CD4 memory T cells and naïve CD4 T cells, while resting mast cells demonstrated the largest negative correlation coefficient with follicular helper T cells in absolute value. We assessed the significance of differences in skin immune cell infiltration between psoriasis patients and controls (Figure 7C). A P value <0.05 in the graph indicates that the corresponding immune cells differ between control and psoriasis groups. Resting mast cells in psoriasis skin lesions were much fewer than in healthy controls, although follicular helper T cells and M0 macrophages were substantially greater, suggesting that these cell types play crucial roles in the immune microenvironment of psoriatic lesions.

Figure 7 Immune cell infiltration in psoriasis. (A) The histogram of immune cells in the training group shows the distribution of immune cells across different sample types. The horizontal axis shows the sample type, whereas the vertical axis reflects the immune cell content. On the right side of the picture, the color of each immune cell is shown, providing a visual representation of the different cell types and their relative abundance in each sample; (B) The heat map displays the relationship between immune cells. The names of immune cells are represented by the horizontal and vertical axes, and the values inside each square show the correlation coefficients between immune cells. Positive correlations are shown as red squares, whereas negative correlations are shown as blue squares. The heat map provides a useful visual representation of the relationships between different immune cell types; (C) The violin plot depicts the difference in immune cell composition between the psoriasis and control groups. The horizontal axis indicates immune cell names, whereas the vertical axis reflects immune cell concentration. The blue violin represents the control group, while the red violin represents the psoriasis group. The violin plot depicts the variations in immune cell content between the two groups in a straightforward and succinct manner.

Correlation Analysis of Characteristic Genes and Immune Cells

The correlation analysis between the characteristic gene SPRR1B and each immune cell is presented in Figure 8. SPRR1B showed positive correlations with M0 macrophages, follicular helper T cells, M1 macrophages, activated dendritic cells, activated CD4 memory T cells, resting NK cells, monocytes, gamma delta T cells, neutrophils, activated mast cells, naive CD4 T cells, resting CD4 memory T cells, naive B cells, and memory B cells. On the other hand, it exhibited negative correlations with resting mast cells, activated NK cells, regulatory T cells (Tregs), plasma cells, resting dendritic cells, M2 macrophages, and eosinophils. The strongest positive correlation was observed with M0 macrophages and follicular helper T cells, while the strongest negative correlation was found with resting mast cells. For a more visually descriptive correlation plot, please refer to Supplementary Figure 1.

Figure 8 The figure presents a scatter plot depicting the relationship between a characteristic gene and immune cells. The horizontal axis reflects the level of expression of the distinctive gene, while the vertical axis represents the content of immune cells. The correlation coefficient (R), which measures the degree and direction of the association between the two variables, is shown. When R > 0, it indicates a positive association between the expression of the characteristic gene and the content of immune cells. Conversely, when R < 0, it denotes a negative correlation. Furthermore, the p value is also presented, which indicates the statistical significance of the correlation between immune cells and the characteristic gene. A p value < 0.05 implies that the observed correlation is unlikely to have occurred by chance alone and is, therefore, significant.

Enriched Hallmark Pathways and Curated Pathways

We identified the enriched hallmark and curated pathways in psoriasis using two distinct methods. GSEA results in Figure 9A demonstrated that the five most significantly enriched hallmark pathways in psoriasis were ALLOGRAFT REJECTION and E2F TARGETS. In Figures 9B-C, the GSVA-enriched hallmark and curated pathways included PROTEIN SECRETION, TGF BETA SIGNALING, NOTCH SIGNALING, FATTY ACID METABOLISM, MYOGENESIS, ANDROGEN RESPONSE, ADIPOGENESIS, BILE ACID METABOLISM, and UV RESPONSE. Additional enriched pathways encompassed those SILENCED BY METHYLATION, related to BREAST CANCER, RESPONSE TO NUTLIN 3A, AGING MIDDLE, PHOSPHORYLATION OF EMI1, FHIT TARGETS, MITOSIS LIN9 TARGETS, E2F4 TARGETS, FREE FATTY ACID RECEPTORS, E2F ENABLED INHIBITION OF PRE REPLICATION COMPLEX FORMATION, IL22 AND IL17A SIGNALING, OAS ANTIVIRAL RESPONSE, MED AND PSEUDOACHONDROPLASIA GENES, NATURAL KILLER CELLS, and INTERLEUKIN 36 PATHWAY. Pathways related to psoriasis include UV RESPONSE, SILENCED BY METHYLATION, IL22 AND IL17A SIGNALING, and INTERLEUKIN 36 PATHWAY.

Figure 9 GSEA and GSVA reveal enriched hallmark pathways in psoriasis. (A) GSEA pinpoints the top five pathways that are considerably enriched in psoriasis skin lesion samples; (B and C) The lollipop plot displays the enriched hallmark pathways and curated pathways in the psoriasis group, as determined by GSVA. GSEA is an abbreviation for gene set enrichment analysis, whereas GSVA is an abbreviation for gene set variation analysis.

Increased Expression of SPRR1B in Skin Lesions by Immunohistochemical Analysis

Furthermore, we used IHC to verify the expression of SPRR1B in skin lesions of psoriasis patients and healthy controls (Figure 10A). We found that SPRR1B was highly in skin lesions of psoriasis patients and confirmed that SPRR1B has the potential as a diagnostic biomarker for psoriasis.

Figure 10 Biological experiments verify the relationship between SPRR1B and psoriasis. (A) Immunohistochemical analysis of SPRR1B expression in psoriasis patients and healthy controls. (B) Expression of SPRR1B in Hacat cells treated with IFN-γ (C) Expression of SPRR1B after HaCat cells were transfected with si-SPRR1B (D–G) Expression of IL17, IL22, KRT6, KRT16 in HaCat cells. Data are expressed as mean ± SD. *P < 0.05; **P < 0.01; ***P < 0.001.

Effect of SPRR1B on Psoriasis-Related Markers

To investigate the impact of SPRR1B on psoriasis-related markers, HaCat cells were treated with IFN-γ to establish an inflammatory cell model. The findings revealed that the expression of SPRR1B increased with the increasing concentration of IFN-γand a significant difference in SPRR1B mRNA expression was observed when the concentration of IFN-γ reached 30ng/mL (P<0.01) (Figure 10B). A specific RNA interference sequence was designed for SPRR1B, resulting in a significant reduction in its expression (Figure 10C). Following treatment with 30ng/mL IFN-γ, the expressions of IL-17, IL-22, Keratins 6 (KRT6), and KRT16 were significantly elevated. However, upon silencing SPRR1B, the elevated levels of IL-17, IL-22, KRT6, and KRT16 were significantly reduced (Figure 10D–G).

Discussion

Psoriasis is a chronic and recurrent inflammatory condition that affects a significant number of individuals. Therefore, timely diagnosis plays a crucial role in managing this condition. However, clinical practice lacks sufficient early diagnostic tools. In this study, we conducted an analysis of four GEO datasets and successfully identified 57 differentially expressed genes (DEGs) between samples from psoriasis patients and healthy individuals. KEGG pathway analysis demonstrated their involvement in the IL-17 signaling pathway, viral protein interaction with cytokines and cytokine receptors, chemokine signaling pathway, cytokine-cytokine receptor interaction, and nucleotide metabolism. These findings suggest that these DEGs play an important role in the inflammatory response process and may influence the development of psoriasis.

To identify potential diagnostic biomarkers for psoriasis, a protein-protein interaction (PPI) network was constructed to study the relationship among 57 DEGs. Two machine learning algorithms were employed to identify two crossover genes, namely LCE3D and SPRR1B, which showed strong associations with psoriasis. Subsequent biological experiments conducted on skin tissue samples confirmed that the RNA and protein levels of SPRR1B were significantly increased in psoriasis lesions.

SPRR1B (small proline-rich protein 1B) is a member of the SPRR family, which encodes proteins that act as envelope components of keratinocytes. This family includes SPRR1 (SPRR1A, SPRR1B), SPRR2 (SPRR2A-SPR2F), SPRR3, and SPRR4. SPRR1 and SPRR2 are found in the skin, conjunctiva, and oral cavity, while SPRR3 is present in the oral cavity and esophageal epithelium. SPRR4 is located in the gut.⁷ Research has demonstrated that SPRR1B can stimulate cell proliferation.⁸ In studies on lung cancer, it was observed that SPRR1B expression reduces lung adenocarcinoma proliferation, induces G2/M phase arrest, and enhances apoptosis.⁹ Moreover, SPRR1B is highly induced during the differentiation of epidermal keratinocytes and is closely associated with KC differentiation.¹⁰ Excessive proliferation of keratinocytes (KCs) is considered a significant factor in the development of psoriasis. KCs, which are a vital component of the skin epidermis, serve as both a physical barrier and a source of local inflammation by attracting and activating immune cells. Numerous studies have demonstrated that KCs can initiate local inflammation by activating antimicrobial peptides and various pro-inflammatory cytokines (such as TNF/IL), thereby triggering the inflammatory response in psoriasis.^11–13 Anne Bowcock et al discovered that the constitutive activation of CARMA2 in KCs contributes to the pathogenesis of psoriasis, indicating that KCs not only induce inflammation but may also serve as a trigger for psoriasis.¹⁴ In addition, it is worth noting that KCs play a crucial role in the maintenance of psoriatic lesions and may also contribute to psoriasis relapse.¹⁵ Based on our analysis and the validation results obtained from skin lesions, it is reasonable to suggest that SPRR1B could potentially be used as a biomarker for understanding the development of psoriasis.

Recent studies have demonstrated the significant role of immune cell infiltration in the development and progression of psoriasis. In 1995, Krueger’s research group discovered a fusion protein called diphtheria toxin fragment (DAB389IL-2) combined with human interleukin 2 (IL-2), which effectively inhibits the growth of activated lymphocytes in psoriasis patients.¹⁶ This finding emphasizes the importance of T cells in the clinical manifestations of psoriasis. There are strong interactions between innate immune cells (such as dendritic cells, macrophages, and neutrophils), adaptive immune cells (B and T cells), and resident skin cells (such as KCs, melanocytes, and endothelial cells), which seem to amplify and enhance immune cell function, leading to chronic inflammation.¹⁷ Therefore, it is crucial from an immunological perspective to assess immune cell infiltration and identify the diverse components of infiltrating immune cells in order to understand the molecular-level causality of psoriasis and develop new immunotherapeutic targets.

In this study, we examined the differences in immune cell patterns between psoriatic lesions and healthy human skin. Our findings revealed significant differences in immune cells between psoriatic lesions and healthy skin. The top three cell types that exhibited the largest differences were resting Mast cells, resting Dendritic cells, and M2 Macrophages. We conducted a correlation analysis between immune cells and observed a strong positive correlation coefficient between activated CD4 memory T cells and naive CD4 T cells. Additionally, we found the largest absolute negative correlation coefficient between resting Mast cells and follicular helper T cells. Compared to healthy controls, psoriasis patients’ skin lesions had significantly lower levels of resting Mast cells, while healthy controls had significantly lower levels of resting T cells. On the other hand, psoriatic lesions had considerably higher levels of follicular helper T cells and M0 macrophages compared to healthy controls. These findings highlight the importance of these specific cell types in the immunological environment of psoriatic lesions.

In psoriasis, macrophages play a crucial role in lymphocyte activation, proliferation, and promoting the inflammatory process. The quantity of macrophages significantly increases in psoriatic lesions,¹⁸ acting as a primary source of TNF-α, a key mediator of chronic inflammation.¹⁹ Macrophages, derived from monocytes, function as phagocytic and antigen-presenting cells within tissues. Upon exposure to different cytokines, unactivated M0 macrophages polarize into M1 or M2 macrophages. The proportion of M0 macrophages corresponds to the severity of psoriasis, with moderate to severe cases having fewer M0 macrophages compared to mild psoriasis.²⁰ M1 macrophages contribute to the activation of the inflammatory response, while M2 macrophages aid in resolving inflammation.²¹ Studies on mice have shown that reducing macrophages alleviates psoriatic inflammation and normalizes Th1 cytokine levels, such as IL-1a, IL-6, IL-23, and TNF-α.^18,22 Our findings align with these observations, highlighting the significance of macrophages in the development and maintenance of psoriatic lesions.

Mast cells play a less studied and still controversial role in psoriasis. It has been suggested that they secrete chemokines CXCL1 and CXCL2, which attract neutrophils during the early stages of inflammation.²³ Additionally, mast cells produce cytokines like IL-1, IL-22, and IL-8, which activate keratin-forming cells and contribute to the chemotaxis of neutrophils to the epidermis, leading to the formation of Munro’s microabscesses.^24,25 In a study by Yongjun Zhang et al mast cell activation was confirmed in psoriatic lesions using immunofluorescence staining. Interestingly, the analysis showed that resting mast cells were almost non-existent in psoriatic skin tissue, which aligns with our findings. The reduced activation of mast cells after etanercept treatment suggests that mast cells play a crucial role in chronic inflammation progression and the maintenance of the skin cell ecological niche.²⁶ However, other studies have proposed that IL-33 produced by mast cells can actually suppress psoriatic inflammation by inhibiting Th1.²⁶ Consequently, further research is warranted to fully understand the involvement of mast cells in the development of psoriasis.

SPRR1B showed a positive correlation with Macrophages M0 and T cells follicular helper, which were significantly increased in psoriatic lesions. On the other hand, it showed a negative correlation with Mast cells resting, which were significantly decreased in psoriatic lesions. These findings align with our observations of cellular content in both psoriatic lesions and healthy skin tissue. SPRR1B, a well-known marker of squamous cell differentiation, has primarily been studied in the context of cell proliferation, differentiation, and apoptosis.^9,10 It has also been investigated as a potential marker for tumors like oral squamous cell carcinoma, lung adenocarcinoma, and gastric cancer.^27,28 This is the first time we propose the involvement of SPRR1B in the immune microenvironment of psoriasis development, suggesting its potential influence on psoriasis progression through the modulation of various immune cells.

Literature data indicates that various cytokines, including IL-17, IL-22, IL-23, TNF-α, and IFN-γ, are associated with the development and occurrence of psoriasis. Among these cytokines, IFN-γ plays a central regulatory role.^17,29 Recent studies have also identified keratin6/16/17 (KRT6/16/17) as a crucial early barrier signal in psoriasis. Upregulation of these keratins leads to changes in KC proliferation, cell adhesion, migration, and inflammatory characteristics. This excessive proliferation of KC and activation of innate immunity occur after disruption of the epidermal barrier, triggering autoimmune activation of T cells and driving psoriasis.³⁰ To investigate the causal relationship between SPRR1B and psoriasis further, we developed an in vitro psoriasis keratinocyte model by stimulating HaCat cells with INF-γ. We observed a significant increase in SPRR1B expression levels in the cell model, consistent with our in vivo findings. Additionally, silencing SPRR1B in the cell model resulted in a significant reduction in mRNA expression levels of inflammatory factors IL-17 and IL-22. Similarly, mRNA levels of KRT6/16 were also significantly reduced upon SPRR1B silencing, reaching levels similar to those of cells in the control group. These findings suggest that SPRR1B not only serves as a marker for psoriasis but may also play a role in influencing the progression of the disease.

Therefore, the highly expressed SPRR1B in psoriasis could potentially serve as an immune marker and even a potential therapeutic target for psoriasis. However, further experiments are required to investigate its molecular mechanisms and its role in various cellular biological functions. This study has some limitations, such as a relatively small sample size, which calls for large-scale clinical trials for validation. Additionally, the data presented here only establish a correlation between psoriasis and immune cells, not causation, emphasizing the need for more comprehensive exploration. Finally, it is necessary to develop diagnostic and treatment solutions for SPRR1B, such as developing SPRR1B rapid detection kits and preparing monoclonal antibodies.

Conclusion

In summary, our investigation has identified and validated a new immune-associated biomarker, SPRR1B, which shows elevated expression levels in patients with psoriasis. The findings of this study shed light on the potential molecular pathways underlying psoriasis and offer promising therapeutic targets for the development of immunotherapy techniques to treat this condition.