Gene Set Enrichment Analysis and Ingenuity Pathway Analysis to verify the impact on Wnt Signaling Pathway in Taodan Granules Treated Psoriasis

Background: Taodan granules (TDGs), traditional Chinese herbals, are effective in treating psoriasis. However, mechanisms of TDGs remain indistinct. The current study aims to indicate the molecular mechanisms of TDGs in treating psoriasis. Methods: Primarily, transcriptional proling was applied to identify differentially expression genes (DEGs). The following was that we used Gene Set Enrichment Analysis (GSEA) and Ingenuity Pathway Analysis (IPA) analysis for functional enrichment analysis. Eventually, RT-PCR was performed to validation. Results: The results revealed that TDGs could regulated the Wnt signaling pathway to ameliorate skin lesions of imiquimod (IMQ)-induced psoriatic mouse models. IPA core network associated with Wnt signaling pathways in TDGs for psoriasis was established. Thereinto zeste homolog (EZH) 2, CTNNB1, TP63, and WD repeat domain (WDR) 5 could be considered as upstream genes in the Wnt signaling pathway. Conclusions: The Wnt signaling pathway with these above upstream genes might be potential therapeutic targets of TDGs for psoriasis. receptors; GSEA: Gene Set Enrichment Analysis; Hdac: histone deacetylase; IL: interleukin; IMQ: imiquimod; IPA: Ingenuity Pathway Analysis; KC: keratinocyte; KEGG: Kyoto Encyclopedia of Genes and Genomes; KLF: Kruppel-like factor; KRT: Keratin; LC-MS/MS: Liquid Chromatograph Mass Spectrometer/Mass Spectrometer; MAPK: mitogen activated protein kinase; NES: normalized enrichment score; PASI: psoriasis area and severity index; PDGF: platelet derived growth factor; PI3K: phosphoinositide 3-kinase; PPI: protein-protein interaction; SFMBT: Sex comb on midleg with four malignant brain tumor domains; SMAD: recombinant mothers against decapentaplegic homolog; SOX: SRY-box transcription factor; SOX: SRY-related HMG-box transcription factor; SPF: Specic pathogen-free; TCL: T-cell leukemia; TDG: Taodan granule; TGF-β: transforming growth factor-β; TP63: tumor protein p63; WDR: WD repeat domain; ZFP: zinc nger protein

thoroughly, and mark it as di200. All the diluted samples were centrifuged at 16000 rpm at 4℃ for 15 min. The supernatant was transferred to liquid vials and stored in a refrigerator at 4℃ for later use.

Animals
Speci c pathogen-free (SPF) -grade male BALB/c mice, weighted 25 ± 3 g, were offered by the Shanghai Medical Experimental Animal Center (SCXK Shanghai 2013-0016, Shanghai, China). The mice were housed in a germfree environment (temperature of 23 ± 2°C). Ethical approvals were obtained by the Ethics Committee of Yueyang Hospital a liated to Shanghai University of Traditional Chinese Medicine (No. YYLAC-2021-107).

Model establishment and interventions
In a nutshell, mice were divided into three groups: control group; psoriatic model group (IMQ group); and psoriatic model with TDGs treated group (IMQ + TDG group) (n = 4). After adaptive feeding, the hair on the back of mice was removed (2 × 2 cm 2 ). Establishment of psoriatic modeling was treated with 62.5 mg of 5% IMQ cream for 6 h on back skin, while mice in the control group were applied isodose petroleum jelly.
Mice in IMQ + TDG group were followed by intragastric administration of 1.8 g/kg TDGs, and mice in IMQ group were conducted with intragastric administration of 1.8 g/kg 0.9% NaCl solution.

