Bacterial Microbiota-derived Extracellular Vesicles in Children With Allergic Airway Diseases: Compositional and Functional Features.

PURPOSE: Bacterial extracellular vesicles (EVs) play crucial roles in bacteria-host interactions. Due to their cargo, EVs are considered fingerprints of the parent cell, which are detectable in body fluids. We studied the composition and function of bacterial microbiota-derived EVs genes in urine to evaluate whether they have specific characteristics concerning allergic airway disease.
METHODS: Subjects were from elementary school surveys and classified into 3 groups according to questionnaires and sensitization to aeroallergens: the allergic airway group (AA, n = 16), atopic controls (AC, n = 7) and healthy controls (HC, n = 26). The bacterial EVs were isolated from voided urine samples, their nucleic acid was extracted for 16S ribosomal RNA pyrosequencing and then characterized using α-diversity, β-diversity, network analysis, intergroup comparison of bacterial composition and predicted functions, and correlation with total immunoglobulin E (IgE), eosinophils% and fractional exhaled NO.
RESULTS: The compositional α-diversity was the highest in AA, while functional α-diversity was the highest in HC. AA had a distinct clustering with the least intersample variation. Klebsiella, Haemophilus, members from Lachnospiraceae and Ruminococcaceae, and the pathways of sphingolipid and glycerolipid metabolism, and biosynthesis of peptidoglycan and lysine were the highest in AA and positively correlated with total IgE or eosinophil%. Genetic information processing function contributed to 48% of the intergroup variance and was the highest in AA. Diaphorobacter, Acinetobacter, and the pathways of short-chain fatty acids and anti-oxidants metabolism, lysine and xenobiotic degradation, and lipopolysaccharide biosynthesis were the lowest in AA and negatively correlated with total IgE or eosinophil%. The bacterial composition and function in AC were closer to those in HC. The bacterial network was remarkably dense in HC.
CONCLUSIONS: The bacterial microbiota-derived EVs in urine possess characteristic features in allergic airway disease with a remarkable correlation with total IgE and eosinophil%. These findings suggest that they may play important roles in allergic airway diseases.

INTRODUCTION

Bacterial extracellular vesicles (EVs) are enriched with bioactive proteins, lipids, nucleic acids and virulence factors; thus, they can act on a variety of signaling pathways and distantly transport their cargo.1 Despite being a significant player in bacteria-bacteria and bacteria-host interactions, only a few studies have addressed the characterization of bacterial microbiota-derived EVs (MEVs) in health and disease.23

Due to its structure, cargo and high stability, EVs can be considered the fingerprint of the parent cell that can be easily detected in body fluids.4 Those characteristics make the EVs an attractive candidate for biomarker studies.5 Human cell-derived EVs in urine have been studied as a biomarker for multiple diseases5; moreover, proteomic analysis of urine EVs has shown a broad representation of EVs from across the body without bias towards the kidney or urine.6 Furthermore, when bacterial EVs were injected intraperitoneally, they spread to the whole mouse body and accumulated in the tissues of the liver, lung, spleen and kidney.7

Allergic airway diseases (AAD) such as asthma and rhinitis are heterogeneous and the implementation of precision medicine in the management of AAD requires the identification of phenotype-specific markers measurable in easily accessible biological fluids, such as urine.

Accordingly, we hypothesized that bacterial EVs in urine might reflect the microbiota consortium dysbiosis in AAD, adding a new piece to the puzzle of the relationship between microbiota and AAD. Therefore, we analyzed the bacterial composition and predicted function of MEV genes in the urine of children suffering from allergic rhinitis with or without concomitant asthma.

MATERIALS AND METHODS

Subject enrolment and sample collection

Subjects were enrolled from a population-based study as part of the Seongnam Atopy Project for Children's Happiness (SAP2017) program.8 They were screened by using the Korean version of the International Study of Asthma and Allergies in Childhood questionnaire9; skin prick test to Dermatophagoides pteronyssinus, Dermatophagoides farinae, cat fur, dog fur, Alternaria, birch tree, orchard grass and Japanese hop8; serum allergen-specific immunoglobulin E (IgE) by ImmunoCAP system (Phadia AB, Uppsala, Sweden) to Dermatophagoides farinae, dog dander, cat dander, birch, Japanese hop, and Alternaria8; body mass index (BMI); and urine analysis.

Rhinitis was considered present if a patient had current rhinitis symptoms and rhinitis treatment. Asthma was considered present if a patient had current asthma symptoms and asthma treatment. Controls were those with neither allergic disease diagnosis ever or treatment ever nor chronic rhinitis symptoms. Sensitization to aeroallergens was diagnosed by a positive skin prick test (wheel diameter > 3 mm) or positive specific IgE (> 0.35 kU/L). The inclusion criteria were an average BMI, normal urine analysis, and the absence of current respiratory illness.

