Microbial communities of the upper respiratory tract and otitis media in children.

Streptococcus pneumoniae asymptomatically colonizes the upper respiratory tract of children and is a frequent cause of otitis media. Patterns of microbial colonization likely influence S. pneumoniae colonization and otitis media susceptibility. This study compared microbial communities in children with and without otitis media. Nasal swabs and clinical and demographic data were collected in a cross-sectional study of Philadelphia, PA, children (6 to 78 months) (n=108) during the 2008-2009 winter respiratory virus season. Swabs were cultured for S. pneumoniae. DNA was extracted from the swabs; 16S rRNA gene hypervariable regions (V1 and V2) were PCR amplified and sequenced by Roche/454 Life Sciences pyrosequencing. Microbial communities were described using the Shannon diversity and evenness indices. Principal component analysis (PCA) was used to group microbial community taxa into four factors representing correlated taxa. Of 108 children, 47 (44%) were colonized by S. pneumoniae, and 25 (23%) were diagnosed with otitis media. Microbial communities with S. pneumoniae were significantly less diverse and less even. Two PCA factors were associated with a decreased risk of pneumococcal colonization and otitis media, as follows: one factor included potentially protective flora (Corynebacterium and Dolosigranulum), and the other factor included Propionibacterium, Lactococcus, and Staphylococcus. The remaining two PCA factors were associated with an increased risk of otitis media. One factor included Haemophilus, and the final factor included Actinomyces, Rothia, Neisseria, and Veillonella. Generally, these taxa are not considered otitis media pathogens but may be important in the causal pathway. Increased understanding of upper respiratory tract microbial communities will contribute to the development of otitis media treatment and prevention strategies.

INTRODUCTION

Otitis media (OM) is associated with over 10 million office visits in the United States each year and is the leading diagnosis for the prescription of antibiotics for children. OM arises in the complex microbial community of the upper respiratory tract. The bacteria most often associated with OM are Streptococcus pneumoniae, Haemophilus influenzae, and Moraxella catarrhalis (1). S. pneumoniae asymptomatically colonizes up to 50% of children (2–7). Colonization of the upper respiratory tract is the first step in infection; even transient colonization provides an opportunity for S. pneumoniae to invade the middle ear space.

Shifts in the composition of microbial communities of the upper respiratory tract are associated with age, vaccination, antibiotic use, and upper respiratory tract infection (7–12). Bacterial interactions can impact microbial community composition as well as disease incidence (13–15). Studies of children’s upper respiratory tract microbial flora often rely on culture to identify taxa of interest (16–19). In our bodies, bacterial cells outnumber human cells (20, 21), and a majority of the bacteria that colonize and infect humans are not routinely cultured or cannot be cultured (22). Therefore, it is likely that the members of the microbial flora that influence S. pneumoniae colonization and subsequent OM are not limited to taxa routinely identified by culture.

The use of 16S rRNA full- and partial-gene sequencing provides an alternative to culture-based methods (23–29). The 16S rRNA gene sequence contains both conserved and hypervariable regions, which can be used for taxonomic classification (30). Hypervariable regions V1 and V2 have been shown to provide a high level of resolution in species discrimination, including identification of the most common OM pathogens (31).

Currently, our knowledge of the nasal microbial flora and its association with OM is limited. In this cross-sectional study of children aged 6 months to 6 years, we describe the complex microbial communities present during upper respiratory tract infection, when children are most susceptible to developing OM (8). We also examine associations among members of microbial communities, S. pneumoniae colonization, and OM.

RESULTS

Study population.

Of 108 children, 20 (18.5%) were 6 to 12 months of age, 29 (26.9%) were 12 to 24 months, 19 (17.6%) were 24 to 36 months, and 40 (37.0%) were over 36 months. A total of 60 (55.6%) children were female, 75 (69.4%) were African-American, and 74 (68.5%) lived in households with other children. Forty-seven (43.5%) of the children were culture positive for S. pneumoniae. Nearly one-fourth (23.2%) of the children were diagnosed with OM. S. pneumoniae was isolated by culture in 15 (60%) of the OM cases. Most children (n = 95; 88.0%) had received the appropriate number of heptavalent pneumococcal conjugate vaccinations for their age. Logistic regression was used to assess the association between OM and demographic risk factors, including age, race, and gender. No significant associations were found, perhaps because all study subjects presented with upper respiratory tract infection, a major OM risk factor (8).

