
==== Front
Medicine (Baltimore)
Medicine (Baltimore)
MD
Medicine
0025-7974
1536-5964
Lippincott Williams & Wilkins Hagerstown, MD

39028999
MD-D-23-07148
00001
10.1097/MD.0000000000039002
3
4000
Research Article
Observational Study
Integrative analysis of gene and microRNA expression profiles reveals candidate biomarkers and regulatory networks in psoriasis
Chen Lu MS 1149492120@qq.com
a
Wang Xiaochen MS 1402243394@qq.com
a
Liu Chang MS 1114060791@qq.com
a
Chen Xiaoqing BS chenxqing@jhun.edu.cn
a
Li Peng MD doglp2005@163.com
b
https://orcid.org/0000-0002-9267-084X
Qiu Wenhong MD a*
Guo Kaiwen MD guowendy@wust.edu.cn
c
a Department of Immunology, Jianghan University, School of Medicine, Wuhan, Hubei, PR China
b Department of Dermatology, Wuhan Central Hospital, Wuhan, Hubei, PR China
c Department of Pathogenic Biology, Wuhan University of Science and Technology, Medical College, Wuhan, Hubei, PR China.
* Correspondence: Wenhong Qiu, Department of Immunology, Jianghan University, School of Medicine, Wuhan, Hubei 430056, PR China (e-mail: qiuwenhong@jhun.edu.cn).
19 7 2024
19 7 2024
103 29 e3900219 8 2023
27 3 2024
28 6 2024
Copyright © 2024 the Author(s). Published by Wolters Kluwer Health, Inc.
2024
https://creativecommons.org/licenses/by-nc/4.0/ This is an open-access article distributed under the terms of the Creative Commons Attribution-Non Commercial License 4.0 (CCBY-NC), where it is permissible to download, share, remix, transform, and buildup the work provided it is properly cited. The work cannot be used commercially without permission from the journal.

Psoriasis (PS) is a chronic inflammatory skin disease with a long course and tendency to recur, the pathogenesis of which is not fully understood. This article aims to identify the key differentially expressed genes (DEGs) and microRNA (miRNAs) of PS, construct the core miRNA-mRNA regulatory network, and investigate the underlying molecular mechanism through integrated bioinformatics approaches. Two gene expression profile datasets and 2 miRNA expression profile datasets were downloaded from the gene expression omnibus (GEO) database and analyzed by GEO2R. Intersection DEGs and intersection differentially expressed miRNAs (DEMs) were each screened. The Metascape database and R software were used to perform enrichment analysis of intersecting DEGs and study their functions. Target genes of DEMs were predicted from the online database miRNet. The protein-protein interaction files of the overlapping target genes were obtained from string and the miRNA-mRNA network was constructed by Cytoscape software. In addition, the online web tool CIBERSORT was used to analyze the immune infiltration of dataset GSE166388, and the relative abundance of 22 immune cells in the diseased and normal control tissues was calculated and assessed. Finally, quantitative reverse transcription polymerase chain reaction (qRT-PCR) was used to verify the relative expression of the screened miRNAs and mRNAs to assess the applicability of DEMs and DEGs as biomarkers in PS. A total of 205 mating DEGs and 6 mating DEMs were screened. 103 dysregulated crossover genes from 205 crossover DEGs and 7878 miRNA target genes were identified. The miRNA-mRNA regulatory network was constructed and the top 10 elements were obtained from CytoHubba, including hsa-miR-146a-5p, hsa-miR-17-5p, hsa-miR-106a-5p, hsa-miR-18a-5p, CDK1, CCNA2, CCNB1, MAD2L1, RRM2, and CCNB2. QRT-PCR revealed significant differences in miRNA and gene expression between inflammatory and normal states. In this study, the miRNA-mRNA core regulator pairs hsa-miR-146a-5p, hsa-miR-17-5p, hsa-miR-106a-5p, hsa-miR-18a-5p, CDK1, CCNA2, CCNB1, MAD2L1, RRM2, and CCNB2 may be involved in the course of PS. This study provides new insights to discover new potential targets and biomarkers to further investigate the molecular mechanism of PS.

biomarkers
miRNAs
psoriasis
regulatory networks
OPEN-ACCESSTRUE
==== Body
pmc1. Introduction

Psoriasis (PS) is a common immune-mediated inflammatory skin disease with a long course and a tendency to recur.[1] The disease can occur in any country of the world and affects people of all ages, but mostly young adults. The disease can have a significant impact on patient’s physical and mental health. The clinical manifestations are mainly erythema and scales, and the whole body can be affected. Scalps and extended limbs are more common, and most of these get worse in winter. At the same time, it is associated with several important diseases, including depression, psoriatic arthritis, and cardiometabolic syndrome.[2] Because it is a chronic relapsing disease and there is currently no specific therapy, many patients need long-term treatment to control the disease recurrence, and all kinds of therapies have certain side effects. At present, there are mainly local drug therapy, phototherapy, oral systemic therapy, biological therapy and so on.[2] In recent years, the rapid development of high-throughput technologies has enabled the identification of many disease markers that can aid in the diagnosis and treatment of patients. It remains the key to fully utilizing effective biomarkers or targets to improve clinical diagnosis and management of patients with PS.

MiRNAs (microRNAs) are small, highly conserved non-coding RNA sequences of 19 to 25 nucleotides. This single-stranded RNA can regulate the expression of protein-coding genes at the post-transcriptional level, participate in the maintenance of normal cell homeostasis, and play an important role in various biological processes.[3] Increasing evidence shows that miRNAs are successfully used as biomarkers of PS for diagnosis, prognosis and monitoring of treatment success.[4]

With the development of bioinformatics, various algorithms and research strategies have been applied to identify the underlying potential regulatory mechanisms of gene networks, which has led to a comprehensive and profound understanding of many diseases and can be applied to the screening of genetic genome and transcriptome alterations, leading to the determination differential expressed genes (DEGs) and their functions.[5,6] At present, most of the methods of disease management and diagnosis of PS are based on clinical evaluation.[7] Due to factors such as individual gene expression differences, lifestyle differences, and environmental influences, the evaluation results are usually inaccurate, and there is still a lack of objective markers that accurately reflect the severity of PS.[7] Transcriptome profiling techniques based on bioinformatics algorithms, such as RNA-sequencing and microarrays, can more accurately measure transcription levels and their subtypes, effectively identifying objective biomarkers of PS.

