Introduction

Chronic rhinosinusitis (CRS) is an inflammatory process affecting the mucosa of the nasal cavity and one or more sinuses. Its clinical manifestations include nasal congestion, nasal leakage, loss or reduction of the sense of smell, and facial or head pain, lasting for more than 12 weeks. These symptoms severely affect the quality of life of patients and impose substantial social and economic burdens.1 The pathogenesis of CRS is complex, involving genetic factors, the sinus microbiota, infections, and environmental influences, and has not yet been fully elucidated.2 Clinically, CRS is generally classified into two subtypes: CRS without nasal polyps (CRSsNP) and CRS with nasal polyps (CRSwNP), with a prevalence ratio of approximately 4:1. This disease is also characterized by its pathophysiological mechanisms. These mechanisms are further classified into type 2 inflammation and non-type 2 inflammation based on the characteristics of the main effector T cells and innate lymphoid cells (ILCs) involved in the inflammatory response.3,4 These distinct inflammatory endotypes reflect the involvement of different immune cell types (eg, T cells, B cells) and cytokines in the immune system.5,6 In diseases such as CRSwNP, variations in inflammatory endotypes may influence disease progression, symptom severity, and treatment responsiveness. The primary pharmacological treatments for CRSwNP include topical and oral glucocorticosteroids, antibiotics, antihistamines, antileukotrienes, and saline nasal rinses. However, conventional therapies often overlook the phenotypic differences and intrinsic pathophysiological mechanisms of CRSwNP, resulting in limited symptomatic relief.7 Currently, the disease is incurable and prone to frequent recurrence. Multiple clinical and biomarker parameters have been reported to be associated with type 2 inflammation and response to biologic therapies, including clinical scores (eg, SNOT-22, LMS) and nasal secretion markers (eg, eosinophil cationic protein, interleukin-5, etc). However, approximately 40–60% of patients show poor response to these biologics, and the predictive validity of these biomarkers still requires further verification.8 Additionally, the complexity of gene expression patterns in CRSwNP represents a crucial factor in fully understanding its underlying pathogenesis. Therefore, the development of clinically applicable biomarkers would enable personalized treatment approaches for CRSwNP patients and could potentially significantly improve therapeutic outcomes.

Whole transcriptome RNA-sequencing (RNA-seq) is a robust technique used to analyze the complete genome of RNA extracted from tissues.9 By comparing the transcriptome of nasal polyps with that of control nasal mucosa, gene mutations associated with the disease can be identified. Moreover, this approach can reveal molecular differences between different types of CRSwNP, which may serve as potential therapeutic targets, thus contributing to the development of personalized treatment options.10 Therefore, RNA-seq serves as a powerful tool for identifying gene signatures and candidate pathways associated with the pathogenesis of CRSwNP, informing novel potential therapeutic targets.

Currently, the recommended treatment for CRSwNP includes pharmacological therapy with systemic (oral) and topical (nasal) glucocorticoids, antimicrobial agents, antihistamines, and leukotriene receptor antagonists.11 However, some patients with severe disease do not respond adequately to pharmacological therapies and may require repeated systemic corticosteroid courses and/or sinus surgery.12 Even after successful surgery, recurrence rate can still be as high as 40% within 18 months, with approximately 80% of patients experiencing recurrence after 12 years. The occurrence of some inflammatory diseases (such as autoimmune diseases, allergic diseases, etc) is related to epigenetic changes.13,14 Epigenetic regulation refers to the modulation of cellular behavior through alterations in gene expression patterns without changing the DNA sequence, typically involved in processes such as cell differentiation, development, and stress responses.15 In these diseases, the epigenetic regulation of genes is abnormal, leading to the overproduction of inflammatory factors and the persistence of the inflammatory response. Inflammatory factors (such as tumor necrosis factor-α, interleukin-1, interleukin-6, etc) regulate histone modification, DNA methylation, and non-coding RNA expression, thereby affecting gene expression and cell function. The inflammatory microenvironment affects the epigenetic characteristics of stem cells, causing them to differentiate into inflammatory cells or fibrocytes.16 In addition, inflammatory factors and inflammatory cells can also regulate the proliferation, invasion and metastasis of tumor cells by influencing the epigenetic characteristics of tumor cells.17 In conclusion, there is a strong link between inflammation and epigenetics. Signaling molecules produced during inflammation can affect epigenetic regulatory mechanisms, and epigenetic changes can also affect the occurrence and development of inflammatory responses. This interaction is of great significance for the occurrence and development of inflammation-related diseases, and also provides new ideas for the treatment of inflammation-related diseases. In recent years, biologic agents targeting type 2 inflammation have emerged as a potential avenue for individualized treatment in difficult-to-control, severe CRSwNP cases.18 Furthermore, although plasma IgE levels are also increased in patients with CRSwNP, there is no association of these IgE levels with the presence of atopy or environmental allergies. It is not yet clear what clinical benefit may be seen if CRSwNP patients are treated with an anti-IgE agent.19 Previous studies have used RNA sequencing to identify transcriptome signature, gene features and candidate pathways associated with CRSwNP by comparing CRSwNP patients and controls will help identify disease targets. The known biological processes implicated in CRSwNP are multifaceted, including airway inflammation, tight junction impairment, pathogen infection and a defective host defence.20 Although eosinophilia reportedly predisposes NP formation.21 Interestingly, in the study by yang et al, suggested that eosinophilia was not a major confounder. The unique DEGs for the aforementioned comparisons revealed that pathways associated with the molecular mechanisms of cancer and neutrophil N-formylmethionyl-leucyl-phenylalanine signaling drive polyp formation. Suppressing exaggerated cell growth might, therefore, be an effective management approach for CRSwNP in the future, given that the rapid recurrence of CRSwNP mimics benign tumor growth.22 These biologics primarily act through several mechanisms, including anti-IgE, anti-interleukin (IL)-4Ra, anti-IL-5, and anti-IL-5Ra, blocking specific inflammatory factors or receptors to modulate the inflammatory response.23 However, the complete remission rate with immunotherapy alone for CRSwNP remains low. Notably, increasing evidence suggests that epigenetic dysregulation is a key factor in the poor response to chemoimmunotherapy in patients with CRSwNP.15 Epigenetics regulates gene expression through DNA methylation, histone modification, and non-coding RNA editing.24 Among them, DNA methylation is integral to gene expression and plays an important role in inflammation. Therefore, normalizing epigenetic regulation is essential to improving the success of immunotherapy for CRSwNP.25

