Identification and Validation of Hub Genes in Hidradenitis Suppurativa

Introduction

Hidradenitis suppurativa (HS), also known as acne inversa, is a recurrent, chronic inflammatory skin disease, which manifests as pimples and inflammatory nodules that progress to abscesses, sinus tracts, and hypertrophic scarring.¹ HS usually begins in adolescence and develops in friction-prone areas, such as the armpits, groin, perianal areas, and buttocks. Characterized by pain, malodor, recurrent attacks, and difficult healing, HS seriously affects patients’ social interactions and work, thus reducing their quality of life.²

The risk factors contributing to HS are primarily genetic and environmental, and are associated with obesity, smoking, and mechanical friction.³ From a pathological perspective, HS is believed to result from hyperkeratosis of follicular openings accompanied by blockage, with consequent damage to the follicular epithelium and the production of a bacterial biofilm. Hair follicle rupture overstimulates the innate immune system, leading to secondary inflammation. The molecular mechanisms underlying the pathogenesis of HS are complex. Recent studies have focused on the role of epigenetic dysregulation, particularly DNA methylation, histone modifications, and non-coding RNAs such as microRNAs (miRNAs) in the pathogenesis of HS.^4,5 The epigenetic alterations give insight into disease pathology and potential therapeutic avenues.⁶ A study by U. Radhakrishna identified 60 CpG sites associated with HS, including 54 sites with low methylation and 6 sites with high methylation.⁷ The following year, 170 cytokine genes related to HS were identified, including 27 highly methylated CpG sites and 143 low-methylated sites. The pathways enriched by these 170 genes were also identified.⁸ A study introduced a novel epigenetic mechanism, showing that specific gene methylation signatures are closely related to pain perception and sensitivity in HS patients. The differentially methylated genes involved in chronic pain in HS may serve as biomarkers and therapeutic targets.⁹ While studies have explored genes and pathways involved in the pathogenesis of HS, most only elucidated a specific aspect of the pathogenetic process. A comprehensive understanding of HS pathogenesis thus requires further study.

In the North American Clinical Management Guidelines for HS, the anti-tumor necrosis factor (TNF)-α antibody adalimumab is the only recommended treatment. This recommendation was based on consistent and high-quality patient-oriented evidence.¹⁰ TNF-α is a key driver of immune response.¹¹ Other therapeutic agents include systemic and topical antibiotics, anti-androgen contraceptives, retinoids, and systemic immunomodulators. However, rigorous and strong evidence for drug efficacy is currently lacking and some even have the potential to exacerbate HS in certain cases (eg, progestogen-only dosing regimens).¹⁰ Taken together, physicians are confronted with a paucity of therapeutic options, and the clinical management of moderate-to-severe HS remains a significant challenge.

Bioinformatic analysis, which has been widely used in the study of multiple diseases,^12,13 can process large amounts of samples within short time and provide valuable information about diseases. The approach facilitates the identification of the gene regulatory networks, hub genes and their associated transcription factors (TFs), microRNAs (miRNAs), and drugs in various diseases. Bioinformatic analysis can be influenced by disease context, tissue type, experimental design and computation approach. Therefore, the identification of hub genes is not universal.

To gain further insight into the pathogenesis of HS, we conducted a bioinformatics analysis of the disease with the aim of identifying disease-associated genes and their corresponding pathways, using cytoHubba plugin and MCODE plugin. The diagnostic significance of hub genes was analyzed, and the network relationships between mRNAs and TFs, miRNAs, and mRNAs, as well as genes and drugs were predicted. Furthermore, we analyzed the immune-infiltrating cells associated with HS and examined the expression levels of relevant genes in the skin tissues of patients and healthy controls. Our findings might contribute to a deeper understanding of the pathogenesis of HS by highlighting specific implicated genes and pathways.

Methods

Data Acquisition and Processing

