In case–control study, 30 children with asthma and 30 controls without allergic diseases were recruited from Nanjing Children’s hospital affiliated to Nanjing Medical University and Jiangsu Province Hospital of Chinese medicine. Covariates were considered in this study including age, sex (male vs. female), BMI (body mass index), passive smoking (yes vs. no). Inclusion criteria: The children with asthma were diagnosed by physicians without taking antibiotics during the relief stage in recent three months. The control group did not have any respiratory system-related diseases or without allergic disease. Exclusion criteria: (1) subjects with recent surgical history were excluded; (2) subjects with recent Oral problems were excluded; (3) subjects with underlying diseases were excluded; (4) subjects with overweight or obesity were excluded. The demographic information of the children with asthma in the study was in Supplementary Table S1. Oral swabs and fecal samples were collected and stored at −80°C for further experiments. Blood samples were collected from control individuals. Peripheral blood mononuclear cells (PBMCs) were isolated using human lymphocyte separation medium (Fcmacs biotech co., ltd.). This study was approved by the ethics Committee of Nanjing Medical University [approval no. Nan med Univ ethical review (2019)02080-1]. Written informed consent was obtained from all guardians.

Total genomic DNA was extracted from oral swabs and fecal samples according to the E.Z.N.A.® soil kit (Omega Biotek, Norcross, GA). The V3-V4 variable region of the bacterial 16S rRNA gene was amplified by PCR, using the forward primer 338F (5’-ACTCCTAC GGGAGGCAGCAG-3′) and reverse primer 806R (5’-GGACTACHVG GGTWTCTAAT-3′). The PCR amplification system included 4 μL 5 × FastPfu buffer, 2 μL 2.5 mM dNTPs, 0.8 μL of each primer (0.5 μM), and 10 ng DNA template. PCR conditions were as following: 95°C denaturation for 3 min, 27 cycles (95°C for 30 s, 55°C for 30 s and 72°C for 45 s), finally 72°C for 10 min. The PCR products were purified and recovered by AxyPrep DNA Gel Extraction Kit (Axygen Biosciences, Union City, CA, USA). Sequencing was performed using the Miseq PE300 platform (Illumina, San Diego, USA). The raw sequences were quality controlled and assembled using Trimmomatic and FLASH software. The operational taxonomic units (OTUs) were clustered based on 97% similarity using UPARSE (version 7.1¹) software, and the OTUs sequences were annotated using the RDP classifier ² software. The free online platform of the Majorbio I-Sanger Cloud Platform (Shanghai Majorbio Bio-pharm Technology Co., Ltd) ³ was used for data analysis.

Metagenomic sequencing was performed on the Illumina NovaSeq platform (Illumina, USA). The raw reads from metagenome sequencing were used to generate clean reads using the fastp ⁴ on the free online platform of Majorbio Cloud Platform. ⁵ These high-quality reads were then assembled to contigs using MEGAHIT, ⁶ which makes use of succinct de Bruijn graphs. Contigs with the length being or over 300 bp were selected as the final assembling result. Open reading frames (ORFs) in contigs were identified using MetaGene. ⁷ Representative sequences of non-redundant gene catalog were annotated based on the NCBI NR database using blastp as implemented in DIAMOND v0.9.19 with e-value cutoff of 1e-5 using Diamond ⁸ for taxonomic annotations. Functional differences were found by LEfSe.

The contigs, which have been assembled, were aligned with the reference gene sequences of NCBI strains using Blast software (version 2.9.0). Sequences with a similarity ≥80% and a length ≥ 200 bp were selected. The selected sequences were then subjected to gene annotation using Prokka software (version 1.13). The annotation results for the sequences identified as 16S rRNA gene were further aligned using the MUSCLE method in MEGA software (version 11.0.13), and the genetic distance was calculated using the p-distance method.