Enhancer of zeste homolog 2 (EZH2) is the catalytic subunit of the polycomb repressive complex 2 (PRC2). It inhibits the transcription of its target genes by catalyzing the trimethylation of lysine 27 on histone H3 (H3K27me3).26 EZH2 plays a crucial role in regulating the division of various immune cells in allergic rhinitis. For example, the exosomal long non-coding RNA GAS5 inhibits Th1 differentiation and promotes Th2 differentiation by downregulating the expression of both EZH2 and T-bet.27 EZH2 deficiency in CD4+ T cells enhances the Th2 immune response in ovalbumin (OVA)-induced allergic airway models.28 Moreover, 3-deazaneplanocin A, an indirect inhibitor of EZH2, inhibits dendritic cell activation in vitro.29 Therefore, EZH2, as a central regulator of immune cell differentiation in allergic rhinitis, plays an important role in chronic inflammation.

In the present study, we aimed to identify biomarkers and pathways associated with the pathogenesis of CRSwNP. To this end, we performed RNA-seq in polyp tissues and normal mucosal tissues to identify transcriptomic features associated with CRSwNP. We also conducted a comprehensive analytical pipeline to determine differentially expressed genes (DEGs) and gene sets enriched in functional pathways. Moreover, we investigated the involvement of key biomarkers in the inflammatory response of CRSwNP through epigenetic regulation in nasal mucosal epithelial cells. This study offers new research directions for the application of immunotherapy in CRSwNP.

Materials and Methods Source of Data

This work was approved by the Institutional Review Board (2022-department-564) of the Ethics Committee of Beijing Chaoyang Hospital. All human studies were conducted according to institutional ethical norms and the guidelines outlined in the Declaration of Helsinki. All participants provided written informed consent before participation. A total of 86 individuals were included between July 2021 and July 2023, comprising 43 patients with CRSwNP who underwent nasal polyp surgery and another 43 patients who had a deviated septum surgically corrected. CRSwNP was diagnosed following the European Position Paper on Rhinosinusitis and Nasal Polyps 2020 guidelines.30 Patients who had taken oral glucocorticoids, antimicrobials, anti-leukotrienes, or antihistamines within four weeks prior to sample collection were excluded from the study. Moreover, patients with concurrent aspirin sensitivity, allergic rhinitis, or asthma, along with a deviated nasal septum, were also excluded. The control group consisted of 43 patients with septal deviation, while the CRSwNP group comprised 43 patients with CRSwNP. Eight pairs of samples (8 chronic rhinosinusitis with nasal polyps and 8 control samples) were collected for transcriptome sequencing to generate the CRSwNP transcriptome sequencing (TS-CRSwNP) dataset. Immunohistochemistry (IHC) analysis was performed on tissue samples obtained from 20 patients with nasal polyps (NPs) in the CRSwNP group and 20 nasal mucosal tissue samples from the control group. Quantitative real-time reverse transcription polymerase chain reaction (qRT-PCR) analysis was performed using tissue samples from 15 patients with NP in the CRSwNP group and 15 nasal mucosal tissue samples from the control group for validation. The GSE136825 dataset, sourced from the Gene Expression Omnibus (GEO) database, was used for gene validation. This dataset (GPL20301) includes RNA-seq data from turbinate samples, including 42 CRSwNP samples, 33 nonpolyp inferior turbinate (CRSwNP-IT) samples and 28 control samples.22 The flowchart of this study is shown in Figure 1.

Figure 1 The flowchart of this study.

Transcriptomic Sequencing

Total RNA from nasal polyp samples was isolated using TRIzol (Invitrogen, Carlsbad, CA, USA) reagent, according to the manufacturer’s instructions. The RNA concentration and purity were quantified using the NanoDrop ND-1000 spectrophotometer (NanoDrop, Wilmington, DE, USA). RNA integrity was evaluated through agarose gel electrophoresis. Oligo (dT) magnetic beads (Dynabeads Oligo [dT]; Thermo Fisher Scientific, Waltham, MA, USA) were used for the specific capture of mRNA with PolyA tails. Subsequently, fragmented RNA was reverse-transcribed into cDNA using SuperScript™ II Reverse Transcriptase (Invitrogen). Two-strand synthesis was performed using Escherichia coli DNA polymerase I (New England Biolabs [NEB], Ipswich, MA, USA) with RNase H (NEB, USA). The RNA-DNA hybrids were converted into double-stranded DNA. The second strand was digested using the UDG enzyme (NEB). PCR amplification was performed to construct a library with a fragment size of 300 ± 50bp. Finally, double-end sequencing was conducted using Illumina novaseq 6000 (LC Bio Technology CO., Ltd. Hangzhou, China) according to standard protocols.

Identification of DEGs, Differentially Expressed microRNAs (DE-miRNAs), Differentially Expressed-Long Non-Coding RNAs (DE-lncRNAs), and Differentially Expressed-Circular RNAs (DE-circRNAs)

DEGs, DE-miRNAs, DE-lncRNAs, and DE-circRNAs between the CRSwNP and control groups were identified using the DESeq2 package (v 1.38.0)31 in the TS-CRSwNP dataset, applying the criteria of adjusted p value < 0.05 and |log2FC| > 0.5.32,33 The Benjamini - Hochberg (BH) method was used for calibration. Analysis of variance was visualized using a volcano map and a heatmap generated using the ggplot2 package (v 3.4.1)34 and pheatmap package (v 1.0.12), respectively. The results for DE-circRNAs were displayed only via a volcano plot.

Functional Annotation of DEGs, DE-miRNAs, and DE-lncRNAs

The target mRNAs of DE-miRNA were predicted using miRTarBase and TarBase, and the DE-miRNA target genes were identified by intersecting the target mRNA and DE-mRNA after removing duplicates. The pairs of DE-lncRNAs and DE-mRNAs with a correlation coefficient (r) > 0.5 and p < 0.05 were screened using the Spearman method, and DE-lncRNAs target genes were obtained after removing duplicates. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of DEGs, DE-miRNAs target genes, and DE-lncRNAs target genes were performed using the clusterProfiler package (v 4.7.1)35 and org.Hs.eg.db (v 3.16.0) (p value < 0.05).