Sequence diversity.

A mean of 1,087.5 sequences per sample (standard deviation [SD], 496.01) was included in the analyses (114,397 total sequences). The average sequence length was 234 bp. The mean Shannon diversity and evenness indices for all samples were 3.0 (SD, 0.83) and 0.62 (SD, 0.12), respectively. Both indices were significantly lower in S. pneumoniae culture-positive samples than in culture-negative samples (2.66 versus 3.25, P = 0.0002, and 0.58 versus 0.66, P = 0.002, respectively). Of 108 children, 7 (6.5%) had used antibiotics within the past 7 days, 15 (13.9%) within the past 14 days, 17 (15.7%) within the past 21 days, and 32 (29.6%) within the past 6 months. Separate analyses of variance (ANOVAs) examined the relationship between the Shannon diversity and evenness indices, pneumococcal colonization, and antibiotic use and the interaction between them. Pneumococcal colonization had a significant effect on diversity (F test, P = 0.0002 for Shannon diversity, P = 0.002 for evenness), but antibiotic use (within the past 7 days, 14 days, 21 days, or 6 months) did not. There was no significant interaction between pneumococcal colonization and antibiotic use. This suggests that antibiotic use within the past 6 months does not confound the relationship between diversity indices and pneumococcal colonization. The Shannon diversity and evenness indices were not significantly different by OM diagnosis (t test, P = 0.6 and 0.64, respectively).

Using the Ribosomal Database Project (RDP) Classifier tool, pyrosequencing provided sufficient information to classify 68.5% of the sequences at the genus and subgenus levels with 90% confidence, 25.7% at the family and subfamily levels, 1.8% at the order and suborder levels, 0.6% at the class and subclass levels, and 3.4% at the phylum level. We identified 389 operational taxonomic units (OTUs). OTU proportions of ≥1% in the overall population (n = 16) are listed in Table 1. OTU proportions were compared by S. pneumoniae colonization (Fig. 1A) and by OM diagnosis (Fig. 1B). The most frequent OTU proportions included potential OM pathogens Streptococcus, Moraxella, and Haemophilus. In the S. pneumoniae-colonized group, we observed fewer “rare” bacteria than in the S. pneumoniae-negative group, as evidenced by the difference in the “all other bacteria” bar, suggesting lower levels of overall diversity in the S. pneumoniae-colonized group (Fig. 1A). The OTU proportion of potential OM pathogens Streptococcus (23.1%) and Haemophilus (10.2%) was greater in the OM group than in the OM-absent group (16.3% and 3.0%, respectively).

TABLE 1

Most frequent nasal swab OTUs

OTU	Frequency (%)b
Unclassified Moraxellaceae	19.00
Streptococcus	17.86
Corynebacterium	7.04
Moraxella	6.46
Haemophilus	4.66
Unclassified Pasteurellaceae	4.09
Staphylococcus	3.84
Acinetobacter	3.44
Dolosigranulum	3.21
Propionibacterium	3.13
Unclassified Proteobacteria	2.59
Lactococcus	2.58
Neisseria	1.45
Actinomyces	1.24
Rothia	1.13
Veillonella	1.05

Frequency of ≥1%.

Percentage of total sequences per nasal microbial community, i.e., per child.

FIG 1

Comparison of the most frequent operational taxonomic unit (OTU) distributions grouped by S. pneumoniae colonization (determined by culture) (A) and by otitis media (OM) diagnosis (B). Percentages of total sequences per nasal microbial community, i.e., per child, that each OTU represented are shown.

PCA.

OTUs representing ≥1% of the microbial flora were used for principal component analysis (PCA), which grouped microbial community members into independent factors. The four resulting PCA factors (factors A to D) represent linear relationships among correlated taxa. Factor A includes Corynebacterium, Dolosigranulum, and Streptococcus (negatively correlated); factor B includes Lactococcus, Propionibacterium, unclassified Proteobacteria, and Staphylococcus; factor C includes unclassified Pasteurellaceae, Moraxella, and Haemophilus; and factor D includes Veillonella, Neisseria, Rothia, and Actinomyces.

