The difference of oropharyngeal microbiome during acute respiratory viral infections in infants and children
Viruses are common causes of acute respiratory infections (ARI) in infants and children. The most frequently viral pathogen was respiratory syncytial virus (RSV) followed by, human rhinovirus (HRV), human parainfluenza virus (HPIV), influenza virus (IFV), human adenvirus (HAdV), and human metapneumovirus (HMPV) in China , accounting for great morbidity and mortality. For example, among children younger than five years globally, it is estimated that there were 643,000 HMPV-associated hospital admissions and 7,700 HMPV-associated in-hospital deaths in 2018 ; and 870,000 IFV associated hospital admissions, and 34,800 overall IFV-associated deaths in 2018 ; furthermore, 3.6 million RSV-associated hospital admissions and 26,300 RSV-associated acute lower respiratory infection in-hospital deaths in 2019 . Although previous studies have identified a number of viral agents, such as IFV, HPIV, HRV, HAdV, and enterovirus (EV) that leads to severe disease requiring hospitalization, in some cases, however, the severe diseases happen among children without any known risk factors such as preterm birth or immunosuppression . The respiratory indigenous microbiota may play an important role in affecting the susceptibility and outcome of acute respiratory infections. For example, the susceptibility to infection of influenza is associated with the upper respiratory tract microbiome at the time of exposure . The composition of microbiota can cause individuals to be associated with the risk of infection . One potential mechanism is that the respiratory microbiota can mediate colonization resistance by preventing the growth and persistence of non-indigenous microorganisms through competing for nutrition and functioning . This microbiota can also stimulate and alter host mucosal immune defenses to prevent invasion of pathogenic microorganisms . Recent studies have shown that there are changes of microbiota during IFV, RSV, or HRV infections – . However, limited study has been focusing on a variety of other common respiratory viruses such as HMPV, HPIV, HAdV, or EV. In this study, we used hospitalization as measure of disease severity to characterize the oropharyngeal microbiome from non-hospitalized and hospitalized infants and children with ARI of one of the seven common respiratory viruses or ARI of no respiratory virus detectable, to determine whether respiratory virus infection associated with the alteration of microbial structure and functional pathways using whole-genome shotgun metagenomic sequencing. 

