Variability and Diversity of Nasopharyngeal Microbiota in Children: A Metagenomic Analysis

The nasopharynx is the ecological niche for many commensal bacteria and for potential respiratory or invasive pathogens like Streptococcus pneumoniae, Haemophilus influenzae, and Neisseria meningitidis. Disturbance of a balanced nasopharyngeal (NP) microbiome might be involved in the onset of symptomatic infections with these pathogens, which occurs primarily in fall and winter. It is unknown whether seasonal infection patterns are associated with concomitant changes in NP microbiota. As young children are generally prone to respiratory and invasive infections, we characterized the NP microbiota of 96 healthy children by barcoded pyrosequencing of the V5–V6 hypervariable region of the 16S-rRNA gene, and compared microbiota composition between children sampled in winter/fall with children sampled in spring. The approximately 1000000 sequences generated represented 13 taxonomic phyla and approximately 250 species-level phyla types (OTUs). The 5 most predominant phyla were Proteobacteria (64%), Firmicutes (21%), Bacteroidetes (11%), Actinobacteria (3%) and Fusobacteria (1,4%) with Moraxella, Haemophilus, Streptococcus, Flavobacteria, Dolosigranulum, Corynebacterium and Neisseria as predominant genera. The inter-individual variability was that high that on OTU level a core microbiome could not be defined. Microbiota profiles varied strongly with season, with in fall/winter a predominance of Proteobacteria (relative abundance (% of all sequences): 75% versus 51% in spring) and Fusobacteria (absolute abundance (% of children): 14% versus 2% in spring), and in spring a predominance of Bacteroidetes (relative abundance: 19% versus 3% in fall/winter, absolute abundance: 91% versus 54% in fall/winter), and Firmicutes. The latter increase is mainly due to (Brevi)bacillus and Lactobacillus species (absolute abundance: 96% versus 10% in fall/winter) which are like Bacteroidetes species generally related to healthy ecosystems. The observed seasonal effects could not be attributed to recent antibiotics or viral co-infection. The NP microbiota of young children is highly diverse and appears different between seasons. These differences seem independent of antibiotic use or viral co-infection.


Introduction
According to the WHO, respiratory tract infections are still among the leading causes of death in children and adults worldwide [1]. The most common pathogens like Streptococcus pneumoniae, Haemophilus influenzae, Neisseria meningitidis and Staphylococcus aureus are normal and transient residents of the nasopharyngeal (NP) niche, where they are embedded in a complex microbiota of generally presumed harmless commensals. The human microbiome in general is assumed beneficial to the host due to stimulation and maturation of immune systems, promotion of mucosal structure and function and providing actual 'colonization resistance' against pathogen invasion [2]. Although colonization by the ''potential pathogens'' of the NP microbiome is mainly asymptomatic, progression towards upper respiratory tract infections, pneumonia or even sepsis and meningitis may occur [3,4]. The exact mechanisms by which this occurs remain largely unknown, although an imbalance in the composition of microbiota, for example by acquisition of new pathogens, viral coinfection or other host or environmental factors have been suggested [5][6][7][8]. In addition, clear correlations between invasive attack rates and season are observed for many of the potential pathogens of the upper respiratory tract [9,10], a phenomenon that cannot be fully explained by concomitant changes in colonization rates of the individual pathogenic bacteria [11,12]. This suggests that local containment of the colonizing pathogenic bacteria by the host and/or the surrounding ecosystem is of major importance in prevention of disease progression. Despite an abundance of data on incidence, prevalence and density of potential pathogens in NP microbiota of children and adults, the detailed composition of the NP microbial community, both during health and disease have not been studied. We, therefore, performed a meta-genomic study on the detailed composition of and variability in NP microbiota in young children sampled during different seasons.