Forty-nine subjects were divided into 3 groups: the allergic airway group (AA, n = 16), the atopic controls group (AC, n = 7) and the healthy control group (HC, n = 26). All the AA subjects had rhinitis, 6 of whom had accompanying asthma with documented sensitization to at least 1 aeroallergen; the rhinitis symptoms were assessed by using rhinitis severity visual analog scale from 0 to 10. The AC subjects were controls sensitized to at least 1 aeroallergen. The HC subjects were controls with no sensitization to the aeroallergens tested.

Subjects further underwent testing for fractional exhaled nitric oxide (FeNO), impulse oscillometry (IOS), acoustic rhinometry, a complete blood analysis and total IgE assessment as described before.9 Fresh voided clean catch mid-stream urine was collected, and 2 mL per sample were stored at −20°C until subjected to metagenomic analysis within 1 month after sampling.

The study protocol was approved by the Institutional Review Board of CHA University (2017-04-049). Written informed consents were obtained from the participant students and their guardians.

Nucleic acid extraction

The urine samples underwent differential centrifugation by microcentrifuge (Labogene 1730R, BMS, Seoul, Korea) at 10,000 g for 10 minutes at 4°C.10 The supernatant was filtered through a 0.22-µm filter to eliminate bacteria and foreign particles. The separated EVs were boiled at 100°C for 40 minutes and centrifuged at 13,000 rpm for 30 minutes at 4°C to eliminate the remaining floating particles and waste, and then the supernatant was taken for DNA extraction. Nucleic acid was extracted by using a DNeasy PowerSoil Kit (QIAGEN, Hilden, Germany) and then quantitated using QIAxpert system (QIAGEN) following the standard protocol in the kit guide (Supplementary Data S1).

16S ribosomal RNA (16S rRNA) sequencing

The gene-specific sequences targeted the 16S V3 and V4 hypervariable region.11 The libraries were prepared using polymerase chain reaction products according to the MiSeq System guide (Illumina, San Diego, CA, USA) and quantified using a QIAxpert (QIAGEN). The sequence reads were clustered into operational taxonomic units (OTUs) using the SILVA 128 database at 97% similarity. In cases the taxon could not be assigned, we assigned it at the higher level and indicated that in parenthesis. The functions of the MEVs genes were predicted using the phylogenetic investigation of communities by reconstruction of unobserved states (PICRUST).12 The OTUs could not be assigned as bacteria or predicted as bacterial functions were excluded before statistical analysis (Supplementary Data S1).

Diversity and clustering

The Shannon diversity index was used for the estimation of α-diversity. A phylogenetic tree was built with FastTree using R version 3.5.3 to compute the Faith phylogenetic diversity. The Bray Curtis dissimilarity index was applied for the analysis of similarity (ANOSIM), principal coordinates analysis (PCoA) and hierarchical using PAST V3 to assess intersample variation (β-diversity).

To identify the variables influencing the samples' clustering, we tested cllustering for association with the age, sex, BMI z-score, family history of allergic diseases, mode of delivery, passive smoking and clinical phenotype using a linear regression model; the sample HC1 was used as the reference point.

Statistical analysis

Kruskal-Wallis 1-way analysis of variance and Wilcoxon rank-sum tests were used for intergroup comparisons. Nonparametric Spearman's rank test was used to test the correlation of bacterial composition and predicted functions with total IgE, eosinophil% and FeNO; a correlation coefficient (r) of > 0.4 or < −0.4 were considered significant. All the tests were 2-sided, and the statistical significance was corrected for the Benjamini-Hochberg false discovery rate (FDR) at a q value < 0.05. We have further assessed the absolute counts of bacterial taxa and functional KEGG orthology groups (KOs) by the fold change, and log2 more than 1 or less than −1 was considered significant. The statistical testing was done with the MatLab 2018 or SPSS 22.

We considered the taxa account for >0.5% relative abundance as core taxa and took them into consideration. Cytoscape version 3.4.0 was used for network visualization of the significant co-occurrences and anti-occurrences of the core genera in the AA and HC. Due to small differences in the mean relative abundances of level 1 functions, we used the similarity of percentage (SIMPER) to assess their contribution to the overall variance between the study groups using PAST V3. We used the iPath tool13 for the projection of the predicted metabolic functions.

RESULTS

Characteristics of the study subjects

The AA subjects had significantly higher total IgE, blood eosinophil count and percentage, and FeNO compared to the HC. The AA subjects'rhinitis symptoms ranged from mild to moderate (3.3–7.8), and their acoustic rhinometry showed average nasal patency. The 6 subjects with concomitant asthma were well-controlled and had average IOS measurements. The AC subjects had significantly higher total IgE, and blood eosinophil count and percentage compared to the HC subjects. The subjects' characteristics are described in Table.

