Integrated analysis reveals that EGR1 promotes epithelial IL33 production in T2 asthma

Abstract

Background

Airway epithelial cells constitute the first line of defense against external noxious stimuli and play crucial roles in the release of epithelial inflammatory cytokines (IL33, IL25 and TSLP), initiating airway allergic inflammatory diseases such as asthma. IL33 plays critical physiological processes in T2-endotype asthma. However, the mechanisms by which allergen exposure triggers IL33 release from airway epithelial cells remain unclear.

Methods

Integrated bioinformatic analysis and transcriptional analysis of bulk RNA-seq and single cell RNA-seq (scRNA-seq) data were used to identify core genes and determine the internal gene network associated with IL33. The expression of EGR1 was subsequently analyzed in vitro in the BEAS-2B cell line and in vivo in a house dust mite (HDM)-induced mouse asthma model. The functional experiments of EGR1 were investigated in vitro via siRNA knockdown and over-expressed plasmid. Chromatin immunoprecipitation (ChIP)-PCR and dual-luciferase reporter assay validation were subsequently performed to investigate the mechanisms by which EGR1 regulates IL33 secretion.

Results

Bulk RNA-seq and scRNA-seq data identified EGR1 as an epithelial cell-derived gene implicated in IL33 expressions in asthma. The comprehensive analysis of multiple datasets indicated that the high EGR1 expression in epithelial cells may suggest a mechanistic basis of T2-endotype childhood asthma. Moreover, we verified that the expressions of EGR1 in airway epithelial cells were elevated both in vitro and in vivo asthma models. EGR1 regulated the production of IL33. Ultimately, ChIP and luciferase reporter assays confirmed that transcription factor EGR1 directly regulate the transcription of IL33 mRNA.

Conclusions

Our integrated bioinformatic analysis elucidated that EGR1 directly regulates the production of IL33 in T2-asthma and provide insights underlying the progression of asthma.

Introduction

Asthma is a chronic respiratory disease with global health implications most common in children [1]. Asthma can be divided into high- and low-T2 subtypes according to its gene expression profile, and the high-T2 inflammatory subtype (T2-endotype) is the predominant subtype [2]. In the initiation of T2 inflammation, airway epithelial cells (AECs) play important roles in recognizing allergens and releasing epithelial cytokines [3]. AECs and their epithelial cytokines are the key drivers of allergen sensitization [4]. After exposure to external stimuli, dysfunctional AECs release epithelial inflammatory cytokines such as IL25, IL33 and TSLP. IL33 activates ST2⁺ mast cells, eosinophils, ILC2s and Th2 cells to produce T2 inflammation cytokines and promote T2 immune response [5, 6]. Although IL4, IL5, and IL13 are classic T2-cytokines, they are expressed at low levels in AECs [2, 7]. Recent studies [2, 8] have demonstrated that POSTN, CLCA1, and SERPINB2 could be used as epithelial genic markers of T2-endotype asthma. Based on the expressions of these genes, patients can be divided into high-T2 and low-T2 subtypes.

The airway epithelium participates in complex interactions with multiple immune cells [7, 9]. Following the stimulation of allergens, AECs produce a wide range of cytokines. IL33 is one of the important epithelial cytokines that can activate downstream signaling pathways, leading to the progression of T2-endotype asthma. A meta-analysis of a genome-wide association study of asthma elucidated susceptibility loci near IL33 and ST2/IL1RL1 [10]. Higher levels of IL33 expressions were detected in severe asthma patients [11]. The evidence supported the predominant contributor roles of IL33.

Epithelial IL33 is a key characteristic cytokine of T2-endotype asthma and is regulated in AECs [12]. However, the mechanisms involved in the upstream regulation of IL33 remain unclear [13]. According to the pathological mechanisms of allergic asthma, we used house dust mite (HDM), the classical inducers of allergic asthma, to stimulate BEAS-2B epithelial cells and detected the production of epithelial inflammatory cytokines reflective of T2-endotype asthma in vitro. In addition, we explored key transcription factors and core genes, which are important for clarifying the pathogenesis of T2-endotype asthma and identifying potential biomarkers for asthma diagnosis and drug development.

