Screening and Identication of Potential Key Biomarkers in Psoriasis : Evidence From Bioinformatic Analysis

Background: Psoriasis is a chronic and prevalent skin condition brought on by various genetic and external factors. Till date, there is no cure for this disease. It is, therefore, crucial to examine the underlying mechanisms leading up to psoriasis. The goal of our study was to explore the mechanistic pathways involved in the molecular pathogenesis of psoriasis. Methods: Using Gene Expression Omnibus (GEO), we performed an extensive analysis of the transcript expression prole in psoriasis patients. The datasets GSE13355, GSE30999, and GSE106992, containing 239 pairs of normal and psoriatic skin samples were arbitrarily assigned to two non-overlapping cohorts for cross-validated differential gene expression analysis. The gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment and gene set enrichment analysis (GSEA) were employed for interpretation, visualization, and unied recognition. The STRING database was used to construct the protein–protein interaction (PPI) network and the hub genes were established using Cytoscape. Additionally, both differentially regulated genes and functional likeness of hub genes were further examined in terms of gene activation and functional outcome. The correlation between normal tissue and inltrating immune cells was analyzed by CIBERSORT. Moreover, ROC analysis was performed to distinguish between skin lesion samples and skin non-lesion samples. In addition, a signaling axis involving lnRNA, miRNA, mRNA, and ceRNA was generated with information from the DEmiRNA, DElncRNA, and DEmiRNA-DEmRNA relationship. Lastly, immunohistochemical evaluations were used to analyze the highest expression of single gene in the whole body within the Human Protein Atlas (HPA) database. Results: The genetic proles of 239 pairs of normal and lesional skin samples were downloaded from three datasets in the GEO database. PPI network revealed a tight interaction among 197 differentially expressed genes (DEGs). Moreover, gene ontology analyses indicated that psoriasis-related DEGs mostly included viral defense genes, type I interferon axis, and its corresponding cellular responses. The Kyoto encyclopedia of genes and the DEGs enrichment analysis showed involvement of the NOD-like receptor signaling network, cytokine-cytokine