Logistic regression analysis was used to examine associations between PCA factors and S. pneumoniae colonization or OM diagnosis. Two of the four factors were significantly associated with a decreased risk of S. pneumoniae colonization (odds ratio [OR] [95% confidence interval {CI}]), as follows: factor A (0.55 [0.35, 0.86]) and factor B (0.34 [0.20, 0.59]) (Table 2). These same two factors were also associated with a decreased risk of OM (OR [95% CI]), as follows: factor A (0.54 [0.30, 0.96]) and factor B (0.44 [0.21, 0.91]) (Table 3). Factors A and B include bacteria that potentially are interfering members of the normal flora (Tables 2 and 3). Furthermore, the presence of Corynebacterium, Dolosigranulum, and Lactococcus was negatively correlated with the genus Streptococcus, with r values of −0.25 (P = 0.01), −0.26 (P = 0.006), and −0.19 (P = 0.05), respectively. Two other factors were significantly associated with an augmented risk for OM (OR [95% CI]), as follows: factor C (1.73 [1.10, 2.73]) and factor D (2.24 [1.26, 3.97]) (Table 3). Factor C includes Haemophilus and unclassified Pasteurellaceae, and factor D includes Veillonella, Rothia, Actinomyces, and Neisseria (Tables 2 and 3).

TABLE 2

Associations between S. pneumoniae colonization determined by culture and factors obtained by PCA as well as individual taxa within factors

Factor	OR (95% CI)b
	PCA factor	Component
Factor A	0.55 (0.35, 0.86)
Corynebacterium		0.99 (0.95, 1.03)
Dolosigranulum		0.99 (0.91, 1.06)
Streptococcusa		1.06 (1.03, 1.09)
Factor B	0.34 (0.20, 0.56)
Lactococcus		0.73 (0.61, 0.89)
Propionibacterium		0.89 (0.79, 1.00)
Unclassified Proteobacteria		0.94 (0.88, 1.02)
Staphylococcus		0.98 (0.93, 1.02)
Factor C	0.94 (0.64, 1.39)
Unclassified Pasteurellaceae		1.08 (1.00, 1.16)
Moraxellaa		1.03 (0.99, 1.07)
Haemophilus		1.02 (0.98, 1.05)
Factor D	0.61 (0.37, 1.02)
Veillonella		0.98 (0.83, 1.15)
Neisseria		0.85 (0.70, 1.04)
Rothia		0.82 (0.62, 1.08)
Actinomyces		0.77 (0.58, 1.02)

 The OTU was negatively correlated with the factor; factor loadings were otherwise positive.

OR and 95% CI values were obtained by separate logistic regression analyses of each factor and each taxa within each factor.

TABLE 3

Associations between OM diagnosis and factors obtained by PCA as well as individual taxa within factors

Factor	OR (95% CI)b
	PCA factor	Component
Factor A	0.54 (0.30, 0.96)
Corynebacterium		0.90 (0.81, 0.99)
Dolosigranulum		0.91 (0.80, 1.04)
Streptococcusa		1.02 (1.00, 1.05)
Factor B	0.44 (0.21, 0.91)
Staphylococcus		0.86 (0.71, 1.06)
Lactococcus		0.92 (0.78, 1.09)
Unclassified Proteobacteria		0.93 (0.83, 1.05)
Propionibacterium		0.97 (0.87, 1.07)
Factor C	1.73 (1.10, 2.73)
Unclassified Pasteurellaceae		1.08 (1.01, 1.15)
Haemophilus		1.04 (1.00, 1.08)
Moraxellaa		0.99 (0.96, 1.06)
Factor D	2.24 (1.26, 3.97)
Veillonella		1.37 (1.07, 1.76)
Rothia		1.35 (1.03, 1.77)
Actinomyces		1.34 (1.07, 1.68)
Neisseria		1.13 (0.99, 1.29)

 The OTU was negatively correlated with the factor; factor loadings were otherwise positive.