Epitopes extracted from the Immune Epitope Database (IEDB) ⁹ were identified as triggers for immune reaction in asthma. These epitopes’ amino acid sequences were then compared to the protein sequences of species. Only epitopes with a minimum identity and coverage of 80% were chosen for subsequent experimental investigation. These microbial mimicry peptides were synthesized by Nanjing Branch of GenScript Biotech Company. PBMCs were isolated from anticoagulated blood using human lymphocyte separation fluid and re-suspended in R1640 complete medium (Zhang et al., 2021). Peptide segments were added to the cells for stimulation at a final concentration of 40 μg/mL, and the cells were cultured at 37°C and 5% CO₂ for 48 h before collection. Total RNA was extracted from the cells using the TRIzol method (Tiangen, China), and the concentration was measured using a nucleic acid protein analyzer (Thermo NanoDrop 2000, US). Reverse transcription was performed using the HiScript II Q RT SuperMix kit (Vazyme Biotech) following the instructions provided. Reverse-transcription quantitative polymerase chain reaction (RT-qPCR) was used to quantify asthma-related inflammatory factors IL-4, IL-5, IL-6, IL-8, IL-13, and IL-17A. The PCR system consisted of cDNA, primer (0.2 μmol/L) and ChamQ Universal SYBR qPCR Master Mix (5 μL), with a total system volume of 10 μL. The primers were synthesized by Nanjing Branch of GenScript Biotech Company, with the following sequences: For IL-4, the forward primer was 5’-CGGCAACTTTGTCCACGGA-3′, and the reverse primer was 5’-TCTGTTACGGTCAACTCGGTG-3′; for IL-5, the forward primer was 5’-GGAATAGGCACACTGGAGAGTC-3′, and the reverse primer was 5’-CTCTCCGTCTTTCTTCTCCACAC-3′; for IL-6, the forward primer was 5’-CCTTCGGTCCAGTTGCCTTCTC CCT-3′, and the reverse primer was 5’-GGGCTGAGATGCCGTCGA GGATGTA-3′; for IL-8, the forward primer was 5’-CCACCGGAAG GAACCATCTC-3′, and the reverse primer was 5’-TTCCTTGGG GTCCAGACAGA-3′; for IL-13, the forward primer was 5′- CCTCATGGCGCTTTTGTTGAC-3′, and the reverse primer was 5′- TCTGGTTCTGGGTGATGTTGA-3′; for IL-17A, the forward primer was 5’-AGATTACTACAACCGATCCACCT-3′, and the reverse primer was 5’-GGGGACAGAGTTCATGTGGTA-3′; for the reference gene β-actin, the forward primer was 5’-CATGTACGTTGC TATCCAGGC-3′, and the reverse primer was 5’-CTCCTTAATGTCA CGCACGAT-3′. According to the instructions of the kit, PCR reactions were performed in a LightCycler® 96 detection system (Roche Diagnostics (Shanghai) Ltd., Switzerland). The amplification parameters were as follows: an initial denaturation step at 95°C for 3 min, followed by 45 cycles of denaturation at 95°C for 10 s, annealing at 60°C for 30 s, and extension at 72°C for 10 s, with a final extension step at 60°C for 60 s. The melting curve was set according to the instrument’s standard. For data analysis, β-actin was used as the reference gene for normalization, and triplicate measurements were performed. The relative expression levels of the target genes were calculated using the 2^-ΔΔCt method.

Fifty milligrams of feces were placed into a 1.5 mL Eppendorf tube, followed by the addition of 1 mL of pre-chilled ultrapure water and magnetic beads. The mixture was then ground at 30 Hz for 3 min at 4°C. The quality control sample was prepared by mixing 50 μL of each individual sample. The sample, along with the quality control mixture, was then centrifuged at 4°C and 12,000 rpm for 15 min, and the supernatant was collected. The remaining pellet was retreated according to the above steps. The above supernatant was mixed well and then centrifuged at 4°C and 12,000 rpm for 15 min. Finally, 200 μL of the supernatant was taken for further analysis.