Table

Characteristics of the study subjects

Variables	AA (n = 16)	AC (n = 7)	HC (n = 26)	P value
				AA-AC	AA-HC	AC-HC
Age, year	11 (11, 12)	10 (10, 12)	10.9 (10, 11.2)	NS	NS	NS
Male, sex (%)	12 (75)	6 (86)	10 (38)	NS	0.023	0.035
BMI, kg/m2	19.5 (17, 21.6)	17.8 (14.8, 20.1)	19.4 (17.3, 22.5)	NS	NS	NS
BMI z-score	0.18 (−0.70, 0.83)	−0.14 (−1.7, 0.19)	0.2 (−0.5, 1)	NS	NS	NS
Blood neutrophil count, cells/µL	45.7 (42.8, 48.6)	47.3 (42.8, 50.2)	46.7 (42.6, 53.2)	NS	NS	NS
Blood eosinophil count, cells/µL	224.5 (148.9, 376.9)	211.2 (122.8, 518.9)	118.2 (80.8, 147.8)	NS	0.001	0.010
Eosinophil%	4.1 (3.1, 7.8)	4 (2, 7)	2 (1.3, 2.2)	NS	0.003	< 0.001
FeNO, ppb	19.5 (13.8, 28.8)	17 (15, 18)	13.5 (10.8, 16.3)	NS	0.003	NS
Acoustic rhinometry*, cm3	1.7 (0.2, 2.6)	0.23 (−0.55, 2.1)	1.9 (1.15, 3.7)	NS	NS	NS
Serum total IgE, kU/L	127 (84, 379.8)	99.1 (36.7, 148)	43.3 (42.6, 53.2)	NS	< 0.001	0.043
Number of positive sIgE†	1–3	1	0	-	-	-
Number of positive SPT‡	1–5	1–2	0	-	-	-
Subjects with positive pollen sensitization	6	3	0	-	-	-
Subjects with positive dust mites sensitization	14	3	0	-	-	-
Subjects with positive cat or dog fur sensitization	7	2	0	-	-	-

Values are presented as median (interquartile range) or number (%).

Statistical significance is determined by using the following tests: Kruskal-Wallis, χ2, or Wilcoxon rank-sum test for pairwise comparison.

AA, allergic airway group; AC, atopic controls group; HC, healthy controls group; NS, not significant; FeNO, fraction of exhaled nitric oxide; BMI, body mass index; IgE, immunoglobulin E; SPT, skin prick test; sIgE, serum immunoglobulin E.

*Acoustic rhinometry: the difference between basal and after nasal decongestion; †Number of positive specific IgE measurements (sIgE > 0.35 kU/L); ‡Number of positive SPT (wheel diameter > 3 mm).

The 16S-rRNA sequencing resulted in 3,462,100 reads clustered into 252,065 OTUs. Bacterial composition Shannon diversity index and Faith phylogenetic diversity were higher in the AA subjects compared to the HC subjects (P < 0.001) and positively correlated with total IgE (P < 0.001, r = 0.45 and 0.50, respectively) and with eosinophil% (P < 0.001, r = 0.41 and 0.46, respectively) (Fig. 1A and B). In contrast, functional diversity was the lowest in the AA subjects compared to the HC and AC subjects (P < 0.001 and P = 0.042, respectively) and negatively correlated with total IgE (P < 0.001 and r = −0.60) (Fig. 1C). In both compositional and functional α-diversity, the AC subjects were between the AA and HC subjects matching its phenotype (Fig. 1).

Fig. 1

Alpha-diversity. (A) Bacterial composition Shannon diversity index and (B) bacterial composition PD were the highest in the AA and differed significantly between the AA and HC. (C) The predicted functions Shannon diversity index was the lowest in the AA and differed significantly from that in the HC and AC.

AA, allergic airway group; AC, atopic controls group; HC, healthy controls group; PD, phylogenetic diversity.

A Wilcoxon rank-sum test was performed.

The AA samples had the least inter-sample variation compared to the HC and AC samples (bacterial composition: P < 0.001 and 0.026, respectively; predicted functions: P < 0.001 and 0.041, respectively) (Supplementary Fig. S1), which explains their distinct clustering in the PCoA space and the hierarchy (Fig. 2). The HC sample hierarchy showed a better clustering in the functional analysis (> 85% similarity) compared to the bacterial composition (< 40% similarity). The AC samples merged the AA and HC samples in the PCoA space and scattered in-between the AA and HC samples without specific clustering in the hierarchy (Fig. 2). The AC samples' bacterial compositions and predicted functions did not differ from those of the HC samples in terms of the Bray Curtice dissimilarity index (P = 0.562 and 0.500, respectively) (Supplementary Fig. S1). This was further confirmed by ANOSIM analysis at 9,999 permutations (corrected P = 1 and 0.480, respectively). It is noteworthy that the Bray Curtis dissimilarity index of the bacterial compositions and the predicted functions had positive correlations with total IgE (P < 0.001 and r = 0.58 and 0.60, respectively).

