Integrated Bioinformatics and Experimental Validation of Inflammatory Cytokine-Related Hub Genes in Pediatric Asthma

Introduction

Asthma is a heterogeneous chronic respiratory disease characterized by airway hyperresponsiveness, mucus hypersecretion, airway remodeling, and variable airflow limitation.^1,2 It remains a major global health challenge, with a steadily increasing prevalence and a substantial burden on patients and healthcare systems worldwide.^3–5 Although inhaled corticosteroids and bronchodilators are effective for many patients, a proportion of individuals still experience poor disease control or persistent symptoms despite treatment. Asthma is also increasingly recognized as a biologically heterogeneous disorder with distinct inflammatory endotypes, broadly including T2-high and T2-low patterns, which differ in immune mechanisms, biomarker profiles, and treatment responsiveness. Therefore, a better understanding of the molecular mechanisms and biomarker profiles underlying asthma is needed to improve disease characterization and management.

Inflammatory cells, particularly eosinophils, mast cells, neutrophils, and specific T lymphocyte subsets, are critical in asthma’s pathophysiology.^6,7 These cells release a variety of inflammatory mediators, such as TNF-α, interleukins IL-4, IL-5, IL-13, and chemokines CCL11 and CCL24, which contribute to airway inflammation, structural remodeling, and disease progression. Although the pathogenic relevance of these mediators has been extensively investigated, an integrated peripheral blood transcriptomic framework focused on inflammatory cytokine-related genes and their immune associations in pediatric asthma remains limited. Compared with airway tissue, peripheral blood provides a relatively accessible and non-invasive source for biomarker discovery and may reflect systemic immune-inflammatory alterations relevant to asthma, although it cannot fully represent the local airway microenvironment. In addition, accumulating evidence suggests that viral triggers, particularly respiratory viruses, can contribute to asthma exacerbation and immune activation, which may be relevant when interpreting infection-related pathway enrichment signals in transcriptomic analyses.⁸

In this study, we integrated two Gene Expression Omnibus (GEO) microarray datasets containing peripheral blood samples from children with asthma and healthy pediatric controls. Differential expression analysis, weighted gene co-expression network analysis (WGCNA), and inflammatory cytokine-related gene intersection were performed to identify candidate genes associated with pediatric asthma. Functional enrichment analyses, including GO, KEGG, and GSEA, were used to explore the underlying biological pathways. In addition, immune cell infiltration analysis using the CIBERSORT algorithm,⁹ transcriptional regulatory analysis, and predictive modeling were conducted to further characterize the immune and molecular features of these candidate genes. Selected genes were then subjected to preliminary validation in clinical samples and an experimental mouse model. Through this integrated strategy, we aimed to identify inflammatory cytokine-related candidate biomarkers in pediatric asthma. The novelty of the present study lies not in proposing definitive therapeutic targets or a completely new asthma mechanism, but in the integrated prioritization and preliminary validation of candidate immune-inflammatory biomarkers specifically in pediatric asthma. 

Materials and Methods

Data Collection