We obtained the GSE137141 and GSE148027 datasets from the GEO database (https://www.ncbi.nlm.nih.gov/geo/)¹⁴ using the R package “GEOquery”.¹⁵ The dataset GSE137141 was obtained from the data platform GPL16699 and screened to obtain eight samples of HS lesion sites in patients diagnosed with Hurley stage II–III HS, and eight samples of normal skin in the healthy group. The dataset GSE148027 was obtained from the data platform GPL570 and screened to obtain 18 samples of HS lesion sites and eight samples of normal skin in the healthy group. In the course of the data annotation process, it was decided to remove one probe corresponding to multiple molecules to avoid an effect on subsequent data processing. When there were multiple probes for the same molecule, only the probe with the largest signal value was retained. Finally, we normalized the data using the “normalizeBetweenArrays” function of the R package “limma”,¹⁶ and visualized the results deploying the “ggplot2” package.

Identification of Differentially Expressed Genes (DEGs)

We used the “limma” package to perform a differential gene expression analysis of the control and HS groups in the GSE137141 and GSE148027 datasets, respectively. The absolute values of logFC > 1 and p.adj < 0.05 were used as thresholds. LogFC > 1 indicated the upregulation of genes, and logFC < −1 indicated the downregulation of genes. Box plots, principal component analysis (PCA) plots, volcano plots, and heat maps were constructed using the “ggplot2” package. The DEGs obtained from each dataset were used as intersections, and Venn diagrams were plotted using the “VennDiagram” package for shared differential genes.

GO and KEGG Enrichment Analyses

The GO knowledge base (https://www.geneontology.org/)¹⁷ is a comprehensive resource for functional genomics. The KEGG database (https://www.kegg.jp)¹⁸ contains information on biochemical pathways. GO annotation analysis and KEGG pathway enrichment analysis of DEGs were performed using the R package “org.Hs.eg.db” and p.adj < 0.05 indicated a statistically significant difference. The “ggplot2” package was employed for visualization.

PPI Network

The STRING database (https://cn.string-db.org/)¹⁹ predicts PPI networks to analyze the functional interactions of proteins and understand the mechanisms of disease development and progression. STRING was used to form a PPI network of DEGs, with a coefficient of 0.4. Cytoscape^20,21 was used to visualize the PPI network.

Identification of Hub Genes and Protein Modules

CytoHubba,²² a Cytoscape plugin, filters genes using different topological analysis methods. We chose the first 10 nodes according to MCC. Significant modules (clusters of genes) were identified using the MCODE plugin with the following filtering criteria: degree cutoff = 2, node score cutoff = 0.2, k-core = 2, and max depth = 100. The top three modules were selected based on their scores.

Receiver Operating Characteristic (ROC) Curves and Area Under the Curve (AUC) Values

ROC curves²³ can be used to evaluate the performance of a binary diagnostic classification method and compare the performance of two or more diagnostic tests to assess whether the identified genes are useful in diagnosing HS. ROC curves were plotted using the “pROC” (ver. 1.18.0) package in R, and the AUC values were calculated. AUC is a quantitative measure of the accuracy of a diagnostic test, and an AUC of >0.9 usually indicates good diagnostic significance.

mRNA–miRNA and mRNA–TF Networks

Next, a prediction of the target miRNAs and TFs of the identified hub genes was conducted. miRNAs were predicted using the miRDB²⁴ and miRWalk²⁵ databases, and intersections were used to identify target miRNAs. Similarly, TFs were predicted using the KnockTF and ChIPBase databases, and those at the intersection were considered as our target TFs. We generated mRNA–miRNA and mRNA–TF interaction networks with the help of Cytoscape.

Target Gene–Drug Networks

The Drug–Gene Interaction Database (DGIdb),²⁶ a web-based resource, houses data pertaining to drugs, genes, and the interactions between them. To provide insight into the relationships between drugs and genes, we predicted potential drugs corresponding to hub genes using DGIdb (filter setting: Regulatory Approval, Approved) and constructed a gene–drug network. Cytoscape was used to visualize the network.

Analysis of Immune Cell Infiltration

The R package “CIBERSORT”²⁷ was employed for immune cell infiltration analysis. We first obtained the built-in data “LM22.txt” in the R package, which recorded the gene expression data of 22 immune cells. Then the gene expression data matrix was prepared, and the data was filtered (p < 0.05) to obtain the immune cell infiltration abundance matrix. Heatmaps of 22 immune cell types in HS samples were generated using the R packages “pplog2” and “LinkET.”

Quantitative Real-Time PCR (qRT-PCR)

From 1 July 2023 to 1 July 2024, we collected 10 healthy control samples and 10 HS patient lesional biopsies from hospital. All samples were from East Asians aged 18 to 60. Specifically, the samples pertaining to patients with HS were collected from individuals who had been clinically diagnosed with the condition and exhibited disease severity classified as levels II to III. The experimental protocol received approval from the Ethics Committee of the Affiliated Zhangjiagang Hospital of Soochow University on 1 July 2023. All participants provided informed consent prior to the study. We collected and preserved skin samples using containers filled with formalin and extracted total RNA from the samples. Reverse transcription was performed using the Applied Biosystems 2720 thermal cycler (Thermo Fisher Scientific, Waltham, MA, USA). A real-time polymerase chain reaction (PCR) instrument (Bio-Rad, Hercules, CA, USA) was used to perform quantitative PCR using MonScript™ ChemoHS qPCR Mix (Monad, Guangzhou, China) according to the manufacturer’s instructions. GAPDH was used as an internal reference for data standardization. Relative expression level was calculated using the 2-ΔΔCt method.

Statistical Analysis

All statistical analyses and visualizations were executed using R software (http://www.r-cret.org; version 4.2.1). Differences were considered statistically significant at p < 0.05. Regarding the comparison of continuous variables between the two groups, the statistical significance of normally distributed variables was evaluated using the independent Student’s t-test, whereas differences between non-normally distributed variables were estimated using the Mann–Whitney U-test (ie, Wilcoxon rank-sum test). All statistical tests were two-sided.

Results

Identification of DEGs in HS

A flowchart illustrating the study’s methodology is presented in Figure 1a. Corrections on the datasets GSE137141 and GSE148207 were performed using the “fdr” method²⁸ to control the false rate (Figure 1b–e). In the GSE137141 dataset, a total of 349 DEGs were identified, consisting of 264 upregulated and 85 downregulated genes (Figure 2a and c, Table S1). In the GSE148027 dataset, 2,303 DEGs were identified, consisting of 1,285 upregulated and 1,018 downregulated genes (Figure 2b and d, Table S2). Volcano plots and heat maps were generated to visualize the results. By separating the intersection of the upregulated and downregulated genes in the two datasets, we obtained 138 upregulated and 42 downregulated genes and presented them in Venn diagrams (Figure 3a and b).

Figure 1 Flowchart of the study, box diagram for data correction, and sample cluster analysis. (a). Flowchart of the study. (b and c).Distribution of gene expression after correction for GSE137141 (b) and GSE148027 (c). Gene expression was relatively consistent between datasets, verifying the validity of normalization. (d and e) PCA plots of group clustering after correction for GSE137141 (d) and GSE148027 (e) data. The results revealed significant differences in gene expression between the control and HS groups.

Figure 2 Differential gene expression analysis. (a and b) Volcano plots of the results of differential gene analysis for GSE137141 (a) and GSE148027 (b). (c and d) Heatmaps of the degree of gene expression for GSE137141 (c) and GSE148027 (d), showing significant differentially expressed genes. Red squares represent upregulated genes, and blue squares represent downregulated genes.

Figure 3 Venn diagrams of upregulated and downregulated genes and bubble and circle plots of GO and KEGG enrichment analysis for the GSE137141 and GSE148027 databases. (a and b) Venn diagrams of upregulated genes (a) and downregulated genes (b) in the GSE137141 and GSE148027 datasets. (c). Bubblediagrams of GO and KEGG enrichment analyses for DEGs, where larger bubbles represent more significant enrichment. (d and e) Circle plots of GO and KEGG enrichment analyses for GSE137141 (d) and GSE148027 (e) datasets.

GO and KEGG Enrichment Analyses

GO annotation analysis was conducted on the DEGs, and the results were visualized in a bubble map, revealing features that were significantly enriched in different categories (Figure 3c). The five most significant terms based on p. adj values were selected for visualization. Similarly, KEGG pathway enrichment analysis was performed, and the top seven pathways are shown according to p. adj values, from the smallest to largest. In the biological process (BP) category, the DEGs were primarily enriched in extracellular matrix organization, extracellular structure organization, and external encapsulating structure organization. In the cellular component (CC) category, DEGs were primarily enriched in external side of plasma membrane and collagen containing extracellular matrix. In the molecular function (MF) category, the enriched terms mainly included serine type endopeptidase activity, serine type peptidase, and serine hydrolase. Enriched KEGG pathways included primary immunodeficiency, chemokine signaling, leukocyte transendothelial migration, IL−17 signaling pathway, PI3K−Akt signaling, cytokine receptor interaction, and cell adhesion molecules. The complete list of enriched GO terms and KEGG pathways were supplied in Table S3. GO and KEGG enrichment results are displayed in Circos maps (Figure 3d and e).

PPI Network and Identification of Hub Target Genes

To evaluate interactions among the 180 DEGs, we constructed a PPI network of DEGs using STRING (Figure 4a). The PPI network comprised 177 nodes and 643 edges. With the implementation of the MCC algorithm of the CytoHubba plugin, 10 hub genes were identified from the PPI network: PTPRC, CD19, IL7R, SELL, ITGAL, CD27, ITGAX, FCGR2A, CD2, and IL2RG (Figure 4b).

Figure 4 PPI network of DEGs (a) and the top 10 hub genes closely related to HS (b).

Protein Module Analysis

Figure 5 illustrates the top three protein modules in terms of significance, with the first significant module containing the 10 hub genes previously identified.

Figure 5 Mcode protein module analysis. (a). Most significant protein-interaction module. Blue nodes indicate hub genes. (b). Second significant protein interaction module. (c). Third significant protein-interaction module.

Diagnostic Significance of Hub Genes

We plotted boxplots of the differential expression of hub genes in the GSE137141 and GSE148027 datasets to determine the expression of hub genes in the control and HS groups (Figure 6a and b). Statistical significance was set at p < 0.05. In both the GSE137141 and GSE148027 datasets, 10 upregulated genes were identified, and no downregulated genes were identified. In other words, all 10 hub genes (PTPRC, CD19, IL7R, SELL, ITGAL, CD27, ITGAX, FCGR2A, CD2, and IL2RG) were upregulated in the HS Group.

Figure 6 Diagnostic significance of hub genes in HS. (a and b) Boxplots of differential expression of hub genes in GSE137141 (a) and GSE148027 (b). (c–f). ROC curves for CD19 and CD2 (c), CD27 and FCGR2A (d), IL2RG and IL7R (e), PTPRC and SELL (f), and ITGAL and ITGAX (g) in the GSE137141 dataset. (h–l). ROC curves for CD19 and CD2 (h), CD27 and FCGR2A (i), IL2RG and IL7R (j), PTPRC and SELL (k), and ITGAL and ITGAX (l) in the GSE148027 dataset.

ROC curves demonstrated the diagnostic significance of hub genes (Figure 6c–l). The closer the curve is to the upper left corner coordinate (0, 1), the greater its diagnostic significance. In the GSE137141 dataset, all hub genes showed good diagnostic significance.

mRNA–TF Network

The KnockTF and ChIPBase databases were used to predict the TFs associated with hub genes in HS (Table S4). Based on TFs intersecting in the Venn diagram of both databases,105 pairs of mRNA–TFs were identified (Figure 7a). We used Cytoscape to produce an mRNA–TF network, in which mRNA–TFs contained nine mRNAs and 69 targeting TFs (Figure 7b).

Figure 7 Prediction and interaction network of target TFs, miRNAs and drugs of hub genes. (a). Venn diagram of TFs obtained from the KnockTF and ChIPBase databases. (b). Networks of hub genes and target TFs. (c). Venn diagram of miRNAs predicted using the miRDB and miRWalk databases. (d). Network of hub genes and targeted drugs from the DGIdb. (e). Network of hub genes and targeted drugs. Blue oval dots represent hub genes, green diamonds represent TFs targeted by hub genes, yellow rectangles represent miRNAs targeted by hub genes, and purple V-shapes represent drugs related to hub genes.

mRNA–miRNA Networks

We used the miRDB and miRWalk databases to predict the target miRNAs of hub genes, and used their intersections as target miRNAs (Figure 7c, Table S5). An mRNA–miRNA interaction network was constructed using Cytoscape (Figure 7d). The network showed that there were 115 miRNAs closely related to 10 hub genes.

Target Gene–Drug Network

We used the DGIdb database to predict drug components that might act on hub genes, and constructed a gene–drug interaction network using Cytoscape. A total of 46 gene–drug pairs were obtained, including eight genes (PTPRC, FCGR2A, IL2RG, ITGAL, CD2, CD19, IL7R, and CD27) and 41 drugs. Among the 10 hub genes, two (SELL and ITGAX) were not predicted to interact with the drugs (Figure 7e).

Immune Cell Infiltration Environment of Hub Genes

We generated heat maps of the correlation between the 22 immune cell types in the GSE13714 and GSE148027 datasets. In the GSE13714 dataset, follicular helper T cells, resting mast cells (r = −0.823707973490796, p = 8.78114794442553e-05), and M2 macrophages (r = −0.673875560341157, p = 0.00420580713529972) were negatively correlated with HS. Activated dendritic cells, memory CD4 T cells (r = 0.751279939228315, p = 0.0007937865770693), regulatory T cells, and neutrophils (r = 0.684991268834017, p = 0.00341137629884325) were positively correlated with the disease (Figure 8a).

Figure 8 Immune cell infiltration in GSE137141 (a) and GSE148027 (b).

In the GSE148027 dataset, resting mast cells and plasma cells (r = −0.721025641025641, p = 3.24526672508263e-05), neutrophils (r = −0.716239316239316, p = 3.87178808889388e-05), and M1 macrophages (r = −0.666202114948614, p = 0.000202882678089322) were negatively correlated with HS. Meanwhile, resting dendritic cells, regulatory T cells (r = 0.6824238626565, p = 0.000122782662943521), activated mast cells, neutrophils (r = 0.681025687196947, p = 0.000128370501825792), and activated memory CD4 T cells (r = 0.670445387322356, p = 0.000178428002754583) were positively correlated (Figure 8b).

Validation of Hub Gene Expression in HS

The expression levels of PTPRC, CD19, IL7R, SELL, ITGAL, CD27, ITGAX, FCGR2A, CD2, and IL2RG in 10 HS and 10 normal skin tissues were examined using qRT-PCR. Compared with that in the control samples, the expression of FCGR2A and IL2RG was significantly upregulated in the HS samples (p < 0.05). PTPRC, CD19, IL7R, SELL, ITGAL, CD27, ITGAX, and CD2 expression was not significantly different between HS and control samples (Figure 9a–j).

Figure 9 qRT-PCR validation of hub gene expression levels within skin lesions in the control and HS groups. (a). IL2RG. (b). FCGR2A. (c). CD19. (d). IL7R. (e). PTPRC. (f). SELL. (g). ITGAL. (h). CD27. (i). ITGAX. (j). CD2. Compared with the control group, *p<0.05, ** p <0.01, *** p<0.001.

Discussion

In the research, we utilized available bioinformatics analysis tools to explore the pathogenesis of HS at the molecular level and thus predict potentially effective therapeutic agents. Functional enrichment analysis of DEGs in HS highlighted the regulation of leukocyte migration and differentiation, serine-type endopeptidase and serine hydrolase activity, as well as interleukin (IL)-17 signaling pathway, cell adhesion molecules, and primary immunodeficiency as implicated in the disease.

HS is a chronic purulent inflammation of the apocrine sweat glands, commonly found in friction areas, such as the axilla, groin, and perianal area. Its clinical manifestations include painful inflammatory nodules, acne, scarring, and sinus tracts. In addition, HS can have serious consequences such as psychological and social disorders, activity limitations, and even skin cancer. The exact mechanism of HS pathogenesis is not well understood, and its treatment in the clinic is not universally effective. Current research has established that HS is a persistent inflammatory disorder, and many pro-inflammatory mediators, including interferon-γ, TNF-α, IL-17, and IL-1, play major roles in its pathogenesis.

The role of serine proteases is to cleave peptide bonds, thus degrading large proteins into small proteins. As described by İlgen et al, autoinflammatory syndromes such as HS are associated with mutations in a gene called proline-serine-threonine phosphatase-interacting protein.^29,30

The IL-17 signaling pathway has been well described in the pathogenesis of HS. The IL-17 family, which includes six members (IL-17A to IL-17F),^31,32 plays a key role in the pathogenesis of HS.^33,34 IL-17 levels are higher in the lesions, peri-lesions, and the unaffected skin of HS patients than in those without HS.^33,35 In addition, the serum IL-17 concentration is increased in HS patients.³⁶ IL-23 is involved in the activation of T helper 17 cells,³⁷ which induces the production of IL-17, thereby activating the pathway and enhancing downstream gene expression.

Currently, the treatment of inflammatory diseases has evolved from general immunosuppression to biological agents that target the IL-23/IL-17 signaling pathway.³⁸ Herein, DEGs associated with HS were also enriched for leukocyte migration and differentiation as well as cell adhesion molecules, which further confirmed the close association of HS with immune and inflammatory responses. Among the hub genes predicted for HS, both SELL and ITGAL have functions in regulating leukocyte migration and adhesion. In addition, immunodeficiency was correlated with HS. A greater susceptibility to HS is observed in HIV-infected patients, with a more aggressive course noted in such cases.³⁹ Furthermore, a case report of a boy with mosaic trisomy 13 and immunodeficiency who developed HS at 18 months of age showed that immunoglobulin therapy was effective.⁴⁰

In the PPI network, we filtered out 10 hub genes (PTPRC, CD19, IL7R, SELL, ITGAL, CD27, ITGAX, FCGR2A, CD2, IL2RG) and three protein modules. All 10 hub genes were presented in the core protein modules. The protein modules further demonstrate the significance of the hub gene, and also suggest that genes such as CCR5, CXCR4, THY1, CCR7 and KLRB1 are closely related to the hub genes. CCR5, as a member of the G protein-coupled receptor protein family, can facilitate communication between immune cells and the environment, thereby playing a crucial role in regulating inflammatory responses.⁴¹ For instance, CCR5 acts as an auxiliary receptor for HIV to enter cells, and it plays an important role in the immune process of HIV patients, and has been extensively studied.^42,43 The second protein module comprises the genes COL4A1, COL5A2, COL4A2, MMP14, LOXL2, and SERPINH1, which facilitate aberrant extracellular matrix remodeling and fibrotic processes.^44,45 The third protein module includes CYTH4, LY86, HCLS1, RAC2, S100A9, S100A8, S100A12, MMP9, and FGR, primarily linked to the immune response.^46–48 These are closely related to the pathogenesis of HS.

Of the 10 aforementioned hub genes, IL2RG and FCGR2A were verified as differentially expressed in HS skin tissues compared with those in normal skin tissues. Both FCGR2A and IL2RG act primarily in immune responses and are targets for a variety of autoimmune diseases.

FCGR2A, also known as immunoglobulin G Fc receptor II, is a member of the immunoglobulin Fc receptor gene family. The proteins encoded by FCGR2A (Fc γ Receptor IIa) are cell-surface receptors found on macrophages and neutrophils, being involved in the process of phagocytosis and immune complex clearance. Macrophages phagocytose and remove immune complexes from tissues and blood, and further facilitate the clearance of antigen–antibody complexes by stimulating the complement system and releasing inflammatory mediators. This series of processes is critical in HS, involving various immune cells and inflammatory factors. Furthermore, serum immunoglobulin G is a marker of disease severity in HS and can be used as an adjunct to clinical severity scoring. Therefore, it may be possible to treat HS with FCGR2A-related drugs (eg, cetuximab, etanercept, adalimumab, infliximab, and trastuzumab). There have been clinical studies of these drugs, such as a 2019 study reviewing TNF-α inhibitors for the treatment of HS, which included adalimumab, etanercept, and infliximab; the study confirmed the reliability of adalimumab for HS.⁴⁹

IL2RG (IL-2 receptor subunit gamma) encodes the IL-2 receptor gamma chain (IL-2Rγ), which is a common receptor subunit for a variety of immune cells and an important component in mediating immune responses. IL-2 can activate downstream signaling pathways, including the PI3K and JAK/STAT signaling pathways, by binding to IL-2R. Furthermore, IL2RG is a necessary component for the activation of downstream signaling pathways. mTOR is an important molecule in the PI3k signaling pathway, and De Vita et al found that mTOR gene expression is statistically correlated with the severity of HS.⁵⁰ The mTORC1 pathway plays an important role in itch signaling.⁵¹ Additionally, melatonin could treat proliferative scars by inhibiting fibroblast function and proliferative scar formation through the PI3K/Akt/mTOR pathway, suggesting that scar formation in HS might also be related to the PI3K/Akt/mTOR pathway. Whether IL2RG plays a role therein needs to be further verified.

Studies have shown that IL2RG can mediate cell proliferation, differentiation, and other signal transduction by participating in various pathways, such as JAK/STAT, essential for the development, differentiation, and function of T, B, and natural killer (NK) cells.⁵² X-linked severe combined immunodeficiency caused by IL2RG mutation is typically characterized by the absence of T cells and NK cells and the presence of some functional B cells, with IL-7 and IL-15, respectively, playing a key role in the maturation of T cells and NK cells. There is no direct evidence of the relationship between IL2RG and HS; however, it can be speculated that IL2RG may be involved in the pathogenesis of HS through modulation of the immune and inflammatory response. Further studies should address this.

Among the predicted hub genes, CD19 and ITGAL are protein-coding genes involved in the endogenous Hippo pathway in leukocytes. Recent studies have preliminarily confirmed that the Hippo pathway can promote extensive fibrosis in HS, which, in turn, interferes with drug penetration and affects the overall therapeutic response. The Hippo pathway is strongly associated with cancer development, and some studies have found that patients with HS are at a higher risk of developing cancer.^53,54 This indicates that it may be possible to interfere with HS pathogenesis by regulating the expression of CD19 and ITGAL. Further, early application of drugs regulating their expression may reduce the incidence of malignant tumors in HS patients. IL7R is a downstream effector of the Notch1 pathway.⁵⁵ Its expression is regulated by dynamic NOTCH1 binding to distal enhancers.⁵⁶ The Notch pathway is able to finely regulate the progression, migration, and apoptosis of keratinocytes in all skin layers, playing a key role in maintaining skin homeostasis^57,58 and sustaining chronic inflammation in HS.⁵⁹ Among the predicted hub genes, multiple genes (CD19, SELL, and ITGAX) were associated with the complement pathway. ITGAX is a complement receptor-related gene that encodes complement receptor 4.⁶⁰ The complement system, as a component of innate immunity, recognizes and triggers an inflammatory response to foreign antigens.⁶¹ In gene expression analyses, complement receptors were higher in HS lesion tissues than in non-lesion skin.⁶² Grand et al suggested that complement mechanism is implicated in HS pathogenesis and therapy.⁶¹ PTPRC is a protein tyrosine phosphatase receptor type C that activates protein tyrosine phosphatases required for T-cell activation through the antigen receptor. In some autoimmune diseases, after B-cell receptor activation, Bruton’s tyrosine kinase (BTK) and spleen tyrosine kinase (SYK) phosphorylation is reduced in memory B cells, while PTPRC expression is upregulated in B cells.⁶³ The expression of tyrosine kinases (BTK, SYK, JAK1, JAK3) is increased in the skin of patients with HS compared to that in healthy populations.⁶⁴ B-cell signaling, encompassing BTK and SYK signaling, is a pivotal event in HS.⁶⁵ It is worth exploring whether PTPRC is related to BTK and SYK, and exactly how it regulates HS.

A series of analyses of the hub genes were carried out. First, ROC curves were used to assess the diagnostic significance of genes. The mRNA–TF and mRNA–miRNA interactions of hub genes were predicted, highlighting the regulatory networks in hub gene transcription. Further, the relationship between genes and drugs was predicted. Our mRNA–TF network contains nine hub gene mRNAs and 69 TFs. Among them, GATA3 is a transcription factor for IL7R, ITGAL, and IL2RG, and its paralog GATA2 is a transcription factor for IL7R. GATA3 is a key factor in the process of Th2 cell differentiation following inflammatory responses. Th2 cells are important regulators of mammalian adaptive immunity. In a previous mouse study, the decrease in GATA3 expression led to alterations in skin innate lymphocytes and promoted hair follicle recirculation.⁶⁶ The transcription factor JUN is a target of ITGAX, IL7R, IL2RG, and CD27. It is highly expressed in various scars, and scar formation can be prevented by blocking JUN-mediated pathological processes.⁶⁷ This suggests that JUN may be involved in scarring and sinus tract formation in patients with HS. Moreover, JUN is a key regulator of the pathogenesis of many skin diseases and may be a target of drug therapy for psoriatic skin inflammation,⁶⁸ photoaging,⁶⁹ systemic sclerosis,⁷⁰ and squamous cell carcinoma⁷¹ of the skin.

MiRNAs play a role in cellular responses and are particularly important in inflammatory regulation and autoimmune diseases.⁷² Hessam et al revealed that the expression of the RNA-induced silencing complex components, which can influence miRNA formation and function, is dysregulated in HS, supporting the hypothesis that miRNAs may be involved in the inflammatory pathogenesis of HS.⁷³ A study in China demonstrated the reduced expression of miR-100-5p in patients having familial HS with nicastrin mutations, which in turn led to the upregulation of p-AKT and promotion of keratin-forming cell growth.⁷⁴ A study by Moltrasio et al showed that miR-24-1-5p, miR-146a-5p, miR26a-5p, miR 206, miR338-3p, and miR-338-5p can be used as biomarkers for the early diagnosis of HS.^75,76 Hessam et al found that the expression of some miRNAs such as miRNA-31-5p was upregulated in lesional HS skin samples compared to healthy controls and speculated that the miRNAs might manipulate the inflammatory pathway in HS. However, the association has not been verified. Further mechanistic studies are needed to illuminate the processes underlying the differential expression of the miRNAs in HS skin and to evaluate their therapeutic potential in treating HS.⁷⁷ In the study conducted by U. Radhakrishna,⁷ a series of miRNAs related to HS were identified and analyzed. The study found that miR-29, miR-200, miR-205, miR-548 and miR-132 were particularly important for skin function. In our research, when searching for miRNAs corresponding to hub genes, all of the above miRNAs were found in the miRDB database or the miRWalk database. MiR-548 mainly acts as a target for the genes IL7R and PTPRC, and is considered a potential biomarker for HS. Previous studies have shown that miR-1238-3p, as a key factor, can exacerbate the malignant progression of skin squamous cell carcinoma.⁷⁸ Moreover, numerous studies have confirmed that patients with HS have a higher risk of developing skin squamous cell carcinoma. Another miRNA, hsa-miR-575, which serves as a target for SELL, has been shown in previous studies to be related to the invasion, proliferation and migration of glioma cells.⁷⁹ Among the newly discovered miRNAs, the miR-3689 family (hsa-miR-3689a-3p, hsa-miR-3689b-3p, hsa-miR-3689c) targets the gene FCGR2A. Among them, miR-3689c has been found to be related to skin tumors.⁸⁰ Furthermore, a study have discovered that hsa-miR-6133 is a common target for the genes IL7R, ITGAL, PTPRC, and SELL. It was found that through the mediation of TGFBR2 and JNK, miR-6133-5p has a strong anti-fibrotic effect.⁸¹ More studies are required to confirm the function of the miRNAs. The 115 miRNAs that we identified provide new insights into the current understanding of HS pathogenesis.

In the predicted gene–drug interaction network, prednisone,⁸² rituximab,⁸³ and cyclosporine⁸⁴ have been used for the treatment of HS.

Efalizumab was used in a trial to treat HS; however, efficacy was not achieved.⁸⁵ Patients with HS are more likely to develop squamous cell carcinoma, and fluorouracil is an effective topical treatment of squamous cell carcinoma.^86,87 However, there have been no cases of HS treated with fluorouracil. Ruxolitinib has not been used directly for the treatment of HS, but an experiment⁸⁸ has shown that ruxolitinib, as a systemic JAK inhibitor, attenuates the expression of inflammatory cytokines and chemokines in the keratinocytes of HS lesions, thus having potential as a topical therapeutic agent.

The following drugs are not used in HS treatment based on the available literature: hydrocortisone butyrate, clobetasol propionate, estradiol valerate, alendronate sodium, hydrogen peroxide, daclizumab, denileukin diftitox (DD), basiliximab, aldesleukin, sirolimus (use of this drug may induce HS),⁸⁹ epoetin alfa, busulfan, cyclophosphamide anhydrous, anti-thymocyte globulin, interferon beta-1a, paclitaxel, doxorubicin hydrochloride, margetuximab-cmkb, tisagenlecleucel, inebilizumab-cdon, brexucabtagene autoleucel, lisocabtagene maraleucel, axicabtagene ciloleucel, tafasitamab, loncastuximab tesirine, benztropine, etoposide, lifitegrast, thrombin, alefacept, blinatumomab, monomethyl fumarate, and octreotide acetatezinc chloride. While hydrocortisone butyrate, estradiol valerate, clobetasol propionate, anhydrous tacrolimus, monomethyl fumarate, and zinc chloride were reported as HS treatments, hydrocortisone,⁹⁰ estradiol,⁹¹ clobetasol,⁹² tacrolimus,⁹³ fumarate,⁹⁴ and zinc⁹⁵ have been used. Of the unused drugs, four that target IL2RG are noteworthy. Daclizumab is a humanized monoclonal antibody against CD25, which is the alpha subunit of the high-affinity IL-2 receptor, which can selectively regulate IL-2 signaling.⁹⁶ Daclizumab treatment inhibits the function of pro-inflammatory activated T-lymphocytes and promotes the expansion of immunoregulatory CD56bright natural killer cells.^97–100 Currently, it is primarily used to treat multiple sclerosis. DD is a recombinant fusion protein of diphtheria toxin fragments and IL-2,¹⁰¹ which targets IL-2 receptor on malignant T lymphocytes,¹⁰² binding to high- and neutral-affinity IL-2Rs through the IL-2 structural domains in human T-cell lymphoma-derived cell lines. Basiliximab is an IL-2 receptor antagonist, which can directionally antagonize the receptor alpha chain of IL-2 (CD25 antigen), blocking the binding of IL-2 to the IL-2 receptor and thus blocking the transmission of proliferative messages to T cells.^103,104 Finally, aldesleukin is a recombinant human IL-2 that promotes effector cell function at high doses. Soluble IL-2 receptor (sIL-2R) levels are elevated in patients with HS.¹⁰⁵ Daclizumab, DD, basilixima, and aldesleuki are IL2RG-targeting drugs that modulate sIL-2R levels; however, no studies have reported their use for the treatment of HS. Therefore, whether daclizumab, denileukin, diftitox, basilixima, and aldesleukin are effective in the treatment of HS deserves further investigation.

We used the “CIBERSORT” package to evaluate immune cell infiltration in patients with HS. M1 macrophages are closely associated with regulatory and activated memory CD4 T cells. Neutrophils are closely associated with resting mast cells, regulatory T cells, and activated mast cells. Resting mast cells are closely associated with follicular helper T cells and plasma cells. However, the mechanisms underlying these associations require further in-depth investigation. Finally, we examined the expression levels of 10 hub genes in HS using qRT-PCR. FCGR2A and IL2RG expression was consistent with the results of bioinformatics analysis.

This study has some limitations. First, the sample size was small, which might have affected the credibility of the data. Second, sex differences are obvious in HS but were not considered in the current study. The proportion of men and women in the dataset sample was not balanced, as all samples in GSE137141 were from female individuals. Third, some studies have reported racial differences in HS. The datasets analyzed herein are from Sweden and Germany, while the experimental qPCR data are from China. Therefore, the racial factor has not been controlled for. In the future, we hope that the hub genes identified will be further confirmed experimentally. Further, the predicted drugs should be experimentally validated to advance the development of HS treatments.

Conclusion

Our bioinformatics analysis led to the identification of 10 hub genes, with roles in signaling pathways. We also identified potential drugs and analyzed their association with immune-infiltrating cells. Among hub genes, FCGR2A is responsible for the phagocytosis and clearance of immune complexes. The dysregulation of these processes may be an important event in the pathogenesis of HS. The PI3K-Akt-mTOR signaling pathway, regulated by the IL2RG gene, showed a close relationship with HS severity and the mechanism of scarring, thus representing a potential therapeutic target. The miR-3689 targeting the gene FCGR2A might also be potential biomarkers.