OR and 95% CI values were obtained by separate logistic regression analyses of each factor and each taxa within each factor.

As factor D taxa are not generally recognized as OM pathogens, we explored additional associations and found that the OTU proportion of each factor D taxon (Veillonella, Rothia, Actinomyces, and Neisseria) in a sample was significantly higher in swabs obtained from children with OM (t test, P = 0.001, 0.02, 0.002, and 0.03, respectively). In addition, samples obtained from children had significantly higher factor D scores if antibiotics had been prescribed within the past 21 days or the past 6 months (t test, P = 0.02 or 0.05, respectively).

DISCUSSION

We describe nasal microbial communities in children with upper respiratory tract infection with and without concurrent OM. We observed lower microbial diversity when S. pneumoniae is present. Using principal component analysis to group microbial community members into factors revealed taxa associated with a decreased risk for S. pneumoniae colonization and OM. These taxa may be useful for the development of multispecies probiotics for OM. We also identified taxa associated with an increased risk for OM. Further study is warranted to define their role in OM pathogenesis.

Vaccine-driven immunological pressures and increasing antibiotic resistance may result in shifts in the upper respiratory tract microbial community that affect the distribution and pathogenic potential of OM-associated bacteria. We used complete-linkage clustering in our diversity calculations through the RDP pyrosequencing pipeline. Average-linkage clustering may yield more accurate estimates of diversity with regard to identifying rare microbiota in short amplicon (60-bp) data sets (32, 33). However, Quince et al. found that diversity inflation was mitigated in data sets of longer sequences (250 bp), and our sequences are of similar lengths (average, 234 bp) (33).

In culture-positive S. pneumoniae samples, diversity and evenness indices were lower than those in culture-negative samples. Children with OM have decreased levels of healthy normal flora bacteria like viridans group streptococci and diphtheroids, while the OTU proportion of OM pathogens increases by 2- to 3-fold (34). The association between diversity indices and OM was not significant here. However, all children in the study were seen for a sick visit, which in this population may have attenuated the differences in microbial diversity by OM status.

The two factors associated with a reduced risk of OM were dominated by the genera Corynebacterium, Dolosigranulum, Propionibacterium, Lactococcus, and Staphylococcus. Together, these two factors (factors A and B) (Tables 2 and 3) identify potential protective members of the normal flora that may interfere with OM pathogens. Streptococcus loaded negatively on factor A and is negatively correlated with both Corynebacterium and Dolosigranulum, suggesting a competitive interaction between this potential OM pathogen and these two genera. Lemon et al. also observed negative correlations between members of the normal flora such as the Corynebacteriaceae and Propionibacteriaceae families and the Staphylococcaceae family, which includes potential OM pathogens, in nasal samples obtained from healthy adults by PhyloChip analysis (35).

Samples obtained from children with high factor C loading are characterized by OM diagnosis as well as by the presence of Haemophilus and the absence of Moraxella. This is consistent with previous literature linking Haemophilus with OM (1). In addition, H. influenzae has been negatively correlated with M. catarrhalis during upper respiratory tract infection (36). It is important to note that as sequence data could not classify at the species level, factor C likely contains other species of Moraxella in addition to the OM pathogen M. catarrhalis.

Several mechanisms may account for the association between factor D and an increased risk of OM. Although the bacteria associated with factor D are not normally recognized as OM pathogens, they may be important to the causal pathway. Higher levels of colonization by anaerobes such as Actinomyces and Veillonella have been identified during acute OM in children up to 2 years of age (37). Neisseria has also been isolated from children with chronic OM (38). In this study, the OTU proportion of Neisseria was significantly higher in swabs obtained from children with OM. Previous studies show that bacterial species can impact the colonization and pathogenicity of other species (39, 40). Therefore, individual species or even specific strains, which cannot cause disease themselves, can interact in a synergistic fashion when in the presence of other coinfecting species. Sibley et al. also found a synergistic effect on pathogenicity between chronic pulmonary infection pathogen Pseudomonas aeruginosa and both Actinomyces and Rothia spp. in the Drosophila model of polymicrobial infections (41). Factor D taxa (Veillonella, Neisseria, Rothia, and Actinomyces) may or may not be pathogenic themselves but may be able to increase the pathogenicity of OM pathogens. Further study is needed to better elucidate the associations among these taxa, known OM pathogens, and the development of OM.

