Identication of hub genes and pathways in psoriasis through bioinformatics and validation by RT-qPCR

Although several studies have attempted to investigate the aetiology and mechanism of psoriasis, the precise molecular mechanism remains unclear. Our study aimed to identify the hub genes and associated pathways that promote its in which interaction and enriched These hub genes were validated in the four aforementioned datasets and M5-induced HaCaT cells using real-time quantitative polymerase chain reaction (RT-qPCR).

In the current study, we merged four independent gene expression datasets (GSE30999 [6], GSE34248 [7], GSE41662 [7] and GSE50790 [8]) from the GEO database using the affy package in R software, and identi ed differentially expressed genes (DEGs) of the merged datasets. Then the DEGs were performed for a series of analyses, such as Gene ontology (GO) enrichment, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis, protein-protein interaction (PPI) network analysis, and hub genes prediction.
In order to verify these predicted hub genes, we compared the expression levels of these genes between the lesional and non-lesional psoriatic skin in the four datasets, and performed RT-qPCR experiments in M5 induced HaCaT cells. Since some studies con rmed M5 (IL-1a, IL-17A, IL-22, TNF-a, and oncostatin M) can induce a better psoriasis cell model than other cytokines, we use M5 to mimic the hyperproliferative pathological state of psoriasis keratinocytes in our experiments [9][10][11][12][13][14] Methods Microarray data source Fig. 1 shows the overall work ow of this study. Using the key words 'psoriasis', 'lesion*', 'Tissue', 'Homo sapiens' and 'Expression pro ling by array' were screened in the National Center of Biotechnology Information (NCBI) -GEO (https://www.ncbi.nlm.nih.gov/geo/). The raw pro les of GSE30999 [6], GSE34248 [7], GSE41662 [7], and GSE50790 [8], were downloaded from the NCBI GEO database. The detailed information of four enrolled microarray datasets is shown in Table 1. In GSE30999, GSE34248, GSE41662, and GSE50790, only psoriatic lesional skin samples and their matched non-lesional skin samples were selected, which resulted in 127 pairs of skin samples from psoriasis patients for subsequent analysis.

DEGs analysis and integration
After the raw data of four datasets were collated, the analyses were performed in the R programming environment (version4.0.2.). The affy package was used to preprocess the collated data by conducting background correction, normalization and expression calculation [15]. Then the ComBat function of sva package was applied for eliminating batch effect between datasets [16]. And the probes were annotated by matrix data combined with HG-U133_Plus_2-na36 annotation le. Lastly we used the limma package of R language to identify the differential expression genes(DEGs). The signi cant differential expression cut-off was set as |logFC| >1.5 and adjusted P<0.05, which were demonstrated as Volcano plots. Furthermore, a heatmap was drawn to show the clustering of samples.
GO and KEGG pathway analysis of DEGs Metascape (http://metascape.org) is a free, powerful and online analysis tool for gene annotation and analysis [17]. It is also an automated meta-analysis tool to annotate genes, conduct enrichment analysis and construct protein-protein interaction networks. In this study, Gene ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were enriched based on Metascape. GO terms included biological process, cellular component, and molecular function. The cutoff criterions were P-value < 0.05, minimum overlap of 3, and minimum enrichment of 1.5. The top 20 enrichment Go terms and KEGG pathways in our study were listed as gures and tables respectively.

Protein-protein interaction(PPI) network construction and hub genes identi cation
The Search Tool for the Retrieval of Interacting Genes/Proteins (STRING, version 11.0, http://www.stringdb.org/) was utilized to construct a PPI network for proteins encoded by DEGs, and the minimum required interaction score was set as 0.7 (referred as high con dence) [18]. In order to further analyze the results from STRING database, Cytoscape 3.8.0 was used as a tool to visualize the PPI network [19]. The modules in the PPI network were extracted using the Cluster Finding algorithm in MCODE with a max depth of 100, a node score cutoff of 0.2, a degree cutoff of 2, and a k-core of 2. The hub genes were identi ed by CytoHubba, which was a useful plugin in Cytoscape. It can provide 12 topological analysis methods to explore top 10 genes. Speci cally, these methods are Maximal Clique Centrality (MCC), Density of Maximum Neighborhood Component (DMNC), Maximum Neighborhood Component (MNC), Degree, Edge Percolated Component (EPC), BottleNeck, EcCentricity, Closeness, Radiality, Betweenness, Stress, and ClusteringCoe cient. A hub gene can be identi ed, if this gene was predicted as one of the top 10 genes in all the 12 methods. Additionally, interactions between hub genes were analyzed back in STRING database. Functional enrichment analysis of hub genes GO and KEGG analyses of 359 DEGs were performed by Metascape, including seven hub genes. Then the functional enrichment analysis of seven hub genes were explored by http: //www. bioinformatics. com.cn, an online platform for data analysis and visualization.

Validation of hub genes in datasets
Expression data of hub genes were extracted from GSE30999, GSE34248, GSE41662, and GSE50790. The correlations between CXCL10 and the other hub genes were con rmed in the four aforementioned datasets. Then the differential expression of these hub genes was validated in these four datasets.  Table 2.

Statistical analysis
The data extracted from GEO datasets were examined by the normality test and homogeneity of variance test. Spearman's correlation analysis was utilized to investigate the correlation between CXCL10 and the other hub genes. T-testing and analysis of variance (ANOVA) testing were used for comparisons between two groups and among three or more groups respectively. GraphPad Prism 8.0.1 was used to perform the above tests.

Identi cation of differential expression genes (DEGs)
A total of 359 DEGs were screened between paired lesional skin and non-lesional skin groups with |logFC|>1.5 and adjusted P<0.05, among which 284 genes were upregulated and 65 genes were downregulated. The volcano plots of DEGs has been shown in Fig. 2, and the heatmap of DEGs has been shown in Fig. 3.

GO and KEGG pathway analysis of DEGs
GO and KEGG pathway enrichment analysis were performed in Metascape. The top 20 GO enrichment items were classi ed into two functional groups: biological process group (15 items), and molecular function group (5 items) (Fig. 4A, B and Table 3). The DEGs were mainly enriched in six clustering groups in the GO enrichment analysis, including response to bacterium, defense response to other organism, antimicrobial humoral response, avonoid glucuronidation, skin development, monocarboxylic acid metabolic process (Fig. 4A). The top 20 KEGG pathways enrichment items are shown in Fig. 4C, D and Table 4. The DEGs were mainly enriched in three clustering groups in the KEGG pathways enrichment analysis, including Steroid hormone biosynthesis, NOD-like receptor signaling pathway, Cytokine-cytokine receptor interaction (Fig. 4C).
Subsequently, the correlation of these seven hub genes were analyzed in STRING database (Fig. 5C). In Fig. 5C, red line means known interaction that has been determined experimentally, and skyblue line represents known interaction from curated databases, these genes including ISG15, STAT1, OASL, IFIT1 and IFIT3. Moreover, yellow green line means text mining evidence and black line represents coexpression, these genes including CXCL10, CXCL1, ISG15, STAT1, OASL, IFIT1 and IFIT3. Light blue line stands for protein homology, these genes including IFIT1 and IFIT3.
Functional enrichment analysis of seven hub genes Seven hub genes were explored by Cytoscape (Table VI). GO and KEGG analyses of seven hub genes were performed by http: //www. bioinformatics. com.cn. The hub genes were enriched in the top ve of these terms, which based on their adjusted P-values were shown in chord plots (Fig 6, Table 7 and Table   8). For GO analysis the top ve terms were response to bacterium, defense response to other organism, antimicrobial humoral response, avonoid glucuronidation, and chemokine receptor binding (Fig. 6A, Table 7). For KEGG pathway analysis, CXCL10, CXCL1, ISG15, STAT1, and IFIT1 were mostly enriched in NOD-like receptor signaling pathway, Cytokine-cytokine receptor interaction, and Amoebiasis, except for OASL and IFIT3 (Fig. 6B). OASL and IFIT3 weren't enriched in the top ve KEGG pathways (Table 8).
Correlation between CXCL10 and the other hub genes Among the correlation of these seven hub genes, CXCL10 was in the central of the network (Fig. 5C), and had a higher rank (Table V), suggesting the correlation of CXCL10 with other hub genes. CXCL10 showed positive correlation with other six hub genes in GSE30999, GSE34248, GSE41662 and GSE50790 (Fig. 7).
Validation of the expression of hub genes in datasets The seven hub genes underwent expression validation in GSE30999, GSE34248, GSE41662 and GSE50790. Except for CXCL1, STAT1, and OASL in GSE30999, all other hub genes were signi cantly upregulated in psoriatic lesional skin samples from the four datasets (P < 0.05, 0.01, 0.001, or 0.0001; Fig.  8). This may be due to small samples in GSE50790(only 4 paired samples).

Validation of hub genes by qRT-PCR
The expression levels of seven predicted hub genes were validated by RT-qPCR(n=3). The results of qRT-PCT showed that the transcription levels of CXCL10, CXCL1, ISG15, STAT1, OASL, IFIT1 and ITIT3 were signi cantly up-regulated in 10 ng/ml M5 induced HaCaT cells (Fig. 9). The expression levels of these genes were consistent with the microarray results.

Discussion
An important nding of this study was the identi cation of hub genes and associated pathways in psoriasis through bioinformatics and RT-qPCR. The predicted seven hub genes contained CXCL10, CXCL1, ISG15, STAT1, OASL, IFIT1, IFIT3, and the mostly enriched pathways were NOD-like receptor signaling pathway, and Cytokine-cytokine receptor interaction. These hue genes and pathways may be vital for the pathogenesis of psoriasis. Some studies reported several genes as hub genes in psoriasis, partly the same with ours [20][21][22].
Among the seven hub genes, CXCL10 was a higher rank gene and had positive correlations with other six hub genes, suggesting it might be the most valuable gene. Furthermore, the RT-qPCR result of CXCL10 con rmed the prediction. CXCL10, also named IFN-γinduced protein 10 (IP-10), is a member of the CXC family of chemokines, which is considered to play an important role in in ammation via its T-cell chemotactic and adhesion properties [23]. Researches indicate CXCL10 is upregulated in the psoriatic skin lesions and serum [23,24]. It has been hypothesized that CXCL10 could be a good marker for psoriasis [25]. Although CXCL10 and CXCL1 belong to the chemokine "CXC" family [26], they play different roles in psoriasis. CXCL10 attracts T helper (Th) 1 cells, and CXCL1 attracts neutrophils [27]. It is reported that CXCL10 production by keratinocytes depends on STAT1 [28]. STAT1 is known as a member of the STAT family, involving in type I and type II interferon (IFN) signalling [29]. The expression of STAT1 is increased in lesional psoriatic skin, and it shows an important role in the pathogenesis of psoriasis [30].
The remaining four genes(ISG15, OASL, IFIT1, IFIT3) are antiviral genes, which may explain the relatively few skin viral infections found in the psoriasis patients. Ubiquitin-like protein ISG15 is an interferonstimulated protein that has a critical role for the control of microbial infections [31,32]. Some studies reported that ISG15 was elevated in lesional psoriatic skin compared to atopic dermatitis skin and healthy skin [33]. IFN-inducible oligoadenylate synthetases-like (OASL) protein belongs to atypical oligoadenylate synthetase (OAS) family, possesses antiviral activity, and boosts the innate immunity [34,35]. Gao LJ et al. also identi ed OASL as a hub gene [21], but it is rarely studied in psoriasis patients. Additionally, OASL wasn't enriched in the top 5 pathways in our study. Although the expression of OASL was up-regulated in M5 induced HaCat cells in the present study, further research is needed to con rm our results, including clinical specimens and animal experiment.
The interferon-induced proteins with tetratricopeptide repeats (IFITs) included IFIT1, IFIT2, IFIT3, and IFIT5 [36]. IFIT1 and IFIT3 are members of IFIT family, which regulate immune responses and function as essential anti-viral proteins [37]. The study indicated that IFIT3 binding to IFIT1 was vital for stabilizing IFIT1 expression, and was indispensable for inhibiting infection of viruses lacking 2'-O methylation [37]. In this study, IFIT1 was enriched in NOD-like receptor signaling pathway, but IFIT3 wasn't. Currently there are few researches on IFITs in psoriasis, while IFIT1 and IFIT3 are overexpressed in oral squamous cell carcinoma that promote tumor growth, regional and distant metastasis [36]. Therefore, this new nding warrants further study.
Nucleotide-binding and oligomerization domain (NOD)-like receptor is a kind of pattern-recognition receptor. It associates with various diseases related to infection and immunity [38]. Some studies show NOD-like receptor signaling pathway is highlighted in psoriatic epidermis [39]. Meanwhile, some researches have found cytokine-cytokine receptor interaction is related to the pathogenesis of psoriasis by combined transcriptomic analysis [40]. These results are consistent with ours.
However, there are some limitations in this study. For example, this study lacks clinical human specimens. Therefore, further experiments should be performed to con rm whether these hub genes can be used as biomarkers in psoriasis.
In conclusion, through an integrative bioinformatic analysis of the GEO data, this study identi ed important hub genes and signaling pathways. These hub genes might be used as potential biomarkers in psoriasis.

Conclusions
The pathogenesis of psoriasis may associate with the seven hub genes (CXCL1, ISG15, CXCL10, STAT1, OASL, IFIT1, IFIT3) and pathways such as NOD-like receptor signaling pathway, and Cytokine-cytokine receptor interaction. These hub genes may be used as potential biomarkers in psoriasis, especially CXCL10.

Declarations
Qichao Jian authored or reviewed drafts of the paper, and approved the nal draft.

Funding
Not applicable.

Availability of data and materials
The following information was supplied regarding data availability: The data is available at NCBI GEO: GSE30999, GSE34248, GSE41667 and GSE50790.
Ethics approval and consent to participate Not applicable.

Consent for publication
Not applicable.

Competing interests
The authors declare that they have no competing interests.

Figure 3
Heatmap of DEGs. The orange and blue colors represent upregulated and downregulated genes, respectively.    Expression of seven hub genes in four datasets. The green bar indicates non-lesional skin groups, and the red bar indicates lesional skin groups. Paired T-testing was performed to compare the means of two groups. *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001. NL, non-lesional skin; LS, lesional skin; ns, no signi cance.