Prediction of crucial epigenetically‑associated, differentially expressed genes by integrated bioinformatics analysis and the identification of S100A9 as a novel biomarker in psoriasis.

Psoriasis is one of the most common immune‑mediated inflammatory diseases of the skin. The identification of the pivotal molecular mechanisms responsible for the disease pathogenesis may lead to the development of novel therapeutic options. The present study aimed to identify pivotal differentially expressed genes (DEGs) and methylated DEGs in psoriasis. The raw data from gene microarrays were obtained from the Gene Expression Omnibus database. The data were processed using packages in Bioconductor. In total, 352 upregulated and 137 downregulated DEGs were identified. The upregulated DEGs were primarily enriched in the 'innate immune defense' response and the 'cell cycle'. The downregulated DEGs were primarily enriched in 'cell adhesion' and 'tight junction pathways'. A total of 95 methylated DEGs were identified, which were significantly enriched in the 'interleukin (IL)‑17 signaling pathway' and the 'response to interferon'. Based on a comprehensive evaluation of all algorithms in cytoHubba, the key epigenetic‑associated hub genes (S100A9, SELL, FCGR3B, MMP9, S100A7, IL7R, IRF7, CCR7, IFI44, CXCL1 and LCN2) were screened out. In order to further validate these genes, the present study constructed a model of imiquimod (IMQ)‑induced psoriasiform dermatitis using mice. The levels of these hub genes were increased in the IMQ group. The knockdown of methylation‑regulating enzyme ten‑eleven translocation (TET) 2 expression in mice attenuated the expression levels of S100A9, SELL, IL7R, MMP9, CXCL1 and LCN2. Furthermore, the hydroxymethylated level of S100A9 was highly expressed in the IMQ group and was significantly decreased by TET2 deficiency in mice. On the whole, using an integrative system bioinformatics approach, the present study identified a series of characteristic enrichment pathways and key genes that may serve as potential biomarkers in psoriasis.

Introduction

Psoriasis is one of the most chronic inflammatory hypertrophic skin diseases, affecting 2-5% of the population worldwide (1,2). Multiple factors may cause psoriasis, including genetics, epigenetic modifications, environmental triggers and immune dysfunction (3,4). Although there are extensive studies available regarding the etiologies and mechanisms underlying psoriasis, the precise molecular mechanisms remain elusive (5). Therefore, crucial molecular markers of psoriasis that are identifiable with more powerful technologies are urgently required.

The methylation and hydroxymethylation of DNA are the complementary processes of epigenetic regulation; these epigenetic alterations are responsible for the modulation of gene expression (6). Ten-eleven translocation (TET)2 is a key methylation-regulating enzyme that can alter the DNA methylation status by converting 5-methylcytosine (5-mC) to 5-hydroxymethylcytosine (5-hmC) and can regulate the expression of certain genes (7). It has previously been reported that TET2 plays an important role in biological processes, involving cell differentiation and proliferation (8). In a recent study, TET2 was also demonstrated to be involved in the progression of psoriatic dermatitis in mice by regulating hydroxymethylation (9). Although some methylation microarrays have been performed to identify differentially methylated genes (DMGs) in psoriasis (10), the comprehensive profile network analysis remains largely elusive.