Therefore, this study obtained multiple gene expression profiles and miRNA expression profiles of clinical PS patients in the gene expression omnibus (GEO) database based on microarrays in transcriptome analysis technology. Co-differentially expressed genes (DEGs) and co-differentially expressed miRNAs (DEMs) in the skin of PS patients were integrated, and key genes and miRNAs that may serve as novel biomarkers were further identified by constructing the miRNA-mRNA regulatory network and finally verified in the inflammatory cell model. Based on clinical databases, transcriptome profiling techniques are used for integrated identification to provide more reliable new insights for PS toward the diagnosis and treatment goals of precision medicine.

2. Materials and methods

2.1. Microarray data

Four data sets with gene expression profiles from PS and healthy controls (GSE166388, GSE153007, GSE175438, and GSE145305) were downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo). GSE166388 microarray data included mRNA expression profiles from 4 PS samples and 4 normal samples, and GSE153007 microarray data included mRNA expression profiles from 14 PS samples and 5 normal samples. GSE175438 microarray data contained miRNA expression profiles from 11 PS samples and 9 normal samples, while GSE145305 microarray data contained miRNA expression profiles from 4 PS samples and 4 normal samples. In these 4 datasets, all PS samples were obtained from diseased skin biopsies from PS patients and control samples were obtained from normal skin biopsies from healthy volunteers.

2.2. Identification of DEGs and DEMs