Fig. 2

PCoA and hierarchy. (A) Bacterial composition and (B) predicted functions, samples as represented by the first 2 principal coordinated based on the Bray Curtis dissimilarity index, where each point represents a single sample. They show the compositional and functional dissimilarity between the AA and HC samples. The AA samples are condensed in a relatively narrow area, while the AC overlapped both AA and HC. (C) Bacterial composition and (D) predicted functions, hierarchical of the study samples; the AA samples arranged adjacently under the same cluster, and the HC samples showed better clustering in the functional profiling. The AC samples were scattered between the AA and HC samples.

AA, allergic airway group; AC, atopic controls group; HC, healthy controls group; PCoA, principal coordinates analysis.

Regarding the variables influencing the clustering, the only factor that had a significant effect on the clustering pattern was the clinical phenotype (P < 0.001); meanwhile, the effects of age, sex, BMI z-score, family history of allergic diseases, mode of delivery and passive smoking were insignificant.

Bacterial composition

Firmicutes dominated the microbiome community in the AA subjects, while Proteobacteria dominated it in the HC subjects. In spite of the wide variability in Firmicutes and Proteobacteria across the study groups, we noticed that their sum was nearly constant (Supplementary Fig. S2A); therefore, we calculated the Firmicutes/ Proteobacteria ratio and compared it among the study groups. The Firmicutes/Proteobacteria ratio was the highest in the AA subjects (5.6) and the lowest in the HC subjects (0.8) (P < 0.001) (Supplementary Fig. S2B), and positively correlated with total IgE (P < 0.001 and r = 0.57).

After initially selecting the core taxa according to the relative abundance, the significant genera were defined according to a significant correlation with total IgE or eosinophil% (q < 0.05 and r >0.40 or <−0.40), the intergroup difference in relative abundance (q < 0.05), and absolute counts (log2 fold change > 1 or < −1) (Fig. 3). The significant genera higher in the HC subjects were members from the phyla Bacteroidetes (Cloacibacterium) and Proteobacteria (Dechloromonas, Comamonas, Diaphorobacter, Comamonadaceae (f), Acinetobacter and Enterobacter). The significant genera enriched in the AA subjects were members from the phyla Actinobacteria (Collinsella), Bacteroidetes (Prevotella 9 and Alistipes), Proteobacteria (Escherichia-Shigella, Klebsiella, and Haemophilus) and Firmicutes (Christensenellaceae R-7 group, 4 Lachnospiraceae, 5 Ruminococcaceae, Catenibacterium, and Dialister). It is remarkable that, despite characteristic depletion in Proteobacteria in the AA subjects, the 3 genera known to be pathogenic (Escherichia-Shigella, Klebsiella and Haemophilus) were significantly higher in the AA subjects and positively correlated with total IgE and eosinophil%. Neither significant enrichment nor depletion was detected in the AC subjects compared to the AA and HC subjects at an FDR q-value < 0.05; however, some differences were detected in comparison to AA subjects at q-value < 0.1, but not to the HC subjects. Prevotella 9, Haemophilus, Acinetobacter, Ruminococcus torques and Lachnospiraceae UCG-008 were significantly depleted in the AC subjects compared to the AA subjects. We also compared subjects with allergic rhinitis alone and those with allergic rhinitis and concomitant asthma, and no statistically significant difference was detected (data are not shown).

Fig. 3

Significant bacterial genera; mean relative abundance, OTU count log2 fold change, and correlation with eosinophil% and total IgE.

Italic font indicates the genera belong to Lachnospiraceae and Ruminococcaceae families. A-, B-, F-, P- denotes the phyla (A-, Actinobacteria; B-, Bacteroidetes; F-, Firmicutes; P-, Proteobacteria).

AA, allergic airway group; AC, atopic controls group; HC, healthy controls group; NS, not significant; IgE, immunoglobulin E.

*Significant difference between the AA and HC groups by using Wilcoxon rank-sum test (q < 0.05); †Significant difference between the AA and AC groups by using the Wilcoxon rank-sum test (q < 0.1); ‡Absolute counts in AA divided by those in HC log2; §Significant correlation by using Spearman's rank correlation coefficient (q < 0.05).