Gene microarray is a high-throughput gene detection technique that has been widely used for the investigation of differentially expressed genes (DEGs), and is an important method for studying the pathogenesis of a variety of diseases (11). In recent years, a number of microarray gene profile studies on psoriasis have been performed. Using microarray data, the gene expression profiles in psoriasis have been disclosed and described (12). Certain studies have integrated analyses of gene expression microarray and methylation microarray by bioinformatics analysis in different types of cancer (13,14). However, to the best of our knowledge, there are currently no systematic studies available focusing on mining expression profiles and methylated DEGs in psoriasis. The integrated bioinformatics analysis will be more reliable and provide valuable new molecular targets to foster the development of accurate diagnosis and novel therapeutic strategies. The present study downloaded original data from the Gene Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/geo/). The integrated gene expression profile and DNA methylation profile microarrays of psoriasis were analyzed using a series of packages in R in order to identify the DEGs. Subsequently, The Database for Annotation, Visualization and Integrated Discovery (DAVID, https://david.ncifcrf.gov/), the Kyoto Encyclopedia of Genes and Genomes (KEGG, http://www.kegg.jp/), the Search Tool for the Retrieval of Interacting Genes (STRING https://string-db.org/) and Cytoscape were used to analyze the functional enrichment, protein-protein interaction (PPI) network, cluster module and hub genes. For further verification of the association of DMGs with epigenetic regulators in vivo, the present study created an animal model of imiquimod (IMQ)-induced psoriasiform dermatitis and knocked down TET2 expression. The results of the present study provide additional potential targets for the diagnosis and therapeutic intervention in patients with psoriasis.

Materials and methods

Data collection

Data on gene expression profiles (GSE13355, GSE14905 and GSE78097) and gene methylation profiles (GSE31835) were downloaded from the GEO database of The National Center for Biotechnology Information (NCBI). The 3 selected gene expression profile microarrays were based on the same platform (GPL570), which can decrease batch differences in the integration analysis. GSE31835 was based on the GPL8490 platform. GSE13355 contained the data from 58 skin lesions of patients with psoriasis and 64 normal healthy controls; GSE14905 contained 33 samples from patients with psoriasis and 21 normal control samples; GSE78097 contained 13 samples from patients with psoriasis and 6 normal control samples; 19 samples from lesions of patients with psoriasis and 10 normal skin samples were contained in GSE31835.

Data preprocessing and the identification of DEGs

The new versions of affymetrix annotation files were download from the official affymetrix website. The raw data of the gene expression profile were pre-processed by the Robust Multi-array Average (RMA) algorithm with the use of the affy package in Bioconductor (https://www.bioconductor.org/) for background correction and normalization (15). Linear models for the microarray data and corrections for the P-values were calculated using the limma package in Bioconductor (16). |log2FC|>1 and P<0.05 were set as the cut-off criteria to indicate statistical significance. The methylation profile microarray GSE31835 was analyzed using GEO2R software. The DMGs were identified with a P-value cut-off of <0.05 and |t| >2-fold. Subsequently, a principal component analysis (PCA) and hierarchical cluster analysis were performed with the scatterplot3d package and gplots package in R software.

Functional enrichment analysis

The DAVID provides gene analysis tools that link genomic information to biological annotations in order to help users identify the most valuable and significant annotations and links (17). DAVID was used for the analysis of DEGs of the KEGG pathway enrichment and Gene Ontology (GO) functional analysis, including cellular components, molecular function and biological processes (18). P<0.05 and counts ≥2 were used as the significance threshold.

PPI network construction and module analysis

STRING is an online tool for searching for known proteins and predicting the interaction between proteins (19). DEGs were mapped into STRING, and a PPI network was established by calculating a combined score of >0.4. The PPI network was constructed with Cytoscape (version 3.5.1) software. Significant modules were selected based on MCODE, with scores >4 and the number of nodes >5 in Cytoscape. The hub genes were identified according to the 12 topological analysis methods of cytoHubba (20).

Animal procedures

A total of 72 female BALB/c mice (weighing 18-20 g) were purchased from Beijing Vital River Laboratory Animal Technology Co., Ltd. The animals were housed in a controlled environment (23±2°C, 12-h light/dark cycle) with sufficient food and water at The Animal Centre of Tongji Medical College. Experimental mice were 8-10 weeks of age. All experimental protocols with animals were conducted following the ARRIVE guidelines. All mouse care and procedures were approved by the Ethics Committee for experimental animals, Tongji Medical College (no. SYXK HUBEI 2016-0057).

The BALB/c mice were randomly divided into 3 groups as follows and the mice in all groups were shaved on their backs (average area of 2×3 cm): i) The control group (n=8), mice were topically treated with vaseline and received an intradermal injection of lentivirus-delivered scrambled RNAs (shRNA-scr) on days 1 and 4; ii) the IMQ group (n=8), skin inflammation was induced by IMQ treatment, as previously reported (21,22). Briefly, mice in the IMQ group received a daily topical dose of 62.5 mg of 5% IMQ cream on their shaved backs for 7 consecutive days and received injections of the same volume of shRNA-scr; the IMQ + sh-TET2 group (n=8), mice received intradermal injections of lentivirus-delivered shRNA against TET2 (sh-TET2) followed by the application of IMQ cream for 7 days. The lentivirus-delivered shRNAs were designed and obtained from Shanghai Genechem Co., LTD. The sequences of the shRNAs were 5′-TTC TCC GAA CGT GTC ACG T-3′ (sh-sc), 5′-GAG CGT TCC TCA GTA TCA TTT-3′ (sh-TET2). On day 8, following the application of 20% urethane (1,500-2,000 mg/kg) (23), the dorsal skin lesions of the mice were taken for use in further experiments. Anesthetic mice were unresponsive to external stimuli, with the loss of righting and palpebral reflexes, but had a constant heart rate and respiratory rate (24). Subsequently, the mice were euthanized by rapid decapitation. The experiments were performed in triplicate.