Methods

RNA-seq dataset analyses

A graphical overview of the integrated bioinformatic analysis study design is shown in Fig. S1. The GEO database was searched with the keywords “children” and “asthma”, and the identified transcriptome sequencing bulk RNA-seq datasets GSE18965 [14], GSE118761 [15], GSE145505 [16], and GSE152004 [17] with airway samples in childhood asthma or wheezing diseases, were downloaded. Cluster analysis was used to categorize each sample in the datasets as high T2 or low T2 according to the expression levels of the T2 inflammatory marker genes POSTN, CLCA1, and SERPINB2 [2]. Heatmaps were generated using ComplexHeatmap (v 2.14.0) and Circlize (v0.4.15).

To obtain the differentially expressed genes (DEGs) associated with high and low T2, the RNA-seq data were analyzed via the R software packages limma (v3.54.2) and DESeq2 (v1.38.3), respectively, and the DEGs were defined as those genes with |logFC|>0.25 and p-value < 0.05. The DEGs were visualized in volcano maps. The intersection of the DEGs of the 4 datasets was analyzed and visualized with the UpSetR package (v1.4.0).

To select the feature genes, two datasets (GSE145505 and GSE152004) were used for feature screening of DEGs using a random forest approach with the randomForest package (v4.7-1.1) and LASSO methods with the glmnet package (v4.1-8) in R. The intersection of the feature-screened genes was selected and visualized via the UpSetR package. The characteristics of the datasets are shown in Table S1.

scRNA-seq dataset analyses

Two single-cell RNA sequencing (scRNA-seq) datasets of GSE164015 and GSE193816 with airway samples of asthma patients were downloaded. We acquired the scRNA-seq GSE164015 [18] comprising more than 35,140 human single cells from bronchoalveolar lavage fluid samples of airway epithelial brushings obtained via bronchoscopy from 4 independent asthma patients after lung segmental allergen challenge with diluent or allergen challenge. Six paired samples of 3 patients with allergen challenge of diluent vs. HDM were selected into analysis, instead of the remaining two samples of 1 patients with allergen challenge of diluent vs. cat. Another scRNA-seq dataset, GSE193816 [19], comprising endobronchial brushing samples via bronchoscopy from 4 allergic asthma patients and 4 allergic controls was analyzed. The participants had an allergic history of allergic rhinitis and/or conjunctivitis to HDM/cat with a positive skin prick test to the same allergen.The characteristics of the datasets are shown in Table S1.

Genes expressed in fewer than three cells in a sample, cells that expressed fewer than 200 genes or more than 5500 genes, and cells with a mitochondrial gene content > 25% or hemoglobin gene content > 20% of the total unique molecular identifier count were excluded. Next, the data were normalized using the default parameters of the log normalization method, and the top 2000 highly variable genes were identified using the FindVariableFeatures function. Z-score transformation of gene expression was performed using the ScaleData function. Subsequent dimensionality reduction analysis was performed using the RunPCA function. Clustering was based on the 20 most significant principal components (PCs). The effect of different resolutions on the clustering was evaluated. Finally, uniform manifold approximation and projection (UMAP) visualizations were constructed using the same number of PCs as the associated clustering.

We focused on epithelial cells. After obtaining the cell clusters, the FindAllMarkers function was used to investigate the DEGs in the clusters. We set the log fold change in average expression between the clusters (avg_logFC) to > 0.25 and the p-value after the Bonferroni correction to < 0.05. Genes were ranked by ascending p-value to evaluate differential expression in each cluster, and cell types were manually identified using airway marker genes from previous studies [20,21,22] and CellMarker [23]. Furthermore, core genes were identified as the genes that intersected DEGs from bulk RNA-seq and critical epithelial DEGs from scRNA-seq.

Transcription factors analysis and ROC analysis

