Distinct Airway Microbiome and Metabolite Profiles in Eosinophilic and Neutrophilic Asthma

Introduction

Asthma is a chronic and heterogeneous disease characterized by airway inflammation and hyperresponsiveness, manifesting in several inflammatory phenotypes.^1,2 Among these, eosinophilic asthma (EA) and neutrophilic asthma (NEA) are two prominent phenotypes, each associated with distinct immune responses and clinical features.^3,4 These phenotypes differ notably in their pathogenesis and response to corticosteroids. While EA typically responds well to corticosteroid therapy, NEA often exhibits steroid resistance and is frequently observed in severe asthma cases, where the use of macrolides may be necessary.^5,6 Despite these well-characterized clinical distinctions, considerable gaps remain in understanding the underlying microbial and metabolic landscapes of these phenotypes. Furthermore, biomarkers to differentiate EA and NEA remain limited, with clinical assessments largely relying on cell counts in blood and sputum. Identifying novel biomarkers for these distinct asthma endotypes could provide deeper insights into disease mechanisms and guide personalized treatment strategies.

Table 1 Clinical Parameters of Discovery Cohort

Table 2 Clinical Parameters of Validation Cohort

Environmental triggers such as allergens and cold air are well-established contributors to asthma exacerbations.⁷ However, emerging evidence supports the hygiene hypothesis, which suggests that increased microbial exposure may reduce asthma risk.⁸ The respiratory microbiota plays a critical role in shaping immune responses, and even small shifts in microbial composition can significantly impact host immunity due to variations in bacterial ligands and metabolite production.⁹ The respiratory tract harbors a diverse bacterial community, with substantial shifts in richness and diversity along its length, from the nasal cavity to the lungs.¹⁰ These microbial changes have been implicated in several respiratory diseases, including lung cancer,¹¹ asthma¹² and chronic obstructive pulmonary disease (COPD).⁹ Growing evidence shows microbiota significantly influence asthma development, with distinct microbial patterns observed across disease severity, treatments, and age groups. Severe asthma exhibits reduced microbial diversity with more Proteobacteria, while mild cases show altered Firmicutes/Bacteroidetes ratios.¹³ Asthma patients exhibit characteristic alterations in their microbiome following corticosteroid treatment. The oral microbiome demonstrated significant reductions in both Absconditabacteriales and Bifidobacterium, nasal microbiota exhibits the most significant alterations, particularly downregulation of Porphyromonas and Prevotella.¹⁴ Importantly, children with asthma have delayed microbiota development, while elderly patients show Bacteroides-dominated profiles.¹⁵ The specific role of microbiota in the pathogenesis of asthma, particularly in relation to distinct inflammatory phenotypes such as EA and NEA, remains poorly understood. Understanding the relationship between airway microbiota and asthma phenotypes could offer new insights into disease mechanisms and help identify prognostic markers for EA and NEA.

Metabolomics, the comprehensive study of small molecules involved in metabolic processes, has become a powerful tool for discovering biomarkers and elucidating disease mechanisms.¹⁶ Recent research has underscored the interplay between lung microbiota metabolites and respiratory health, showing that microbial-derived metabolites can influence disease progression. Yuan et al highlighted the correlation between the respiratory microbiome and serum metabolomics, suggesting microbial influence on metabolic activity.¹⁷ However, research integrating microbiome and metabolomic data in the context of distinct asthma phenotypes remains limited.¹⁸ Integrating microbiome and metabolomics data, particularly within the lung, could reveal novel pathways involved in EA and NEA pathogenesis and enhance the development of diagnostic and therapeutic strategies.

In this study, we used 16S rRNA sequencing for bacterial profiling and non-targeted metabolomics on sputum from EA, NEA patients, and HC. By comparing microbial and metabolite profiles, we aimed to explore distinct microbial-metabolite correlations. This approach sought unique pathways in asthma subtypes, potentially identifying biomarkers and targets. Our findings offer insights for personalized asthma management.

Materials and Methods

Study Population

Patients with diagnosed asthma were recruited from the First Affiliated Hospital of Guangzhou Medical University between April 2023 and June 2024. Asthma diagnosis followed the Global Initiative for Asthma (GINA) guidelines,¹⁹ according to fulfilment of 1 or more of the following criteria:(1) an increase > 12% in forced expiratory volume in 1 second (FEV1) and an increase > 200 mL after inhalation of 400 μg of albuterol; (2) a 20% reduction in FEV1 in response to a provocative concentration of inhaled methacholine (PC20 < 10 mg/mL); (3) and/or an increase > 20% in FEV1 after 2 weeks of treatment with systemic oral corticosteroids (OCS) or inhaled corticosteroids (ICS).²⁰ (4) Participants were aged between 18 and 70 years. (5) All study participants attended a single visit that encompassed assessments of lung function, asthma symptoms, asthma-specific quality of life, and medication use. We enrolled patients and Healthy control (HC) with no history of smoking. (6) Patients were initially diagnosed with mild-to-moderate asthma exacerbation or persistent asthma according to the American Thoracic Society guidelines and GINA guidelines.²¹ (7) Patients had been routinely treated with fluticasone or an equivalent dose of inhaled corticosteroids for over 6 months but had experienced more than 1 acute exacerbations in the past half year. Exclusion criteria: History of chronic lung diseases or any other diagnosed conditions, ongoing immunotherapy, malignancy diagnosis, current/recent ex-smokers (<1 year), recent respiratory infection, or antibiotic use in past 4 weeks.