Western blot analysis

The samples were lysed with cell lysis buffer (RIPA Lysis Buffer; cat. no. P0013B; Beyotime), and the supernatant was collected following centrifugation at 13,000 × g for 5 min at 4°C. The concentrations of proteins were measured using a BCA protein assay kit (cat. no. G2026; Servicebio). Equal quantities of protein were separated via SDS-PAGE (stacking gel, 5%; separating gel, 8%-TET2, 12%-S100A9) and then transferred onto polyvinylidene difluoride membranes. After blocking with 5% skim milk at room temperature for 1 h, the membranes were incubated with primary antibodies against GAPDH (1:1,000; cat. no. AF7021; Affinity Biosciences), TET2 (1:1,000; cat. no. ab94580; Abcam) or S100A9 (1:1,000; cat. no. NB110-89726SS; Novus Biologicals) at 4°C overnight, followed by the appropriate horseradish peroxidase-conjugated secondary antibodies (anti-rabbit IgG, HRP-linked antibody/anti-mouse IgG, HRP-linked antibody 1:2,000; cat. no. 7074/7076; Cell Signaling Technology) for 1 h at room temperature. Signals were detected using DNR bio-imaging systems according to the manufacturer's protocol. The densitometric quantification analysis of proteins was performed using ImageJ software (Ver. 1.45s; National Institutes of Health).

Reverse transcription-quantitative PCR (RT-qPCR)

Total RNA was extracted from the whole skin samples using TRIzol® reagent (Invitrogen; Thermo Fisher Scientific, Inc.). cDNA was synthesized by a reverse-transcriptase reaction with a PrimeScript RT Reagent kit (Takara Biotechnology). The specific cDNA fragments were amplified using a Real-time PCR Detection System and Power SYBR-Green (Takara Biotechnology) with the relevant primers. The data are expressed as the relative expression to the GAPDH gene. The relative expression level of each target gene was analyzed using the 2−ΔΔCq method and normalized to GAPDH (25). The PCR primer sequences are presented in Table SI.

Hydroxymethylated DNA immunoprecipitation (hmeDIP)

The extraction of genomic DNA was performed using the SimpleDIP™ Hydroxymethylated DNA IP kit (Cell Signaling Technology) according to the manufacturer's protocol. Genomic DNA was fragmented to sub 100-500 bp using an Ultrasonic Homogenizer. The 5-hmC mouse monoclonal antibody (1:1; cat. no. 51660, Cell Signaling Technology) was added to the single stranded DNA samples. The samples were incubated overnight at 4°C with rotation. After washing and eluting, the immunoprecipitated DNA was purified using spin columns. RT-qPCR was then utilized for the detection of specific hydroxymethylated DNA levels. Quantitative analysis was performed using the following formula: 10% × 2(C[T] 10% Input Sample-C[T] IP Sample).

Statistical analysis

Comparisons between groups were performed using an unpaired two-tailed Student's t-test. In some analyses, one-way analysis of variance (ANOVA) was used, as indicated. Tukey's post hoc test was used for multiple comparisons among diverse groups. Data are expressed as the means ± standard deviation, and P<0.05 was considered to indicate a statistically significant difference. GraphPad Prism software (version 6.0; GraphPad Software) was used for statistical analysis.

Results

Identification of DEGs in psoriasis

Following quantification and normalization, DEGs were identified based on a threshold of P<0.05 and fold change >2. The present study screened 951, 1,530 and 2,836 DEGs in the psoriasis samples compared with the normal samples in the GSE13355, GSE14905 and GSE78097 datasets, respectively. Using GSE13355 as a representative, the results of the PCA and hierarchical cluster analysis of the DEGs are presented in Fig. 1A and B. The other two PCAs and heatmaps are presented in Fig. S1. In total, 498 common DEGs were screened out, including 352 upregulated, 137 downregulated and 9 non-specific genes (Fig. 1C; Table SII).