Results and Discussion
We studied the NP microbiota composition of 96 healthy 18months old children. Their characteristics are depicted in Table  S1. Being aware of the current discussions on the artefacts that may be introduced by pyrosequencing [13,14], we applied a stringent protocol for filtering and clustering of sequences. The approx. 1 100 000 generated sequences (on average 11000 sequences per sample) yielded about 92 000 unique sequences, representing 13 taxonomic phyla and 243 species-level phyla types (OTUs). The data were normalized for equal numbers of reads per sample. The 5 most predominant phyla were Proteobacteria (64%), Firmicutes (21%), Bacteroidetes (11%), Actinobacteria (3%) and Fusobacteria (1.4%) ( Figure 1). In addition, we found representatives of Cyanobacteria, probably reflecting plant chloroplasts obtained through inhalation. Sporadically and/or in low abundance we found sequences for the candidate divisions OD1, TM7 and BRC1 and the phyla Deinococcus-Thermus, Nitrospira, Planctomycetes and Chloroflexi. On a lower taxonomic level, the most prevalent genera were Moraxella (40%), Haemophilus (20%), Streptococcus (12%), and Flavobacterium (10%). Other fairly common genera were Dolosigranulum (5%), Corynebacterium (2%), Neisseria (2%) and Fusobacterium (1%). The 30 most common OTUs representing almost 98% of all reads, and their relative and absolute presence are shown in Table 1 (For the complete list of OTUs; see Table S2). Although the top 6 predominant phyla are identical to those of neighbouring microbiota, the composition, i.e. relative contribution of each phyla to those microbiota seems fairly different. In the oral cavity, microbiota are dominated by Firmicutes followed by Proteobacteria and Bacteroidetes (overall 50% Gram-positive bacteria), whereas the microbiome of the nostril contains more than 80% gram-positive bacteria, mostly Actinobacteria and Firmicutes [15]. These data, therefore, suggest different dynamics (i.e., different biological equilibria) in the NP microbiome.
There was a high inter-individual variability in the composition of the microbiota up to phyla level, and in the relative abundance of the individual bacterial inhabitants (Table 1). This resulted in a limited core microbiome (as representing .0.1% of sequenes and being present in all 96 children) of specific phyla only, namely Proteobacteria and Firmicutes, however no single OTU was found universally. Because of the observed high inter-individual variation, we applied a less strict definition of core microbiome, i.e. OTUs present in more than 50% of all samples and representing .0.1% of the sequences. With this definition we observed a core microbiome of Moraxella, Haemophilus influenzae, Enhydrobacter (Proteobacteria), Streptococcus, Dolosigranulum (Firmicutes), and Corynebacterium (Actinobacteria) ( Table 1).
Principal component analysis identified 3 distinct clusters of microbiota profiles correlating strongly with a predominance (.50% of sequences per sample) of single OTU's, i.e. Moraxella (OTU 1), Haemophilus influenzae (OTU 2), and Streptococci (OTU 3), respectively, connected by a group of community profiles representing mixed microbiota ( Figure 2). Additionally, we observed transition zones for microbiota profiles between the Haemophilus-and Moraxella-dominated clusters, but not between Haemophilus-or Moraxella-dominated clusters and the Streptococci-dominated cluster, which might implicate potential interactions between microbiota profiles.
Since respiratory and invasive infections are associated with fall/winter season, we analysed our samples concordantly. When distinguished by the time of sampling (fall/winter versus spring) groups of children did not differ significantly in demographics or life style characteristics, infectious symptoms, medical history or environmental parameters (Table S1), reducing the likelihood of internal confounders as a cause of potential seasonal correlations with microbiota profiles. However, with respect to microbiota profiles, we observed marked differences between samples obtained in fall/winter versus samples from spring ( Table 2). In samples obtained in late fall and winter, we observed a predominance of Proteobacteria (relative abundance (% of all sequences): 75% versus 51% in spring), Fusobacteria (absolute abundance (% of children): 14% versus 2% in spring), and Cyanobacteria (absolute abundance: 64% versus 30% in spring; relative abundance: 0.27% versus 0.09% in spring) were significantly more abundant compared to spring, whereas Bacteroidetes were more frequently present in samples obtained in spring (relative abundance: 19% versus 3% in fall/winter, absolute abundance: 91% versus 54% in fall/winter) (Figure 3a). On OTU level we observed amongst others more Bacillus, Brevibacillus and Lactobacillus species, and Flavobacterium and B. fragilis (both Bacteroidetes) in samples from spring compared to fall/winter. In addition we found less a-Proteobacteria, Oxalobacteriaceae, Microbacteriaceae, Ralstonia, Pseudomonas and Acidovorax (all Proteobacteria), Cyanobacteria, and Porphyromonas catoniae (Bacteroidetes) in samples from spring compared to fall/ winter ( Figure 3b). When re-evaluating the core microbiome per individual season (i.e. OTUs present in more than 50% of samples of a certain season), we observed an additional core of Proprionibacterium and Cyanobacteria for fall/winter and an additional core of Flavobacteria, Brevibacillus and Bacillus (almost exclusively) for spring. The latter groups of bacteria, i.e. Bacteroidetes and (Brevi)bacillus and Lactobacillus species, are generally related to protection against overgrowth of pathogenic species due to the production of bacteriocins and other inhibitory substances [7,16]. In other microbiota like the gastrointestinal and vaginal tract they are highly related to maintenance of a balanced microbiome as well [17][18][19]. Since infections with respiratory pathogens, especially pneumonia, are strongly related to fall and winter season [9,10], the presence and abundance of these bacteria in respiratory microbiota in spring might therefore suggest in general a more balanced respiratory microbiome in this specific season as well protecting against onset of respiratory or invasive infections.
We tested all samples for the presence of respiratory viruses by q-PCR methods and detected one or more viruses in 67% of samples (Table S3). We found no evidence for associations between the observed seasonal shift in microbiota and the overall presence of respiratory viruses, nor for any of the individual viruses when tested by SAM analysis (Figure 3b). Although these data do not exclude an effect of respiratory viruses on microbiota composition, they do suggest other season-related factors like environmental factors (temperature, humidity, smoke exposure, crowding), or nutrient-or vitamine-related effects, or a combina-tion of factors might be important for the observed shifts in microbiota profiles [3,16,20]. Interestingly, day-care attendance or smoke exposure could not be related to the observed shifts in microbiota, although the latter was encountered very rarely. This further underlines that different or more complex effect may be responsible for the observed phenomenon. In addition, no association was observed between seasonal changes in microbiota and recent antibiotic use, which was rather limited in this population (Table 1, Figure 3b). Because of the explicit correlation between season and microbiota profiles, correlations between other environmental determinants and microbiota could not be accurately tested, which underlines that in future studies one needs to control and power for seasonal effects.
As previously mentioned, the majority of the children (87%) had a predominant Gram-negative NP community profile (.50% Gram-negative bacteria). This is probably due to predominance of Gram-negative Moraxella and Haemophilus species known to reside specifically at this body site [3,4]. On average, 76% of the overall NP microbiome in children was composed of Gramnegative bacteria, however with a wide range of 9-99% Gramnegative bacteria per sample. In addition, there was a higher contribution of Gram-negative bacteria to microbiota obtained in fall/winter (81%) compared to spring (72%) (Independent samples t-test: p = 0.044). This could potentially explain some of the observed differences between gram-negative ratios at the NP microbiome and at other human microbiota, where seasonal changes in composition have so far not been studied.
Finally, we studied the inter-individual diversity in NP microbiota overall and in relation to season. We observed highly diverse microbiota, with on average 40 OTUs per sample, and a high inter-individual diversity with 20-87 OTUs per individual. With respect to season, there was no significant difference in diversity between fall/winter (average: 38 OTUs, range: 20-77) and spring (average: 43 OTUs, range: 20-87) (independent samples T-test: p = 0.083).
As internal control, we compared conventional culturing results for the potential pathogens S. pneumoniae (71% positive), H. influenzae (69% positive) and M. catarrhalis (88% positive) with sequencing results for the OTUs Streptococcus (OTU 3, OTU 23, OTU 31), H. influenzae (OTU 2, OTU 10), and Moraxella (OTU 1, OTU 6) respectively ( Figure S1) and found strong correlations between S. pneumoniae and Streptococcus OTU 3 and 31 (p,0.0001) but not OTU 23, which is probably another Streptococcus species, and H. influenzae and H. influenzae OTU 2 and OTU 10 (p,0.0001). For M. catarrhalis, we were only able to find a positive correlation between M. catarrhalis and Moraxella OTU 6 with independent samples t-test (p = 0.003) but not with Spearman's Correlation, which may be explained by the low number of Moraxella negative individuals, making a comparison between binary and quantitative data difficult. Also the presence of other Moraxella species with high sequence homology might interfere with a strong correlation between these results.
In conclusion, to our knowledge this is the first report describing in detail the composition of and the variability within the human NP microbiota assessed at the depth of next generation sequencing. In line with other human body habitats, we found a complex, diverse and highly variable microbiota with a relatively limited core microbiome. There is considerable seasonal variation in NP microbiota. This implies the time of sampling should be considered when describing or comparing NP microbiomes, and preferably controlled for when other potential determinants like the impact of viruses or antibiotics on microbiota profiles or the correlation between microbiota profiles and diseases will be tested. Whether these seasonal changes in composition of the NP microbome are causally related to seasonal occurrence of respiratory tract infections remains to be determined, though seems relevant for further understanding of pathogenesis of infectious diseases and in the long run potentially for understanding of effects of current and design of future preventive measures.