Overview of the study population and oropharyngeal microbiome This study originally selected 154 infants and children aged less than 12 years from six hospitals in China. Forty-five samples were excluded due to no valid sequencing data ( N = 33), rarefaction ( N = 10), and potential contamination ( N = 2), and finally 109 samples were included. The characteristics of participants are shown in Table and . The median number of microbial reads per sample was 355,052 (IQR 18,888–3,495,916). All 49 inpatients and 60 outpatients were less than 12 years old, with the majority aged 0-5 years (61.2% of inpatients and 58.3% of outpatients) and a slightly higher proportion of male (57.10%, 53.3%). As of the disease-related symptoms, more than half of the participants had fever (59.2%, 43.3%), cough (83.7%, 76.7%), and sore throat (61.2%, 56.7%), whereas only a few participants reported headache (8.2%, 8.3%) and fatigue (4.1%, 5.0%). The overall taxonomic alpha diversity was not statistically different between inpatients and outpatients as measured by Shannon index and observed feature ( Fig. ), and also not differed by the types of infected virus (Shannon index P = 0.092, observed feature P = 0.461; Figs. and ). However, the difference was not significant within each respiratory virus infected group when compared inpatients and outpatients (Figs. and ), or patients with and without specific disease-related symptoms (Figs. and ). In the beta diversity analyses, the matrices based on Bray Curtis and Jaccard were significantly different between inpatients and outpatients as indicated by the first two principal coordinates analysis (PCoA) vectors (Figs. and ) and the Microbiome Regression-Based Kernel Association Test (MiRKAT) models (both Bray Crutis and Jaccard P < 0.001; Tab. ). Significant difference was also found between patients with disease-related symptoms, including fever (Bray Crutis P = 0.022; Jaccard P = 0.011), congestion (Bray Crutis P = 0.017; Jaccard P = 0.037), and oropharyngeal hyperaemia (Bray Crutis P = 0.070; Jaccard P = 0.021; Table ), compared to those without. The overall beta diversity also differed across infected virus groups (Figs. and ). For the analyses within each virus-infected group, the beta diversity was significant different between HRV-infected inpatients and HRV-infected outpatients (Bray Crutis P = 0.040; Jaccard P = 0.089), as well as between HMPV infected inpatients and HMPV infected outpatients but only observed in Jaccard matrix (Bray Crutis P = 0.241; Jaccard P = 0.008; Figs. – and Tab. ). As measured by the distance-based coefficient of determination ( R 2 ) using permutational multivariate analysis of variance (PERMANOVA) test, disease severity (inpatient or outpatient) explained most of the overall microbial community composition (i.e., beta diversity), followed by several disease-related symptoms including oropharyngeal hyperaemia, fever, and congestion. Most respiratory virus infection explained less than 1%, except HAdV and HMPV significantly correlated with Bray-Curtis matrix and HRV significantly correlated with Jaccard (Fig. ). Oropharyngeal microbes associated with disease severity and disease-related symptoms The relative abundance of oropharyngeal phyla differed by disease severity and infected virus (Fig. ). The most abundance phyla in our study population were Proteobacteria (42.4% for inpatient, 59.4% for outpatient), Firmicutes (30.2%, 26.7%), Actinobacteria (17.6%, 12.5%), Bacteroidetes (7.6%, 0.9%), and Fusobacteria (1.3%, 0%). The presence of Proteobacteria was significantly associated with disease severity. The odds of being inpatients were 0.27 (95% CI = 0.10–0.73, P = 0.010) with the presence of Proteobacteria (Fig. ). The relative abundance of bacterial phyla was associated with respiratory virus infections. For every increase in one percent of the relative abundance, Proteobacteria (OR = 7.62, 95% CI = 1.92–30.32, P = 0.004) increased, and Firmicutes (OR = 0.13, 95% CI = 0.02–0.78, P = 0.026) and Actinobacteria (OR = 0.03, 95% CI = 0.001–0.62, P = 0.024) decreased, the odds of HMPV infection; Bacteroidetes (OR = 9.91, 95% CI = 1.44–68.36, P = 0.020) increased the odds of HAdV infection; Actinobacteria increased the odds of HPIV infection (OR = 6.59, 95% CI = 1.28–33.89, P = 0.024) and decreased the odds of HRV infection (OR = 0.01, 95% CI = 0.0001–0.90, P = 0.045; Fig. ). The 15 most abundance species in our study population were Stenotrophomonas maltophilia (14.4% for inpatient, 9.9% for outpatient), Rothia mucilaginosa (14.0%, 8.0%), Brucella anthropi (0.9%, 17.3%), Achromobacter aegrifaciens (4.6%, 11.1%), Achromobacter marplatensis (2.8%, 8.7%) Streptococcus infantis (5.1%, 3.2%), Streptococcus parasanguinis (2.4%, 5.2%), Streptococcus mitis (2.1%, 3.9%), Streptococcus salivarius (2.6%, 3.5%), Pseudomonas oleovorans (2.3%, 2.9%), Neisseria subflava (3.7%, 0.9%), Phyllobacterium myrsinacearum (4.1%, 0%), Dolosigranulum pigrum (2.0%, 1.6%), Cutibacterium acnes (3.2%, 0.1%), and Haemophilus influenzae (2.1%, 0.9%; Fig. ). Among all the microbial species detected in our study population, sixteen species had more than twofold higher abundance in outpatients than inpatients (Fig. and Tab. ), where Enterobacter kobei , A. marplatensis , A. aegrifaciens , Pseudomonas putida , Achromobacter xylosoxidans , B. anthropi , Achromobacter pestifer , Ochrobactrum sp.CGA5 , and Achromobacter insolitus had an false discover rate (FDR) less than 0.05. Five species had higher abundance in inpatients compared to outpatients, but not statistically significant after FDR adjustment. When restricted this analyses within each virus infected group, A. marplatensis , A. xylosoxidans , and A. aegrifaciens among HMVP infection patients, and A. xylosoxidans and Delftia acidovorans among respiratory virus undetectable patients were twofold higher abundance in inpatients than in outpatients (Fig. ). However, no species were significantly different at an FDR of 5% for the infection group stratified analyses. According to the LDA score by linear discriminant analysis and effective size comparisons (LEfSe) analysis, fifteen taxa from order to species levels (order Burkholderiales and Enterobacterales ; family Brucellaceae , Alcaligenaceae , and Enterobacteriaceae ; genus Brucella , Achromobacter , and Enterobacter ; species B. anthropi , A. aegrifaciens , A. marplatensis , P. putida , A. xylosoxidans , E. kobei , and A. pestifer ) were enriched among outpatients, but no taxa were found to be enriched among inpatients (Fig. , ). When analyzed the taxonomic abundance of disease-related symptoms, statistical significance was only found between patients with fatigue and those without, with six taxa enriched in patients with fatigue at a FDR of 5% (Table ). Oropharyngeal community types and participant characteristics Dirichlet multinomial mixture (DMM) model was used to unsupervised clustering the sample to three community types, with 60 samples assigned to community type 1 (CT1), 27 samples assigned to community type 2 (CT2), and 22 samples assigned to community type 3 (CT3). The alpha diversity of CT1 and CT2 was significantly lower than CT3 (Fig. ), and the overall microbial composition was significant different between the three groups (Fig. , ). Among the top driver species (Fig. ), the relative abundance of R. mucilaginosa was the highest within each community, but was higher in CT1 and CT3 than in CT2 when compared across different communities. S. maltophilia and A. aegrifaciens were enriched in CT1; A. xylosoxidans and E. kobei were enriched in CT2; and S. parasanguinis and S. salivarius were enriched in CT3 (Fig. ). The distribution of CT1 and CT3 was similar between inpatients and outpatients, whereas the proportion of CT2 was higher in outpatients than in inpatients (Fig. ). The prevalence of CT1 was higher among patients with IFV, EV, and HAdV infections; CT2 was higher among patients with HMPV infections; and CT3 was higher among patients with RSV and HAdV infections. Participants of CT1 and CT2 were more likely to have cough (81.7% of CT1, 85.2% of CT2), and participants with CT3 were more likely to have fever (77.3%). Microbial pathways associated with disease severity and disease-related symptoms A total of 431 microbial pathways were identified, with 334 pathways present in 5–95% of the population. Among them, 244 pathways had more than twofold higher abundance in outpatients compared to inpatients, and 231 were statistically significant at an FDR of 5% (Fig. and Table ). The quantile-quantile plot of P value showed that most of the pathways had a higher observed P value than the expected P value (Fig. ). Many of these pathways are linked to metabolic degradation, biosynthesis, or cellular respiration. For example, the top 25 most distinctive enriched pathways in outpatients include 13 degradation pathways 14DICHLORBENZDEG-PWY (1,4-dichlorobenzene degradation), CATECHOL-ORTHO-CLEAVAGE-PWY (catechol degradation to β-ketoadipate), GALACTARDEG-PWY (D-galactarate degradation I), GLUCARGALACTSUPER-PWY (superpathway of D-glucarate and D-galactarate degradation), PWY-1361 (benzoyl-CoA degradation I), PWY-5417 (catechol degradation III), PWY-5431 (aromatic compounds degradation via β-ketoadipate), PWY-561 (superpathway of glyoxylate cycle and fatty acid degradation), PWY-6084 (3,5-dichlorocatechol degradation), PWY-6182 (superpathway of salicylate degradation), PWY-6185 (4-methylcatechol degradation), PWY-8187 ( l -arginine degradation XIII), and TYRFUMCAT-PWY ( l -tyrosine degradation I), four biosynthesis pathways ARGSYNBSUB-PWY ( l -arginine biosynthesis II), PWY-5855 (ubiquinol-7 biosynthesis), PWY-6292 (superpathway of l -cysteine biosynthesis), and PWY-7400 ( l -arginine biosynthesis IV), as well as a six other metabolic pathways P105-PWY (TCA cycle IV), P124-PWY (Bifidobacterium shunt), PWY-5138 (fatty acid β-oxidation IV), PWY-5676 (acetyl-CoA fermentation to butanoate II), PWY-6969 (TCA cycle V), and PWY66-391 (fatty acid β-oxidation VI; Fig. ). When analyzed the pathway abundance of disease-related symptoms, no significant difference was found at an FDR of 5% (Fig. ). Sensitivity analysis As we had only one HPIV-infected inpatient with a valid sequencing result, the asymmetry of inpatients and outpatients in this group may bias our results. Therefore, we conducted a sensitivity analysis that excluded all HPIV-infected patients. Results were similar when this group of patients was excluded (Fig. ). 

