TGM1 as a novel signature gene in psoriasis identified by integrative bioinformatics and experimental validation

Introduction

Psoriasis is an immune-mediated dermatological condition characterized by erythematous and scaly plaques affecting the entire integumentary system (1). Empirical observations indicate that it severely impacts 2–4% of the population, leading to serious psychological anguish and diminished quality of life (2,3). Numerous contributing factors have been identified in this context (4). Genetic underpinnings and external influences, including smoking cessation, anxiety and alcohol consumption, have been shown to contribute to the development of psoriasis (4).

Systematic exploration has revealed that psoriasis is predominantly regulated by interactions among numerous cytokines and multifaceted signaling pathways (5). Notably, studies have demonstrated that several cytokines, including IL-2, IL-6, IL-17, IL-23, TNF-α and IFN-γ, are important for the advancement of psoriasis (6,7). Consequently, biological therapies targeting these cytokines, particularly the TNF-α, IL-12/IL-23, IL-17 and IL-23/IL-39 pathways, have been recognized as promising treatments for psoriasis (8,9). The Janus kinase (JAK)/STAT pathway is a primary inflammatory mechanism involved in the development of psoriasis. Specific cytokines, particularly IL-17 and IL-23, are important in this process by relaying signals and regulating the transcriptional expression of targets within the JAK/STAT pathway. This process is important for the advancement of psoriatic disease (10). Signaling pathways, including MAPK and PI3K/AKT signaling, also contribute to the development of psoriasis by influencing the disease state (11,12). Notably, while existing research has yielded some insights into psoriasis mechanisms, current findings remain insufficient. Therefore, it is key to identify novel diagnostic markers to clarify the pathogenesis of psoriasis.

The present study utilized the GSE30999, GSE53552 and GSE13355 datasets from the Gene Expression Omnibus (GEO) database containing lesional skin (LS) and non-lesional skin (NLS) tissues from individuals with psoriasis. The datasets were divided into a training set, GSE30999, and two validation sets, GSE53552 and GSE13355. The analysis of these datasets employed the following methods: Variable selection, model training, protein-protein interaction (PPI) analysis, Gene Set Enrichment Analysis (GSEA) and single-gene immune infiltration analysis. To specifically address the course of psoriatic disease, 163 distinct differentially expressed genes (DEGs) were identified. The DEGs were examined utilizing the following machine learning algorithms: Least absolute shrinkage and selection operator (LASSO) regression, support vector machine-recursive feature elimination (SVM-RFE) and random forest (RF). The findings from these investigations yielded notable insights into the mechanism of transglutaminase 1 (TGM1). To clarify the role of TGM1 during psoriasis, an investigation of PPIs was conducted. Furthermore, single-gene GSEA and single-gene immune infiltration analysis were implemented to elucidate the underlying activities and biological pathways of TGM1.

Materials and methods

Data collection

Gene expression profiling data of psoriasis were retrieved from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). The GSE30999, GSE13355 and GSE53552 datasets were used in the present study. GSE30999 (13) includes 85 patients with psoriasis including 85 biopsy samples of LS and 85 matched biopsy samples of NLS. GSE13355 (14) includes 58 patients with psoriasis including 58 biopsy samples of LS and 58 biopsy samples of NLS. GSE53552 (15) includes 25 patients with psoriasis including 25 biopsy samples of LS and 24 biopsy samples of NLS. GSE30999 was selected as training dataset and both the GSE13355 and GSE53552 datasets were selected as validation datasets.

Identification of DEGs