Massive significant bacterial interrelationships were detected in the HC subjects. All of the core genera, except Staphylococcus, were interrelated in the HC subjects with 530 co-occurrences or anti-occurrences, giving a complex interaction network unlike that of the AA subjects, which included only 14 genera with 11 correlations (Fig. 4).

Fig. 4

Bacterial network analysis. Each edge represents a significant correlation colored to indicate either positivity (red) or negativity (blue). Nodes colored to indicate the enriched side; brown if higher in HC, green if higher in AA. Node size is proportionate to the count of significant correlations, either positive or negative; the largest is Bacteroides (48 correlation). Node border indicates the presence or absence of a correlation with total IgE and eosinophil%; yellow if correlated and white if not. Correlations were tested by Spearman's rank correlation coefficient with a correlation cut-off r-value of greater than 0.4 or less than −0.4, false discovery rate q < 0.05. (A) Bacterial network in the AA: 14 genera had a significant interrelationship, with 11 either co-occurrences or anti-occurrences. 13 of those genera are from the genera higher in AA and one is from the genera higher in the HC and had an anti-occurrence relationship. Only 5 of the 14 are from the genera correlating with total IgE and eosinophil%. (B) Bacterial network in the HC: 48 out of the 49 core genera were interrelated with 530 co-occurrences or anti-occurrences relationships, giving a complex interaction network unlike that of the AA. The genera that were found higher in AA had a co-occurrence relationship with each other and an anti-occurrence with the genera that were higher in the HC, which also had a co-occurrence relationship with each other except the Enterococcus. Enterococcus was higher in the HC but with no statistical significance.

AA, allergic airway group; HC, healthy controls group; IgE, immunoglobulin E.

Predicted functions

The genes of level 1 predicted functions: metabolism, genetic information processing, environmental information processing and cellular processes differed significantly among the 3 groups with significant correlations with total IgE and eosinophil%. Genetic information processing was the superior contributor in the intergroup variance (48%), followed by metabolism (21%) (Fig. 5).

Fig. 5

Level 1 functions relative abundance, SIMPER analysis, and correlation with total IgE. (A) The level 1 functions differ significantly between the 3 groups (Wilcoxon rank-sum test, q < 0.05); metabolism is the most abundant function. (B) Genetic information processing is the major contributor to the overall variance among the study groups. (C) All level 1 functions have a moderately significant correlation with total IgE.

AA, allergic airway group; AC, atopic controls group; HC, healthy controls group; SIMPER, similarity of percentage; IgE, immunoglobulin E.

*Wilcoxon rank-sum test (q < 0.05); †Spearman's rank correlation coefficient.

The abundance of many level 3 metabolic pathway genes differed significantly between the AA and HC subjects (q < 0.05), a few differed between the AA and AC subjects, and none differed between the HC and AC subjects. The relative abundances and 95% confidence interval of the relevant metabolic pathways are represented in Fig. 6, and the significant KOs fold change (log2 > 1 or < −1) are shown in Supplementary Fig. S3. It is essential to illustrate that the metabolic pathway comprises both anabolism and catabolism, and it was not feasible to elucidate which direction overshoots the other.

Fig. 6

Level 3 metabolic pathways relevant to AAD: The data are shown as mean relative abundance and a 95% confidence interval.

AA, allergic airway group; AC, atopic controls group; HC, healthy controls group; AAD, allergic airway diseases.

*Functions differed significantly between the AA and the other 2 groups (Wilcoxon rank-sum test, q < 0.05); †Functions differed between the AA and HC group (Wilcoxon rank-sum test, q < 0.05). No statistically significant difference detected between the AC and HC groups.

The metabolic pathway genes enriched in the HC subjects comprised the short-chain fatty acids (SCFAs) propanoate and butanoate metabolism, retinol metabolism, ascorbate and aldarate metabolism, glutathione metabolism, ubiquinone and other terpenoid-quinone biosynthesis, LPS biosynthesis, flagellar assembly, and lysine degradation (q < 0.05). These functions had negative correlations with total IgE (q < 0.05 and r > 0.45), and most of their significant KOs were higher in the HC subjects. The lipoic acid metabolism was also higher in the HC subjects, but none of its KOs was significant. The genes of cytochrome P450 xenobiotics biodegradation were extremely depleted in the AA subjects (q < 0.05), negatively correlated with total IgE (q < 0.05 and r > 0.45), and most of its KOs were higher in the HC subjects; nevertheless, the steroid degradation pathway could not be detected in the MEV genes in any of the study subjects.