This study originally selected 154 infants and children aged less than 12 years from six hospitals in China. Forty-five samples were excluded due to no valid sequencing data ( N = 33), rarefaction ( N = 10), and potential contamination ( N = 2), and finally 109 samples were included. The characteristics of participants are shown in Table and . The median number of microbial reads per sample was 355,052 (IQR 18,888–3,495,916). All 49 inpatients and 60 outpatients were less than 12 years old, with the majority aged 0-5 years (61.2% of inpatients and 58.3% of outpatients) and a slightly higher proportion of male (57.10%, 53.3%). As of the disease-related symptoms, more than half of the participants had fever (59.2%, 43.3%), cough (83.7%, 76.7%), and sore throat (61.2%, 56.7%), whereas only a few participants reported headache (8.2%, 8.3%) and fatigue (4.1%, 5.0%). The overall taxonomic alpha diversity was not statistically different between inpatients and outpatients as measured by Shannon index and observed feature ( Fig. ), and also not differed by the types of infected virus (Shannon index P = 0.092, observed feature P = 0.461; Figs. and ). However, the difference was not significant within each respiratory virus infected group when compared inpatients and outpatients (Figs. and ), or patients with and without specific disease-related symptoms (Figs. and ). In the beta diversity analyses, the matrices based on Bray Curtis and Jaccard were significantly different between inpatients and outpatients as indicated by the first two principal coordinates analysis (PCoA) vectors (Figs. and ) and the Microbiome Regression-Based Kernel Association Test (MiRKAT) models (both Bray Crutis and Jaccard P < 0.001; Tab. ). Significant difference was also found between patients with disease-related symptoms, including fever (Bray Crutis P = 0.022; Jaccard P = 0.011), congestion (Bray Crutis P = 0.017; Jaccard P = 0.037), and oropharyngeal hyperaemia (Bray Crutis P = 0.070; Jaccard P = 0.021; Table ), compared to those without. The overall beta diversity also differed across infected virus groups (Figs. and ). For the analyses within each virus-infected group, the beta diversity was significant different between HRV-infected inpatients and HRV-infected outpatients (Bray Crutis P = 0.040; Jaccard P = 0.089), as well as between HMPV infected inpatients and HMPV infected outpatients but only observed in Jaccard matrix (Bray Crutis P = 0.241; Jaccard P = 0.008; Figs. – and Tab. ). As measured by the distance-based coefficient of determination ( R 2 ) using permutational multivariate analysis of variance (PERMANOVA) test, disease severity (inpatient or outpatient) explained most of the overall microbial community composition (i.e., beta diversity), followed by several disease-related symptoms including oropharyngeal hyperaemia, fever, and congestion. Most respiratory virus infection explained less than 1%, except HAdV and HMPV significantly correlated with Bray-Curtis matrix and HRV significantly correlated with Jaccard (Fig. ). 

The relative abundance of oropharyngeal phyla differed by disease severity and infected virus (Fig. ). The most abundance phyla in our study population were Proteobacteria (42.4% for inpatient, 59.4% for outpatient), Firmicutes (30.2%, 26.7%), Actinobacteria (17.6%, 12.5%), Bacteroidetes (7.6%, 0.9%), and Fusobacteria (1.3%, 0%). The presence of Proteobacteria was significantly associated with disease severity. The odds of being inpatients were 0.27 (95% CI = 0.10–0.73, P = 0.010) with the presence of Proteobacteria (Fig. ). The relative abundance of bacterial phyla was associated with respiratory virus infections. For every increase in one percent of the relative abundance, Proteobacteria (OR = 7.62, 95% CI = 1.92–30.32, P = 0.004) increased, and Firmicutes (OR = 0.13, 95% CI = 0.02–0.78, P = 0.026) and Actinobacteria (OR = 0.03, 95% CI = 0.001–0.62, P = 0.024) decreased, the odds of HMPV infection; Bacteroidetes (OR = 9.91, 95% CI = 1.44–68.36, P = 0.020) increased the odds of HAdV infection; Actinobacteria increased the odds of HPIV infection (OR = 6.59, 95% CI = 1.28–33.89, P = 0.024) and decreased the odds of HRV infection (OR = 0.01, 95% CI = 0.0001–0.90, P = 0.045; Fig. ). The 15 most abundance species in our study population were Stenotrophomonas maltophilia (14.4% for inpatient, 9.9% for outpatient), Rothia mucilaginosa (14.0%, 8.0%), Brucella anthropi (0.9%, 17.3%), Achromobacter aegrifaciens (4.6%, 11.1%), Achromobacter marplatensis (2.8%, 8.7%) Streptococcus infantis (5.1%, 3.2%), Streptococcus parasanguinis (2.4%, 5.2%), Streptococcus mitis (2.1%, 3.9%), Streptococcus salivarius (2.6%, 3.5%), Pseudomonas oleovorans (2.3%, 2.9%), Neisseria subflava (3.7%, 0.9%), Phyllobacterium myrsinacearum (4.1%, 0%), Dolosigranulum pigrum (2.0%, 1.6%), Cutibacterium acnes (3.2%, 0.1%), and Haemophilus influenzae (2.1%, 0.9%; Fig. ). Among all the microbial species detected in our study population, sixteen species had more than twofold higher abundance in outpatients than inpatients (Fig. and Tab. ), where Enterobacter kobei , A. marplatensis , A. aegrifaciens , Pseudomonas putida , Achromobacter xylosoxidans , B. anthropi , Achromobacter pestifer , Ochrobactrum sp.CGA5 , and Achromobacter insolitus had an false discover rate (FDR) less than 0.05. Five species had higher abundance in inpatients compared to outpatients, but not statistically significant after FDR adjustment. When restricted this analyses within each virus infected group, A. marplatensis , A. xylosoxidans , and A. aegrifaciens among HMVP infection patients, and A. xylosoxidans and Delftia acidovorans among respiratory virus undetectable patients were twofold higher abundance in inpatients than in outpatients (Fig. ). However, no species were significantly different at an FDR of 5% for the infection group stratified analyses. According to the LDA score by linear discriminant analysis and effective size comparisons (LEfSe) analysis, fifteen taxa from order to species levels (order Burkholderiales and Enterobacterales ; family Brucellaceae , Alcaligenaceae , and Enterobacteriaceae ; genus Brucella , Achromobacter , and Enterobacter ; species B. anthropi , A. aegrifaciens , A. marplatensis , P. putida , A. xylosoxidans , E. kobei , and A. pestifer ) were enriched among outpatients, but no taxa were found to be enriched among inpatients (Fig. , ). When analyzed the taxonomic abundance of disease-related symptoms, statistical significance was only found between patients with fatigue and those without, with six taxa enriched in patients with fatigue at a FDR of 5% (Table ). 

