Microbiome-Transcriptome Interactions Related to Severity of Respiratory Syncytial Virus Infection

Abstract

Respiratory syncytial virus (RSV) is a major cause of lower respiratory tract infections and hospital visits during infancy and childhood. Although risk factors for RSV infection have been identified, the role of microbial species in the respiratory tract is only partially known. We aimed to understand the impact of interactions between the nasal microbiome and host transcriptome on the severity and clinical outcomes of RSV infection. We used 16 S rRNA sequencing to characterize the nasal microbiome of infants with RSV infection. We used RNA sequencing to interrogate the transcriptome of CD4+ T cells obtained from the same set of infants. After dimension reduction through principal component (PC) analysis, we performed an integrative analysis to identify significant co-variation between microbial clade and gene expression PCs. We then employed LIONESS (Linear Interpolation to Obtain Network Estimates for Single Samples) to estimate the clade-gene association patterns for each infant. Our network-based integrative analysis identified several clade-gene associations significantly related to the severity of RSV infection. The microbial taxa with the highest loadings in the implicated clade PCs included Moraxella, Corynebacterium, Streptococcus, Haemophilus influenzae, and Staphylococcus. Interestingly, many of the genes with the highest loadings in the implicated gene PCs are encoded in mitochondrial DNA, while others are involved in the host immune response. This study on microbiome-transcriptome interactions provides insights into how the host immune system mounts a response against RSV and specific infectious agents in nasal microbiota.

Conflict of interest statement

The authors declare no competing interests.

Figures

Figure 1

An overview of the data and analyses performed in this study.

Figure 2

(A) Overview of study timelines, including when each type of data was collected. All 58 infants included in this study had both transcriptomic and microbiome data collected for at least one time point. (B) Venn diagrams showing the overlap in the infants that had either transcriptomic or microbiome data collected during the acute visit (57 of the 58 infants), the overlap in the transcriptomic data collected at either the initial (Visit 1) or follow-up visit (Visit 3), and the overlap in the infants microbiome data collected at either the initial (Visit 1) or follow-up visit (Visit 2). Numbers are shown separately for mild and severe groups, with the total number in each category indicated below the Venn diagram. (C) Differential gene expression analysis comparing Visit 1 versus Visit 2 RNA-Seq samples yields 27 significantly differentially expressed genes (FDR < 0.05). The log2 of the expression levels of these genes are shown as boxplots combining mild and severe samples. The genes on the left side of the panel have higher mean expression during the acute visit and those on the right side of the panel have higher expression levels in the follow-up visit (D,E) The log2 expression levels for EZH1 comparing (D) data from infants with severe versus mild RSV infection and (E) data from the acute and follow-up visit. The reported FDR p-values are based on limma analysis. (F,G) Abundance profiles of nasal microbiota comparing groups of mild and severe infants at each of the two visits, including (F) observed diversity and (G) Shannon diversity index.

Figure 3

(A) Principal component analysis was performed on Visit 1 samples for which we had both transcriptome and microbiome data to reduce the transcriptome into gPCs (top) and the microbiome into cPCs (bottom). (B) Cumulative distribution of the amount of variance explained by the top gPCs and cPCs. The vertical and horizontal lines indicate the number of PCs which explain 95% of variance and which were included in our network analysis. (C) The coordinates of the gene and clade principal components across the 40 analyzed samples from Visit 1. (D) Spearman correlation across the loadings of the top 13 gPCs and the top 10 cPCs.

Figure 4

(A) The weights of edges predicted for each of the 40 sample-specific networks obtained from the LIONESS analysis. Columns are grouped based on each infant’s infection severity, sex, and race. Rows are sorted based on the significance of edges. (B) The significance of each edge as defined from multivariate linear analysis comparing edge-weights between two groups, mild and severe, and corrected for sex and race. The top significant edges with nominal p-value less than 0.05 are labeled.

Figure 5

GO term enrichment analysis for the top 100 genes based on the loadings of (A) gPC1, (B) gPC8, and (C) gPC3. These bubble plots include the 25 most significant GO terms identified in each analysis based on the FDR significance. All terms shown are significant at an FDR < 0.05. Size of each bubble indicates the number of genes annotated to the respective GO term and the color indicates the percentage of the top 100 genes annotated to that term. The top 25 GO terms enriched in the other identified gPCs are shown in Supplemental Fig. 6A,B and all GO terms enriched at an FDR < 0.05 are shown in Supplemental Fig. 6C–G.