Transcription factor analysis via single-cell regulatory network inference and clustering (SCENIC) was performed via the R software Scenic package (v1.3.1) and Python software Pyscenic v0.12.0 (VIB Center for Brain Disease Research, Laboratory of Computational Biology, Leuven, Belgium) [24]. To analyze transcription factor regulation, SCENIC consists of three steps: GENIE3 (v1.24.0) inference of gene coexpression networks, RcisTarget (v1.22.0) analysis of transcription factor-binding motifs and construction of transcription factor-target gene networks, and AUCell (v1.24.0) identification of cells with active gene sets and calculation of regulatory activity.

The ssGSEA algorithm of the R software GSVA package (v1.46.0) was used for immune infiltration analysis of the core genes and immune cells. To evaluate the predictive power of genes associated with childhood asthma, a receiver operating characteristic (ROC) curve was plotted using the R software pROC package (v1.18.5), and the area under the curve (AUC) was calculated to assess the prediction accuracy of childhood asthma. AUC values > 0.8, 0.6–0.8, and 0.5–0.6 represent high, medium, and low differentiation degrees, respectively, and closer to 1 indicates a better predictive value.

Asthma model

Eight-week-old female C57BL/6JGpt mice were purchased from the Laboratory Animal Center of Chongqing Medical University in China and maintained in a controlled environment at ∼23 °C with 40–50% humidity under a 12/12 h dark/light cycle. The mice were randomly divided into a control group and an asthma group, with five mice per group. Each mouse in the asthma group was intranasally administered each time a dose of 20 µg/40 µl HDM (D. pteronyssinus, Greerlabs, USA) in the sensitization stage (day 0 and day 14) and in the challenge stage (day 21, 23, 25, 27 and 29). The control group received the same volume of 0.9% saline at the same time points. Lung function examination was performed on day 30. The mice were sacrificed on day 31, and lung tissue samples were collected. In vitro, a human airway epithelial cell line (BEAS-2B) was incubated with or without HDM extract (30 µg/mL) for 3 h to construct asthma and control cell models.

Cell culture and transfection

The BEAS-2B cells were cultured in DMEM (Gibco, USA) containing 10% fetal bovine serum (Gibco, USA) and 1% penicillin‒streptomycin solution (NCM Biotech, China) at 37°C with 5% CO₂. Cells were seeded at a density of 1.5 × 10⁵ cells per well in 12-well plates. Transient transfection was performed once cells reached 60–70% confluence using Lipofectamine™ 2000 (Invitrogen, USA) with small interfering RNAs (siRNAs) (Sangon Biotech, Shanghai, China) for the EGR1 knock-down and pcDNA3.1 (+) plasmid (Genecreate, Wuhan, China) for the EGR1 over-expression vector, following the manufacturers’ standard protocols. The following siRNA sequences were used: EGR1 sense 5’-GCAGCAGCAGCACCUUCAATT-3’ and antisense 5’-UUGAAGGUGCUGCUGCUGCTT-3’. After transfection with siRNA (siRNA 50 nM and Lipo2000 2 µl per well) for 48 h or with the over-expression plasmid (plasmid 1 µg and Lipo2000 2 µl per well) for 48 h, the cells were treated with or without HDM (30 µg/mL) for 3 h. All experiments were repeated as least three times.

Real-time quantitative PCR

The total RNA of the cells and tissues was extracted with Total RNA Extraction Reagent (Bioflux, China). An Evo M-MLV RT Reaction Mix Kit (Accurate Biotechnology, China) was used to reverse transcribe the extracted RNA into cDNA, according to the manufacturers’ protocols. Then, real-time quantitative polymerase chain reaction (qPCR) was conducted with a SYBR Green Premix Pro Taq HS qPCR Tracking Kit (Accurate Biotechnology, China) on a Bio-Rad CFX PCR System (Bio-Rad, USA). Reactions were carried out in a total volume of 10 µL, including 5 µL of SYBR Green Pro Taq HS, 0.2 µL of each specific primer, 2.6 µL of RNase-free water, and 2 µL of cDNA. The primer sequences are listed in Table S2. The PCR reactions were performed for 40 cycles at 95 °C for 10 s and 61.4 °C for 30 s.The relative expression of EGR1 and IL33 was normalized to that of GAPDH (internal reference) via the 2^−ΔΔCt method.