Samples
We randomly selected 150 NP samples from a cohort of 330 healthy children 18 months of age who had participated in a randomised controlled trial studying the effect of reduced-dose schedules of 7-valent pneumococcal conjugate vaccine (PCV-7) performed in a general community in the Western part of The Netherlands where the control children received PCV-7 only after the trial was finished at the age of 24 months [21]. An  Table 1. Thirty most common OTUs or 'species-level' phylotypes (ranked by predominance, i.e. absolute presence among the approx. 1 100 000 reads).

Microbial cultivation
The nasopharyngeal swabs were plated onto a 5% sheep blood agar plate, a 5% sheep blood agar plate with 5 mg/L gentamicin, a chocolate agar plate and a Haemophilus chocolate agar plate. Agar plates were incubated at 35uC for 48 h; the blood agar plate aerobically, the blood agar plate with gentamicin and the chocolate agar plates with raised CO2. Identification of S. pneumoniae, H. influenzae, M. catarrhalis and S. aureus was based on colony morphology and conventional methods of determination.

Real time PCR for bacterial DNA
The total bacterial load of the samples was established by quantitative PCR.  Real-time PCR for viruses One 200 ml aliquot of swab ''rinse'' solution was used to extract viral nucleic acids using the MagNA Pure LC total nucleic acid isolation kit (Roche Diagnostics, catalogue 03 038 505 001, Basel, Switzerland) as described previously [22]. Detection of viral pathogens was performed in parallel, using real-time PCR assays for bocavirus (HBoV), polyomaviruses (WUPyV and KIPyV), respiratory syncytial virus (RSV) A and B, influenzavirus (IV) A and B, para-influenzavirus (PIV) 1-4, human rhinoviruses (HRV), adenoviruses, human coronavirus OC43, NL63, HKU and 229E, and human metapneumovirus (hMPV). Real-time PCR procedures were performed as described previously [22] Briefly, samples were assayed in duplicate in a 25 ml reaction mixture containing 10 ml (c)DNA, 12.

