Integrated nasopharyngeal airway metagenome and asthma genetic risk endotyping of severe bronchiolitis in infancy and risk of childhood asthma

Abstract

Background: Infants with bronchiolitis are at increased risk of developing asthma. Growing evidence suggests bronchiolitis is a heterogeneous condition. However, little is known about its biologically distinct subgroups based on the integrated metagenome and asthma genetic risk signature and their longitudinal relationships with asthma development.

Methods: In a multicentre prospective cohort study of infants with severe bronchiolitis (i.e. bronchiolitis requiring hospitalisation), we profiled nasopharyngeal airway metagenome and virus at hospitalisation, and calculated the polygenic risk score of asthma. Using similarity network fusion clustering approach, we identified integrated metagenome-asthma genetic risk endotypes. In addition, we examined their longitudinal association with the risk of developing asthma by the age of 6 years.

Results: Out of 450 infants with bronchiolitis (median age 3 months), we identified five distinct endotypes, characterised by their nasopharyngeal metagenome, virus and asthma genetic risk profiles. Compared with endotype A infants (who clinically resembled "classic" bronchiolitis), endotype E infants (characterised by a high abundance of Haemophilus influenzae, high proportion of rhinovirus (RV)-A and RV-C infections and high asthma genetic risk) had a significantly higher risk of developing asthma (16.7% versus 35.9%; adjusted OR 2.24, 95% CI 1.02-4.97; p=0.046). The pathway analysis showed that endotype E had enriched microbial pathways (e.g. glycolysis, l-lysine, arginine metabolism) and host pathways (e.g. interferons, interleukin-6/Janus kinase/signal transducers and activators of transcription-3, fatty acids, major histocompatibility complex and immunoglobin-related) (false discovery rate (FDR)<0.05). Additionally, endotype E had a significantly higher proportion of neutrophils (FDR<0.05).

Conclusion: In this multicentre prospective cohort study of infant bronchiolitis, the clustering analysis of integrated-omics data identified biologically distinct endotypes with differential risks of developing asthma.

Figure 1.. Study design and analytic workflow

The analytical cohort consists of 450 infants with severe bronchiolitis in a multi-center prospective cohort study—the 35th Multicenter Airway Research Collaboration (MARC-35). First, we computed a distance matrix for each of four datasets: metagenome taxonomy (species-level) and functional data, virus (including binary variable of RSV, RV, and RV-C, and genomic load of RSV and RV) data, and asthma PRS data. Then, we computed an affinity matrix of each dataset separately, and generated a fused affinity matrix by similarity network fusion before clustering analysis. To identify integrated endotypes of severe bronchiolitis, we used spectral clustering. Second, we characterized the clinical features and asthma genetic risk of the integrated endotypes. We subsequently investigated the association of the endotypes and the risk of developing childhood asthma and its subtypes (i.e., atopic asthma and non-atopic asthma). Lastly, to examine the distinct function of each endotype, we performed functional class scoring analysis with the hallmark gene sets and GO terms. We also investigated the association of the endotypes with seven blood immune cells. Abbreviations: GO, gene ontology; IgE, immunoglobulin E; PRS, polygenic risk score; RSV, respiratory syncytial virus; RV, rhinovirus.

Figure 2.. Taxonomy profile, clinical feature, and asthma genetic risk of integrated endotypes

A. The t-distributed stochastic neighbor embedding (tSNE) of the integrated endotypes. We used the five eigenvectors in the spectral clustering in the tSNE calculation. Each dot represents a single infant with bronchiolitis, who cluster together according to their integrated nasopharyngeal metagenome, virus, and asthma genetic risk profiles. B. The relative abundance of the five most abundant nasopharyngeal microbial species among infants with bronchiolitis. The box plots show the distribution (median and interquartile range) of the bacterial species according to integrated endotypes. C. The distribution of integrated endotypes with major clinical and virus features. Ribbons connect individual integrated endotypes to the major clinical and virus characteristics. For the age variable (continuous), the ribbon width represents mean age (i.e., the narrower the ribbon width, the younger infants mean age in the corresponding endotype). For other variables (binary), ribbon width represents the proportion of infants within the endotype who have the corresponding clinical or virus features, which was scaled to a total of 100%. D. The distribution of the integrated endotypes based on their asthma genetic risk as measured by asthma PRS. Asthma PRSs were categorized into deciles. Deciles 1-2 denote low genetic risk group, deciles 3-8 denote medium genetic risk group, deciles 9-10 denote high genetic risk group.

Figure 3.. Association of integrated endotypes with risk of developing asthma and its subtypes

The forest plot shows the analysis for the association of the integrated endotypes (endotype A as the reference) with the risk of developing asthma, atopic asthma, and non-atopic asthma. The multivariable logistic regression models adjusted for potential confounders, including age, sex, race/ethnicity, parent history of asthma, number of previous breathing problems.

Figure 4.. Association of integrated endotypes with microbial functional pathways

A. Association of the integrated endotypes with microbial metabolic (MetaCyc) pathways. The heatmap shows the top 25 most significant pathways based on five group overall comparison by using the Kruskall-Wallis test (all FDR<0.05, green color bar on the left side). The asterisks show the significance of the association between endotypes (endotype A as the reference) by using Wilcoxon rank-sum test. B. Association of the integrated endotypes with microbial KEGG pathways. The heatmap shows the top 25 most significant pathways based on five group overall comparison by using Kruskal-Wallis test (all FDR<0.05, green color bar on the left side). The asterisks show the significance of the association between endotypes (endotype A as the reference) by using Wilcoxon rank-sum test.

Figure 5.. Association of integrated endotypes with nasopharyngeal airway host transcriptome pathways and blood immune cells

A. Functional class scoring analysis using Molecular Signatures Database(MSigDB) hallmark genesets. We identified 11 differentially enriched host pathways (FDR <0.05) with the highest absolute value of normalized enriched score for the endotype A (reference) versus E comparison. The MSigDB hallmark collection contains 50 genesets, which summarize and represent specific well-defined biological states or processes and display coherent expression. B. Functional class scoring analysis using Gene Ontology (GO) terms. We identified seven differentially enriched host pathways (FDR <0.05) associated with endotype E. C. Association of the integrated endotypes with seven host blood immune cells. The y-axis denotes proportion of the blood immune cells. Abbreviations: FDR, false discovery rate; GO, gene ontology; NES, normalized enrichment score.