Introduction
Psoriasis is a chronic autoimmune skin disorder manifested by scaly erythema and regulated by genetic and environmental agents [1]. Although psoriasis poses no threat to life, it can present serious mental and economic burden to patients. Fortunately, with better understanding of psoriasis [2][3], an increasing number of therapies have emerged. The mainstream psoriatic therapy includes systemic application of retinoic acid, immunosuppressant, cyclosporine, biological agents, external use of glucocorticoid, retinoic acid drugs, and so on. However, the adverse reactions are serious. Due to the unclear etiology and complex pathogenesis of psoriasis, it still remains incurable. Hence, there is great urgency in the comprehension of molecular networks involved in psoriasis to enhance the quality of psoriatic care.
Psoriasis is mediated by polygenic inheritance T cells. The excessive activation of immune cells is caused by a variety of pathogenic factors. At present, multiple new drugs are available for the management of psoriasis. The biological agents used for psoriasis therapy in China are anti-TNF-α, anti-IL17A, anti-IL12/23, and anti-IL23. Anti-TNF-αwas the rst commercially available biological agent for psoriasis, and it includes Etanercept, Adalimumab, and In iximab. Among them, Etanercept sequesters serum TNF-α, which limits its activity, and produces an anti-in ammatory response [4]. Likewise, Adalimumab also binds to TNF -α and blocks its interaction with the TNF receptor on the surface of p55 and p75 cells, thereby reducing in ammation. Lastly, In iximab is a human mouse chimeric monoclonal antibody, which strongly sequesters both soluble and transmembrane forms of TNF -α, thereby making them inactive [5]. According to reports, Adalimumab and In iximab fared better than Etanercept in controlling psoriatic skin lesions. Moreover, the onset speed of Adalimumab and In iximab remained slower than that of anti-IL17A [6].The anti-IL17A drugs prevalent in China are Scuchizumab and Iqizumab. Both drugs have similar e cacy in psoriatic lesion resolution and have a better rate of improvement in PASI than monoclonal antibodies [7][8][9]. Despite the presence of multiple psoriatic drugs, there is still potential for advancements, especially in terms of making the drugs affordable and in reducing double infection risk. Hence, understanding the intricate details of psoriatic pathogenesis will add value to the future of psoriatic care.
Multiple studies have revealed that a variety of chemokines and receptors mediate pathogenesis of psoriasis [10][11], and the CXC family plays an important role. Some chemokines stimulate the chemotactic induction of in ammatory T cells (such as CXCL10, CXCL16) and neutrophils (such as CXCL1 and CXCL5) during psoriasis. The chemokine receptors involved during this process include: T cells: cxcri, CXCR2, CXCL3, CXCL10, and CXCL16; and neutrophils: CXCL1, CXCL5, CXCR3, and CXCR6 [12]. There is emerging evidence that CXCL10 and its receptor CXCR3 are expressed in keratinocytes and dermis of plaque psoriasis patients. Moreover, upon effective treatment, CXCL10 levels diminish, whereas CXCL10 levels remain high in the serum of early psoriatic patients, but reduce in serum of patient with longer duration of the disease [13]. Multiple studies have revealed that the in ltrating leukocytes in psoriatic lesions strongly express CCR4 and CXCR3, whereas, in the dermis, the ligands CCL17, CCL22, and CXCL9, CXCL10, and CXCL8 are highly expressed. Alternately, others have shown that the CCL2-and CXCL10-related antibodies diminish the expression of these ligands in psoriatic lesions, thereby reducing severity of psoriasis, as evidenced by the PASI score [14]. CXCL8 is known to regulate neutrophil in ltration, angiogenesis, and keratinocyte proliferation in psoriatic lesions. It is worth mentioning that serum CXCL8 in psoriatic patients remains signi cantly elevated and is positively correlated with the severity of the disease [15][16]. Nevertheless, the speci c function of CXCL8 in psoriatic progression remains unknown and requires further investigation.
Over the past decade, advances in microarray screening and bioinformatics investigations have enabled the genome-wide search for differentially expressed genes (DEGs) and physiological networks linked to speci c diseases, including psoriasis. Yet, there have been limitations. For example, a high false-positive incidence was noted in an isolated microarray study. Additionally, differences in microarray procedures, variations in statistical analysis, and limited sample populations may have, in unison, introduced discrepancies in conclusions across studies. Hence, additional extensive examinations are warranted to reliably recognize psoriasis-related DEGs in the search of effective targeted therapy and prognostic biomarkers for psoriasis.
The goal of our work was to analyze DEGs, collected from three microarray datasets from the Gene Expression Omnibus (GEO) database, comparing psoriatic skin to healthy skin. Additionally, we conducted network analyses to further evaluate underlying moledular mechanisms involved in psoriasis.
Using the results of our analysis, we established an extensive network of genes and functional pathways involved in psoriasis. Moreover, we demonstrated CXCL8 to be a promising psoriatic biomarker.

Screening of DEGs
DEG analysis was done with R 3.4.1's (https://www.r-project.org/) Affy and Limma package [20][21]. The corresponding plots were made by R and the genes that followed the following criteria were regarded as DEGs: P 0.05 and log |FC| 1.0. Finally, we generated Venn charts [22], including an overlap of total, upregulated, and downregulated DEGs from all 3 datasets.

GO and KEGG enrichment analyses of DEGs
Gene Ontology (GO) [23] enrichment analysis and Kyoto encyclopedia of Genes and Genomes (KEGG) [24] functional network analysis were conducted using "clusterPro ler" [25] within R. This allowed us to recognize potential DEG roles and relationships with corresponding diseases. Gene sets carrying a p value 0.05 was regarded as signi cant enrichment.