The genes of the sphingolipid (SL) metabolism, glycerolipid metabolism, lysine biosynthesis and peptidoglycan biosynthesis pathways were the highest in the AA subjects, positively correlated with total IgE (q < 0.05 and r > 0.45), and most of their significant KOs were higher in the AA subjects. The SL metabolism also significantly correlated with eosinophil% (q < 0.05 and r = 0.55) and did not differ between the AA and AC subjects. The oxidative phosphorylation pathway did not differ among the study groups, but most of its significant KOs were higher in the AA subjects (Supplementary Fig. S3).

The projection of the metabolic pathways detected spread along most of the global metabolic pathway map (Supplementary Fig. S4) where the denser projections were in carbohydrate metabolism, energy metabolism, lipid metabolism and amino acid metabolism pathways, and the less dense projections were in glycan biosynthesis and metabolism as well as in xenobiotics biodegradation and metabolism pathways.

DISCUSSION

Our results demonstrated significant differences between the AA and HC subjects in diversity metrics, bacterial composition, bacterial interrelationships and predicted functions with evident correlations between total IgE and eosinophil%, supporting that MEVs may play roles in the implication of microbiota in AAD.

The AA subjects in this study had the highest compositional Shannon diversity index and Faith phylogenetic diversity, consistent with results of previous studies on airway microbiota in adults with atopic and eosinophilic asthma.1415 In contrast, functional α-diversity metrics were the highest in the HC subjects, which may indicate a higher capability for competition in regard to invasion or overgrowth of harmful bacteria. We observed significant correlations between α-diversity metrics and total IgE or eosinophil%, which suggests that the atopic background is an important factor for the differences among the study groups.

The homogeneity of the AA samples in both compositional and functional analysis was striking. They showed distinct clustering in the PCoA space and the hierarchy with significant differences compared with the HC and AC samples, which may indicate the ability of the urine MEVs to characterize AAD. The HC samples clustered at a higher similarity in the functional analysis compared to the compositional analysis, suggesting that the functional analysis may be superior to the compositional analysis in defining symbiosis.

The Firmicutes/Bacteroidetes ratio has commonly been used to report gut microbiome symbiosis and dysbiosis16; meanwhile, our study samples showed a characteristic trend in the Firmicutes/Proteobacteria ratio as it was the highest in the HC subjects and the lowest in the AA subjects. This can be simply attributed to the differences in the sampling site since the MEVs in urine might express a broader portion of the microbiome community. In addition, we suggest that enrichment with the phylum Firmicutes at the expense of the rest of the community might be a sign of dysbiosis.

Members from phylum Proteobacteria abundant in our HC subjects were reported to play protective roles against atopy. Experimental exposure to Acinetobacter species protected against allergic sensitization and lung Inflammation via induction of strong Th1 and anti-inflammatory responses by immune cells and skin cells,17 and protected against parental asthma through Toll-like receptor (TLR) signaling.18 In a relatively large study on environmental and skin microbiota, Gammaproteobacteria and its genus Acinetobacter were significantly depleted in atopic subjects' skin19; furthermore, in healthy but not in atopic individuals, Gammaproteobacteria was associated with interleukin (IL)-10 expression, and Acinetobacter species was associated with a robust balance in anti-inflammatory response and Th1/Th2 gene expression.1719 We found 3 members from the Comamonadaceae family (Diaphorobacter, Comamonas testosterone and undefined Comamonadaceae genus) that were significantly higher in the HC subjects. Diaphorobacter which represents 79% of the Comamonadaceae family played a role in the T cells balance and then the response to allergens via direct links with TLR2, which is a potent inducer of IL-10 that further bridges into a network of anti-inflammatory transcripts including forkhead box P3, transforming growth factor-β, TLR2, and aryl hydrocarbon receptor gene.17 However, Comamonas testosteroni has been shown to possess steroid-responsive degradation pathways; hence, it was suggested to play a role in steroid-resistant asthma, especially after the Comamonadaceae family in bronchial brushings were reported to correlate with bronchial hyperresponsiveness in suboptimally controlled asthmatics.20 Nevertheless, we did not detect any steroid degradation activity in any of our samples; in addition, the cytochrome P450 xenobiotics biodegradation pathways were enriched in the HC subjects, not in the AA subjects.

Despite the enrichment with Proteobacteria in the HC subjects, the commonly known pathogens (namely the genus Haemophilus the family Enterobacteriaceae, and its genera Klebsiella, and Escherichia-Shigella) were enriched in the AA subjects. Haemophilus has been frequently reported in association with neutrophilic and atopic asthma in adults and a higher risk of asthma in children,141521 which could be attributed to its ability to induce Th2 and Th17 responses.21 Enterobacteriaceae was reported in the bronchial brushings of patients with severe asthma.22 Klebsiella, especially Klebsiella pneumonia, was found to be abundant in patients with severe asthma compared to those with moderated asthma and controls.22 In this study, K. pneumonia formed 90% of the total Klebsiella abundance and was the highest in the AA subjects.

