Weighted Gene Co-Expression Network Analysis of Oxymatrine in Psoriasis Treatment

Introduction

Psoriasis is a chronic inflammatory skin disease, affecting 2–4% of the world’s population.¹ Moreover, obesity, fatty liver, metabolic syndrome, cardiovascular disease,² and diabetes are often associated with moderate and severe psoriasis.³ Conventional long-term systemic treatment of patients with psoriasis is limited by poor tolerance and toxicity.⁴ In patients with psoriasis, T cells produce high levels of cytokines, such as interleukin (IL)-17, which stimulate keratinocyte proliferation.^5,6 Certain immunosuppressants (methotrexate and cyclosporine) and retinoic acid have been shown to control the disease and relieve symptoms quickly and effectively.⁷ While these drugs are primarily used in the treatment of patients with moderate and severe psoriasis, their adverse effects are often more serious than those of topical drugs.⁸ In contrast, biological agents greatly improve psoriasis treatment in moderate to severe psoriasis cases.⁹ However, prolonged treatment might lead to adverse reactions, such as infections, allergic reactions, autoimmune reactions, and malignant tumors.¹⁰ Hence, new treatment strategies are required for patients who cannot tolerate the current ones.

Oxymatrine is a natural alkaloid compound that primarily exists in Sophora flavescens.¹¹ It possesses a wide range of biological activities, including antioxidant, anti-inflammatory,¹² antibacterial,¹³ antiviral¹⁴ and antitumor¹⁵ activities, while effectively protecting the heart,¹⁶ liver,¹⁵ lungs,¹⁷ kidneys,¹⁸ and blood vessels.¹⁹ In our previous clinical trial, oxymatrine was found to effectively treat severe plaque psoriasis and inhibit recurrence.²⁰ In addition, oxymatrine regulates lipid metabolism and substantially improves psoriasis-induced comorbidities, such as fatty liver and metabolic syndrome.²¹

In our previous study, we utilized transcriptomic analysis to identify differentially expressed genes before and after oxymatrine treatment and identified potential targets for the treatment of psoriasis.²² Hence, in the current study, we aimed to characterize the target genes and regulated pathways of oxymatrine in the treatment of psoriasis. To this end, weighted gene co-expression network analysis (WGCNA) was used to construct a network of complex genes and pathways that may be affected by oxymatrine during the treatment period. WGCNA uses information of thousands of the most varied genes, or all genes, to identify gene clusters of interest (referred to as gene modules) and performs critical association analyses with traits.²³ Hence, it is used to identify hub genes that drive vital cellular signaling pathways from large-scale gene expression profiles.²⁴ Thus, it serves as a convenient and effective means to screen core genes that may be potential biomarkers for clinical prognosis and treatment.²⁵ In this study, we constructed a novel co-expression network in psoriatic lesions before and after oxymatrine treatment using WGCNA and psoriasis area and severity index (PASI) with full-length transcriptional sequencing to explore changes in skin lesions upon oxymatrine treatment. Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), and protein–protein interaction (PPI) analyses were used to explore potential biomarkers and possible action mechanisms of oxymatrine in the treatment of psoriasis. Therefore, the existing drugs are limited to the treatment of this chronic recurrent psoriasis. It is urgent to explore new drugs for the treatment of psoriasis. Our clinical observation found that oxymatrine can effectively treat moderate and severe psoriasis and inhibit recurrence. Therefore, it is very important to explore the mechanism of its treatment of psoriasis.

Methods

Subject Enrollment, Treatment, and Sampling