Samples Collection

Sputum induction involved isotonic saline with a bronchodilator, processed as in.¹⁴ ICS and LABAs were stopped for 24 hours.²² Asthma patients were classified as EA or NEA based on sputum cell counts:²³ EA had ≥3% eosinophils; NEA had <3% eosinophils and neutrophil thresholds of ≥75.57% (<20 years), ≥61.61% (20–40 years), ≥63.25% (40–60 years), or ≥67.25% (>60 years). Samples with >10% squamous cells were excluded. Precipitates were used for microbiome analysis, supernatants stored at −80°C for metabolism analysis.²⁴ All participants underwent spirometry per ERS guidelines to assess lung function and airway obstruction severity. Included patients had mild-to-moderate, clinically stable exacerbations safe for bronchoscopy.

Ethical approval was obtained from the Ethics Committee of The First Affiliated Hospital of Guangzhou Medical University, with approval codes 2022 No.121 and 2024 No. G-007. Written informed consent was obtained from all participants prior to inclusion in the study. This study was conducted in accordance with the ethical principles of the Declaration of Helsinki. A schematic representation of the study can be found in Figure 1.

Figure 1 Visual representation of the study design. A cohort study combining microbiome analysis, metabolomics analysis, and machine learning applications to identify and validate biomarkers and pathways associated with asthma phenotypes and disease mechanisms.

16S rRNA Amplicon Sequencing

Microbial DNA from sputum was extracted using the FastDNA Spin Kit and purity was checked with a Nanodrop. Integrity was assessed via agarose gel electrophoresis. PCR amplified the V3-V4 regions of the 16S rRNA gene according to universal primers (515F: 5′-GGACTACHVGGGTWTCTAAT-3′, 806R: 5′- GGACTACHVGGGTWTCTAAT- 3′). PCR products were purified, quantified, and sequenced on an Illumina NovaSeq6000 at the Novogene Bioinformatics Technology Co., Ltd. The paired-end reads were assigned using barcode and then were merged with FLASH (Fast Length Adjustment of SHort reads) (V1.2.11) software.²⁵ High-quality filtering of the raw tags was performed to obtain clean tags using QIIME (Quantitative Insights into Microbial Ecology) (Version QIIME2-202202). After that, chimeric sequences were compared with Silva Database and using the UCHIME algorithm to detect and remove chimera sequences. A final effective sequence analysis was performed using Uparse software (Uparse V7.0.1001).²⁶

Microbiomics Study

Sequences with a similarity index greater than 97% were assigned to the same operational taxonomic units (OTUs). Alpha diversity metrics (Chao1, Shannon Diversity, Simpson, and Richness indice) were calculated for the OTU table using alpha rarefaction. The beta diversity metric was calculated on the OTU table and visualized using principal coordinate analysis (PCoA) and non-metric multidimensional scaling (NMDS) analysis. The alpha diversity and beta diversity analysis were calculated with QIIME software (version 1.9.1) and displayed with R software (version 2.15.3). Metagenomeseq test was used to calculate significant differences between two groups. Metabolic function was predicated with Phylogenetic Investigation of Communities by Reconstruction of Unobserved States (PICRUSt) based on the high-quality tags of the microbiota.²⁷ Random forest analysis was the classification of samples based on a machine learning model for the purpose of filtering variables.

Non-Targeted Metabolomics Profiling

Metabolites were extracted from sputum and tested with a liquid chromatography-tandem mass spectrometry (LC-MS/MS). The metabolomics analysis was performed using a Vanquish UHPLC system (ThermoFisher, Germany) in Novogene Co., Ltd. (Beijing, China). Samples were injected onto a Hypersil Gold column (100×2.1 mm, 1.9μm) using a 12-min linear gradient at a flow rate of 0.2 mL/min. The eluents for the positive and negative polarity modes were eluent A (0.1% FA in Water) and eluent B (Methanol). The solvent gradient was set as follows: 2% B, 1.5 min; 2–85% B, 3 min; 85–100% B, 10 min; 100–2% B, 10.1 min; 2% B, 12 min. Q ExactiveTM HF mass spectrometer was operated in positive/negative polarity mode with spray voltage of 3.5 kV, capillary temperature of 320°C, sheath gas flow rate of 35 psi and aux gas flow rate of 10 L/min, S-lens RF level of 60, Aux gas heater temperature of 350°C. The raw data files generated by UHPLC-MS/MS were processed using the Compound Discoverer 3.3 (CD3.3, ThermoFisher) to perform peak alignment, peak picking, and quantitation for each metabolite.

