Identification of potential key molecules and signaling pathways for psoriasis based on weighted gene co-expression network analysis

Abstract

BACKGROUND

Psoriasis is a chronic inflammatory skin disease, the pathogenesis of which is more complicated and often requires long-term treatment. In particular, moderate to severe psoriasis usually requires systemic treatment. Psoriasis is also associated with many diseases, such as cardiometabolic diseases, malignant tumors, infections, and mood disorders. Psoriasis can appear at any age, and lead to a substantial burden for individuals and society. At present, psoriasis is still a treatable, but incurable, disease. Previous studies have found that microRNAs (miRNAs) play an important regulatory role in the progression of various diseases. Currently, miRNAs studies in psoriasis and dermatology are relatively new. Therefore, the identification of key miRNAs in psoriasis is helpful to elucidate the molecular mechanism of psoriasis.

AIM

To identify key molecular markers and signaling pathways to provide potential basis for the treatment and management of psoriasis.

METHODS

The miRNA and mRNA data were obtained from the Gene Expression Omnibus database. Then, differentially expressed mRNAs (DEmRNAs) and differentially expressed miRNAs (DEmiRNAs) were screened out by limma R package. Subsequently, DEmRNAs were analyzed for Gene Ontology and Kyoto Encyclopedia of Genes and Genomics functional enrichment. The “WGCNA” R package was used to analyze the co-expression network of all miRNAs. In addition, we constructed miRNA-mRNA regulatory networks based on identified hub miRNAs. Finally, in vitro validation was performed. All experimental procedures were approved by the ethics committee of Chinese PLA General Hospital (S2021-012-01).

RESULTS

A total of 639 DEmRNAs and 84 DEmiRNAs were identified. DEmRNAs screening criteria were adjusted P (adj. P) value < 0.01 and |logFoldChange| (|logFC|) > 1. DEmiRNAs screening criteria were adj. P value < 0.01 and |logFC| > 1.5. KEGG functional analysis demonstrated that DEmRNAs were significantly enriched in immune-related biological functions, for example, toll-like receptor signaling pathway, cytokine-cytokine receptor interaction, and chemokine signaling pathway. In weighted gene co-expression network analysis, turquoise module was the hub module. Moreover, 10 hub miRNAs were identified. Among these 10 hub miRNAs, only 8 hub miRNAs predicted the corresponding target mRNAs. 97 negatively regulated miRNA-mRNA pairs were involved in the miRNA-mRNA regulatory network, for example, hsa-miR-21-5p-claudin 8 (CLDN8), hsa-miR-30a-3p-interleukin-1B (IL-1B), and hsa-miR-181a-5p/hsa-miR-30c-2-3p-C-X-C motif chemokine ligand 9 (CXCL9). Real-time polymerase chain reaction results showed that IL-1B and CXCL9 were up-regulated and CLDN8 was down-regulated in psoriasis with statistically significant differences.

CONCLUSION

The identification of potential key molecular markers and signaling pathways provides potential research directions for further understanding the molecular mechanisms of psoriasis. This may also provide new research ideas for the prevention and treatment of psoriasis in the future.

Key Words: Psoriasis; MicroRNAs; Weighted gene co-expression network analysis; Functional enrichment; MicroRNA-mRNA regulatory network

Core Tip: Psoriasis is a common chronic, recurrent, immune-regulatory skin and joint disease. Although psoriasis is widespread and has significant negative impact on patients’ life quality, it has not yet been fully diagnosed and treated. Moreover, it is also associated with many other diseases. So far, psoriasis is still a treatable, but incurable, disease. We use weighted gene co-expression network analysis to identify key modules and microRNAs (miRNAs) related to psoriasis and explore potential key pathways related to psoriasis through the targeting relationship of miRNA-mRNA. This provides new research ideas for the prevention and treatment of psoriasis in the future.

INTRODUCTION

Psoriasis is a common chronic, recurrent, immune-regulatory skin and joint disease. It has a variety of clinical skin manifestations, but the most common clinical manifestations are erythematous, scaling papules, and plaques[1]. Previous studies have found that psoriasis is affected by family history and age[2]. Mechanical stress, air pollutants, sun exposure, infection, lifestyle, obesity, dyslipidemia, and mental stress are associated with the progression of psoriasis[3]. The pathogenesis of psoriasis is more complicated and often requires long-term treatment. Mild to moderate psoriasis can be treated topically with a combination of glucocorticoids, vitamin D analogues, and phototherapy. Moderate to severe psoriasis usually requires systemic treatment[4]. Psoriasis is also associated with many diseases, for example cardiometabolic diseases, malignant tumors, infections and mood disorders[5]. So far, psoriasis is still a treatable, but incurable, disease. Therefore, identification of key genes and signaling pathways are of great significance for understanding the molecular mechanism of psoriasis and provides potential basis for the treatment and management of psoriasis.