Five of the 16 oxymatrine-treated psoriasis patients with skin lesions were randomly selected for full-length transcriptome sequencing. Sixteen patients were given 0.6 g/100 mL oxymatrine (Tianqingfuxin; Zhengdatianqing Company Ltd., Jiangsu, China) intravenously once a day for 8 weeks. Further details can be found in our previous report.²² Previously, we analyzed the differentially expressed genes from our transcriptome data and concluded that oxymatrine-mediated treatment of psoriasis involves multiple genes and pathways. In this study, we combined the changes in clinical PASI scores of patients with psoriasis before and after oxymatrine treatment with the WGCNA. In total, 15 PASI scores were obtained; five each for patients before and after oxymatrine treatment and five for normal subjects whose PASI was calculated to be zero. The patients were admitted to the dermatology clinic of Ningxia Medical University General Hospital. The inclusion criteria, exclusion criteria, criteria for withdrawal from the trial, oxymatrine treatment, and acquisition of normal skin were according to a previous report.²² This research was approved by the Ethics Committee of the General Hospital of Ningxia Medical University (clinical trial registration number: CHICCTR-TRC-14004301, China). All patients signed an informed consent form, and all procedures involving human participants were conducted following the tenets of Helsinki Declaration.

WGCNA

The WGCNA and good Samples Genes functions used the Bioconductor software package in the R programming language to determine whether the gene data met the requirements for network analysis.²³ In this method, highly correlated genes are clustered to form a module and their overall expression level is represented by the module characteristic gene (MES), which is associated with the external phenotype. The main steps of network construction are as follows: (1) define the similarity matrix; (2) select the weight coefficient β to transform the similarity matrix into an adjacency matrix;²⁶ (3) transform the adjacency matrix into a topological overlap matrix (TOM); (4) perform hierarchical clustering based on TOM to obtain a hierarchical clustering tree; (5) apply the dynamic tree cutting method to identify modules from the hierarchical clustering tree; (6) calculate the module feature element (ME) of each module. The Pearson correlation coefficient of the MES between modules was calculated, and the 1-Pearson correlation coefficient was defined as the average distance between the MES of each module. The average linkage hierarchical clustering method based on the minimum size (genome) 30 was used to cluster the MES of all modules, and those with high similarity were merged to form a co-expression network. Finally, a dynamic tree was used to divide the modules of hierarchical clustering results and merge the modules with genes <30 and cutting height <0.25.²⁷ In addition, we carried out WGCNA analysis of the corresponding mRNAs of all genes and further analyzed their expression levels.

Identification of Clinically Significant Modules in WGCNA

The discovery of modules and genes with biological or clinical significance is the primary goal of many co-expression analyses. The relationship between modulus and clinical characteristics was evaluated based on the correlation between modular characteristic genes and clinical characteristics. When dealing with sample characteristics, a measure of statistical significance between module feature genes and traits, can be defined, for example, using correlation (r) or p values. Modules with higher trait significance are considered to be related to clinical traits.²³

Identification of Hub Genes in Key Modules

First, the module eigengenes of the first principal component were calculated to determine the expression level of genes in the module. Strong correlation module was determined according to the Pearson correlation coefficient between the module and psoriasis. Hub genes were then determined. Gene significance (GS) was defined as the correlation between gene expression and various traits. In addition, module member (MM) was defined as the association between gene expression and characteristic genes of each module. Subsequently, the correlation between GS and MM was tested to verify module-character associations.

Functional Enrichment Analysis

GO analysis was performed to reveal the differentially expressed genes that have unique biological significance.²⁸ The KEGG database was used to identify this important approach.²⁹ Cluster Profiler packages³⁰ in R were used to analyze and demonstrate GO annotations and KEGG approaches.²⁸

PPI Analysis