For the metabolomic analysis, the method was well established in our laboratory. Briefly, a Dionex Ultimate 3,000 UPLC (Dionex, Germering, Germany) combined with a Q-Orbitrap-MS system equipped with a heated electrospray ionization source (Q-Exactive plus MS, Thermo Scientific, MA) was used. The column was a C₁₈ column (100 mm × 2.1 mm, 1.9 μm), and the column temperature was set at 40°C. The mobile phase solvents, solvent A (0.1% formic acid-water) and solvent B (0.1% formic acid-acetonitrile), were used for sample elution. The column was washed with 100% solvent B and equilibrated with 100% solvent A. The flow rate was 0.2 mL/min, and the injection volume was 5 μL. The MS parameters were as follows: sheath gas flow rate of 10 arbitrary units, spray voltage of 3.0 kV, capillary temperature of 320°C, and auxiliary gas heater temperature of 30°C. A heated electrospray ionization source was used, and the samples were detected in full scan mode. Data for samples in the positive ion mode and negative ion mode were collected in the range of 30–800 m/z.

The raw data was processed using the MSConvert software to convert it to the RAW format. The R software was used for normalization, data transformation, centering, and standardization. Then, the data was imported into the online database XCMS Online Metabolomics (Scripps Research Institute, USA). After exporting the data, a pre-selection was performed based on the criteria of metabolite features >50% and relative standard deviation (RSD) ≤30%. Subsequently, the R software (version 3.5.1) was used for principal component analysis (PCA) to simplify and reduce the dimensionality of the data. The overall distribution between sample groups was analyzed using distance algorithms. The variable importance in the projection (VIP) was used as an important parameter to indicate disruptions in the body’s metabolic network. Candidate variables were selected based on VIP ≥1.0 and statistical significance (adjusted p-value <0.05). Metabolites were compared with the HMDB Library ¹⁰ based on positive and negative ion modes (M + H and M-H).

Human bronchial epithelial (HBE) cells were generously provided by Dr. Haiyan, Chu from the School of Public Health, Nanjing Medical University. HBE cells were cultured in DMEM (Gibco, USA) containing 10% FBS (TransGen Biotech Co., Ltd., Beijing, China) supplemented with 1% penicillin–streptomycin at 37°C with 5% CO₂. Prevotella bivia (ATCC 29303) and modified reinforced clostridial broth (ATCC Medium 2,107) were purchased from Mingzhou Biotechnology Co., LTD. The strain was cultured in Medium 2,107 at 37°C in the Bugbox anaerobic system (Bug Box, Ruskinn Technology, Leeds, UK) for 48 h. The cultures were centrifuged at 5,000 rpm for 5 min at 4°C. The supernatant was filtered to remove bacteria. HBE cells were treated for 24 h either 10% of bacterial supernatant or ATCC Medium 2,107.

HBE cells were washed twice with PBS and resuspended in 100 μL ultrapure water extraction solution (containing protease inhibitors, PMSF, and EDTA). Five hundred microliter of lipid extraction solution and 50 μL cell suspension were mixed by vortexing for 2 min. The mixed samples were sonicated for 5 min, added with 100 μL water, vortexed for 1 min, and then centrifuged at 4°C and 12,000 rpm for 10 min. The 200 μL supernatant was collected and dried with nitrogen. The dried supernatant was re-dissolved in 200 μL of isopropanol/acetonitrile (50%/50%) for further analysis. Metabolic data were acquired by two data acquisition systems, i.e., ultra-performance liquid chromatography (UPLC, ExionLC™ AD) and tandem mass spectrometry (MS/MS, QTRAP® 6,500+). The samples were separated on a Thermo Accucore™ C30 column (2.6 μM, 2.1 mm × 100 mm i.d.) at a flow rate of 0.35 mL/min, and the temperature of 45°C. The solvent system was as follows: A: 0.1% formic acid and 10 mmol/L ammonium formate in acetonitrile/water (60%/40%); B: 0.1% formic acid and 10 mmol/L ammonium formate in acetonitrile/isopropanol (10%/90%). The samples were detected in full scan mode. The characteristic ions of each lipid metabolite were processed by multiple reaction monitoring. Lipid metabolites were identified by MWLDB (Metware lipidomics database v 3.0) provided by MetWare, Co., Ltd. ¹¹ (Wuhan, China). The semiquantitative analysis of the metabolites were calculated based on the peak area ratio of the analyte to the internal standard (Supplementary Table S2).

