The paeonol target gene autophagy-related 5 has a potential therapeutic value in psoriasis treatment

Data source and pre-processing

The Gene Expression Omnibus (GEO; http://www.ncbi.nlm.nih.gov/geo/) (Clough, 2016) repository was retrieved using the keywords of “psoriasis”, “lesion”, and “homo sapiens” to obtain microarray datasets from psoriasis patients. The eligible datasets (GSE30999, GSE13355, GSE14905, and GSE41662; Affymetrix platform; Table S1) were selected following screening criteria: (1) microarray expression profile from skin tissue of patients undergoing psoriasis; (2) ≥ 40 samples of lesion and non-lesion cases. The flow chart of data analysis in this study is shown in Fig. 1.

DEGs identification by meta-analysis

We employed the MetaQC package (https://cran.r-project.org/web/packages/MetaQC/index.html) in R 3.4.1 to implement quality control (QC). The principal component analysis and standardized mean rank score were also used to evaluate and screen data information. The DEGs were extracted using MetaDE.ES in MetaDE (Chang, Sibille & Tseng, 2013) package (https://cran.r-project.org/web/packages/MetaDE). The cutoffs of consistently DEGs identification were set as tau² = 0 and Q pval > 0.05 of heterogeneity test, false discovery rate (FDR) < 0.05 and —log₂fold change (FC)— > 0.263. Functional analyses of these DEGs, including the Gene Ontology (GO) biological process and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway, were performed by DAVID (Huangda, Sherman & Lempicki, 2009) (version 6.8; https://david.ncifcrf.gov/). P value < 0.05 was used as the threshold of significant enrichment.

Screening of module DEGs related to psoriasis

Weighted gene co-expression network analysis (WGCNA) (Chen et al., 2012) has widely applied to identify the gene module associated with diseases. GSE30999 was set as the training dataset while other thee datasets (GSE13355, GSE14905, and GSE41662) served as the validation datasets. The stable gene module linked with disease features was screened using WGCNA software (version 1.61; https://cran.r-project.org/web/packages/WGCNA/index.html) with the thresholds of gene number ≥ 25, cutHeight = 0.995, preservation z score >5, and correlation >0.6.

Protein-protein interaction (PPI) network construction

The STRING (Szklarczyk et al., 2017) database (version 10.5; https://string-db.org/) was used to examine the protein product interactions of genes in the module. The PPI network was constructed and visualized with Cytoscape (Shannon et al., 2003) (version 3.6.1; http://www.cytoscape.org/).

Selecting in database

The genes and pathways related to psoriasis and paeonol were selected from the Comparative Toxicogenomics Database (CTD, 2020 update; http://ctd.mdibl.org/). The immunity-related genes were selected from the AmiGO 2 (http://amigo.geneontology.org/amigo). The overlapping genes between the DEGs and the genes in the related pathways were selected and regarded as the target of further experiments. The MCODE plugin (http://apps.cytoscape.org/apps/mcode) was used to identify modules in the PPI network.

Cell line, culture conditions, treatment, and psoriatic model

Human immortalized keratinocyte HaCaT cell line (Central Laboratory of Peking University Third Hospital, Beijing, China) were maintained in complete Dulbecco’s modified Eagle’s medium (DMEM; Hyclone, Thermo Scientific, Epsom, UK). To mimic psoriatic conditions, HaCaT cells (1 ×10⁶cells/ml) were cultured in complete DMEM for 12 h in 6-well plates, followed with incubation in serum-free DMEM supplemented with 100 ng/ml IL-22 and 10 ng/ml TNF-α (PeproTech House, London, UK) for 24 h (Wu et al., 2020), with or without paeonol (National Institutes for Food and Drug Control, Beijing, China; 40 µg/ml) (Meng et al., 2017; Quan, 2019). All cells were incubated at 37 °C 5% CO₂.

Plasmid construction and cell transfection

The full-length coding sequence (CDS) of the human autophagy-related 5 (ATG5) gene was multiplied by PCR using the ATG5 specific primers with Bam HI/Xho I restriction enzyme sites (Table 1). The plasmid pc-ATG5 was constructed by cloning the PCR products into the pcDNA3.1 vectors (Genechem Co. Ltd, Shanghai, China). The short interfering RNAs (siRNA) targeting ATG5 (si-ATG5) and scramble sequences were purchased from the Genechem Co. Ltd. HaCaT cells (1 × 10⁵ cells/well) were transfected with si-ATG5, pc-ATG5, scramble sequences, and pcDNA3.1 vector for 24 h. Cell transfection was performed in triplicate using Lipofectamine 2000 regents (Invitrogen, Carlsbad, CA, USA). Psoriatic HaCaT cells (1 × 10⁵ cells/well) were stimulated with 100 ng/ml IL-22 and 10 ng/ml TNF-α for 24 h, with or without paeonol treatment.

Proliferation assay

Psoriatic HaCaT cells were treated with trypsin at 12 h and 24 h post-treatment. Cell Counting Kit-8 (CCK8; Dojindo, Japan) solution was added into cell culture and incubated for 2 h. Cell viability at 450 nm absorbance was analyzed using a microplate reader (Thermo Labsystems, Helsinki, Finland).

Cell apoptosis assay

Cell apoptosis was analyzed using flow cytometry and an annexin V-Cy5-labeled Apoptosis Detection Kit (Beyotime Institute of Biotechnology, Nanjing, China). Briefly, HaCaT cells were harvested and were then suspended in 5 µl of Annexin V-Cy5 and 5 µl of PI. A FACS Calibur flow cytometer (BD Biosciences, Franklin Lakes, NJ, USA) was used for apoptosis analysis.

Measurement of cytokines

The profiles of IL-1β and IL-6 in cell culture were detected using the enzyme-linked immunosorbent assay (ELISA) and commercial ELISA kits (Cusabio Biotech Corporation, USA). A microplate reader (Thermo Labsystems) was used for data analysis.

Protein extraction and Western blot analysis

After treatment for 24 h, HaCaT cells were lysed in triplicate using RIPA buffer (Beyotime Institute of Biotechnology, Beijing, China). Protein concentration was determined using a BCA protein assay kit (Thermo Fisher Scientific, Inc., Waltham, USA). Proteins (30 µg) were separated using 10% SDS-PAGE and were transferred onto the polyvinylidene fluoride (PVDF) membranes (Millipore, Billerica, MA, USA). Membranes were blocked using 5% skim milk (Beyotime) and incubated with the first primary antibodies including anti-ATG5 (1: 1000, Abcam), anti-Beclin 1 (1: 1500, Abcam), anti-p62 (1: 1000, Abcam), and anti-GAPDH (1: 10000, Abcam) at 4 °C overnight. The secondary incubation with HRP goat anti-rabbit IgG (1: 20000, Boster Biotechnology, Wuhan, China) was conducted at room temperature for 1 h. Protein bands were analyzed using the enhanced chemiluminescence (ECL) system.

RNA isolation and quantitative real-time PCR

Total RNA was extracted from psoriatic HaCaT cells using TRIzol Reagent (Invitrogen) at 24 h after stimulation or transfection. The RNA was reversely transcribed to cDNA using a PrimeScript RT-polymerase kit (TaKaRa, Dalian, China) to the manufacturer’s protocol.

The expression levels of genes were detected using an SYBR ExScript qRT-PCR Kit (TaKaRa). GAPDH was used as the internal control. Gene-specific PCR primer pairs were synthesized by Sangon (Shanghai, China; Table 1). PCR amplification was conducted according to the following reaction conditions: 95 °C for 5 min; 38 cycles of 95 °C for 30 s, 60 °C for 40 s, and 72 °C for 45 s. The relative expression levels of genes were analyzed using the 2-^△△Ct methods.

Statistical analysis

Each experiment was performed in triplicate in this study. Statistical analysis was performed using the GraphPad Prism 8.0 software (GraphPad Prism Inc., La Jolla, CA, USA). All data were expressed as mean ± standard deviation. Data differences across groups were analyzed using the ordinary one-way analysis of variance (ANOVA) test followed by the Tukey test as post hoc. P values <  0.05 were considered statistically significant.

Identification of consistently DEGs and functional analyses

Data analysis showed that four datasets exhibited even distribution and good quality (Fig. 2A). A total of 779 DEGs were identified (Table S1). The heatmap clustering of the DEGs in samples is shown in Fig. 2B. Functional enrichment analysis indicated that these genes were enriched with GO biological processes associated with regulation of phosphorus metabolic process (GO:0051174), negative regulation of cell proliferation (GO:0008285), and lipid biosynthetic process (GO:0008610; Fig. 2C), and KEGG pathways including PPAR signaling pathway (hsa03320), complement and coagulation cascades (hsa04610), and ECM-receptor interaction pathway (hsa04512; Fig. 2C).

WGCNA extraction of gene modules related to psoriasis

Before WGCNA, we found that the gene expression profiles in four datasets had high and positive correlations and connectivity (Figs. S1A and S1B). The soft-thresholding power = 10 when the square of the correlation coefficient = 0.9 and mean connectivity =1 (Figs. 3A–3B). Ten WGCNA modules were identified in the training dataset (Fig. 3C) and the validation datasets (GSE13355, GSE14905, and GSE41662; Figs. 3D–3F). Module preservation analysis showed that eight modules related to psoriasis traits. Four modules, including the black, white, red, and turquoise, had negative correlations with the psoriatic phenotype (preservation z score >5, correlation <−0.6, and p < 0.05, Fig. 4) and two modules, including blue (98 upregulated DEGs) and brown (79 upregulated DEGs) modules, had positive correlations with psoriasis (preservation z score >5 and correlation >0.6, Fig. 4).

PPI network construction

The PPI network was constructed using the 177 upregulated DEGs in the blue and brown modules. Accordingly, the PPI network consisted of 236 edges (interaction pairs) and 133 nodes (upregulated gene products; Fig. S2). These genes enriched in 31 biological processes, including cell cycle (GO:0007049) and apoptosis (GO:0006915; Table S2), and 12 KEGG pathways including Cell cycle (hsa04110), Natural killer cell mediated cytotoxicity (hsa04650), Cytokine-cytokine receptor interaction (hsa04060), and T cell receptor signaling pathway (hsa04660; Table 2). The 133 upregulated genes are listed in Table S3.

Selection of psoriasis-related and paeonol-targeted genes

To select the genes that have an important role in psoriasis, genes associated with immunity and psoriasis were identified from the 133 DEGs. A total of 119 DEGs (Table S3) were included in the lists of immunity-related genes (n = 3,324) and psoriasis-related genes (17,246). Also, nine common pathways related to psoriasis, paeonol, and DEGs were identified from the CTD, including the Jak-STAT signaling pathway (hsa04630), Cytokine-cytokine receptor interaction (hsa04060), and RIG-I-like receptor signaling pathway (hsa04622; Table 3). Only two genes (ATG5 and IL12B) were overlapped between the DEGs and genes in the CTD database.

Also, a total of 35 paeonol-targeted genes were identified from the CTD database and the network of paeonol-targeted genes is shown in Fig. 5A. ATG5 was the only overlapped gene between paeonol-targeted genes and the DEGs (Fig. 5A). ATG5 interplayed with IL-6, CASP3, TNF, Akt1, and mTOR (Fig. 5A). Besides, we identified an MCODE module (score=15.73) consisting of 19 genes including IL-6, CASP3, TNF, Akt1, and mTOR. This module was related to 15 pathways including Apoptosis (hsa04210), Toll-like receptor signaling pathway (hsa04620), MAPK signaling pathway (hsa04010), and RIG-I-like receptor signaling pathway (hsa04622; Fig. 5B). These results indicated that the ATG5 gene might have an important role in psoriasis and might be a target of paeonol.

Paeonol suppresses psoriasis

To investigate the influence of ATG5 and paeonol intervention on psoriasis, we established an in vitro psoriatic model using IL-22/TNF-α stimuli and treated it with paeonol. IL-22/TNF-α stimuli increased HaCaT cell viability (p = 0.0110; Fig. 6A), ATG5 expression (p < 0.0001, Fig. 6B), and the expression of IL-6 and IL-1β (p < 0.01; Figs. 6C–6F). However, paeonol significantly reversed IL-22/TNF-α-induced changes in HaCaT cells. Paeonol decreased ATG5 expression (p = 0.0026) and the expression of IL-6 and IL-1β in IL-22/TNF- α-treated HaCaT cells (p < 0.05; Figs. 6C–6F). We also found that IL-22/TNF-α stimuli induced less obvious reduction in the apoptotic percentage of HaCaT cells, from 4.58 ± 0.47% in control cells to 3.36 ± 0.07% in model cells (IL-22/TNF- α-treated HaCaT cells; p = 0.0550, Figs. 6G and 6H). The effect of paeonol on apoptosis was not evident (p = 0.3482; Fig. 6G and 6H). These results might indicate that paeonol has an inhibitory effect on psoriasis.

Si-ATG5 increases IL-6 and IL-1β and promotes apoptosis

We transfected si-ATG5 and pc-ATG5 into HaCaT cells and detected the expression of IL-6 and IL-1β to investigate whether the ATG5 gene could be used as a therapeutic target for the management of psoriasis. PCR analysis confirmed that si-ATG5 and pc-ATG5 transfection significantly decreased and increased the expression level of ATG5 in model cells, respectively (p < 0.0001; Fig. 7A). Also, si-ATG5 dramatically increased the contents of IL-6 and IL-1 β in model cells compared with negative controls (NC; p = 0.0011 for IL-6, and p = 0.0393 for IL-1β; Figs. 7B and 7C). However, we observed that pc-ATG5 transfection significantly decreased the contents of IL-6 and IL-1β in model cells compared with NC (Fig. 7B and 7C).

Besides, we found si-ATG5 significantly decreased HaCaT cell viability (p < 0.0001; Fig. 8A) and increased the apoptotic percentage of model cells (from 3.84 ± 0.15% to 14.45 ± 1.07%, p < 0.0001; Fig. 8B). However, pc-ATG5 transfection prevented the effects of IL-22/TNF-α stimuli on HaCaT cell viability and apoptosis (Figs. 8A and 8B). These results showed that ATG5 expression played crucial roles in controlling inflammation and cell viability in HaCaT cell apoptosis.

Expression of Beclin 1 and p62

To investigate the effect of ATG5 expression on autophagy in HaCaT cells, we detected the expression levels of two autophagy-related proteins, including Beclin 1 and p62. Western blot analysis showed that IL-22/TNF-α stimuli significantly increased the fold changes of ATG5, Beclin 1, and p62 proteins compared with control cells (p < 0.001; Figs. 9A–9D). Paeonol significantly decreased the expression of ATG5 and Beclin 1 proteins in HaCaT cells (p < 0.05; Figs. 9A–9C) but not p62 ( p = 0.3334; Fig. 9D). Besides, si-ATG5 transfection decreased the expression of Beclin 1, ATG5, and p62 proteins compared with NC (p < 0.0001; Figs. 9A–9D) and pc-ATG5 transfection increased Beclin 1, ATG5, and p62 proteins (p < 0.01; Figs. 9A–9D), respectively. However, the effects of paeonol on the expression of Beclin 1, ATG5, and p62 proteins in HaCaT cells transfected with si-ATG5 or pc-ATG5 were no evident.