This study was approved by the Ethics Committee of the Children’s Hospital of Chongqing Medical University. Written informed consent was obtained from the legal guardians of the study participants before enrollment.

A total of 32 preschoolers with RW (defined as ≥ 3 physician-diagnosed wheezing episodes)  (21, 22), 15 children with SA  (23), and 18 control individuals (patients without current respiratory tract infection, current wheezing, or history of allergy or wheezing) undergoing bronchoscopy were enrolled in the study. Collection of bronchoalveolar lavage fluid (BALF) from controls and patients with RW was carried out using standard procedures  (24). BALF was gently aspirated and centrifuged at 2 500 rpm for 5 min at 4°C after collection. The bronchoalveolar lavage (BAL) cells were resuspended in phosphate-buffered saline (PBS) and stored at −80°C. Details on subject characteristics are included in Supplementary File 1 ; Table S1 .

RNA-seq data (raw count matrix) of RW, SA, and healthy controls were obtained from the Gene Expression Omnibus (GEO) database at the National Center for Biotechnology Information (NCBI; accession number GSE118761). This dataset included nasal and tracheal brushing samples from RW (n = 14), SA (n = 13), and healthy (n = 14) groups ( Supplementary File 1 ; Table S2 ). To prevent sampling noise caused by lowly expressed genes, genes with low expression (< 30 counts in total across all nasal or tracheal samples) were excluded. For subsequent analysis and visualization of count data, variance stabilizing transformation (VST) normalized expression values were calculated using the “DESeq2”  (25) R package (v1.34.0).

The DEGs between groups (RW or SA vs. control) were identified using the “DESeq2” R package (v1.34.0), with P < 0.05 and |log2fold-change| > 0.5 considered significant. To identify common DEGs (co-DEGs) between RW and SA in the upper and lower airways, the “UpSetR” (26) R package (v1.4.0) was used to construct an UpSet diagram.

The scRNA-seq data (count matrix) of airway samples collected from 18 healthy children were obtained from the FigShare repository ( Supplementary File 1 ; Table S2 ). The R package Seurat  (27) (v4.1.2) was used to process the scRNA-seq data. To ensure high-quality single cells were used for downstream analysis, cells expressing fewer than 500 genes and cells with more than 15% mitochondrial reads were filtered out. In total, 38 399 filtered cells were used for further analysis.

Raw data were normalized using the “NormalizeData” function, and 2 500 highly variable genes were identified using “FindVariableFeatures”. Principal component analysis (PCA) was then performed for dimensionality reduction after data scaling. The top 40 principal components were selected for downstream analysis. The Uniform Manifold Approximation and Projection (UMAP) algorithm was used for cell visualization. The “FindClusters” function was used for cell clustering. To annotate cell clusters, the DEGs for each cell cluster were identified by comparing each cluster to all other clusters with the “FindAllMarkers” function. Genes with adjusted P < 0.05 were considered DEGs. The cell subsets were annotated based on the DEGs and known markers  (28).