Metabolomics Study

Principal component analysis (PCA) and orthogonal partial least squares-discriminant analysis (OPLS-DA) were used to explore whether all samples can be significantly clustered in different groups. We applied univariate analysis (t-test) to calculate the statistical significance (P-value). The metabolites with VIP > 1 and P-value< 0.05 and fold change≥1.2 or FC≤0.833 were considered to be differential metabolites.^28–30 Volcano plots were used to filter metabolites of interest which based on log2 (Fold Change), and -log10(P-value) of metabolites by ggplot2 in R language P-value < 0.05 was considered as statistically significant and correlation plots were plotted by corrplot package in R language. The functions of these metabolites and metabolic pathways were studied using the KEGG database.

Statistical Analysis

Statistical significance of the difference was firstly estimated using the F-test, followed by the Student’s or Welch’s t-test in two-group analysis according to the homogeneity of variance. Data normality was assessed by the Kolmogorov–Smirnov method. Comparisons between the two groups were performed by t-test with homogeneity of variance conversely. Mann–Whitney-Wilcoxon test was performed for the dataset, which does not follow the normality. Differences in airway bacteria diversity and relative abundance among three groups were analyzed by the nonparametric Kruskal–Wallis test. P-value of < 0.05 were considered statistically significant. All statistical analysis about the identification of distinct metabolites was completed with R (v4.0.2) basic statistical packages.

To evaluate the discriminatory power of classifiers in distinguishing between healthy controls and patients, the Receiver Operating Characteristic (ROC) Curve was employed. Based on the differential metabolites and microbiota obtained by the above analysis, ROC analysis was used to obtain curve and calculate the area under the curve (AUC), providing a quantitative measure of the model’s performance. Furthermore, we ventured to enhance the predictive capability of our model by integrating the identified biomarkers.

Results

Characteristics of All Participants

Two cohorts were established to explore asthma heterogeneity based on induced sputum cytology: a discovery cohort with 44 asthma patients (24 with eosinophilic asthma, EA, and 20 with neutrophils asthma, NEA) and 20 healthy controls (HC), and a validation cohort consisting of 26 asthma patients (10 EA, 8 NEA) and 8 hC.

Baseline characteristics revealed no statistically significant differences in age, gender distribution, body mass index (BMI) across the groups. Serum total IgE levels were significantly higher in asthmatic patients compared to HC, with EA showing higher levels than NEA. Lung function analysis demonstrated significant intergroup differences: the percentages of FEV1 (P = 0.018), FEV1/FVC ratio (P = 0.026), and FVC (P = 0.013) were higher in the HC group than in both EA and NEA groups. Furthermore, these lung function indicators were higher in EA compared to NEA. EA exhibited significant eosinophilia in both blood and sputum, accompanied by lower neutrophil counts compared to NEA. NEA also had the longest hospital stay compared to the other groups. Clinical demographics of the discovery cohort are presented (Table 1), and the validation cohort showed high consistency with the discovery cohort (Table 2).

Microbiota in the Respiratory Tract of Different Asthma Patients versus HCs

Estimation of Sequencing Depth

To investigate these dynamics, 16S rRNA amplicon sequencing of sputum samples was performed to explore shifts in lung microbiota and uncover alterations in the structure and function of the lower respiratory tract microbiota. Based on Research statistics, raw 16S rRNA sequencing data using primers 515F and 806R were filtered using QIIME software, followed by sequence splicing with FLASH software to generate tag sequences. A total of 64 high-quality samples were obtained from sputum sequencing, with an average of 87,495 valid sequences per sample. At a 97% similarity threshold, 2696 Operational Taxonomic Units (OTUs) were identified. Rarefaction curves for each group indicated that sequencing data reached saturation, suggesting sufficient depth of coverage (Figure S2A). Rank abundance curves provided detailed insights into the richness and evenness of microbial diversity across the sample sets (Figure S2B). Additionally, the species accumulation boxplot demonstrated a steady increase in species diversity as sample size increased, with curves plateauing at the maximum sample size of 64, further supporting the adequacy of sequencing depth (Figure S2C).

Alpha-Diversity and Beta-Diversity