Western blotting

The cells were harvested, and the denatured protein samples were resolved on 10% SDS‒polyacrylamide gels, electrophoresed and then wet-transferred onto polyvinylidene difluoride membranes. The membranes were blocked with 5% bovine serum albumin at room temperature and incubated with primary antibodies overnight at 4 °C, followed by incubation with a horseradish peroxidase-conjugated secondary antibody at room temperature for 1 h. Proteins were visualized using ECL reagents and a Bio–Rad gel imaging system. The primary antibodies used in the study included an anti-EGR1 rabbit monoclonal antibody (1:1000 dilution; Cell Signaling Technology, USA) and an anti-beta-actin mouse polyclonal antibody (1:5000 dilution; ZenBio, China). The secondary antibodies used were an anti-goat rabbit antibody (1:8000 dilution; ZenBio, China) or an anti-goat mouse antibody (1:8000 dilution; ZenBio, China).

Immunofluorescence

Tissue sections fixed in paraffin (4 μm in thickness) were rehydrated and dewaxed. Antigen repair was performed, and the tissues were subsequently incubated at 4 °C overnight with anti-EGR1 antibodies (1:200 dilution, Cell Signaling Technology, USA). On the next day, the sections were incubated with fluorescent secondary antibodies (Donkey anti-Rabbit IgG H + L, Alexa Fluor™ 488, 1:500 dilution, Thermo Scientific, USA) and sealed with antifade mounting medium with DAPI (Beyotime, China). Images were then taken (magnification: ×20 or ×60) with a confocal microscope (Nikon C2⁺).

Enzyme-linked immunosorbent assay (ELISA)

The levels of IL33 in the cell culture supernatants were measured via ELISA with commercial IL33 (MeiMian, China) kits according to the manufacturer’s instructions.

Chromatin immunoprecipitation (ChIP) and dual luciferase assay

In accordance with the instructions of the Chromatin IP Kit (Waals™, China), the cells were fasten in 1% formaldehyde (Thermo Scientific, USA) at 37 °C for 6 min, and the reaction was terminated by treatment with glycine solution for 5 min. Then, cell cross-linking was completed, nuclei were isolated, and the chromatin was fragmented using a Bioruptor instrument (Diagenode). Chromatin was precipitated with normal Rabbit IgG (10uL) (Cell Signaling Technology, USA) or anti-EGR1 (10 µL) with protein A/G magnetic beads. Each sample was evaluated by triplicate qPCR. The Ct value of each ChIP DNA fragment was standardized to that of the input DNA fragment detected by the same qPCR assay (ΔCt) to illustrate the difference in chromatin sample preparation. The %input for each ChIP fraction was calculated as 2^[Ct(input)-Ct (ChIP)]*100%, and anti-IgG fold enrichment was evaluated. The sequences of the ChIP primers used are presented in the supplemental material (Table S3).

For the dual-luciferase reporter assay, cells were seeded at a density of 0.8 × 10⁵ cells per well in 24-well plates, and a human IL33 full-length promoter (-2000 ~ -1) luciferase reporter plasmid (Genecreate, Wuhan, China) and EGR1 expression plasmid were transiently cotransfected into the cells. The luciferase activity was measured using the dual-luciferase reporter assay system and normalized to that of Renilla luciferase TK-pRL.

Statistical analyses

Transcriptomic data from bulk RNA-seq and scRNA-seq were analyzed via R version 3.5.1. Unless otherwise stated, statistical analyses of the experiments were performed in GraphPad Prism v8.0 (GraphPad Software, San Diego, CA, USA). Data were presented as the mean ± SD for the number of samples. The unpaired Student’s t-test, one-way ANOVA, or repeated-measures ANOVA as appropriate for each set of experimental conditions was used to evaluate the differences between groups. Differences were considered statistically significant at *p < 0.05, **p < 0.01, ***p < 0.001, or ****p < 0.0001.