Perturbations of the bacterial composition in the current AA subjects were mainly represented as enrichment with members from the Firmicutes phylum. Half of the significant genera higher in our AA subjects belonged to the Lachnospiraceae and Ruminococcaceae families, which were reported to be associated with atopy in multiple studies. Enrichment with Ruminococcus and Subdoligranulum in the gut were linked to a higher risk of food sensitization, and the latter was associated with a lower butyrate concentration and SCFA in feces.2324 Blautia and Ruminococcus were enriched in the gut of children at high risk of asthma compared to those with average risk.25 Analysis of gut microbiota in infant twins associated the development of allergies, especially respiratory allergies, with a higher abundance of Lachnospiraceae and Ruminococcus gnavus.26 R. gnavus-fed mice developed airway inflammation characterized by expansion of Th2 cells in the colon and lung with infiltration of colon and lung parenchyma by eosinophils and mast cells.26 Taken together, the enrichment with those families could be an indicator of allergy. The family Christensenellaceae and its genus Christensenellaceae R-7 were higher in the AA subjects. Christensenellaceae was suggested as evidence for the effect of the human genome on microbiome community composition after it was found to be the most heritable taxon in twin pairs.27 We have tested Christensenellaceae and Christensenellaceae R-7 for the maternal history of allergic diseases, as the infant microbiota is mainly obtained from the mother and they were significantly associated with a positive maternal history of allergic disease (P = 0.038 and 0.031, respectively). These results may raise a possible role for the microbiome in disease inheritance.

We observed enrichment with a member of Bacteroidetes (Prevotella 9) in the AA subjects compared to the HC and AC subjects. There is increasing evidence linking elevated Prevotella to chronic inflammatory diseases through Th17 mediated mucosal inflation28; in addition, it was reported to be higher in children with sensitization to food compared to controls,23 which is consistent with our results. However, the Prevotella abundance in bronchial brushings was increased in atopic non-asthmatics compared to atopic asthmatics and controls in adults.14 It is necessary to clarify whether this discrepancy is related to differences in age groups.

We used network analysis to visualize the potential cooperation and competition between the microbiota in our study groups. The wide gap in the complexity of the bacterial associations in the HC subjects compared to that in the AA subjects reflects the cohesiveness of the microbiome community under healthy conditions compared to disease status. This may propose a broader definition for dysbiosis to include not only disturbed bacterial composition but also a disrupted bacteria-bacteria relationship. As most of the genera in the AA network are those without a correlation with eosinophil% nor total IgE, we proposed that the association with allergy status may have substituted or dispersed the inter-bacterial relationship, and the bidirectional relationship between host eosinophil and microbiome reinforced our proposal.29

There have been few studies describing predicted functions in AAD. This study reported significant results in metabolism, genetic information processing, environmental information processing and cellular processes pathways showing the ability of MEVs in urine to express the functional derangement in AAD. Enrichment with SCFAs metabolic pathway genes among the HC subjects was consistent with their protective role against lung inflammation in mice and their association with a lower risk of asthma in children.3031 Furthermore, multiple studies have demonstrated their anti-inflammatory properties by modifying the T cells function via immunologic and epigenetic mechanisms.3233

Relative enrichment among our HC subjects with LPS biosynthesis genes is consistent with a previous report on the lung microbiome.14 Exposure to a bacterial LPS-rich environment ‘farm dust’ induced tolerance in lung epithelial cells via suppression of Th2 activity, and the administration of LPS to germ-free animals was sufficient to restore oral tolerance in mice.3435 Therefore, we suggest that higher LPS biosynthesis activity of the human microbiome could potentiate the propensity to tolerance. Peptidoglycan and flagella are thought to participate in the gut lung axis by the same mechanism as LPS; however, neither was reported in prior AAD microbiome studies.

The abundance of SL and glycerolipid metabolism genes were significantly higher in the AA subjects and were positively correlated with total IgE. Both are bioactive lipids involved in cell signaling leading to bronchoconstriction and inflammation.3637 Also, altered SL metabolism with an augmented synthesis of sphingosine-1-phosphate was demonstrated in response to allergen challenge in young adults with allergic rhinitis with or without asthma and was suggested to participate in the development of the asthma phenotype in house dust mite allergic patients.38 Moreover, increases in glycerolipid metabolites were reported in global plasma metabolomic profiling of asthmatics.39 Bacterial cell-membrane glycosphingolipids are recognized by iNKT cells and can activate their induction of IL-4 and IL-13.40 A small group of microbes can produce SL, including Bacteroides and Sphingomonas.41 In this study, the relative abundance of Sphingomonas was extremely low (< 0.5%), unlike Bacteroides, which ranged from 1% to 11% and was significantly the highest in the AA (q < 0.05). Consistently, higher abundances of gene function related to the SL metabolic pathways have been reported in the fecal microbiome of infants with bronchiolitis compared to healthy controls,42 and we noticed that the taxonomic profiling of the same subject had a Bacteroides-dominant profile.43 Bacteroides has been considered a healthy bacteria due to its ability to produce SCFAs16; however, its ability to produce SL41 has made its abundance double-edged. This demonstrates the necessity of functional analysis in explaining the relative enrichment or depletion of bacteria.