We searched the online STRING database (http://www.string-db.org/) to identify and predict interactions between genes or proteins. The connectivity of the nodes in the network was very high, and the connectivity of each node in the network diagram of the light-green module was 149. Therefore, according to the weight value of the gene in the high-green module, the first 100 related pairs were selected to construct the network graph, and the degree of the network graph was recalculated and presented using different color gradients with red indicating the strongest interaction. We then used the target hub genes to construct a PPI network using the STRING database (https://string-db.org/) to identify potential functional proteins and important genes.

Histological Analysis and Immunohistochemistry (IHC)

Samples were embedded in paraffin, cut to 5-mm-thick cross-sections and placed on glass slides. The slides were then stained with hematoxylin and eosin (H&E). The spinous layer (epidermal thickness) of the skin was measured at the midpoints of the left and right sides of the slice, using image analysis software (Image-Pro Plus 6.0; Media Cybernetics, Silver Spring, USA). The average values of the three sites were taken as the epidermis thickness of the samples to be measured. The slides were deparaffinized and rehydrated, endogenous peroxidase was inactivated with methanol containing 0.3% hydrogen peroxidase for 30 min at room temperature (25°C), and the antigen retrieval method was induced by heat. The sections were then incubated with rabbit anti-human antibody at a 1:350 dilution of CNFN (Abcam, Cambridge, UK) and S100A8 (Abcam) for 24 h at 25°C and were rinsed with phosphate-buffered saline (PBS) three times. The reaction enhancement solution was added, incubated at room temperature for 15 min, and rinsed with PBS for three times. The sections were then incubated for 30 min with secondary antibodies at room temperature and rinsed with PBS three times. Finally, DAB was used to develop color, re-dye, and seal the film. Images were observed and photographed under a microscope (Motic Incorporation Ltd., Canada). Image-Pro Plus 6.0 immunohistochemical image analysis system was used to analyze the collected images. Using Image-Pro Plus, the epidermis thickness was measured at the midpoints of the left and right sides of the slice, and the average values of the three sites were taken as the epidermis thickness of the samples to be measured.

Quantitative Real-Time Polymerase Chain Reaction (qRT-PCR)

A total of 15 RNA samples were subjected to qRT-PCR with β-actin (ACTB) as the internal reference gene, using SYBR Green QPCR Mix (DF Biotech., CHENGDU, China). Subsequently, 1000 μL of TRIzol was added to the samples, mixed thoroughly, and incubated at room temperature for 30 min. Next, 200 μL of chloroform was added to the samples, which were then vortexed and incubated at room temperature for 10 min. Gene expression was quantified using the 2^–ΔΔCt method. These sequences are listed in Table 1.

Table 1 List of Primers (5′ - 3′)

Western Blotting

Protein was extracted from skin lesions using 400 μL of radioimmune precipitation assay (RIPA) buffer (Epizyme, Shanghai, China). Protein concentrations were measured using the BCA Protein Assay Kit (KeyGen Biotech, Nanjing, China). Equal amounts of protein were separated using 10% sodium dodecyl sulfate (SDS) electrophoresis and then transferred onto polyvinylidene difluoride membranes. The membranes were labeled with the following primary antibodies: β-actin (1:6000 dilution; Abcam, Cambridge, UK), S100A8 rabbit mAb (1:1000 dilution; Affinity Biosciences, USA), and SPRR2A mouse mAb (1:2000 dilution; Novus Biologicals, USA). Image density analysis was performed using ImageJ software (NIH Image, Bethesda, MD, USA).

Statistical Analyses

The data were analyzed using GraphPad Prism (version 8.0; GraphPad Software, San Diego, CA, USA) and are presented as mean ± standard deviation. Student’s t-test was used for comparisons between pre- and post-treatment, and the threshold for statistical significance was set at p < 0.05.

Results

Clinical Therapeutic Effect of Oxymatrine

The PASI score of five patients decreased from the baseline level of 42.06 ± 8.45 to 7.10 ± 5.74 after treatment (p < 0.001, Figure 1A).²² Eight weeks of treatment caused the skin lesions to subside, with some pigmentation spots noted in the regression sites.

Figure 1 Images of clinical samples. (A) Images of psoriasis lesions before and after treatment. (B) H&E staining of tissues from patients and normal group. Scale bar = 50 µm. (C) Images of H&E-stained sections showing epidermal hyperplasia. (D) Cluster tree and heat map of samples according to the expression data of 15 samples. This is a cluster analysis of the 15 samples we have obtained. The average expression level of genes was analyzed using WGCNA. *vs normal group, **p < 0.01, #vs before treatment, ##p < 0.01. Abbreviations: H&E, hematoxylin and eosin; WGCNA, weighted gene co-expression network analysis.

H&E staining showed loss of the granular layer, hyperplasia of the epidermis, thickening of the spinous layer, incomplete keratosis, and infiltration of inflammatory cells in patients with psoriasis. Before and after oxymatrine treatment, the thickness of the spinous layer of the skin decreased (p < 0.01, Figure 1B). After treatment, the skin lesions markedly improved, and erythema scales were appreciably reduced. Compared to the normal group, all layers of the skin in the patients tended to be normal after treatment. In addition, we measured the thickness of epidermal hyperplasia and found that considerably reduced after treatment and recovered to a level that was not noticeably different from that of normal controls (Figure 1C).

Data Preprocessing

The 15 samples were divided into two clusters, with five samples in the first cluster and 10 samples in the second cluster (Figure 1D). The first group of five samples comprised all patients with psoriasis and the second cluster comprised five patients after oxymatrine treatment and five normal controls. This arrangement showed that the chosen skin tissue had good consistency. After oxymatrine treatment, patients did not show prominent differences compared to the normal control group.

Weighting Coefficient β Selection

The co-expression network has the characteristics of a scaling network, that is, P (k)∼k⁻¹, where k represents the connectivity of the nodes. Therefore, according to the scale-free network rule, weighting coefficient β must satisfy the condition of a negative correlation between log (k) and log [P (k)]. The larger the correlation coefficient, the more significant the characteristics of scale-free networks. According to the correlation coefficient between log (k) and log [P (k)], the best weighting coefficient β was selected between different weighting coefficients β. Considering the complex calculations and running time, it is reasonable to select β = 6 to construct a co-expression network when the correlation coefficient is >0.85 (Figure 2A). When β = 6, there is a negative linear correlation between log (k) and log [P (k)]. Based on the results presented in Figure 2A, we conclude that when β = 6, the co-expression network is scale-free. A co-expression network was constructed based on the optimal soft threshold. Clinically relevant gene expression modules were screened. Only by finding the optimal soft threshold could we build a co-expression network based on this value, divide the genes into different modules, and draw the gene clustering tree.

Figure 2 Determination of soft-thresholding power and network heatmap plot of psoriasis samples. (A) Chart showing the correlation coefficients of log(k) and log(p(k)) corresponding to different soft thresholds; the second chart shows the mean values of gene adjacency coefficients corresponding to different soft thresholds, reflecting the average connectivity level of the network. (B) Clustering dendrograms of all genes, with dissimilarity based on topological overlap, together with assigned module colors; 15 co-expression modules were constructed and displayed in different colors. Modules are assembled according to gene similarity. (C) Heatmap of topological overlap in the gene network. The gene dendrogram and module assignment are shown on the left and at the top. In the heatmap, each row and column correspond to a gene, the light color denotes low topological overlap, whereas a progressively darker red shade denotes a higher topological overlap. (D) Cluster analysis.

Weighted Gene Correlation Network Analysis

According to the WGCNA algorithm, the correlation and adjacency matrices of the sample gene expression profile were calculated and transformed into the TOM. Simultaneously, we obtained a genetic system cluster tree. No significant difference was observed in the interactions among the different modules, indicating a high degree of independence (Figure 2C). As shown in Figure 2B, the clustering tree contained 15 co-expression modules, represented by black, blue, brown, dark-green, dark-gray, dark-magenta, dark-olive green, dark-orange, dark-turquoise, green, light-green, light-yellow, pale-turquoise, red, and steel-blue. The color of the module to which the corresponding branch is assigned, is listed below. Merged dynamic is the color after merging similar modules. Grey indicates that the genes in the module do not belong to any module. The number of genes in each module was as follows: 218 (black), 626 (blue), 351 (brown), 127 (dark-green), 173 (dark-gray), 239 (dark-magenta), 40 (dark-olive green), 358 (dark-orange), 186 (dark-turquoise), 272 (green), 2136 (light-green), 286 (light-yellow), 50 (pale-turquoise), 369 (red), and 51 (steel-blue). Hierarchical clustering of the samples (Figure 2D) showed that genes clustered in the same module had similar expression and function according to certain clinical features. These findings indicate that oxymatrine treatment of psoriasis involved a complex network. In summary, we first calculated the co-expression correlation coefficient among genes according to the measured gene expression level. Based on the calculation results of Figure 2A, the genes were clustered, and the gene tree was generated. Different branches of the cluster tree represent different gene modules, and different colors represent different modules. Based on the weighted correlation coefficient of genes, genes were classified according to expression patterns, and genes with similar patterns were classified as a module. The module genes most related to clinical symptoms were selected for follow-up analysis.

Correlation Between Modules and Identifying Clinically Significant Modules

Identifying genes related to a clinical trait is essential to explore the underlying molecular mechanisms. In this study, the four clinically related factors were age, sex, PASI score, and duration of psoriasis. We primarily observed the clinical symptom genes related to psoriasis, and therefore, mainly focused on the modules related to PASI. The correlation coefficient between ME and Pearson’s correlation coefficient of clinical features was calculated, a significant module was obtained, and the co-expression module was further screened. As shown in Figure 3A, we identified five modules that were significantly related to PASI (ME light-green [correlation coefficient (r) = 0.9, p = 4e-06], ME black [r = −0.68, p = 0.005], ME dark-green [r = −0.66, p = 0.008], ME green [r = −0.57, p = 0.03], and ME dark-turquoise[r = −0.56, p = 0.03]). The correlations among the remaining modules were not significant. By generating a feature gene adjacency heat map, we found that the light-green module exhibited the strongest correlation with the other modules; hence, it was identified as the central module (Figure 3B).

Figure 3 Module-trait associations and analysis of three important modules. (A) Heatmap of adjacencies in the eigengenes network. We can see which modules are closely related to phenotypic PASI. Light-green module exhibited the strongest correlation with PASI. Red represents high adjacency (positive correlation), whereas blue represents low adjacency (negative correlation). (B) Each row corresponds to a module eigengene, column to a trait. Each cell contains the corresponding correlation and p values. The cells are color-coded by correlation. (C) Scatterplots of GS (y-axis) versus MM (x-axis) for three modules (light-green, dark-green, and black modules). (D) MES of three modules. The above line shows the expression levels of all genes in the module (y-axis) and the same sequence sample (x-axis). The following line shows the corresponding ME expression values for the module (y-axis) and the same sorted sample (x-axis). Abbreviations: PASI, Psoriasis Area and Severity Index; GS, gene significance; MM, module member; MES, module eigengenes.

Target Module Analysis

MM value of each gene was calculated to identify hub genes in the module. The correlation between GS and MM of each module gene was analyzed to assess whether the MM value was closely related to the PASI degree. The results showed that the correlation coefficients between the GS and MM value of PASI were the highest in the high-green module (r = 0.92, p < 1E-200), dark-green module (r = 0.71, p = 9.1E-21), and black module (r = 0.7, p =2E-33, Figure 3C). The high-green, dark-green, and black modules had the highest positive correlation with the PASI, which is consistent with the above conclusion (Figure 3D). ME expression was highly correlated with that of all genes in the module.

GO and KEGG Analyses of the Light-Green Module

In the GO analysis, the genes associated with PASI were significantly enriched for biological processes (BP), cellular components (CC), and molecular functions (MF). There were 374, 60, and 87 significant pathways in BP, CC, and MF, respectively (p < 0.05). The top four BP terms were keratinocyte differentiation (p = 1.44E-07), skin barrier establishment (p = 3.83E-07), intermediate filament cytoskeleton organization (p = 3.22E-06), and cell division (p = 1.12E-05). The top four CC terms were cornified envelope (p = 5.24E-10), desmosome (p = 2.65E-07), extracellular space (p = 1.50E-05), and cytosol (p = 7.63E-05). The top four MF terms were serine-type endopeptidase inhibitor activity (p = 1.85E-06), chemoattractant activity (p = 3.39E-05), chemokine activity (p = 0.000138126), and receptor of advanced glycation end-product (RAGE) binding (p = 0.000428888). Thus, oxymatrine treatment considerably improved the PASI in patients with psoriasis, and altered many related CC terms. In particular, the abnormal barrier function of the skin improved and was regulated after oxymatrine treatment. While the lesions and PASI score were improved, cellular structures, such as the cytoskeleton and keratinized capsules, were regulated. Hence, oxymatrine maintained the dynamic balance of cells (Figure 4A–C).

Figure 4 Enrichment analysis. (A) GO enrichment of biological processes. (B) GO enrichment of cellular components. (C) GO enrichment of molecular function. (D) KEGG pathway enrichment. Abbreviations: GO, gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

In the KEGG analysis, 37 pathways were significantly enriched (p < 0.05). The top four enriched pathways associated with improved PASI scores were the IL-17 signaling pathway (p = 5.65E-06), NOD-like receptor signaling pathway (p = 7.97E-06), cytokine-cytokine receptor interaction (p = 9.48E-06), and cell cycle (p = 1.47E-05, Figure 4D). With the improvement in skin lesions in patients with psoriasis upon oxymatrine treatment, the abnormal pathways related to inflammation and proliferation improved. Hence, these pathways represent probable targets of oxymatrine.

Identification of Hub Genes

The light-green module contained 2136 genes. The PPI results for the genes in this module revealed 11,175 interacting gene pairs. We selected the top 100 gene pairs with the most significant weight coefficients, in which 38 gene pairs were determined to exhibit strong interactions after oxymatrine treatment.

The top five genes (cornifelin (CNFN), S100A8, SPRR2A, SPRR2D, and SPRR2E) were identified as hub genes that contribute to the effective treatment of psoriasis with oxymatrine (Figure 5A). CNFN and S100A8 levels were verified using immunohistochemistry. In the normal epidermis, CNFN was primarily distributed in the upper granular layer and stratum corneum (Figure 5B1). In contrast, in patients with psoriasis, keratin was widely distributed in the spinous layer (Figure 5B2). After oxymatrine treatment, keratin expression markedly decreased (Figure 5B3). In patients with psoriasis, S100A8 was expressed in all layers of the epidermis with the cytoplasm and nuclei of cells in each layer appearing stained. S100A8 expression in the lesions of patients before treatment was higher than that in healthy control skin (p < 0.01, Figure 5C1). Following oxymatrine treatment, S100A8 abundance decreased in patients (p < 0.01, Figure 5C2) with no significant difference between the after-treatment and normal control groups (p > 0.05, Figure 5C3). Hence, oxymatrine treatment restored S100A8 expression in skin lesions to normal levels.

Figure 5 Hub gene screening and verification. (A) Protein–protein interaction network of light-green module genes. The edge represents interactions between genes. A degree is used to describe the importance of the protein nodes in the network. (B) Immunohistochemical staining for CNFN and IOD values of CNFN. (C) Immunohistochemical staining for S100A8 IOD values of S100A8. (D) qRT-PCR of CNFN, S100A8, SPRR2A, SPRR2D, and SPRR2E mRNA expression. (E) Western blotting of S100A8, SPRR2A. (F) Schematic diagram of the relationship between EDC and psoriasis and oxymatrine. *vs normal group, *p < 0.05, **p < 0.01, #vs before treatment, #p < 0.05, ##p < 0.01. Abbreviations: CNFN, cornifelin; IOD, integrated optical density; qRT-PCR, quantitative real-time polymerase chain reaction.

The mRNA expression levels of the screened hub genes (CNFN, S100A8, SPRR2A, SPRR2D, and SPRR2E) was analyzed using qRT-PCR. Although the CNFN results were not significant, a trend was evident; CNFN mRNA expression increased in patients with psoriasis and decreased after oxymatrine treatment. Moreover, the expression of S100A8 (ONT.1) in patients with psoriasis before treatment was higher than that in normal subjects (p < 0.05). Meanwhile, a significant decrease in its expression was observed following oxymatrine treatment (p < 0.05) with no statistical difference between the post-treatment and normal groups (p > 0.05), indicating that the expression level of S100A8 mRNA returned to the normal level after oxymatrine treatment. (Figure 5D). This trend was also observed for S100A8 (ONT.2), SPRR2A, SPRR2D, SPRR2E (ONT.1), and SPRR2E (ONT.2). Hence, the ameliorative profile of oxymatrine for psoriasis may be altered at the mRNA level, thus regulating the abundance of abnormal proteins.

Western blotting analysis further revealed that the abundance of small proline-rich proteins (SPRR)2A was increased in patients with psoriasis before treatment (vs control group, p < 0.01); however, it decreased after oxymatrine treatment (vs before treatment, p < 0.05). This trend was also observed for S100A8 abundance (Figure 5E).

We also screened modules highly related to the clinical PASI scores at the mRNA level. PPI analysis results of the red module was consistent with the results of the above core genes (Figures S1 and S2). A schematic diagram of the relationships among CNFN, S100A8, SPRR, psoriasis, and oxymatrine have been presented in Figure 5F.

Discussion

Oxymatrine successfully achieved certain clinical outcomes in the treatment of psoriasis.³¹ In our study, skin lesions of patients with psoriasis before and after treatment with oxymatrine were sequenced in the full-length transcriptional group, followed by WGCNA analysis in combination with the PASI score. The light-green module (top 1) contained 2136 genes, with the most significant number of genes. In addition, the black, dark-green, green, and dark-turquoise modules were related to PASI changes, indicating that oxymatrine treatment of psoriasis involves a complex network.

We then sought to analyze the possible genes and pathways regulated by oxymatrine in treating psoriasis according to the light-green module genes. The KEGG analysis primarily identified enrichment of IL-17, as well as the NOD-like receptor (NLR) signaling pathway, cytokine-cytokine receptor interaction, and cell cycle pathway. IL-17 is a key factor in psoriasis pathogenesis.^32,33 In particular, Bunte et al found that T cells produce high levels of IL-17, promoting a self-magnifying feedforward inflammatory response in keratinocytes.³⁴ T cells also produce large amounts of other types of cytokines that interact with each other, resulting in skin thickening and inflammatory cell infiltration.^35,36 Meanwhile, using RNA sequencing analysis, Tervaniemi et al identified NLR signal transduction in immune-related pathways of psoriasis, with NOD2 upregulated in psoriatic epidermis.³³ Indeed, the epidermal cell cycle is altered in patients with psoriasis compared with that in normal subjects.³⁷ However, our results indicate that the marked decrease in PASI after oxymatrine treatment depends on the improvement in pathways associated with IL-17, NLR, cell cycle, and cytokines, characterized by a reduction in the PASI score and recovery of clinical lesions.

The GO term analysis primarily revealed enrichment of keratinocyte differentiation, skin barrier establishment, cell division, chemoattractant activity, chemokine activity, and (RAGE) receptor binding, among others. Psoriasis is caused by abnormalities between infiltrating immune cells and activated keratinocytes.³⁸ In psoriasis, chemokine receptors are expressed in immune cells to produce a myriad of cytokines. Cytokines and chemokines play a key role in the production of psoriatic plaques.³⁹ In addition, RAGE contributes to the regulation of various pro-inflammatory genes and is the main factor promoting psoriasis, which is an immune-mediated disease. RAGE plays an important role in inflammatory response and participates in the regulation of various pro-inflammatory genes. It is a potential regulator in the pathogenesis of psoriasis.^40,41 Our findings suggest that, with a prominent improvement in skin lesions and PASI score, oxymatrine not only regulates skin barrier integrity, but also affects the complex regulatory network of chemokines and cytokines.

The epidermal differentiation complex (EDC) is a gene cassette on human chromosome 1q21 that contains many genes associated with epidermal differentiation.^42,43 The cross-linking of epidermal proteins eventually leads to the establishment of a cornified envelope (CE), which is a thick peripheral protein envelope that stabilizes each keratinocyte and is essential for epidermal barrier function.⁴⁴ In healthy individuals, there is a precise balance between the proliferation of basal cells and desquamation of the keratinized layer, whereas, in some skin diseases, the balance is disrupted. For instance, psoriatic keratinocytes proliferate abnormally, and CE formation is incomplete. Mutation of EDC plays an important role in the pathogenesis of psoriasis.⁴⁵ Overexpression of CNFN—a protein component of epidermal keratinocytes⁴⁶—alters the protein composition of CE.⁴² HUB genes of the light-green module (CNFN, S100A8, SPRR2A, SPRR2D, and SPRR2E) were screened using PPI analysis. S100A8, SPRR2A, SPRR2D, and SPRR2E belong to the EDC.⁴⁷ CNFN expression is considerably increased in psoriatic skin.⁴⁸ S100A8—also known as MRP8—binds to keratin filaments in keratinocytes, leading to skin inflammation and cell differentiation.⁴⁹ Meanwhile, SPRR2A, SPRR2D, and SPRR2E belong to the SPRR family of genes. SPRRs encode the precursor proteins of the keratinocyte envelope and are specifically expressed during keratinocyte terminal differentiation.⁵⁰ Indeed, the SPRR protein is overexpressed in psoriatic skin;⁵¹ its expression is closely related to epidermal differentiation, and its role in keratinocyte differentiation makes it a candidate gene for psoriasis.^50,52 Our results showed that before treatment, the expression of EDC-related genes and CNFN was upregulated compared to the normal levels; however, the expression decreased notably after oxymatrine treatment. Moreover, the skin lesions of patients with psoriasis noticeably improved upon treatment, indicating that oxymatrine regulated the abnormal expression of EDC and CE. This is also consistent with the GO and KEGG analysis results related to keratinocyte differentiation. Additionally, oxymatrine has been suggested to improve abnormal CE in patients with psoriasis. Hence, CNFN, S100A8, and SPRRs are closely associated with excessive proliferation of keratinocytes in patients with psoriasis. Indeed, oxymatrine considerably reduced the abnormal expression of CNFN, S100A8, and SPRRs. We hypothesize that oxymatrine may improve skin lesions in patients with psoriasis by affecting EDC on 1q21.

The protective barrier is formed by complex genetic networks comprising cell differentiation, enzyme activity, and cell connection.⁴⁷ The abnormal expression of EDC in psoriasis affects the balance between keratinocyte proliferation and differentiation, leading to barrier dysfunction. Oxymatrine regulates keratinocyte differentiation by altering EDC expression and coordinating various signal transduction pathways. Hence, oxymatrine not only improves the development and repair of epithelial tissue but also maintains the homeostasis of keratosis. According to the EDC treatment strategy, gene networks cooperate with each other to maintain the stability of epidermal differentiation after oxymatrine treatment.

This study has certain limitations. First, it is necessary to further clarify the molecular mechanisms underlying the effects of the hub genes on psoriasis treatment both in vivo and in vitro. Second, similar to other analysis methods, a larger sample will produce more reliable and refined results. Therefore, more clinical samples are needed to improve our findings and clarify the potential mechanism by which these hub genes affect the efficacy of oxymatrine in treating psoriasis.

Conclusions

Oxymatrine can significantly improve the PASI score by regulating multiple pathways and inflammatory factors. Oxymatrine affects the cell cycle, inhibits abnormal proliferation of keratinocytes, and can improve the skin lesions of patients with psoriasis and make erythema scales disappear. More specifically, the EDC locus may be a target of oxymatrine action. Thus, oxymatrine is a promising drug and provides a new treatment strategy for psoriasis.