Nasopharyngeal microbiota in hospitalized children with Bordetella pertussis and Rhinovirus infection

Despite great advances in describing Bordetella pertussis infection, the role of the host microbiota in pertussis pathogenesis remains unexplored. Indeed, the microbiota plays important role in defending against bacterial and viral respiratory infections. We investigated the nasopharyngeal microbiota in infants infected by B. pertussis (Bp), Rhinovirus (Rv) and simultaneously by both infectious agents (Bp + Rv). We demonstrated a specific nasopharyngeal microbiome profiles for Bp group, compared to Rv and Bp + Rv groups, and a reduction of microbial richness during coinfection compared to the single infections. The comparison amongst the three groups showed the increase of Alcaligenaceae and Achromobacter in Bp and Moraxellaceae and Moraxella in Rv group. Furthermore, correlation analysis between patients’ features and nasopharyngeal microbiota profile highlighted a link between delivery and feeding modality, antibiotic administration and B. pertussis infection. A model classification demonstrated a microbiota fingerprinting specific of Bp and Rv infections. In conclusion, external factors since the first moments of life contribute to the alteration of nasopharyngeal microbiota, indeed increasing the susceptibility of the host to the pathogens' infections. When the infection is triggered, the presence of infectious agents modifies the microbiota favoring the overgrowth of commensal bacteria that turn in pathobionts, hence contributing to the disease severity.

Modulation of the nasopharyngeal microbiota in the presence of B. pertussis infection has been described in the mouse models only 15,16 . These studies evidenced a preventive action of the nasopharyngeal resident microbiota against B. pertussis colonization 15 , while the use of antibiotics produces the opposite effect, favoring B. pertussis colonization 16 .
Understanding the modulation of ecology of nasopharynx in humans with B. pertussis infection may help to understand determinants and potential diagnostic, preventive and therapeutic implications for this disease. In addition, no data are available on the changes that may derive from the simultaneous infection by B. pertussis and Rhinovirus.
In this study, we investigated the characteristics of nasopharyngeal microbiome in infants with B. pertussis and Rhinovirus simultaneous infection with those of either infection alone.

