ABSTRACT Respiratory syncytial virus (RSV) is the most frequently detected respiratory virus in children with acute lower respiratory tract infection. Previous transcriptome studies have focused on systemic transcriptional profiles in blood and have not compared the expression of multiple viral transcriptomes. Here, we sought to compare transcriptome responses to infection with four common respiratory viruses for children (respiratory syncytial virus, adenovirus, influenza virus, and human metapneumovirus) in respiratory samples. Transcriptomic analysis showed that cilium organization and assembly were common pathways related to viral infection. Compared with other virus infections, collagen generation pathways were distinctively enriched in RSV infection. We identified two interferon-stimulated genes (ISGs), CXCL11 and IDO1, which were upregulated to a greater extent in the RSV group. In addition, a deconvolution algorithm was used to analyze the composition of immune cells in respiratory tract samples. The proportions of dendritic cells and neutrophils in the RSV group were significantly higher than those in the other virus groups. The RSV group exhibited a higher richness of Streptococcus than the other virus groups. The concordant and discordant responses mapped out here provide a window to explore the pathophysiology of the host response to RSV. Last, according to host-microbe network interference, RSV may disrupt respiratory microbial composition by changing the immune microenvironment. IMPORTANCE In the present study, we demonstrated the comparative results of host responses to infection between RSV and other three common respiratory viruses for children. The comparative transcriptomics study of respiratory samples sheds light on the significant roles that ciliary organization and assembly, extracellular matrix changes, and microbial interactions play in the pathogenesis of RSV infection. Additionally, it was demonstrated that the recruitment of neutrophils and dendritic cells (DCs) in the respiratory tract is more substantial in RSV infection than in other viral infections. Finally, we discovered that RSV infection dramatically increased the expression of two ISGs (CXCL11 and IDO1) and the abundance of Streptococcus. KEYWORDS: host expression profile, respiratory syncytial virus, adenovirus, influenza virus, human metapneumovirus, microbiome INTRODUCTION Respiratory syncytial virus (RSV) is the most important viral pathogen causing acute lower respiratory tract infection (ALRTI) in children less than 5 years old worldwide ([44]1), posing huge medical and economic burdens globally ([45]2). However, there is currently no effective anti-RSV medicine or vaccine, and clinical practice is restricted to supportive care ([46]3). Therefore, aiming to reveal new therapeutic targets and help with vaccine research and development, insight into the pathophysiology of RSV infection is needed. The interaction between viruses and hosts is crucial in the pathogenesis of viral infection. The molecular mechanisms associated with specific diseases and individual responses to infection can be understood and described through transcriptomics ([47]4). There are currently a few studies about host expression profiles for RSV infection, but most of them are limited to systemic transcriptional profiles in blood ([48]5[49]–[50]7). Although peripheral blood mononuclear cells are involved in innate immune responses to infection, they are not directly infected by RSV and are unable to reflect immune cell composition changes in the airway, which is the primary infection site of RSV ([51]8). Therefore, analysis of the transcriptome of respiratory specimens has the potential to gain valuable insights into the pathogenesis of RSV. In addition to host-virus interactions, the host microbiome composition at the infection site frequently influences the host response against viruses indirectly and the overall risk of respiratory tract infection development ([52]9). There are reports about commensal and pathogenic microbial disruption after RSV infection ([53]7, [54]10), but few studies have reported comparisons of microbial composition between respiratory samples from patients with RSV and those with other common respiratory viruses. In this work, transcriptome sequencing (RNA-seq) was applied to profile the transcriptome and microbiome of respiratory samples from 25 patients with RSV infection, 52 patients with other virus infections, and 15 healthy children. Here, we report comparative results of host responses to infection between RSV and other viruses. We mapped their similarities and differences at the gene level, pathway level, and cell proportion level to gain insight into the specific immune signatures of RSV infection. In addition, we analyzed the active microbiome of three different groups and identified specific biomarkers of different groups. Host-microbe network interference was built to explore the relationship between the host immune response and microorganisms during RSV and other virus infections. RESULTS Host gene expression profiles among the RSV group, non-RSV group, and control group. First, we analyzed whether the transcriptome results of different sample types were highly heterogeneous before comparing the transcriptome of different viruses. Principal-coordinate analysis (PCoA) of RNA expression between 25 sputum samples (SPs) and 53 nasopharyngeal swabs (NPs) showed there was no significant heterogeneity between the nasopharyngeal swab and the airway aspirate (see Fig. S1 in the supplemental material). We enrolled and sequenced 25 and 15 airway aspirations from pediatric patients with RSV infection (defined as the RSV group) and healthy children (defined as the control group), respectively. RNA-seq identified 6,215 differentially expressed genes (DEGs; 6,117 overexpressed, 98 underexpressed) with |log[2]FC| of ≥2 and adjusted P value of <0.05, referred to as the “RSV signature” ([55]Fig. 1A). In addition, the differential expression analysis of the metatranscriptome from 52 patients with other viruses (18 patients with influenza virus [flu], 13 patients with adenovirus [AdV], and 21 patients with human metapneumovirus [hMPV], defined as the non-RSV group) identified 2,141 DEGs (947 overexpressed, 1,194 underexpressed) compared to the control group, referred to as the “non-RSV signature” ([56]Fig. 1A). The PCoA of the global expression profile demonstrated that there were substantial differences among the groups ([57]Fig. 1B). In addition, to identify a statistically significant set of genes differentially expressed in patients with RSV compared with those with other viral infections, we compared 25 patients with RSV versus 52 patients with non-RSV viral infections. A total of 6,793 DEGs were screened with |log[2]FC| of ≥2 and adjusted P value of <0.05 ([58]Fig. 1A). By comparing these DEGs with the RSV signature, 2,063 genes overlapped and were referred to as the RSV-specific gene signature ([59]Fig. 1C). For pathway analysis, the above-described RSV-specific gene signature DEGs were subjected to gene ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), and Reactome pathway enrichment analysis. GO enrichment analysis identified a set of biological processes (BPs) associated with cilium biological processes, such as cilium organization and assembly, and with extracellular component activity, such as extracellular structure and extracellular matrix (ECM) organization ([60]Fig. 1D); a set of cellular components (CCs), such as collagen-containing extracellular matrix, ion channel complex, and cation channel complex ([61]Fig. 1D); and a set of molecular functions (MFs), such as ATPase activity, active transmembrane transporter activity, and extracellular matrix structural constituents ([62]Fig. 1D). Meanwhile, KEGG analysis highlighted the involvement of protein digestion and absorption, metabolism of xenobiotics by cytochrome P450, drug metabolism-cytochrome P450, etc. ([63]Fig. 1E). Moreover, Reactome pathway analysis of this gene set identified collagen generation and assembly along with extracellular matrix as the top Reactome gene ontology terms ([64]Fig. 1F). FIG 1. [65]FIG 1 [66]Open in a new tab RNA-seq data for patients with RSV versus controls and pathway analysis of the RSV signature. (A) Volcano plots showing DEGs in control versus RSV, control versus non-RSV, and non-RSV versus RSV. Upregulated genes (red) and downregulated genes (blue) are marked. (B) PCoA of global gene expression among different viral infection groups and control groups. (C) Venn diagram showing DEGs in RSV versus non-RSV (left circle) and RSV versus control (right circle). A total of 2,063 genes differentially expressed in RSV versus control and RSV versus non-RSV represented RSV-specific gene signatures. (D to F) GO term enrichment analysis of biological processes (BP), cellular component (CC), molecular function (MF) (D), KEGG (E), and Reactome database (F). The gene ratio (x axis) is the ratio of the number of genes in our data enriched in a given gene set (pathway) to the total number of genes in that pathway. Differential immune gene expression and immune cell infiltration between the RSV and non-RSV groups. Next, to explore the difference in immune responses against viruses between the RSV group and non-RSV group, immune-related genes (IRGs) were screened by P values of <0.05 in the RSV group compared with the control and non-RSV groups. A total of 14 genes were found to be significantly differentially expressed in the RSV group compared with both the non-RSV and control groups ([67]Table 1). The heatmap showed that the expression of CCL8, IL36G, FGF2, IL17C, IFNL2, IFNL1, IFNL3, and IFNB1 was low in the non-RSV group, while their expression decreased in only a few samples in the RSV group ([68]Fig. 2A; Fig. S2). In addition, the expression levels of CXCL11 and IL36G in the RSV group were remarkably higher than those in the non-RSV and control groups; however, the expression of some IRGs, such as TNFSF12, TNFSF13, and IFNL1, in the RSV group showed the opposite trend ([69]Fig. 2A, Fig. S2, and [70]Table 1). TABLE 1. Immune-related genes with significantly different expression (P < 0.05) comparing RSV with non-RSV or controls Gene Control vs RSV P value Non-RSV vs RSV P value IFNL3 0.00034 0.00003 IFNB1 0.013 0.021 CCL18 0.0117 0.04926 CXCL11 0.000024 0.0048 TNFSF13 0.0136 0.0336 CCL8 0.0016 0.0315 IFNL2 0.0021 0.01 TNFSF12 0.00563 0.01841 IFNL1 0.00231 0.00043 PDGFRB 0.034 0.017 IL36G 0.0016 0.0014 FGF2 0.02 0.042 IL17C 0.044 0.014 CLCF1 0.0008 0.0381 [71]Open in a new tab FIG 2. [72]FIG 2 [73]Open in a new tab (A) Heatmap for differentially expressed immune-related genes and interferon-stimulated genes (ISGs) with P values of <0.05 in the RSV group versus the control group and non-RSV group. (B) Changes in cell proportions when comparing patients with RSV to non-RSV and controls. Note the trends of increased dendritic cells and neutrophils in RSV infection. The expression of interferon-stimulated genes (ISGs) between the RSV group and non-RSV group or control group was also compared, and a total of 19 ISGs were screened out ([74]Table 2) with P values of <0.05. S100A8 and DUSP5 exhibited high expression levels in all three groups ([75]Fig. 2A; Fig. S2). The expression levels of CCL8, DEFB1, ANXA2R, GPATCH11, SPSB1, MTX1, SSBP3, CDK18, and ALYREF were downregulated in the non-RSV group, while their levels showed an upregulated pattern in the RSV group. The expression levels of some ISGs, such as ABTB2, ABLIM3, CXCL11, SOCS1, and MYOF, were partially downregulated in the non-RSV group but upregulated in almost all samples of the RSV group ([76]Fig. 2A; Fig. S2). In addition, the expression levels of CXCL11 and IDO1 showed considerably increased patterns in the RSV group compared with both the control and non-RSV groups ([77]Table 1; Fig. S2). Further validation by quantitative PCR (qPCR) showed that transcription levels of CXCL11 and IDO1 were significantly increased after RSV and flu infection of A549 cells ([78]Fig. 3), but the CXCL11 level in supernatant was only significantly increased after RSV infection ([79]Fig. 4). TABLE 2. Interferon-stimulated genes with significantly different expression (P < 0.05) comparing RSV with non-RSV or controls Gene Control vs RSV P value Non-RSV vs RSV P value SSBP3 0.0008 0.02292 SPSB1 0.000038 0.039 CDK18 0.00345 0.04714 GPATCH11 0.000016 0.047 MYOF 0.02414 0.01266 CXCL11 0.000024 0.0048 CCL8 0.0016 0.0315 ABTB2 4.6E-06 0.017 SMAD3 0.0138 0.0165 ABLIM3 0.0028 0.0206 CD9 0.00022 0.00789 S100A8 0.02808 0.00012 DUSP5 0.0241 0.0043 ALYREF 0.00304 0.02158 IDO1 0.0023 0.0357 MT1X 0.044 0.031 DEFB1 0.0351 0.00026 SOCS1 0.01 0.025 ANXA2R 0.037 0.024 [80]Open in a new tab FIG 3. [81]FIG 3 [82]Open in a new tab Real-time PCR quantification of selected immune-related genes in infected and control cells at 24, 48, and 72 h postinfection. RSV strains A2, HAdV-7, and PR8 were propagated in A549 cells with a multiplicity of infection (MOI) of 0.1. Statistical analysis was performed with unpaired t test. P values of <0.05 (*), <0.01 (**), and <0.001 (***) were considered significant. Bars in the figure panels show means and standard errors of the means (SEMs). FIG 4. [83]FIG 4 [84]Open in a new tab Cell-free supernatants were analyzed for CXCL11 (left) and IL-36G (right) secretion by ELISA. Statistical analysis was performed unpaired t test. P values of <0.05 (*), <0.01 (**), and <0.001 (***) were considered significant. Bars in the figure panels show means and standard errors of the means (SEM). Furthermore, we used CIBERSORT to estimate the proportion of 22 immune cell types in bulk gene expression in respiratory samples from the RSV group, non-RSV group, and control group. In patients with RSV infection, we found that the proportion of immune cells from the innate immune system, including activated dendritic cells (DCs) and neutrophils, increased, while CD8-positive (CD8^+) T cells and follicular helper T cells in the adaptive immune system significantly decreased compared with the non-RSV group and control group ([85]Fig. 2B and Fig. S3). These results were in line with previous reports demonstrating increased DC and neutrophil counts and decreased T-cell counts in patients with RSV infection ([86]11[87]–[88]13). Microbial diversity of the active microbiome among the RSV group, non-RSV group, and control group. Next, we characterized the microbiome diversity from collected respiratory samples obtained by RNA-seq. According to our data, the Shannon diversity was not significantly different between the RSV group and non-RSV group (P > 0.05) or between the RSV group and control group (P > 0.05) ([89]Fig. 5A). The non-RSV group had a significantly greater richness index (richness, Chao1, and ACE) than the RSV group (P < 0.001) ([90]Fig. 5A), and the RSV group exhibited a relatively low richness among these three groups ([91]Fig. 5A). Streptococcus was considerably more prevalent in the RSV group than in the non-RSV group or control group at the level of the top 10 genera (median, 41.5%; false-discovery rate [q] = 0.05) ([92]Fig. 5B). Prevotella was, however, enriched in the non-RSV group (median, 11.5%; q = 0.05) ([93]Fig. 5B). These two genera occupy the oral cavity and upper respiratory tract. Among Streptococcus genera, Streptococcus parasanguinis exhibited considerably higher abundance in the RSV group than in the control group, while Streptococcus salivarius in the RSV group had the highest abundance among the three groups ([94]Fig. 5C). The total microbial composition (β-diversity) analysis revealed significant differences among the RSV group, non-RSV group, and control group (P = 0.001; permutational multivariate analysis of variance [PERMANOVA] test) ([95]Fig. 5D). Then, we utilized linear discriminant analysis (LDA) to investigate specific species within distinct groups. Streptococcus parasanguinis, Streptococcus salivarius, Streptococcus vestibularis, and Gemella sanguinis were significantly enriched in the RSV group, whereas Capnocytophaga gingivalis, Capnocytophaga sputigena, and the other 10 species were more prevalent in the non-RSV group (LDA score > 3, P = 0.05) ([96]Fig. 5E). To validated prevalence of Streptococcus both in the RSV and non-RSV groups, we did 16S rRNA sequencing to confirm the relative abundance of Streptococcus of metatranscriptome. We found a significant correlation of relative abundance between the two methods (R = 0.7, P = 3.3E-09) (Fig. S4A and B). FIG 5. [97]FIG 5 [98]Open in a new tab (A to E) Microbial community diversity within the RSV, non-RSV, and control groups. (A) Alpha diversity of the lung microbiota in cases of RSV, non-RSV, and control, measured by Shannon diversity index, richness index, Chao1 index, and ACE index. (B and C) Relative abundance of the top 10 most significantly enriched genera (B) and species (C) in the RSV, non-RSV, and control groups. (D) Composition of the upper respiratory microbiota in the RSV, non-RSV, and control groups, measured by principal-coordinate analysis via Jensen-Shannon divergence (JSD). (E) Specific species in the RSV, non-RSV, and control groups using linear discriminant analysis (LDA) score of >3; P = 0.05. Interaction between the active lung microbiome and host gene expression. To evaluate the relationship between the microbiome in the upper respiratory tract and systemic host responses induced by diverse viral infections, we examined host genes and species that were differentially expressed and active in the three groups. We eliminated correlation scores above 0.3 and P[adj] values below 0.05. The results revealed strong interactions (including correlation) between IRGs, ISGs, respiratory syncytial viruses, Streptococcus salivarius, and Streptococcus vestibularis (the top two species of airway microbiota enriched in RSV samples) ([99]Fig. 6), enabling the inference of gene-microbe connections within an integrated network. We concentrated on the differential IRGs and ISGs since they are a crucial component of the core host signature and are a group of genes that respond to RSV infections. RSV was strongly linked to numerous IRGs and ISGs, such as CXCL11, SOCS1, IFNB1, IFNL2, DEFB1, and FGF2 ([100]Fig. 6). Intriguingly, Streptococcus salivarius and Streptococcus vestibularis, two species significantly upregulated in the RSV group, were found to interact with RSV indirectly via IL36G and DUSP5, respectively, rather than directly ([101]Fig. 6). Since IL36G and DUSP5 showed significant increases in our transcriptome analysis, we further tested the results via qPCR. Transcription and secretion levels of IL36G were significantly upregulated in RSV infection in vitro samples ([102]Fig. 3 and [103]4). This implied that RSV may alter the microenvironment of the upper respiratory tract by altering the expression level of immune genes in patients, or otherwise, the presence of specific microbes colonizing the upper respiratory tract might affect the expression of specific immune genes to favor RSV infection. FIG 6. FIG 6 [104]Open in a new tab Network inference through co-occurrence analysis reveals gene-microbe associations of the immune-related and interferon-stimulated genes. Gray lines denote interactions between nodes (representing both microbes and host genes), with line thickness reflecting their observed interaction strength. Microbes within the co-occurrence network are color coded with respect to their groups. DISCUSSION In this article, we identified 2,063 RSV-specific gene signature DEGs in the RSV group compared with the non-RSV and control groups. Most of these RSV-specific DEGs exhibited downregulating patterns. While these genes were subjected to GO, KEGG, and Reactome pathway enrichment analysis, the majority of them were enriched in cilium function-related pathways. The mucociliary clearance system is the first line of defense of the respiratory tract, which initially removes pathogens through the mechanical mucus secretion of the respiratory epithelium and the pulsation of motor cilia ([105]14[106]–[107]16). RSV preferentially infects airway ciliated epithelial cells ([108]17), causing loss of cilium, overproduction of mucus, and necrosis of ciliated cells. Bronchiolitis is frequently associated with RSV infections and is characterized by mucus plug formation in the bronchiole lumens, partially due to loss of function of the mucociliary clearance system ([109]18). RSV infection promotes significant ciliary dyskinesia in the early phase ([110]19), and a sharp decrease in the ciliary beat frequency a few hours post-RSV infection has been reported ([111]20). Ultrastructural studies revealed abnormalities in RSV-infected cells, characterized by nondirectional tubular systems and atypical axonemal basal bodies ([112]21). The expression of β-tubulin and two ciliogenesis-related genes, FOXJ1 and DNAI2, also decreased ([113]21). Overall, RSV infection could lead to functional changes in the cilium and impair cilium organization and assembly. In addition to cilium-related pathways, most genes were enriched in extracellular structure organization and extracellular matrix organization pathways in GO BP enrichment analysis. According to GO CC and GO MF enrichment analysis, most of the enriched genes also referred to extracellular matrix components and extracellular matrix structural constituents. Interestingly, enrichment of extracellular matrix-related genes was specific to RSV infection compared with both the non-RSV infection and control groups. Studies have indicated that RSV infection could potentiate changes in extracellular matrix dynamics, which might be linked to asthma development ([114]22[115]–[116]24). RSV-infected human lung fibroblasts could induce hyaluronan-enriched extracellular matrix production, promoting the adhesion of mast cells, which have long been recognized as a type of critical cell in the development of allergic airway inflammation and asthma. The binding between mast cells and the ECM was accompanied by augmented expression of proteases, chymase, and tryptase by mast cells, which led to a proinflammatory milieu and tissue damage during RSV infections ([117]22). An in vivo study on juvenile and adult RSV-infected mice suggested that sustained myeloid infiltration observed in juvenile mice was related to extracellular matrix remodeling, facilitating the conditions for the development of asthma later in life ([118]25). Xu et al. indicated that RSV infection could change the topographical map of chromatin, forming nucleosome-free regions on genes responsible for ECM production and remodeling, including TGFB1, MMP9, and GFPT2, which encode the rate-limiting enzyme in the hexosamine biosynthetic pathway ([119]26), leading to epigenetic reprogramming and tissue remodeling ([120]26). Our study provided evidence at the transcript level that RSV infection could alter ECM-related gene expression patterns, leading to extracellular matrix production and airway remodeling, and might be associated with long-term airway hyperresponsiveness, which makes those infected with RSV in early life predisposed to asthma. We also found that CXCL11 and IDO1 genes exhibited high expression levels in the RSV group. CXCL11, also known as interferon gamma (IFN-γ)-inducible protein 9 (IP-9) or interferon-inducible T-cell alpha chemoattractant (I-TAC), is an immune chemokine belonging to the CXC family of chemokines and mainly regulates immune cell recruitment, differentiation, migration, and activation ([121]27). It is induced by IFN-γ and IFN-β and weakly by IFN-α and has the highest affinity for CXCR3 ([122]28). Studies have indicated that CXCR3 signaling pathways play a crucial role in the priming and activation of CD4^+ T cells and in the balance between effector CD4^+ T cells (Th1 and Th17) and regulatory CD4^+ T cells (Th2 and Tr1) ([123]29). The CXCL11-CXCR3 axis has been linked to inflammation, autoimmune disease, and cancer. The binding of CXCL11 to CXCR3 leads to internalization of CXCR3 ([124]30), making it less accessible to CXCL9/CXCL10 and shifting the Th1/Th2 balance into Th2 ([125]29). CXCL11/CXCR3 binding induces an immunotolerizing state via p70 kinase/mTOR in STAT3- and STAT6-dependent pathways, characterized by high levels of interferon 10 (IL-10^hi) Tr1 and IL-4^hi Th2 cells, to restrain inflammation ([126]31). Many studies have indicated that CXCL11 levels increase to various extents in response to viral infections, such as influenza virus ([127]32), SARS-CoV-2 ([128]33), human immunodeficiency virus ([129]34), and hepatitis B virus ([130]35). During RSV infection, Th2-biased polarization, which might be a consequence of high levels of CXCL11, may be associated with severe clinical manifestations. Severe/hypoxic bronchiolitis exhibited a higher level of IL-4 (a Th2 marker) than that of a group with nonhypoxic bronchiolitis ([131]36). A higher IL-4/IFN-γ ratio in nasal lavage fluid was detected in RSV-infected infants with acute bronchiolitis than in their counterparts with only RSV-related upper respiratory tract infection alone ([132]37). A cohort study also indicated that children with hypoxic RSV low respiratory tract infection showed a predominance of Th2 cytokines ([133]38). Thus, our results indicated that CXCL11 might play an important role in the pathogenesis of RSV infection by inducing a Th2 polarization immune response, although further studies are needed to investigate the specific mechanism. Indoleamine 2,3-dioxygenase (IDO) is a rate-limiting enzyme involved in tryptophan-related metabolites ([134]39). IDO regulates the immune response in vitro and in vivo, mainly shifting the Th balance toward Th2 and in favor of the Th2-related cytokine profile ([135]40). It was also indicated that IDO was involved in the selective apoptosis of proliferating Th1 cells ([136]40). Only a few studies have focused on the relationship between RSV and IDO, and these results are debatable. A study indicated that IDO played a protective role during RSV infection. The antiviral ability of IFN-γ in RSV pathogenesis was probably mediated by IDO since IDO was potently induced by IFN-γ and IDO knockdown in vitro could diminish the antiviral effect of IFN-γ ([137]41). However, other studies mainly focused on the adverse effects of IDO during RSV infection. The upregulation of IDO in RSV-infected mesenchymal stem cells might affect immune cell proliferation, which might be responsible for the lack of protective immunity against RSV and for the long-term outcomes of RSV-induced lung diseases such as asthma and chronic obstructive pulmonary disease ([138]42). IDO bioactivity and expression could be induced by RSV in monocyte-derived DCs in a virus replication-dependent manner, which inhibited T-bet (a Th1-cell marker) expression on T cells and thus led to Th2 polarization. A lower Th1 response to RSV might hinder viral clearance and probably promote the development of allergic diseases ([139]43). Jin et al. showed that RSV infection reduced kynurenic acid production and inhibited the transformation from Th17 to Foxp3^+ regulatory T cells (Tregs) (Th17/Treg balance) by altering IDO activity in plasmacytoid DCs (pDCs) ([140]44), which might disrupt asthma tolerance later in life ([141]45). High expression of IDO1 in bat neutrophils was observed in a single-cell transcriptome analysis study of viral infection response in bats, which may play an important role in limiting inflammation in the context of neutrophil-activated conditions ([142]46). In addition, El-Zaatari, M. et al. showed that neutrophil recruitment to infected cecum was reduced in IDO1-deficient mice ([143]47), while the increase of local concentration of l-kynurenine, a reaction product of IDO1, inhibited neutrophil entry into infected tissues ([144]48). Whether RSV utilizes IDO1 to attenuate neutrophil chemotaxis in lung tissue and thus reduce lung immune response remains to be investigated. Our results are in accordance with previous studies that indicated that neutrophils constitute the major part of infiltrating cells in the airways of infants with RSV bronchiolitis ([145]11). Neutrophils exert antimicrobial effects by phagocytosis, degranulation, secretion of broad antimicrobial proteins, and formation of neutrophil extracellular traps (NETs) ([146]49). Neutrophils can help eliminate infected cells and reduce the amount of infectious virus in vitro ([147]50). Human cathelicidin LL-37, which is primarily expressed by neutrophils and epithelial cells, has antiviral activity against RSV by directly affecting RSV viral particles and exerts protective effects against cell death ([148]51). RSV infection could induce classical reactive oxygen species (ROS)-dependent NET release, which was able to bind and isolate RSV virion particles, inhibiting them from reaching the target cells ([149]52). Serine proteases released from NETs can also interact with and inactivate RSV ([150]53). However, deteriorating effects such as tissue damage and systemic symptoms come with the influx and infiltration of neutrophils in response to RSV infections. Studies have revealed that the extent of neutrophil infiltration is positively associated with clinical severity ([151]54, [152]55). Greater neutrophil transepithelial migration, activation, and degranulation were associated with poorer ciliary function and more severe epithelial cell damage ([153]56). Infiltration of neutrophils accumulating in the small airway impairs gas exchange ([154]57) and promotes mucus production, which also leads to airway obstruction ([155]58). Activated neutrophils can release antimicrobial mediators, such as elastase ([156]59) and neutrophil-derived ROS, exerting antiviral effects ([157]60), but are also cytotoxic to host cells at the same time. Thus, although neutrophils could help us eliminate the infection, inappropriate activation of neutrophils during RSV infection would lead to adverse outcomes. In addition to neutrophils, we found that DCs in the RSV group were significantly higher than those in the non-RSV group and control group. DCs can initiate and regulate adaptive immune responses according to peripheral cues, forming a remarkable cellular network against foreign and self-antigens ([158]61). Although DCs may not be optimal permissive cells for RSV, RSV infects DCs probably as a strategy to interfere with the activation of antigen-specific immune responses and assist their infectious process indirectly. The interaction between RSV-infected DCs and T cells could evade adaptive immunity and inhibit T-cell-mediated virus clearance by interfering with immunological synapse assembly ([159]62). RSV-exposed monocyte-derived DCs upregulate IFN-α and IFN-γ, suppressing the proliferation of CD4^+ T cells ([160]63). RSV could also alter the balance in the proportion of DC phenotypes in pathogenesis. pDCs play a protective role during RSV, characterized by controlling RSV replication and enhancing their elimination, suppressing pulmonary inflammation, and inhibiting the development of airway hyperresponsiveness ([161]64), while lung myeloid DCs (mDCs) may have detrimental effects on the host during RSV infection ([162]13). The number of pDCs was lower and that of mDCs was higher in the peripheral blood of children with RSV bronchiolitis than in healthy controls, and they still exhibited significantly lower quantities of pDCs even 3 months after discharge ([163]65). Allergen sensitization and severe lung pathology caused by RSV infection are likely responsible for the accumulation and activation of mDCs, which might skew Th cell polarization toward a pathogenic direction. The turbulence caused by RSV in DC balances might have a long-term effect on the host and be responsible for the development of asthma later in life ([164]66). The differences in microbial diversity and composition among the RSV group, non-RSV group, and control group were also analyzed. We demonstrated that patients with RSV infection exhibited a relatively lower microbiological richness in their respiratory tract than patients in the non-RSV group or control group. Compared with the other groups, the RSV group was associated with a change in the abundance of specific microbiota, distinguished by a significant increase in Streptococcus, which was in accordance with previous studies ([165]7, [166]10). Among the species of Streptococcus, the abundance of Streptococcus salivarius and Streptococcus vestibularis was significantly higher than that in either the non-RSV or control groups, indicating that these two species of Streptococcus might be involved in the pathogenesis of RSV infection. Integrative analysis of the active lung microbiome and transcriptome data could help us gain insight into the relationship between microbial composition and host genes. Thus, using gene and microbe coexpression network analysis, we characterized associations between host responses to RSV infection and the microbiome. RSV was positively correlated with some IRGs and ISGs, including CXCL11, IDO1, DUSP5, and IL36G. Among these genes, IL36G/IL-36γ was positively associated with Streptococcus salivarius, while DUSP5 was positively correlated with Streptococcus vestibularis. IL-36 belongs to the IL-1 family, comprising IL-36α, IL-36β, IL-36γ, and IL-36Ra. IL-36γ can be expressed by macrophages, monocytes, and bronchial epithelial cells under stimulation by proinflammatory cytokines, including IL-1, tumor necrosis factor (TNF), IL-17, Pseudomonas aeruginosa, and Toll-like receptor 3 (TLR3) ligands ([167]67). IL-36 is a potent regulator of DCs, characterized by upregulation of markers of activated DCs and increases in IL-6 and IL-12 production ([168]68). IL-36γ might play a role in the development of allergen-induced inflammation in the lung. A study indicated that IL-36γ was upregulated in the lungs of asthma-susceptible A/J mice compared to asthma-resistant C3H/HeJ mice ([169]69). Another study demonstrated that IL-36γ promoted airway hyperresponsiveness, induced neutrophil influx, enhanced NF-κB activity, and increased chemokine production ([170]70). In addition, IL-36γ was significantly upregulated at the mRNA level in RSV-infected A549 cells. Thus, IL-36γ might contribute to the pathogenesis of RSV infection and potentially favor secondary Streptococcus salivarius infection. The expression of DUSP5 was upregulated in clinical samples and also elevated in in vitro samples. We presumed that DUSP5 was indirectly involved in RSV infection in the presence of Streptococcus vestibularis, which was absent in cell cultures. DUSPs are key regulators in the MAPK pathway that dephosphorylate key signaling molecules ([171]71). DUSP5 shows phosphatase activity toward ERK (one of the best-known members of the MAPK family), leading to both activation and nuclear translocation of ERK ([172]72). ERK-1/2 activity is required for efficient RSV infection and replication ([173]73). Additionally, RSV-induced ERK-1/2 expression was involved in the posttranscriptional upregulation of RANTES, which was linked to severe RSV disease ([174]74, [175]75). Thus, DUSP5 probably promotes RSV infection via the ERK pathway. However, our study only indicated that RSV might interact with Streptococcus salivarius and Streptococcus vestibularis via IL36G and DUSP5, respectively. If RSV infection leads to a change in the abundance of these two bacteria via these two genes, further studies will be needed to elucidate the specific mechanisms. One limitation of the present work was the relatively small sample size due to the difficulty encountered in collecting respiratory specimens of children, thus restricting the statistical power of the analysis used for identifying differences between groups. As a result, we only concentrated on the characteristics that were most closely associated with upper respiratory virus, microbiota, and host immune response. This may have made it impossible for the inquiry to identify more than a small number of potential microbiome types. The true respiratory microbiomes may be more diverse, and they may have more complex interactions with the human immune system. Nevertheless, the information pertaining to the stratified organization of the host immune signature and its correlation with different traits is statistically significant and supports the notion of an ongoing host-microbe interaction. Based on previous studies on metagenome data, it is alluring to anticipate that the variance in the microbiome is stratified in several dimensions in distinct viral infections. This stratification may represent changes in the homeostasis of the host and microbiome caused by various viruses. Collectively, the comparative transcriptomics study of respiratory samples sheds light on the significant roles that ciliary organization and assembly, extracellular matrix changes, and microbial interactions play in the pathogenesis of RSV infection. Additionally, it was demonstrated that the recruitment of neutrophils and DCs in the respiratory tract is more substantial in RSV infection than in other viral infections. Finally, we discovered that RSV infection dramatically increased the expression of two ISGs (CXCL11 and IDO1) and the abundance of Streptococcus. MATERIALS AND METHODS Subjects and clinical samples. Nasopharyngeal swabs (NPs) or sputum samples (SPs) were collected from children diagnosed with pneumonia in the National Multicenter Pneumonia Pathogen Surveillance Program organized by Beijing Children’s Hospital. Total white blood cell counts less than 12 × 10^9/L, C-reactive protein concentration in serum less than 20 mg/L, procalcitonin concentration in serum less than 0.5 ng/mL, and negative sputum culture were considered to rule out bacterial infection on the first routine blood test after admission ([176]76). Samples with only one respiratory virus pathogen identified by PCR analysis were included in the study. Other inclusion criteria were age of less than 18 years and the presence of a parent or guardian capable of providing informed consent. Exclusion criteria included any underlying medical condition that required regular medical care, receipt of immunosuppressive medications, including corticosteroids, within the preceding 30 days and receipt of antibiotics within 7 days. Children hospitalized for selected surgery (such as hypospadias or hydronephrosis) and with no signs of any infection were enrolled as a control group. Demographic and laboratory parameters were obtained for each enrolled patient (see Table S1 in the supplemental material). Metatranscriptome sequencing and data processing. According to the manufacturer's recommendations, a 1-μL sample was processed with Turbo DNase (Life Technologies, USA) to deplete the host DNA background. RNA was extracted using a QIAamp UCP pathogen minikit (Qiagen, Valencia, CA, USA), reverse transcribed, and amplified using an Ovation RNA-Seq system (NuGen, CA, USA). Following fragmentation, the library was constructed using Ovation Ultralow system v2 (NuGen, CA, USA) and sequenced on an Illumina NextSeq 550 (single-end 75 bp) ([177]77). Raw sequencing data were processed using fastp ([178]78) to remove reads containing adapters or ambiguous “N” nucleotides and low-quality reads. All clean data were mapped to the human genome hg19 using HISAT2 ([179]79) with default parameters. FeatureCounts ([180]80) was used to quantify expression at the gene level. Read count normalization and differentially expressed analyses were performed using the DESeq2 package ([181]81). Differentially expressed genes (DEGs) were determined with an adjusted P value of <0.05 and absolute log fold change (log[2]FC) of ≥2. The clusterProfiler package ([182]82) was used for gene ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), and Reactome pathway enrichment analysis of DEGs. Benjamini-Hochberg-adjusted P values of <0.05 showed significant enrichment. The immune-related gene (IRG) list was downloaded from the Immport database ([183]83), and interferon (IFN)-stimulated genes (ISGs) were retrieved from the study of Wu et al. ([184]84). According to the gene expression profile of IRGs and ISGs, the Kruskal-Wallis test was used to select genes that were significantly different from the other two groups in the RSV group and visualized by the pheatmap package ([185]85). To infer the composition of immune cells, the CIBERSORT algorithm ([186]86) with the original CIBERSORT gene signature file LM22 and 1,000 permutations were used to examine the relative proportions of 22 invasive immune cell types in each sample. Microbial taxonomy assignment and host-microbial network interaction. Taxonomy assignment of the upper respiratory microbiome was performed by Kraken 2 ([187]87) ([188]https://ccb.jhu.edu/software/kraken2/) with the Standard database ([189]https://benlangmead.github.io/aws-indices/k2). The taxonomic relative abundance was calculated using Bracken (version 2.5). The alpha diversity was calculated using the Shannon index, richness, and Chao1 and ACE estimators. Wilcoxon rank sum tests were used to evaluate the differences in α-diversity and species richness among the three different groups ([190]77). For β-diversity analysis, the principal-coordinate analysis scaling (PCoA) plot was based on Jensen-Shannon divergence (JSD) ([191]88). The PERMANOVA test was used to determine the overall community composition difference between the elderly and young samples. These calculations and visualizations were implemented using the R package vegan. To elucidate the microbial and host co-occurrence patterns, species with frequencies greater than 40% were kept for correlation calculations using the SparCC algorithm ([192]89). Only significant (P < 0.05) and robust (SparCC |r| > 0.3) correlations were retained for further co-occurrence network analyses. Visualization of networks was then implemented in Cytoscape ([193]90). 16S rRNA sequencing and data processing. DNA was extracted using a QIAamp UCP pathogen minikit (Qiagen, Valencia, CA, USA). The V3-V4 hypervariable region of the bacterial 16S rRNA gene was amplified with barcoded primer set 341F (CCTAYGGGRBGCASCAG) and 806R (GGACTACNNGGGTATCTAAT) with an expected amplicon length of 450 to 475 bp. Sequencing of the amplicons was performed using an Illumina NovaSeq 6000 instrument (Illumina, USA) (250-bp read length, paired-end protocol). Reads were analyzed by QIIME2 using the SILVA database ([194]91, [195]92). Cells and viruses. Type II alveolar epithelial cells (A549) (CCL-185) were obtained from ATCC (Manassas, VA, USA). The cells were maintained in Dulbecco’s modified Eagle medium (DMEM; Gibco) supplemented with 10% fetal bovine serum (Gibco) and 1% penicillin-streptomycin solution (Gibco) at 37°C and 5% CO[2]. RSV strain A2, human adenovirus type 7 (HAdV-7), and influenza A strain PR8 were propagated in A549 cells before quantification of viral titer by 50% tissue culture infective dose (TCID[50]) assay. Cell infections, qPCR, and ELISA. All in vitro infection experiments were performed with a multiplicity of infection (MOI) of 0.1. In vitro samples were collected at 24, 48, and 72 h postinfection, cells were subsequently lysed with RLT buffer, and RNA was extracted using a Qiagen RNeasy kit (Qiagen, Valencia, CA, USA). cDNA was reverse transcribed from RNA extracts (500 ng) using a PrimeScript one-step reverse transcriptase PCR (RT-PCR) kit version 2 (TaKaRa, Japan). qPCR was carried out using TB Green premix Ex Taq (Tli RNase H Plus) (TaKaRa, Japan) with the primers listed in [196]Table 3. Fold changes of target genes in infected samples were calculated by using the threshold cycle (ΔΔC[T]) method and normalized to a GAPDH (glyceraldehyde-3-phosphate dehydrogenase) mRNA level ([197]93). Cell-free supernatants were frozen and thawed before quantification of CXCL11 and IL36G was performed using a human CXCL11 enzyme-linked immunosorbent assay (ELISA) kit (Abcam) and human IL36G ELISA kit (Signalway Antibody), respectively. Statistics were performed with GraphPad Prism 8 software using unpaired t test. P values of <0.05 (*), <0.01 (**), and <0.001 (***) were considered significant. Bars in the figure panels show means and standard errors of the means (SEMs). TABLE 3. Primer sequences used in this study Gene name Direction Sequence (5′–3′) CXCL11 Forward GACGCTGTCTTTGCATAGGC Reverse GGATTTAGGCATCGTTGTCCTTT IDO1 Forward TCTCATTTCGTGATGGAGACTGC Reverse GTGTCCCGTTCTTGCATTTGC IL36G Forward AGGAAGGGCCGTCTATCAATC Reverse CACTGTCACTTCGTGGAACTG DUSP5 Forward GCGACCCACCTACACTACAAA Reverse CTTCATAAGGTAAGCCATGCAGA [198]Open in a new tab Ethics statement. This study was performed in strict accordance with the human subject protection guidance and was approved by the Ethical Review Committee of Beijing Children’s Hospital. Written informed consent was obtained from the participants’ parents or guardians. Data availability. Meta-transcriptomic sequencing data were deposited in the Genome Warehouse in the National Genomics Data Center (National Genomics Data Center Members and Partners, 2021) under project [199]PRJCA014022, which is publicly accessible at [200]https://ngdc.cncb.ac.cn/. ACKNOWLEDGMENTS