Dirichlet multinomial mixture (DMM) model was used to unsupervised clustering the sample to three community types, with 60 samples assigned to community type 1 (CT1), 27 samples assigned to community type 2 (CT2), and 22 samples assigned to community type 3 (CT3). The alpha diversity of CT1 and CT2 was significantly lower than CT3 (Fig. ), and the overall microbial composition was significant different between the three groups (Fig. , ). Among the top driver species (Fig. ), the relative abundance of R. mucilaginosa was the highest within each community, but was higher in CT1 and CT3 than in CT2 when compared across different communities. S. maltophilia and A. aegrifaciens were enriched in CT1; A. xylosoxidans and E. kobei were enriched in CT2; and S. parasanguinis and S. salivarius were enriched in CT3 (Fig. ). The distribution of CT1 and CT3 was similar between inpatients and outpatients, whereas the proportion of CT2 was higher in outpatients than in inpatients (Fig. ). The prevalence of CT1 was higher among patients with IFV, EV, and HAdV infections; CT2 was higher among patients with HMPV infections; and CT3 was higher among patients with RSV and HAdV infections. Participants of CT1 and CT2 were more likely to have cough (81.7% of CT1, 85.2% of CT2), and participants with CT3 were more likely to have fever (77.3%). 

A total of 431 microbial pathways were identified, with 334 pathways present in 5–95% of the population. Among them, 244 pathways had more than twofold higher abundance in outpatients compared to inpatients, and 231 were statistically significant at an FDR of 5% (Fig. and Table ). The quantile-quantile plot of P value showed that most of the pathways had a higher observed P value than the expected P value (Fig. ). Many of these pathways are linked to metabolic degradation, biosynthesis, or cellular respiration. For example, the top 25 most distinctive enriched pathways in outpatients include 13 degradation pathways 14DICHLORBENZDEG-PWY (1,4-dichlorobenzene degradation), CATECHOL-ORTHO-CLEAVAGE-PWY (catechol degradation to β-ketoadipate), GALACTARDEG-PWY (D-galactarate degradation I), GLUCARGALACTSUPER-PWY (superpathway of D-glucarate and D-galactarate degradation), PWY-1361 (benzoyl-CoA degradation I), PWY-5417 (catechol degradation III), PWY-5431 (aromatic compounds degradation via β-ketoadipate), PWY-561 (superpathway of glyoxylate cycle and fatty acid degradation), PWY-6084 (3,5-dichlorocatechol degradation), PWY-6182 (superpathway of salicylate degradation), PWY-6185 (4-methylcatechol degradation), PWY-8187 ( l -arginine degradation XIII), and TYRFUMCAT-PWY ( l -tyrosine degradation I), four biosynthesis pathways ARGSYNBSUB-PWY ( l -arginine biosynthesis II), PWY-5855 (ubiquinol-7 biosynthesis), PWY-6292 (superpathway of l -cysteine biosynthesis), and PWY-7400 ( l -arginine biosynthesis IV), as well as a six other metabolic pathways P105-PWY (TCA cycle IV), P124-PWY (Bifidobacterium shunt), PWY-5138 (fatty acid β-oxidation IV), PWY-5676 (acetyl-CoA fermentation to butanoate II), PWY-6969 (TCA cycle V), and PWY66-391 (fatty acid β-oxidation VI; Fig. ). When analyzed the pathway abundance of disease-related symptoms, no significant difference was found at an FDR of 5% (Fig. ). 

As we had only one HPIV-infected inpatient with a valid sequencing result, the asymmetry of inpatients and outpatients in this group may bias our results. Therefore, we conducted a sensitivity analysis that excluded all HPIV-infected patients. Results were similar when this group of patients was excluded (Fig. ). 

