Nasopharyngeal Microbiota, Host Transcriptome, and Disease Severity in Children with Respiratory Syncytial Virus Infection

Abstract

Rationale: Respiratory syncytial virus (RSV) is the leading cause of acute lower respiratory tract infections and hospitalizations in infants worldwide. Known risk factors, however, incompletely explain the variability of RSV disease severity, especially among healthy children. We postulate that the severity of RSV infection is influenced by modulation of the host immune response by the local bacterial ecosystem.

Objectives: To assess whether specific nasopharyngeal microbiota (clusters) are associated with distinct host transcriptome profiles and disease severity in children less than 2 years of age with RSV infection.

Methods: We characterized the nasopharyngeal microbiota profiles of young children with mild and severe RSV disease and healthy children by 16S-rRNA sequencing. In parallel, using multivariable models, we analyzed whole-blood transcriptome profiles to study the relationship between microbial community composition, the RSV-induced host transcriptional response, and clinical disease severity.

Measurements and main results: We identified five nasopharyngeal microbiota clusters characterized by enrichment of either Haemophilus influenzae, Streptococcus, Corynebacterium, Moraxella, or Staphylococcus aureus. RSV infection and RSV hospitalization were positively associated with H. influenzae and Streptococcus and negatively associated with S. aureus abundance, independent of age. Children with RSV showed overexpression of IFN-related genes, independent of the microbiota cluster. In addition, transcriptome profiles of children with RSV infection and H. influenzae- and Streptococcus-dominated microbiota were characterized by greater overexpression of genes linked to Toll-like receptor and by neutrophil and macrophage activation and signaling.

Conclusions: Our data suggest that interactions between RSV and nasopharyngeal microbiota might modulate the host immune response, potentially affecting clinical disease severity.

Keywords: disease severity; microbiota; nasopharynx; respiratory syncytial virus; transcriptome profiling.

Figure 1.

Nasopharyngeal microbiome composition in young children with respiratory syncytial virus (RSV) infection and healthy control (HC) subjects and associations with host characteristics. (A) Dendrogram visualizing an average linkage hierarchical clustering of individuals on the basis of the Bray–Curtis dissimilarity matrix. The branches of the tree structure were colored according to health status (HC subjects, green; patients with RSV, yellow). Information on age and disease severity is depicted adjacent to the branch ends. Stacked bar charts show the relative abundance of the 15 highest-ranked operational taxonomic units (OTUs) and of residual bacteria specified by phylum. OTUs are color coded according to phylum: Firmicutes, red; Proteobacteria, blue; Actinobacteria, yellow; and Bacteroidetes, green. On the basis of clustering indices, an optimal number of 11 clusters was identified, 5 of which comprised more than three study participant samples. Classifier taxa for these five clusters (color-coded horizontal panels) were: Streptococcus (STR), light red; Haemophilus influenzae (HPH), dark blue; Corynebacterium (COR), yellow; Moraxella (MOR), light blue; and Staphylococcus aureus (STA), dark red. Gray panels mark individuals not included in any of these five clusters. (B) A nonmetric multidimensional scaling (NMDS) plot was used to visualize the associations between nasopharyngeal microbiota clusters (see color coding in A) and host characteristics: age group (round arrowheads), health status (arrowheads, bold text), and antibiotic treatment (arrowheads, bold italic text), depicting both the individual nasopharyngeal microbiota composition (data points [n = 132] and ellipses [SD of data points] colored according to cluster) and the 10 highest-ranked OTUs. We observed that RSV hospitalization was related to STR- and HPH-enriched profiles. Additionally, we observed an association between outpatients, older age, and MOR enrichment, whereas STA-dominated profiles were observed primarily in young, healthy infants. The orthogonal orientation (∼90°angle) of the vectors associated with health status and age imply that these host characteristics are not highly correlated, which was verified by permutational multivariate analysis of variance (unadjusted R² severity, 5.7%; age, 6.7%; and interaction severity–age, 5.9%). AB = antibiotic treatment; Inp = inpatients; Outp = outpatients.