‘Limma (version 3.58.1)’, a package of R software (version 4.3.0, http://www.r-project.org/), was used to identify DEGs in LS and NLS of patients with psoriasis (16). The present study used |log2[fold change (FC)]|>1 and P<0.05 as cutoff thresholds to extract statistically significant DEGs. However, when looking for signature genes, not all DEGs were included in later calculations. For example, Guan et al (17) created three algorithms to choose potential genes for lung cancer from a total of 51 PPI-related DEGs. Wei et al (18) adopted three methods to analyze the DEGs of the most important modules in the results of a weighted gene co-expression network analysis to screen significant variables. Therefore, the present study only used the significant DEGs for subsequent analysis. DEGs with |FC|>2 were initially isolated, but the top DEGs with |log2FC|>3 were then screened for the following analysis in order to prevent interference of genes with little variation and to reduce the workload of the calculations. Eventually, the R packages ‘pheatmap’ and ‘ggplot2’ were used to depict the DEGs in a heatmap and a volcano plot.

GSEA

GSEA is a computer tool that assesses the type of gene expression in a certain functional gene set and reveals the underlying biochemical pathways driving complex disease (19). In the present study, the ‘clusterProfiler’ package (4.10.0) was employed to perform GSEA (20). The reference gene set was ‘c2.cp.kegg_medicus.v2024.1.Hs.symbols.gmt’ (gsea-msigdb.org/gsea/msigdb/human/collections.jsp#C2), downloaded from the MsigDB database (https://www.gsea-msigdb.org/gsea/msigdb). To avoid interference from unchanged genes, the top 10,000 genes from all probes in the GSE30999 dataset were selected for GSEA based on |logFC| descending order. Kyoto Encyclopedia of Genes and Genomes (KEGG, genome.jp/kegg/) pathway terms with a P-value <0.05 were selected as significant.

In addition, the present study investigated the potential function of signature genes in psoriasis via single-gene GSEA. The ‘clusterProfiler’ R package was utilized to calculate the correlation between signature and other genes and the genes were then ordered according to their correlation scores. The ranking genes were chosen as the test gene set for KEGG enrichment analysis. P<0.05 was used as a criterion for significant items.

Functional and pathway enrichment analysis

The present study used Gene Ontology (GO) and KEGG analyses to unravel the functional implications and potential biological processes (BPs) of the significant DEGs with |log2FC|>3. GO analysis, a widely recognized bioinformatics tool, makes it possible to systematically classify and annotate genes according to BP, cellular component (CC) and molecular function (MF) (21). By subjecting genes to GO and KEGG pathway analysis, the present study identified the notable BPs and pathways that were substantially enriched in these genes. The present study used the bioinformatics website (bioinformatics.com.cn/) to perform GO and KEGG pathway analysis. Only GO and KEGG terms with P<0.05 were kept as significantly enriched. The top 10 results in GO and KEGG pathway analyses were selected and plotted as bubble plots and gene-pathway association network diagrams.

PPI network analysis

The identified genes were imported into the STRING database (version 12.0) (https://cn.string-db.org/) (22) to obtain the PPI network, with a minimum required interaction score of ≥0.4. The PPI network was visualized using Cytoscape v3.10.1 (23).

Screening of signature genes using machine learning methods

To further identify the signature genes for psoriasis, the present study adopted LASSO, SVM-RFE and RF analyses. A 10-fold cross-verification of LASSO was performed using the ‘glmnet’ package (version 4.1–8, cran.r-project.org/web/packages/glmnet/index.html) of R software to identify significant genes, with the minimal λ value considered optimal. The RF algorithm was executed by the ‘randomForest’ package (version 4.7–1.1, http://cran.r-project.org/web/packages/randomForest/index.html) of R software, and the top 15 genes were selected as potential candidates. The ‘e1071’ R package (version 1.7–16, http://cran.r-project.org/web/packages/e1071/index.html) was used to implement the SVM-RFE algorithm to detect the classifier with the least possible cross-validation error. The genes identified by all three methods were considered to be the signature genes for psoriasis, as shown in a Venn diagram. Finally, the area under the curve (AUC) of the receiver operating characteristic curve (ROC) was analyzed using the R package ‘pROC’ (version 1.18.5, search.r-project.org/CRAN/refmans/pROC/html/pROC-package.html) to assess the diagnostic efficacy of signature genes in the GSE30999, GSE13355 and GSE53552 datasets.

Immune infiltration analysis

To explore the distribution of immune cells in LS and NLS of patients with psoriasis in the GSE30999 dataset, the present study employed the CIBERSORT algorithm of the R software based on the ‘CIBERSORT’ R package (version 0.1.0, http://github.com/Moonerss/CIBERSORT/blob/main/R/CIBERSORT.R) (24), which contains a gene expression matrix of 22 immune cell types. The difference in infiltrated immune cells between the LS and NLS groups was analyzed using the Wilcoxon rank-sum test and visualized with a boxplot using the ‘ggpubr’ R package (version 0.6.0, cran.r-project.org/web/packages/ggpubr/index.html). P<0.05 was considered to indicate a statistically significant difference. The Spearman correlation between immune cell types was calculated and illustrated using the ‘corrplot’ tool (version 0.92, cran.r-project.org/web/packages/corrplot/index.html) in R.

For single-gene immune infiltration analysis, the Spearman correlation between the expression of signature genes and the gene expression matrix of 22 infiltrating immune cell types was calculated and displayed in a lollipop graph visualized using the ‘ggpubr’ R package (version 0.6.0, cran.r-project.org/web/packages/ggpubr/index.html). The immune cells significantly associated with signature genes (P<0.05) were extracted and scatter plots were generated.

Patients and tissue samples

The present study included 10 patients with psoriasis (4 female and 6 male patients; age, 20–49; mean age, 36.3±6.3 years) and 10 healthy controls (5 female and 5 male individuals; age, 21–50; mean age, 34.3±6.7 years). The control samples were obtained from the Plastic Surgery Department of Tianjin Academy of Traditional Chinese Medicine Affiliated Hospital (Tianjin, China). Patients with psoriasis were recruited from the Department of Dermatology of Tianjin Academy of Traditional Chinese Medicine Affiliated Hospital (Tianjin, China) between October 2022 and March 2023, and were selected based on objective criteria such as age, sex and health status to avoid self-selection bias. The inclusion criteria for patients with psoriasis were: i) Clinical presentation consistent with typical features of plaque psoriasis; ii) histological evidence of abnormal epidermal proliferation and keratinocyte differentiation (25); iii) disease duration ≥2 years, with severity during the active phase assessed by experienced dermatologists using the Psoriasis Area and Severity Index (26); and iv) no systemic treatment within 4 weeks prior to sampling and no topical treatment within 2 weeks. A total of four doctors were involved in the diagnosis of psoriasis and samples with disagreements were excluded from the study. Exclusion criteria included comorbid diabetes, renal insufficiency, history of malignancy, severe cardiovascular or cerebrovascular disease, and other autoimmune or immunodeficiency disorders. Patient samples were obtained from skin lesions on the upper arm or leg (1.0×1.0 cm), and control samples were collected from normal skin during plastic or reconstructive surgery. All samples were fixed in 4% paraformaldehyde at 4°C for 12 h and subsequently processed into 5 µm paraffin sections.

H&E staining and immunohistochemistry (IHC)

H&E staining and IHC were conducted as described previously (27). Sections were incubated with a primary antibody against TGM1 (1: 100 cat. no. 12912-3-AP; Proteintech Group, Inc.) overnight at 4°C. The slides were rinsed with PBS with 0.1% Tween-20 and treated with the secondary goat anti-rabbit HRP antibody (1:200; cat. no. HS101-01; TransGen Biotech Co., Ltd.) for 1 h at room temperature. Images were captured using a light microscope (Leica DM2000; Leica Microsystems GmbH).

TGM1 overexpression

HaCaT cells (immortalized human keratinocyte cell line; cat. no. CLS300493; Cell Line Service) were cultured in DMEM (high glucose) (cat. no. PM150210; Pricella®; Elabscience Bionovation Inc.) supplemented with 10% TransSerum® EQ Fetal Bovine Serum (cat. no. FS201-02; TransGen Biotech Co., Ltd.) and 1% penicillin-streptomycin (cat. no. FG101-01; TransGen Biotech Co., Ltd.) at 37°C. The overexpression plasmid pLV3-CMV-TGM1-human and the control plasmid pLV3-CMV-Empty-human (2 µg/ml, MiaoLing Plasmid Platform,) were transfected into HaCaT cells using Lipofectamine® 3000 reagent (cat. no. L3000015; Invitrogen; Thermo Fisher Scientific, Inc.) according to the manufacturer's protocol at 37°C for 72 h. Reverse transcription-quantitative PCR (RT-qPCR) was used to assess the mRNA expression levels of TGM1 in the HaCaT cells 72 h after transfection.

TGM1 knockdown

HaCaT cells were transfected with small interfering RNA (siRNA) targeting TGM1 (si-TGM1) using Lipofectamine® 2000 (cat. no. 11668030; Invitrogen; Thermo Fisher Scientific, Inc.) according to the manufacturer's protocol. The sequence for si-TGM1 was as follows: Sense (5′-3′), GGUGAAUAGUGACAAGGUGUAdTdT; and antisense (5′-3′), UACACCUUGUCACUAUUCACCdTdT (75 nM, Ruilai Biotechnology (Tianjin) Co., Ltd.). A non-targeting siRNA was used as the negative control [sense, 5′-UUCUCCGAACGUGUCACGUdTdT-3′ and antisense [5′-ACGUGACACGUUCGGAGAAdTdT-3′; 75 nM, Ruilai Biotechnology (Tianjin) Co., Ltd.). Following 48 h of transfection at 37°C, HaCaT cells were stimulated with a combination of recombinant human IL-17A (10 ng/ml; PeproTech, Inc.; Thermo Fisher Scientific, Inc.), IL-22 (10 ng/ml; PeproTech, Inc.; Thermo Fisher Scientific, Inc.), Oncostatin M (10 ng/ml; PeproTech, Inc.; Thermo Fisher Scientific, Inc.), TNF-α (10 ng/ml; PeproTech, Inc.; Thermo Fisher Scientific, Inc.) and IL-1α (10 ng/ml; PeproTech, Inc.; Thermo Fisher Scientific, Inc.), collectively referred to as M5. The cells were exposed to M5 for 0, 3, 6 or 12 h at 37°C and harvested for subsequent analysis.

RT-qPCR

Total RNA collected from skin tissue and HaCat cells was extracted using TransZol-Up reagent (cat. no. ET111-01; TransGen Biotech Co., Ltd.) according to the manufacturer's instructions. cDNA was synthesized using the All-in-One FirstStrand cDNA Synthesis SuperMix for RT-PCR (cat. no. AE341-02; TransGen Biotech Co., Ltd.). The RT-PCR reaction procedure is as follows: 37°C 15 min, 85°C 1 min, 4°C holding. qPCR was performed on a 7500 Fast Real-Time PCR System (Applied Biosystems; Thermo Fisher Scientific, Inc.) with PerfectStart Green qPCR SuperMix (cat. no. AQ602-24; TransGen Biotech Co., Ltd.). The RT-qPCR reaction procedure is as follows: (95°C 30 sec, 95°C 5 sec, 60°C 30 sec) for 40 cycles, 95°C 15 sec, 60°C 1 min, 95°C 15 sec. Relative mRNA expression levels were calculated using the 2−ΔΔCq method (28), with GAPDH as the internal control. The primer sequences used for RT-qPCR are listed in Table I.

Table I. Primer sequences for reverse transcription-quantitative PCR.

Table I.

Primer sequences for reverse transcription-quantitative PCR.

Gene	Forward (5′-3′)	Reverse (5′-3′)
Human-GAPDH	GGAGCGAGATCCCTCCAAAAT	GGCTGTTGTCATACTTCTCATGG
Human-IL-1α	TGGTAGTAGCAACCAACGGGA	ACTTTGATTGAGGGCGTCATTC
Human-IL-1β	CAGAAGTACCTGAGCTCGCC	GGTCGGAGATTCGTAGCTGG
Human-IL-6	CCTGAACCTTCCAAAGATGGC	TTCACCAGGCAAGTCTCCTCA
Human-IL-23	ATGTTCCCCATATCCAGTGTG	GCTCCCCTGTGAAAATATCCG
Human-S100A8	AGTGTCCTCAGTATATCA	CATCTTTATCACCAGAATG
Human-S100A9	CAACACCTTCCACCAATAC	TCATTCTTATTCTCCTTCTTGAG
Human-K1	AGAGTGGACCAACTGAAGAGT	ATTCTCTGCATTTGTCCGCTT
Human-K6	CTGGAGGCATCCAAGAGGTCA	GCAGGGTCCACTTTGTTTCCA
Human-K10	TCCTACTTGGACAAAGTTCGGG	CCCCTGATGTGAGTTGCCA
Human-K16	CAGCGAACTGGTACAGAGCA	GTTCTCCAGGGATGCTTTCA
Human-K17	GGTGGGTGGTGAGATCAATGT	CGCGGTTCAGTTCCTCTGTC
Human-TGM1	GCACCACACAGACGAGTATGA	GGTGATGCGATCAGAGGATTC

[i] S100, S100 calcium binding protein; K, keratin; TGM1, transglutaminase 1.

Western blot analysis

HaCaT Cells were lysed using RIPA lysis buffer (TransGen Biotech Co., Ltd.) containing phenylmethylsulfonyl fluoride as a protease inhibitor. Protein concentrations were quantified with the Omni-Easy™ Ready-to-Use BCA Protein Assay kit (Epizyme; Ipsen Pharma) according to the manufacturer's protocol. Equal amounts of total protein (15–20 µg/lane) were separated by 10–15% SDS-PAGE and transferred to PVDF membranes. The membranes were blocked with blocking buffer (5% skimmed milk powder in TBST) for 1 h at room temperature and subsequently incubated overnight at 4°C with primary antibodies, including anti-TGM1 (1:1,000; cat. no. 12912-3-AP; Proteintech Group, Inc.) and anti-β-actin (1:1,000, cat. no. YT0099; ImmunoWay Biotechnology Company).

Following three 10-min washes with TBS with 0.1% Tween-20, membranes were incubated with HRP-conjugated secondary antibodies (1:2,000) at room temperature for 1 h. Protein bands were then visualized with ECL reagent (cat. no. SQ201; EpiZyme, using an enhanced chemiluminescence detection system (Bioworld Technology, Inc.), and images were acquired using a Tanon 5200 chemiluminescence imaging system (Tanon Science and Technology Co., Ltd.). The gray values were analyzed using ImageJ (version 1.54, imagej.net/ij/).

Statistical analysis

Statistical analyses were performed using GraphPad Prism software 8.0 (Dotmatics). Data are presented as the mean ± SD for at least three individual experiments. The statistical significance of differences was determined using the unpaired, two-tailed Student's t-test. Other statistical analyses in the present study were performed using R software 4.3.0 (r-project.org/). P<0.05 was considered to indicate a statistically significant difference.

Results

Identification of DEGs

The ‘limma’ R package retrieved 1,845 DEGs from the GSE30999 dataset based on P<0.05 and FC>2 or FC<-2, including 1,103 upregulated and 742 downregulated genes (Fig. 1A). Most of the significant DEGs were upregulated in LS of psoriasis compared with NLS. In addition, 163 significant DEGs were selected based on |log2FC|>3 and used for subsequent analysis (Fig. 1B).

Figure 1. DEG analysis in the GSE30999 dataset. (A) Volcano plot showing the 1,845 DEGs with |FC|>2 and P<0.05. (B) Results of screening of the top DEGs with |log2FC|>3 for the following analysis shown as a heatmap. DEGs, differentially expressed genes; FC, fold change; NLS, non-lesional skin; LS, lesional skin; FDR, false discovery rate.

Functional enrichment analysis of the gene set in the GSE30999 dataset

To identify pathways enriched in the GSE30999 dataset within the context of psoriasis, GSEA was conducted. To avoid disturbance from unperturbed genes, the present study selected the top 10,000 genes from a total of 19,099 probes in the GSE30999 dataset based on descending order of |log(FC)| values for GSEA. GSEA showed that a number of upregulated pathways, including ‘Cytokine-JAK-STAT signaling pathway’, ‘DNA replication licensing’, ‘Origin unwinding and elongation’ and ‘Pre-IC formation’, and downregulated pathways, including ‘ACTH/cortisol signaling pathway’, ‘Keap1-Nrf2 signaling pathway’, ‘RTK-PLCG-ITPR signaling pathway’ and ‘TSH-TG signaling pathway’, were enriched in the GSE30999 dataset (Fig. 2).

Figure 2. Gene Set Enrichment Analysis of the top 10,000 differentially expressed genes in the GSE30999 dataset. (A) A total of four upregulated pathways were enriched in the GSE30999 dataset. (B) A total of four significant downregulated pathways were enriched in the GSE30999 dataset. (C) Ridgeline plots of the four upregulated pathways and four downregulated pathways enriched in the GSE30999 dataset. JAK, Janus kinase; ACTH, adrenocorticotropic hormone; Keap1, Kelch-like ECH-associated protein 1; Nrf2, nuclear factor erythroid 2-related factor 2; PLCG, phospholipase Cγ; ITPR, inositol 1,4,5-trisphosphate receptor; TSH, thyroid-stimulating hormone; TG, thyroglobulin. RTK, receptor tyrosine kinase; Pre-IC: Pre-Initiation Complex.

Functional annotation of 163 DEGs

The present study used GO enrichment analysis to obtain a more precise functional overview of the top 163 DEGs. A total of 284 BP terms, 17 CC terms and 50 MF terms were identified, and the top 10 GO terms of each category are shown in Fig. 3. PBH was obtained via Benjamini and Hochberg (BH) correction. The top 10 GO BP terms were ‘antimicrobial humoral response’ (PBH=2.11×10−6), ‘antimicrobial humoral immune response mediated by antimicrobial peptide’ (PBH=8.87×10−6), ‘response to virus’ (PBH=1.04×10−5), ‘defense response to virus’(PBH=1.04×10−5), ‘negative regulation of viral genome replication’ (PBH=1.22×10−5), ‘type I interferon signaling pathway’ (PBH=2.07×10−5), ‘cellular response to type I interferon’ (PBH=2.07×10−5), ‘response to type I interferon’ (PBH=2.52×10−5), ‘regulation of viral genome replication’ (PBH=2.52×10−5) and ‘neutrophil chemotaxis’ (PBH=2.94×10−5; Fig. 3A and D). The most significant GO CC and MF terms were ‘specific granule lumen’ (PBH=1.29×10−4) and ‘chemokine activity’ (PBH=5.97×10−5), respectively (Fig. 3B, C, E and -F).

Figure 3. Functional analysis of the top 163 differentially expressed genes. (A) Top 10 GO biological process terms. (B) Top 10 GO cellular component terms. (C) Top 10 GO molecular function terms. (D) Cnetplot of the top 10 GO biological process terms. (E) Cnetplot of the top 10 GO cellular component terms. (F) Cnetplot of the top 10 GO molecular function terms. (G) Dotplot of the top 10 KEGG pathways. (H) Cnetplot of KEGG pathways. (I) Emapplot of interacting KEGG pathways. KEGG, Kyoto Encyclopedia of Genes and Genomes pathways. GO, Gene Ontology.

A total of 21 pathways were identified using KEGG analysis and the top 10 pathways are shown in Fig. 3. These were ‘Viral protein interaction with cytokine and cytokine receptor’ (PBH=9.67×10−5), ‘IL-17 signaling pathway’ (PBH=3.53×10−4), ‘Influenza A’ (PBH=1.67×10−2), ‘Chemokine signaling pathway’ (PBH=2.41×10−2), ‘Cytokine-cytokine receptor interaction’ (PBH=2.41×10−2), ‘NOD-like receptor signaling pathway’ (PBH=7.56×10−2), ‘Epithelial cell signaling in Helicobacter pylori infection’ (PBH=1.29×10−1), ‘Coronavirus disease-COVID-19’ (PBH=1.97×10−1), ‘Transcriptional misregulation in cancer’ (PBH=2.31×10−1) and ‘Glycosaminoglycan degradation’ (PBH=2.53×10−1) (Fig. 3G-I).

Selection of signature DEGs using machine learning techniques

The present study utilized three machine learning methods (LASSO, RF and SVM-RFE) to extract the possible diagnostic DEGs from a total of 163 DEGs. The LASSO regression algorithm identified 13 core genes with optimal λ values (λ.min, 0.01264872) in the GSE30999 dataset (Fig. 4A and B), which were SYNCRIP, VNN1, GBP1, TGM1, TMPRSS11D, CYP2C18, PRKCQ, RSAD2, FUT3, HAO2, NLRP2, TMEM86A and LYPD5. Using a RF algorithm, 35 genes with an importance score >1.0 were identified, and the top 15 in importance were chosen as signature DEGs (Fig. 4C and D). Through SVM-RFE analysis, a total of 29 DEGs were identified with an accuracy of 0.975 and an error rate of 0.0252 (Fig. 4E and F). Ultimately, these prediction methods identified TGM1. Therefore, TGM1 was employed as a diagnostic marker of psoriasis in later investigations (Fig. 4G).

Figure 4. Identification of biomarkers of psoriasis based on machine learning algorithms. (A) Optimal penalization coefficient λ in the LASSO model and (B) average deviance values for each model with a given λ. (C) Association between the total number of trees in the RF algorithm and the error rates, and (D) the top 15 signature genes. (E) Error rate of the curve after 5-fold cross-validation of SVM-RFE algorithm to select the signature genes. Error rate predicted by the 29 genes is 0.0252. (F) Accuracy rate predicted by 29 genes was 0.975. (G) Venn diagram of the biomarkers identified using LASSO, SVM-RFE and RF algorithms. RF, random forest; SVM-RFE, support vector machine-recursive feature elimination; LASSO, least absolute shrinkage and selection operator; TGM1, transglutaminase 1, CV, cross validation.

Identification of diagnostic markers in psoriasis

PPIs are important in numerous biological pathways. The majority of proteins serve their functions through interactions with a large number of other proteins. Therefore, the 163 identified DEGs were submitted to the STRING database to acquire their interaction data with TGM1. Of the 163 DEGs, a total of 52 genes were disconnected, and thus, removed from the PPI network. After removing the 52 disconnected genes, the final PPI network contained 111 nodes and 388 edges. Among the 111 connected proteins, CNFN, SPRR3, SDR9C7, PI3, KRT16, DSG3, PKP1 and KLK6 were directly connected with TGM1 (Fig. 5A). Additionally, in comparison with the NLS tissue, the LS tissues of psoriasis exhibited a significant increase in TGM1 expression (Fig. 5B). ROC curves highlighted the robust diagnostic potential of TGM1 as biomarkers for psoriasis (Fig. 5C).

Figure 5. Identification of biomarkers of psoriasis. (A) TGM1 and its connected genes in the protein-protein interaction network constructed using 163 differentially expressed genes with |log2(FC)|>3. (B) Box plot of TGM1 expression in LS and NLS tissues of patients with psoriasis. (C) A receiver operating characteristic curve was employed to assess the diagnosis relevance of TGM1 in psoriasis. LS, lesional skin; NLS, non-lesional skin; AUC, area under the curve; TGM1, transglutaminase 1.

Identification of potential pathways of TGM1 by single-gene GSEA

Single-gene GSEA is a bioinformatics tool used to investigate the biological pathways enriched by genes associated with a specific gene, and thus, to reveal its particular biological activity. Single-gene GSEA of TGM1 was conducted on LS samples, and revealed that TGM1 was primarily enriched in upregulated pathways, including the ‘ARNO-ARF-ACTB_G signaling pathway’ [false discovery rate (FDR), 0.0062], ‘Cytokine-JAK-STAT signaling pathway’ (FDR, 0.0140), ‘Tight junction-Actin signaling pathway’ (FDR, 0.0320), ‘TLR2/4-MAPK signaling pathway’ (FDR, 0.0064) and ‘TNF-NFκB signaling pathway’ (FDR, 0.0325) (Fig. 6A). Downregulated pathways included the ‘EP/NE-ADRB-cAMP signaling pathway’ (FDR, 1.71×10−3), ‘GF-RTK-PI3K signaling pathway’ (FDR, 9.21×10−3), ‘RTK-PLCG-ITPR signaling pathway’ (FDR, 1.57×10−3), ‘Translation initiation’ (FDR, 2.29×10−5) and ‘Wnt signaling modulation, LGR/RSPO’ (FDR, 8.13×10−3) (Fig. 6B).

Figure 6. Functional enrichment for TGM1 via single-gene Gene Set Enrichment Analysis. (A) Up- and (B) Downregulated signaling pathways involving TGM1. TGM1, transglutaminase 1.

Validation of signature genes

The present study first examined TGM1 expression in the GSE53552 and GSE13355 datasets. The results showed that TGM1 expression was significantly increased in LS tissues of psoriasis compared with the NLS tissues in the GSE53552 (Fig. 7A and B) and GSE13355 (Fig. 7D and E) datasets. Additionally, the diagnostic efficacy of TGM1 was also validated using the GSE53552 and GSE13355 datasets, which showed a high associated value with AUC values of 1.000 and 0.893, respectively (Fig. 7C and F). Finally, IHC was performed to evaluate TGM1 protein expression levels in affected areas of patients with psoriasis. IHC showed that TGM1 expression exceeded normal levels in psoriatic lesions (Fig. 8).

Figure 7. Validation of TGM1 expression in validation datasets. Volcano plots showing the differentially expressed genes in the (A) GSE53552 and (D) GSE13355 datasets based on |FC|>2 and P<0.05. The dark red dot indicates TGM1. Boxplots showing TGM1 expression in the LS and NLS of patients with psoriasis in the (B) GSE53552 and (E) GSE13355 datasets. Receiver operating characteristic curves of TGM1 in the (C) GSE53552 and (F) GSE13355 datasets. NLS, non-lesional skin; LS, lesional skin; FC, fold change; AUC, area under the curve; TGM1, transglutaminase 1.

Figure 8. TGM1 expression in psoriatic skin validated by IHC. (A) H&E staining of skin tissues of healthy controls and patients with psoriasis. (B) IHC of TGM1 protein in skin tissues of healthy controls and patients with psoriasis. (C) TGM1-positive cells in skin tissues of healthy controls and patients with psoriasis. ***P<0.001 vs. healthy control. IHC, immunohistochemistry; TGM1, transglutaminase 1.

Immune cell infiltration in psoriasis

The immune cell infiltration results based on GSE30999 showed that plasma, naïve CD4 T cells, activated memory CD4 T cells, follicular T helper cells, γδT cells, activated NK cells, M1 macrophages, M2 macrophages, resting dendritic cells, activated dendritic cells, resting Mast cells, eosinophils and neutrophils were significantly increased in the LS group compared with the NLS group, while activated natural killer cells and resting mast cells were significantly decreased (Fig. 9A). The correlation study revealed a significant positive correlation between activated dendritic and follicular helper T cells as well as between memory B cells and naïve CD4 T cells (Fig. 9B). A significant negative correlation was observed between activated dendritic cells and resting mast cells as well as between follicular T helper cells and resting mast cells (Fig. 9B). In order to further validate the association between TGM1 expression and the immune component, the present study employed single-gene immune infiltration analysis. The results showed association between TGM1 expression and eosinophils, activated dendritic cells and follicular T-helper cells (Fig. 9C-F). There was an association between TGM1 expression and resting mast cells (Fig. 9C and G).

Figure 9. Correlation of TGM1 expression with infiltrating immune cells in psoriasis. (A) Expression of immune cells between LS and NLS. (B) Correlation between 22 types of infiltrating immune cells between LS and NLS. (C) Lollipop graph showing the correlation between TGM1 expression and 22 types of infiltrating immune cells in psoriasis. Scatter plots showing that (D) eosinophils, (E) activated dendritic cells, (F) follicular T helper cells and (G) resting mast cells were significantly correlated with TGM1 expression. *P<0.05, **P<0.01, ***P<0.001 and ****P<0.0001 vs. NLS. NLS, non-lesional skin; ns, not significant; TGM1, transglutaminase 1; NK, natural killer; abs, absolute; cor, correlation coefficient.

TGM1 aggravates the inflammatory response and keratinocyte differentiation

To validate the prediction results of machine learning techniques, transfection of HaCaT cells was performed to support the effect of TGM1. The present study first performed an overexpression experiment. The mRNA and protein expression levels of TGM1 were significantly increased after transfection with the overexpression plasmid compared with those in the negative control group (Fig. 10A and B). Additionally, mRNA expression levels of some inflammatory cytokines, such as IL-1α, IL-1β, IL-6 and IL-23, were significantly increased after TGM1 overexpression (Fig. 10C). Finally, the expression levels of several markers of keratinocyte differentiation were also increased following TGM1 overexpression, of which S100 calcium binding protein (S100)A8, keratin (K)6, K10 and K17 showed significant increases in expression compared with the control group (Fig. 10D and E).

Figure 10. mRNA expression levels of inflammatory cytokines and differentiation markers of keratinocytes after TGM1 overexpression. TGM1 (A) mRNA and (B) protein expression following sh-TGM1 treatment. (C) mRNA expression levels of inflammatory factors IL-1α, IL-1β, IL-6 and IL-23. mRNA expression levels of (D) S100A8 and S100A9, as well as (E) K1, K6, K10, K16 and K17 after sh-TGM1 treatment. *P<0.05, **P<0.01, ***P<0.001 and ****P<0.0001. Neg, negative control; S100, S100 calcium binding protein; K, keratin; TGM1, transglutaminase 1.

To further verify the effect of TGM1, knockdown experiments were subsequently performed using siRNA. As shown in Fig. 11A and B, TGM1 mRNA and protein expression levels were significantly reduced by si-TGM1 transfection compared with those in the negative control group. mRNA expression levels of inflammatory cytokines, such as IL-1β, was downregulated following TGM1 downregulation at 12 h (Fig. 11C). Similarly, the expression levels of the keratinocyte differentiation markers S100A8, S100A9, K1 were also decreased at 12 h (Fig. 11D and E).

Figure 11. mRNA expression levels of inflammatory cytokines and differentiation markers of keratinocytes after TGM1 knockdown. TGM1 (A) mRNA and (B) protein expression after si-TGM1 transfection. (C) mRNA expression levels of inflammatory factors IL-1α, IL-1β and IL-6. mRNA expression levels of (D) S100A8 and S100A9, as well as (E) K1, K6, K10, K16 and K17 following M5 stimulation after si-TGM1 transfection. *P<0.05, **P<0.01, ***P<0.001 and ****P<0.0001 vs. si-Neg. si, small interfering RNA; Neg, negative control; TGM1, transglutaminase 1; S100, S100 calcium binding protein; K, keratin.

Discussion

Psoriasis is a notable public health issue with complex pathophysiology, with clinically relevant biomarkers (4). Despite advancements, the etiological mechanisms of psoriasis remain yet to be fully elucidated, necessitating systematic exploration of disease pathways to establish novel theoretical frameworks for early diagnosis and targeted interventions.

The present study examined psoriasis-related microarray data from the GEO database. The present study initially performed GSEA on the GSE30999 dataset, uncovering significant enrichment of the ‘Cytokine-JAK-STAT signaling pathway’ and ‘Keap1-Nrf2 signaling pathway’. The JAK/STAT pathway, a recognized regulator of inflammatory and autoimmune disorders (29), facilitates the transcriptional activation of psoriasis-associated cytokines such as IL-17 and IL-23, driving keratinocyte hyperproliferation and disease symptoms (25). Clinical evidence additionally indicates that pharmacological inhibition of this pathway reduces cutaneous cytokine production (30,31). The Keap1/Nrf2 axis, an important regulator of redox homeostasis and cutaneous barrier integrity (32), exhibits therapeutic significance in psoriasis, with Nrf2 activators proving effective in moderate-to-severe cases (33,34). The enhancement of these pathways supports the biological validity of the results of the present study.

Subsequent differential expression analysis identified 163 significant DEGs with |log2FC|>3 from 1,845 candidates in the GSE30999 dataset. GO and KEGG functional enrichment analyses identified two important pathways (‘Viral protein interaction with cytokine and cytokine receptor’ and ‘IL-17 signaling pathway’). The etiology of psoriasis entails T cell-mediated cytokine cascades, including cytokines such as TNF-α, IL-17A, IL-23 and IL-36, which dysregulate keratinocyte dynamics via receptor-mediated signaling (35). IL-17 receptor activation markedly stimulates STAT3-dependent keratinocyte proliferation, intensifying inflammatory responses (36). The IL-17 family, mostly produced by T helper 17 cells, displays isoform-specific functions in chronic inflammation and psoriasis progression, with IL-17A exhibiting the most notable pathophysiological impact (37–40). Targeted biologics, such as secukinumab and ixekizumab, effectively neutralize IL-17A, thus modulating cytokine networks and enhancing clinical outcomes in moderate-to-severe psoriasis (41).

The present study utilized three machine learning methods, LASSO, RF and SVM-RFE, to identify diagnostic biomarkers, ultimately identifying TGM1 as a possible biomarker for psoriasis. TGM1, a calcium-dependent enzyme in the epidermal differentiation complex, catalyzes ε-(g-glutamyl) lysine crosslinks during cornified envelope formation, a key process for epidermal barrier function (42). Barrier dysfunction, a characteristic of psoriasis and related dermatoses (43,44), results from cytokine-driven keratinocyte hyperproliferation (45) and differentiation imbalances associated with disease severity (46). TGM1 mutations are associated with laminar ichthyosis (47). However, the upregulation of TGM1 in psoriatic lesions indicates a poorly understood pathogenic function (48). TGM1 exhibits diagnostic utility in multiple malignancies, such as adrenocortical carcinoma and bladder cancer (49,50), underscoring its potential as a dual-purpose biomarker in psoriasis.

Single-gene GSEA and immune infiltration analyses further clarified the molecular roles of TGM1. Upregulated pathways, including the ‘cytokine-JAK-STAT signaling pathway’, ‘TLR2/4-MAPK signaling pathway’ and ‘TNF-NFκB signaling pathway’, and downregulated pathways, such as ‘GF-RTK-PI3K signaling pathway’ and ‘Wnt signaling modulation, LGR/RSPO’, correspond with established psoriatic pathobiology (5). Immune profiling revealed positive correlations between TGM1 expression and eosinophils, activated dendritic cells and follicular T-helper cells, and a negative correlation with resting mast cells, suggesting that TGM1 serves a role in innate immune dysregulation. Cross-dataset validation using both the GSE53552 (AUC, 1.000) and GSE13355 (AUC, 0.893) datasets supported TGM1 upregulation, further validated by IHC revealing increased protein levels in psoriatic lesions compared with controls. Finally, cell experiments were conducted to verify the potential effects of TGM1 on psoriasis and revealed that TGM1 may aggravate the inflammation response and keratinocyte differentiation. Collectively, these findings established TGM1 as an important mediator and therapeutic target in psoriasis.

The present multimodal investigation integrating bioinformatics, machine learning and experimental validation identified TGM1 as a functionally significant biomarker in the pathogenesis of psoriasis. The observed upregulation in affected skin, pathway linkage and immunological correlations offer mechanistic insights into disease progression. The findings of the present study enhance the understanding of psoriatic pathophysiology and suggest TGM1-targeted strategies for diagnostic and therapeutic innovation.

Acknowledgements

The authors would like to thank Professor Litao Zhang and Dr Lin Li (Tianjin Academy of Traditional Chinese Medicine Affiliated Hospital, Tianjin, China) for providing the wax blocks of the patients with psoriasis.

Funding

The present work was supported by The Science and Technology Development Fund of Tianjin Education Commission for Higher Education (grant no. 2022ZD047).

Availability of data and materials

The data generated in the present study may be requested from the corresponding author.

Authors' contributions

PG and MS contributed to the literature search and study design. JZ, QY, HC and JH analyzed and interpreted data. JL and QZ assisted with the experiments and wrote the manuscript. BZ was responsible for project administration and secured funding for the study. BZ and LH conceived and designed the experiments, revised the manuscript and confirmed the authenticity of all the raw data. All authors read and approved the final version of the manuscript.

Ethics approval and consent to participate

For patient samples, written informed consent was obtained from each patient and the study was approved by the Ethics Committee of Tianjin Academy of Traditional Chinese Medicine Affiliated Hospital (approval no. LLKY2022-39; Tianjin, China). The study was performed in accordance with The Declaration of Helsinki.

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

References