Figure 1

Identification of DEGs in psoriasis. (A) Principal component analysis of the samples in the GSE13355 dataset. The results revealed that psoriasis samples and normal samples could be well distinguished according to these DEGs. Orange dots represent normal samples and gray dots represent psoriasis samples. (B) Hierarchical cluster analysis of DEGs in GSE13355. Gene expression levels were indicated by colors: Red represents high expression level, green represents low expression level. The orange area represents normal samples and the gray area represents psoriasis samples. (C) Venn diagram of overlapping DEGs from 3 microarray datasets (GSE13355, GSE14095 and GSE78097). A |log2FC|>1 and P-value <0.05 were set as the cut-off criteria. DEGs, differentially expressed genes; FC, fold change.

DEGs GO analysis and pathway enrichment analysis

In order to investigate the biological functions of the DEGs, the present study performed a functional pathway enrichment analysis of the DEGs through GO analysis and KEGG pathway enrichment in DAVID. A P<0.05 was considered to indicate a statistically significant difference. As presented in Fig. 2A, in the biological processes group, the upregulated DEGs were primarily involved in 'the type I interferon signaling pathway', 'the defense response to virus', 'mitotic nuclear division', 'the immune response' and 'the inflammatory response'; in the molecular function group, the upregulated DEGs were significantly involved in 'chemokine activity', 'serine-type endopeptidase activity', 'CXC motif chemokine receptor (CXCR) 3 chemokine receptor binding', 'advanced glycosylation end-product specific receptor (RAGE) receptor binding' and 'regulatory region DNA binding'; in the cellular component group, the upregulated DEGs were primarily involved in 'cytosol', 'condensed chromosome kinetochore', 'spindle microtubule' and 'cornified envelope'. The downregulated DEGs were primarily involved in 'cell adhesion', 'positive regulation of protein kinase B signaling', 'RNA polymerase II transcription factor binding', 'transcription factor binding' and 'extracellular matrix' (Fig. 2B).

Figure 2

Enriched GO terms and KEGG pathways of upregulated and downregulated DEGs identified in psoriasis. (A) Significant enriched GO terms of upregulated genes. Blue, biological processes; green, molecular function; yellow, cellular components. The top 10 significant terms of each category in upregulated DEGs are displayed. (B) Significant enriched GO terms of downregulated genes. Blue, biological processes; green, molecular function; yellow, cellular component. The top 5 significant terms of each category in downregulated DEGs are displayed. (C) Significant enriched KEGG pathways of upregulated and downregulated DEGs associated with psoriasis. P<0.05 was considered to indicate a statistically significant difference. GO, gene ontology; DEGs, differentially expressed genes; KEGG, the Kyoto Encyclopedia of Genes and Genomes.

The significantly enriched KEGG pathways of upregulated DEGs associated with psoriasis were the 'influenza A pathway', 'chemokine signaling pathway', 'cell cycle pathway', 'RIG-I-like receptor signaling pathway' and 'p53 signaling pathway'. The downregulated genes were primarily enriched in 'tight junction pathways' (Fig. 2C).

Module analysis and hub gene selection from the PPI network of DEGs

In order to screen out key modules and hub genes, a PPI network of 475 DEGs was established with a combined score of >0.4, including 475 nodes and 2,556 edges (Fig. S2), which was based on the STRING database. A module analysis was performed using MCODE in Cytoscape. MCODE scores >4 and nodes >5 were set as the cut-off criteria. The 4 significant modules are presented in Fig. 3. Subsequently, the present study analyzed the pathways involved in modules 1 and 2 using ClueGO, the plugin of Cytoscape. The 'Cell cycle', 'NOD-like receptor signaling pathway' and 'RIG-I-like receptor signaling pathway' were the significantly enriched pathways (Fig. S3). The hub genes were screened by cytoHubba in Cytoscape (data not shown) and it was revealed that TOP2A, CDK1, IL8, KIAA0101, CCNB1, BIRC5, STAT3, TYMS, UBE2C, ABCA12, TMPRSS11D and S100A9 may all play a pivotal role in the molecular mechanisms underlying psoriasis.

Figure 3