The CIBERSORTx online platform (https://cibersortx.stanford.edu/) was applied to infer the cellular composition of the bulk airway transcriptomes. We first prepared and uploaded the single-cell expression matrix from our scRNA-seq analysis according to the CIBERSORTx instructions using default parameters. We ran “CIBERSORTx” and obtained a signature matrix of 18 cell types, then uploaded the gene expression matrix data for the RW, SA, and control groups. The previously obtained signature matrix was used for deconvolution analysis, with all parameters set to default. After running “CIBERSORTx”, we obtained the relative proportions of the 18 cell subsets in each sample.

WGCNA was performed using the “WGCNA”  (15) R package (v4.1.2) and the VST-normalized gene expression profile data were used for network construction. We first constructed nasal and tracheal co-expression gene networks for RW and SA in parallel using default parameters. The RW modules were then correlated with the SA modules ( Supplementary File 1 ; Figures S1A, B ). We calculated the overlap of each pair of RW-SA modules and used Fisher’s exact test to assign a P-value to each pairwise overlap. Next, we built nasal and tracheal consensus networks for RW and SA, respectively. Based on scale-free topology criteria, we selected appropriate soft-threshold power values for network construction ( Supplementary File 1 ; Figures S2 , S3 ). The Dynamic Tree Cut method was used to identify different modules, with modules showing similar expression patterns then merged ( Supplementary File 1 ; Figures S4A, B ). The minimum module size was set to 40 and the DeepSplit parameter was set to 2.5.

In order to identify consensus modules significantly correlated with clinical features of RW and SA (and in the same positive or negative direction), we performed Pearson correlation analysis of consensus module eigengenes and clinical traits (atopy and cellular compositions inferred from CIBERSORTx). Modules showing significant positive correlations with both RW and SA, as well as important phenotypic traits, were considered as key consensus modules.

To identify the biological functions of the co-DEGs and genes within the key consensus modules, Gene Ontology (GO) (29) and Kyoto Encyclopedia of Genes and Genomes (KEGG)  (30) pathway analyses were performed using the “clusterProfiler” R package (v4.2.2), with significant enrichment considered at P < 0.05.

Two machine learning algorithms, Random Forest (RF)  (31) and support vector machine-recursive feature elimination (SVM-RFE)  (32), were used to screen hub genes (i.e., genes within key consensus modules) highly correlated with RW and SA. Specifically, RF and SVM-RFE were used to estimate how well each candidate hub gene correctly classified RW or SA compared to the controls, and to detect the number of feature genes required to separate groups with maximum accuracy. The RF and SVM-RFE modeling procedures were based on three-fold cross-validation (CV) using the “randomForest” and “e1071” R packages, respectively. The classification accuracies of different numbers of feature genes were determined for RF and SVM-RFE, with those showing the highest classification accuracy retained to determine the final hub genes. Hub genes were selected by intersecting the co-DEGs and common feature genes identified by RF and SVM-RFE.

For testing hub gene efficacy, receiver operating characteristic (ROC) curves and corresponding areas under the ROC curves (AUC) were calculated for each hub gene based on their standardized expression levels.

To validate the expression profiles of the hub genes in RW and SA, external validation was performed with two independent SA nasal gene expression datasets (GSE19187 and GSE65204) and one RW nasal gene expression dataset (GSE103166). The GSE103166 and GSE65204 datasets were used to validate shared mechanisms between RW and SA. The GSE103166 dataset contained 56 RW and 21 healthy individuals, the GSE65204 dataset contained 36 SA and 33 healthy individuals, and the GSE19187 dataset contained 13 SA and 11 healthy individuals. Detailed information on these datasets is provided in Supplementary File 1 ; Table S2 .

Total RNA was extracted from BAL cells using TRIzol reagent (Invitrogen, USA), and purified using a Micro Total RNA Extraction Kit (Tianmo Biotech, China). cDNA was synthesized using a PrimeScript RT Kit (TaKaRa, Japan) according to the manufacturer’s instructions. Reactions were carried out in a total volume of 10 μL, including 5 μL of TB Green^®Premix Ex Taq™ II (TaKaRa, Japan), 0.2 μL of each specific primer, 2.6 μL of ddH₂O, and 2 μL of cDNA. The relative expression levels of the hub genes were calculated using the2^-ΔΔCt method. GAPDH was used as an internal reference. The specific primers for each gene are provided in Supplementary File 1 ; Table S3 .

To identify the potential regulatory transcription factors (TFs) for the hub genes, we performed enrichment analysis of TF binding motifs (TFBMs) and TFs surrounding the transcription start site (TSS) of genes using the “RcisTarget”  (33) R package (v1.14.0). Significantly enriched TFBMs (normalized enrichment score (NES) > 3.0) were annotated to TFs using the provided annotation database. The TF-target network was visualized with Cytoscape (v3.8.2)  (34).

All continuous data are expressed as mean ± standard deviation (SD). All statistical analyses were conducted using R (v4.1.2; https://www.r-project.org/). Wilcoxon’s rank-sum test was used to compare gene expression levels ( Figures 1F , 2C, D and Supplementary File 1 ; Figures S9E–G ) and cellular compositions inferred by CIBERSORTx between groups ( Figure 3C ). Pearson correlation analysis ( Figures 2E, F ) was performed to obtain correlation coefficients (r) and P-values. A P-value of < 0.05 was considered statistically significant.

Compared with the healthy controls, we identified 567 DEGs (335 up-regulated and 232 down-regulated) and 549 DEGs (201 up-regulated and 348 down-regulated) in the RW nasal and tracheal samples, respectively, and 838 DEGs (152 up-regulated and 686 down-regulated) and 412 DEGs (189 up-regulated and 223 down-regulated) in the SA nasal and tracheal samples, respectively ( Figure 4A ). The UpSet diagram in Figure 4B shows the number of overlapping DEGs between RW and SA across nasal and tracheal samples. In total, 16 co-DEGs were shared between RW and SA across the upper and lower airways ( Figures 4A, C ). Compared to the healthy controls, these genes were overrepresented in RW and SA. The co-DEGs were retained for identifying hub genes shared by RW and SA.

To investigate the biological behaviors of the co-DEGs, GO functional annotation and KEGG enrichment analysis were performed. Results indicated that the co-DEGs were mainly involved in sensory perception of taste, wound healing, endopeptidase inhibitor activity, and regulation of cell-matrix adhesion ( Figure 4D ). Salivary secretion and pancreatic secretion were the major biological pathways involved ( Figure 4E ).

To understand changes in the RW and SA airway microenvironments as well as potential critical cell populations shared by RW and SA, we performed cell deconvolution of the upper and lower bulk transcriptome data to deduce changes in cell-type frequency within the airways of RW and SA.

First, we analyzed the scRNA-seq data using cells from airway samples of 18 healthy children. After stringent quality control, 38 399 cells were retained for further analysis. Unsupervised clustering by UMAP identified 18 cell subsets ( Figure 3A ) based on representative marker genes from previous research   (28). The cell clusters identified by cell lineage-specific marker gene expression are shown in Supplementary File 1 ; Figure S5 . Based on the CIBERSORTx algorithm, a signature matrix of 2 360 genes in 18 cell clusters was created ( Supplementary File 2 ; Table S4 ).

Next, using the signature matrix derived from scRNA-seq analysis, we calculated the relative proportions of the 18 cell subsets for the RW, SA, and control samples using CIBERSORTx. The histogram in Figure 3B shows the proportions of each cell type in the upper and lower airways of preschool and school-age children. We then examined fractional differences in each cell type between RW and SA and the control group. Interestingly, a higher fraction of mast cells and lower fraction of club cells were observed in both the upper and lower airways of the RW and SA patients compared to the healthy controls ( Figure 3C ).

To explore the possibility that RW and SA share certain mechanisms and biological pathways, we analyzed similarities in their gene expression networks across nasal and tracheal transcriptomes using WGCNA. WGCNA can identify gene modules with similar gene expression patterns, thus providing biological insight based on the principal that genes with highly correlated expression patterns are likely to participate in the same biological processes. By constructing co-expression gene networks for RW and SA in parallel, we identified 25 and 25 network modules for RW and SA from their nasal transcriptomes, respectively ( Supplementary File 1 ; Figure S1A ), and 26 and 25 network modules for RW and SA from their tracheal transcriptomes, respectively ( Supplementary File 1 ; Figure S1B ). To determine whether network modules identified in RW could also be identified in SA, we explored the correlations between RW-specific and SA-specific modules. Results showed that most RW-specific modules had one or more SA module counterparts in both the upper and lower airways ( Supplementary File 1 ; Figures S1A, B ), thus suggesting possible similarities in the underlying transcriptional networks between RW and SA.

Next, we built consensus networks using the nasal and tracheal transcriptome data, with the resulting network modules representing robust gene co-expression patterns shared by RW and SA. The nasal and tracheal consensus networks contained 23 and 25 modules of co-expressed genes, respectively ( Supplementary File 1 ; Figures S4A, B ). To determine key consensus network modules highly correlated with both RW and SA, we performed Pearson correlation analysis between consensus module eigengenes and clinical traits of RW and SA, including the fractions of the 18 cell subsets inferred by CIBERSORTx. Heatmaps were constructed to show the relationship between nasal ( Supplementary File 1 ; Figures S6A, B ) and tracheal ( Supplementary File 1 ; Figures S7A, B ) consensus modules and clinical trait information for RW and SA. To obtain consensus module-trait relationships across both RW and SA, we integrated the separate consensus module-trait relationships specific to RW and SA. We retained the lower absolute value in the two sets with the same correlation sign, and “NA” for those with the opposite trend. Figures 5A, B shows the integrated consensus module-trait relationships in both RW and SA. Of note, in the upper and lower airway consensus networks, only two consensus modules (orange module in the nasal consensus network and plum module in the tracheal consensus network) were significantly positively correlated with RW and SA. Interestingly, these two consensus modules were also significantly positively correlated with atopy and mast cells. Therefore, the nasal orange module and tracheal plum module were selected as key consensus gene modules shared by RW and SA.

Next, we explored whether key consensus module gene signatures could be distinguished between children with RW and SA and healthy subjects. We performed hierarchical clustering analyses using the 66 ( Supplementary File 3 ; Table S5 ) and 98 ( Supplementary File 4 ; Table S6 ) genes in the nasal orange module and tracheal plum module, respectively. Heatmap visualization of the genes in the nasal orange module ( Figure 6A ) and tracheal plum module ( Figure 6B ) indicated that: (1) genes were more highly expressed in the RW and SA individuals; (2) hierarchical clustering grouped samples into two major clusters, with RW and SA clustered together based on similar expression patterns; and (3) hierarchical clustering using key consensus module gene signatures clearly distinguished RW and SA patients from the controls. Furthermore, PCA of the genes in the nasal orange module and tracheal plum module also distinguished RW and SA patients from the controls in both the nasal and tracheal samples ( Figures 6C, D ). Notably, the RW and SA individuals tended to cluster together. These results suggest the potential biological significance of the co-expressed genes in the key consensus modules and highlight the possibility of shared underlying molecular mechanisms between RW and SA.

We further explored shared mechanisms between RW and SA by identifying significantly overrepresented biological pathways and GO terms for the two key consensus modules. Overall, the nasal and tracheal consensus modules were associated with several similar GO terms, including fatty acid oxidation, fatty acid metabolic process, and sensory perception of bitter taste ( Figures 6E, F ). Fatty acid metabolism-related pathways were also significantly enriched in both the nasal orange and tracheal plum modules ( Figures 6E, F ).

To validate the convergent mechanisms in RW and SA, we performed similar analysis using two larger validation datasets, GSE103166 (56 RW patients and 21 controls) and GSE65204 (36 SA patients and 33 controls). As expected, fatty acid metabolism-related biological processes and pathways were significantly enriched for both RW and SA ( Supplementary File 1 ; Figure S8 ).

To further screen hub genes most relevant to RW and SA within the key consensus modules, we applied the RF and SVM-RFE machine learning algorithms. Gene expression profiles in the nasal orange module (66 genes) and tracheal plum module (98 genes) were extracted, and the RF and SVM-RFE algorithms were run. Twelve and 14 feature genes in the nasal orange module were identified with the best classification accuracy based on RF and SVM-RFE, respectively ( Figures 1A, C ), while 21 and 15 feature genes in the tracheal plum module were identified with the best classification accuracy based on RF and SVM-RFE, respectively ( Figures 1B, D ). As shown in the UpSet diagram in Figure 1E , five genes (CST1, CST2, CST4, POSTN, and NTRK2) overlapped between the two algorithms. These five genes were also co-DEGs between RW and SA ( Figure 4C ) and were therefore selected as hub genes.

To verify the efficacy of the hub genes, ROC analyses were performed to evaluate the power of each hub gene to distinguish RW and SA patients from healthy controls. Overall, the results showed that each hub gene displayed a moderate ability to distinguish children with RW and SA from healthy controls across both nasal and tracheal samples ( Supplementary File 1 ; Figures S9A–D ).

To further validate transcriptome data accuracy, we analyzed the expression levels of the hub genes using two independent SA gene expression datasets (GSE19187 and GSE65204) and one RW gene expression dataset (GSE103166). As expected, all five genes were significantly up-regulated in both RW and SA patients compared to the controls ( Supplementary File 1 ; Figures S9E–G ). Next, we examined the mRNA expression levels of hub genes using qRT-PCR of BAL samples from 18 controls, 15 children with SA, and 32 preschoolers with RW. Similar results were obtained after comparing the expression levels of hub genes, revealing that CST1, CST2, CST4, POSTN, and NTRK2 were highly overrepresented in preschoolers with RW and patients with SA ( Figure 1F ). Furthermore, CST1, CST2, CST4, and NTRK2 displayed moderate-to-high discriminatory power (AUC = 0.91, 0.77, 0.90, and 0.89, respectively) and POSTN showed high discriminatory power (AUC = 0.92) for distinguishing RW from controls ( Figure 1G ).

We next performed enrichment analysis of TFs to identify potential upstream regulators of the hub genes. All enriched TFs and their corresponding AUC and NES values are shown in Supplementary File 5 ; Table S7 . The TF-gene interaction networks of the top 10 enriched TFs for each hub gene and corresponding hub genes are shown in Supplementary File 1 ; Figure S10 . Results indicated that the master TF regulator of CST1, CST2, and CST4 was TRIM28, which had the highest NES, while the main TF regulator for POSTN and NTRK2 was TBPL2.

Previous studies in adults have shown that POSTN and several other genes can be used to define subsets with T2-high and T2-low asthma  (35–37). Thus, we explored whether the combined expression signatures of the five hub genes could be used to classify RW subjects into different subsets. Using the nasal and tracheal transcriptomes, we performed hierarchical clustering analysis of all subjects with RW based on the expression levels of POSTN, CST1, CST2, CST4, and NTRK2. Approximately 30% of RW patients with high hub gene expression were grouped together in one major branch of the dendrogram in both the nasal ( Figure 2A ) and tracheal samples ( Figure 2B ). We then compared the expression levels of 12 T2 inflammatory markers  (35, 38) between RW patients with high and low hub gene expression. Interestingly, most T2 marker genes (e.g., interleukin-5 (IL-5), IL-13, and IL1RL1) showed significantly higher expression in subjects with high hub gene expression than in those with low hub gene expression ( Figures 2C, D ). Furthermore, we found high positive correlations between the expression levels of hub genes and T2 inflammatory genes ( Figures 2E, F ). Taken together, these results suggest that the identified hub genes could potentially be used for defining subsets (T2-high vs. T2-low) of preschoolers with RW.