Results

Identifying the key genes associated with high-T2 inflammation in childhood asthma via bulk RNA-seq

To identify the key genes associated with high-T2 inflammation in childhood asthma, four bulk RNA-seq datasets were analyzed. The samples were clustered into high-T2 and low-T2 subgroups on the basis of the expression of POSTN, CLCA1, and SERPINB2 and visualized in a circular heatmap (Fig. 1A-D). The high-T2 and low-T2 subgroups were compared to obtain DEGs, which were visualized via a volcano map (Fig. 1A-D). 43 co-DEGs were obtained from the intersection of the DEGs from the four datasets (Fig. 1E).

Fig. 1

Key genes were identified to be associated with high-T2 inflammation in childhood asthma by bulk RNA-seq: (A-D) DEGs identification of T2-high vs. T2-low children: GSE18965, GSE118761, GSE145505, GSE152004; (E) Upset plot of 43 co-DEGs by the intersection of four datasets. (F) The feature genes selection by Lasso in GSE145505; (G) The feature genes selection by Lasso in GSE152004; (H) Upset plot of co-DEGs from the gene selection by the randomForest and Lasso in GSE145505 and GSE152004

To select the feature genes associated with high-T2 inflammation in childhood asthma patients among the 43 co-DEGs, we conducted machine learning with the Lasso and random forest methods in the GSE145505 and GSE152004 datasets. Lasso regression with lambda.min was used to screen feature genes (Fig. 1F‒G). The random forest algorithm ranked the importance of the genes according to the mean decrease in the Gini coefficient (Fig. S2 A-D). The intersection of the above feature screening process was subsequently performed, and the results identified ten genes (EGR1, SERPINB2, STEAP4, MAGT1, POSTN, ACE2, CFH, PTHLH, PRKD1, and SCNN1B) (Fig. 1H) as key genes closely related to high-T2 inflammatory asthma in airway specimens from children. However, airway specimens contain numerous cell types, including epithelial cells and various immune cells. Bulk RNA-seq cannot distinguish precisely the exact cell subtype from which these genes are derived. Consequently, we considered that scRNA-seq might provide additional useful information.

Identifying the epithelial cell-derived genes associated with epithelial inflammation in the T2 endotype via scRNA-seq

Two sets of scRNA-seq data were analyzed to determine the exact cellular origin of these epithelial cell-derived genes. In the scRNA-seq of GSE164015, a total of 35,140 cells and 23,294 genes were harvested (Fig. S3 A-B), and a total of 12 cell types were identified (Fig. S4), among which the largest number of cell types were epithelial cells. The fine classification (12 cell types) and the rough classification (epithelial cells and immune cells) of the cell types were visualized using UMAP (Fig. 2A-B). The marker genes of the cell types are shown in Fig. 2C. The epithelial cell subtype-specific DEGs between the HDM-incubated group and control group were identified via FindMarkers (Fig. S5), revealing a total of 1978 upregulated and downregulated genes; the volcano map of these DEGs is shown in Fig. 2D.

Fig. 2

Epithelial cell-derived genes were identified by scRNA-seq: (A-B) UMAP of clusters of rough label and final label in GSE164015; (C) The visualization of marker genes by clusters of GSE164015; (D) DEGs of HDM and CON groups of epithelial cells in GSE164015; (E-F) UMAP of clusters of rough label and final label in GSE193816; (G) The visualization of marker genes by clusters of GSE193816; (H) DEGs of asthma and non-asthma groups of epithelial cells in GSE193816

The GSE193816 scRNA-seq dataset yielded 28,604 cells and 52,152 genes (Fig. S6 A). These cells were categorized into 19 cell types (Fig. S6 B), among which the largest number of cell types were epithelial cells (Fig. 2E-G). The same method was used to analyze the DEGs between different cell types (Fig. S7). Analysis of the epithelial cell subtype-specific DEGs between the allergic asthma group and the allergic nonasthma group revealed a total of 1308 upregulated and downregulated genes (Fig. 2H).