The top 4 modules from the PPI network. (A) Module 1 contained 48 nodes and 1,106 edges. (B) Module 2 contained 29 nodes and 383 edges. (C) Module 3 contained 13 nodes and 61 edges. (D) Module 4 contained 6 nodes and 12 edges. MCODE scores >4 and nodes >5 were set as the cut-off criteria. PPI, protein-protein interaction; DEGs, differentially expressed genes.

Identification of abnormal methylated-DEGs in psoriasis

The gene methylation microarray GSE31835 was obtained from the public GEO database. Based on the cut-off standards, a total of 2,730 DMGs was primitively identified between the psoriasis samples and healthy controls. By the construction of a Venn diagram (data not shown), the consistent differential genes between DEGs and DMGs were obtained (Table I).

Table I

Abnormal methylated-differentially expressed genes in psoriasis.

Term

Gene name

Methylated-differentially expressed genes

TNIP3, XDH, HSD17B2, HMGCS2, SELL, FOXM1, C10orf99, IFI27, PPP1R1B, ADRB1, PI15, IRF7, CCND1, ZNF677, PI3, PDK4, CCR7, PDZK1IP1, IL19, ATP11B, CRABP2, GALNT6, CHI3L2, FUT3, CLEC7A, GBP1, NFE2L3, OAS2, KRT16, GZMB, FCGR3B, SMPD3, SLAMF7, LRRC8B, ZDHHC21, STS, ALDH1A3, LRG1, CD47, RORC, S100A7, S100A12, PRELP, SERPINB3, IVL, APOC1, CYP4B1, ARSF, TRIM22, SLC23A2, TIMP3, GJB6, APOL6, SLC2A13, FOXC1, DNASE1L3, KLHDC7B, FADS2, ARNTL2, CCL27, ZBTB16, KYNU, CCL20, CRISPLD1, TMEM40, CHAC1, LYN, BTC, TCN1, ID4, SLC7A5, GDPD3, GJB2, GBP6, WNT2B, ISG20, CXCL1, CDK5R1, TMEM116, LYPD5, IL7R, PTAFR, KLK10, ATP12A, SPRR3, CFB, AIM2, LCN2, KRT6B, UHRF1, MAMDC2, KCNJ15, IFI44, S100A9, MMP9

Identification of key candidate genes and pathways with the methylated-DEGs PPI network

A PPI network of 95 methylated-DEGs was constructed using Cytoscape software (Fig. 4A). According to the ClueGO analysis, these genes were primarily involved in the 'interleukin (IL)-17 signaling pathway', 'Staphylococcus aureus infection', 'response to interferon' and 'granulocyte chemotaxis' (Figs. 4B and S4). GBP1, IRF7, OAS2, IFI44, GBP6, IFI27 and ISG20 were selected as significant module genes using MCODE (Fig. 4C). The top 10 genes, calculated by Betweenness algorithms, are presented in Fig. 4D. S100A9, SELL, FCGR3B, MMP9, S100A7, IL7R, IRF7, CCR7, IFI44, CXCL1 and LCN2 were the top outstanding key genes based on a comprehensive evaluation of 12 algorithms in cytoHubba, which were selected when the genes contained in five or more algorithms.

Figure 4

The PPI network for methylated-differentially expressed genes. (A) PPI network of 95 methylated-differentially expressed genes was constructed by STRING and reconstructed by Cytoscape. (B) KEGG pathway enrichment analysis of methylated-differentially expressed genes by ClueGO software. P<0.05 was considered to indicate a statistically significant difference. (C) Module analysis of methylated-differentially expressed genes by the MCODE in Cytoscape, including 7 nodes and 17 edges. Scores >4 and nodes >5 were set as the cut off criteria. (D) The top 10 hub genes were calculated by Betweenness algorithms in Cytoscape. The deeper the color, the more important the gene. PPI, protein-protein interaction; KEGG, the Kyoto Encyclopedia of Genes and Genomes; STRING, Search Tool for the Retrieval of Interacting Genes/Proteins.

Evaluation of the hub DEGs in IMQ-induced psoriasiform dermatitis

