key: cord-0004058-zgurb8gb
authors: Sonawane, Abhijeet R.; Tian, Liang; Chu, Chin-Yi; Qiu, Xing; Wang, Lu; Holden-Wiltse, Jeanne; Grier, Alex; Gill, Steven R.; Caserta, Mary T.; Falsey, Ann R.; Topham, David J.; Walsh, Edward E.; Mariani, Thomas J.; Weiss, Scott T.; Silverman, Edwin K.; Glass, Kimberly; Liu, Yang-Yu
title: Microbiome-Transcriptome Interactions Related to Severity of Respiratory Syncytial Virus Infection
date: 2019-09-25
journal: Sci Rep
DOI: 10.1038/s41598-019-50217-w
sha: 74af665d4b129e0ed662f50fa0d215f4ff2c2c52
doc_id: 4058
cord_uid: zgurb8gb

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.

In addition to major risk factors such as congenital or chronic cardiopulmonary disease, meta-analyses have identified additional risk factors for RSV, including preterm birth, low birth weight, siblings at home, day-care attendance, and maternal smoking [6] [7] [8] . Infants with RSV infection in early childhood are also at a higher risk of developing asthma and wheezing [9] [10] [11] . Although prophylactic treatments are available to prevent RSV infection in the most at-risk infants, attempts are ongoing to identify an effective and safe vaccine or small molecule drug to reduce the health burden of RSV.

Biomarker discovery through expression analysis is an important step in assessing disease severity and in distinguishing RSV from other common respiratory viruses [12] [13] [14] . Transcriptional profiles often reflect a host's immune response to the virus, helping to explain disease progression and characterize its severity. Several studies have reported that T cell mediated response is crucial in clearing the viral load during RSV infection 15, 16 . Recent reports also indicate that RSV infects CD4 + and CD8 + T cells and affects T cell function 12, 17, 18 , implicating T cells as potential biomarkers for RSV severity.

The microbiome in the respiratory tract is known to influence the course of acute infectious diseases 19 . The succession pattern of the nasal microbial community might influence host responses to RSV, thereby modulating inflammation and possibly disease severity. Indeed, several studies indicate that the nasal microbial composition affects the overall risk of developing respiratory tract infections 20, 21 and is associated with the severity of acute respiratory symptoms 22 .

Joint characterization of the nasal microbiome and host transcriptome may provide valuable insights into how viral infections influence the host response. However, the impact of the interactions between the nasal microbiome and the host transcriptome on the severity and clinical outcomes of RSV infection has not been fully understood. Recent studies have reported associations between nasal microbial compositions and whole blood gene expression 22 . In this work, we perform an integrative analysis to study associations between nasal microbial compositions and transcriptional profiles of bloodstream CD4 + T cells.

Network approaches provide an important framework for understanding complex relationships that influence human health 23 . In our study, we construct a network model that correlates microbial taxonomic profiles and host transcriptomic profiles. Examining this network in the context of RSV severity highlights important patterns of transcriptomic response related to immune processes and viral infection. An overview of our analysis approach is summarized in Fig. 1 .

Baseline characteristics of participants and data source. We analyzed data for 58 infants with RSV infection, each of whom had both CD4 + gene expression and nasal microbiome data. The RSV disease severity was measured by the Global Respiratory Severity Score (GRSS). This allowed us to divide infants into two groups based on severity: mild and severe (see Methods section). We emphasize that GRSS is a composite score that reflects the worst values during the entire illness, rather than a single time point, allowing us to address the problem of clinical variability among subjects. GRSS also recapitulates more frequently used factors such as viral load (which has been associated with RSV disease severity in several studies 24, 25 ) and need for hospitalization (Supplemental Fig. 1 ), thus acting as excellent surrogate for disease severity 26 . Among the infants, 23 had mild illness and 35 had severe illness. There was no difference in any demographic characteristics between mildly and severely ill subjects ( Table 1) . The data collection occurred in three stages ( Fig. 2A) . transcriptome data analysis. We first characterized the gene expression data for these infants, as measured by RNA-sequencing (RNA-Seq) of CD4 + T cells collected at two distinct time points, with 46 samples collected during the acute visit (Visit 1) and 34 samples collected at follow-up day 12-16 (Visit 2). Out of 22 subjects who had RNA-Seq samples collected at both the acute and follow-up visit, 7 were mild and 15 were severe (Fig. 2B) .