Figure 2.

Independent associations between disease severity, age, and the relative abundance of classifier taxa. Multivariable linear regression modeling was used to assess the independent association between disease severity, defined as the need for hospitalization, and age (independent variables), and the relative abundance of classifier taxa (Haemophilus influenzae, Streptococcus, Corynebacterium, Moraxella, and Staphylococcus aureus; outcome variable). Adjusted effect sizes (aESs) are visualized by a forest plot (outpatients, gray squares; inpatients, white squares; age, black squares) and 95% confidence intervals (CIs) (bars). Benjamini-Hochberg–corrected P values (q values) were calculated. Background colors indicate the classifier species (Streptococcus, light red; H. influenzae, dark blue; Corynebacterium, yellow; Moraxella, light blue; and S. aureus, dark red). We observed that higher abundance of H. influenzae and Streptococcus was positively associated (i.e., aES > 0) with hospitalization, whereas S. aureus enrichment was negatively associated (i.e., aES < 0) with being an inpatient and was observed more often in young infants.

Figure 3.

Transcriptional modular repertoire analysis stratified by nasopharyngeal microbiome cluster. (A) Heat map depicting the difference between the proportion of over- and underexpressed genes (modular expression; Figure E3) per cluster compared with healthy control (HC) subjects on a per-module basis. The color gradient indicates the proportion of differentially expressed genes per module (modular expression) ranging from −100% (blue) to 100% (red). Modular expression between −10% and 10% is shown in white and reflects no differences compared with HC subjects. Although IFN-related modules (M1.2, M3.4, M5.12) were similarly overexpressed in all clusters, we observed increased expression of inflammation (M3.2, M4.2, M4.6, M4.13, M5.1, M5.7), neutrophil (M.5.15), and monocyte/macrophage modules (M4.14) in the Haemophilus influenzae (HPH) and Streptococcus (STR) clusters. Contrarily, the shared adaptive immune response to respiratory syncytial virus was underexpressed and distinguished by particularly strong suppression of B cell– (M4.10) and T/natural killer (NK) cell–related modules (M3.6, M3.1, M4.15) in the STR cluster. (B) Box plots showing the per-cluster median fold-change gene expression relative to the HC baseline. A representative selection of modules is shown. Median fold-change expression values of inflammation (M3.2 and M4.6) and neutrophil modules (M5.15) in the HPH-enriched cluster were significantly higher than that in the Moraxella (and Staphylococcus aureus) clusters. Statistically significant differences in median fold change between clusters was assessed using a Kruskal–Wallis test (P values shown) with Nemenyi post hoc test (bars; *P < 0.05). COR = Corynebacterium; MOR = Moraxella; STA = Staphylococcus aureus.

Figure 4.