Further integrated analysis intersecting the bulk RNA-seq and scRNA-seq DEGs revealed that EGR1, SERPINB2, STEAP4, and SCNN1B were core epithelial cell-derived genes associated with T2 inflammation in asthma (Fig. S8).

To explore the associations between the four core genes and T2-related inflammatory cytokines, we performed correlation analysis on different datasets of childhood asthma patients. Since SERPINB2 has been shown to be a marker of T2 inflammation in asthma [8, 25], the remaining three genes (EGR1, STEAP4, and SCNN1B) were analyzed in the GSE18965 and GSE118761 datasets (Fig. 3A). Relatively consistent results were obtained in the different childhood asthma datasets, indicating that EGR1 levels were significantly positively correlated with IL33 production (Fig. 3A-C). There was no significant correlation between EGR1 and IL25 or TSLP expression.

Fig. 3

Transcription factor EGR1 may be associated with epithelial inflammatory IL33 production in T2 asthma: (A) The correlation of the key gene (EGR1) and T2 markers; (B-C) The correlation between EGR1 and IL33 in GSE18965 and GSE118761; (D-E) Transcription factor SCENIC analysis by scRNA-seq; (F-H) ROC curve of the predictive accuracy for EGR1 in childhood asthma in GSE18965, GSE118761 and GSE145505

Additionally, to analyze the associations between genes and the immune microenvironment in high-T2 asthma patients, immune infiltration gene set variation analysis (GSVA) was conducted. These results revealed that EGR1 levels were positively correlated with CD4⁺ T cells, eosinophils, and mast cells in all the datasets except GSE118761 (Fig. S9). These findings further indicate that EGR1 may play an important role in immunity in T2 asthma. Therefore, among the four core genes, we subsequently focused on EGR1 and especially its association with IL33 production.

Transcription factor analysis reveals that EGR1 is associated with epithelial inflammatory IL33 production in T2 asthma

To further explore the upstream transcription factors in the epithelial cell type associated with T2-related inflammatory asthma, we performed transcription factor regulatory network analysis of scRNA-seq via the R package SCENIC and used the HDM stimulation group in the GSE164015 dataset to predict potential regulatory factors in cell types. The prediction results for SCENIC revealed that EGR1 is a potential transcription factor that specifically regulates the core genes of epithelial cells (Fig. 3D-E). In addition, the JASPAR database shows predicted EGR1 binding sites within the promoter of the target gene IL33, suggesting that EGR1 could directly regulate the transcription of IL33 mRNA.

To further investigate the role of EGR1 in childhood asthma, we evaluated datasets GSE18965, GSE118761, and GSE145505 to determine the predictive power of this gene for childhood asthma diagnosis. EGR1 had moderate predictive value for childhood asthma diagnosis, with AUCs of 0.78, 0.61, and 0.62 in these respective datasets (Fig. 3F-H), which supported the predictive efficacy of the core EGR1 gene in the diagnosis of childhood asthma.

EGR1 expressions are upregulated in the airway epithelium in an HDM-induced asthma model in vivo and in vitro

To explore the expressions of EGR1 after allergen stimulation, we conducted in vitro experiments in BEAS-2B and in vivo experiments in a mouse asthma model. In vitro, BEAS-2B cells were stimulated with HDM extract (30 µg/mL) for 3 h to construct a cellular asthma model. Results indicated that the EGR1 mRNA and EGR1 protein levels were markedly increased in the HDM-stimulated group (Fig. 4A-C), followed by the hypersecretion of IL33 (Fig. 4D-E). These results suggested that HDM extract exposure promoted the expression of EGR1 and IL33 in vitro.

Fig. 4