We used principal component analysis (PCA) to investigate the data structure in light of various measured clinical variables (Supplemental Fig. 2 ). This analysis indicated that enrollment season was a significant contributor to structure in our expression data. Therefore, we applied batch-correction to remove this signal from the data. We then performed differential expression analysis, comparing the expression levels of genes between the acute and follow-up visit. We identified 27 genes that were differentially expressed at a Benjamini-Hotchberg false discovery rate (FDR) less than 0.05 (Fig. 2C) . Many of the genes that had higher expression levels in the acute visit are known to be important in mediating host immune response. This included IFITM1 and IFITM2, which are known to inhibit the infection and replication of respiratory syncytial virus 27,28 , IFI27, a known biomarker for RSV 29 , and IRF7, a gene associated with suppression of innate immunity response 30 . Other genes with known associations to RSV and higher expression in the acute visit included SOCS3 31 , MX1 32 , and the ISGylation pathway genes USP18 and ISG15 33 . Our analysis also identified several genes with higher expression in the follow-up visit, including MS4A1, which encodes the CD20 protein, and MCOLN1. CD20 + B cells are prominent in the lung tissue of infants with fatal RSV infection 34, 35 and MS4A1 is upregulated in infants after the administration of live attenuated influenza vaccine, indicating an association with immune system processes 36 . RSV activates innate immunity through the toll like receptor (TLR) pathway. MCOLN1 has been associated with TLR signaling through modulating viral pathogen-associated molecular pattern (PAMP) along with trafficking of single-stranded RNA (ssRNA) into lysosomes 37, 38 . MCOLN1 also regulates autophagy 39 , as expected in the convalescence phase of an infection.

We also evaluated if any genes changed their expression levels in response to RSV infection severity. We found 600 genes differentially expressed at a nominal p-value significance of less than 0.05. Out of these, the most significantly differentially-expressed gene was EZH1 (p-value = 2.3 × 10 −6 ; FDR = 0.026) (Fig. 2D,E) . EZH1 was previously included in a biosignature proposed as a molecular diagnosis tool for RSV infection 40 . Interestingly, both EZH1 and IFITM have been implicated in immune response signaling indicating resolving infection 27 .

Because we were only able to identify a handful of genes, primarily related to immune response, as significantly differentially-expressed in these same infants, we also investigated whether differences between severe and mild RSV infection could be captured using other methods. In particular, we examined whether gene expression differences between mild and severe RSV infection might be better captured through changes in variability rather than in mean expression levels. We evaluated the differential variability in gene expression levels between mild versus severe samples using the F-test and identified 641 (5.5%) genes that were differentially variable at an FDR significance threshold of 0.05. Gene Ontology enrichment analysis 41 identified a number of biological processes nominally associated with these genes (Supplemental Fig. 3 ). Many of these were associated with mitochondrial activity such as mitochondrial gene expression and translation. Mitochondria are important for viral suppression of the innate immunity 30 . www.nature.com/scientificreports www.nature.com/scientificreports/ Microbiome data analysis. We next characterized microbiome data from nasal swabs obtained by 16 S rRNA sequencing. This included data collected at two time points, with 51 samples collected during the acute visit (Visit 1) and 40 collected during a one-month follow-up visit (Visit 3); 23 of these samples were from infants with mild RSV infection and 35 were from infants with severe RSV infection (Fig. 2B) .

We used both PCA and Principal Coordinate Analysis (PCoA) with various dissimilarity metrics, including rooted-Jensen-Shannon divergence (rJSD) 42 as well as weighted and unweighted Unifrac distance 43 , to reduce dimensionality and visualize the microbiome data (Supplemental Fig. 4) 44 . We also applied MaAsLin (Multivariate association Analysis with Linear modeling) to test for significant relationships between microbial taxa and clinical outcome, after adjusting for sex, race, and enrollment season 45 . At a nominal p-value of 0.05, this analysis identified one operational taxonomic unit (OTU) -H. influenza -as positively associated with severity, and two OTUs -Ralstonia and Streptococcus -as negatively associated with severity. We also assessed the influence of visit and infection severity on the microbiome composition by using two measures to quantify the α-diversity: (1) observed number of operational taxonomic units (OTUs), i.e., the OTU richness and (2) the Shannon index 46 . We found that both α-diversity measures have higher values for infants with severe infection compared to those with mild infection during the acute visit (Visit 1), as well as higher values in the acute as compared to the post-acute visit (Visit 3) for infants with severe infection (Fig. 2F ,G). In particular, the Shannon diversity is significantly different between severe versus mild samples in Visit 1 (p-value = 0.011), and between the Visit 1 and Visit 3 samples among infants with severe RSV infection (p-value = 0.0057). The infection score was estimated during Visit 1 (the acute phase). Severity of infection generally decreases over time, potentially explaining why we observe that the α-diversity decreases between the Visit 1 and Visit 3 samples. integrative analysis of transcriptome and microbiome data. Although it has been suggested that host gene expression is influenced by the microbiome, the biological mechanisms that may facilitate these types of interactions are largely unknown 47 . Integrative analysis of microbiome and transcriptome data could help us understand the relationships between host gene expression, microbial composition, and disease pathogenesis. Therefore, we applied a network-based approach, which systematically evaluates the interdependence of multiple biological entities instead of looking at each one independently, to identify microbiome-transcriptomic relationships important for RSV severity.