Wilcoxon’s rank sum tests and student’s t tests were used to determine the significance of differences in alpha diversity, composition, genetic distance and inflammatory factor expression between two groups. The aforementioned tests were carried out by utilizing either R version 3.5.1 or GraphPad Prism version 8.0. A two-sided p-value of less than 0.05 was considered statistically significant.

In this study, the mean (± SD) of age of children with asthma and the control was 6.52 ± 2.60 years and 6.06 ± 2.78 years, respectively; the number of males and females in the group with asthma was 14:16 and 18:12 in the control group. There was no significant difference in the distribution of age and sex between two groups. BMI of all subjects was all in the normal range. All of subjects reported without passive smoking.

The oral microbiota and gut microbiota of both the control and the children with asthma were analyzed through 16S rRNA gene sequencing. Compared to the control group, the Chao, Sobs and Shannon indices of oral microbiota and gut microbiota in the children with asthma were significantly decreased, and the Simpson index increased (Figure 1A). Beta (β) diversity analysis using PCoA (Principal co-ordinates analysis) was used to visualize the differences among groups. Bray-Curtis distance demonstrated statistically significant differences in β diversity of oral microbiota and gut microbiota between two groups (Figure 1B).

The dominant bacterial phyla in oral microbiota were Firmicutes, Bacteroidetes, Proteobacteria, Actinobacteria, and Fusobacteria, while the top eight dominant bacterial genera were Streptococcus, Prevotella, Veillonella, Neisseria, Haemophilus, Actinomyces, Leptotrichia and Granulicatella (Figure 2A); the dominant phyla of gut microbiota were Firmicutes, Bacteroidetes, Actinobacteria, Proteobacteria, and Verrucomicrobiota; the top eight dominant genera were Bacteroides, Blautia, Faecalibacterium, Bifidobacterium, Enterococcus, Allobaculum, Romboutsia, and Clostridium_sensu_stricto_1 (Figure 2B), which were consistent with previous studies (Hu et al., 2022; Zheng et al., 2022).

In oral microbiota, the relative abundance of Proteobacteria was significantly decreased, and although the relative abundance of Firmicutes and Bacteroidetes was increased, there was no statistical significance in the children with asthma. At the genus level (relative abundance >5%), the relative abundance of Prevotella and Veillonella was significantly increased, and Streptococcus showed a decreased trend in the children with asthma (Figure 2C). In gut microbiota, there was an increase in the relative abundance of Firmicutes and Bacteroidetes without significance, and the relative abundance of Proteobacteria was significantly increased and Actinobacterioat was significantly decreased in the children with asthma. At the genus level (relative abundance >5%), the relative abundance of Enterococcus and Allobaculum was significantly increased, while the relative abundance of Blautia, Bifidobacterium, Faecalibacterium, Eubacterium_hallii_group, and Fusicatenibacter was significantly decreased in the children with asthma (Figure 2D).

One hundred and twenty-one shared OTUs coexisted in oral and gut microbiota shown in the Venn diagram (Supplementary Figure S1). The relative abundance of the top 20 shared genera showed the relationship between oral and gut microbiota (high in oral and low in gut microbiota, vice versa) (Figure 3A). Furthermore, we conducted a co-occurrence network analysis on the distribution patterns between samples and species. In the control group, only Streptococcus was present in both the oral cavity and gut, while Prevotella and Veillonella were exclusively found in the oral cavity. In the group with asthma, Prevotella, Veillonella, and Streptococcus were observed in both oral cavity and gut (Figure 3B).

Further analysis of microbiota metagenomic sequencing results revealed a significant increase in the abundance of P. bivia, P. disiens, and P. oris in both oral cavity and gut. Additionally, there was a significant increase in the relative abundance of Bacteroides fragilis in the oral and gut microbiota of children with asthma, while the relative abundance of Streptococcus thermophilus was significantly decreased (Figures 3C,D). To confirm the oral-gut migration of these species in the children asthma, we conducted a strain-level analysis based on 16S rRNA gene using oral and fecal metagenome of the control and fecal metagenome of the group with asthma. Genetic distance of 16S rRNA gene was calculated by p-distance. The strains of Prevotella bivia, P. disiens, P. oris and S. thermophilus in the gut of the group with asthma had a closer relationship to those in the oral cavity of the control, indicating these species originated from the oral microbiota (Figures 4A–E).