TET2 is a key DNA methylated regulatory enzyme that converts methylation to hydroxymethylation, thereby regulating gene expression. In order to verify the hub DMGs screened in the present study, a mouse model of IMQ-induced psoriasis with or without the knockdown of TET2 expression was established. Compared with the control group, the IMQ group exhibited significant psoriasiform dermatitis of the erythema, thickening and scaling, whereas the skin lesions of lentivirus-delivered shTET2-injected mice were significantly decreased (Fig. 5A). The mRNA levels and protein concentration of TET2 were decreased in the sh-TET2-injected mice compared with the IMQ group (Fig. 5B, D and E). Due to the significantly high expression of S100A9 in the IMQ group, the western blot bands were slightly conjoined, which may affect the accurate quantification of the protein. The expression levels of the hub genes were detected by RT-qPCR. The mRNA expression levels of SELL, FCGR3B, MMP9, S100A7, IL7R, IRF7, CCR7, IFI44, CXCL1 and LCN2 were increased in the IMQ group compared with the control group. The expression levels of SELL, IL7R, MMP9, CXCL1 and LCN2 were suppressed by sh-TET2 treatment (Fig. S5). Furthermore, compared with the control group, the IMQ group exhibited a significantly increased expression of S100A9, which was evidently suppressed by sh-TET2 treatment at both the mRNA and protein levels (Fig. 5C, F and G). When investigating the molecular mechanisms underlying the effects of TET2 on S100A9 expression, it was revealed that the hydroxymethylation level of S100A9 was significantly decreased in the IMQ + sh-TET2 group (Fig. 5H). These data thus indicate that S100A9 may be regulated by epigenetic approaches and may play an important role in the pathogenesis of psoriasis.

Figure 5

TET2 regulates the expression of S100A9 in a model of psoriasiform dermatitis. (A) Representative lesions of BALB/c mice following continuous treatment for 8 days. (B and C) TET2 and S100A9 mRNA levels in the skin lesions were evaluated by RT-qPCR. The protein expression levels of (D and E) TET2 and (F and G) S100A9 in the skin lesions were measured by western blot analysis. (H) The DNA hydroxymethylation level of S100A9 was analyzed by hMeDIP-qPCR. The results are presented as the means ± standard deviation. All experiments were performed 3 times, and the representative results are presented. *P<0.05, **P<0.01 and ***P<0.001. S100A9, S100 calcium binding protein A9; TET2, ten-eleven translocation 2; IMQ, imiquimod; hMeDIP-qPCR, Hydroxymethylated DNA Immunoprecipitation-quantitative PCR; shRNA-scr, lentivirus-delivered scrambled RNAs; shRNA-scr, control group; IMQ+shRNA-scr, IMQ group; IMQ+sh-TET2, IMQ+sh-TET2 group.

Discussion

The present study analyzed three gene expression profiles and one gene methylation profile. A total of 352 upregulated DEGs and 137 downregulated DEGs were identified. In a comparison between the DEGs and DMGs, 95 overlapping genes were selected. The IL-17 signaling pathway and Staphylococcus aureus infection pathway were the main enriched pathways. The biological processes involved were the response to interferon and granulocyte chemotaxis. S100A9, SELL, FCGR3B, MMP9, S100A7, IL7R, IRF7, CCR7, IFI44, CXCL1 and LCN2 were the top outstanding genes. Furthermore, S100A9 was regulated by epigenetic regulator TET2 in a mouse model of psoriasis.

As the three raw gene expression microarrays had the same normalization and quality control, the present study obtained highly standardized and comparable data. Gene expression datasets combined with comprehensive bioinformatics analysis can provide reliable and new insights into potential crucial diagnostic and therapeutic biomarkers. Researchers have previously focused on key genes and pathways in different diseases (26). According to the analysis of pathways and GO terms in the present study, DEGs were enriched in the 'immune response', 'inflammatory response', 'cell adhesion', 'chemokine signaling pathway' and 'cell cycle', which is consistent with extensive previous studies on psoriasis (5,27,28). Furthermore, all the top 3 significant terms of the biological process analysis and one of the main pathways of the upregulated DEGs were associated with responses to virus, suggesting that the pathogenesis of psoriasis may be closely associated with viral infection. A recent study revealed that activated RIG-I induced the production of IL-23 and triggers psoriasiform dermatitis in mice, highlighting that the dysregulation of antiviral immune responses may initiate the inflammation of psoriasis (29). Extracellular vesicles transport cytosolic proteins, lipids and RNA, mediating intercellular communication between different cells (30). Circulating endothelial- and platelet-derived extracellular vesicles are increased in psoriasis (31) and are decreased following effective anti-TNF-α treatment (32). The present study demonstrated that downregulated DEGs were evidently involved in extracellular vesicles (OGN, APOE, EDIL3 and PRELP) in the GO cell component analysis. Thus, these molecules may be potential biomarkers of psoriasis.