The RNA-seq data for pediatric asthma were obtained from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). The GSE106230 dataset contained nine pediatric asthma cases and three healthy pediatric control samples, the GSE117038 dataset included nine pediatric asthma cases and three healthy pediatric control samples, and the GSE195599 dataset contained four pediatric asthma cases and three healthy pediatric control samples. All samples were human peripheral blood samples. In addition, 2,001 inflammatory cytokine-related genes were retrieved from the GeneCards database (http://www.genecards.org/). All resources were accessed on 24/08/2024.

Identification of Differentially Expressed Genes (DEGs) in Pediatric Asthma

Principal component analysis (PCA) was performed on the combined dataset (GSE117038 and GSE106230) to visualize sample distribution and assess overall between-group separation. Differential expression analysis between pediatric asthma samples and healthy pediatric control samples was conducted using the “limma” R package (version 3.52.2), with the thresholds set at |logFC| > 0.5 and p < 0.05. The “ggplot2” R package (version 3.3.2) and the “pheatmap” R package (version 1.0.12) was used to visualize the expression of DEGs.¹⁰

Weighted Gene Co-Expression Network Analysis (WGCNA)

The co-expression network was constructed in the combined dataset using “WGCNA” R package (version 1.7–3).¹¹ Pediatric asthma and healthy pediatric control status were used as the clinical trait for module screening. First, genes with variance in the top 25% were selected for network construction. Samples were then clustered to detect potential outliers and ensure the reliability of subsequent analyses. A sample dendrogram and trait heatmap were generated, and the soft-thresholding power was determined according to the scale-free topology criterion. The adjacency matrix was then calculated and transformed into a topological overlap matrix to assess gene connectivity. Gene modules were identified using the dynamic tree-cutting algorithm, with the minimum module size set to 150. Finally, correlations between each module and the clinical trait were evaluated, and the module showing the strongest association with pediatric asthma was selected for subsequent analysis.

Function Analysis of Characteristic Genes

The characteristic genes were retrieved by intersecting the DEGs, model genes and inflammatory cytokines related genes. The GO function including biological process (BP), cellular components (CC) and molecular functions (MF) enrichment analysis, and KEGG pathway enrichment analysis were conducted by “clusterprofiler” R package (version 4.4.4), and the enrichment results were drawn by “enrichplot” R package (version 1.10.2).^12,13

Screened for the Key Genes of Asthma

The protein-protein interaction (PPI) network was constructed to investigate the relationship of characteristic genes by the “STRING” online website (https://string-db.org), and “Cytoscape” (version 3.8.2) (Confidence = 0.4) was used to visualize the results.^14,15 Finally, the genes of the top 5 maximum mass centrality (MCC) obtained by PPI network were defined as hub genes.

In order to study the ability of hub genes to distinguish the asthma, the “pROC” R package (version 1.18.0) was used to draw the receiver operating characteristic (ROC) curves of each target genes in combined dataset and GSE195599 datasets. Target genes with the area under the ROC curve (AUC) value ≥ 0.7 were defined as key genes.

Establishment and Assessment of Nomogram

A nomogram was developed in the training set utilizing the rms package (version 6.5.0) to evaluate the probability of biomarkers in predicting the onset of asthma. Additionally, a calibration curve was generated to assess the accuracy of the nomogram, utilizing the rms package (version 6.5.0). Simultaneously, the receiver operating characteristic (ROC) curve was plotted utilizing the pROC package (version 1.18.0), and the area under the curve (AUC) value was computed to evaluate the diagnostic effectiveness of the nomogram model for asthma. An AUC value exceeding 0.7 indicated favourable predictive performance for the nomogram model. Moreover, Furthermore, decision curve analysis (DCA) was performed utilizing the rmda package (version 1.6) (https://CRAN.R-project.org/package=rmda) to evaluate the clinical utility of the nomogram.

Construction of Regulatory Network

Molecular regulatory networks offered insights into the fundamental mechanisms underlying gene regulation in disease processes. The NetworkAnalyst online platform (https://www.networkanalyst.ca/) provided access to ENCODE database (https://www.encodeproject.org/) for the purpose of predicting the transcription factors (TFs) of targeted key genes. The TF-mRNA regulatory network was visualized utilizing the Cytoscape software (version 3.8.2).

Drug/Compound Prediction and Molecular Docking

The Comparative Toxicogenomics Database (CTD) (http://ctdbase.org/) identified drugs or compounds potentially targeting key genes, with those having an interaction count below 2 being excluded. The drug/compoundg-key gene network map was visualized utilizing the Cytoscape software (version 3.8.2). Drugs or compounds with the highest interaction count and 3D structures were selected as drug/compound candidates. The candidate drug/compound was uploaded to the public chemistry database (PubChem) database (https://pubchem.ncbi.nlm.nih.gov/) to retrieve their respective 3D structures. Simultaneously, the 3D structures of proteins corresponding to key genes were obtained from the Research Collaboratory for Universal Protein Resource (UniProt) (https://www.uniprot.org/). Molecular docking analysis was conducted utilizing the CB-Dock2 online platform (https://cadd.labshare.cn/cb-dock2/php/index.php).

GSEA of Key Genes

The correlation coefficient between each key gene and all genes in the combined dataset was calculated, and the GSEA of key genes was performed using the “clusterProfiler” R package (version 4.4.4)⁹ with “org.Hs.egdb” as the background gene set.

Analysis of Immune Micro-Environment

To evaluate the immune cell infiltration landscape, the single-sample Gene Set Enrichment Analysis (ssGSEA) algorithm was employed to calculate the relative abundances of 28 distinct immune cell types in both asthma and normal control samples within the combined dataset. The specific signature gene sets used to define these 28 immune cell subpopulations were obtained from Charoentong et al, and the detailed gene lists are provided in Table S1.¹⁶ Differences in the infiltration levels of these immune cells between the two groups were compared using the Wilcoxon rank-sum test. Furthermore, Spearman’s rank correlation analysis was conducted to investigate the associations between the expression of key genes and the differentially infiltrated immune cells.

Real-Time Quantitative Polymerase Chain Reaction (RT-qPCR)

Based on 10 pairs of asthma blood samples, RT-qPCR was performed to validate the expression of key genes. These samples were obtained with patient consent and approved by the Ethics Committee of The Second Hospital of Hebei Medical University. Total RNA was extracted using TRIzol (Thermo Fisher, Shanghai, CN), and mRNA was reverse-transcribed into cDNA for qPCR using the SureScript First-strand cDNA synthesis kit (Servicebio, Wuhan, CN). For the normalization strategy, GAPDH was utilized as the internal reference, and relative expression levels were determined via the 2^−ΔΔCT method. The primer sequences for the key genes are listed in Table 1.

Table 1 The Primer Sequences for the Key Genes

Animal Model of Asthma

All procedures were approved by the Ethics Committee/IACUC of the Second Hospital of Hebei Medical University (approval No. 2021-AE012). Male BALB/c mice (6–8 weeks) were randomly assigned to two groups (n = 6 mice per group) and housed under SPF conditions. Mice were sensitized i.p. with 20 μg OVA + 2 mg alum on days 0 and 14, then challenged with 1% OVA aerosol for 30 min/day on days 21–27; controls received PBS. Pulmonary function testing was performed under isoflurane inhalation anesthesia (2–3% in oxygen via a precision vaporizer). All procedures involving anesthesia or euthanasia complied with the American Veterinary Medical Association (AVMA) Guidelines for the Euthanasia of Animals (2020 Edition). To minimize distress at the experimental endpoint (day 28), mice were humanely euthanized by isoflurane overdose (≥5% in oxygen in an induction chamber) followed by cervical dislocation to ensure death, as recommended by AVMA guidelines. The collected lung tissues (n = 6 per group) were subsequently processed for downstream RT-qPCR analysis and pathological staining. Animals were monitored at least daily with predefined humane endpoints observed.

Histological Examination

After harvesting, the lung tissues were fixed in 4% paraformaldehyde for 24 hours, embedded in paraffin, and sectioned at 4 μm thickness for histological analysis:

Hematoxylin and Eosin (HE) Staining: Sections were deparaffinized, rehydrated, and stained with hematoxylin for 5 minutes followed by eosin for 3 minutes. HE staining was used to assess general morphology and inflammatory cell infiltration in lung tissue.

Periodic Acid-Schiff (PAS) Staining: PAS staining was used to evaluate mucus production in the airways. Sections were oxidized in periodic acid solution for 5 minutes, followed by Schiff reagent for 15 minutes. The intensity of mucus staining was evaluated to assess goblet cell hyperplasia.

Statistical Analysis

All data are presented as mean ± standard deviation (SD). For immune infiltration analysis, differences between the two groups were compared using the Wilcoxon rank-sum test, and correlations between candidate hub genes and differentially infiltrated immune cell types were assessed using Spearman correlation analysis. For experimental validation, statistical significance between two groups was evaluated using Student’s t-test, whereas one-way ANOVA was used for multiple-group comparisons. For enrichment analyses and GSEA, multiple testing correction was performed using the Benjamini–Hochberg method where applicable. Statistical analyses were performed using R software (version 4.1.0), and P < 0.05 was considered statistically significant.

Results

Identification of Characteristic Genes in Asthma

The PCA results showed that the combined dataset without obvious batch effect could be used for further analysis Figure 1A. There were 213 DEGs (56 up-regulated and 157 down-regulated) between 18 asthma and 6 normal samples Figure 1B. The WGCNA was used to identify the module genes significantly associated with asthma. Sample clustering analysis was implemented, and the results showed that there were no outlier samples Figure 1C. The soft threshold value = 4 was identified as the optimal soft threshold. When the ordinate scale-free fit index, signed R² approached the threshold value of 0.85 (red line), the network was close to scale-free distribution and mean connectivity was close to 0 (Figure 1D). A total of 14 modules were obtained Figure 1E. Then the correlation between modules and traits was evaluated, and the results showed that the yellow module had a significant correlation with asthma (|cor| = 0.5 and p = 0.01) (Figure 1F). Then, 1,058 genes in the key module were screened for subsequent analysis. Finally, 12 characteristic genes were obtained by overlapping the 213 DEGs, 1,058 module genes and 2,001 inflammatory cytokines related genes (Figure 1G). Besides, these 12 characteristic genes have high correlations with each other (Figure 1H, Table S2).

Figure 1 Identification of characteristic genes. (A) Disease and normal samples expression boxplots and PCA box plot. (B) Heat map for DEGs. Higher the expression the redder the colour, the lower the expression the bluer the colour. (C) Merged data sample clustering and phenotypic information. (D) Soft threshold screening. (E) Cluster diagram of module genes. (F) Heat map of correlation between modules and clinical traits. (G) Venn diagram of characteristic genes. (H) Characteristic genes Spearman correlation analysis.

Enrichment Analysis of 12 Characteristic Genes

GO function enrichment and KEGG pathway analyses were performed to assess the function of 12 characteristic genes. These genes were enriched to 346 GO BP, 22 GO CC, 24 GO MF, and 7 KEGG pathway, including respiratory burst, positive regulation of cell adhesion, leukocyte migration, osteoclast differentiation, NOD-like receptor signaling pathway, NOD-like receptor signaling pathway and etc. (Figure 2A and B, Tables S3 and S4).

Figure 2 The enrichment analysis of characterisatic genes. (A) GO enrichment histogram (TOP10). Green for BP, Orange for CC, blue for MF. (B) KEGG enriched bubble diagram.

CD19, RELB, IL12RB1 and RETN Genes Identified as Key Genes

The PPI network was constructed with 11 characteristic genes and 15 reciprocal relationships (Confidence = 0.15), and the top 5 genes of MCC (CD19, RELB, IL12RB1, TNFRSF14, RETN) were defined as hub genes Figure 3A and B. In addition, the AUC value of CD19, RELB, IL12RB1, RETN were greater than 0.7 in both combined dataset and GSE195599 datasets (Figure 3C and D). Therefore, CD19, RELB, IL12RB1, and RETN were retained as candidate hub genes for subsequent analyses. (Table S5).

Figure 3 Identification of key genes. (A) Networks diagram of interactions of characterised genes. The greater the centrality of the largest clusters, the darker the color. (B) PPI network node maximum clique centrality ranking. (C and D) The ROC curve of GSE117038 (C) and GSE195599 (D) datasets.

Nomogram Exhibited Promising Performance in Evaluating the Diagnosis of Asthma

The construction of a nomogram was undertaken for the purpose of further evaluation of the diagnostic value of CD19, RELB, IL12RB1 and RETN for asthma. The nomogram model indicated an association between increased total points from CD19, RELB, IL12RB1 and RETN and a higher likelihood of asthma occurrence (Figure 4A). Furthermore, as demonstrated by the calibration curves in Figure 4B, the probabilities predicted by the nomogram closely adhered to the reference line, showed acceptable predictive performance. The DCA results indicated that the nomogram model offered a greater net benefit Figure 4C, suggesting potential diagnostic utility. Besides, the AUC value of this nomogram model was 0.852 (Figure 4D), suggesting that the model’s predictive accuracy for asthma had attained an acceptable level.

Figure 4 Development and validation of a diagnostic nomogram for asthma risk based on CD19, RELB, IL12RB1, and RETN. (A) Nomogram construction: Higher total points from combined CD19, RELB, IL12RB1, and RETN expression correlate with increased asthma likelihood. (B) Calibration curve: Predicted probabilities closely align with the 45-degree reference line.(C) DCA: The nomogram demonstrates superior net benefit across threshold probabilities.(D) ROC curve: The model achieves an AUC of 0.852.

Multiple Transcriptional Regulators Coordinate the Expression of CD19, RELB, IL12RB1 and RETN

The regulatory factors of the CD19, RELB, IL12RB1 and RETN were subjected to further investigation. Specifically, CD19 was regulated by SP1 and NR2F1; RELB by SP1 and SMARCA5; IL12RB1 by NR2F1, SMARCA5, and REST; and RETN by REST (Figure 5). The findings indicated that these pivotal asthma-related genes are governed by intricate transcriptional regulatory mechanisms and likely participate in the disease’s multifaceted regulatory process.

Figure 5 CD19 was regulated by SP1 and NR2F1; RELB by SP1 and SMARCA5; IL12RB1 by NR2F1, SMARCA5, and REST; and RETN by REST.

Compound Screening and Molecular Docking Analysis of Asthma-Associated Candidate Biomarkers

Compound prediction analysis identified 23, 61, 16, and 55 candidate compounds for CD19, RELB, IL12RB1, and RETN, respectively (Figure 6A). Among the top-ranked compounds with available 3D structures, bisphenol A, tetrachlorodibenzodioxin, benzo(a)pyrene, and tretinoin were selected for CD19, RELB, IL12RB1, and RETN, respectively (Table S6). The predicted binding energies for the docking of CD19 with bisphenol A, RELB with tetrachlorodibenzodioxin, IL12RB1 with benzo(a)pyrene, and RETN with tretinoin were −6.4, −6.2, −6.0, and −5.7 kcal/mol, respectively, and the corresponding docking poses are shown in Figure 6B–E. Using a docking score threshold of ≤ −5 kcal/mol, these results suggested potential interactions between the selected compounds and the candidate genes.

Figure 6 Drug prediction and molecular docking analysis of asthma-associated biomarkers. (A) Predicted potential drug counts targeting CD19 (23 drugs), RELB (61), IL12RB1 (16), and RETN (55).(B–E) Molecular docking results and binding free energies for top ligand-receptor pairs: (B) CD19-bisphenol A complex (−6.4 kcal/mol), (C) RELB-tetrachlorodibenzodioxin complex (−6.2 kcal/mol), (D) IL12RB1-tretinoin complex (−6.0 kcal/mol). (E) RETN-tretinoin complex (−5.7 kcal/mol). All binding energies exceed the robustness threshold (≤ −5.0 kcal/mol), confirming strong ligand-receptor affinity.

GSEA Revealed Enrichment of HSV-1 Infection–Related Gene Sets in the Candidate Genes

The GSEA results for the four candidate genes are shown in Figure 7A–D. GSEA suggested that herpes simplex virus 1 infection–related pathways were enriched for each of the four genes, and the corresponding leading-edge genes are summarized in Table S7. In addition, these leading-edge genes were mainly involved in immune-related, ribonucleoprotein complex biogenesis, ribosome biogenesis, rRNA processing, and spliceosome-associated functions. Gene sets associated with ribonucleoprotein complex biogenesis, ribosome biogenesis, rRNA processing, and spliceosome were mainly enriched in CD19, RELB, and IL12RB1.

Figure 7 The GSEA enrichment analysis of key genes. (A–D) GSEA enrichment analysis of for (A) CD19, (B) RELB, (C) IL12RB1, (D) RETN. The left side of each subgraph denotes KEGG and the right side denotes GO.

Key Genes May Influence Asthma Through Immune Cells

MDSC (84.86%, 83.43%), Neutrophil (79.6%, 79.38%), Monocyte (670.63%, 9.44%) were the top three immune cells in both asthma and normal groups Figure 8A and B. In addition, the Activated B cell, Central memory CD8 T cell and T follicular helper cell were significantly different between asthma and normal samples Figure 8C. The correlation analysis results showed that there was a significant strongest positive correlation between CD19 and Activated B cell, and a significant strongest negative correlation between CD19 and Type 1 T helper cell (p < 0.01). Besides, there was also a significant positive correlation between IL12RB1 and Central memory CD8 T cell Figure 8D and E.

Figure 8 Immune infiltration analysis of key genes. (A) Immune infiltration score heat map. (B) Rose chart of immune cell proportion. (C) Box plots of immune infiltration scores. (D) Heat map of key gene and immune cell correlation(E) Scatter map of correlation between key genes and immune cells. *, P < 0.05; **, P < 0.01.

Expression Validation of Key Genes

RT-qPCR showed that CD19, RELB, and IL12RB1 were significantly upregulated in asthma samples compared with normal samples (p < 0.01), whereas RETN did not show a statistically significant difference (Figure 9A–D).

Figure 9 Relative mRNA expression levels of key genes in the Normal and Case groups. The relative mRNA expression levels of IL12RB1 (A), RETN (B), CD19 (C), and RELB (D) were measured and normalized to GAPDH. Data are presented as mean ± SD. Statistical significance was determined using an unpaired Student’s t-test. ns, not significant; *, P < 0.05; ***, P < 0.001; ****, P < 0.0001 compared with the Normal group.

Validation in Animal Model

To further validate the differential gene expression observed in asthma, we conducted experiments in a mouse model. An OVA-induced asthma model was established in mice. Histological analysis of lung tissues using hematoxylin and eosin (H&E) staining confirmed significant airway inflammation, characterized by the infiltration of inflammatory cells, thickening of the airway epithelium, and hyperplasia of goblet cells. Additionally, periodic acid-Schiff (PAS) staining demonstrated an increased presence of mucus-producing goblet cells, further validating the successful induction of the asthma phenotype in the model Figure 10A.

Figure 10 Lung histopathological changes and related gene expression levels in the OVA-induced airway inflammation model. (A) Representative images of Hematoxylin-Eosin (HE) and Periodic Acid-Schiff (PAS) staining of lung tissues from the Control and OVA groups. The two left columns show low-magnification views (10×), and the two right columns show high-magnification views (40×). (B–D) The relative mRNA expression levels of CD19 (B), RELB (C), and IL12RB1 (D) in lung tissues of the Control and OVA groups were measured by RT-qPCR. Data are presented as mean ± SD [or SEM]. * P < 0.05 compared with the Control group.

Subsequently, lung tissues were harvested from both OVA-induced asthmatic mice and control mice for RT-qPCR analysis to measure the expression levels of CD19, RELB, and IL12RB1. Consistent with the human asthma samples, the expression of these genes was significantly elevated in the lung tissues of asthmatic mice compared to controls (p < 0.05). These findings further support the role of these genes in asthma pathogenesis Figure 10B–D.

Discussion

Asthma is characterized by persistent chronic airway inflammation, involving the activation and interplay of various inflammatory cells and cytokines. Inflammatory cytokines, which are endogenous peptides with potent biological effects, are primarily produced by immune system cells. Upon release into the extracellular environment, these molecules initiate inflammatory cascades that contribute to a spectrum of pathological conditions, including asthma. A murine model study demonstrated that pro-inflammatory cytokines can target the murine respiratory system, ultimately inducing asthma-like phenotypes.¹⁷ Furthermore, research has shown that upregulated NF-κB, IL-4, IL-9, and FOXP3 play pivotal roles in sustaining and exacerbating the inflammatory milieu associated with asthma.^18,19 Consequently, comprehensive analysis of inflammatory gene expression may facilitate biomarker discovery and improve understanding of immune-inflammatory dysregulation underlying pediatric asthma.^20,21

In the present work, through comprehensive analysis of gene expression data from pediatric asthma patients and healthy pediatric control samples, we identified 12 DEGs closely associated with pediatric asthma. Further gene enrichment analysis highlighted a central role for cell adhesion and leukocyte migration in the pathogenesis of asthma. Research indicates that cell adhesion factors, including ICAM-1, TGFβ, and L-selectin, play a pivotal role in the migration of inflammatory cells and their interaction with airway tissues. They significantly promote the migration of inflammatory cells like monocytes, eosinophils, and T cells towards the airway tissues, thereby exacerbating pulmonary inflammation and airway hyper-responsiveness.²² Notably, it has been demonstrated that treatment approaches that target adhesion molecules, like inhaling integrin antagonists, can effectively limit the recruitment of inflammatory cells in lung tissue.²³ This approach highlights the potential relevance of adhesion-related pathways to asthma-associated inflammatory responses and supports their further investigation in asthma biology. These findings support the potential involvement of cell adhesion and leukocyte migration in asthma-related immune dysregulation. Because pediatric asthma differs from adult asthma in immune maturation, allergen sensitization patterns, and developmental context, integrated biomarker investigation in childhood asthma may provide complementary information that cannot be directly extrapolated from adult studies.

CD19, RELB, IL12RB1, TNFRSF14, and RETN were identified as hub genes by PPI analysis, and four of them (CD19, RELB, IL12RB1, and RETN) were subsequently retained as candidate hub genes for further analyses. Importantly, these genes should be interpreted as candidate biomarkers rather than definitive mechanistic drivers or therapeutic targets at the current stage. Of which, CD19, as a B-cell surface molecules, is involved in B-cell activation and regulation. Research by Brosseau et al revealed that CD19 can regulate the secretion of IL10 by CD9+ B cells, thereby inducing the apoptosis of T cells, potentially playing a significant role in the pathophysiology of asthma.²⁴ Further, Chai et al’s study showed that inhalation of high concentrations of SO₂ significantly reduces CD19 expression, leading to nasal epithelial injury and increased airway hyperreactivity in rats, thus emphasizing the importance of CD19 in asthma pathogenesis RELB, a member of the NF-κB family, has been implicated in the function and maturation of dendritic cells.²⁵ Nair et al developed a RELB-deficient mouse model and discovered that the absence of RELB exacerbates airway inflammatory responses, elevates levels of chemokines, IL4, and IL5 in lung tissue, leading to collagen deposition and epithelial thickening, causing spontaneous allergic airway inflammation.²⁶ Pareek S et al’s study also demonstrated that downregulation of RELB is associated with an increase in neutrophils in lung tissue, further supporting its possible involvement in asthma-related immune regulation.²⁷ IL12RB1, a component of the IL12 and IL23 receptor complex, is commonly considered an inducer of delayed-type hypersensitivity reactions. Recent research through genome-wide association analysis in European populations found that IL12RB1 expression is related to changes in lung function in asthma patients.²⁸ Additionally, Dara et al’s research found that gene mutations in IL12RB1 are closely associated with the onset of asthma in Europeans.²⁹ RETN, a molecule involved in inflammatory responses and immune regulation, may be relevant to airway inflammatory responses and the pathological process of asthma. M T Kasaia et al, through proteomic analysis of serum and sputum secretions from asthma patients, found low levels of RETN in the sputum of asthma patients, supporting its possible involvement in the onset of asthma.³⁰ Collectively, these findings suggest that these four genes may have potential biomarker value in pediatric asthma; however, the current evidence is insufficient to define them as therapeutic targets, and stronger in vivo, in vitro, and clinical evidence is required before any translational significance can be established. From a broader biological perspective, the prioritization of CD19, RELB, IL12RB1, and RETN suggests that pediatric asthma may involve coordinated alterations in B-cell regulation, NF-κB-related immune signaling, cytokine receptor-mediated responses, and systemic inflammatory regulation. Although the present study does not establish direct mechanisms, these findings are consistent with the concept that pediatric asthma is shaped by complex immune-inflammatory dysregulation rather than isolated single-gene effects. Given the differences between pediatric and adult asthma in immune maturation, sensitization patterns, and developmental context, these candidate biomarkers may provide pediatric-specific clues that warrant further validation rather than direct extrapolation from adult studies.

  Emerging evidence underscores the contribution of innate immune cells to diverse asthma phenotypes. B cells, as central mediators of allergic reactions and tolerance, exert both pathogenic and regulatory effects.³¹ Activated B cells can produce immunoglobulin E (IgE), when bound to allergens, triggers an inflammatory cell response including mast cells and eosinophils, leading to airway inflammation and asthma symptoms.³² However, recent research has also revealed that B cells can exhibit anti-inflammatory properties by secreting cytokines such as IL10, which reduce the activation of T cells and suppress the production of pro-inflammatory cytokines. Moreover, Carole et al’s study found that regulatory B cells expressing CD19 can induce the apoptosis of T cells by secreting IL10, thereby effectively improving airway inflammatory responses.²⁴ Additionally, CD19+ regulatory B cells interact with Type 1 T helper (Th1) cells, inhibiting Th1 cell differentiation and function through the secretion of anti-inflammatory cytokines like IL-10, thereby ameliorating Th1-mediated inflammatory responses and tissue damage.³³ Our data reveal reduced peripheral-blood B-cell levels in asthma, a positive correlation between CD19 expression and B-cell abundance, and an inverse correlation with Th1 cells, indicating a potential CD19+ B-cell–Th1 axis in disease pathogenesis. Our study highlights the potential role of IL12RB1 in central-memory CD8+ T cells, during asthma development. Given that these cells provide rapid and durable immune responses.³⁴ IL12RB1 may modulate asthma via memory-T-cell regulation. However, limited patient-level metadata precluded mechanistic dissection in the present analysis.  Importantly, the present study should be interpreted as an exploratory biomarker-discovery analysis rather than a definitive mechanistic or therapeutic target identification study. Although the highlighted genes were supported by integrated transcriptomic screening and preliminary validation, the current evidence remains insufficient to establish direct causal roles or justify these molecules as clinically actionable therapeutic targets in pediatric asthma. In the context of current precision medicine for severe asthma, including biologics targeting IgE, IL-4/IL-13, IL-5, and TSLP, our findings should be viewed as a complementary transcriptomic resource for candidate prioritization rather than an immediate translational framework.

Notably, inconsistency was observed for RETN across analyses. Although RETN was prioritized as a candidate hub gene based on the integrated bioinformatic screening strategy, its differential expression was not statistically significant in the human RT-qPCR validation cohort. In addition, subsequent animal validation was focused on CD19, RELB, and IL12RB1. This discrepancy may be attributable to the limited sample size of the validation cohort, biological heterogeneity between public transcriptomic datasets and clinical samples, and differences in tissue source or disease status. Similarly, although RELB was validated at the expression level, the downstream evidence supporting RELB appeared less extensive than that for CD19 and IL12RB1. Therefore, RELB and especially RETN should be interpreted as candidate biomarkers requiring further validation rather than definitive markers.

Although this study identified four candidate biomarkers potentially related to asthma, the exact biological mechanisms underlying these associations remain to be further elucidated. Several limitations should be acknowledged. First, the integrated analysis of public datasets may have been influenced by platform-related variation, batch effects, and differences in sample processing across cohorts, which could affect the stability and generalizability of the identified biomarkers. Second, the transcriptomic analyses were based primarily on peripheral blood samples and may not fully reflect the local airway microenvironment and tissue-specific inflammatory processes in asthma. Third, the validation cohort was relatively small, which limited statistical power and reduced the strength of external validation, particularly for RETN. Therefore, the current findings should be interpreted as preliminary and require confirmation in larger pediatric cohorts. Fourth, protein-level validation was lacking, limiting biological and clinical interpretability. In addition, the molecular docking results provide only in silico support for potential interactions and cannot substitute for direct biochemical or functional validation. Finally, although the nomogram showed favorable predictive performance in the available datasets, its clinical applicability remains to be confirmed in larger, independent, and prospective cohorts. Overall, the main contribution of this study lies in the integrated identification and preliminary validation of candidate inflammatory cytokine-related biomarkers in pediatric asthma rather than in establishing definitive mechanisms or therapeutic targets.

Conclusions

In this study, we used multiple bioinformatics tools to identify and validate four candidate hub genes that may be associated with the development of asthma. Our findings emphasize the critical importance of investigating immune cell subpopulations comprehensively to understand the complex biology of asthma. Moreover, our findings provide substantial support for the formulation of more efficacious treatment strategies aimed at addressing the complexities of asthma pathology.