Lysine metabolism showed a characteristic pattern in this study since the AA subjects were depleted with lysine degradation and enriched with lysine biosynthesis genes. This finding is consistent with 2 whole metagenomic sequencing studies of the upper airway microbiome in children, adolescents, and the elderly that showed depleted lysine degradation in the asthma group.44 Lysine modifications are a critical step in collagen post-translational modifications,45 which has a pivotal role in airway remodeling in asthma.46 Importantly, although the pulmonary functions of our AA subjects were within the normal range, the functional analysis was able to detect this difference.

We detected enrichment with genes of the anti-oxidants metabolic pathways, glutathione, lipoic acid, retinol, and ascorbate and aldarate metabolism, and the ubiquinone and other terpenoid-quinone biosynthesis pathways in the HC subjects. In contrast, the abundance of oxidative phosphorylation pathway genes did not differ between the AA and HC subjects, with most of the oxidative phosphorylation enzymes were higher in the AA subjects. These findings suggested a favorable reactive oxygen compounds/anti-oxidants balance in the microbiome community in the HC subjects. The reactive oxygen compounds released by bronchial epithelial cells, eosinophils, and neutrophils are believed to be a part of the inflammatory response detected in AAD47; in addition, reactive oxygen compounds signaling is an essential mediator of the immune system and inflammatory processes that modulate the gut epithelial barrier.48

The AA subjects in this study were significantly enriched in genetic information processing function genes at the 3 hierarchical levels. Information processing functions were mostly overlooked in previous AAD microbiome studies. However, exploring the MEVs brought attention to those functions. The bacterial EVs play roles in the bacteria-bacterial interaction, including cell to cell signaling and horizontal gene transfer.50 Genes shared within a specific site may likely reflect adaption to a particular habitat or niche; hence, intense exchange of the genetic processing genes may enhance replication and biofilm formation of the abundant taxa inducing and sustaining dysbiosis.

There is a growing orientation towards machine learning for microbiome data that aims at a more in-depth understanding microbiome community paving the way for feasible diagnostic and therapeutic applications. A bioinformatics team has studied microbiome data from daily stool samples using machine learning aiming at the prediction of dysbiosis.49 They concluded that the microbiome functions and metabolic byproducts are superior to taxonomy in the prediction of dysbiosis.49 Looking through the results of the functional analysis studies, including ours, we noticed that they were more consistent and well-framed compared to those of the taxonomy; therefore, we propose that the microbiome function may surpass the microbiome taxonomy in the biomarkers' field.

The values of the bacterial composition and predicted functions measured in the AC were always between the AA and HC subjects. However, in terms of statistical significance, the AC differed from the AA subjects, but not the HC subjects. This might be interpreted as the features of MEV genes in urine in atopy without allergic symptoms are not yet shifted from normal, or the small number of the AC subjects (n = 7) is not enough to make a statistical difference from the HC subjects.

The limitations of our study are the small number of the AC subjects and undetermined allergy medications. However, strengths include a strict selection of the study subjects from the same age group and average BMI to minimize the effect of confounders. The AA subjects MEVs were well-characterized with significant correlations with total IgE and eosinophil% in spite of being a population-based study using a non-invasive technique. We suggest that concomitant sampling from the gut, airway and urine will better demonstrate the ability of MEVs in urine to reflect the microbiome consortium and reveals its biomarker potentials.

In conclusion, the MEVs' diversity and clustering differed between the AA and HC subjects. Bacteria repeatedly detected in relation to AAD and metabolic functions with evident relevance to immune signaling and airway remodeling were significantly higher in the AA; while, bacteria and metabolic functions reported to play protective roles against atopy were significantly higher in the HC subjects. We found an association between the features of MEVs and levels of total IgE and eosinophil%. The characteristics of MEVs in the AC subjects were between the AA and HC subjects with a close similarity to the HC subjects. The bacterial interrelationships were markedly diminished in the AA subjects. The genetic information processing genes were notably enriched in the AA subjects, which may be a cause of dysbiosis. Our results demonstrate the importance of MEVs in the understanding of the role of bacterial microbiota in AAD.