Aberrant DNA methylation can influence the expression of hub genes, and can induce a variety of pathophysiological alterations (33). Therefore, the present study identified 95 abnormally methylated DEGs by overlapping the DEGs and DMGs. These genes were primarily involved in the 'IL-17 signaling pathway', 'response to interferon' and 'granulocyte chemotaxis', which is reasonable, as these biological processes have been demonstrated to increase the risk of developing psoriasis (34). The top 11 hub genes were selected by cyto-Hubba analysis. S100A7 is a calcium-binding protein that participates in a number of biological functions, such as proliferation and inflammation (35). S100A7 has been demonstrated to be upregulated in psoriasis and is also named psoriasin, as it was originally discovered in psoriasis (36). The knockdown of the histone demethylase JMJD3 in keratinocytes has been shown to result in the significant downregulation of the expression of S100A7 (37). MMP9 is an enzyme from the gelatinase family that regulates cell matrix composition and proliferation (38). Elevated concentrations of MMP9 have been described in patients with psoriasis (39). The aberrant methylation of the promoter region of MMP9 has been shown to decrease the levels of RNA and protein secretion in lymphoma cells (40). FCGR3B, IL7R, IRF7, CCR7 and LCN2 have also been shown to be highly expressed in patients with psoriasis (41-45). IRF7 is a key transcriptional regulator of type I interferon (IFN)-dependent immune responses and plays a vital role in the response against viruses (46). TET2 can maintain the hypomethylation of IRF7 to regulate the TLR7/9-elicited IFN response in plasmacytoid dendritic cells (pDCs) (47). Hence, these genes may be potential candidate genes regulated by epigenetics that are involved in the pathogenesis of psoriasis.

S100A9 is a calcium-binding protein of the S100 family, which plays a prominent role during the amplification of the inflammatory response in autoimmunity, including chemotaxis, adhesion and the degranulation of neutrophils (48). Through a comprehensive bioinformatics analysis, S100A9 was one of the identified hub genes of methylated DEGs. Furthermore, the research in the present study demonstrated that the knockdown of epigenetic regulator, TET2, in mice significantly decreased the expression of S100A9, both at the mRNA and protein level. Compared with the IMQ group, the hydroxymethylation level of S100A9 was also significantly decreased in the IMQ + sh-TET2 group. These findings highlight the importance of the regulation of S100A9 by an abnormal methylation status.

To the best of our knowledge, there are limited studies that have used integrated bioinformatics analyses to gather valuable information and provide important clues for future research. Using computational systems biology, Sevimoglu and Arga reported that the JAK/STAT signaling pathway was significant for psoriasis. They screened out several biomarkers and therapeutic candidates by analyzing the raw data of gene expression datasets (49). Certain feature algorithms filtered gene expression profiles in order to acquire classification models of psoriasis (50). Unlike previous studies, the present study comprehensively analyzed DEGs and DMGs to identify active modules and key epigenetically-associated DEGs.

The integrated bioinformatics analysis of both the gene expression profile and the methylation profile is more reliable and provide valuable, novel molecular targets to foster the development of accurate diagnosis and novel therapeutic strategies. However, one limitation of this study was the lack of clinical human specimens. Additional studies with clinical materials and experiment validation are required in order to verify the key candidate genes. Furthermore, the present study analyzed only a single methylation profile microarray; thus, more comparable microarray data are required so that more valuable data can be obtained. The present study selected S100A9 as a representative gene, and preliminarily evaluated the expression and epigenetic regulatory mechanism of S100A9 in psoriasis. Further experiments are required to clarify the exact association between TET2 and S100A9, and the molecular mechanisms of action.

In conclusion, the results from the present study not only identify a series of DEGs and DEMGs, but also analyze the significant modules and key candidate genes. In addition, in order to further verify the bioinformatics analysis data, the present study detected the expression levels of hub genes (S100A9, SELL, FCGR3B, MMP9, S100A7, IL7R, IRF7, CCR7, IFI44, CXCL1 and LCN2) in a mouse model of psoriasiform dermatitis by interfering with the expression of methylation-regulating enzyme TET2. These key genes may serve as potential biomarkers in psoriasis.