Alpha-diversity was assessed to analyze microbial community diversity within individual groups. Metrics such as Chao1, Shannon Diversity, Simpson, and Richness indices did not reveal significant differences in species richness or diversity, indicating that α-diversity remained consistent across the three study groups (P > 0.05, Kruskal–Wallis tests) (Figure 2A). In contrast, beta-diversity, which compares microbial community composition across different groups, demonstrated significant differences between EA and NEA (P = 0.03, Wilcoxon test, unweighted UniFrac distance), and between NEA and HC (P = 9e-04, Wilcoxon test, unweighted UniFrac distance) (Figure 2B). These findings were supported by the Principal Coordinates Analysis (PCoA) plot, which illustrated shifts in microbial community composition among the groups (Figure S1A). The PCoA, based on Bray–Curtis dissimilarities, highlighted statistically significant differences in overall microbiota structure between HC and NEA, as well as between EA and NEA. Non-metric multidimensional scaling (NMDS) analysis, based on Bray-Curtis similarity distance, further demonstrated distinct group separation (stress = 0.089, Figure S1B). A Venn diagram (Figure 2C) visually depicted the shared and unique OTUs among the groups, providing insights into their compositional differences.

Figure 2 Characteristics of the airway microbiome composition. (A) Alpha diversity analysis. (B) Microbial community beta-diversity. (C) Venn diagram of OTUs. (D) Relative abundance of top 10 phylum and genus. (E) Heatmap of microbial relative abundances of top 10 phylum. (F) Heatmap of microbial relative abundances of top 10 genus. *P <0.05, **P<0.01.

Distribution of Microbiota Taxonomic Composition in Different Asthma Patients

All sequences were classified into 30 phyla. The relative abundance of the top 10 phyla across the three groups is displayed in Figure 2D. The most abundant phyla in the airway microbiota were Proteobacteria, Firmicutes, Bacteroidetes, and Fusobacteriota. EA patients exhibited a higher abundance of the Proteobacteria phylum, while NEA patients showed higher abundances of Bacteroidetes, Aenigmarchaeota, and Firmicutes. However, there were no significant differences in phylum-level distribution between HC and asthma patients (P > 0.05, MetagenomeSeq, Figure 2E). Despite similar phylum-level distributions, genus-level analysis revealed distinct patterns. Both healthy controls and asthma patients shared four major bacterial general: Stenotrophomonas, Streptococcus, Achromobacter, and Neisseria. Notably, Veillonella was more abundant in NEA patients compared to HC, while Achromobacter was more prevalent in NEA patients than in EA patients (P < 0.05, MetagenomeSeq, Figure 2F). These results underscore the nuanced differences in microbial composition at the genus level, despite similarities at the broader phylum level.

Altered Composition of the Respiratory Tract in Different Asthma Patients

We utilized the MetagenomeSeq method to identify genus-level taxa with significant statistical and biological differences among the three groups. In the comparison between EA and HC, Caulobacter and Acidovorax were more abundant in EA than in HC (P < 0.05, Figure 3A). In the comparison between NEA and HC, a marked enrichment of Enterococcus, Pseudomonas, Acinetobacter, Rothia, Porphyromonas, and Veillonella was observed in the NEA group (P < 0.05, Figure 3B). Further comparison between EA and NEA revealed that Achromobacter, Abiotrophia, and Rothia were more abundant in NEA compared to EA (P < 0.05, Figure 3C). This differential abundance pattern highlights the unique microbial landscapes characterizing EA and NEA, which may contribute to their distinct clinical manifestations and disease mechanisms.

Figure 3 Prediction model of airway microbiota for EA and NEA based on genus-level abundances using random forests. (A–C) Genus-level taxa differed between groups. (D) Variable importance ranking chart. (E) ROC curve of the model. (F) COG functional analysis of microbiota between NEA and HC. (G) COG functional analysis of microbiota between EA and NEA. *P <0.05, **P<0.01.

With the confirmation of airway microbial dysbiosis among the patients, we next explored the potential of the airway microbiota as diagnostic indicators. Using random forest analyses, we sought to identify taxonomic markers that could distinguish between EA and HC, as well as between NEA and the other two groups. This analysis constructed a random forest predictive model based on OTU feature selection. By selecting general with the highest area under the curve (AUC) values, we identified key markers (Figure S3). For distinguishing EA from HC, five general were selected based on MeanDecreaseAccuracy (Figure 3D), and a 10-fold cross-validation was applied. ROC curves were plotted to evaluate predictive performance, yielding an AUC of 66.78%, indicating moderate discriminatory ability (Figure 3E). However, the model demonstrated robust performance in distinguishing NEA from both HC and EA. Remarkably, only three general were required to achieve the highest AUC value of 83.81% when distinguishing NEA from HC, underscoring strong discriminatory potential. These findings suggest that the airway microbiota holds significant promise as a diagnostic tool, particularly in identifying NEA patients from healthy controls. Additionally, Haemophilus, Methylobacterium, and Rothia emerged as potential biomarkers for identifying individuals predisposed to developing neutrophils asthma.