MicroRNAs (miRNAs) are a class of small non-coding single-stranded RNAs that regulate gene expression[6]. MiRNAs play a vital role in the progression of diseases, for example cancer, infectious diseases, and immune diseases[7]. MiRNAs can regulate cell proliferation, keratinocyte differentiation, apoptosis, and atypical immune activation in psoriasis[8]. Previous studies have found that miR-187 can inhibit the proliferation of keratinocytes by targeting B7 homolog 3 protein. In addition, in the mouse model of psoriasis, overexpression of miR-187 can reduce acanthosis and reduce the severity of the disease[9]. MiR-183-3p can inhibit the proliferation and migration of keratinocytes in psoriasis by inhibiting growth factor receptor binding 2-associated binding protein 1[10]. MiR-320b negatively regulates keratinocyte proliferation in psoriasis by targeting AKT serine/threonine kinase 3[11]. The low expression of miR-31 in dermal mesenchymal stem cells (DMSCs) in patients with psoriasis results in the increased expression of some of its target genes, which in turn promotes the activation of T lymphocyte by inhibiting the proliferation of DMSCs, thereby contributing to the pathogenesis of psoriasis[12]. Previous studies have found that miR-146 regulates inflammatory responses in keratinocytes and skin fibroblasts, which may affect the pathogenesis of psoriasis[13]. Furthermore, in human keratinocytes, the expression of miR-146a is induced by pro-inflammatory cytokines tumor necrosis factor-α (TNF-α), interleukin-1B (IL-1B), and IL-17a, and the expression of miR-146b is induced by interferon-γ (IFN-γ) and IL-22. MiR-203 promotes keratinocyte proliferation by targeting liver X receptor-α and peroxisome proliferator-activated receptor-γ in psoriasis[14]. MiR-125 regulates keratinocytes proliferation by regulating targeted genes[15,16]. MiR-197 over expression inhibits keratinocytes proliferation induced by IL-22 and keratinocytes migration[17]. In addition, miR-99 also plays an important role in regulating the abnormal proliferation and differentiation of keratinocytes in psoriasis[18]. These studies demonstrate that miRNAs play a vital role in the pathogenesis of psoriasis. Currently, miRNAs studies in psoriasis and dermatology are relatively new. Therefore, the identification of key miRNAs in psoriasis is helpful to elucidate the molecular mechanism of psoriasis.

Weighted gene co-expression network analysis (WGCNA) can be used for finding clusters (modules) of highly correlated genes[19]. WGCNA has been used to identify key genes for many diseases, including cancer[20], cardiovascular disease[21], and immunological diseases[22]. In addition, WGCNA has also been used to identify key mRNAs and long noncoding RNAs in psoriasis. So far, we have not found studies that use WGCNA to identify key miRNAs in psoriasis. Therefore, we use WGCNA to identify key modules and miRNAs related to psoriasis and explore potential key pathways related to psoriasis through the targeting relationship of miRNA-mRNA.

MATERIALS AND METHODS

Data sources and processing

The miRNA and mRNA data were obtained from the Gene Expression Omnibus (GEO) database[23]. The keyword we searched in the GEO database was “psoriasis”. Cell line or animal level studies and single-sample studies were excluded. Five datasets GSE13355, GSE66511, GSE145054, GSE142582, and GSE129373 were involved in this study (Table 1). All data were derived from skin samples. GSE13355 and GSE66511 were mRNA data from GPL570 and GPL16288 platform, respectively. GSE145054, GSE142582, and GSE129373 were miRNA data from GPL19117, GPL20301, and GPL11154, respectively. The gene expression matrix files in GSE13355 and GSE145054 were downloaded. The GPL annotation file was used to annotate the expression matrix. The probe ID was converted to the gene symbol. Multiple probes correspond to the same gene, taking the average value. The probe that corresponds to multiple genes was removed. The original data of gene expression profile in GSE66511, GSE142582, and GSE129373 were downloaded, and performed by logarithm processing. Then, the combat function in “sva” R package was used to remove batch effects.

Table 1 Details of five datasets.

GEO ID	Samples, Normal:Case	Platforms	Author	Type	Source
GSE13355	64:58	GPL570	Gudjonsson JE	mRNA	Skin
GSE66511	12:12	GPL12688	Keermann M	mRNA	Skin
GSE145054	4:4	GPL19117	Mildner M	miRNA	Skin
GSE142582	5:5	GPL20301	Yu Z	miRNA	Skin
GSE129373	9:9	GPL11154	Srivastava A	miRNA	Skin

Differential expression analysis

After the above pretreatment, differentially expressed mRNAs (DEmRNAs) and differentially expressed miRNAs (DEmiRNAs) were screened out by limma R package[24]. DEmRNAs screening criteria were adjusted P (adj. P) value < 0.01 and |logFoldChange| (|logFC|) > 1. DEmiRNAs screening criteria were adj. P value < 0.01 and |logFC| > 1.5.

Functional enrichment analysis

To investigate the possible involvement of DEmRNAs in the biological processes related to psoriasis, we performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomics (KEGG) functional enrichment analyses. David database was used for the enrichment analysis of DEmRNAs[25], and the screening criterion was P < 0.05. GO enrichment analysis, which is comprised of biological processes (BP), cellular components (CC), and molecular functions (MF). The “GOplot” R package was used to visualize the enrichment results.

Construction of weighted gene co-expression network

The “WGCNA” R package was used to analyze the co-expression network of all miRNAs. To detect outliers, the “hclust” function was used to cluster the sample data. Subsequently, to construct a scale-free topology, the “pickSoftThreshold” function was used to select an appropriate soft threshold power regulator (height 0.90, β value 2). Calculate the adjacency matrix based on the kernel value. The adjacency matrix was transformed into topological overlap matrix (TOM) and corresponding dissimilarity matrix (1-TOM). Genes with similar expression patterns were gathered together, and modules were divided according to the function of “cutreeDynamic” with default parameters. Since the modules identified by the dynamic tree cutting algorithm may be similar, they are combined with a truncation of 0.25 height[26].

Hub modules and miRNAs