Antibiotic resistance may also play a role in the observed association between factor D taxa and OM. Samples in this study obtained from children who were prescribed antibiotics had significantly higher factor D scores. Both Neisseria and Veillonella have been shown to have widespread antibiotic resistance (42–44). Fluoroquinolone resistance has emerged recently in Neisseria (44), and Veillonella spp. are among the most prevalent tetracycline-resistant bacteria in the oral cavity (42). The survival of Streptococcus mutans increases during antibiotic treatment when cocultured with Veillonella parvula in a biofilm (45). As a result of repeated OM-related courses of antibiotics, children with a history of recurrent OM certainly have the potential to develop a more resistant microbial flora.

One limitation of our study is that data on viral pathogens, although important in OM, were not collected. In addition, as this was a cross-sectional study, we could not characterize within-subject variability. All children presented with upper respiratory tract infection symptoms. Additional studies comparing healthy and sick children are necessary. The Roche/454 Life Sciences pyrosequencing platform is one of several methods available for high-throughput sequencing and circumvents the need for creating clone libraries. However, a limitation of the short-read sequencing technology is that it does not provide sufficient sequence data for the classification of species. Thus, the majority of sequences could be classified only at the genus level. In addition, due to potential PCR and sequencing biases, the relationship between the number of sequences per OTU and the bacteria in the population may not be linear.

The present study is based on nasal sampling, which has the advantage of being easier and is more comfortable for children than nasopharyngeal sampling (46). Other studies have used data obtained from nasopharyngeal sampling to characterize pneumococcal colonization patterns in children (6, 47). However, Rapola et al. have shown that the rates of OM pathogens S. pneumoniae and H. influenzae isolation do not differ by site (48). Although regions of the upper respiratory tract harbor distinct microbial communities, they are physically connected, and thus, a high level of overlap is observed between sites (35).

Our culture-independent classification of bacteria present in children’s nasal mucosa during upper respiratory tract infection contributes to the understanding of nasal microbial communities, development of OM, and polymicrobial interactions. Increased understanding of complex microbial communities will advance the development of efficacious prevention and treatment protocols for OM, including probiotic therapies that target high-risk individuals, which in turn can reduce the risk of permanent damage caused by OM.

MATERIALS AND METHODS

Study design and participants.

Nasal swabs were collected as part of a cross-sectional study of Philadelphia, PA, children presenting with upper respiratory tract infection symptoms at one of six primary care clinics within the Pediatric Research Consortium of the Children’s Hospital of Philadelphia during the 2008-2009 winter respiratory virus season. Samples obtained for this study (n = 108) were randomly selected from swabs collected between 9 December 2008 and 2 January 2009 from participants aged 6 months to 6 years (n = 237). Demographic data collected from each participant through a case report form at the clinic visit included age, gender, race, ethnicity, and number of children in household. Clinical data, including diagnosis at visit, comorbidities, vaccination history, and recent antibiotic use, were obtained from electronic medical record extraction, including International Classification of Diseases, 9th Edition (ICD-9) codes. Trained research assistants obtained informed consent for participation in the study during the clinic visit. The Institutional Review Board (IRB) of the University of Pennsylvania approved the study protocol.

Nasal swab collection and processing.

Anterior nasal swabs (49) were collected from each child by inserting a rayon-tipped swab moistened with saline approximately 1 cm into the anterior nares. Swabs were placed in STGG (skim milk, Oxoid tryptone soya broth, glucose, and glycerol) transport media until processed (19). Pneumococcal colonization was determined by standard microbiological methods. DNA was extracted using the QIAamp DNA minikit (Qiagen, Valencia, CA) protocol for extraction of Gram-positive bacteria from nasal swabs, with one modification: DNA was eluted into 50 µl.