We focused on the 40 infants from Visit 1 that had paired transcriptome and microbiome data (see Fig. 2B ) and performed dimension reduction on these data using PCA (Fig. 3A ). This analysis identified 13 gene principal components (gPCs) and 10 clade principal components (cPCs) that explained 95% of the variance in each of their respective data (Fig. 3B) .

The loadings of each cPC and gPC represent a pattern of highly correlated microbial and transcriptional abundances, respectively (Fig. 3C ). In order to relate these patterns and identify associations between highly varying genes and clades, we calculated the Spearman correlation between the loadings of the top gPCs and cPCs (Fig. 3D ). We find several gPC-cPC pairs that are highly correlated, such as gPC5 and cPC4 (most positively correlated; ρ = 0.48), and gPC4 and cPC7 (most negatively correlated; ρ = −0.43). Interestingly, these relationships are not limited to the top few PCs, which are normally the primary focus of dimension reduction analysis. Instead, these relationships highlight that prominent patterns in microbiome data may be associated with more subtle patterns in gene expression, and vice versa.

integration of microbiome/transcriptome relationships with clinical characteristics. To interpret these results, we considered this correlation matrix as edges in a bipartite graph, where nodes are gPCs and cPCs. This network framework can help us identify important transcriptomic-microbiomic (gPC-cPC) relationships. However, since this network was derived using information from all samples, by itself it is unable to shed light on which of these relationships might be associated with differences in the various phenotypic or clinical properties of the input samples, including disease severity.

To overcome this limitation, we applied LIONESS 48, 49 , which employs a jackknife approach to reverse engineer a set of sample-specific networks, to our gPC/cPC correlation network. This allowed us to construct separate correlation networks between gPCs and cPCs for each of the infants and to analyze gPC-cPC relationships in light of clinical information for these infants.

The weights of edges across these networks are shown in Fig. 4A . We compared the distribution of the sample-specific edge-weights between the mild and severe groups and identified six edges which were nominally www.nature.com/scientificreports www.nature.com/scientificreports/ significant with a p-value less than 0.05: gPC1 with cPC2, gPC8 with cPC8 and cPC4, gPC3 with cPC10, gPC13 with cPC9, and gPC12 with cPC1 (Fig. 4B) .

To assess the robustness of this result, we performed a sensitivity analysis by randomly selecting 10 samples each from mild and severe groups and repeating this analysis 1000 times. This analysis allows us to assess the robustness of our results given our relatively low sample number as well as the disparity between the number of mild and severe samples. If an edge is robustly different between the severe and mild groups, we will recover that association repeatedly across the randomizations. We found that the top edges obtained above were frequently (at least 20% of the time compared to 5% expected by chance) identified as significant with the same direction of effect (Supplemental Fig. 5 ).

Microbiome/transcriptome relationships identify clades that may impact the host immune response. Six clade principal components, cPC2, cPC8, cPC10, cPC9, cPC1, and cPC4 (ordered by significance) were identified as differentially-associated with gPCs in the context of RSV severity based on our network analysis. To understand why these cPCs might have a different relationship with the host transcriptome in the context of RSV severity, we identified their associated top OTUs (Table 2) . We find that genus like Streptococcus, Corynebacterium, Alloiococcus, Haemophilus influenzae, and Staphylococcus are among the top OTUs in cPC2, Ralstonia, Corynebacterium, Neisseriaceae (family, genus not known), and Pseudomonas are among the top OTUs in cPC8, and Corynebacterium, Ralstonia, Alloiococcus, and Staphylococcus are among the top OTUs in cPC10. Some of these same OTUs were also identified as OTUs that discriminate between severe and mild RSV infection in our MaAsLin analysis.

The gene principal components that were identified as differentially-associated with cPCs in the context of RSV severity included gPC1, gPC8, gPC3, gPC13, and gPC12. We identified the top genes from these gPCs based on their loadings (Table 3) . Interestingly, we found that 10 of top 20 genes in gPC1 are mitochondrial genes. We also found IFITM1 or IFITM2 in the top loadings of gPC3. These two genes are also significantly differentially-expressed between the acute and post-acute visit, but not between infants with severe versus mild RSV infection (see Fig. 2C ). CCR7 (Chemokine receptor type 7), a top gene in gPC1, gPC8, and gPC3, and ILR7 (interleukin 7 receptor) from gPC8 and gPC12 have been found to be downregulated in RSV patients 50 . Other genes present in the loadings of multiple gene PCs included L-ribosomal proteins (RPL family), which are involved in various pathophysiological process, and SELL (L-Selectin), a cell surface lectin mostly expressed in leukocytes, which is down-regulated in RSV infection. SELL also plays a key role in the recruitment of neutrophils to roll along the endothelium to the infected tissue 51, 52 .