Results
Patients' clinical features. Socio-demographic characteristics of patients included in the study are shown in Table 1.
The majority of patients were males (57.4%) with a median age of 1.5 months; patients with B. pertussis infection were slightly older than those with Rhinovirus infection or coinfection.
Children with coinfection were more likely breastfed compared to those with Rhinovirus even if the difference is small (57% vs 50% respectively, p = 0.029). On the other hand, children with B. pertussis infection were more likely formula-fed compared to those with coinfection (64% vs 14%, p = 0.029). Table 2 shows the clinical features among the three groups. Cough, paroxysmal cough, whooping and apnea were more frequent in patients with B. pertussis infection. Cyanosis was occurred more frequently in the coinfection group by comparison to the Bp or Rv groups. In addition, the coinfection group more likely uses antibiotics before admission.
Dendrogram analysis, based on clinical features' correlation, revealed that C-reactive protein (CRP), fever and antibiotic treatment at the admission (others' than macrolides) were correlated together and constituted a separated cluster from the other clinical variables. The other features were subgrouped into: (i) apnea, cyanosis, macrolides' assumption and antibiotics treatment before the admission, and (ii) delivery, feeding cough, paroxysmal cough and inspiratory stridor, complications, vomit and leukocytes plus lymphocytes count (Fig. S1).
Nasopharyngeal microbiota profiling. A total of 2,393,694 sequencing reads were obtained from 54 nasopharyngeal aspirates, with a mean value of 44,328 sequences per sample. We identified an overall of 245 operational taxonomic units (OTUs), belonging to 15 phyla and 95 families.
To assess the overall differences of microbial community structures for the three groups, ecological parameters were based on alpha-diversity indexes (i.e., Shannon, Chao1, and Observed species indexes).
Despite no statistical differences were resulted, samples belonging to Rv group showed higher values of all indexes followed by Bp samples and then Bp + Rv samples (Fig. 1, panels a-c).
Regarding beta diversity analysis, only Bray Curtis algorithm evidenced a statistically significant separation amongst groups (PERMANOVA p value ≤ 0.05), while phylogenetic unweighted and weighted UniFrac algorithms not (PERMANOVA p values > 0.05) (Fig. 1, panels d-f). Analyses of the intra-group distances revealed www.nature.com/scientificreports/ statistically significant differences between Bp and Rv for Bray-Curtis and Unweighted UniFrac metrics ( Fig. 1, panels g-i).
The OTU distribution was investigated at the phylum and genus levels. At phylum level, Proteobacteria, Firmicutes, Actinobacteria and Bacteroidetes were the most representative phyla in the three groups (Fig. S2, Panel A). At genus level, Enterobacteriaceae, Veillonella, Staphylococcus, Gemellaceae, Alcaligenaceae and Achromobacter were increased in Bp ecosystem, while Streptococcus and Moraxella in Rv group (Fig. S2, Panel B).
Kruskal-Wallis test confirmed the higher abundance of Achromobacter and Alcaligenaceae in Bp group compared to the other groups (raw p value < 0.05) (Fig. 1, panels j and k), and the higher abundance of Moraxella and Moraxellaceae in Rv compared with others (raw p value < 0.05) (Fig. 1, panels l and m).
In the comparison amongst the three groups the principal component analysis (PCA) (Fig. S3, Panel a) and the partial least square (PLS) ( Correlation between clinical features and nasopharyngeal microbiota global distribution. At family level, Pearson's correlation revealed that all considered clinical features showed a correlation with a least one OTU. In particular, Alcaligenaceae was positively correlated with cough, paroxysmal cough and inspiratory stridor, while negatively with CRP. Moraxellaceae was negatively correlated with paroxysmal cough and inspiratory stridor (Fig. 2, panel a).
At genus level, antibiotic treatment before the admission was positively correlated with Paracoccus and Haemophilus and negatively with Rothia. Cough was negatively correlated with Clostridium, Blautia, Staphylococcus and Ruminococcus and positively with Achromobacter.
Paroxysmal cough was positively correlated with Neisseria and Achromobacter and negatively with Moraxella. Inspiratory stridor was positively correlated with Campylobacter and Achromobacter and negatively with Blautia, Enterococcus and Moraxella. Cyanosis was negatively correlated with Aggregatibacter. Apnea was positively correlated with Prevotella, Neisseria, Fusobacterium and Porphyromonas and negatively with Bacteroides, Oscillospira, Bifidobacterium and Rothia. CRP had positive correlation with Haemophilus and negative with Achromobacter. Finally, complications were positively related with Methanobrevibacter and Streptococcus (Fig. 2, panel b).

Comparison between Bp-and Rv-related microbiota and patients' clinical profiles.
To focus onto microbiota specific features associated to B. pertussis or Rhinovirus infection and to highlight their correlations with patients' clinical traits, we focalized on Bp-and Rv-related microbiota, excluding the coinfection group (Bp + Rv).
The Pearson's correlation performed on OTU profiles at genus level revealed two wide clusters. The first (a) was associated to Rv patients (9/16) with a nasopharyngeal microbiota profile rich in genera belonging to Actinobacteria, Bacteroidetes, Firmicutes and Proteobacteria, and related to low levels of leukocytes + lymphocytes (14/16) (Fig. 3). The second (b) was composed by two sub-clusters, namely b 1 and b 2 . The cluster b 1 was prevalently associated to Bp patients (6/7), characterized by high levels of leukocytes and lymphocytes, Caesarian section delivery and formula milk or mixed milk feeding and by a nasopharyngeal microbiota principally composed by Proteobacteria (Fig. 3). The patients in cluster b 2 , mostly Rv infected, was characterized by low level www.nature.com/scientificreports/ of leukocytes and lymphocytes, no antibiotic at the admission and not specific feeding and delivery modality and with a nasopharyngeal microbiota profile prevalently composed by Firmicutes and Proteobacteria (Fig. 3).

Model classifications analysis.
To investigate if nasopharyngeal microbiota profile could be predictive of a specific type of infection, we used the model classification analysis. This analysis, at genus level, revealed that the microbiota had capability to classify the 90% of the Bp-and Rv-related patients by the models Extra Trees Classifier, Gradient Boosting Classifier, Bagging Classifier and Decision Tree Classifier. A lower level of classification for Bp-(80%) was obtained compared to the classification of Rv-related patients (100%) (Table S1). Moreover, applying Random Forest Classifier and XGBRF Classifier, the microbiota profile reached the 100% of capability in correctly classifying both groups of patients. In Fig. 4 were reported the OTUs selected by the model analysis and the respective contribution value (importance score) in the model prediction.

Discussion
Our study describes for the first time the nasopharyngeal microbiota in children with pertussis.
We showed a different microbiota profiles consisting in the increase of Enterobacteriaceae, Veillonella, Staphylococcus Gemellaceae, Alcaligenaceae and Achromobacter in Bp group and of Streptococcus and Moraxella in Rv group.
The absence of comparison between B. pertussis infected patients with healthy subjects, due to the invasiveness of the nasopharyngeal aspirate in healthy infants, represents a limitation of this study. For this reason and in consideration of the high prevalence of cases of B. pertussis-Rhinovirus coinfection, we compared the nasopharyngeal microbiota profile of infants affected by single infection of B. pertussis or Rhinovirus with that of coinfected patients, who routinely underwent nasopharyngeal aspirate.
The nasopharyngeal microbiota of children is densely colonized by commensal bacteria, such as Moraxella, Haemophilus, Streptococcus, Flavobacterium, Dolosigranulum, Corynebacterium, Neisseria and Fusobacterium. However, Streptococcus pneumoniae, Haemophilus influenzae and Moraxella catarrhalis are considered potential airway pathogens 17 . When the nasopharyngeal eubiotic status is perturbed, these potentially pathogenic microorganisms can invade adjacent sites and cause different diseases [17][18][19][20][21][22] . It remains not clear which is the mechanism favoring pathogens' colonization; however the interaction amongst multiple pathogens including viral coinfections, the progressive shift from health to disbiotic status and the effect of antibiotic administration could represent triggering factors 17,20,23 .
Based on our results, a reduction of nasopharyngeal microbial richness was evident in coinfection related microbiota compared to that related to single infections. The reduction of the diversity of the nasopharyngeal microbiota has been actually associated with the increase frequency of upper respiratory infections 4,24 . Moreover, considering the comparison between coinfected and Rv groups, this effect could also depend to macrolide intake. In fact, 90.5% of coinfected patients took macrolides respect to only 16.7% of Rv. On the contrary, also  Analyzing the distance matrices and the composition of the microbiota, the co-infected microbiota profile was similar to both single infections' microbiota, even if with a higher similarity to the Rv group.
The comparison amongst the three groups revealed the increase of Alcaligenaceae and Achromobacter in the Bp infection group. Achromobacter spp., belonging to Alcaligenaceae family, are environmental microorganisms, causing a series of opportunistic illnesses, including respiratory diseases 25,26 . In our study, Achromobacter resulted related to cough, paroxysmal cough and inspiratory stridor, suggesting a role of this microorganism in the disease prognosis. Moreover, Moraxellaceae and Moraxella were increased in Rv infection group, as confirmed also by literature [27][28][29] . During Rv infection it has been described a colonization by as S. pneumoniae, H. influenzae and M. catarrhalis 3,13,30 . In particular, it seems that the nasopharyngeal colonization by Streptococcus and Moraxella precedes the Rv infection, probably due to the decreasing in ciliary function of epithelial cells and increasing of pro-inflammatory molecule release 30 . The simultaneous presence of Rv and Moraxella seems to be associated with severe lower respiratory infections, suggesting a correlation between specific viral/bacterial interactions and increase of disease severity 3 .
Evaluating the impact of delivery and feeding modality, count of leukocytes, administration of antibiotics it was possible to identify specific taxa as representative of Rv and Bp groups. Particularly, patients with vaginal delivery, maternal feeding, low count of leukocytes plus lymphocytes, low administration of antibiotics and rich microbiota profile resulted those infected by Rhinovirus. On the contrary, patients infected by B. pertussis were characterized by caesarian section delivery, formula milk feeding, high counts of leukocytes, antibiotic administration and high level of Proteobacteria. Thanks to these evidences we speculate that these conditions could shape the microbiota to be more prone to the airway infection by B. pertussis.
In conclusion, external factors since the first moments of life contribute to the alteration of nasopharyngeal microbiota, and this event could increase the susceptibility of the host to the pathogens' infections. In our opinion, when the infection is triggered, the presence of infectious agents could further modify the microbiota ecological niche favoring the overgrowth of other commensal bacteria that turn in pathobionts, hence contributing to the disease severity.

Methods
Study design and population. This study was performed at Bambino Gesù Children's Hospital (OPBG), of Rome, Italy. From January 2016 to December 2019, all infants younger than 1 year of age, accessing the emergency room with symptoms compatible with pertussis, according to European Centre for Disease Prevention and Control (ECDC) case definition, were routinely screened for pertussis and respiratory viruses in the context of an enhanced surveillance program. A subset of patients with either B. pertussis (Bp) or Rhinovirus (Rv) or both infections (Bp + Rv), for whom a leftover of the nasopharyngeal aspirate was available for the analysis of the pharynx microbiota, was selected for this study. The study was conducted in accordance with the Declaration www.nature.com/scientificreports/ of Helsinki, and the protocol was approved by Bambino Gesù Children's Hospital Ethics Committee (1355_ OPBG2017). Informed consent was obtained from the parents or legal guardians of all participants.
Data collection and laboratory confirmation of pertussis. Sociodemographic and clinical data were collected for each patient. Nasopharyngeal aspirate were collected within 24 h of hospital admission and processed immediately, or stored at − 80 °C until analysis. Two-hundred μl of nasopharyngeal sample were used for nucleic acids extraction by the EZ1 Virus Mini Kit v. 2.0 on the EZ1 Advanced XL platform (Qiagen, GmbH, Hilden, Germany). Bordetella Real Time (RT) PCR kits were used to investigate the presence of B. pertussis, amplifying the inter-space (IS) 481 target (Bordetella R-gene™ assay; Argene, Biomerieux, Marcy l'Etoile, France). To prevent misdiagnosis of Bordetella holmesii as B. pertussis, all samples positive for B. pertussis were confirmed by the amplification of ptxP (promoter of pertussis toxin) locus by a RT-PCR assay specific for B. pertussis. Moreover, culture for B. pertussis was also performed, using Regan-Lowe and Bordet-Gengou selective culture media at 37 °C under aerobic conditions. Nasopharyngeal aspirates were also processed by a multiplex RT-PCR using a specific panel detecting the 16   Biocomputational and statistical analyses. Illumina Miseq reads were checked for quality, length and chimera presence using the Qiime v1.9 pipeline 31 , and USEARCH algorithms 32 . Then, sequences were organized into Operational Taxonomic Units (OTUs) with a 97% of clustering threshold of pairwise identity by UCLUST with open reference method 32 . PyNAST v.0.1 software 31 was used to carry out a multiple sequence alignment against Greengenes 13_08 database with a 97% similarity for bacterial sequences and to build a phylogenetic tree 33 . Alpha diversity was performed by scikit-bio (http:// scikit-bio. org/) of Python 3 package using Shannon index, Chao1 estimator, and observed species indices, and the p value for group comparisons was determined by analysis of variance (ANOVA). Beta diversity analysis of samples, based on phylogenetically informed weighted and unweighted Unifrac distance matrices 34 and Bray Curtis matrices, was used and graphed by principal coordinate analysis (PCoA) plots. The association between the covariates and beta diversity measures was assessed by permutational analysis of variance (PERMANOVA) 35 .
Categorical variables were tabulated as frequencies and percentages and compared using Pearson's Chi square test.
OTUs were filtered to select only the features of interest by retaining OTUs with interquartile range (IQR) ≠ 0 and OTUs present in at least 25% of the samples.
The non-parametric Mann-Whitney U-test, Wilcoxon signed-rank, Kruskal-Wallis test, corrected for FDR p value ≤ 0.05, were used to compare OTUs relative abundance amongst groups.
Both the OTUs and clinical variables (considered both raw values and categorical variables) were clustered using the Pearson's correlation coefficient as the distance measure. In order to correlate clinical variables and OTUs, the Pearson's coefficient with relative p value was calculated and represented by heat map.
Multiple machine learning (ML) was trained for the classification tasks. The pipeline consisted of a tenfold cross-validation with a train-test split of 70-30%. To evaluate the model, the global and the single class accuracies were considered. The tested models were: Logistic Regression, SGD Classifier, Logistic Regression CV, Hist Gradient Boosting Classifier, Random Forest Classifier, Extra Trees Classifier, Gradient Boosting Classifier, Bagging Classifier, Ada Boost Classifier, XGB Classifier, XGBRF Classifier, MLP Classifier, Linear SVC, SVC, Gaussian NB, Decision Tree Classifier, Quadratic Discriminant Analysis, K Neighbors Classifier, Gaussian Process Classifier.