Identification and functional annotations of respiratory syncytial virus (RSV) and microbiota-specific differentially expressed genes between patients with RSV and healthy control (HC) subjects. (A) Heat maps depict differentially expressed genes between HC subjects and each microbiota cluster (limma, stringent filter: P < 0.01, log₂-fold change > 1.25, Benjamini–Hochberg multiple test correction), adjusted for age and sex. Normalized expression is indicated as overexpressed (red) or underexpressed (blue) compared with the median expression of HC subjects (yellow). Heat maps were scaled according to number of samples and transcripts per comparison. The highest number of differentially expressed transcripts was observed in the Streptococcus (STR; 1,435 genes) and Haemophilus influenzae (HPH) clusters (1,315 genes) compared with the other microbiota clusters. (B) Venn diagram showing the intersection between differentially expressed genes derived from pairwise comparisons between patients with RSV stratified by microbiome profiles and HC subjects (A). We observed a considerable overlap between the STR and HPH clusters (464 genes). In contrast, the Corynebacterium (COR, 323 genes), Moraxella (MOR, 651 genes), and Staphylococcus aureus clusters (STA, 207 genes) showed a small number of differentially expressed genes versus HC subjects. Using ingenuity pathway analysis, shared and cluster-specific genes were extracted and functionally annotated: (C) genes for each of the clusters were intersected and those shared between all or all but one clusters (209 genes; common RSV signature); (D) genes shared between or unique to the STR and HPH clusters (1,257 genes; shared STR/HPH signature); and (E) genes shared between COR, MOR, and STA, but not the STR and HPH clusters (258 genes; shared COR/MOR/STA signature), were extracted. For each gene list of interest, the 10 strongest associated over-/underexpressed pathways were visualized; Benjamini-Hochberg–corrected P values (q values) associated with a specific pathway were calculated using Fisher’s exact test and visualized on a log₁₀ scale. (C) We hypothesized that the genes shared between all or four out of five clusters would represent the common RSV signature, as this feature is shared between all clusters. Indeed, we observed these genes were strongly associated with IFN signaling (q value = 8.1 × 10⁻⁸). (D) Genes shared between or unique to the HPH and STR clusters were involved (q value ranging from 6.5 × 10⁻⁴ to 5.9 × 10⁻⁶) in the activation of neutrophils (N-formylmethionyl-leucyl-phenylalanine [fMLP], triggering receptor expressed on myeloid cells [TREM]-1, IL-8, and IL-17A signaling) and macrophages (inducible nitric oxide synthase [iNOS] signaling and production of nitric oxide and reactive oxygen species [ROS]). In addition, genes related to the pattern recognition of bacteria and viruses were up-regulated, including Toll-like receptor (TLR) genes TLR4, TLR6, and TLR8. (E) Genes shared between COR, MOR, and STA, but not the STR and HPH clusters, were not significantly associated with any of the 10 highest-ranking canonical pathways (q = 0.48). EIF = eukaryotic initiation factor; mTOR = mechanistic target of rapamycin; PRRs = pattern recognition receptors; trx = transcripts.

Figure 5.

Nasopharyngeal microbiome composition in infants with respiratory syncytial virus (RSV) infection and associations with host characteristics and gene principal components (gPCs). (A) A nonmetric multidimensional scaling (NMDS) plot was used to jointly visualize individual microbiota profiles (data points [n = 74] and ellipses [SD of data points] are colored according to cluster), the 10 most-abundant nasopharyngeal microbiota, host characteristics (age and hospitalization status), and gPCs extracted by principal component analysis from the host transcriptome dataset. (B) Multivariable linear regression models were fitted to investigate the association between classifier taxa (explanatory variables) and gPCs (outcome variables) within the RSV cohort (n = 74), corrected for age and multiple testing (Benjamini–Hochberg correction). We used a forest plot to visualize adjusted effect size (aES) (colored squares) and 95% confidence intervals (CIs) (bars). Only gPCs with significant results (defined as q < 0.15, depicted in bold) are shown. We verified the positive associations between Streptococcus (STR) and Haemophilus influenzae (HPH) and gPC1 and between Staphylococcus aureus (STA) enrichment and gPC6 and gPC7, and the negative association between Corynebacterium (COR) abundance and gPC3 observed using NMDS (A). Color coding: STR, light red; HPH, dark blue; COR, yellow; Moraxella, light blue; and STA, dark red. (C) Heat map depicts the log₂-fold change gene expression of gPC1 genes identified using DAVID pathway enrichment analyses for each individual, stratified by health status/microbiota cluster. Pathways were selected based on a cluster enrichment score of greater than 1 and a Benjamini-Hochberg–corrected P value < 0.05. Genes covered by more than one probe are shown as [gene abbreviation].1 or [gene abbreviation].2 (e.g., MAPK14.1). Normalized expression is indicated as overexpressed (red) or underexpressed (blue) compared with the median expression of healthy control subjects (yellow). Explanations of the abbreviations of the depicted genes are given in Table E5. HC = healthy controls; Inp = inpatients; MOR = Moraxella; Outp = outpatients; TLR = Toll-like receptor.