Transcriptome sequencing
On the 12th day, mice were euthanized and the back skin was extracted for Transcriptome sequencing.
Illumina HiSeqTM 2500 was used for sequencing and the transcriptome analysis was performed by Shanghai OE Biotechnology Co., Ltd. Standardization disposal was performed using DESeq software (version 1.18.0) was used to standardize the gene count for each sample. The differentially expressed genes (DEGs) were screened according to the results of |log2FoldChange| > 1 and p-value < 0.05.

Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis
After obtaining DEGs, DAVID (https://david.ncifcrf.gov/) was used to analyse the signi cance of GO and KEGG analysis. A p-value of < 0.05 was considered as differentially statistical signi cance.

IPA
DEGs with |log2FoldChange| > 2 as well as p-value < 0.05, pathway analysis and the construction of protein-protein interaction (PPI) network of the expression data were executed with IPA (QIAGEN, Redwood City, CA, USA). Results with |z-score| ≥ 2 and/or overlap p-value < 0.05 were considered as differentially statistical signi cance. Upstream networks were described on the IPA database. 2.7. GSEA GSEA (https://www.gsea-msigdb.org/gsea/index.jsp) with probes ranked by signal-to-noise ratio and statistical signi cance determined by 1,000 gene set permutations was used to investigate differences among groups for exploring the potential molecular mechanisms and functions of DEGs [16,17]. A pvalue of < 0.01 and false discovery rate (FDR) with a q-value of < 0.05 were considered as differentially statistical signi cance.

RT-PCR
On day 12, the mice were euthanized with CO 2 inhalation. The back skin was extracted and preserved in TRIzol reagent kit. The Reverse Transcription System First Strand cDNA Synthesis Kit was used with 20.0 a reaction volume. Real-time uorescent PCR was performed for RT-PCR. The speci c experimental method was the same as the previous study [18]. The primer sequences were revealed in Additional le 2.
2.9. Statistical methods SPSS 24.0 (IBM Corp., Armonk International Business Machines, New York, USA) was used for analysing the data, described as mean ± standard deviation (SD). A t-test was used to compare the two groups, while a p-value of < 0.05 was considered as statistical signi cance. ug/ml Danshensu possessed the highest content in TDGs. A total ion ow chromatogram of TDGs was shown in Additional le 3. The structure of the major components was displayed in Additional le 4.
Afterwards, to further determine the effect of TDG treatment for psoriasis, the transcriptional pro les of skin lesions treated or nor treated with TDGs of psoriatic mouse models were detected and compared based on RNA sequencing [13]. The DEGs following TDG treatment were identi ed by RNA sequencing, and the data were stored in the SRA database (SRP292449). A total of 1233 DEGs were identi ed, with 539 upregulation as well as 694 downregulation (Additional le 5A). Biological characteristics of the potential targets for TDGs were estimated by KEGG and GO analysis. KEGG analysis indicated the most upregulated gene categories (Top 3) were for the neuroactive ligand-receptor interaction, estrogen signaling pathway, biosynthesis of unsaturated fatty acids, and glycine, serine, and threonine (Gly-Ser-Thr) metabolism; and the most downregulated gene categories (Top 3) were for cytokine-cytokine receptor interaction, Staphylococcus aureus infection, osteoclast differentiation, and chemokine signaling pathway (Additional le 5B). GO analysis manifested that the most upregulated gene categories (Top 3) were for the intermediate laments, keratin laments, cellular components, and structural molecule activity; and the most downregulated gene categories (Top 3) were in ammatory response, immune system process, extracellular space, and innate immune response (Additional le 5C). Therefore, the above results revealed the effect of TDG treatment for psoriasis was closely related to metabolism, in ammation as well as immune regulations.

Upstream analysis of DEGs following TDG treatment
To further reveal the therapeutic pathway of TDG in the treatment of psoriasis, we analysed the potential upstream regulatory functions of DEGs following TDG intervention. We selected DEGs with |log2FoldChange| > 2 and p-value < 0.05 based on the sequencing results to construct the core regulatory network in the IPA database (  (Fig. 3A). Given the above, it was con rmed that the therapeutic effect of TDGs on psoriasis might be closely related to the Wnt signaling pathway. Wnt signaling pathway has proved to be activated in the pathogenesis of in ammatory diseases including psoriasis vulgaris, atherosclerosis, rheumatoid arthritis, and sepsis [19]. The differentially expressions of the core enrichment mRNA of three relevant pathways of Wnt were displayed in Fig. 3B.
For further analyse the possible regulatory relationships in the three enriched Wnt related signaling pathways, we constructed relevant IPA core network. These genes associated with Wnt signaling pathway (GO: 0030178), consisted of complex, cytokines, enzymes, group/complex, growth factor, kinase, peptidase, transcription regulator, transmembrane receptor, and others, were used for constructing the core regulatory network on the IPA database. Platelet derived growth factor (PDGF) BB, Wnt, extracellular signal-regulated kinase (ERK)1/2, and phosphoinositide 3-kinase (PI3K) were shown to have the highest correlation in the regulations of other proteins in this pathway following TDG treatment (Fig. 3C). The proteins encoded by these genes associated with Wnt-activated receptor activity (GO: 0042813) including complex, G-protein coupled receptor, group, growth factor, kinase, transmembrane receptor, together with others were integrated into IPA core network. G protein-coupled receptors (Gpcr), epidermal growth factor (EGF), as well as protein kinase B (Akt) indicated higher levels of activities (Fig. 3D). Furthermore, Genes of complex, cytokine, enzyme, G-protein coupled receptor, group, growth factor, transcription regulator, transmembrane receptor and others were consolidated into IPA core network associated with Wnt signaling pathway (GO: 0090090). Axis inhibition protein (AXIN) 1, Wnt, along with histone deacetylase (Hdac) showed signi cant activities in regulating other proteins (Fig. 3E).
Similarly, we constructed IPA core network associated with Wnt signaling pathway (mmu04310), while genes consisted of complex, enzyme, group, kinase, ligand-dependent nuclear receptor, transcription regulator and others. In this pathway, CTNNB1, mitogen activated protein kinase (MAPK) 8, together with casein kinase (Ck) 2 played key regulatory roles in the Wnt signaling pathway following TDG treatment in psoriasis (Fig. 4C).
For verifying whether the six transcription regulators of upstream had signi cant corresponding functions in the TDG treatment for psoriasis, we conducted RT-PCR in the back skin lesions of IMQ-induced psoriatic mouse models. The results con rmed that compared with normal skin, the EZH2, CTNNB1 and WDR5 mRNA levels of psoriatic back lesions were speci cally elevated, while TP63 levels was dropped. Following TDG treatment, the decline in mRNA expression levels of EZH2, CTNNB1 together with WDR5, and the rise in mRNA expression levels of as TP63 was emerged (Fig. 5B).

Discussion
TDGs have favorable effects on relieving psoriasis. Our previous studies have also proved that TDGs are effective for psoriatic patients [8,9], while the in vivo results [13] con rmed that TDGs affected KC proliferation and in ammatory responses to alleviate IMQ-induced psoriatic symptoms in animal models. To integrally control the quality of TDGs, LC-MS/MS analysis was for application. The results indicated that TDGs contained Rutin, Caffeic acid, Cryptotanshinone, Tanshinone IIA, Formononetin, Liquiritin, Danshensu, and other active ingredients. Tanshinone IIA, Cryptotanshinone, and Danshensu are the three main bioactive components of Salvia miltiorrhiza Bunge, while all of them has proved to alleviate psoriasis [12,[20][21][22]. Formononetin proved to has strong anti-proliferation properties [23]. Previously, Rutin, Caffeic acid, along with Liquiritin were reported to decrease in ammation, while Liquiritin had a certain protective effect on skin injury, and inhibited angiogenesis [24][25][26][27]. However, the mechanisms of TDGs in the treatment of psoriasis remains vague, while the speci c action pathways of TDGs are worth further exploration. Therefore, the current study aims to use GSEA and IPA core network to investigate the transcriptional regulation mechanism, and systematically obtain molecular functions in TDGs.
Based on transcriptome sequencing, we constructed the core regulatory network in the IPA database for DEGs following TDG application. In order to further explore the possible relationships between regulatory factors and downstream molecules, upstream analysis in IPA predicted the active or inhibited upstream regulators following TDGs. Results revealed the most remarkable inhibited transcription regulators were KLF4, TCL1A and CEBPE. KLF4, a transcription factor, regulates a diverse array of cellular processes, including development, differentiation, proliferation, and apoptosis. Compared with non-psoriatic skin, KLF4 protein levels were signi cantly increased in psoriatic lesions in patients [28]. Excessive KLF4 can increase the level of histone H3 acetylation in Keratin (KRT) 17 promoter region by synergistic EP300, and mediate the over-expression of KRT17 in psoriatic lesions [29]. TCL1A, a stem cell marker, is abundantly expressed in embryonic stem cells, identi ed as an oncogene in various hematological malignancies, beyond that also speci cally expressed in proliferative solid tumors. TCL1A expression proved to be higher in colorectal cancer tissues than that in adjacent normal tissues, and signi cantly correlated with tumor differentiation, TNM stage and Ki67 positive rate [30]. CEBPE, one of the CCAAT/enhancer binding proteins, plays critical role in multiple physiological and pathological processes, and is highly expressed in tumor diseases [31]. The transcriptional activities of CEBPE are regulated by interactions with other transcription factors and/or post-translational modi cation (such as acetylation) [32]. For another, the most distinct activated transcription regulators following TDGs were PAX1, ZFP36 and SOX7. PAX1, a pivotal tumor suppressor gene, is responsible for regulating cell differentiation and maturation. PAX1 gene methylation detection can accurately monitor cervical cancer precancerous lesions, thus has high research value [33]. ZFP36 is a type of transcription factor with nger-like domains, associated with autoimmune diseases, arthritis and other syndromes in mice and humans, and has regulatory functions on gene expression, cell differentiation, and embryo development. Mice with conditional deletion of Tristetraprolin (TTP) (encoded by the ZFP36 gene) in KCs (Zfp36 / K14-Cre mice) pullulated exacerbated in ammation in the IMQ-induced psoriatic models. Furthermore, the Zfp36 / K14-Cre mice developed progressively spontaneous pathology, including systemic in ammation, psoriatic-like lesions, and dactylitis [34]. TTP was downregulated in broblasts deriving from psoriatic patients, when compared to those deriving from healthy individuals, while psoriatic broblasts exhibited abnormal in ammasome activities. The above phenomena proved to be related to ZFP36 promoter methylation [35]. SOX7 is a member of the SOX family of transcription factors. Studies have indicated that SOX7 is a tumor suppressor gene, down-regulated in most tumors. The regulatory mechanisms may be through regulating the transcription process mediated by Wnt-β-catenin signaling pathway, inhibiting tumor proliferation, migration, and invasion [36][37][38]. Given all this, these upstream transcription factors played central roles in the improvement of psoriasis following TDGs.
In order to further explore and analyse the changes of the gene sets following TDG intervention, we wielded GSEA to critical GO and KEGG. GSEA results con rmed that the treatment of TDGs was closely relative with the Wnt signaling pathway. Previous studies have con rmed that Wnt signaling pathway, as a crucial pathway of proliferative signal transduction, plays a negative regulatory role in psoriasis. βcatenin, an important transcription factor of Wnt signaling pathway, is translocated in the skin granular layer of psoriasis [39]. Inhibition of the β-catenin encoding gene CTNNB1 in HaCaT cells leads to the decrease of Cyclin D1 and Axin2 expressions, thereby inhibiting IL-22-induced cell proliferation [40]. After stimulating HaCaT cells, IWP-2 (Wnt inhibitor) can inhibit cell proliferation and secretion of proin ammatory factors, as well as promote cell differentiation [41]. Next, the core network analysis was performed using IPA, allowed us to analyse the coordinate expression changes at pathway levels [42]. The results indicated that PDGFBB, Wnt, ERK1/2, PI3K, Gpcr, EGF, Akt, AXIN1, Hdac, CTNNB1, MAPK8, as well as Ck2 in Wnt signaling pathway were considered as critical targets of TDGs.
For exploring the upstream regulation of Wnt signaling pathway following the intervention of TDGs, we applied the upstream analysis of IPA to forecast upstream transcription factors, and carried out experimental veri cation. Combined with results of veri cation via animal experiments, critical transcription factors with prominent enrichment included EZH2, CTNNB1 and WDR5, were remarkably down-regulated by TDGs. The methyltransferase EZH2 as a valid target for psoriasis therapy, consistent with our results, overexpresses in skin lesions of psoriatic mouse models and HaCaT cells. In vivo, GSK126, the inhibitor of EZH2, GSK126 can ameliorate the IMQ-induced psoriatic lesions. EZH2 contributes to the development of psoriasis by inhibiting the transforming growth factor-β (TGFβ)/recombinant mothers against decapentaplegic homolog (SMAD) pathway impairment of miR-125a-5p-mediated Sex comb on midleg with four malignant brain tumor domains (SFMBT) 1 inhibition [43][44][45]. CTNNB1 gene is located on chromosome 3p21, and the mutation of the gene exon 3 causes nuclear accumulation of β-catenin. Targeting CTNNB1 and subsequently affecting the downstream factors, CyclinD1 and Axin2, can inhibit IL-22-induced proliferation of HaCaT and HKC cells, provided diagnostic markers and novel targets for psoriatic treatment [40]. The WD40 protein family member WDR5 in the Wnt signaling pathway has several functions on tumorigenesis and development of multiple organ tumors. It has been demonstrated that overexpression of WDR5 is associated with poor prognosis in patients with esophageal squamous cell carcinoma (ESCC), while WDR5 may act as a potential novel prognostic biomarker for ESCC [46]. However, the impact on psoriasis of WDR5 has not been elucidated. At the same time, the expression of TP63 mRNA was signi cantly increased following TDGs by verifying. TP63 can act as transcription factors, activating or repressing expression from a variety of gene promoters, and is believed to be crucial for normal development of ectodermal derived structures such as skin and oral mucosa. It has been indicated the downregulation of TP63 mRNA in psoriatic lesions compared to both clinically normal skin from patients and matched healthy controls [47]. Gao et al. reported that the expression of TP63 in psoriatic lesions was increased after treatment, which was consistent with our veri cation results [48]. Notably, SOX11 and FOS of the signi cant upstream transcription factors enriched by IPA analysis revealed no obvious difference in mRNA levels in psoriatic lesions, normal skin, and lesions with TDGs treatment. Although differences were not statistically signi cant, but upstream subtle changes might lead to mutative expressions of downstream Wnt signaling pathways. In summary, we predicted and veri ed that TDGs alleviated psoriasis via regulating integrant upstream transcription factors in the Wnt signaling pathway.
Our data provide evidence that TDGs may improve psoriasis by regulating the Wnt signaling pathway. In the future, further experiments on this pathway can be conducted in vivo and in vitro to further explore the mechanism of TDGs.

Conclusion
The present study used modular pharmacology analysis to prove that TDGs regulated the Wnt signaling pathway to ameliorate skin lesions of IMQ-induced psoriatic mouse models, and identi ed the upstream regulators in the Wnt signaling pathway following TDG treatment.