Construction of a Competing Endogenous RNAs (ceRNAs) Regulatory Network

The StarBase and miRNet databases were utilized to predict lncRNAs targeting miRNAs. The predicted lncRNAs were intersected with DE-lncRNAs to obtain co-lncRNAs. Combined with the prediction results from the previous step, the lncRNA-miRNA-mRNA network was constructed using Cytoscape software.36

The StarBase and circBank databases were utilized to predict circRNAs targeting miRNAs. The relationship pairs between circRNAs and DEGs with r > 0.5 and p < 0.05 were identified using the Spearman method. Finally, a circRNA-miRNA-mRNA network was constructed.

Integration of ceRNA Networks and Functional Enrichment

The lncRNA-miRNA-mRNA and circRNA-miRNA-mRNA networks were integrated, and relationships linked by the same miRNA were used to construct ceRNA networks. Genes within the ceRNA network were considered candidate genes. GO and the KEGG enrichment analysis of these candidate genes were performed using the clusterProfiler package (v 4.7.1)35 and org.Hs.eg.db (v 3.16.0), with a significance threshold of p < 0.05.

Machine Learning Screening for Biomarkers

First, to explore the interactions between candidate genes, a protein–protein interaction (PPI) network was established using the STRING website (https://string-db.org) (confidence = 0.4). Subsequently, the sub-networks were analyzed using the plugin CytoHubba, which is capable of assigning values to each node in the network by applying multiple topological algorithms to identify key genes and sub-networks. Based on this, the five gene networks with the highest degree values were selected as the core networks, and the genes in the networks were identified as hub genes. Feature genes were first screened based on the hub genes using the XGBoost model, which was constructed via the xgboost package (v 1.7.3.1) (https://CRAN.R-project.org/package=xgboost). XGBoost is an ensemble learning algorithm based on gradient boosting decision trees, which can optimize model performance by adjusting various hyperparameters. In constructing the XGBoost model, max_depth was set to 5, eta to 0.3, objective to “binary:logistic”, and nround to 25. Additionally, feature genes were screened using the Support Vector Machines-Recursive Feature Elimination (SVM-RFE) model, constructed using the caret package (v 6.0–93) (https://CRAN.R-project.org/package=caret). SVM-RFE identified the most critical feature subset for classification by progressively eliminating features. A five-fold cross-validation strategy was employed to more robustly evaluate model performance. Biomarkers were selected by overlapping SVM-RFE-feature-genes and XGBoost-feature-genes, and their diagnostic performance was assessed using receiver operating characteristic (ROC) curves. Finally, the correlation between the biomarkers was calculated using Spearman correlation analysis using the corrplot package (v 0.92) (https://github.com/taiyun/corrplot).

Construction of the Artificial Neural Network (ANN) Model

The ANN model was developed using a neural network software package based on biomarkers that minimizes the error between the model’s predicted and actual values by adjusting the weights and biases.37 The importance of each biomarker for the predictive outcome of the model was computed using the ANN network. Subsequently, the performance of the model was evaluated using a confusion matrix and ROC curve analysis.

Single-Gene Gene Set Enrichment Analysis (GSEA) Analysis

Single-gene GSEA was performed to identify the enriched regulatory pathways and biological functions of the biomarkers using clusterProfiler (v 4.7.1),35 with adjusted p < 0.05 and normalized enrichment score |NES| ≥ 1. The results for the top 5 pathways with the largest absolute values of NES, both greater than 0 and less than 0, were visualized.

Relationship Between Biomarkers and Immune Factors

The CIBERSORT algorithm, which is specifically designed to estimate the abundance of different cell types in mixed cell populations, was used in this study to calculate the proportions of 15 immune cell subtypes in each sample (excluding immune cells with zero abundance in more than 75% of the samples).38 To investigate the association between biomarkers and different immune features, the immune features (chemokines and major histocompatibility complex [MHC]) associated with CRSwNP were retrieved from the TISIDB database. Differential profiles of these immune features between the CRSwNP and control groups were then assessed using the Wilcoxon test. Finally, the correlation between biomarkers with differential immune cells, chemokines, and MHC was analyzed using the Spearman method.

Construction of a mRNA-Transcription Factor (TF) Network

To understand the transcriptional regulatory mechanisms of biomarkers, TFs targeting the biomarkers were predicted using the ChEA3 and NetworkAnalyst databases, and the mRNA-TF regulatory network was obtained using Cytoscape.36

Construction of Biomarker-Drug Interaction Network

To explore new therapeutic targets for CRSwNP treatment, drug predictions for biomarkers were performed. Initially, the DGIDB database (https://dgidb.org/) was used to predict potential drugs targeting the biomarkers. Subsequently, a biomarker-drug network was constructed based on the predicted results.

Analysis of Molecular Regulatory Axes

The sub-network containing biomarkers was extracted from the ceRNA network and displayed using the Sankey diagram.

Relationship Between Biomarkers and Diseases

To analyze the role of biomarkers in other ear-nose-throat (ENT) diseases, the Comparative Toxicogenomics (CTD) database (https://ctdbase.org/) was used to assess their association with various ENT conditions. The results are presented using bar charts.

Expression Analysis

To validate the expression of the biomarkers in CRSwNP and control tissues, expression analysis was performed using both the TS-CRSwNP and GSE136825 datasets.

Cell Lines and Culture

Human normal nasal mucosal epithelial cell line (HNEpCs) was obtained from WHELAB (Shanghai, China). Cells were cultured in epithelial cell media (M1005A; WHELAB) at 37°C in a 5% CO2 incubator.

siRNA and Plasmid Transfection

The siRNA used in this study was developed by SyngenTech (Beijing, China). The plasmid for overexpression of EZH2 was synthesized by Thermo Fisher Scientific. Lipofectamine 3000 (Life Technologies, Thermo Fisher Scientific) was used for siRNA and plasmid transfections into HNEpCs, following the manufacturer’s instructions. The siRNA sequences are as follows:siEZH2-1:5′-GACUCUGAAUGCAGUUGCUTT-3′;

siEZH2-2:5′-CGGCUCCUCUAACCAUGUUTT-3′; siEZH2-3:5′-CAGAAGAACUAAAGGAAAATT-3′; and control:5′-UUCUCCGAACGUGUCACGUTT-3′.

The sequence information on EZH2 overexpression plasmids is available in the Supplementary Information 1. The empty vector pCDNA3_Zsgreen-T2A-puro was used as a control.

RNA Extraction, Reverse Transcription, and Quantitative Real-Time Reverse Transcription Polymerase Chain Reaction (qRT-PCR)

qRT-PCR is a technique used to measure the quantity of RNA in biological samples by converting RNA into DNA (cDNA) through reverse transcription and then amplifying and quantifying the cDNA in real-time using fluorescent markers. The First-Strand cDNA Synthesis Kit (Transgene, Beijing, China) was used to reverse-transcribe 1 g of the extracted total RNA. The resulting cDNAs were then assessed using the QuantStudioTM 5 System for qRT-PCR with SYBR mix (Transgene, Beijing, China).

The following primers were used:

EZH2: F, 5′-AATCAGAGTACATGCGACTGAGA-3′;

R,5′-GCTGTATCCTTCGCTGTTTCC-3′;

Insulin-like growth factor 2 mRNA-binding protein 1 (IGF2BP1): F, 5′-TCGTTGCAAGACCTTACCCTT-3′;

R, 5′-GGCAGCCACATCATTCTCATAG-3′;

High-mobility gene group A2 (HMGA2): F, 5′-CTCAAAAGAAAGCAGAAGCCACTG-3′;

R, 5′-TGAGCAGGCTTCTTCCGAACAACT-3′

GAPDH: F, 5′-GGAAGCTTGTCATCAATGGAAATC-3′;

R, 5′-TGATGACCCTTTTGGCTCCC-3′.

Immunohistochemistry (IHC) Analysis

Immunohistochemistry (IHC) analysis is a staining technique that utilizes the specificity of antibodies to target and visualize antigens in tissue sections, aiding in the understanding of cell types and their functions within a tissue. Tissue samples were heated at 65 °C for 2 h, then rehydrated and dewaxed using xylene and gradient ethanol. Samples were then boiled in an EDTA solution (pH = 8.0) for 20 min for antigen retrieval. To block non-specific binding, 2% goat serum was used. Tissue samples were then incubated with primary antibodies overnight at 4°C (anti-IGF2BP1, dilution 1:200; anti-EZH2, dilution 1:100; anti-HMGA2, dilution 1:100; all obtained from Abcam, Cambridge, UK).

Immunofluorescence (IF) Staining

IF staining is a technique used to visualize specific antigens in cells or tissues by labeling them with fluorescent antibodies, allowing for the detection and localization of these antigens under a microscope. The expression of H3k27me3 and HMGA2 was detected using immunofluorescence staining. Following a 20-min exposure to 4% paraformaldehyde, cells were blocked for 2 h at room temperature using 1% BSA (Gibco; Thermo Fisher Scientific) and 0.1% Triton-X. After overnight incubation at 4°C with the primary antibody against H3k27me3 (anti-H3k27me3, 1:1000 dilution; Cell Signaling Technology, Danvers, MA, USA) and HMGA2 (anti-HMGA2, 1:100 dilution; Abcam), the cells were incubated for 2 h at room temperature with the goat anti-rabbit IgG-H&L secondary antibody (cat. no. ab150077; 1:500; Abcam). Nuclei were stained with DAPI (Beyotime Institute of Biotechnology, Suzhou, China) for 5 min at room temperature, and stained cells were examined under a fluorescence microscope.

Statistical Analysis

All bioinformatics analyses were performed using R software (v 4.2.2) (R Foundation, Vienna, Austria). The Benjamini - Hochberg method was used for calibration. Spearman correlation analysis was conducted to assess relationships between variables. Data from different groups were compared using the Wilcoxon test, with p < 0.05 considered statistically significant.

Results Acquisition of DEGs, DE-miRNAs, and DE-lncRNAs and Construction of lncRNA-miRNA-mRNA and CircRNA-miRNA-mRNA Networks

In total, 3298 DEGs (2683 upregulated and 615 downregulated), 201 DE-miRNAs (105 upregulated and 96 downregulated), 204 DE-lncRNAs (170 upregulated and 34 downregulated), and 16 DE-circRNAs (16 downregulated) between the CRSwNP and control groups were obtained (Figure 2A–D, Supplementary Tables 1–4). As presented in the heatmap, the CRSwNP group exhibited decreased expression levels of LTF, BPIFB2, ZG16B, AZGP1, DMBT1, STATH, MUC5B, PPP1R1B, LPO, and LYZ, but elevated expression levels of CTAG2, CST1, CCL18, IGFBP3, HCRTR2, GOLGA6L1, UNC93A, SLITRK3, FAM71B, and SHH compared to the control group (Figure 2E–G).

Figure 2 Identification of differentially expressed genes (DEGs), differentially expressed microRNAs (DE-miRNAs), differentially expressed-long non-coding RNAs (DE-lncRNAs), and differentially expressed-circular RNAs (DE-circRNAs). (A–D) Volcano plots of DEGs (A), DE-miRNAs (B), DE-lncRNAs (C), and DE-circRNAs (D). (E–G) Heatmaps of the DEGs (E), DE-miRNAs (F), and DE-lncRNAs (G). (H) Construction of competing endogenous RNAs (ceRNA) networks. Orange ovals, purple diamonds, dark green triangles, and light green triangles indicate mRNAs, miRNAs, circRNAs, and lncRNAs, respectively.

Through prediction, we identified 5525 pairs of DE-miRNA-DE-lncRNA regulatory relationships. The lncRNA-miRNA-mRNA network contained 579 DEGs (including CSNK1A1L, CYP2U1, ATP13A1, ATP13A1, and ASPM), 57 miRNAs (including hsa-mir-519a-3p, hsa-mir-106b-5p, hsa-mir-20b-5p, and hsa-mir-92b-3p), and 29 lncRNAs (including SYNJ2-IT1, CBR3-AS1, INTS6-AS1, and FER1L6-AS2) (Supplementary Figure 1A). The hsa-let-7a-5p-AGO4 and other similar pairs were specific miRNA-mRNA pairs, while the hsa-let-7c-5p-CDKN2B-AS1 and others were miRNA-lncRNA pairs (Supplementary Table 5). Similarly, the circRNA-miRNA-mRNA network contained 16 DEGs (LRIG3 SFSWAP, STX3, ONECUT2, DNAJC6, ATP6V1B2, AGO4, CFL2, CRY2, LTA4H, BEND4, XBP1 SRCAP, ID1, PDGFB, and MTUS1), 3 miRNAs (hsa-let-7a-5p, hsa-let-7c-5p, and hsa-let-7g-5p), and 7 circRNAs (hsa_circ_0001009, hsa_circ_0000489, hsa_circ_0000799, hsa_circ_0000095, hsa_circ_0001725, hsa_circ_0000160, and hsa_circ_0000384) (Supplementary Figure 1B, Supplementary Table 6). After integrating the lncRNA-miRNA-mRNA and circRNA-miRNA-mRNA networks, we obtained the ceRNA network. This network comprised 46 DEGs (including AGO4, USP14, DNA2, PDGFB, and GP5), 3 miRNAs, 2 lncRNAs (LINC00885 and CDKN2B-AS1), and 7 circRNAs (hsa_circ_0001009, hsa_circ_0000489, hsa_circ_0000799, hsa_circ_0000095, hsa_circ_0001725, hsa_circ_0000160, and hsa_circ_0000384) (Figure 2H, Supplementary Table 7).

Enrichment Analysis of DEGs, DE-miRNAs, and DE-lncRNAs

Enrichment analysis indicated that the DEGs were associated with 613 GO terms and 42 KEGG pathways. The target genes of DE-miRNAs were associated with 566 GO terms and 53 KEGG pathways, while those of DE-lncRNAs target were associated with 601 GO terms and 42 KEGG pathways. The DEGs were mainly enriched in GO terms such as metal ion transmembrane transporter activity and actin-based cell projection (Figure 3A, Supplementary Table 8), and KEGG pathways such as protein digestion and absorption and synaptic vesicle cycle (Figure 3B, Supplementary Table 9). The target genes for DE-miRNAs were mainly enriched in GO terms such as mitotic cell cycle phase transition and nuclear division (Figure 3C, Supplementary Table 10), and in KEGG pathways such as the MAPK signaling pathway and cell cycle (Figure 3D, Supplementary Table 11). The target genes for DE-lncRNAs target genes were mainly enriched in GO terms such as metal ion transmembrane transporter activity and synaptic membrane (Figure 3E, Supplementary Table 12), and KEGG pathways such as glutamatergic synapse and protein digestion and absorption (Figure 3F, Supplementary Table 13).

Figure 3 Results of the enrichment analysis of the differentially expressed genes (DEGs), differentially expressed microRNAs (DE-miRNAs), and differentially expressed-long non-coding RNAs (DE-lncRNAs). The relationship between DEGs (A), DE-miRNAs (C), DE-lncRNAs (E), and major terms annotated by Gene Ontology (GO) is indicated by the tree diagram and circos graph. The column plot and circle diagram of the top 10 Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment results for the DEGs (B), DE-miRNAs (D), and DE-lncRNAs (F).

Functional Enrichment Analysis of Candidate Genes

We used the 46 mRNAs in the ceRNA network as candidate genes for subsequent analyses. To further investigate the functions performed by these candidate genes, we performed functional enrichment analyses. The results indicated that the candidate genes were enriched in 419 GO terms and 5 KEGG pathways. The enriched GO terms included CRD-mediated mRNA stabilization, regulation of muscle cell differentiation, and meiotic spindle organization (Supplementary Figure 2, Supplementary Table 14). The GO biological process (BP) enrichment analysis identified several functional categories, primarily including the regulation of the miRNA metabolic process, RNA catabolic process, meiotic cell cycle, and myoblast differentiation, among others (Supplementary Figure 3A). KEGG enrichment highlighted key pathways, such as collecting duct acid secretion, microRNAs in cancer, and the p53 signaling pathway (Supplementary Figure 3B, Supplementary Table 15). The KEGG pathways were predominantly distributed across five major categories: human diseases, organismal systems, cellular processes, environmental information processing, and metabolism (Supplementary Figure 3C).

Screening and Performance Evaluation of Biomarkers

We constructed a PPI network based on 46 candidate genes, resulting in a network with 17 nodes and 16 edges, after excluding isolated proteins (Figure 4A). Subsequently, we identified five hub genes (EZH2, ASPM, IGF2BP3, HMGA2, and IGF2BP11) through further analysis using the plugin CytoHubba (Figure 4B). Moreover, using machine learning screening, we selected four SVM-RFE-feature genes (EZH2, IGF2BP3, HMGA2, and IGF2BP11) and four XGBoost-feature genes (IGF2BP1, EZH2, HMGA2, and ASPM) (Figure 4C and D). Subsequently, we overlapped these eight feature genes and identified three biomarkers (IGF2BP1, EZH2, and HMGA2) (Figure 4E).

Figure 4 Screening of the three biomarkers. (A) Construction of the protein–protein interaction (PPI) network of 17 proteins. Each green circle represents a protein, while the gray lines represent the interaction relationships between proteins. (B) Construction of the core network with genes of TOP 5 degree value. The intensity of the red color indicates a higher degree value of the protein node. (C) Construction of Support vector machines-recursive feature elimination (SVM-RFE) models for the five candidate genes. (D) Column plot of genetic coefficient for the 4 XGBoost-feature genes. (E) Venn diagram of the four SVM-RFE-feature genes and four XGBoost-feature genes.

The results of the ROC curves revealed that the area under the curve (AUC) values of the three biomarkers all exceeded 0.7, indicating that these genes had good diagnostic values (Figure 5A–C). Correlation analysis revealed a significantly positive correlation between EZH2 and HMGA2 (Figure 5D). Moreover, the ANN analysis revealed that the optimal number of hidden neuron layers was three (Figure 5E). Notably, HMGA2 was the most influential variable affecting the predictive outcomes of the model (Figure 5F). The confusion matrix indicated that the model achieved 100% accuracy, while the ROC curve further indicated good stability (AUC = 0.719) (Figure 5G–H).

Figure 5 Performance evaluation of the three biomarkers. (A–C) Receiver operating characteristic (ROC) curves of the enhancer of zeste homolog 2 (EZH2) (A), the high-mobility gene group A2 (HMGA2) (B), and Insulin-like growth factor 2 mRNA-binding protein 1 (IGF2BP1) (C). The horizontal axis represents specificity, and the vertical axis represents sensitivity. The area under the curve (AUC) represents the accuracy (prediction performance). (D) Correlation analysis of the three biomarkers. (E) Construction of the artificial neural network (ANN) for 3 biomarkers. The colors of the lines represent positive and negative contributions, with red representing a positive contribution and gray a negative contribution. (F) Importance of the independent variables on the predicted results of the model. The horizontal coordinate is the hub gene, and the vertical coordinate represents the importance of the gene to the model outcome. (G and H) Confusion matrix (G) and ROC curve (H) for the three biomarkers.

Single-Gene GSEA and Immune-Related Analyses of Biomarkers

After constructing the ANN network, we analyzed the function of the biomarkers. The enrichment results indicated that the biomarkers were mainly enriched in KEGG terms such as neuroactive ligand-receptor interaction and olfactory transduction (Figure 6A–C, Supplementary Table 16–18). We further analyzed the differences in the proportions of immune factors between the CRSwNP and control groups. In total, one immune cell type (monocytes), three chemokines (CCL18, CCL28, and CXCL17), and one MHC (HLA-DOA) were significantly different between the CRSwNP and control groups (Figure 6D–F). Moreover, correlation analysis of biomarkers with differential immune factors revealed that HMGA2 had the strongest negative association with CXCL17, whereas IGF2BP1 exhibited the strongest positive relationship with HLA-DOA (Figure 6G).

Figure 6 Single-gene Gene Set Enrichment Analysis (GSEA) and immune-related analyses of the three biomarkers. (A–C) Single-gene GSEA of EZH2 (A) HMGA2 (B), and IGF2BP1 (C). (D) Differences in the abundance of immune cells between the CRSwNP and control groups. *p < 0.05. (E) Differences in chemokine expression between the CRSwNP and control groups. *p < 0.05, **p < 0.01. (F) Differences in MHC expression between the CRSwNP and control groups. *p < 0.05. (G) Correlation analysis of biomarkers with differential immune factors. (H) Expression of the three biomarkers in the training dataset. *p < 0.05, **p < 0.01. (I) Expression of three biomarkers in the GSE136825 dataset. **p < 0.01, ****p < 0.0001.

Abbreviation: ns, not significant.

Establishment of TF-mRNA Regulatory Networks, Regulatory Axes, and Diseases Associated with Biomarkers

To identify mRNA-TF networks of biomarkers in CRSwNP, we used the online databases ChEA3 and NetworkAnalyst, identifying a total of 106 TFs associated with the biomarkers. Subsequently, we established an mRNA-TF regulatory network containing 109 nodes (3 biomarkers) and 106 TFs (including CCND1, E2F6, OLIG2, and TEAD4) (Supplementary Figure 4A). In the DGIDB database, we identified one biomarker (EZH2) targeted by seven therapeutic drugs (doxorubicin, vorinostat, selumetinib, MK-2206, CPI-1205, tazemetostat, and dabrafenib) (Supplementary Figure 4B).

The ceRNA sub-network of biomarkers contained 3 mRNAs (EZH, HMGA2, and IGF2BP1), 12 miRNAs (including hsa-mir-485-5p, hsa-mir-98-5p, and hsa-let-7c-5p) and 9 lncRNAs (including CASC2, CDKN2B-AS1, and DUXAP8) (Supplementary Figure 4C). Correlation analysis of biomarkers with ENT diseases revealed that both IGF2BP1 and EZH2 were highly correlated with allergic rhinitis, while HMGA2 was highly associated with nose neoplasms (Supplementary Figure 4D).

Validation of the Expression of the Biomarkers

Expression validation results from the TS-CRSwNP dataset indicated that the three biomarkers significantly differed between control and CRSwNP tissues, with high expression levels in the CRSwNP group (Figure 6H). The results obtained using the validation cohort were consistent with those observed in the TS-CRSwNP dataset, except for IGF2BP1 (Figure 6I).

Experimental Verification

We further assessed the protein expression levels of the three biomarkers (IGF2BP1, EZH2, and HMGA2) in samples obtained from 20 patients with CRSwNP and 20 control participants using IHC analysis. The results indicated that IGF2BP1 expression was predominantly located in the cytoplasm, EZH2 was mostly expressed in the cytoplasm, while HMGA2 was expressed both in the cytoplasm and the nucleus. A subsequent quantitative analysis revealed that the number of cells positive for EZH2 (p < 0.01; Figure 7A), HMGA2 (p < 0.0001; Figure 7B), and IGF2BP1 (p < 0.001; Figure 7C) expression was significantly higher in the CRSwNP group than in the control group.

Figure 7 Quantification of EZH2, HMGA2, and IGF2BP1 expression. (A) Representative immunohistochemical (IHC) staining images and quantification of EZH2 levels in the control and CRSwNP group (Student’s t-test, n =20 per group). The boxes show high-power images (×60) obtained from low-power images (×20). The black arrows indicate EZH2. The scale bar on the left image is 200 μm, and the scale bar on the right image is 50 μm. **p < 0.01. (B) Representative IHC staining images and quantification of HMGA2 levels in the control and CRSwNP group (Student’s t-test, n =20 per group). The boxes show high-power images (×60) obtained from low-power images (×20). The black arrows indicate HMGA2. The scale bar on the left image is 200 μm, and the scale bar on the right image is 50 μm. ****p < 0.0001. (C) Representative IHC staining images and quantification of IGF2BP1 levels in the control and CRSwNP groups (Student’s t-test, n =20 per group). The boxes show high-power images (×60) obtained from low-power images (×20). The black arrows indicate EZH2. The scale bar on the left image is 200 μm, and the scale bar on the right image is 50 μm. **p < 0.01. (D) Expression of EZH2 in the control and CRSwNP groups (Student’s t-test, n =15 per group). *p < 0.05. (E) Expression of HMGA2 in the control and CRSwNP groups (Student’s t-test, n =15 per group). **p < 0.01. (F) Expression of IGF2BP1 in the control and CRSwNP groups (Student’s t-test, n =15 per group). *p < 0.05.

We then performed qRT-PCR to validate our findings by assessing the mRNA expression levels of IGF2BP1, EZH2, and HMGA2 in 15 CRSwNP samples and 15 control samples. The results revealed that the mRNA expression levels of EZH2 (p < 0.05; Figure 7D), HMGA2 (p < 0.01; Figure 7E), and IGF2BP1 (p < 0.01; Figure 7F) were significantly higher in the CRSwNP group than in the control group.

To further elucidate the molecular mechanism by which EZH2 regulates HMGA2 expression, we downregulated EZH2 expression. Specifically, we synthesized three siRNAs targeting EZH2 mRNA sequences and transfected them into HNEpCs. As presented in Figure 8A, siEZH2-1 exhibited more potent inhibition, resulting in a more pronounced downregulation of EZH2 expression. We measured the mRNA levels of H3k27me3 and HMGA2 using qRT-PCR. Following the downregulation of EZH2 expression, the expression level of H3k27me3 increased, whereas that of HMGA2 decreased in HNEpCs (Figure 8B and C). We further used IF to detect changes in the protein expression profile of H3k27me3 and HMGA2 after the downregulation of EZH2. We semi-quantitatively scored the immunofluorescence results to assess the correlation between reduced EZH2 expression levels and the expression of H3k27me3 and HMGA2. The figure presents representative images of H3k27me3 and HMGA2 expression after the downregulation of EZH2 expression in HNEpCs, revealing that this downregulation resulted in increased H3k27me3 expression levels (Figure 8D) and decreased HMGA2 expression levels (Figure 8E).

Figure 8 Immunofluorescence analysis of EZH2, HMGA2, and IGF2BP1 expression based on three different siRNAs targeting EZH2 mRNA sequence. (A) Expression of EZH2 in the control group and HNEpCs transfected with three different siRNAs targeting EZH2 mRNA sequences. *p < 0.05, **p < 0.01. (B) Expression of H3k27me3 in the control group and HNEpCs transfected with siEZH2. ***p < 0.001. (C) Expression of HMGA2 in the control group and HNEpCs transfected with siEZH2. ****p < 0.0001. (D) Representative immunofluorescence images and quantified data showing the expression of H3k27me3 in the control group and HNEpCs transfected with siEZH2. Nuclei were stained with DAPI. Scale bar = 50 μm. *p < 0.05. (E) Representative immunofluorescence images and quantified data showing the expression of H3k27me3 in the control cell and HNEpCs transfected with siEZH2. Nuclei were stained with DAPI. Scale bar = 50 μm. **p < 0.01.

To assess the effects of the upregulation of EZH2 expression, we synthesized three plasmids overexpressing EZH2 and transfected them into HNEpCs. Figure 9A demonstrates that the OE-EZH2 plasmid exhibited a stronger effect, resulting in more significant EZH2 overexpression than that induced by the two other plasmids. We used qRT-PCR to detect the mRNA levels of H3k27me3 and HMGA2 in response to EZH2 upregulation. Following the upregulation of EZH2 expression, the expression level of H3k27me3 was significantly decreased, whereas that of HMGA2 was increased (Figure 9B and C). We used IF to identify changes in the protein expression patterns of H3k27me3 and HMGA2 following the upregulation of EZH2. In HNEpCs with upregulated EZH2 expression, H3k27me3 expression levels were decreased (Figure 9D) whereas those of HMGA2 were increased (Figure 9E). These results confirm that EZH2 promotes HMGA2 expression by downregulating H3k27me3 expression and inhibits HMGA2 expression by upregulating H3k27me3 expression in HNEpCs.

Figure 9 Immunofluorescence analysis of EZH2, HMGA2, and IGF2BP1 expression based on three different plasmids overexpressing (OE) EZH2. (A) Expression of EZH2 in the control group and HNEpCs transfected with three different plasmids overexpressing (OE) EZH2. **p < 0.01. (B) Expression of H3k27me3 in the control group and HNEpCs transfected with OE-EZH2. ***p < 0.001. (C) Expression of HMGA2 in the control group and HNEpCs transfected with OE-EZH2. ****p < 0.0001. (D) Representative immunofluorescence images and quantified data showing the expression of H3k27me3 in the control group and HNEpCs transfected with OE-EZH2. Nuclei were stained with DAPI. Scale bar = 50 μm. ***p < 0.001. (E) Representative immunofluorescence images and quantified data showing the expression of H3k27me3 in the control group and HNEpCs transfected with OE-EZH2. Nuclei were stained with DAPI. Scale bar = 50 μm. *p < 0.05.

Discussion

In the present study, we performed RNA-seq using CRSwNP and control tissues, enabling the identification of candidate genes and pathways associated with the disease. Pathway analysis revealed key signaling pathways involved in the pathogenesis of CRSwNP. Comparison of CRSwNP tissues with control ones revealed alterations in key signaling pathways such as metal ion transmembrane transporter activity and MAPK signaling pathway, which are associated with chronic airway inflammation.39,40 The identified DEGs, target genes of DE-miRNAs, and DE-target genes of lncRNAs were also enriched in pathways related to metal ion transmembrane transporter activity and MAPK signaling. Notably, we identified three biomarkers (IGF2BP1, EZH2, and HMGA2) that are closely associated with CRSwNP.

IGF2BP1 is a key regulator of mRNA metabolism and transport during organismal development. It is primarily expressed in embryos, with relatively low expression in adults and healthy tissues.41 Mechanistically, IGF2BP1 binds to the methyl group of RNA to stabilize it and prevent its degradation, thereby promoting protein translation. Functionally, IGF2BP1 plays an important role in regulating development and cell proliferation.42 Consistent with previous studies, we observed that IGFBP1 was highly expressed in the CRSwNP group compared to the control group. The observed differences in IGF2BP1 expression between the training and validation sets in this study may stem from multiple factors. Variations in sample collection, processing, and storage protocols, as well as differences in experimental equipment, could all contribute to expression level discrepancies. Additionally, individual variations such as age, sex, disease status, and lifestyle factors may influence IGF2BP1 expression, and the demographic differences between subjects in the training and validation sets might have amplified this effect.

The EZH2 gene is the human homolog of the Drosophila zeste gene enhancer and belongs to the family of epigenetic regulators.43 EZH2 performs two main biological functions, depending on polycomb group genes (PcGs). Specifically, it plays a role in various biological processes, and its expression is associated with cancer onset, progression, metastasis, metabolism, drug resistance, and immune modulation. Therefore, EZH2 has become a potential target for immunotherapy.44 In allergic rhinitis, EZH2 can promote Th1 cell differentiation or inhibit Th2 cell differentiation, thereby exacerbating the pathological process of the disease.27 Furthermore, EZH2 may play a pro-inflammatory role in allergic airway inflammation, as inhibiting its expression can alleviate allergic airway inflammation.26 Consistent with findings in other immune-related diseases, in our study, EZH2 was highly expressed in nasal polyp tissues, suggesting that EZH2 plays an important role in the pathogenesis of CRSwNP. In general, EZH2 expression is suppressed through epigenetic modifications, leading to the inhibition of the transcription of target genes, a mechanism that has received considerable attention in tumor research.45 In this study, we focused on the role of EZH2 in NPs. EZH2 is involved in the regulation of inflammation,46 potentially playing a role in the pathophysiology of chronic sinusitis with NPs. Notably, inhibiting EZH2 significantly reduces inflammation and lung injury in acute lung injury (AKI) induced by sepsis.47 EZH2 inhibition can alleviate AKI induced by ischemia-reperfusion and reduce the associated inflammatory responses.48

HMGA2 is a transcriptional regulator essential for embryonic development, with its expression decreasing during postembryonic stages.49 However, it is often re-expressed during tumorigenesis.50 RNA-sequencing analyses of clinical and normal samples indicate that HMGA2 expression is upregulated in many cancers.51 HMGA2 induces cancer cell proliferation by promoting the transition into the S phase of the cell cycle, inhibiting apoptosis, and inducing DNA repair. In addition, HMGA2 regulates miRNA expression by activating multiple signaling pathways, including the MAPK/ERK, TGFβ/Smad, PI3K/AKT/mTOR/NFKB, and RKIP pathways.49 In our study, EZH2 and HMGA2 were highly expressed and synergistic, suggesting that these two key genes may collaboratively play an important role in the pathogenesis of CRSwNP.

Epigenetic abnormalities are important factors contributing to the development of various diseases.52 Histone modifications in chromosomes, which serve as key markers of epigenetic inheritance, are critical in regulating the transcription of disease-associated genes by controlling their expression.24 EZH2, a prominent epigenetic target, is a histone methyltransferase that can affect cellular functions by catalyzing the trimethylation of histone H3 lysine 27 (H3K27). Its expression has been detected in HNEpCs.53 In our study, we found that, similar to findings in tumor tissues, EZH2 overexpression in HNEpCs decreased the methylation level of H3K27, whereas EZH2 downregulation increased the methylation level of H3K27. Recent advances in immunotherapy have highlighted the importance of the immune system in disease control, but only a minority of patients currently benefit. Studies have shown that targeting epigenetic factors that inhibit immune cell function, such as EZH2 inhibition, can improve the efficacy of immunotherapy by reshaping the disease microenvironment and enhancing the coordination of the immune system.54,55 In addition, lncRNA PiHL inhibits its localization in HMGA2 promoter region by binding to EZH2, thereby reducing H3K27me3 methylation level and promoting HMGA2 expression.56 In our study, we found that HMGA2 exhibited a strong synergistic effect with EZH2. Specifically, when EZH2 was overexpressed, H3K27me3 expression was downregulated whereas HMGA2 expression was upregulated. In contrast, downregulated EZH2 expression resulted in increased H3K27me3 levels and decreased HMGA2 expression. These results suggest that targeting EZH2 may pave the way for enhancing immune surveillance in CRSwNP. We also found that EZH2 and HMGA2 were co-regulated by hsa-let-7c-5p. We speculated that hsa-let-7c-5p may promote the high expression of EZH2 through histone modification, and then activate the expression of HMGA2 in CRSwNP through a more complex regulatory mechanism, revealing the potential interaction between EZH2 and HMGA2. It is now increasingly understood that CRS phenotypes do not necessarily match these classic endotypic presumptions, and the underlying pathogenetic processes in CRS may encompass type 1, type 2, or type 3 inflammation, or a mix. However, from the perspective of epigenetics, genetic predisposition is stronger in patients with CRS with nasal polyps compared with those without nasal polyps (CRSsNP).57 Among the many epigenetic regulatory modalities, methylation modification of histone H3 at its K27 position mediates the silencing of gene expression, and EZH2 is a key enzyme molecule that induces methylation modification of H3K27.58 Furthermore, the three biomarkers (IGF2BP1, EZH2, and HMGA2) demonstrate significant correlations with the differential chemokine CCL18. Studies have revealed that CCL18 is highly expressed in M2 macrophages and aligns with the immune characteristics of type 2 inflammatory endotypes, suggesting its potential critical role in type 2 inflammation.59 Dysregulation of EZH2 is closely associated with CRSwNP. HMGA2 could synergize with EZH2 to induce epithelial cell destruction; therefore, the study of the immunomodulatory role of EZH2 and the related epigenetic regulatory mechanism is of great significance for the prevention and diagnosis of related immune disease. We further hypothesize that IGF2BP1, EZH2, and HMGA2 may play pivotal roles in type 2 inflammatory endotypes of CRSwNP through epigenetic regulatory mechanisms. These genes likely interact and modulate chemokine expression (eg, CCL18) via their coordinated actions in specific immune cells, thereby promoting the initiation and perpetuation of type 2 inflammatory responses.

The EZH2 signaling pathway, known as a driver of cancer progression, can potentially interfere with anti-tumor immune responses, making it an important therapeutic target for cancer therapy.60 In this study, seven small molecule drugs were selected based on the prediction of EZH2. VORINOSTAT, as a histone deacetylase inhibitor, affects gene expression by regulating chromatin accessibility and has been widely used in the treatment of hematological malignancies.61 TAZEMETOSTAT is a potent targeted epigenetic regulator that specifically inhibits EZH2.