We performed Gene Ontology (GO) enrichment analysis on the top 100 genes associated with each of these gPCs. For gPC1, top significant terms included "immune system process", "immune response", "viral process", Continued "T cell activation", and "leukocyte activation". Terms associated with gPC8 included "viral gene expression", "viral transcription", and "multi-organism process". Terms associated with gPC3 were associated with signaling pathways like "cytokine-mediated signaling", "type I interferon", "innate immune response", "defense response", "response to other organism", "viral genome replication", "defense response", and "regulation of viral genome replication" (Fig. 5) . Similar functions were identified for gPC12 and gPC13 (Supplemental Fig. 6 ). All of these GO terms are consistent with viral infection and immune response in infected infants. The association of these viral and immune processes with the taxa identified in the associated cPCs allows us to hypothesize that certain microbial species may modulate the immune response in the host CD4 + T cells.

Our microbiome data analysis indicated that taxonomic diversity is highest in infants with severe infection during the acute visit, and decreases over time. Moreover, using MaAsLin and at a nominal p-value of 0.05, we identified one OTU (H. influenza) that is positively associated with RSV severity, and two OTUs (Ralstonia and Streptococcus) that are negatively associated with RSV severity. These results are partially consistent with previous observations that infants within H. influenzae-enriched clusters mount a distinct host inflammatory response characterized by the overexpression of genes related to toll-like receptor signaling and neutrophil recruitment and activation 22 . Interestingly, our subsequent network analysis implicated many of the same OTUs as MaAsLin did. Top OTUs in the network-identified clade PCs included Streptococcus, Corynebacterium, Alloiococcus, H. influenzae, Staphylococcus, Moraxella, Ralstonia, and Pseudomonas. These are consistent with a similar study analyzing the microbiome from nasopharyngeal bacterial swabs alongside whole blood transcriptomic data collected from RSV infected infants 22 . In addition, among these identified OTUs, incidence of co-infection of Streptococcus, Haemophilus influenza, and Moraxella with RSV has been reported in studies of nasopharyngeal aspirate samples from RSV infected infants [53] [54] [55] . During acute respiratory illness most infants show stable colonization of these microbes with Alloiococcus or Moraxella 21 . Our transcriptome analysis found that genes involved in immune response, such as IFITM1 and IFITM2, decreased in expression after the acute visit, while EZH1, which has also been implicated in immune response signaling 27 , was expressed at lower levels in infants with severe RSV infection. Interferons are widely expressed and are among the key genes responsible for mediating immune response 27, 56 . In particular, the IFN-inducible protein is localized in the cell membrane and endocytic vesicle and is important in restricting viral entry to the cell 57-59 . In 27 , IFITM genes have been shown to potentially inhibit RSV infection by interfering with virus entry and subsequent viral multiplication. Despite these promising findings, we note that standard differential expression only identified one gene as significantly (FDR < 0.05) differentially expressed based on RSV infection severity -EZH1 with an FDR of 0.026.

In this paper, we have proposed and applied a novel approach for integrating data from multiple omics platforms, in our case the host transcriptome and nasal microbiome, in order to extract meaningful associations. Our method constructs infant-specific correlation networks between the transcriptome and microbiome in reduced dimensions using PCA. Linear modeling of these interactions, accounting for various confounding factors, allowed us to identify associations between top gene PCs and clade PCs that differ between infants with mild versus severe RSV infection. This allowed us to better understand the etiology of RSV infection and its impact on disease severity by highlighting key associations between active components in the transcriptome and dominant constituents of the microbial composition.

The functional enrichment analysis of the top genes in the network-identified gene PCs indicated that the top-loading genes are highly enriched in pathways related to immune response and viral infection. These top genes included IFITM1 and IFITM2, which were significantly expressed between the acute and post-acute visit, but not between infants with severe versus mild infection, and CCR7, which has been found to be downregulated in RSV patients 50 . Interestingly, many of the top genes are encoded in mitochondrial DNA and included members of NADH dehydrogenase (MT-ND4L and MT-ND5), members of cytochrome c oxidase (MT-CO1, MT-CO2, and MT-CO3), and members of ATP synthase (MT-ATP6 and MT-ATP8). Mitochondrial function can be modulated by viruses 60, 61 and mitochondrial genes play a key role in the host immune response 62, 63 . These genes are also important for viral suppression of natural immunity in RSV infection 30 .