Functional predictions of microbial communities based on 16S rRNA gene sequences were conducted using PICRUSt2 with the COG database. Distinct functional categories within the COG database were identified between NEA and HC groups, reaching statistical significance (P < 0.05, Figure 3F), including DNA-binding response regulators and glycosidases. A comparative analysis between EA and NEA microbial communities also revealed statistically significant differences (P < 0.05, Figure 3G) in five key functional domains, including DNA-binding transcriptional regulators, ADP-ribose pyrophosphatase YjhB, predicted dehydrogenase, signal transduction histidine kinase, and DNA-binding response regulators. These functional differences provide valuable insights into the structure and functionality of microbial communities and their interactions under distinct environmental conditions.

Metabolomics Analysis of Different Asthma Patients versus HC

Multivariate Analysis of Metabolomic Data

To further elucidate the pathogenesis of asthma in different patient subtypes, residual sputum samples were analyzed for metabolomics. In positive ion mode (ESI+), 807 metabolites were identified, with 252 annotated against the KEGG database. In negative ion mode (ESI-), 532 metabolites were identified, with 193 annotated. The PCA score plots in both 2D and 3D for the sputum metabolome in ESI+ and ESI- modes (Figure S4A and B) demonstrated strong stability across samples in both ionization modes, confirming the reliability of the metabolomic system. A clear segregation was observed between EA, NEA, and HC groups, indicating significant metabolic differences. Each sample was represented as a point on the score plots, and permutation tests (n = 200) were used to validate the OPLS-DA model. The permutation plots showed that all blue Q2-values to the left were lower than the original points on the right, confirming the validity of the models. OPLS-DA clearly separated metabolites between the groups under investigation (Figure S4C–H), highlighting distinct metabolic profiles across asthma subtypes.

Identification of the Differential Metabolites in Sputum

A combined analysis of positive and negative ions was conducted to comprehensively understand metabolite functions. Metabolites that distinguished asthma patients from controls were selected based on fold change (cutoff >1.2 or <0.833) and VIP scores from OPLS-DA.

A total of 269 metabolites showed significant differences between EA patients and HC (Figure 4A, D, G), with 180 metabolites upregulated and 89 downregulated in EA patients (VIP > 1, P < 0.05). Notably, lipids and lipid-like molecules such as 15(S)-HpETE, Lipoxin B4, and 15-Epiprostaglandin E1 were elevated in EA patients, while glycerophospholipids such as Thymidine 5′-monophosphate and Uridine monophosphate (UMP) were downregulated.

Figure 4 Metabolic profiles among three groups. (A–C) PLS-DA plot for pairwise comparison between groups (D–F) Volcano plot of metabolites of pairwise comparison between groups. (G–I) Matchstick chart of the top 20 different metabolites between EA and HC, NEA and HC, EA and NEA.

For NEA patients, 224 metabolites exhibited significant differences compared to HC (Figure 4B, E, H), with 149 metabolites upregulated and 75 downregulated (VIP > 1, P < 0.05,). Notable metabolites in NEA included elevated levels of PC 16:0_16:1 and D-Raffinose, while methadone-d9 and 2-Methylglutaric acid were downregulated. In the comparison between EA and NEA groups, 166 metabolites showed significant differences in relative abundance (Figure 4C, F, I), with 88 upregulated and 78 downregulated in EA patients (VIP > 1, P < 0.05). Key metabolites such as 5-Hydroxymethyluracil, methadone-d9, and 5-Hydroxyindole-3-acetic acid were elevated in EA patients, while Quinolinic acid and PC 16:0_16:1 was downregulated in NEA patients.

Enrichment Analysis and Correlation of Metabolites in Sputum

The Kyoto Encyclopedia of Genes and Genomes (KEGG) database provides a comprehensive resource for the functional interpretation of genomes and high-throughput data generated from large-scale experimental technologies. The altered metabolites in this study were primarily associated with amino acid and lipid metabolism. In the comparison between EA and HC, the significantly different metabolites were enriched in 51 metabolic pathways, with three pathways showing significant differences: pyrimidine metabolism, phenylalanine metabolism, and bile secretion (P < 0.05) (Figure 5A). For NEA versus HC, 46 pathways were identified, with four showing significant differences, including tryptophan metabolism, protein digestion and absorption, and phenylalanine metabolism (P < 0.05) (Figure 5B). In the comparison between EA and NEA, 38 pathways were enriched, with four pathways showing significant differences, including arachidonic acid metabolism, pyrimidine metabolism, serotonergic synapse, and leishmaniasis (P < 0.05) (Figure 5C).

Figure 5 Enrichment Analysis and correlation of metabolites in sputum. (A–C) KEGG pathway enrichment analysis of the top 20 differential metabolites between EA and HC, NEA and HC, EA and NEA. (D–F) Correlation analysis of the top 20 differential metabolites between EA and HC, NEA and HC, EA and NEA.