Through aligning the amino acid sequences of P. bivia, P. disiens, P. oris, and B. fragilis with the IEDB database, microbial mimicry peptides that exhibited more than 80% similarity to asthma antigenic epitopes were selected (Supplementary Table S3). These peptides were used to stimulate PBMCs, and the expression of relevant inflammatory factors was subsequently evaluated. The results showed that Ppb from P. bivia induced IL-4, IL-5, IL-6, and IL-8 expression; Rpd from P. disiens induced IL-6 and IL-8 expression. Fbf was derived from B. fragilis, and induced IL-4, IL-5, IL-6, IL-8, and IL-17A expression (Figure 4F).

Based on gut metagenomic data, we identified KEGG pathways with an LDA (Linear discriminant analysis) score greater than 2 by LEfSe analysis (Supplementary Figure S2). The asthma group was characterized by amino acid metabolism and lipid metabolism. Spearman correlation analysis was performed between the above 5 different species and selected pathways. B. fragilis, P. bivia, P. disiens and P. oris had associations with these two types of metabolic pathways (Figure 5A).

In order to elucidate the functional changes in gut microbiota and their interaction with the host, we conducted a comprehensive analysis of the metabolome in fecal samples using an untargeted ultrahigh-performance liquid chromatography-mass spectrometry (UPLC-MS) approach. The score plots generated from PCA depicted distinct separation between the control and the group with asthma (Figure 5B). The differential metabolism were determined based on a variable importance in projection (VIP) score greater than 1 and a false discovery rate lower than 0.05. The certain metabolism such as malic acid, SDA, ALA, and EPA (fatty acyl); DG and MG (glycerolipids); LPC, PC, and PA (glycerophospholipids) were decreased, while L-carnitine (fatty acyl) and Cer (sphingolipids) were increased compared to the control (Figure 5C; Supplementary Table S4).

Next, a redundancy analysis (RDA) using metagenomic and metabolomic data revealed a complex co-occurrence pattern and correlation between bacteria associated with childhood asthma and characterized metabolites from the essential metabolic pathways. In the gut of the children with asthma, P. bivia, P. disiens, and P. oris showed a negative correlation with glycerophospholipid (LPC, PC, and PA), and B. fragilis showed a positive correlation with L-carnitine and Cer. As a probiotic, S. thermophilus showed an opposite trend compared to Prevotella and Bacteroides (Figure 5D).

Human bronchial epithelial (HBE) cells were cultured with the supernatant of P. bivia and were detected with lipid metabolome based on UPLC-MS/MS. PCA performed on all the samples revealed that a clear clustering of QC (quality control) samples, and lipid components in the control group and treatment group were well separated (Supplementary Figure S3). A total of 1,221 lipid metabolites were detected in the samples, categorized into 45 lipid subclasses, and the types and quantities of lipid subclasses detected were presented in Supplementary Figure S4A. The distribution of various lipid subclasses between the two groups was shown in Supplementary Figure S4B. Based on volcano plot analysis, 220 differential lipid metabolites were selected with VIP > 1 and adjusted p-value <0.05. Among them, 55 showed a significant increase, while 165 showed a significant decrease in the treatment group (Figure 6A). Through visualizing the distribution of differential lipids in different lipid subclasses (Supplementary Figure S5) and filtering the top 20 based on fold change, in the treatment group, it revealed that a significant decrease in LPC, LPG, PA, PC, PG, PI, and PS, and a significant increase in CE, Cer, and TG (Figures 6B,C). These results were partially consistent with findings from previous non-targeted metabolomics analysis. Most of lipid metabolites (87.72%) were enriched in metabolic pathways (Supplementary Figure S6). Based on the differential lipid metabolites, KEGG pathway enrichment analysis revealed pathway disruptions in the treated group, primarily enriched in glycerophospholipid metabolism (p < 0.05) (Figure 6D).