Importantly, our analysis indicates that the association of some clade PCs with gene PCs is context-dependent and a function of RSV infection severity. Further, the enrichment of immune and viral pathways in the network-identified gene PCs implies that specific microbial taxa in the nasal microbiome may impact host immune response, potentially mediating RSV infection severity. In this context, the identification of mitochondrial genes in our integrative analysis is an especially intriguing finding. The relationship between mitochondria and the microbiome is only beginning to emerge 64, 65 , and has not been previously described in the context of RSV severity. However, interestingly, the role of microbiota has been noted in the immune response to influenza 66 .

We also point out that the sign of the identified associations (negative versus positive correlations) may imply that both adverse as well as cooperative relationships exist between certain bacterial species and immune and www.nature.com/scientificreports www.nature.com/scientificreports/ viral defense response. For example, clade PC4 and clade PC8 both have positive associations with the gene PC8, indicated a potentially multifactorial host response. However, in the context of this analysis, we note that it is impossible to establish whether the microbiome-transcriptome associations we identified are indeed causal, and, if so, their direction of effect (whether the microbiome influences the transcriptome, vice versa, or both). Without data-collection in infants prior to RSV infection, we can only hypothesize as to whether targeting these associations might help to prevent or minimize the impact of RSV infection.

Finally, it is also important to note that, although highly intriguing, the associations identified in our network analysis were only nominally significant. This likely reflects heterogeneity in the disease as well as a limitation Table 3 . Top 20 genes with positive loadings in the gPCs (with associated cPCs in parentheses) identified as associated with significant edges in the sample-specific network analysis. www.nature.com/scientificreports www.nature.com/scientificreports/ in statistical power due to sample size. However, we point out that the network analysis uncovered very similar microbial taxa as MaAsLin did, lending confidence to our results. Furthermore, highly plausible biological pathways related to viral processes and immune response were implicated in RSV severity using this integrative analysis. These pathways were not found using more standard approaches, such as differential-expression analysis in the context of infection severity, and despite the fact that we did observe some immune genes differentially-expressed between the acute and post-acute visits. We believe these results illustrate the strength of using an integrative analysis approach to bring novel insights into the disease.

This study has several limitations. First, the number of samples with simultaneous collection of transcriptome and microbiome data was small. Hence most of the statistical analysis with interesting results, although nominally significant, did not show significance after adjusted for multiple hypothesis testing. Second, there may have been other host factors such as maternal or household smoking history, breast-feeding status, number of siblings, and antibiotic use, that could affect the severity of RSV infection but were not considered in our study. Third, this study was only able to provide associations between genes and microbiota without exploring the potential causal relationship.

In this study, we performed dimension reduction and differential-association analyses to quantify patterns in the nasal microbiome and host transcriptome that are associated with RSV infection severity. We also used a network-based approach to integrate these two data types and identify higher-order associations among genes and clades. This integrative analysis allowed us to systematically quantify relationships among the nasal microbiome, host transcriptome, and disease severity for RSV infection. Our results suggest that certain associations between the microbiome and transcriptome are modulated based on RSV infection severity, and that particular microbial taxa impact host immune response, with a key role for mitochondrial genes. Overall, these findings on microbiome-transcriptome associations provides novel insights into how the immune system mounts a response against RSV. Based on our findings, a key future direction for our group is to study the potential mechanisms by which the nasal microbiome impacts host immune response, such as how gene regulation is altered or otherwise affected by microbial composition, or how nasal bacteria influence the severity of RSV from a pathophysiologic perspective.

Study participants. RSV-infected infants were enrolled from three cohorts in order to capture the full spectrum of disease severity. A birth cohort was enrolled at the University of Rochester Medical Center's (URMC), Strong Memorial Hospital and Highland Hospital, and Rochester General Hospital (RGH) for two winter seasons, extending from August 15 to February 1 for 2012-2013 and 2013-2014, and followed by passive and active surveillance for development of RSV infection during the winter months (November 1-April 1). A second cohort was enrolled in pediatric offices or the emergency room at URMC's Golisano Children's Hospital or RGH when respiratory symptoms were present. The third cohort was enrolled on admission to the hospital with documented RSV infection. Eligible infants were full-term (>36 weeks gestation), healthy infants born after May 1 and less than 10 months of age at infection.