Data Analysis
GS-FLX sequencing data were processed as previously described [23]. In brief, we trimmed data by removing primer sequences and low-quality data, sequences that did not have an exact match to the reverse primer, that had an ambiguous base call (N) in the sequence, or that were shorter than 50 nt after trimming. We then used the GAST algorithm [24] to calculate the percent difference between each unique sequence and its closest match in a database of 69816 unique Eubacterial and 2779 unique Archaeal V5-V6 sequences, representing 323499 SSU rRNA sequences from the SILVA database [25]. Taxa were assigned to each full-length reference sequence using several sources including Entrez Genome entries, cultured strain identities, SILVA, and the Ribosomal Database Project Classifier [26]. In cases where reads were equidistant to multiple V5-V6 reference sequences, and/or where identical V5-V6 sequences were derived from longer sequences mapping to different taxa, reads were assigned to the lowest common taxon of at least two-thirds of the sequences. The operational taxonomic units (OTUs) were created by aligning unique sequences and calculating distance matrices as previously described [23] and using DOTUR [27] to create clusters at the 3% level.
Only those sequences that were found at least 5 times were included in the analyses. This strict and conservative approach was chosen to preclude inclusion of sequences from potential contamination or sequencing artefacts. To compare the relative abundance of OTUs among samples, the data were normalized for number of sequenced reads obtained for each sample. To reduce the influence of abundant taxa on principal component analyses, the normalized abundance data were log2 transformed. Unsupervised data analysis, Principle Component Analysis, and hierarchical clustering was performed using MeV software package as part of TM4 microarray software suite [28].
Seasonal differences in phyla distribution were studied by using Mann Whitney U test (SPSS software Version 15.0). Seasonal differences in OTU patterns as well as potential correlations between respiratory viruses and OTU patterns were studied by the Significant Analysis of Microarrays (SAM analysis) -a nonparametric statistical technique for finding significant differences between microarray data of groups based on experimental conditions [29], implemented in the MeV software package [28]. To determine significant differences between microbiota profiles, we used Pearson's correlation with average linkage clustering method and a FDR significance criterion of ,0.05.
Independent samples T-test, and Spearman correlation coefficients (SPSS software Version 15.0) were used for testing correlations between conventional cultures and pyrosequencing data. Independent samples t-test was also used to compare the contribution of gram negative and positive bacteria to the microbiota in different seasons and to test for differences in diversity between seasons. Figure S1 Graphs showing correlation between conventional culture results for S. pneumoniae, H. influenzae, and M. catarrhalis on the x-axis (absent/present) and sequencing results (as % of total microbiota profile) for Streptococcus (OTU 3), H. influenzae (OTU 2 plus 10), and Moraxella (OTU 6) on the y-axis. (TIFF)