GEO2R (http://www.ncbi.nlm.nih.gov/geo/geo2r/) is an online, interactive tool based on the Limma package to identify DEGs by comparing samples from the GEO series. We used GEO2R to identify mRNAs (DEGs) and miRNAs (DEMs). DEG screening criteria: |log2FC|≥1 and adjusted P value < .05; DEMs screening criteria: |log2FC|≥1 and adjusted P value < .05. Results were visualized through heat maps and volcano plots using online bioinformatics tools (http://www.sangerbox.com/tool).

The PS-related mRNAs were obtained by section DEGs of 2 datasets (GSE166388 and GSE153007). The PS-related miRNAs were obtained by section DEMs of 2 datasets (GSE175438 and GSE145305). Intersection DEGs and DEMs were obtained via the online bioinformatics tool VENNY 2.1.0. (https://bioinfogp.cnb.csic.es/tools/venny/).

2.3. GO and KEGG enrichment analysis

Gene ontology (GO) is a bioinformatics tool for gene annotation in 3 parts: molecular function (MF), biological process (BP), and cell component (CC). The Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis can perform functional analysis and important pathway analysis of DEGs.[8–10] GO function annotation and KEGG pathway enrichment analysis were performed for all overlapping DEGs via an online bioinformatics database called Metascape database (https://metascape.org/gp/index.html#/main/step1).

For overlapping upregulated DEGs and downregulated DEGs, GO enrichment analysis was performed with R package clusterProfiler (version 3.14.3). Functional GO enrichment results were obtained. P < .05 and FDR < 0.1 were considered statistically significant.

For overlapping upregulated DEGs and downregulated DEGs, the clusterProfiler R package (version 3.14.3) was used for enrichment analysis to obtain KEGG enrichment results. P < .05 and FDR < 0.1 were considered statistically significant.

2.4. Prediction of miRNAs targets

The downstream target genes of overlapping DEMs in datasets GSE175438 and GSE145305 were predicted by miRNet (https://www.mirnet.ca/miRNet/home.xhtml), an online bioinformatics tool. The miRNet website predicts target genes based on 2 algorithms: miRanda and TarPmiR. The MiRanda algorithm mainly emphasizes the evolutionary conservation of miRNA and target gene link sites and uses RNAFold to calculate thermodynamic stability, with score = 140 and energy = 1 as the default settings. The TarPmiR algorithm sets the probability cut > 0.5. We compared the predicted target genes to cut difference mRNAs (DEGs) and used the online bioinformatics tool VENNY 2.1.0 to obtain the cut points, which were used as cut point-dysregulated target genes for further assembly of regulatory networks.

2.5. MiRNA-mRNA network construction and module analysis

Target gene protein-protein interaction (PPI) files were obtained from string online database (https://cn.string-db.org/), high confidence > 0.7. PPI files, miRNA, and downstream target gene correspondence files were sequentially merged and imported into Cytoscape 3.9.1 bioinformatics software for visual analysis to construct the interaction network between differential miRNAs and differential mRNAs. We used the CytoHubba plugin in Cytoscape to calculate and select the top 10 elements in the network based on the Degree algorithm.

2.6. Immuneinfiltration analysis

Immunoinfiltration analysis was performed on dataset GSE166388 using the online web tool CIBERSORT (https://cibersortx.stanford.edu/). CIBERSORT is an immunoinfiltration analysis tool developed by a team of researchers at Stanford University. Based on transcriptome data, CIBERSORT uses a deconvolution algorithm to estimate the composition and relative abundance of immune cells in a mixed-cell population with high precision and confidence. It can fully reflect the infiltration of immune cells in the disease samples and is conducive to the study of the immune mechanism of the occurrence and development of diseases. At the same time, CIBERSORT has certain limitations: his analysis results can be affected by sample quality, data sources, and technical biases. The dataset GSE166388 included 4 disease samples (GSM5070321, GSM5070322, GSM5070323, and GSM5070324) and 4 healthy controls (GSM5070317, GSM5070318, GSM5070319, and GSM5070320). The relative abundances of 22 immune cells in diseased and normal control tissues were calculated and evaluated. Run mode is relative, permutations are set to 1000.

2.7. Cell culture and inflammatory cell model construction

The human keratinocyte lineage immortalized HaCaT cells (Procell, Wuhan, China) were cultured in a complete medium supplemented with 10% FBS (Procell, Wuhan, China) + 90% MEM (Procell, Wuhan, China), 37.5% CO2 was configured. The cells were seeded in 6-well plates and cultured overnight. The next day, the medium was discarded, washed with PBS, and a fresh medium containing 10 ng/mL TNF-α (PEPROTECH, USA) and IFN-γ (PEPROTECH, USA) was added. The inflammatory cell model was obtained by continuing the cultivation for 24 hours and the expression levels of IL-6 and IL-1β were detected to determine whether the inflammatory cell model was successfully constructed.

2.8. Quantitative reverse transcription polymerase chain reaction (qRT-PCR) analysis

Total RNA in cells was extracted and reverse transcribed with Trizol (Ambion, USA) according to the manufacturer instructions. QRT-PCR was performed on the CFX Connect Real-time PCR Detection System (Bio-Rad, USA). U6 and GAPDH were selected to normalize miRNA and gene expression levels. Relative expression was calculated using the 2−ΔΔCt method. The primer sequences used in the qRT-PCR analysis are listed in Table 1.

Table 1 Primers used in qRT-PCR.*

Primers	Sequences(5，→ 3，)	
hsa-miR-146a-5p-F	CGCGTGAGAACTGAATTCCA	
hsa-miR-146a-5p-R	ATCCAGTGCAGGGTCCGAGG	
hsa-miR-18a-5p-F	CCTGTGCATAAGGTGCATCTAGTG	
hsa-miR-18a-5p-R	ATCCAGTGCAGGGTCCGAGG	
hsa-miR-106a-5p-F	GGTCCGGAAAAGTGCTTACAGTG	
hsa-miR-106a-5p-R	ATCCAGTGCAGGGTCCGAGG	
hsa-miR-17-5p-F	CCTCTGCCAAAGTGCTTACAGTG	
hsa-miR-17-5p-R	ATCCAGTGCAGGGTCCGAGG	
U6-F	AGAGAAGATTAGCATGGCCCCTG	
U6-R	AGTGCAGGGTCCGAGGTATT	
CDK1-F	AAACTGGCTGATTTTGGCCTTG	
CDK1-R	GTTGAGTAACGAGCTGACCC	
CCNA2-F	ACCGTTCCTCCTTGGAAAGC	
CCNA2-R	CAGGGCATCTTCACGCTCTAT	
CCNB1-F	AAGGCTGTGGCAAAGGTGTAA	
CCNB1-R	ACAACTTTTCCGCTACCCTACA	
MAD2L1-F	AATCGTGGCCGAGTTCTTCT	
MAD2L1-R	TTACAAGCAAGGTGAGTCCGT	
RRM2-F	GCCACACCATGAATTGTCCG	
RRM2-R	ATGGTAAGTCACAGCCAGCC	
CCNB2-F	GGAAGTCATGCAGCACATGG	
CCNB2-R	TTCCTATCAGTGGGGAGGCA	
IL-6-F	TTCGGTCCAGTTGCCTTCTC	
IL-6-R	CTGAGATGCCGTCGAGGATG	
IL-1β-F	AAGTACCTGAGCTCGCCAGT	
IL-1β-R	CTTGCTGTAGTGGTGGTCGG	
GAPDH-F	AATTCCATGGCACCGTCAAG	
GAPDH-R	AGCATCGCCCCACTTGATTT	
* qRT-PCR = quantitative reverse transcription polymerase chain reaction.

2.9. Statistical analysis

All data were expressed as mean ± standard deviation (SD). Statistical analyses were calculated using GraphPad Prism (version 8.0.1). The t test was used for statistical test. P < .05 was considered statistically significant.

3. Results

3.1. Identification of DEGs in PS

The original data of 2 independent datasets (GSE166388 and GSE153007) were downloaded from the GEO database; DEGs were identified by GEO2R. A total of 547 DEGs (366 upregulated and 181 downregulated mRNAs) and 2232 DEGs (1108 upregulated and 1124 downregulated mRNAs) were screened from the GSE166388 and GSE153007 data sets, respectively. The expression of different mRNAs between PS samples and healthy controls was visualized by heatmaps (Fig. 1A and B). Simultaneously, these DEGs were visualized in the volcano plots (Fig. 1C and D). In addition, the online bioinformatics tool VENNY 2.1.0 was used to obtain 205 slice DEGs, including 159 upregulated mRNAs and 46 downregulated mRNAs (Fig. 1E).

Figure 1. DEGs from the GSE166388 and GSE153007. (A) Heat map of GSE166388 differentially expressed genes. (B) Heat map of GSE153007 differentially expressed genes. (C) Volcano plot of GSE166388 differentially expressed genes. (D) Volcano plot of GSE153007 differentially expressed genes. (E) Intersection DEGs of GSE166388 and GSE153007. DEGs = differentially expressed genes.

3.2. GO and KEGG enrichment analysis

First, Metascape was used for GO and KEGG analysis to examine the functions and paths of all 205 identified intersection DEGs. The results of the GO functional annotation of the Top 20 were obtained by analysis (Fig. 2A) and mainly include chromosome segregation, innate immune response, cyclin-dependent protein kinase holoenzyme complex, DNA metabolic process, meiotic chromosome segregation, regulation of DNA metabolic process, response to type I interferon, GTP binding, positive regulation of cell death, protein homodimerization activity, DNA replication, cytokine-mediated signaling pathway, interleukin-27-mediated signaling pathway, melanosome, positive regulation of cell cycle process, ubiquitin-like protein ligase binding, proteasome core complex, alpha-subunit complex regulation of peptidase activity, regulation of deoxyribonuclease activity and protein domain specific binding. The top 7 KEGG pathways mainly include cell cycle, hepatitis C, Epstein-Barr virus infection, proteasome, apoptosis, viral lifecycle-HIV-1, IL-17 signaling pathway (Fig. 2B).

Figure 2. GO and KEGG enrichment analysis by Metascape database. (A) Top20 results of 205 intersection DEGs GO analysis. (B) Top7 results of 205 intersection DEGs KEGG analysis.[8–10] DEGs = differentially expressed genes, GO = gene ontology, KEGG = Kyoto encyclopedia of genes and genomes.

Then the top 5 GO annotations including BP, CC and MF and the Top 10 KEGG pathways were performed on 159 overlapping upregulated genes and 46 overlapping downregulated genes. Results showed that the most enriched terms for upregulated genes were innate immune response in BP, cytosol in CC and GTP binding in MF (Fig. 3A). KEGG pathways of upregulated genes were mainly enriched in cell cycle, hepatitis C and influenza A. The visualization of KEGG enrichment analysis of upregulated genes is shown in Figure 3B to D. At the same time, the most enriched terms for downregulated genes were positive regulation of carbohydrate metabolic process in BP, beta-catenin destruction complex in CC and alpha-actinin binding in MF (Fig. 4A). The downregulated genes were mainly enriched in insulin resistance, insulin signaling pathway and Wnt signaling pathway. KEGG enrichment analysis of downregulated genes was visualized in Figure 4B to D.

Figure 3. GO and KEGG analysis of upregulated DEGs, via the R package. (A) Top 5 significant enrichment terms of upregulated DEGs in GO annotation. Blue, red and green bars represent BP, CC and MF, respectively. The length of the bars represents the gene counts. (B) Top 10 significant enrichment terms of upregulated DEGs in KEGG pathways.[8–10] The depth of colors indicates the P value, and the size of circles shows the gene counts in the pathway. (C) Lollipop chart for Top 10 KEGG pathway. (D) Circos plot for Top 10 KEGG pathway. BP = biological process, CC = cellular component, DEGs = differentially expressed genes, GO = gene ontology, KEGG = Kyoto encyclopedia of genes and genomes, MF = molecular function.

Figure 4. GO and KEGG analysis of downregulated DEGs, via the R package. (A) Top 5 significant enrichment terms of downregulated DEGs in GO annotation. Green, red and blue bars represent BP, CC and MF, respectively. The length of the bars represents the gene counts. (B) Top 10 significant enrichment terms of downregulated DEGs in KEGG pathways.[8–10] The depth of colors indicates the P value, and the size of circles shows the gene counts in the pathway. (C) Lollipop chart for Top 10 KEGG pathway. (D) Circos plot for Top 10 KEGG pathway. BP = biological process, CC = cellular component, DEGs = differentially expressed genes, GO = gene ontology, KEGG = Kyoto encyclopedia of genes and genomes, MF = molecular function.

In summary, the results showed that the upregulated differential genes in PS play an important role in the regulation of the innate immune response, the interferon signaling pathway, measles, and the cell cycle pathway. These upregulated genes are involved in the inflammatory and immune responses of PS and regulate the proliferation, differentiation, and apoptosis of cells. The downregulated differential genes in PS play important roles in metabolism, the insulin pathway, and the Wnt signaling pathway. These downregulated genes are involved in mediating intracellular signaling and transcription factor regulation, regulating glucose transport, protein synthesis, cell proliferation, and related cell and tissue survival.

3.3. Identification of DEMs in PS

The original data from 2 independent data sets (GSE175438 and GSE145305) were downloaded from the GEO database; DEMs were identified by GEO2R. A total of 79 DEMs (all upregulated miRNAs) and 66 DEMs (35 upregulated and 31 downregulated miRNAs) were screened from the GSE175438 and GSE145305 datasets, respectively. The expression of different miRNAs between PS samples and healthy controls was visualized by heatmaps (Fig. 5A and B). Meanwhile, these DEMs have been visualized in the volcano plots (Fig. 5C and D). Furthermore, after using the online bioinformatics tool VENNY 2.1.0 to obtain the intersection, 6 DEMs were selected from the 2 datasets that were both upregulated (Fig. 5E, Table 2).

Table 2 The details of the overlapping DE-miRNAs.*

The DE-miRNAs		
Symbol	Up/Down	
hsa-miR-31-5p	up	
hsa-miR-146a-5p	up	
hsa-miR-18a-5p	up	
hsa-miR-106a-5p	up	
hsa-miR-1307-3p	up	
hsa-miR-17-5p	up	
* DE-miRNAs = differentially expressed miRNAs.

Figure 5. DEMs from the GSE175438 and GSE145305. (A) Heat map of GSE175438 differentially expressed miRNAs. (B) Heat map of GSE145305 differentially expressed miRNAs. (C) Volcano plot of GSE175438 differentially expressed miRNAs. (D) Volcano plot of GSE145305 differentially expressed miRNAs. (E) Intersection DEMs of GSE115293 and GSE145305. DEMs = differentially expressed miRNAs, miRNA = MicroRNA.

3.4. Prediction of miRNAs targets

The downstream target genes of 6 intersection DEMs were predicted from the online bioinformatics tool miRNet database, and 7878 target genes were predicted. Simultaneously, 205 crossover DEGs and miRNA predictive target genes were cut with VENNY 2.1.0 to yield 103 dysregulated crossover target genes for further analysis and construction of the miRNA-mRNA network (Fig. 6).

Figure 6. Intersection results of intersection DEGs and predicted target genes. DEGs = differentially expressed genes.

3.5. MiRNA-mRNA network construction and Hub gene analysis

PPI files obtained from string, miRNA, and downstream target gene correspondence files were merged and imported into Cytoscape 3.9.1 to construct the interaction network between PS-cut differential miRNAs and 103 target mRNAs (Fig. 7A). It can be seen that hsa-miR-146a-5p is associated with the most target genes, while hsa-miR-1307-3p is associated with the fewest target genes. The CytoHubba plugin in Cytoscape was used for the calculation and the Degree algorithm was used to extract the core network (Fig. 7B). The results showed that the top 10 elements including upregulated hsa-miR-146a-5p, upregulated hsa-miR-17-5p, upregulated hsa-miR-106a-5p and upregulated hsa-miR-18a-5p in PS and 6 hub genes with upregulated expression, including CCNB1, CCNA2, CDK1, MAD2L1, RRM2 and CCNB2 (Table 3).

Table 3 The details of the Top 10 elements.

node_name	Degree	Up/Down	
hsa-miR-146a-5p	58	up	
hsa-miR-17-5p	39	up	
CCNB1	25	up	
hsa-miR-106a-5p	24	up	
CCNA2	24	up	
CDK1	24	up	
hsa-miR-18a-5p	23	up	
MAD2L1	22	up	
RRM2	21	up	
CCNB2	21	up	

Figure 7. MiRNA-mRNA network construction and module analysis. (A) miRNA-mRNA regulatory network, triangle symbol represents miRNA, ellipse symbol represents mRNA, orange is upregulated, green is downregulated. (B) CytoHubba plugin extracted the core network. The red icon had the highest score, followed by the orange icon and the yellow icon had the lowest score. miRNA = microRNA.

3.6. Immune infiltration analysis

The results of the immune infiltration analysis of the dataset GSE166388 showed that the frequency of several immune cell groups was different in PS patients. These include naïve B cells, memory B cells, CD4 memory-activated T cells, follicular helper T cells, regulatory T cells, resting and activated natural killer cells, monocytes, M0, M1 and M2 macrophages, resting and activated dendritic cells, dormant and activated mast cells and neutrophils. Immunoinfiltration between healthy control samples and PS samples was visualized by stacking a bar graph, a heat map, and a boxplot (Fig. 8A, C, and D). The correlation heat map of 22 immune cells is shown in Figure 8B. In general, the infiltration levels of T cells, CD4 memory quiescence, Tregs, and dendritic cells that are quiescent were higher in the control group. In the PS group, CD8 T cells, follicular helper T cells, and activated dendritic cells had higher infiltration rates, while resting mast cells and Tregs had lower infiltration rates.

Figure 8. Immune infiltration analysis of GSE166388. (A) Stacked bar chart of GSE166388 immune infiltration results. (B) Correlation heatmap of 22 immune cells. (C) Heat map of immune infiltration results of GSE166388. (D) Box plot of normal samples and PS samples. PS = psoriasis.

3.7. Validation of identified miRNAs and genes

Through the interaction network of miRNAs and mRNAs described above, we degraded 4 miRNAs and 6 genes in the core module. To verify the expression of these miRNAs and genes in the inflammatory cell model and in normal human immortalized keratinocytes (HaCaT), qRT-PCR experiments were performed after HaCaT cells were induced with 10 ng/mL TNF-α and IFN-γ (TI). The untreated group was the control group and the TI group was treated with 10 ng/mL TI. The results of qRT-PCR showed that the relative expression levels of IL-6 (P < .001) and IL-1β (P < .001) in the TI group were significantly higher than those in the control group, indicating that the inflammatory cell model was successfully constructed.

The expression levels of miRNAs, IL-6 and IL-1β in TI and control groups are shown in Figure 9, the expression levels of genes in TI and control groups are shown in Figure 10. In the inflammatory cell model, the relative expression level of hsa-miR-146a-5p (P = .0013) was higher than that of the control group. The relative expression of hsa-miR-18a-5p (P = .0267) was significantly lower than that of the control group. No statistically significant difference in hsa-miR-17-5p and hsa-miR-106a-5p between inflammatory cell models and healthy controls was observed. QRT-PCR results showed that the relative expression levels of CDK1 (P < .001), CCNB1 (P = .001) and RRM2 (P = .0341) were significantly higher in the inflammatory cell model than in the control group. The relative expression levels of CCNA2 (P = .001), CCNB2 ((P < .001) and MAD2L1 (P = .0043) were significantly lower than that of the control group. These results suggest that these DEMs and DEGs in psoriatic lesions might be involved in the pathogenesis of PS and differentially expressed in inflammatory diseases.

Figure 9. The relative expression of IL-6, IL-1β, miRNAs between the control group and TI group. miRNA = microRNA, TI = TNF-α and IFN-γ.

Figure 10. The relative expression of genes between the control group and TI group. TI = TNF-α and IFN-γ.

4. Discussion

PS is a common chronic relapsing inflammatory skin disease with complex pathogenesis. The origin and recurrence of the disease are not clear and are influenced by the interaction of several factors such as the immune system, genetic factors and environmental factors. Although there is currently no complete cure, the physical and psychological damage to patients can be minimized by identifying patients in the early stage of the disease and preventing secondary diseases through incremental lifestyle changes and personalized treatment methods.[2] The specific clinical detection index system for early and recurrence of the disease needs to be improved, so it is still necessary to explore the biomarkers and mechanisms of the disease. Gene expression can provide information about changes in the expression of mRNA transcripts and provide new insights into disease pathogenesis. Transcriptome analysis techniques such as microarrays and RNA-sequencing (RNA-seq), are relevant tools for the discovery of new biomarkers and therapeutic targets.[11] Through data mining and information integration, disease-related biomarkers can be screened effectively and correctly, which is beneficial for research into disease mechanisms.

miRNAs are small endogenous non-coding RNAs that have emerged as key regulators of gene expression and powerful biomarkers of disease. miRNAs can bind to the 3’-untranslated region (3’-UTR) of protein-coding mRNA in full or incomplete complementation to regulate degradation or translational repression.[12,13] MiRNAs have the potential to regulate keratinocyte proliferation, differentiation, apoptosis and inflammation, as well as the activation of subsets of T cell.[14] Evidence for the multiple roles of miRNAs in inflammatory skin diseases is rapidly accumulating.[15] Studies have confirmed that miR-21, miR-31, miR-146a, miR-155, and miR-203 are significantly upregulated in skin lesions in patients with PS.[16] MiRNAs are important regulators in the development of PS, and it is of interest to screen and explore miRNAs that could be potential biomarkers of PS.

Through a bioinformatics approach, we fully utilized published PS databases to screen miRNAs and mRNAs that could be used as potential biomarkers for PS. In our study, we used microarray datasets and the online bioinformatics tool GEO2R to screen the DEGs and DEMs between PS diseased skin tissue and normal skin tissue, and successfully obtained 205 slice DEGs from datasets GSE166388 and GSE153007, 6 section DEMs were screened from datasets GSE175438 and GSE145305. Then miRNet is used to predict the target genes, and the intersection of the target genes and cut DEGs is taken to obtain the dysregulated target genes for further analysis. Four key miRNAs and 6 hub genes in the regulatory pair were screened by constructing a miRNA-mRNA regulatory network. The qRT-PCR results showed that the selected miRNAs and mRNAs were statistically significant.

For the 10 miRNA-mRNA regulatory pairs composed of the selected core elements, the targeting relationship between miR-146a-5p and the predicted CCNA2 has been verified experimentally in existing literature.[17] Meanwhile, Ma X et al confirmed the targeting-binding relationship between miR-17-5p and RRM2 through the dual-luciferase reporter assay.[18] The relationship between other miRNAs and mRNAs was verified by crosslinking immunoprecipitation and high-throughput sequencing, photoactivatable ribonucleoside-enhanced cross-linking and immunoprecipitation, and other methods. More details are shown in Table 4.

Table 4 The details of the miRNA-mRNA core regulation pairs.

Verification of miRNA-mRNA regulatory pairs composed of core elements	
miRNA ID	Target	Experiment	PMID of literature	
hsa-miR-106a-5p	CCNB1	HITS-CLIP/PAR-CLIP	22100165/21572407/22012620/23024010/31101765/	
			23446348/23592263/23313552/24906430/	
			26061048/27292025/30455455/22927820/	
			22291592/23824327/24668909/24389009	
hsa-miR-146a-5p	CCNB1	NONE	NONE	
hsa-miR-17-5p	CCNB1	HITS-CLIP/PAR-CLIP	22100165/21572407/22012620/23024010/22291592/	
			23446348/23592263/23313552/24668909/24038734/	
			24389009/26061048/25768906/27150721/27292025/	
			29386283/30455455/22927820/31101765/34914716/	
			23824327/24906430/26701625/33406413	
hsa-miR-18a-5p	CCNB1	HITS-CLIP/PAR-CLIP	22100165/23024010/23592263/24906430/24389009/	
			26061048/27150721/26701625/30455455/22927820	
hsa-miR-146a-5p	CCNA2	qRT-PCR/Western blot/HITS-CLIP	19944095/23313552/26061048	
hsa-miR-18a-5p	CCNA2	Chimeric fragments/HITS-CLIP	24857550/26061048	
hsa-miR-17-5p	MAD2L1	Chimeric fragments/PAR-CLIP	22291592/24857550/27292025	
hsa-miR-18a-5p	RRM2	PAR-CLIP	22100165/21572407/22012620/22291592/23446348/	
			23592263/24668909/27292025/26701625	
hsa-miR-17-5p	RRM2	Dual-luciferase activity assay	37115505/22473208/22100165/21572407/22012620/	
		RNA-Seq	22291592/23446348/23592263/24668909/	
		HITS-CLIP/PAR-CLIP	27292025/26701625/30670076/35269913	
hsa-miR-106a-5p	RRM2	HITS-CLIP/PAR-CLIP	22100165/21572407/22012620/22291592/23446348/	
			23592263/24668909/26061048/27292025/22927820	

Cytokines have a wide range of biological activities that help coordinate the body response to infection and play important roles in the immune response, tissue repair, and inflammatory response. The occurrence of inflammation is often accompanied by the secretion of a large number of cytokines such as IL-6, IL-1β, IFN-γ, and the secretion of such cytokines is an important factor affecting the onset of inflammation in some diseases. Therefore, cytokines can be used as biomarkers to play a crucial role in the study of inflammatory diseases. Studies have shown that serum expression of cytokines such as IL-6, IFN-γ and TNF-α is significantly higher in patients with Ps than in healthy controls.[19–21] In our study, we used 10 ng/mL TI to induce HaCaT cells to construct an inflammatory cell model. We first demonstrated the relative expression levels of inflammatory cytokines IL-6 and IL-1β. The qRT-PCR results showed that the relative expression levels of IL-6 (P < .001) and IL-1β (P < .001) were significantly higher in the TI group than in the control group, indicating that the inflammatory cell model was constructed successfully.

QRT-PCR results showed that hsa-miR-146a-5p and hsa-miR-18a-5p were differentially expressed in inflammatory cell models and healthy controls. Furthermore, no statistically significant difference in hsa-miR-17-5p and hsa-miR-106a-5p was observed between inflammatory cell models and healthy controls. There was a significant difference in the high expression of hsa-miR-146a-5p (P = .0013) and the relative expression of hsa-miR-18a-5p (P = .0267) was significantly lower than that of the control group. The results showed that some miRNAs were differentially expressed in inflamed and normal HaCaT cells.

Sonkoly E et al reported that miR-146a-5p is overexpressed in the skin of patients with PS and may be involved in the regulation of the innate immune response and the tumor necrosis factor (TNF-α) pathway.[22] According to the ROC analysis in the literature, the combination of miR-146a-5p and miR-203a-3p levels in serum can discriminate PS patients from healthy controls more reliably than either miRNA alone, proving that multiple miRNAs play a synergistic role play pathogenesis of PS, and the combination of miRNAs may be more reliable in the diagnosis of PS.[23]

Studies have shown that the expression levels of TNF-α, IL-6, IL-8 and other pro-inflammatory factors in the intestinal epithelium of mice are upregulated after LPS stimulation, and the expression of miR-146a-5p is significantly upregulated. miR-146a-5p can promote the regeneration of intestinal epithelium by promoting cell proliferation.[24] In THP-1 cells, the expression of miR-146a could be upregulated by LPS stimulation.[25] The expression of miR-146a-5p and ox40l is significantly increased in the condylar chondrocyte inflammation model induced by IL-1β and TNF-α. miR-146a-5p can regulate T cell-mediated immunity by targeting ox40l in osteoarthritis.[26] Studies have reported the differential expression of small and medium RNAs in endothelial cells under inflammatory conditions, which confirmed that TNF-α can regulate the expression of small and medium RNAs in endothelial cells and the expression has certain tissue specificity.[27] The above reports indicated that the expression of miR-146a-5p was upregulated in an inflammatory state, which was consistent with the results in this paper. The pathogenesis of PS is related to the abnormal proliferation of keratinocytes. Under normal circumstances, keratinocyte proliferation and differentiation have a good balance to maintain the homeostasis of the epidermis, but this balance is broken in the inflammatory state. The abnormal up-regulation of miR-146a-5p expression in inflammatory states may be involved in the regulation of the biological functions of keratinocytes, resulting in abnormal proliferation and differentiation of keratinocytes, thus affecting the occurrence and development of PS.

Torri A et al found that miR-106, which is part of the in vitro CD4 T cell-derived miRNA signature, was consistently upregulated in serum from patients with PS compared to healthy controls and that miR-106a-5p was specific in Th1/Th17 cell-derived extracellular vesicles (EVs) in vitro. The expression of Treg-derived miR-146a-5p is increased in the serum of patients with PS.[28] Delic D et al stimulated keratinocytes with 10 ng/mL IFN-γ for 24 hours, and expression of hsa-miR-106a-5p and miR-17-5p increased, which differs from our results. The reason may lie in the difference caused by the action of a single inflammatory cytokine and the combined action of multiple inflammatory cytokines.[29]

In our study, the qRT-PCR results showed that the relative expression levels of CDK1 (P < .001), CCNB1 (P = .001), and RRM2 (P = .0341) were significantly higher in the inflammatory cell model than in the control group. The relative expression levels of CCNA2 (P = .001), CCNB2 (P < .001), and MAD2L1 (P = .0043) were significantly lower than that of the control group.

In the skin lesions of PS patients, CDK1 expression is widespread in all layers of the epidermis.[30] CCNB1 is expressed at low levels in the nucleus and cytoplasm of healthy skin, mainly in the basal cell layer of the epidermis and some hair follicles, while CCNB1 expression is significantly higher in psoriatic lesions and all layers of the epidermis than in healthy skin.[30]

Related studies have shown that the cell cycle of PS patients is significantly shortened and the high expression of CDK1 and CCNB1 promotes the proliferation of keratinocytes in PS.[31] Meanwhile, Wilms tumor 1-binding protein promotes keratinocyte proliferation by regulating CCNA2 and CDK2, thereby promoting the development of PS.[32] In the study by Li AH et al, it was confirmed that the expression of CCNB1 and CCNB2 in the serum of patients with PS was higher than that of healthy controls. CCNB1 and CCNB2 may not only regulate the cell cycle of keratinocytes through a variety of pathways, but also regulate immune cells such as macrophages and mast cells, support the release of important signaling molecules, and are involved in and promote the pathogenesis and progression of PS.[33] These results suggest that abnormally high expressions of CDK1 and CCNB1 may lead to changes in the cell cycle process, significantly shortening the cell cycle, promoting the abnormal proliferation of keratinocytes, and participating in the regulation of immune cells such as macrophages and mast cells. It is involved in the pathophysiological process of PS by influencing the biological functions of keratinocytes and immune cells. The abnormally high expression of CDK1 and CCNB1 in this study was consistent with the literature.

Interestingly, in our study, we used psoriatic lesions for bioinformatic analysis, and an in vitro inflammatory cell model was used to verify the relative expression levels, which showed consistent and divergent expression results. In HaCaT cells stimulated with 10 ng/mL TI, expression of miR-18a-5p, CCNA2, CCNB2, and MAD2L1 was downregulated, contradicting the results of the bioinformatic analysis. This indicates that some specific miRNAs and genes have different expression activities in different inflammatory environments and in different cells and tissues.

The mRNAs regulated by DEMs in psoriatic skin can elicit various functional effects such as skin barrier/keratinocyte function, lipid metabolism, cell signaling/cell cycle and immune response mechanism. The relationship between miRNAs and mRNA is often based on predictive models, but the expression of genes in the human body is regulated by a variety of factors.

This study has several limitations. We have validated only the most important miRNAs and genes in cell models, and more detailed studies are needed to elucidate the specific mechanisms and roles of miRNA-mRNA regulation in PS. In the future, we plan to fully utilize clinical tissue samples and blood samples for cross-validation to more comprehensively explore the related mechanisms of miRNA and genes in PS.

5. Conclusion

Our study attempted to identify key regulatory miRNA-mRNA pairs, significantly DEMs and genes in PS lesion skin through bioinformatic analysis. This study showed that expression of hsa-miR-146a-5p, CDK1, CCNB1 and RRM2 was upregulated under inflammatory conditions. However, the expression of hsa-miR-18a-5p, CCNA2, CCNB2 and MAD2L1 was downregulated in an inflammatory state. These DEMs and DEGs in skin lesions of patients with PS may be involved in PS pathogenesis, suggesting that they have potential biomarker value in PS and provide information for PS diagnosis and treatment. Their role in PS could provide new clues to further investigate the relevant molecular mechanisms of PS in the future.

Author contributions

Conceptualization: Wenhong Qiu, Kaiwen Guo.

Data curation: Lu Chen, Xiaochen Wang, Chang Liu.

Formal analysis: Lu Chen, Xiaoqing Chen, Peng Li.

Funding acquisition: Wenhong Qiu.

Software: Xiaochen Wang, Chang Liu.

Validation: Lu Chen.

Writing – original draft: Lu Chen, Wenhong Qiu.

Writing – review & editing: Wenhong Qiu, Kaiwen Guo.

Abbreviations:

BP biological process

CC cellular component

DEGs differentially expressed genes

DEMs differentially expressed miRNAs

GEO gene expression omnibus

GO gene ontology

KEGG Kyoto encyclopedia of genes and genomes

MF molecular function

miRNA microRNA

PPI protein-protein interaction

PS psoriasis

qRT-PCR quantitative reverse transcription polymerase chain reaction

TI TNF-α and IFN-γ

This study was funded by National Natural Science Foundation of China (Grant number 31671092), Cooperative Education between Industry and Education (Construction of New Engineering, New Medical, New Agricultural and New Liberal Arts) of the Department of Higher Education of the Ministry of Education (Grant number 202102585001) and Hubei Education Department (Grant number 2020435). The research performed in this study was in compliance with the laws of China and the authors’ respective institutions.

This is a study involving only bioinformatics. The Medical ethics committee of Jianghan University has confirmed that no ethical approval is required.

The authors have no conflicts of interest to disclose.

The datasets generated during and/or analyzed during the current study are publicly available.

How to cite this article: Chen L, Wang X, Liu C, Chen X, Li P, Qiu W, Guo K. Integrative analysis of gene and microRNA expression profiles reveals candidate biomarkers and regulatory networks in psoriasis. Medicine 2024;103:29(e39002).
==== Refs
References

[1] Lowes MA Bowcock AM Krueger JG . Pathogenesis and therapy of psoriasis. Nature. 2007;445 :866–73.17314973
[2] Griffiths CEM Armstrong AW Gudjonsson JE . Psoriasis. Lancet (London, England). 2021;397 :1301–15.33812489
[3] O’Brien J Hayder H Zayed Y . Overview of MicroRNA Biogenesis, mechanisms of actions, and circulation. Front Endocrinol. 2018;9 :402.
[4] Yang SC Alalaiwe A Lin ZC . Anti-Inflammatory microRNAs for treating inflammatory skin diseases. Biomolecules. 2022;12 :1072.36008966
[5] Huang da W Sherman BT Lempicki RA . Bioinformatics enrichment tools: paths toward the comprehensive functional analysis of large gene lists. Nucleic Acids Res. 2009;37 :1–13.19033363
[6] Barrett T Wilhite SE Ledoux P . NCBI GEO: archive for functional genomics data sets--update. Nucleic Acids Res. 2013;41 :D991–995.23193258
[7] Krishnan VS Kõks S . Transcriptional basis of psoriasis from large scale gene expression studies: the importance of moving towards a precision medicine approach. Int J Mol Sci . 2022;23 :6130.35682804
[8] Kanehisa M Goto S . KEGG: Kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28 :27–30.10592173
[9] Kanehisa M . Toward understanding the origin and evolution of cellular organisms. Protein Sci. 2019;28 :1947–51.31441146
[10] Kanehisa M Furumichi M Sato Y . KEGG for taxonomy-based analysis of pathways and genomes. Nucleic Acids Res. 2023;51 :D587–92.36300620
[11] Rioux G Ridha Z Simard M . Transcriptome profiling analyses in psoriasis: a dynamic contribution of keratinocytes to the pathogenesis. Genes. 2020;11 :1155.33007857
[12] Bartel DP . MicroRNAs: genomics, biogenesis, mechanism, and function. Cell. 2004;116 :281–97.14744438
[13] Guo H Ingolia NT Weissman JS . Mammalian microRNAs predominantly act to decrease target mRNA levels. Nature. 2010;466 :835–40.20703300
[14] Wu R Zeng J Yuan J . MicroRNA-210 overexpression promotes psoriasis-like inflammation by inducing Th1 and Th17 cell differentiation. J Clin Invest. 2018;128 :2551–68.29757188
[15] Jinnin M . Various applications of microRNAs in skin diseases. J Dermatol Sci. 2014;74 :3–8.24530178
[16] Raaby L Langkilde A Kjellerup RB . Changes in mRNA expression precede changes in microRNA expression in lesional psoriatic skin during treatment with adalimumab. Br J Dermatol. 2015;173 :436–47.25662483
[17] Hsieh CH Rau CS Jeng SF . Identification of the potential target genes of microRNA-146a induced by PMA treatment in human microvascular endothelial cells. Exp Cell Res. 2010;316 :1119–26.19944095
[18] Ma X Fu T Ke ZY . MiR-17- 5p/RRM2 regulated gemcitabine resistance in lung cancer A549 cells. Cell cycle (Georgetown, Tex.). 2023;22 :1367–79.37115505
[19] Arican O Aral M Sasmaz S . Serum levels of TNF-alpha, IFN-gamma, IL-6, IL-8, IL-12, IL-17, and IL-18 in patients with active psoriasis and correlation with disease severity. Mediators Inflamm. 2005;2005 :273–9.16258194
[20] Michalak-Stoma A Bartosińska J Raczkiewicz D . Multiple cytokine analysis of Th1/Th2/Th9/Th17/Th22/Treg cytokine pathway for individual immune profile assessment in patients with psoriasis. Med Sci Monitor. 2022;28 :e938277.
[21] Andersen CSB Kvist-Hansen A Siewertsen M . Blood cell biomarkers of inflammation and cytokine levels as predictors of response to biologics in patients with psoriasis. Int J Mol Sci . 2023;24 :6111.37047086
[22] Sonkoly E Wei T Janson PC . MicroRNAs: novel regulators involved in the pathogenesis of psoriasis? PLoS One. 2007;2 :e610.17622355
[23] Jinnin M . Recent progress in studies of miRNA and skin diseases. J Dermatol. 2015;42 :551–8.25917002
[24] Chen X Li W Chen T . miR-146a-5p promotes epithelium regeneration against LPS-induced inflammatory injury via targeting TAB1/TAK1/NF-κB signaling pathway. Int J Biol Macromol. 2022;221 :1031–40.36096257
[25] Nahid MA Pauley KM Satoh M . miR-146a is critical for endotoxin-induced tolerance: implication in innate immunity. J Biol Chem. 2009;284 :34590–9.19840932
[26] Yu D Wei W Hefeng Y . Upregulated ox40l Can Be Inhibited by miR-146a-5p in Condylar Chondrocytes Induced by IL-1β and TNF-α: a possible regulatory mechanism in osteoarthritis. Int Arch Allergy Immunol. 2021;182 :408–16.33147588
[27] Liu P Hu L Shi Y . Changes in the Small RNA expression in endothelial cells in response to inflammatory stimulation. Oxid Med Cell Longevity. 2021;2021 :8845520.
[28] Torri A Carpi D Bulgheroni E . Extracellular microRNA signature of human helper T cell subsets in health and autoimmunity. J Biol Chem. 2017;292 :2903–15.28077577
[29] Delić D Wolk K Schmid R . Integrated microRNA/mRNA expression profiling of the skin of psoriasis patients. J Dermatol Sci. 2020;97 :9–20.31843230
[30] Xue X Yu J Li C . Full-length transcriptome sequencing analysis of differentially expressed genes and pathways after treatment of psoriasis with oxymatrine. Front Pharmacol. 2022;13 :889493.35721124
[31] Ni X Lai Y . Keratinocyte: a trigger or an executor of psoriasis? J Leukoc Biol. 2020;108 :485–91.32170886
[32] Kong Y Wu R Zhang S . Wilms’ tumor 1-associating protein contributes to psoriasis by promoting keratinocytes proliferation via regulating cyclinA2 and CDK2. Int Immunopharmacol. 2020;88 :106918.32866786
[33] Li AH Chen YQ Chen YQ . CCNB1 and CCNB2 involvement in the pathogenesis of psoriasis: a bioinformatics study. J Int Med Res. 2022;50 :3000605221117138.35949173