To identify the important modules associated with psoriasis, the module eigengene (ME) of each module was calculated using the “moduleEigengenes” function. Then, the correlation between ME and psoriasis was analyzed using Pearson method. Subsequently, the module with the highest correlation with psoriasis was identified as the hub module. According to module connectivity (MM) and clinical trait relationship (GS) of each gene in the hub module, the candidate hub miRNAs were screened out[27]. The screening criteria were MM > 0.8 and GS > 0.5. Finally, the intersection of candidate hub miRNAs and DEmiRNAs were selected as real hub miRNAs.

MiRNA-mRNA network

To explore the correlation between miRNA and mRNA, we constructed a miRNA-mRNA regulatory network. Target mRNAs of miRNAs were predicted using miRDB (http://mirdb.org) database. MiRNA-mRNA pairs involved in common DEmRNA and negatively regulated were selected to construct the network. Cytoscape (www.cytoscape.org/) was used to visualize the miRNA-mRNA regulatory network.

In vitro validation

We collected skin tissue samples from healthy control individuals and patients with psoriasis. Basic information (including age, sex, etc) of patients and healthy controls were recorded in detail during sample collection (Table 2). The specific inclusion criteria of patients with psoriasis: (1) Patients had at least one well-circumscribed erythematous and squamous lesion, which had been confirmed by at least two dermatologists; (2) Patient’s pathological tissues were confirmed by clinical histopathology; (3) Patients had no systemic anti-psoriatic treatment 2 wk before the skin biopsy; and (4) Patients had no topical anti-psoriatic treatment 1 wk before the skin biopsy. Patients with psoriasis who were treated with immune-related drugs before sampling and whose clinical information was incomplete were excluded. The individuals in the normal control group were sex and age matched with the psoriasis group, and without history of psoriasis or other autoimmune diseases.

Table 2 Clinical information of healthy control individuals and patients with psoriasis in vitro validation.

Group	Sample number	Sampling location	Sex	Age	Height, cm	Body weight, kg	BMI	Smoking	Drinking	History of ASCVD	History of hypertension	History of hyperlipidemia	History of diabetes	Psoriatic arthritis	Severity index, PASI	Whether to collect non-lesion parts
Normal control group	C1	Right forearm extension side	Female	47	158	54	21.6	No	No	No	No	No	No
	C2	Right sole of the foot	Male	60	178	70	22.1	Yes	Yes	No	No	No	No
	C3	Left sole of the foot	Male	59	183	81	24.2	No	No	No	No	No	No
	C4	Left armpit	Female	27	160	52	20.3	No	Yes	No	No	No	No
	C7	Left forearm extension side	Female	18	160	50	19.5	No	No	No	No	No	No
	C9	Right shoulder	Female	27	161	62	23.9	No	No	No	No	No	No
	C11	Right calf	Female	60	156	55	22.6	No	No	No	No	No	No
Psoriasis group	E2-2	Right forearm	Male	19	177	58	18.5	No	No	No	No	No	No	No	10.5	Yes
	E3-2	Back	Female	58	157	58	23.5	No	No	No	No	Yes	No	No	8	Yes
	E5-2	Left calf stretched side	Male	55	176	68	22	No	No	No	Yes	No	Yes	No	13.4	Yes
	E6-2	Outer left thigh	Female	44	166	67	24.3	No	No	No	No	No	No	No	11.9	Yes
	E7-2	Right thigh internal test	Male	62	170	66	22.8	Yes	No	No	No	No	No	No	8.6	Yes
	E8-2	Right forearm extension side	Male	48	170	58	20.1	Yes	No	No	No	No	No	No	15.4	Yes
	E10-2	Back	Male	18	172	60	20.3	Yes	No	No	No	No	No	No	22.2	Yes

ASCVD: Arteriosclerotic cardiovascular disease; BMI: Body mass index.

According to screening criteria, 7 normal and 7 patient skin tissue samples were obtained. Firstly, total RNA of the samples was extracted using TRIzol. Then, fastKing cDNA first strand synthesis kit and miRNA first strand cDNA synthesis kit (Tailing Reaction) were used for mRNA and miRNA reverse transcription, respectively. Subsequently, SuperReal PreMix Plus (SYBR Green) and miRNA quantitative PCR kit (SYBR Green Method) were used to perform real-time polymerase chain reaction (RT-PCR) validation of mRNA and miRNA, respectively. GAPDH, ACTB, and hsa-U6 were internal reference. ABI7300 fluorescence quantitative PCR instrument was used for detection. Finally, the 2^-ΔΔCt method was used for relative quantitative analysis of the data[28]. This study complied with the Declaration of Helsinki. Informed consent was obtained from the participants. All experimental procedures were approved by the ethics committee of Chinese PLA General Hospital, No. S2021-012-01.

Statistical analysis

All statistical analyses were performed using R software. Limma R package was used to screen for DEmRNAs and DEmiRNAs. Pearson’s method was used to analyze the correlation between ME and psoriasis. RT-PCR validation data were statistically analyzed by one-way ANOVA. P < 0.05 was statistically significant.

RESULTS

Analysis of DEmRNAs and DEmiRNAs

After pretreatment of the raw data, 15780 mRNAs and 837 miRNAs were screened out (Figure 1). Compared with the control group, 639 DEmRNAs were identified in the psoriasis group (adj. P value < 0.01 and |logFC| > 1). Among them, 497 were up-regulated and 142 were down-regulated. Compared with the control group, 84 DEmiRNAs were identified in the psoriasis group (adj. P value < 0.01 and |logFC| > 1.5). Among them, 70 were up-regulated and 14 were down-regulated. Heat maps showed that there were significant differences of mRNA (Figure 2A) and miRNA (Figure 2B) expression between psoriasis group and control group. The volcano maps of DEmRNAs and DEmiRNAs were shown in Figure 2C and D.

Figure 1 Removal of data batch effect. A: Ven diagram of batch effect of mRNA; B: Boxplot of batch effect of mRNA; C: Ven diagram of batch effect of miRNA; D: Boxplot of batch effect of miRNA.

Figure 2 Identification of differentially expressed mRNAs and differentially expressed miRNAs. A: Heat map of top 100 differentially expressed mRNAs (DEmRNAs). Complete‑linkage method combined with Euclidean distance was used to construct clustering; B: Heat map of differentially expressed miRNAs (DEmiRNAs); C: Volcano map of DEmRNAs. Blue, red, and gray points represent down-expressed, up-expressed, and not DEmRNAs, respectively; D: Volcano map of DEmiRNAs.

Analysis of biological functional

Functional analysis of 639 DEmRNAs was performed using David database (P < 0.05). GO enrichment results showed that most of the DEmRNAs were distributed in inflammatory response (GO:BP), immune response (GO:BP), cytosol (GO:CC), extracellular exosome (GO:CC), protein homodimerization activity (GO:MF), and other biological functions (Figure 3A-C). KEGG enrichment results demonstrated that DEmRNAs were significantly enriched in metabolic pathways, influenza A, chemokine signaling pathway, and cytokine-cytokine receptor interaction (Figure 3D and Table 3).

Figure 3 Gene Ontology and Kyoto Encyclopedia of Genes and Genomics analysis of differentially expressed mRNAs. A: Biological processes terms were enriched of differentially expressed mRNAs (DEmRNAs) by Gene Ontology (GO) function analysis; B: Cellular components terms were enriched of DEmRNAs by GO function analysis; C: Molecular functions terms were enriched of DEmRNAs by GO function analysis; D: Kyoto Encyclopedia of Genes and Genomics enrichment of DEmRNAs.

Table 3 Kyoto Encyclopedia of Genes and Genomics function enrichment of differentially expressed mRNAs.

Category	hsa	Term	Count	P value	Genes
KEGG_PATHWAY	hsa03320	PPAR signaling pathway	14	3.07E-06	GK, MMP1, ADIPOQ, LPL, FADS2, ACADL, ACOX2, FABP5, FABP7, ACSBG1, PLIN1, ANGPTL4, SLC27A4, PPARD
KEGG_PATHWAY	hsa05164	Influenza A	22	1.05E-05	NLRX1, RSAD2, DDX58, STAT1, TMPRSS4, MX1, MAPK13, PYCARD, CXCL10, SOCS3, OAS1, OAS2, IL-1B, OAS3, TNFSF10, IRF7, CCL2, PRSS3, KPNA2, PRSS2, MYD88, IRF9
KEGG_PATHWAY	hsa05133	Pertussis	12	2.56E-04	PYCARD, C1QB, CASP7, NOS2, CALML5, IL-1B, IRF1, IRF8, CALML3, FOS, MYD88, MAPK13
KEGG_PATHWAY	hsa05162	Measles	16	4.14E-04	DDX58, STAT1, MX1, CD3G, CD3D, OAS1, CCNE2, CCNE1, OAS2, IL-1B, OAS3, TNFSF10, IRF7, PRKCQ, MYD88, IRF9
KEGG_PATHWAY	hsa04668	TNF signaling pathway	14	4.87E-04	MLKL, CCL20, CXCL1, FOS, NOD2, SELE, MMP9, MAPK13, CXCL10, SOCS3, CASP7, IL-1B, CCL2, JUNB
KEGG_PATHWAY	hsa04062	Chemokine signaling pathway	19	7.44E-04	LYN, CXCL9, CCL22, STAT1, CCL20, CXCR4, CXCL1, CXCL13, CXCL10, HCK, CCL8, PLCB4, CXCR2, RAC2, CCL2, CCR7, CCL19, CCL18, CCL27
KEGG_PATHWAY	hsa04060	Cytokine-cytokine receptor interaction	21	0.002838	CXCL9, IL20, IL4R, CCL22, CCL20, CXCR4, IL19, CXCL1, CXCL13, CXCL10, CCL8, IL-1B, CXCR2, TNFSF10, CCL2, CCR7, CCL19, IL7R, CCL18, CCL27, TNFRSF21
KEGG_PATHWAY	hsa05160	Hepatitis C	14	0.003635	DDX58, STAT1, IFIT1, CLDN1, MAPK13, SOCS3, OAS1, OAS2, IRF1, OAS3, CLDN8, IRF7, CLDN17, IRF9
KEGG_PATHWAY	hsa00590	Arachidonic acid metabolism	9	0.003655	PLA2G2F, HPGDS, GPX2, PLA2G4D, PLA2G2A, PLA2G3, ALOX12B, CYP2E1, CYP4F8
KEGG_PATHWAY	hsa04110	Cell cycle	13	0.00555	BUB1B, TTK, CDC6, CDC25B, CCNA2, CDC20, CCNB2, CCNB1, CCNE2, PTTG1, CCNE1, CDK1, MAD2L1
KEGG_PATHWAY	hsa05142	Chagas disease (American trypanosomiasis)	11	0.011418	C1QB, GNA15, PLCB4, NOS2, IL-1B, CCL2, CD3G, FOS, CD3D, MYD88, MAPK13
KEGG_PATHWAY	hsa04114	Oocyte meiosis	11	0.017478	CDC20, CCNB2, CCNB1, PTTG1, CALML5, CCNE2, CCNE1, CDK1, CALML3, MAD2L1, AURKA
KEGG_PATHWAY	hsa00592	Alpha-Linolenic acid metabolism	5	0.018996	PLA2G2F, FADS2, PLA2G4D, PLA2G2A, PLA2G3
KEGG_PATHWAY	hsa05168	Herpes simplex infection	15	0.020979	DDX58, STAT1, TAP1, FOS, IFIT1, SOCS3, OAS1, OAS2, IL-1B, OAS3, IRF7, CDK1, CCL2, MYD88, IRF9
KEGG_PATHWAY	hsa04380	Osteoclast differentiation	12	0.021031	FOSL1, SOCS3, FCGR3B, STAT1, LCK, IL-1B, BLNK, ACP5, FOS, JUNB, IRF9, MAPK13
KEGG_PATHWAY	hsa04064	NF-kappa B signaling pathway	9	0.028614	LYN, BCL2A1, DDX58, LCK, IL-1B, BLNK, PRKCQ, CCL19, MYD88
KEGG_PATHWAY	hsa04621	NOD-like receptor signaling pathway	7	0.028679	PYCARD, IL-1B, CARD6, CCL2, CXCL1, NOD2, MAPK13
KEGG_PATHWAY	hsa00591	Linoleic acid metabolism	5	0.031316	PLA2G2F, PLA2G4D, PLA2G2A, PLA2G3, CYP2E1
KEGG_PATHWAY	hsa00120	Primary bile acid biosynthesis	4	0.031945	CYP39A1, CH25H, ACOX2, CYP7B1
KEGG_PATHWAY	hsa00140	Steroid hormone biosynthesis	7	0.033377	SULT2B1, HSD11B1, STS, HSD3B1, HSD17B2, CYP2E1, CYP7B1
KEGG_PATHWAY	hsa01100	Metabolic pathways	64	0.037857	ST6GALNAC1, GDA, GLDC, PIK3C2G, PYGL, HK2, SMPD3, IL4I1, HPGDS, PNP, ACADL, SPTLC2, KYNU, NAMPT, LIPG, HYAL4, PGM2, TK1, UPP1, AASS, PLA2G4D, ARG1, GPT2, AMPD3, NME1, ALDH3A1, PLCB4, ACOX2, CMPK2, ACSBG1, CYP2E1, IDO1, PLA2G3, ALOX12B, FUT2, CYP2C18, FUT1, FUT3, TYMP, HSD11B1, UGCG, RDH12, HSD17B2, RDH16, ATP6V0A4, XDH, HAO2, PLA2G2F, GALNT6, RRM2, GK, GCH1, NOS2, PLA2G2A, HSD3B1, CYP4F8, ALDH4A1, SQLE, HAL, AKR1B10, DHRS9, POLE2, POLR3G, HPSE
KEGG_PATHWAY	hsa05143	African trypanosomiasis	5	0.047416	PLCB4, IL-1B, SELE, MYD88, IDO1

WGCNA

WGCNA was performed on 837 miRNAs. The samples were clustered to remove abnormal samples (GSM3711960 and GSM3711959) (Figure 4A and B). After calculation, when the soft threshold β was 2, it was approximately a scale-free topology (Figure 4C). After determining the soft threshold, the cluster tree graph was constructed. Subsequently, with the minimum number of genes for modules set to 20, the modules with dissimilarity < 25% merged using the dynamic tree cutting method. Finally, 6 modules were determined (Figure 4D).

Figure 4 Construction of weighted gene co-expression network. A: Sample screening dendrogram; B: Sample dendrogram and trait heat map; C: Scale-free fit index of different soft threshold power (β) and mean connectivity of various soft threshold power; D: The cluster dendrogram of the miRNAs. Each branch represents one gene, and every color below represents one module.

Uncovering of hub modules and miRNAs

The correlation between modules and psoriasis was analyzed by Pearson’s method. The results showed that turquoise module had the highest correlation with psoriasis (r = 0.96) (Figure 5A). Therefore, turquoise module was considered the hub module. Subsequently, 21 miRNAs were screened out from the turquoise module as candidate hub miRNAs (MM > 0.8 and GS > 0.5) (Figure 5B). Intersection of DEmiRNAs and candidate hub miRNAs was obtained (Figure 5C). Ten intersecting miRNAs were identified as real hub miRNA. Among them, 5 hub miRNAs were up-regulated (hsa-miR-21-3p, hsa-miR-21-5p, hsa-miR-31-5p, hsa-miR-18a-5p, and hsa-miR-33b-3p) and 5 hub miRNAs were down-regulated (hsa-miR-181a-2-3p, hsa-miR-181a-5p, hsa-miR-6510-3p, hsa-miR-30a-3p, and hsa-miR-30c-2-3p).

Figure 5  Hub module and hub miRNAs. A: Heat map of the correlation between module eigengene and psoriasis. The correlation between turquoise module and psoriasis was highest; B: Scatter plot of miRNAs in the turquoise module. MiRNAs with module connectivity > 0.8 and clinical trait relationship > 0.5 were selected as candidate hub miRNAs; C: The 10 intersecting miRNAs were identified as real hub miRNAs.

MiRNA-mRNA regulatory network

Target mRNAs of 10 hub miRNAs were predicted using miRDB database, but only 8 hub miRNAs (4 up-regulated and 4 down-regulated) were predicted to the corresponding target mRNAs. Eleven target mRNAs of up-regulated hub miRNAs (hsa-miR-21-3p, hsa-miR-21-5p, hsa-miR-18a-5p, and hsa-miR-33b-3p) were matched with DEmRNAs, and 72 target mRNAs of down-regulated hub miRNAs (hsa-miR-181a-2-3p, hsa-miR-181a-5p, hsa-miR-30a-3p, and hsa-miR-30c-2-3p) were matched with DEmRNAs. Ninety-seven negatively regulated miRNA-mRNAs were involved in the miRNA-mRNA regulatory network (Figure 6); for example, hsa-miR-21-3p/hsa-miR-18a-5p-F3, hsa-miR-21-5p-claudin 8 (CLDN8), hsa-miR-33b-3p-PLCB4, hsa-miR-30a-3p-IL-1B, hsa-miR-181a-5p-C-C motif chemokine ligand 8 (CCL8), hsa-miR-181a-5p/hsa-miR-30c-2-3p-C-X-C motif chemokine ligand 9 (CXCL9), and hsa-miR-30c-2-3p-KYNU.

Figure 6 Regulatory network of miRNA-mRNA. Red, yellow, oval, and V-shaped represent up-regulated, down-regulated, differentially expressed mRNAs, and hub differentially expressed miRNAs, respectively.

Functional enrichment of target DEmRNAs

Functional analysis of 83 target DEmRNAs was performed using David database (P < 0.05). GO enrichment results demonstrated that most of the target DEmRNAs were distributed in signal transduction (GO:BP), immune response (GO:BP), integral component of membrane (GO:CC), plasma membrane (GO:CC), ATP binding (GO:MF), and other biological functions (Figure 7A-C). KEGG enrichment results demonstrated that the target DEmRNAs were significantly enriched in influenza A, hepatitis C, measles, and chemokine signaling pathway (Figure 7D and Table 4).

Figure 7 Gene Ontology and Kyoto Encyclopedia of Genes and Genomics analysis of target differentially expressed mRNAs. A: biological processes terms were enriched of target differentially expressed mRNAs (DEmRNAs) by Gene Ontology (GO) function analysis; B: Cellular components terms were enriched of target DEmRNAs by GO function analysis; C: Molecular functions terms were enriched of target DEmRNAs by GO function analysis; D: Kyoto Encyclopedia of Genes and Genomics enrichment of target DEmRNAs.

Table 4 Kyoto Encyclopedia of Genes and Genomics function enrichment of target differentially expressed mRNAs.

Category	hsa	Term	Count	P value	Genes
KEGG_PATHWAY	hsa05164	Influenza A	6	0.002804	STAT1, OAS2, IL-1B, OAS3, TNFSF10, MAPK13
KEGG_PATHWAY	hsa05160	Hepatitis C	5	0.006506	STAT1, OAS2, CLDN8, OAS3, MAPK13
KEGG_PATHWAY	hsa05162	Measles	5	0.006506	STAT1, OAS2, IL-1B, OAS3, TNFSF10
KEGG_PATHWAY	hsa04062	Chemokine signaling pathway	5	0.020355	CXCL9, CCL8, PLCB4, CCL22, STAT1
KEGG_PATHWAY	hsa04620	Toll-like receptor signaling pathway	4	0.021751	CXCL9, STAT1, IL-1B, MAPK13
KEGG_PATHWAY	hsa04060	Cytokine-cytokine receptor interaction	5	0.047526	CXCL9, CCL8, CCL22, IL-1B, TNFSF10

RT-PCR validation

IL-1B, CXCL9, CLDN8, CCL8, hsa-miR-181a-5p, hsa-miR-30a-3p, and hsa-miR-21-5p related to psoriasis were selected from hub miRNAs and target DEmRNAs for in vitro validation. Primers for each mRNAs and miRNAs were shown in Table 5. IL-1B, CXCL9 and CCL8 were up-regulated and CLDN8 and hsa-miR-181a-5p were down-regulated in psoriasis tissues (Figure 8). Among them, IL-1B, CXCL9, and CLDN8 showed significant difference in expression levels. In addition, we found that expression trends of hsa-miR-30a-3p and hsa-miR-21-5p were contrary to the results of bioinformatics analysis. The reason for the inconsistency between RT-PCR and bioinformatics analysis results may be the small sample size. Further research is needed.

Figure 8 Real-time polymerase chain reaction validation. A: Interleukin-1β; B: C-X-C motif chemokine ligand 9; C: C-C motif chemokine ligand 8; D: Claudin 8; E: hsa-miR-181a-5p; F: hsa-miR-21-5p; G: hsa-miR-30a-3p. ^aP < 0.05, P < 0.05 was considered significant.

Table 5 Primer sequence in the real-time polymerase chain reaction.

Primer name	Primer sequence, 5’ to 3’
GAPDH-F (internal reference)	5-GGAGCGAGATCCCTCCAAAAT-3
GAPDH-R (internal reference)	5-GGCTGTTGTCATACTTCTCATGG-3
ACTB-F (internal reference)	5-CATGTACGTTGCTATCCAGGC-3
ACTB-R (internal reference)	5-CTCCTTAATGTCACGCACGAT-3
IL-1B-F	5-TTCGACACATGGGATAACGAGG-3
IL-1B-R	5-TTTTTGCTGTGAGTCCCGGAG-3
CXCL9-F	5-CCAGTAGTGAGAAAGGGTCGC-3
CXCL9-R	5-AGGGCTTGGGGCAAATTGTT-3
CCL8-F	5-TGGAGAGCTACACAAGAATCACC-3
CCL8-R	5-TGGTCCAGATGCTTCATGGAA-3
CLDN8-F	5-CTACAGGCAGCCAGAGGACT-3
CLDN8-R	5-ACAGGGATGAGCACCACCAT-3
hsa-U6 (internal reference)
hsa-miR-30a-3p	5-CTTTCAGTCGGATGTTTGCAGC-3
hsa-miR-181a-5p	5-AACATTCAACGCTGTCGGTGAGT-3
hsa-miR-21-5p	5-TAGCTTATCAGACTGATGTTGA-3

DISCUSSION

Although psoriasis is widespread and has significant negative impact on patients’ quality of life, it has not yet been fully diagnosed and treated. In the research, we used the WGCNA method to identify 10 hub miRNAs that may be related to the pathogenesis of psoriasis. Then, target mRNAs of hub miRNAs were predicted by using miRDB database. Only 8 hub miRNAs were predicted to the corresponding target mRNAs. Subsequently, to understand the key biological functions involved in DEmRNAs and target DEmRNAs, we performed GO and KEGG functional analysis. Results demonstrated that they were significantly enriched in immune-related biological functions, for example, immune response, cell chemotaxis, and chemokine signaling pathway.

Hsa-miR-21-5p is abnormally expressed in psoriasis, but the specific molecular mechanism remains unclear[29]. In this study, results demonstrated that hsa-miR-21-5p expression was up-regulated and negatively regulated with CLDN8. CLDN8 is down-regulated in psoriasis and has an important molecular regulatory role[30]. CLDN8 has also been identified as a key downstream component of the IL-9 and IL-23 inflammatory cascade[31,32]. In addition, KEGG functional enrichment analysis showed that CLDN8 was enriched in hepatitis C. Cohen et al[33] found that psoriasis is associated with hepatitis C[33]. Cathelicidin, toll like receptor 9 (TLR9), IFN-γ, and TNF-α are inflammatory cytokines. Hepatitis C may increase susceptibility to psoriasis by up-regulating these inflammatory factors[34,35]. Thus, we hypothesized that hsa-miR-21-5p may play a vital regulatory role in hepatitis C through regulating CLDN8 and, thus, affect the pathogenesis of psoriasis.

Hsa-miR-30a-3p can affect the migration and proliferation of cancer cells by targeting related genes[36,37]. Hsa-miR-30a-3p is also contacted with platelet apoptosis and adhesion in immune thrombocytopenia[38]. However, we have not found any report about hsa-miR-30a-3p in psoriasis. In this study, we found that hsa-miR-30a-3p was up-regulated and negatively regulated with multiple DEmRNAs in psoriasis. It may play a vital regulatory role in psoriasis by targeting these DEmRNAs. We found the pro-inflammatory factor IL-1B in target DEmRNAs of hsa-miR-30a-3p. IL-1B has been found to play an important role in autoimmune or autoinflammatory conditions[39]. IL-1B is abundant in the tissue fluid of psoriasis and participates in the disease progression through dual effects[40]. First, it induces insulin resistance through p38 mitogen-activated protein kinase (p38MAPK), preventing insulin-dependent differentiation of keratinocytes, while IL-1B promotes keratinocyte proliferation, both hallmarks of psoriasis. Through KEGG functional enrichment analysis, we also found that IL-1B participates in multiple signal pathways, for example, measles, toll-like receptor signaling pathway, and cytokine-cytokine receptor interaction. Measles can relieve psoriasis through an immunosuppressive effect[41]. Toll-like receptor signaling pathway play an important role in psoriasis by affecting keratinocyte production[42]. The cytokine-cytokine receptor interaction is related to the occurrence and progression of psoriasis through combined transcriptomic analysis[43], which is consistent with our analysis. Thus, we hypothesized that hsa-miR-30a-3p may play a vital molecular role in the progression of psoriasis by targeting DEmRNAs to regulate multiple biological signaling pathways.

Hsa-miR-181a-5p is involved in the catabolic pathway of chondrocytes and oxidative stress in osteoarthritis[44]. Hsa-miR-181a-5p can also regulate the pathogenesis of sepsis-related inflammation through CRNDE/hsa-miR-181a-5p/TLR4 pathway[45]. In a case-control study, hsa-miR-181a-5p was significantly down-regulated in psoriasis[29]. The miRNA-mRNA regulatory network results demonstrated that hsa-miR-181a-5p was negatively regulated with multiple DEmRNAs. Among these DEmRNAs, we found inflammatory mediators CXCL9 and CCL8. CXCL9 is an important chemokine involved in T cell recruitment and is up-regulated in the plasma of patients with psoriasis[46,47]. Increased CXCL9 expression can aggravate the progression of psoriasis[48]. CXCL9 and hsa-miR-30c-2-3p were also negatively regulated in the miRNA-mRNA regulatory network. Oncology studies have shown that hsa-miR-30c-2-3p play a vital role in tumor pathogenesis by regulating the proliferation, apoptosis, migration, and invasion of cancer cells[49,50]. So far, we have not found any studies on hsa-miR-30c-2-3p in psoriasis. This article may be the first to report that hsa-miR-30c-2-3p plays a role in the pathogenesis of psoriasis. As a chemokine, CCL8 is involved in immune regulation and inflammatory processes in a variety of diseases[51-53]. Although no relevant studies on CCL8 have been found in psoriasis, the expression of CCL8 is up-regulated in atopic dermatitis[54]. In addition, CXCL9 and CCL8 were found to be enriched in chemokine signaling pathway and cytokine-cytokine receptor interaction by functional analysis. Therefore, this further suggests that hsa-miR-181a-5p and hsa-miR-30c-2-3p may play a regulatory role in psoriasis by targeting DEmRNAs to mediate multiple biological signaling pathways.

Current research results highlighted that silencing hsa-miR-181a-2-3p could enhance cadmium-induced inflammatory response and activation of inflammasome[55]. In the research, we found that hsa-miR-181a-2-3p was negatively regulated with multiple DEmRNAs. Among these DEmRNAs, OAS3 was found to be involved in hepatitis C and measles in KEGG functional enrichment. Some researchers have found that the OAS3 is concerned with the occurrence and progression of psoriasis through transcriptomic analysis, which is consistent with our analysis[56,57]. Thus, we hypothesized that hsa-miR-181a-2-3p may play a regulatory role in the progression of psoriasis by targeting OAS3 to mediate hepatitis C and measles. This provides potential molecular research directions for further research on the pathogenesis of psoriasis.

The expression of hsa-miR-21-3p in psoriasis was significantly increased[58]. Moreover, hsa-miR-21-3p plays a pro-inflammatory role in keratinocytes, and high expression in the skin is concerned with psoriasis[59]. So far, we have not found relevant reports about hsa-miR-18a-5p in psoriasis. However, hsa-miR-18a-5p promotes the proliferation and migration of pulmonary smooth muscle cells by targeting notch receptor 2[60]. Oncology studies have shown that hsa-miR-18a-5p promotes melanoma cell proliferation by targeting EPH receptor A7 signaling[61]. Hsa-miR-18a-5p can affect keratinocytes apoptosis by targeting B-cell lymphoma/leukemia-2-like protein 10 in patients with toxic epidermal necrolysis and is related to the skin erythema or erosion area of drug eruptions[62]. The miRNA-mRNA regulatory network regulatory network result demonstrates that hsa-miR-21-3p and hsa-miR-18a-5p jointly negatively regulate coagulation factor III, tissue factor (F3). This finding provides new research ideas for the pathogenesis of psoriasis in the future.

Results of previous studies demonstrate that down-regulation of interleukin 1 receptor associated kinase 3 by hsa-miR-33b-3p can alleviate the inflammation and apoptosis induced by IL-1B in CHON-001 cells[63]. As a key miRNA of major depression disorder and Kawasaki disease, hsa-miR-33b-3p play an important role in their pathogenesis[64,65]. Hsa-miR-33b-3p has also been reported in cancer[66]. In the miRNA-mRNA regulatory network, phospholipase C beta 4 (PLCB4) is the only negatively regulated target DEmRNA of hsa-miR-33b-3p. Neutrophils are an important part of the innate immune system and an early line of defense against microbial invasion. PLCB4 can regulate the number of neutrophils[67]. In addition, KEGG functional analysis result demonstrated that PLCB4 was enriched in chemokine signaling pathway. Thus, we hypothesized that hsa-miR-33b-3p may play a vital molecular role in psoriasis by targeting PLCB4 to regulate chemokine signaling pathway.

This study has some limitations. First, sample size of in vitro validation was small, leading to a certain degree of error between RT-PCR validation results and bioinformatics analysis results. Further studies with a larger sample size are needed. Secondly, the specific regulatory mechanism of the identified genes and signaling pathways in psoriasis remain unclear, so further research is needed.

CONCLUSION

In conclusion, identification of potential key molecular markers and signaling pathways provides potential molecular research directions for further understanding the pathological mechanisms of psoriasis. This may also provide new research ideas for the prevention and treatment of psoriasis in the future.

ARTICLE HIGHLIGHTS

Research background

Previous studies have found that microRNAs (miRNAs) play an important regulatory role in the progression of various diseases. Currently, miRNAs studies in psoriasis and dermatology are relatively new.

Research motivation

Although psoriasis is widespread and has significant negative impact on patients’ life quality, it has not yet been fully diagnosed and treated.

Research objectives

Identification of key miRNAs in psoriasis is helpful to elucidate the molecular mechanism of psoriasis.

Research methods

Differentially expressed mRNAs (DEmRNAs) and differentially expressed miRNAs were screened out by limma R package. DEmRNAs were analyzed for Gene Ontology and Kyoto Encyclopedia of Genes and Genomics functional enrichment. The “Weighted gene co-expression network analysis (WGCNA)” R package was used to analyze the co-expression network of all miRNAs. We constructed miRNA-mRNA regulatory networks based on identified hub miRNAs.

Research results

We identified a large number of DEmRNAs and screened possible signaling pathways related to psoriasis, for example, toll-like receptor signaling pathway, cytokine-cytokine receptor interaction, and chemokine signaling pathway. Ten hub miRNAs were identified by WGCNA. Eight hub miRNAs predicted the corresponding target mRNAs. Ninety-seven negatively regulated miRNA-mRNA pairs were involved in the miRNA-mRNA regulatory network, for example, hsa-miR-21-5p-CLDN8, hsa-miR-30a-3p-IL-1B and hsa-miR-181a-5p/hsa-miR-30c-2-3p-CXCL9.

Research conclusions

The identification of potential key molecular markers and signaling pathways provides potential research directions for further understanding the molecular mechanisms of psoriasis.

Research perspectives

This study provide new research ideas for the prevention and treatment of psoriasis in the future.