Correlation analysis helps to elucidate the synergistic or mutually exclusive relationships between different metabolites, providing insight into their regulatory interactions during biological state changes. A heat map of the top 20 differential metabolites between EA and HC shows that most lipids and lipid-like molecules were positively correlated (Figure 5D). In the comparison between NEA and HC, 5-hydroxyindole-3-acetic acid, derived from indoles and derivatives, showed a positive correlation with 5-hydroxymethyluracil, derived from diazines, while no strong correlations were observed among other metabolites (Figure 5E). In the EA versus NEA comparison, fatty acyls such as 5(S),15(S)-DiHETE exhibited positive correlations with Protectin D1, a lipid mediator with anti-inflammatory and immunoregulatory effects (Figure 5F). The positive correlation between these two metabolites may suggest a shared inflammatory response or immunoregulatory mechanism between the EA and NEA groups.

Random Forest Analysis Discovered Potential Metabolic Biomarkers for Asthma Patients’ Stratification

In the comparison between EA and HC, KEGG enrichment analysis identified three significantly different pathways—pyrimidine metabolism, phenylalanine metabolism, and bile secretion (P < 0.05). Among these, pyrimidine metabolism had the smallest P-value, indicating the highest reliability. Consequently, we established a random forest model based on the eight differential metabolites enriched in the pyrimidine metabolism pathway to identify metabolic biomarkers capable of distinguishing EA from HC. The mean decrease Gini score highlighted the prominent metabolites differentiating EA from HC (Figure 6A). The top metabolites associated with EA included methylmalonate, 5-methylcytosine, dUDP, and 2-deoxyuridine (Figure 6B). To enhance the accuracy of biomarker identification, these four metabolites were combined in a ROC analysis, achieving an AUC of over 0.90, significantly higher than that of individual metabolites (Figure 6C).

Figure 6 Metabolic markers of EA and NEA by random forest models. (A–C) Candidate metabolites classifying EA and HC. (D–F) Candidate metabolites classifying NEA and HC. (G–I) Candidate metabolites classifying EA and NEA. (J) ROC curve analysis of the candidate biomarkers.

In the comparison between NEA and HC, tryptophan metabolism exhibited the smallest P-value. We developed a random forest model using the seven differential metabolites enriched in the tryptophan metabolism pathway to identify biomarkers distinguishing NEA from HC (Figure 6D). The top metabolites associated with NEA included quinolinic acid, 5-hydroxyindole-3-acetic acid, acetyl-N-formyl-5-methoxykynurenamine, and 2-ketoadipic acid (Figure 6E), with a combined AUC of 0.90 in the test set (Figure 6F).

In the comparison between EA and NEA, the arachidonic acid metabolism pathway showed the most significant P-value. A random forest model was established with nine differential metabolites from this pathway to identify metabolic biomarkers distinguishing EA from NEA (Figure 6G). The top metabolites associated with NEA included prostaglandin E2, 15(S)-HpETE, PC 34:1, PC 32:0, and 11,12-epoxy-(5Z,8Z,11Z)-icosatrienoic acid. The expression of these metabolites in each group is shown in Figure 6H, with a combined AUC of 0.93 in the test set (Figure 6I).

To further validate the metabolites identified through non-targeted metabolomics, an external validation cohort was analyzed. Consistent with the discovery cohort results, significant differences were observed, particularly in the comparison between EA and HC. The four metabolites involved in pyrimidine metabolism—methylmalonate, 5-methylcytosine, dUDP, and 2-deoxyuridine—showed consistent high AUC results in the validation cohort. Additionally, using these four metabolites to differentiate NEA from HC achieved an AUC of 0.91 in the validation set, demonstrating excellent diagnostic performance. Similar findings were confirmed in comparisons between EA and NEA (Figure 6J).

Correlation Between Airway Microbiome and Metabolites

To investigate the relationship between the altered airway microbiome and metabolites involved in asthma pathogenesis, we calculated Spearman correlation between differential bacteria and the top 40 differential metabolites across the different groups. In the comparison between EA and HC, we found that the Caulobacter genus positively correlated with the metabolite 6-(3-hydroxybutan-2-yl)-5-(hydroxymethyl)-4-methoxy-2H-pyran-2-one, though no specific functional annotation was available in KEGG (Figure 7A).

Figure 7 Spearman correlation heatmap for significantly altered metabolites and general among the three groups. (A) Correlation heatmap between EA and HC. (B) Correlation heatmap between NEA and HC. (C) Correlation heatmap between EA and NEA. *P < 0.05, **P < 0.01.