This study described the microbiome of the oropharynx during acute respiratory infection caused by a variety type of viruses in infants and children. We found that the overall taxonomic alpha diversity differed by the types of infected virus, and the beta diversity differed by disease severity, disease-related symptoms, and the types of infected virus. Nine species had significantly higher abundance in outpatients than in inpatients, with five of them in the genus Achromobacter . Three microbial community types were identified. The prevalence of CT1 was higher among patients with IFV, EV, and HAdV; CT2 was higher among patients with HMPV; and CT3 was higher among patients with RSV and HAdV infections. Participants of CT1 and CT2 were more likely to have cough, and CT3 were more likely to have fever. Changes in the upper respiratory tract microbiome during virus infection have been reported in several previous studies. One study conducted in the US described the nasopharyngeal microbial profile in infants with ARI due to HRV and RSV, and they found that the overall taxonomic composition was significantly different between infants infected with HRV and those infected with RSV . However, in their study, HRV ARI was found to have higher Shannon index compared to RSV ARI, which is inconsistent with our study that the Shannon index is higher in the RSV group. The difference may be explained by the sample types, as they collected the nasopharyngeal samples and HRV is more likely to colonize in the nasal mucosa , while the samples collected in our study were from the oropharynx. Other studies compared the upper respiratory tract microbiome to healthy individuals. A study from southern China characterized the oropharyngeal microbiome among outpatients with IFV A, IFV B, HRV, and RSV infections, as well as healthy controls . Results showed that IFV A, IFV B, and HRV, but not RSV, exhibited significant changes in the microbial composition compared to the healthy controls. However, different microbial composition between healthy and RSV-infected individuals was observed in a study of infant nasal microbiome in the US , but the result of IFV was confirmed by another US study that compared the nasopharyngeal microbiome among individuals infected with IFV A/H3N2, IFV B, and healthy controls . Although no significant difference in the Shannon index was detected between the influenza and control groups in that US study, the difference of microbial compositions was statistically significant. The core bacteria in our study were not exactly the same as the others. The top five most abundance phyla of our study were Proteobacteria , Firmicutes , Actinobacteria , Bacteroidetes , and Fusobacteria . In a study from UK described the oropharyngeal microbiome among adults with or without IFV A infection, the top five most abundance phyla were the same, but the order of abundance in descending was Firmicutes , Fusobacteria , Proteobacteria , Bacteroidetes , and Actinobacteria . High abundance of Streptococcus , Moraxella , Corynebacterium , and Haemophilus have been found in nasopharyngeal samples in previous studies , . Specifically, Haemophilus influenza is a bacterium that commonly found in ARI . In the US study of infant RSV infection and nasal microbiome, H. influenzae was associated with disease severity . Although this association was not found in our study, we did observe very different proportion of Haemophilus influenza between inpatients (3.2%) and outpatients (0.1%). Our study observed five species in genus Achromobacter have higher abundance in outpatients compared to inpatients. Achromobacter is an obligately aerobic Gram-negative bacterium and can be isolated from different kinds of environmental habitats. Certain species such as A. xylosoxidans has been reported to be able to colonize in the human respiratory tract and cause disease in certain populations such as individuals with cystic fibrosis . However, limited evidence is available for the role of Achromobacter in acute respiratory infections, and the underlying mechanism for the different abundance between inpatients and outpatients needs additional research. In the analyses restricted within each infected group, D. acidovorans were found to be twofold higher abundance in inpatients than in outpatients among respiratory virus undetectable patients. Delftia was reported as a potential reagent and laboratory contaminant genera in a previous study . However, none of our negative controls contained this species, but it appeared in 19.3% (21/109) of the real samples. There is no clue that Delftia is a contaminant in the current study. In the analysis of community types, we found CT1 and CT2 was driven mostly by Achromobacter spp ., and CT1 and CT3 was driven mostly by Streptococcus spp . Streptococcus includes Gram-positive bacteria that are commensals and opportunistic pathogens for humans . More than 60 Streptococcus species have been recognized so far. Among them, Streptococcus pneumoniae is considered as one of the most pathogenic types of Streptococcus for humans and often finds in ARI, but we did not see many Streptococcus pneumoniae infected cases in our study. The Streptococcus species enriched in CT1 and CT3, including S. infantis , S. mitis , S. parasanguinis , and S. salivarius , were defined as viridans, which have been reported to be associated with dental caries, arthritis, and infective endocarditis but may not be implicated in ARI development . Pathways related to l -arginine metabolism, including ARGSYNBSUB-PWY ( l -arginine biosynthesis II), PWY-7400 (L-arginine biosynthesis IV), and PWY-8187 (L-arginine degradation XIII), have the highest geometric mean abundance among outpatients, but were significantly lower among inpatients. L-arginine is the substrate NOS uses to generate NO, which can directly lead to inhibition of virus and indirectly regulate inflammation and immune responses . In a study of COVID-19, adults and children infected SARS-CoV-2 had significantly lower levels of L-arginine compared to controls . In another study, levels of L-arginine were shown to inversely be associated with the severity of COVID-19, which was consistent with our findings that the bioactivity of l -arginine had an impact of disease severity . Although many studies attempting to characterize the microbial profiles in the upper respiratory tract among individuals infected with IFV, RSV, and HRV, to our knowledge, this is the first study to describe the oropharyngeal microbiome of infants with HMPV, HPIV, HAdV, or EV ARI by shotgun sequencing. Limitations of our study should also be noted. First, this study is cross-sectional, and thus we were unable to evaluate the dynamic change of microbiome during ARI. We cannot determine whether the difference of microbiome was due to the infection, or the microbiome predisposed to infection. Second, the sample size of our study is relatively small and therefore limited our statistical power to identify small differences in associations, especially for the analysis within each virus group. Larger study is needed in the further. Third, the distribution of the reads is relatively wide, which may have an impact on the microbial composition. Additionally, although patients taken any antiviral drugs before visiting the hospital were excluded, information of antibiotics or steroids use was not collected, which might have potentially confounded the results. However, the oral microbiome has been found to recover quickly and was robust against antibiotic-induced disturbances . Finally, our study restricted to participants with single virus infection, but virus co-infections in ARI were very common , especially for HAdV and EV . Our findings might not be generalizable to ARI patients with more than one virus infections. In conclusion, we found associations of oropharyngeal microbiome with ARI disease severity, disease-related symptoms, and the types of infected virus. Further research to understand the dynamic change of microbiome during ARI and how these affect the development of the immune response in early life is warranted. 