EGR1 expressions are upregulated in HDM-induced asthma model: (A-C) EGR1 expression levels of BEAS2B stimulated by HDM by qPCR, WB and IF (×60 magnification); (D) IL33 levels in cultural supernatant of BEAS2B; (E) IL33 mRNA expression of BEAS2B; (F) Mechacholine-induced airway resistance in different groups of asthma model mice in vivo; (G) Lung tissues with control and asthma by HE staining in vivo; (H-J) EGR1 expression levels of the asthma model mice in vivo by RT-PCR, WB and IHC (×20 magnification). Results were analyzed using unpaired Student’s t-test, except that data of (F) were analyzed using repeated-measures ANOVA. Differences were considered statistically significant at p < 0.05 (*p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001). All experiments were conducted with three replicates

In vivo, an asthma mouse model was constructed via treatment with HDM extract. The airway resistance in the asthma group was significantly greater than that in the control group (Fig. 4F). The pathological features of the asthma model mice included thickened bronchial walls, inflammatory infiltration of the submucosa, inflammatory cell infiltration of the lung interstitium and thickening of the lung interstitium (Fig. 4G). The levels of EGR1 mRNA and EGR1 protein in the asthma group were significantly greater than those in the control group (Fig. 4H-J), which supported the increased expression of EGR1 in the asthma model.

EGR1 regulates the levels of the epithelial inflammatory cytokine IL33

To experimentally verify whether EGR1 can regulate IL33 production in the HDM-stimulated asthma model, we performed plasmid-mediated overexpression (oe-EGR1) and siRNA-mediated knockdown (si-EGR1) experiments in vitro (Fig. S10 A-B). The results indicated that oe-EGR1 promoted the productions of IL33 (Fig. 5A-B) and that si-EGR1 suppressed the expressions of IL33 (Fig. 5C). These results suggest that EGR1 can regulate the expressions of epithelial IL33 in the HDM-stimulated cellular model.

Fig. 5

EGR1 regulates IL33 expressions: (A) The IL33 mRNA level in BEAS2B with the overexpression of EGR1 simulated by HDM; (B) The IL33 protein level in BEAS2B by ELISA with the overexpression of EGR1 simulated by HDM; (C) The IL33 mRNA level in BEAS2B with the knockdown of EGR1 simulated by HDM; (D) The binding sites predicted by JASPER database; (E-G) ChIP for the mechanism of EGR1 transcriptional regulation on IL33 promotor; (H) The dual-luciferase reporter assays for the mechanism of EGR1 transcriptional regulation on IL33 promotor. Results of (A, B, C, H) were analyzed using one-way ANOVA, and results of (F, G) were analyzed using Student’s t-test. Differences were considered statistically significant at p < 0.05 (*p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001). The experiments were conducted with three technical replicates except ChIP and dual-luciferase reporter assay

Furthermore, to delve deeper into the molecular mechanism by which EGR1 regulates IL33, we conducted JASPAR prediction, ChIP‒PCR and dual-luciferase reporter assays. The JASPAR database showed direct binding sites for EGR1 within the IL33 transcription initiation region, which indicated that EGR1 can increase gene transcription by binding to specific DNA sequences in the IL33 promoter region. Next, we designed primers for three of the predicted binding sites (-1548 bp to -1535 bp, -1125 bp to -1111 bp, and − 602 bp to -583 bp) for ChIP–PCR (Fig. 5D). The ChIP‒PCR results revealed dramatic enrichment of the IL33 promoter upon immunoprecipitation with the anti-EGR1 antibody relative to the control IgG, indicating the direct binding relationship of EGR1/IL33 (Fig. 5E‒G). Moreover, a dual-luciferase reporter assay (Fig. 5H) revealed that upregulation of EGR1 led to increased luciferase activity in cells transfected with the IL33 reporter construct. These results confirm the direct transcriptional regulation of epithelial EGR1 on IL33 promoter sequences.

Discussions