In the comparison between NEA and HC, all six different general exhibited positive correlations with metabolites (Figure 7B). Notably, Veillonella, Acinetobacter, Pseudomonas, and Enterococcus were positively associated with LPE 18:0, 2-linoleoyl glycerol, and cis-aconitic acid. LPE 18:0 and 2-linoleoyl glycerol are lipids and lipid-like molecules, while cis-aconitic acid is an organic acid involved in the 2-oxocarboxylic acid metabolism pathway. The Porphyromonas genus correlated positively with salicylic acid, a benzenoid involved in the phenylalanine metabolism pathway. The Rothia genus correlated positively with metabolites such as N“, N“′-di[1-(4-nitrophenyl)ethylidene]carbonic dihydrazide, 1-(4-fluorophenyl)-2-(4-methoxyphenyl)-4-(2-naphthyl)butane-1,4-dione, and dihydrothymine, while showing negative correlations with Lypressin acetate and 2-methoxyestrone.

In the comparison between EA and NEA, interesting correlations were observed between specific bacterial general and distinct metabolites. Abiotrophia and Achromobacter general exhibited a positive correlation with thymidine 5′-monophosphate, a nucleoside involved in the antifolate resistance pathway (Figure 7C). The Abiotrophia genus also showed a negative correlation with 3-methyluridine, another nucleoside. Achromobacter positively correlated with (±)9-HpODE, suggesting a potential link between this genus and specific lipid peroxidation products. Additionally, Rothia correlated positively with N1-cyclohexyl-4-(4-fluorobenzyl)-1,4-diazepane-1-carbothioamide. These findings highlight the intricate interplay between airway microbial communities and their metabolic byproducts in different asthma phenotypes, potentially contributing to the heterogeneity observed in asthma pathogenesis. Additionally, we compared the correlations between differentially microbial/metabolites and granulocyte percentages in EA and NEA. The correlation analysis revealed that Abiotrophia positively correlated with sputum neutrophils. Achromobacter showed strong associations with blood neutrophil counts (Figure S5A). These findings support EA and NEA as distinct endotypes, with Abiotrophia and Achromobacter linked to neutrophilic inflammation. The differential metabolites 15(S)-HpETE and 11,12-Epoxy-(5Z,8Z,11Z)-icosatrienoic acid exhibited significant positive correlations with eosinophils in both sputum and blood (Figure S5B). The compartment-specific (sputum vs blood) and phenotype-dependent patterns highlight potential biomarkers for targeted therapies.

Discussion

This study highlights the pivotal role of the lower airway microbiota in asthma pathogenesis, aligning with the growing body of research on the microbiome’s influence on asthma outcomes.^12,31,32 By integrating 16S rRNA sequencing with untargeted metabolomics, we sought to uncover distinct microbial and metabolic signatures across EA, NEA, and HC, advancing our understanding of asthma’s heterogeneity.

At the phylogenetic level, our study provides robust evidence for the existence of a core respiratory microbiome, comprising Proteobacteria, Firmicutes, Bacteroidota, and Fusobacteriota, which serves as a foundational component of respiratory health across both healthy individuals and those with asthma. This core microbiome underscores the critical role of these microbial groups in maintaining respiratory homeostasis, regardless of asthma status. However, genus-level differences revealed significant shifts in microbial composition, particularly between EA and NEA. The enrichment of Caulobacter and Acidovorax in EA suggests these general may contribute to eosinophilic inflammation.^33,34 The presence of these bacteria may influence the local immune response, promoting the recruitment and activation of eosinophils, which are central to the inflammatory process in EA. While the increased abundance of Enterococcus, Pseudomonas, Acinetobacter, Rothia, Porphyromonas, and Veillonella in NEA reflects a microbial signature commonly associated with infection and inflammation. This microbial signature is commonly associated with infection and inflammation, indicating a potential role for these bacteria in driving the inflammatory pathways characteristic of NEA.

Functionally, the identification of two significantly altered classes-DNA-binding response regulators and glycosidases-between NEA and HC underscores the profound impact of NEA on cellular processes. DNA-binding proteins, as transcriptional regulators, influence gene expression related to inflammation, immune responses, and cell proliferation.³⁵ Similarly, glycosidases, which are key enzymes involved in glycan processing, also exhibit significant alterations in NEA. This suggests that NEA may disrupt the normal glycan metabolism, further illustrating the complexity of cellular changes associated with this asthma subtype.³⁶ The interplay between glycosidases and other cellular components may be crucial in understanding the pathogenesis of NEA and developing targeted therapeutic strategies.³⁷ In addition, the significant differences in five functional categories, including transcriptional regulators, ADP-ribose pyrophosphatase YjhB, predicted dehydrogenases, signal transduction histidine kinases, and DNA-binding response regulators, between EA and NEA, span a wide range of metabolic and regulatory processes, reflecting the sophisticated nature of microbial interactions within each asthma phenotype. These findings reveal the sophisticated nature of the microbial interactions within each phenotype and the need for targeted therapeutic strategies to address their distinct pathogenic mechanisms.