RSV-infected infants underwent evaluation by two members of the study team (a physician and a nurse). Demographic data, illness symptoms, findings on physical examination, results of laboratory and radiograph results were recorded. A nasal swab was obtained using a medium sized flocked swab (Copan Diagnostics Inc. Murrieta, CA, cat. no. 501CS01) and placed in 2 ml of sterile UV-inactivated water for quantitative RSV reverse-transcriptase polymerase chain reaction (RT-PCR). A 2-3 ml sample of heparinized blood was collected for CD4 gene expression studies as previously described 12 . RSV severity score. RSV disease severity was measured using the Global Respiratory Severity Score (GRSS) on a scale from 0-10, with higher scores representing more severe infection 67 . The GRSS was developed using an unbiased data-driven approach with nine clinical variables (including the infant's general appearance, the presence of wheezing, rales, retractions, cyanosis, lethargy, and poor air movement). The maximal age-adjusted respiratory rate, as well as the worst room air oxygen saturation, are also included in the score. The GRSS is highly predictive of other potential parameters of severity, such as need for hospitalization and duration of hospitalization 67 . The GRSS has an excellent ability to discriminate between mild and severe disease with AUC (area under the receiver operating characteristic curve) ~0.961. An infant with a GRSS ≤ 3.5 (or > 3.5) is classified as having mild (or severe) infection, respectively. For more information see 67 . We note that using continuous or discrete severity scores in our analyses resulted in similar results, so for simplicity we treated RSV infection severity as a binary variable.

CD4 + t cell extraction. CD3 + , CD4 + , CD8 − T cells were isolated from freshly collected peripheral blood samples as previously described 12, 68 . Briefly, within 24 hours of collection, Ficoll-purified peripheral blood mononuclear cells (PBMCs) were stained and sorted into major lymphocyte populations. Sorted cells were immediately lysed and homogenized in RNA extraction buffer and stored for later purification.

transcriptomic data collection and analysis. RNA extraction, processing and normalization. RNA-Seq data for CD4 + T cells were collected at two distinct time points, with 46 samples collected during the acute visit (Visit 1) and 34 samples collected on follow-up day 12-16 (Visit 2). RNA purification, library preparation, sequencing and data processing were essentially the same as previously described 12, 68, 69 . Briefly, library preparation was performed using the NexteraXT library kit (Illumina, San Diego, CA) following the SMARter Ultra Low amplification kit (Clontech, Mountain View, CA). Libraries were sequenced using the Illumina HiSeq 2500 at a target depth of ~20 million 100-bp single end reads per sample. Raw reads were mapped to Human Genome GRCh38 (annotation of GENCODE 23) and normalized by FPKM (fragments per kilobase of transcript per million reads). The transcriptome data contained reads for 11,576 genes with unique Gene Symbol annotations across 80 total samples. www.nature.com/scientificreports www.nature.com/scientificreports/ Differential expression analysis. Differential gene expression analysis was performed on FPKM normalized RNA-Seq data using the Bioconductor 70 limma package (version 3.30.13) 71 . We note that principal component (PC) analysis indicated an association of the leading PC with enrollment season. Several methods have been developed to directly regress out signals related to batch in expression data [72] [73] [74] . We used ComBat function from sva (version 3.26.0) to correct for differences from enrollment season 75 . In our differential-expression analysis we compared groups based on RSV infection severity (mild vs. severe) while correcting for sex and race. RSV viral load from nasal swab and nasal wash did not show any association with severity of infection 26 and were not included as covariates. Genes were considered significantly differentially expressed (DE) if their Benjamini-Hochberg corrected p-value was less than 0.05.

Differential variability analysis. We used var.test() in R (3.4.3 (2017-11-30) ) to perform an F-test on log2-transformed gene expression data in order to statistically quantify differences in variance in gene expression levels between the mild and severe RSV-infected groups. Genes were considered significantly differentially variable if their Benjamini-Hochberg corrected p-value was less than 0.05.

Gene Ontology pathway analysis. We used the topGO package (version 2.30.1) 41 in R to assess the enrichment of Gene Ontology (GO) pathways in a given gene list. Gene lists analyzed included (1) genes that were identified as significantly differentially-variable (see above), as well as (2) the top 100 genes associated with a gene principal component (PC) from principal component analysis on transcriptome data, based on the associated PC-loadings (see below). In the topGO (version 2.30.1) analysis, we used the 11,576 genes with measured reads in our RNA-Seq data as a background.

Microbiome data collection and analysis. 16S rRNA sequencing. Nasal swab specimens were collected during acute RSV-related illness (Visit 1) and one month later (Visit 3). As described in 69 , bacterial 16 S rRNA from these samples was extracted, amplified, and sequenced, and the resulting data were used to determine the taxonomic compositions, in terms of the relative abundances of those present operational taxonomic units (OTUs). Briefly, the V3-V4 hypervariable regions were targeted for amplification and sequenced using an Illumina MiSeq platform according to a paired end 2 × 300 bp read protocol. Preliminary read processing and quality control were performed using the Quantitative Insights into Microbial Ecology (QIIME) software package 76 , and a closed-reference OTU picking was done with USEARCH and the GreenGenes reference database 76 . The final microbiome data contained information for 1,022 distinct OTUs across 91 samples.

Multivariate association Analysis with Linear modeling (MaAsLin). We used MaAsLin (version 0.0.5) to test for significant relationships between microbial clusters and clinical outcome (severe or mild infection). MaAsLin is a multivariate statistical framework that finds associations between clinical metadata and potentially high-dimensional experimental data 76 . In contrast to transcriptomic data, the application of batch correction methods to microbiome data is still nascent 77 . Therefore, when we applied MaAsLin we adjusted for potential confounding factors, including sex and race, as well as enrollment season, which was removed by batch-correction in our transcriptomic data analysis. All microbiome samples from both Visit 1 and Visit 3 were used when running MaAsLin. integrative analysis of transcriptomic and microbiomic data. Principal Component Analysis (PCA) on transcriptomic and microbiomic data. To maximize statistical power in the integrative analysis, we first performed dimension reduction 47 . In particular, we applied principal component analysis (PCA) to host transcriptomics data to create a handful of gene principal components (gPCs), and to the nasal microbiomics data to create a few clade principal components (cPCs) 76 . To achieve that, we used the prcomp() function in the stats (version 3.4.3) package in R to perform PCA on the OTU relative abundance profiles and the FPKM-normalized RNA-Seq gene expression profiles (after applying batch-effect correction for enrollment season). In this integrative analysis, we restricted ourselves to samples from 40 subjects with both microbiome and transcriptome data collected during Visit 1. The top 13 gPCs and top 10 cPCs, which explained 95% of the variance in the transcriptomic and microbiomic data, respectively, were selected for further analysis.

Correlation between gPCs and cPCs. We constructed a Spearman correlation matrix comparing the top 13 gPCs and the top 10 cPCs, using cor() function from the stats (version 3.4.3) package in R. We treated this global 13×10 correlation matrix as the weighted adjacency matrix of a complete bipartite graph (G α , where the subscript α denotes the fact that G α is derived using all the 40 input samples) that contains two types of nodes, gPCs and cPCs.

Linear Interpolation to Obtain Network Estimates for Single Samples (LIONESS). In order to relate gPC/cPC associations to clinical variables, we applied the LIONESS method to G α to construct sample-specific (infant-specific) correlation networks 49 . LIONESS works under the assumption that the global correlation network represents a linear combination of N different networks, one from each of the N input samples. Therefore, to construct the network for sample q, we first exclude sample q and calculate the Spearman correlation matrix using the remaining samples (G (α−q) ). We then use the LIONESS equation, = − +

) , to find the network estimate for that sample (G q ). We applied LIONESS to the nasal microbiome and host transcriptome samples collected from the same set of 40 infants during Visit 1. The end result of this analysis was 40 sample-specific networks, i.e., bipartite graphs relating gPCs and cPCs.