Asthma is chronic inflammatory airway disease with remarkable heterogeneity, and its development is influenced by both genetic and environmental factors [3]. Asthma is usually divided into high- and low-T2 endotypes according to the expression levels of molecules involved in T2 inflammation [9]. The most common endotype of childhood asthma is the high-T2 endotype. The airway epithelium, as the first line of interaction between the respiratory tract and the environment, contributes to the pathogenesis of asthma upon contact with various environmental factors, including allergens, infectious viruses, and air pollutants [4,5,6]. The immunologic functions of the airway epithelium are also thought to play crucial roles in asthma development [3]. Three epithelium-derived cytokines, TSLP, IL25, and IL33 have been recognized as key mediators of T2-immunity [13]. Genetic polymorphisms in these cytokines and their receptors are associated with the pathogenesis of asthma, and agents that block TSLP have shown clinical efficacy in certain patients [4]. However, the molecular and cellular mechanisms of allergens exposure triggering the release of epithelial IL33 remain unclear [12]. In this study, we explored the function of EGR1 on epithelial IL33 production and its role in high-T2 asthma development. The levels of epithelial EGR1 are elevated in HDM-induced asthma. The expressions of EGR1 were shown to be concentrated in epithelial cell clusters via scRNA-seq. The increased expressions of EGR1 were consistent with the hypersecretion of IL33. Functional experiments verified that EGR1 could directly regulate the production of IL33. These results highlight the critical role of EGR1 in IL33 production in asthma models.

This comprehensive analysis of multiple bulk RNA-seq and scRNA-seq datasets revealed that EGR1 is an epithelial cell-derived transcription factor associated with the T2 immune response, especially in the context of IL33 production. EGR1 (early growth response 1) is a Cys2-His2-type zinc finger transcription factor that can bind to specific DNA promoter sequences [26]. Previous studies [27,28,29] have reported that EGR1 can be activated by external stimuli such as oxidative stress, endotoxins, and cytokines and participates in a variety of biological processes, including cell proliferation, cell differentiation, apoptosis, and immune inflammation. EGR1 expression is upregulated in various acute and chronic respiratory inflammatory diseases, amplifying pathological signals from the extracellular environment and driving disease progression [27, 30, 31]. A previous study revealed that EGR1 is a gatekeeper of inflammatory enhancers in human macrophages [28]. In our study, we innovatively investigated the upregulated expression of EGR1 in a HDM-induced asthma model. Consistent with our results, a previous study demonstrated that polymorphisms in the EGR1 gene are associated with total IgE levels and atopy in children with asthma [32]. In adults with allergic rhinitis [33], higher levels of EGR1 were positively associated with atopy and increased total IgE production in the plasma. In addition, EGR1 was induced concomitant with increased MUC5AC levels and mucus production in human bronchial epithelial cells incubated with cigarette extracts [34].

Transcription factors have attracted increasing attention as critical regulators of various physiological and pathological processes. EGR1 regulates thousands of downstream target genes [27, 30]. As a key transcription factor, EGR1 is involved in signal transduction and the activation of inflammatory immune cells and plays a crucial role in promoting inflammatory respiratory disease progression [27, 28, 35]. As mentioned above, the upstream regulation of epithelial inflammatory IL33 productions in asthma still needs to be explored further [36]. As there have few studies investigated the role of EGR1 in epithelial IL33 expressions in asthma, we explored the underlying mechanisms in vitro. Using a combination of online tools and experimental assays, such as SCENIC analysis, ChIP and dual-luciferase reporter assays, we predicted and verified that EGR1 could directly regulate IL33 productions. These findings demonstrate that EGR1 could serve as a potential target for early IL33 intervention in allergic asthma.

The potential limitations of this study should be acknowledged. As lower airway specimens are inaccessible in most pediatric cohorts, we used adult airway specimens for scRNA-seq, which may not fully recapitulate the conditions in children. However, we predominantly used scRNA-seq for analyses of epithelial cell clusters. In the future, scRNA-seq of pediatric airway samples should be performed. Moreover, further exploration of the role of EGR1 in asthma and extension of this research to other signaling pathways or other cell types should be performed to more fully understand the mechanisms underlying the development of childhood asthma.

Taken together, the transcription factor EGR1 is upregulated in asthma and induces airway inflammation by promoting epithelial IL33 production. Our findings indicate that EGR1 may serve as a potential target for asthma management strategies, and additional research is encouraged to identify the exact mechanism by which EGR1 regulates T2 asthma.