Our metabolomics analysis revealed significant metabolic differences between EA and NEA, providing new molecular-level evidence for understanding the heterogeneity of asthma. In EA, the upregulation of pro-inflammatory lipids such as 15(S)-HpETE and lipoxin B4 is highly consistent with the previously reported Th2-type inflammatory characteristics in the airways.³⁸ The downregulation of nucleic acid precursors such as thymidine 5′-monophosphate and UMP reflects metabolic reprogramming in airway epithelial cells under chronic inflammatory conditions, aligning with the mitochondrial dysfunction observed in EA patients.³⁹ NEA exhibits distinct metabolic features: the increase in PC 16:0_16:1 may affect airway smooth muscle reactivity by altering cell membrane fluidity, a mechanism validated in the lipidomics study.⁴⁰ Meanwhile, the downregulation of methadone-d9 and 2-methylglutaric acid may indicate NEA-specific tryptophan metabolism disorders. These findings deepen our understanding of asthma endotyping: EA is dominated by arachidonic acid metabolism, whereas NEA involves phospholipid remodeling and microbial co-metabolism, establishing a robust biochemical basis for the clinical and immunological differences between these two major asthma endotypes.

In EA, altered pyrimidine metabolism was among the most significantly affected pathways, indicating its role in regulating inflammatory responses and cellular proliferation.⁴¹ Altered pyrimidine metabolism played a key role, with four metabolites (Methylmalonate, 5-Methylcytosine, dUDP, 2-Deoxyuridine) identified as robust biomarkers. In NEA, alterations in tryptophan metabolism were the most significant, suggesting that modulation of this pathway may be key to understanding NEA-associated metabolic changes. Tryptophan breakdown products play crucial roles in the synthesis of bioactive compounds like serotonin and melatonin, which influence immune function, mood, and sleep regulation.^42,43 Additionally, the involvement of protein digestion and absorption pathways in NEA indicates possible disruptions in nutrient processing, further highlighting its metabolic complexity. Biomarkers for NEA, including quinolinic acid, 5-hydroxyindole-3-acetic acid, acetyl-N-formyl-5-methoxykynurenamine, and 2-ketoadipic acid, demonstrated strong predictive value. Arachidonic acid metabolism distinguished EA and NEA, with potential therapeutic targets identified.^44,45 Altered arachidonic acid metabolites in allergic rhinitis serum suggested a metabolic association.¹⁶ The consistency of these metabolic markers across validation cohorts strengthened their reliability and clinical significance. Overall, these findings provided insights into metabolic pathways and biomarkers for targeted therapies and diagnostics.

In exploring the link between airway microbiota alterations and metabolite changes in disease pathogenesis, our Spearman correlation analysis revealed intriguing findings. The airway microbiota may impact the host’s metabolic state, influencing asthma onset and inflammation development. Comparisons showed broader microbe-metabolite correlations in EA vs HC, with specific general correlated with lipids. In asthmatics with EA vs NEA, different general were linked to nucleoside metabolism and other pathways. These results suggest microbial communities influence disease progression by modulating metabolic pathways, potentially affected by microbiota structure and nutritional homeostasis. Additionally, we also investigated how the differential microbial communities and metabolites correlate with granulocyte levels in EA and NEA. Although no statistically significant correlations were observed between the differential microbial communities and neutrophil levels in sputum or blood, likely due to the small sample size and the low abundance of the differentially abundant microbial taxa, a positive correlation trend was still noticeable. This suggests that airway microbial dysbiosis may primarily be involved in neutrophil-dominated inflammatory pathways in asthma, whereas key metabolites are mainly positively correlated with eosinophils in both blood and sputum, with lipid mediators orchestrating systemic recruitment of eosinophils. For instance, the 15(S)-HpETE and epoxy fatty acid metabolic pathways, known to modulate type 2 inflammation,⁴⁶ have been implicated, thereby revealing important mechanistic links between microbial/metabolic features and granulocyte-driven inflammation in asthma phenotypes.

Our multi-omics approach has offered novel insights into complexities of asthma, particularly by highlighting differences in airway microbial communities, functional characteristics, and metabolic pathways between EA and NEA. However, the journey toward identifying robust biomarkers for asthma stratification and personalized treatment is fraught with challenges, such as sample size limitations and methodological variability persist, necessitating further validation in larger cohorts to confirm reproducibility and clinical utility. Despite these limitations, our study represents a crucial step toward personalized asthma treatment, laying the foundation for targeted strategies. Future research should focus on elucidating the causal mechanisms underlying microbiome-metabolite interactions in asthma to identify robust biomarkers for more precise and effective management.

Conclusion

This study highlights the profound differences in airway microbial community structure, functional characteristics, and metabolite levels between eosinophilic asthma, neutrophils asthma, and healthy controls. By integrating 16S rRNA sequencing, metabolomics analysis, and the random forest algorithm, we have gained significant insights into the pathogenesis of asthma. This comprehensive approach not only enhances our understanding of the disease’s complexity but also lays a foundation for advancing personalized approaches in asthma care, facilitating more precise and targeted disease characterization and monitoring.