GSEA analysis
The relationship between CXCL8 and functional network genes was analyzed with gene set enrichment analysis (GSEA, v2.2,http://www.broad.mit.edu/gsea/) [26]. Furthermore, its association with phenotype was explored with TCGA gene sets that included 33 forms of cancer from MSigDB. The gene sets that showed marked enrichment by genes, based on elevated CXCL8levels, [false discovery rate (FDR) 0.05] were de ned as enriched gene sets.

Integration of the protein-protein interaction (PPI) network
The PPI network was estimated with the online Search Tool for the Retrieval of Interacting Genes (STRING; http://string-db.org)(version 10.0) [27] database and it was formed with the STRING database, with a combined interaction score 0.4 deemed as signi cant. Cytoscape (version 3.4.0) was employed for the visualization of potential underlying signaling pathways [28]. Finally, Ggplot2 package [29] was used to generate a bar chart illustrating distribution of relevant genes.

Immune in ltration analysis
CIBERSORT was employed for immune cell in ltration analysis [30] and facilitated the identi cation of alterations in in ltrating immune cell between psoriatic and healthy skin specimen. The corresponding bubble graph was visualized with ggplot2 package.
Receiver operating characteristic (ROC) curve analysis ROC curve analysis, conducted with the "pROC" package of R [31], was chosen for determining speci city and sensitivity of the diagnosis [32]. The maximum value of the sum of speci city and sensitivity was set as the threshold for individual hub genes.

Construction of the ceRNA network
The relationship between DElncRNAs and DEmiRNAs were estimated with miRcode [33]. The DEmiRNAtargeted mRNAs were obtained based on starBase [34], miRDB [35], and miRWalk[36]. To be deemed as candidate mRNAs, the mRNAs had to be identi ed by all data sets and required intersection with DEmRNAs to remove DEmiRNAs-targeted DEmRNAs. Next the lncRNA-miRNA-mRNA ceRNA axis was generated using the DEmiRNA-DElncRNA and DEmiRNA-DEmRNA relationships, and was visualized with the Cytoscape software [28]. Please refer to Figure 11A for a detailed depiction of the lncRNA-miRNA-mRNA ceRNA network analysis. Potential associations between ceRNA levels were analyzed with linear regression. Immunohistochemical (HIC) analysis of lymph nodes through HPA database To explore the diagnostic signi cance of CXCL8, whole-body distribution of this gene (both at the transcript and protein level) was examined using the Human Protein Atlas (HPA, https://www.proteinatlas.org/), an open access website providing protein maps within organs, tissues, and cells [37].

Analysis of DEGs
We utilized the Affy and Limma packages to analyze the three GEO data sets. The Venn diagram, consisting of 197 DEGs, from the three data sets, is illustrated in Figure 1A. The DEGs were further analyzed using the STRING database to identify potential links with other DEGs. An interaction network diagram for these genes was additionally generated ( Fig. 2A). Figure 2B illustrates each DEG as a node and presents the degrees of each node. Based on this analysis, the degrees of the CXCL8 gene were found to be the largest, suggesting an essential role of CXCL8 in psoriasis. Hence, CXCL8 was selected for further analysis.
Analyzing single gene differences Next, we separated all samples into a high-or low-expressing group, based on their level of CXCL8 expression, and analyzed the differences between the groups (Figure 3). The GSE13355 data set had 62 high-expressing genes and 62 low-expressing genes. Likewise, the GSE30999 data set had 82 highexpressing genes and 82 low-expressing genes. Lastly, the GSE106992 data set had 96 high-and 96 lowexpressing genes. Moreover, the heat map showed a marked difference in CXCL8 levels between the highand low-expressing groups. Finally, based on the volcano plot, most genes were down-regulated in the GSE13355 and GSE30999 data sets, while the remaining genes were up-regulated in GSE106992.

CXCL8 enrichment analysis
The Venn diagram of 254 low-and high-expressing co-DEGs from the three data sets is shown in Figure  2B. In addition, we performed GO enrichment analysis and KEGG enrichment analysis on these genes ( Figure 4). The GO analysis showed that these genes were markedly enriched in the "defense response to virus", "type I interferon signaling pathway", and "cell response to type I interferon" categories. Moreover, based on the KEGG analysis, these genes are markedly enriched in the "NOD-like receptor signaling pathway", "interaction between cytokine and cytokine receptor", and "IL−17 signaling pathway" categories (Table 1).

GSEA and analysis of pan-cancer pathway
The biological pathways that were signi cantly altered in the CXCL8 high-expressing verses lowexpressing group were determined using GSEA. These pathways were mainly related to the complement-like network ( Figure 5A). To examine CXCL8 pathway enrichment in each cancer type, pan-cancer pathway analysis of CXCL8 was performed ( Figure 5B) to identify the distribution pattern of CXCL8 gene in pan-cancer.

Immune in ltration analysis
The immune in ltration analyses of all three data sets were performed separately and the immune cell composition of each sample are displayed in a bar chart. We also analyzed correlation between immune cell expression and differences within them, which are displayed in a heat map and volcano plot respectively. The analytical conclusions from data sets GSE13355, GSE30999, and GSE106992 are provided in Figures 6, 7, and 8 respectively.

ROC analysis
To better elucidate the relationship between CXCL8 expression and classi cation of diseased verses healthy tissues, we performed ROC analysis of the CXCL8 gene among the 3 data sets (Figure 9). The area under the ROC curves was 0.941 for GSE13355, as shown in Fig. 9A. The area under the ROC curve was 0.935 for GSE30999, as shown in Fig. 9B. The area under the ROC curve was 0.794 for GSE106992, as shown in Fig. 9C.

HIC analysis
The CXCL8 gene was further analyzed, using the HPA database, to obtain CXCL8 gene expression and distribution patterns throughout the human body. We discovered that CXCL8 was most highly expressed within the lymph node. We, thus, selected the lymph node for HIC analysis, and downloaded the HIC map of the CXCL8 gene in the lymph node ( Figure 11).

Discussion
Psoriasis is a highly prevalent and chronic in ammatory skin disorder. According to the latest epidemiological survey, the annual incidence rate of psoriasis in European countries is about 0.96 per thousand [38]. This disease is often chronic and recurrent in nature, which offers different degrees of negative impact on the life of patients, involving physiology, psychology, daily life, social contact, family and other aspects. A foreign survey on psoriatic patients revealed that 79% of the patients agreed that the disease had a negative impact on their lives, and 10% of the patients entertained the idea of suicide [39].
Unfortunately, despite the availability of multiple therapeutics for psoriasis, this disease remains incurable and a vast majority of treated patients experience adverse reactions to psoriatic drugs. As a result, there is great urgency in the development of highly speci c and e cient psoriatic biomarkers. Here, we utilized bioinformatics analysis of three separate microarray data sets, involving a total of 265 psoriatic vs. 213 healthy samples, to identify genes speci c to the pathogenesis of psoriasis. Among the DEGs studied, CXCL8 showed the largest fold change and is intimately linked to viral defense, type I interferon axis, and its corresponding cellular response. We further delineated the CXCL8 role in psoriasis by assessing its overall expression pro le along with evaluating its underlying mechanism of action.
Multiple studies have demonstrated CXCL8 to be an essential cytokine. It often regulates tumor growth, invasion, and neovascularization and can be found in excess in tissues from skin squamous cell carcinoma patients [40]. Earlier reports have also indicated that CXCL8 modulates the formation of psoriasis Munro micro abscess. We, therefore, assessed CXCL8 expression in patients with psoriasis. Given that the role of CXCL8 in psoriasis has not been fully elucidated, we also explored the underlying function of CXCL8 in psoriasis.
Here, we conducted a comprehensive analysis of three microarray data sets to retrieve DEGs involved with psoriasis verses healthy skin specimen. We discovered 187 DEGs in total among the 3 data sets. We further explored interrelationships between DEGs, using GO and KEGG enrichment analyses. Moreover, we conducted analysis with the STRING database to verify the links among the DEGs. 30 DEGs were chosen as hub genes carrying degrees ≥20. Among them, CXCL8 exhibited the largest node degrees with 76, indicating its strong importance in the development and/or progression of psoriasis. Furthermore, using GO enrichment analysis, we demonstrated that alterations in the markedly obvious modules occurred in response to virus, type I interferon axis and its corresponding cellular response. Conversely, using KEGG, it was shown that the DEGs were mainly in the NOD-like receptor axis, interaction between cytokine and cytokine receptor and the IL−17 axis. Multiple studies have reported a strong relation between the IL−17 axis and psoriasis. In fact, IL−17 expression is markedly elevated in the psoriatic skin, compared to healthy skin [41]. It was also shown in keratinocytes that IL-17 stimulated the production of proin ammatory cytokines like IL-6 and IL-8, which exacerbated psoriasis [42]. Moreover, topical administration of imiquimod, a Toll-like receptor (TLR) 7/8 ligand and powerful immunologic stimulator, resulted in psoriasis-resembling dermatitis in mice, with simultaneous increases in IL-17A and IL-17F levels [43]. Alternately, both CsA and anti-TNF-a agents downregulated IL-17A, IL-23p19, and chemokine (C-C motif) ligand 20 in psoriatic skin, thereby improving eruptions [42.44-46]. Collectively, these evidences suggest a strong involvement of IL-17A in the progression of psoriasis.
Interleukin-22 (IL-22) belongs to the IL-10 cytokine family [47][48] and is intimately linked to psoriatic in ammation (Fig. 1). It is primarily made by Th17 cells. However, T22 have also been shown to exclusively produce IL-22. The IL-22 receptor can be found on keratinocytes, wherein IL-22 stimulates epidermal hyperplasia via a rise in keratinocyte proliferation. Clinically, serum IL-22 is markedly elevated in psoriatic patients, as opposed to healthy individuals [49]. Additionally, IL-22 is known to be highly expressed in psoriatic lesions and the levels fall with anti-psoriatic therapy [39]. Based on these data, IL-22 not only plays a crucial role in keratinocyte proliferation, but also regulates pasoriasis development and progression.
Dendritic cells produce interleukin (IL)-23, which regulates Th17 proliferation. The TH17 cells, in turn, are responsible for stimulating a variety of cytokines, including IL-17A, IL-17F, and IL-22. Among them, IL-17A and IL-22 was reported to be involved in keratinocyte proliferation and positively regulate tumor necrosis factor (TNF)-a, chemokine (C-X-C motif) ligand (CXCL)1, and CXCL8. Taken together, these conclusions suggest the involvement of these genes in multiple mmune-provoked disease states.
Using single gene ROC analysis, we revealed that CXCL8 has great potential to serve as a bio-marker for psoriatic diagnosis (AUC values 0.7). Moreover, based on CXCL8 predictions from starBase, miRDB and miRWalk, we generated a CXCL8-ceRNA network that included hsa-miR-1294, hsa-miR-140-3p, hsa-miR-185-5p, hsa-miR-4306, hsa-miR-4644, and hsa-miR-493-5p. lncRNAs modulate immune and in amma tory pathways via multi-mechanistic regulation of gene expression [50]. Among these mechanisms, the ceRNA axis has been widely investigated. Based on the ceRNA theory, lncRNA activates target gene expression by sequestering miRNA that would otherwise suppress gene activation. [51]. Qiao et. al. [52], for instance, revealed a lncRNA-MSX2P1 axis that stimulates S100A71 and facilitates IL-22-induced keratinocytes production via depressing miR-6731-5p action. Similarly, Li et al [53] reported that lncRNA H19 modulates keratinocyte differentiation by sequesteting miR-130b-3p and augmenting desmoglein 1 expression. Consistent with these studies, our work also demonstrated a strong involvement of ceRNA in the development and progression of psoriasis via sequestering multiple potential miRNAs.
This study highlighted a promising bio-marker for psoriasis. Nevertheless, our work had some notable limitations. Firstly, the strati cation analysis and interaction between factors did not receive detailed investigation. Moreover, due to its retrospective nature, we also couldn't rule out selective and recall biases. Hence, we propose prospective future studies to increase the statistical power and obtain more reliable conclusions. Secondly, apart from bioinformatics being a powerful tool for genome-wide studies, inadequate knowledge of underlying mechanisms may have limited our identi cation of effective psoriatic bio-markers. Hence, CXCL8 must be examined extensively at the molecular, cellular, and organismal level to con rm its signi cance in psoriatic diagnosis.

Conclusions
In summary, this study analyzed 265 psoriatic skin lesion samples to identify DEGs responsible for psoriasis progression. We revealed 197 total DEGs associated with psoriasis and we recognized CXCL8 as a key bio-marker for psoriasis. Further investigations are warranted to explore all aspects of CXCL8 function in psoriasis.