Global, regional, and national disease burden estimates of acute lower respiratory infections due to respiratory syncytial virus in young children in 2015: a systematic review and modelling study

Determining the burden of respiratory syncytial virus disease: the known and the unknown

Incidence, risk factors and hospital burden in children under five years of age hospitalised with respiratory syncytial virus infections

Risk factors for mortality from acute lower respiratory infections (ALRI) in children under five years of age in low and middle-income countries: a systematic review and meta-analysis of observational studies

Trends in bronchiolitis hospitalizations in the United States

Risk factors for respiratory syncytial virus associated with acute lower respiratory infection in children under five years: Systematic review and meta-analysis

Fatality rates in published reports of RSV hospitalizations among high-risk and otherwise healthy children

Viral Bronchiolitis in Children

Respiratory syncytial virus and recurrent wheeze in healthy preterm infants

Asthma and allergy patterns over 18 years after severe RSV bronchiolitis in the first year of life

Adults face increased asthma risk after infant RSV bronchiolitis and reduced respiratory health-related quality of life after RSV pneumonia

Association of Dynamic Changes in the CD4 T-Cell Transcriptome With Disease Severity During Primary Respiratory Syncytial Virus Infection in Young Infants

Whole blood gene expression profiles to assess pathogenesis and disease severity in infants with respiratory syncytial virus infection

Whole blood gene expression in infants with respiratory syncytial virus bronchiolitis

The CD4 T cell response to respiratory syncytial virus infection

Respiratory syncytial virus-viral biology and the host response

Protective and dysregulated T cell immunity in RSV infection

Respiratory Syncytial Virus (RSV) Infects CD4+ T Cells: Frequency of Circulating CD4+ RSV+ T Cells as a Marker of Disease Severity in Young Children

Associations of Nasopharyngeal Metabolome and Microbiome with Severity among Infants with Bronchiolitis. A Multiomic Analysis. American journal of respiratory and critical care medicine 196

The infant nasopharyngeal microbiome impacts severity of lower respiratory infection and risk of asthma development

Nasopharyngeal Microbiota, Host Transcriptome, and Disease Severity in Children with Respiratory Syncytial Virus Infection. American journal of respiratory and critical care medicine 194

Network medicine approaches to the genetics of complex diseases