Roche/454 Life Sciences pyrosequencing of a 16S rRNA gene fragment targeting hypervariable regions V1 and V2 was used to describe the taxa of bacteria colonizing the nasal mucosa. 16S rRNA gene amplification was performed as described (23). The forward primer sequence included the Roche/454 Life Sciences adaptor primer B, TC-linker, and the 27F conserved bacterial 16S rRNA forward primer: 5′ GCCTTGCCAGCCCGCTCAGTCAGAGTTTGATCCTGGCTCAG 3′ (23). The reverse primers carried the Roche/454 Life Sciences adaptor primer A, a unique eight-base barcode (denoted here as 8 N), CA-linker, and the 338R conserved bacterial primer: 5′ GCCTCCCTCGCGCCATCAGNNNNNNNNCACTGCTGCCTCCCGTAGGAGT 3′ (23). PCR reactions were done in triplicate, and products were cleaned by QIAquick PCR purification (Qiagen). Samples were pooled in equimolar amounts and submitted for pyrosequencing.

Pyrosequencing results analysis.

Amplicons were sequenced at the Yale University Center for Genomics and Proteomics using the SLR70 Roche/454 Life Sciences pyrosequencing platform. The sequences were first searched for the linker, primers, and their reverse complements using in-house programs (http://graphics.med.yale.edu/trim/). Two errors were allowed, and the identified primer sequences were trimmed from each sequence read. Sequence reads that did not contain the 5′-end primer were removed from the data set. The same program was also used for barcode identification. Barcodes were identified within the first 15 bases of the reads, with one error allowed. Sequence reads were binned into separate FASTA files based on the individual barcodes, and the barcode sequences were trimmed. Current methods for chimera detection are designed for longer sequences (~500 bp or longer), such as ChimeraSlayer and Bellerophon (50, 51; http://microbiomeutil.sourceforge.net/). In an effort to minimize the number of chimeras, we set strict alignment criteria; all sequences had to align to at least 100 bp of the 16S rRNA gene, and any sequence aligning outside the 27 and/or 338 positions of the 16S rRNA gene was discarded.

Sequences were analyzed using tools available from the Ribosomal Database Project (RDP). Individual sequences were aligned using the Aligner tool. Aligned sequence files for each child were processed by complete-linkage clustering using distance criteria. We used a maximum cluster distance cutoff of 3% (97% identity). These data were used to calculate the Shannon diversity and evenness indices (25, 26, 52, 53). Shannon’s index measures the degree of nonredundancy within a microbial community, and the evenness index measures the distribution of species or groups within a sample.

Sequences were classified at 90% by the RDP Classifier tool (52). The RDP website recommends using a 50% bootstrap value for short pyrosequences to increase the number of sequences classified (54). Claesson et al. (54) and Wang et al. (52) note that the proportion of sequences classified correctly increases as the bootstrap cutoff is increased. Therefore, we selected a 90% cutoff in order to classify the sequences with greater accuracy. The sequence data could not be classified at the species level; OTUs (operational taxonomic units) were defined by grouping together all sequences belonging to the same genus. Sequences that were unclassified at the genus level were classified and grouped at the next lowest taxonomic level. In total, there were 389 OTUs. We normalized the data by calculating the OTU proportion—the proportion of total sequences per nasal microbial community, i.e., per child, that each OTU represented. These data were used for principal component analysis (PCA) and all subsequent analyses.

Data analysis.

Associations among community diversity indices, the presence of S. pneumoniae by culture, OM diagnosis, and antibiotic use within the past 7 days, 14 days, 21 days, or 6 months were explored using Student’s t test, Fisher’s exact chi-square test, and analysis of variance (ANOVA). PCA was used to group microbial community members into four factors representing linear relationships among selected component taxa (SAS 9.1.3; SAS Institute, Cary, NC). OTUs that were more frequent than 1% of the microbial community were included as component taxa (n = 16) in the PCA. An eigenvalue of 1 and an orthogonal rotation were specified for PCA. PCA factor components had significant loadings of at least ±0.4. Associations among PCA factors, PCA factor scores, and the presence of S. pneumoniae by culture and OM diagnosis were explored using Student’s t test, correlation, and logistic regression.