Study population This study included participants from six hospitals in Chongqing Municipality, China, between October 2021 and April 2022. Participants were inpatients or outpatients with ARI, which were defined as a sudden onset of at least one of the following respiratory symptoms: cough, sore throat, shortness of breath, or nasal congestion. Outpatients were enrolled from departments of pediatrics, fever clinic, respiratory, and infections of the selected hospital. Inpatients were those who met the ARIs case definition and required hospitalized for more than 24 h, as determined by the doctors based on the results of auscultation, blood routine examination, blood gases and pulse oxygen saturations, C-reactive protein, and CT scans. Inpatients were hospitalization for intravenous infusion and nebulized medication treatment, and got continuous monitoring of their conditions. None of the inpatients in this study required supplemental oxygen treatment. Hospitalization was used as measure of disease severity. The infants and children required for hospitalized in this cohort were healthy term-born and did not experience vomiting or diarrhea. Subjects who have taken any antiviral drugs before visiting the hospital were excluded. A total of 3373 patients were recruited during this period. Demographic and disease-related symptom data were collected by a standardized interview-based questionnaire. Oropharyngeal swabs were collected prior to therapeutic measures, sorted under 2–8 °C immediately, and shipped to the local Center of Disease Prevention and Control daily. All samples were aliquoted into two sterile cryovials. The first aliquot was used to test for IFV within one week for surveillance purpose required by China CDC according to their standard operating protocol . The remaining samples were aliquoted into two sterile cryovials and shipped to the central laboratory in Beijing by dry ice. All cryovials were stored at −70°C until further tested. All ethical regulations relevant to human research participants were followed. PCR was conducted for all samples to detect EV, HAdV, HMPV, HRV, IFV, HPIV, RSV, as well as human bocavirus and coronavirus (described in detail below), and 1212 samples with PCR-confirmed single virus infection was selected as the sampling frame. R package “matchit” was used to match inpatients and outpatients with age and sex as covariates, in a ratio of one to one. Eight groups were generated: seven groups were single infection of EV, HAdV, HMPV, HRV, IFV, HPIV, or RSV, respectively, and the last group was individuals without any of these viruses as well as human bocavirus and coronavirus infection (i.e., virus undetectable group). After that, infants and children aged less than 12 years were selected, and “runif” function was used to generate random numbers for all matched pairs. For the groups with more than ten inpatients (EV, HAdV, HMPV, HPIV), ten pairs were selected from small to large random numbers; for groups with less than ten inpatients (IFV N = 8, HRV N = 8, RSV N = 9), all pairs were selected. To fill up the batch, 12 pairs of the virus undetectable group were selected from small to large random numbers. Finally, 154 participants were selected for microbiome detection. Laboratory procedures The remaining samples of the first aliquot were used for respiratory virus detection. Viral nucleic acid was extracted using the nucleic acid extraction and purification kit (Guangzhou Bay Area Biotechnology Co., Ltd,) according to the manufacturer instructions. In brief, 300 μl sample was put into a deep-well multi-well plate, and then 20 μl proteinase K as well as 20 μl nucleic acid protection solution added to the same plate. The materials were then transferred to the automated nucleic acid extraction instrument for extraction. This procedure was completed in 30 min. Then a multiplex real-time quantitative PCR assay (Biomax Biological technology Co, Ltd) was performed to detect EV, HAdV, HMPV, HRV, HPIV, and RSV. Each PCR reaction system included 5 μl template and amplified for 45 cycles, with a cycle threshold value less than 40 to be defined as positive. The second aliquot of the oropharyngeal sample was used for whole-genome shotgun metagenomic sequencing. The extraction batches were created by randomly selecting study participants to have an even distribution of both inpatients and outpatients and of infected viruses within each batch. The sample was thawed and 300 μl of material was transferred into a new tube. After transferring, 10 µl IC, 150 µl NPB, and 7 µl HLE were added to the tube and thoroughly mixed. After incubated for 15 min at 35 °C, 20 µl PK was added and thoroughly mixed. Samples were then incubated for another 10 minutes at 37 °C. The processed samples were transferred to the PathoXtract TM system (Willingmed Corp.) for automatic DNA extraction and library construction according to the manufacturer instructions. DNA quantity was evaluated using Qubit dsDNA HS (Thermo Fisher). One positive and one negative control were included in each run. Sequencing of the libraries was performed on NextSeq CN550 (Illumina) using single-end 75 bp reads and at least 20 million sequencing reads were acquired for each sample. Bioinformatics Thirty-three samples had no valid sequencing data. The raw shotgun sequencing reads of the remaining 121 samples were quality checked, quality trimmed, and filtered by KneadData (v0.12.0, https://huttenhower.sph.harvard.edu/kneaddata ), which integrates the tools FastQC, Trimmomatic, and Bowtie2 . Host filtering was performed by aligning the reads to a human reference genome (GRCh38). The non-human reads were then processed based on taxonomic profiling using MetaPhlAn4 (v4.0.6) and functional profiling using HUMAnN3 (v3.8) . The functional information was obtained by aligning the sequences based on the UniRef90 diamond full reference database (version 201901b) . Relative abundances of taxa from phylum to species level were created. To reduce bias for relative abundance calculation, samples with less than 1000 reads were removed from downstream analysis, and each sample was rarefied to 1000 microbial sequences for taxonomic alpha and beta diversity analyses based on the rarefaction curves for alpha diversity (Fig. ). After rarefaction, a total of 111 samples were retained. Potential contaminants were identified using the ‘decontam’ package. C. acnes was identified as a contaminant and removed from the taxonomic table. After removal, two more samples did not reach the 1000 reads threshold and excluded for further analyses. Finally, 109 samples were included in this study. Distribution of reads is shown in Fig. . Alpha diversity measures, including observed features and the Shannon index, and beta diversity matrices, including Bray–Curtis and Jaccard, were calculated using QIIME 2 (v2023.5) . Principal coordinates analysis (PCoA) was calculated from the two beta-diversity distance matrices. For the quality control samples, the positive controls appeared to cluster similarly based on visual inspection of the PCoA plots of the two beta diversity measures (Fig. ). In the negative controls, only one negative control sample had 1065 reads and other negative controls did not find any reads after decontamination. Statistical and reproducibility Characteristics of study participants were summarized by disease severity. Differences between categorical groups in alpha diversity were determined by Kruskal–Wallis test and Wilcoxon rank-sum test. PCoA plots were generated using the first five PCoA vectors, which explained a total of 22.7% and 22.8% of the variability in Bray Curtis and Jaccard, respectively. The difference of the overall beta-diversity distance matrices by disease severity (inpatient, outpatient), disease-related symptoms (yes, no), or infected virus (positive, netative) was tested using Microbiome Regression-Based Kernel Association Test (MiRKAT) . The proportion of variance in the microbiome that can be explained by host and virus factors was calculated by the distance-based coefficient of determination ( R 2 ) using the PERMANOVA test, and the 95% CI was calculated using 1000 bootstrap samples. The relative abundances of microbial taxa were arcsine square root transformed for analyses. Logistic regression models were used to calculate prevalence odds ratios (ORs) and 95% confidence intervals (CI) for the association of taxa presence/absence with disease severity and taxa relative abundance with infected virus. The fold changes of the microbial species and pathways between outpatients and inpatients as well as between patients with and without specific disease-related symptom were calculated by log 2 -transformed of the ratio of feature (taxa or pathways) relative abundance multiplying prevalence in outpatient to inpatients as: \documentclass[12pt]{minimal}
				\usepackage{amsmath}
				\usepackage{wasysym} 
				\usepackage{amsfonts} 
				\usepackage{amssymb} 
				\usepackage{amsbsy}
				\usepackage{mathrsfs}
				\usepackage{upgreek}
				\setlength{\oddsidemargin}{-69pt}
				\begin{document}$${{Fold}-{change}}_{i}={{{\mathrm{log}}}}_{2}\left(\frac{{{GeoM}}_{{Out},i}\times {{Prev}}_{{Out},i}}{{{GeoM}}_{{In},i}\times {{Prev}}_{{In},i}}\right)$$\end{document} F o l d − c h a n g e i = log 2 G e o M O u t , i × P r e v O u t , i G e o M I n , i × P r e v I n , i where GeoM Out ,i and GeoM In ,i represent the geometric mean of feature i in the subset of outpatients or inpatients where the feature is present, and Prev Out ,i and Prev In ,i represent the prevalence of feature i among this subset of patients. A stratified analysis was then performed to assess the fold change of the taxa in outpatients versus inpatients by infected virus. Signature taxonomic groups were also identified by linear discriminant analysis and effective size comparisons (LEfSe) . Microbial community types were defined by the Dirichlet multinomial mixture (DMM) model . The number of community types was determined by identifying the point of the Laplace approximation where Dirichlet components increase led to minor reductions in model fit (Fig. ). For the functional pathways analyses, we restricted to pathways present in 5% to 95% of the population. Quantile-quantile plots were generated by two-sided Wilcoxon rank-sum P values for all genes and pathways. Reporting summary Further information on research design is available in the linked to this article. 

This study included participants from six hospitals in Chongqing Municipality, China, between October 2021 and April 2022. Participants were inpatients or outpatients with ARI, which were defined as a sudden onset of at least one of the following respiratory symptoms: cough, sore throat, shortness of breath, or nasal congestion. Outpatients were enrolled from departments of pediatrics, fever clinic, respiratory, and infections of the selected hospital. Inpatients were those who met the ARIs case definition and required hospitalized for more than 24 h, as determined by the doctors based on the results of auscultation, blood routine examination, blood gases and pulse oxygen saturations, C-reactive protein, and CT scans. Inpatients were hospitalization for intravenous infusion and nebulized medication treatment, and got continuous monitoring of their conditions. None of the inpatients in this study required supplemental oxygen treatment. Hospitalization was used as measure of disease severity. The infants and children required for hospitalized in this cohort were healthy term-born and did not experience vomiting or diarrhea. Subjects who have taken any antiviral drugs before visiting the hospital were excluded. A total of 3373 patients were recruited during this period. Demographic and disease-related symptom data were collected by a standardized interview-based questionnaire. Oropharyngeal swabs were collected prior to therapeutic measures, sorted under 2–8 °C immediately, and shipped to the local Center of Disease Prevention and Control daily. All samples were aliquoted into two sterile cryovials. The first aliquot was used to test for IFV within one week for surveillance purpose required by China CDC according to their standard operating protocol . The remaining samples were aliquoted into two sterile cryovials and shipped to the central laboratory in Beijing by dry ice. All cryovials were stored at −70°C until further tested. All ethical regulations relevant to human research participants were followed. PCR was conducted for all samples to detect EV, HAdV, HMPV, HRV, IFV, HPIV, RSV, as well as human bocavirus and coronavirus (described in detail below), and 1212 samples with PCR-confirmed single virus infection was selected as the sampling frame. R package “matchit” was used to match inpatients and outpatients with age and sex as covariates, in a ratio of one to one. Eight groups were generated: seven groups were single infection of EV, HAdV, HMPV, HRV, IFV, HPIV, or RSV, respectively, and the last group was individuals without any of these viruses as well as human bocavirus and coronavirus infection (i.e., virus undetectable group). After that, infants and children aged less than 12 years were selected, and “runif” function was used to generate random numbers for all matched pairs. For the groups with more than ten inpatients (EV, HAdV, HMPV, HPIV), ten pairs were selected from small to large random numbers; for groups with less than ten inpatients (IFV N = 8, HRV N = 8, RSV N = 9), all pairs were selected. To fill up the batch, 12 pairs of the virus undetectable group were selected from small to large random numbers. Finally, 154 participants were selected for microbiome detection. 

The remaining samples of the first aliquot were used for respiratory virus detection. Viral nucleic acid was extracted using the nucleic acid extraction and purification kit (Guangzhou Bay Area Biotechnology Co., Ltd,) according to the manufacturer instructions. In brief, 300 μl sample was put into a deep-well multi-well plate, and then 20 μl proteinase K as well as 20 μl nucleic acid protection solution added to the same plate. The materials were then transferred to the automated nucleic acid extraction instrument for extraction. This procedure was completed in 30 min. Then a multiplex real-time quantitative PCR assay (Biomax Biological technology Co, Ltd) was performed to detect EV, HAdV, HMPV, HRV, HPIV, and RSV. Each PCR reaction system included 5 μl template and amplified for 45 cycles, with a cycle threshold value less than 40 to be defined as positive. The second aliquot of the oropharyngeal sample was used for whole-genome shotgun metagenomic sequencing. The extraction batches were created by randomly selecting study participants to have an even distribution of both inpatients and outpatients and of infected viruses within each batch. The sample was thawed and 300 μl of material was transferred into a new tube. After transferring, 10 µl IC, 150 µl NPB, and 7 µl HLE were added to the tube and thoroughly mixed. After incubated for 15 min at 35 °C, 20 µl PK was added and thoroughly mixed. Samples were then incubated for another 10 minutes at 37 °C. The processed samples were transferred to the PathoXtract TM system (Willingmed Corp.) for automatic DNA extraction and library construction according to the manufacturer instructions. DNA quantity was evaluated using Qubit dsDNA HS (Thermo Fisher). One positive and one negative control were included in each run. Sequencing of the libraries was performed on NextSeq CN550 (Illumina) using single-end 75 bp reads and at least 20 million sequencing reads were acquired for each sample. 

Thirty-three samples had no valid sequencing data. The raw shotgun sequencing reads of the remaining 121 samples were quality checked, quality trimmed, and filtered by KneadData (v0.12.0, https://huttenhower.sph.harvard.edu/kneaddata ), which integrates the tools FastQC, Trimmomatic, and Bowtie2 . Host filtering was performed by aligning the reads to a human reference genome (GRCh38). The non-human reads were then processed based on taxonomic profiling using MetaPhlAn4 (v4.0.6) and functional profiling using HUMAnN3 (v3.8) . The functional information was obtained by aligning the sequences based on the UniRef90 diamond full reference database (version 201901b) . Relative abundances of taxa from phylum to species level were created. To reduce bias for relative abundance calculation, samples with less than 1000 reads were removed from downstream analysis, and each sample was rarefied to 1000 microbial sequences for taxonomic alpha and beta diversity analyses based on the rarefaction curves for alpha diversity (Fig. ). After rarefaction, a total of 111 samples were retained. Potential contaminants were identified using the ‘decontam’ package. C. acnes was identified as a contaminant and removed from the taxonomic table. After removal, two more samples did not reach the 1000 reads threshold and excluded for further analyses. Finally, 109 samples were included in this study. Distribution of reads is shown in Fig. . Alpha diversity measures, including observed features and the Shannon index, and beta diversity matrices, including Bray–Curtis and Jaccard, were calculated using QIIME 2 (v2023.5) . Principal coordinates analysis (PCoA) was calculated from the two beta-diversity distance matrices. For the quality control samples, the positive controls appeared to cluster similarly based on visual inspection of the PCoA plots of the two beta diversity measures (Fig. ). In the negative controls, only one negative control sample had 1065 reads and other negative controls did not find any reads after decontamination. 

Characteristics of study participants were summarized by disease severity. Differences between categorical groups in alpha diversity were determined by Kruskal–Wallis test and Wilcoxon rank-sum test. PCoA plots were generated using the first five PCoA vectors, which explained a total of 22.7% and 22.8% of the variability in Bray Curtis and Jaccard, respectively. The difference of the overall beta-diversity distance matrices by disease severity (inpatient, outpatient), disease-related symptoms (yes, no), or infected virus (positive, netative) was tested using Microbiome Regression-Based Kernel Association Test (MiRKAT) . The proportion of variance in the microbiome that can be explained by host and virus factors was calculated by the distance-based coefficient of determination ( R 2 ) using the PERMANOVA test, and the 95% CI was calculated using 1000 bootstrap samples. The relative abundances of microbial taxa were arcsine square root transformed for analyses. Logistic regression models were used to calculate prevalence odds ratios (ORs) and 95% confidence intervals (CI) for the association of taxa presence/absence with disease severity and taxa relative abundance with infected virus. The fold changes of the microbial species and pathways between outpatients and inpatients as well as between patients with and without specific disease-related symptom were calculated by log 2 -transformed of the ratio of feature (taxa or pathways) relative abundance multiplying prevalence in outpatient to inpatients as: \documentclass[12pt]{minimal}
				\usepackage{amsmath}
				\usepackage{wasysym} 
				\usepackage{amsfonts} 
				\usepackage{amssymb} 
				\usepackage{amsbsy}
				\usepackage{mathrsfs}
				\usepackage{upgreek}
				\setlength{\oddsidemargin}{-69pt}
				\begin{document}$${{Fold}-{change}}_{i}={{{\mathrm{log}}}}_{2}\left(\frac{{{GeoM}}_{{Out},i}\times {{Prev}}_{{Out},i}}{{{GeoM}}_{{In},i}\times {{Prev}}_{{In},i}}\right)$$\end{document} F o l d − c h a n g e i = log 2 G e o M O u t , i × P r e v O u t , i G e o M I n , i × P r e v I n , i where GeoM Out ,i and GeoM In ,i represent the geometric mean of feature i in the subset of outpatients or inpatients where the feature is present, and Prev Out ,i and Prev In ,i represent the prevalence of feature i among this subset of patients. A stratified analysis was then performed to assess the fold change of the taxa in outpatients versus inpatients by infected virus. Signature taxonomic groups were also identified by linear discriminant analysis and effective size comparisons (LEfSe) . Microbial community types were defined by the Dirichlet multinomial mixture (DMM) model . The number of community types was determined by identifying the point of the Laplace approximation where Dirichlet components increase led to minor reductions in model fit (Fig. ). For the functional pathways analyses, we restricted to pathways present in 5% to 95% of the population. Quantile-quantile plots were generated by two-sided Wilcoxon rank-sum P values for all genes and pathways. 

Further information on research design is available in the linked to this article. 

Supplementary information Description of Additional Supplementary File Supplementary Data Reporting summary Transparent Peer Review file 