Respiratory syncytial virus infections in hospitalized infants: association between viral load, virus subgroup, and disease severity

Viral Load Dynamics and Clinical Disease Severity in Infants With Respiratory Syncytial Virus Infection

Virus-Specific Antibody, Viral Load, and Disease Severity in Respiratory Syncytial Virus Infection

Human respiratory syncytial virus infection is inhibited by IFN-induced transmembrane proteins

Interferon-induced Transmembrane Protein 1 restricts replication of virus that enter cells via the plasma membrane

Plasticity and virus specificity of the airway epithelial cell immune response during respiratory virus infection

Viral degradasome hijacks mitochondria to suppress innate immunity

Respiratory syncytial virus modifies microRNAs regulating host genes that affect virus replication

Time-course of transcriptome response to respiratory syncytial virus infection in lung epithelium cells

ISG15 Is Upregulated in Respiratory Syncytial Virus Infection and Reduces Virus Growth through Protein ISGylation

Innate immune signals modulate antiviral and polyreactive antibody responses during severe respiratory syncytial virus infection

The Human Immune Response to Respiratory Syncytial Virus Infection

The expression of B & T cell activation markers in children's tonsils following live attenuated influenza vaccine

Mucolipin 1 positively regulates TLR7 responses in dendritic cells by facilitating RNA transportation to lysosomes

Mucolipin-2 Cation Channel Increases Trafficking Efficiency of Endocytosed Viruses

MCOLN1 is a ROS sensor in lysosomes that regulates autophagy

A Cross-Study Biomarker Signature of Human Bronchial Epithelial Cells Infected with Respiratory Syncytial Virus

Improved scoring of functional groups from gene expression data by decorrelating GO graph structure

{A new class of metric divergences on probability spaces and its applicability in statistics}

UniFrac: a new phylogenetic method for comparing microbial communities

A guide to enterotypes across the human body: meta-analysis of microbial community structures in human microbiome datasets

Dysfunction of the intestinal microbiome in inflammatory bowel disease and treatment

Human microbiome analysis

Associations between host gene expression, the mucosal microbiome, and clinical outcome in the pelvic pouch of patients with inflammatory bowel disease

Sexual dimorphism in gene expression and regulatory networks across human tissues

Estimating Sample-Specific Regulatory Networks

Downregulation of IL7R, CCR7, and TLR4 in the cord blood of children with respiratory syncytial virus disease

The pathogenesis of respiratory syncytial virus disease in childhood

Shedding of L-selectin and PECAM-1 and upregulation of Mac-1 and ICAM-1 on neutrophils in RSV bronchiolitis

Incidence of bacterial coinfection with respiratory syncytial virus bronchopulmonary infection in pediatric inpatients

High incidence of pulmonary bacterial co-infection in children with severe respiratory syncytial virus (RSV) bronchiolitis

Risk of concurrent bacterial infection in preterm infants hospitalized due to respiratory syncytial virus infection

The role of IFN in respiratory syncytial virus pathogenesis

The broad-spectrum antiviral functions of IFIT and IFITM proteins

Distinct patterns of IFITM-mediated restriction of filoviruses, SARS coronavirus, and influenza A virus

The IFITM proteins mediate cellular resistance to influenza A H1N1 virus, West Nile virus, and dengue virus

Viruses as modulators of mitochondrial functions

Mitochondria Apply the Brake to Viral Immunity

Mitochondrial respiratory-chain adaptations in macrophages contribute to antibacterial host defense

Mitochondria in the regulation of innate and adaptive immunity

The Microbiome-Mitochondrion Connection: Common Ancestries, Common Mechanisms

Microbiota regulates immune defense against respiratory tract influenza A virus infection

Development of a Global Respiratory Severity Score for Respiratory Syncytial Virus Infection in Infants

Flow-based sorting of neonatal lymphocyte populations for transcriptomics analysis

The Healthy Infant Nasal Transcriptome: A Benchmark Study

Orchestrating high-throughput genomic analysis with Bioconductor

limma powers differential expression analyses for RNA-sequencing and microarray studies

Adjusting batch effects in microarray expression data using empirical Bayes methods

A general framework for multiple testing dependence

Capturing heterogeneity in gene expression studies by surrogate variable analysis

The sva package for removing batch effects and other unwanted variation in high-throughput experiments

QIIME allows analysis of high-throughput community sequencing data

Correcting for batch effects in case-control microbiome studies

This project has been funded in whole or in part with Federal funds from the National Institute of Allergy and Infectious Diseases, National Institutes of Health, Department of Health and Human Services, under Contract No. HHSN272201200005C. The authors declare no competing interests.

Supplementary information accompanies this paper at https://doi.org/10.1038/s41598-019-50217-w.

The authors declare no competing interests.Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.