The temporal dynamics of the tracheal microbiome in tracheostomised patients with and without lower respiratory infections.

BACKGROUND: Airway microbiota dynamics during lower respiratory infection (LRI) are still poorly understood due, in part, to insufficient longitudinal studies and lack of uncontaminated lower airways samples. Furthermore, the similarity between upper and lower airway microbiomes is still under debate. Here we compare the diversity and temporal dynamics of microbiotas directly sampled from the trachea via tracheostomy in patients with (YLRI) and without (NLRI) lower respiratory infections.
METHODS: We prospectively collected 127 tracheal aspirates across four consecutive meteorological seasons (quarters) from 40 patients, of whom 20 developed LRIs and 20 remained healthy. All aspirates were collected when patients had no LRI. We generated 16S rRNA-based microbial profiles (~250 bp) in a MiSeq platform and analyzed them using Mothur and the SILVAv123 database. Differences in microbial diversity and taxon normalized (via negative binomial distribution) abundances were assessed using linear mixed effects models and multivariate analysis of variance. RESULTS AND DISCUSSION: Alpha-diversity (ACE, Fisher and phylogenetic diversity) and beta-diversity (Bray-Curtis, Jaccard and Unifrac distances) indices varied significantly (P<0.05) between NLRI and YLRI microbiotas from tracheostomised patients. Additionally, Haemophilus was significantly (P = 0.009) more abundant in YLRI patients than in NLRI patients, while Acinetobacter, Corynebacterium and Pseudomonas (P<0.05) showed the inverse relationship. We did not detect significant differences in diversity and bacterial abundance among seasons. This result disagrees with previous evidence suggesting seasonal variation in airway microbiotas. Further study is needed to address the interaction between microbes and LRI during times of health and disease.

Introduction

The population of children with a tracheostomy is increasing [1] and these children are at high risk of developing lower respiratory infections (LRI) [2] that require intensive care [3]. A retrospective analysis of 917 children aged 0–18 years from 36 children's hospitals who received a tracheostomy in 2002 demonstrated that in the 5-year follow-up period, the mean number of hospitalizations experienced per child was 3.8 (SD: 4.4; range: 0–34) and 46% of the hospitalizations were for “respiratory” diagnoses [2].

Despite this risk of LRI, there is no consensus about how providers should treat common acute respiratory infections in tracheostomised children [3, 4]. While common respiratory viruses and bacteria in the airway tract are identifiable utilizing standard laboratory techniques and have provided interesting avenues for research (e.g., [5]), these common organisms may only represent a small fraction of the commensal and pathogenic microbes living in the airways of an individual [6]. Indeed, the traditional conceptual model of these potentially life-threatening LRIs (i.e., pathogen causing disease in a sterile lung) may be too simplistic [7]. More current conceptual models of LRIs account for the highly functional bacterial communities (i.e., microbiota) that are present in the upper and lower airways (nose to lung) [8-13] and influence both immune [14-18] and inflammatory responses [18, 19].

However, the role that the airway microbiota plays during LRI is still poorly understood since: i) airway microbiome research focuses mainly on the nasopharynx and oral cavity and far less on the lower airways due to easier access and ethical considerations [20], and ii) most microbiome studies during both health and disease are cross-sectional [21-26] and neglect the temporal dynamics of the microbial communities [11]. Furthermore, upper airway samples are frequently used as a proxy for the lower airway, but the similarity of lower and upper airway microbiotas during disease is still under debate [27, 28]. Hence, since studying the dynamics of lower respiratory microbiotas during both health and disease (e.g., infection) is a challenging endeavor, new approaches or models [29] are needed to continue their investigation.

We have addressed these knowledge gaps by directly examining the dynamics of tracheal microbiotas (i.e., bypassing the nose and mouth) in a longitudinal study of 40 tracheostomised children and young adults with and without LRIs. We aimed to determine the differences in the microbiotas of these patients and their variation over one calendar year. We hypothesized that tracheal microbiotas differ in composition and structure according to whether or not individuals developed LRIs.

Materials and methods

Ethics

This study was approved by the Boston Children's Hospital Institutional Review Board (IRB), which requires that consent is obtained and documented prior to conducting study procedures and collection of samples for research. Written consent was obtained from all independent participants or their legal guardians using the Boston Children's Hospital IRB approved informed consent documents (IRB No P00007853).

Cohort

The Critical Care, Anesthesia, Perioperative Extension (CAPE) and Home Ventilation Program began in June 2007 for children and young adults with respiratory technology dependence in an effort to enhance outpatient and transitional-care services [30-32]. The participants’ diagnoses include, but are not limited to, spinal muscular atrophy, spastic quadriplegia with respiratory insufficiency, muscular dystrophy, mitochondrial disorders, spinal cord injuries, idiopathic hypotonia, and a range of complex chronic lung disease. All participants in this study resided in New England (USA) at the time of sampling.

Sample collection and storage

Samples were collected every meteorological season or quarter (Q) over the course of one calendar year starting in the fall of 2013 and ending in the summer of 2014: Q1 (fall) = 10/2013 to 11/2013, Q2 (winter) = 12/2013 to 02/2014, Q3 (spring) = 03/2014 to 05/2014 and Q4 (summer) = 06/2014 to 08/2014. All samples were collected when patients were healthy (i.e., no symptoms of LRI for previous four weeks), since our goal was to determine whether microbiotas from individuals who develop a LRI are different from those who do not develop a LRI. A LRI was defined as any illness causing increased mucus production and requiring increased oxygen delivery or higher ventilator settings over baseline. In order to ensure standardized sample collection, we observed parents or visiting nurses collecting the first sample in person during a home or clinic visit. The aspirate was stored at 0°C within 15 minutes of being collected. The samples were picked up by the study team, transported on ice to Boston Children’s Hospital, and stored at -80°C.

High-throughput sequencing

We aimed to sequence one tracheal aspirate sample per quarter (four quarters total) from both participants who acquired a LRI (yes LRI = YLRI) during the study and those who remained healthy (no LRI = NLRI), hence rendering four samples per participant. Q1 to Q4 samples in both YLRI and NLRI patients are isochronous. Total DNA was extracted using the QIAGEN QIAamp DNA Kit (Catalog # 51304). All samples were incubated in 200 μL of lysozyme-TE buffer pH = 8.0 for 30 minutes at 37°C. All extractions yielded >5 ng/μL of total DNA (as indicated by NanoDrop 2000 UV-Vis Spectrophotometer measuring). DNA extractions were prepared for sequencing using the Schloss’ MiSeq_WetLab_SOP protocol [33]. Each DNA sample was amplified for the V4 region (~250 bp) of the 16S rRNA gene and all libraries were sequenced together in a single run of the Illumina MiSeq sequencing platform at University of Michigan Medical School. Negative controls processed as above showed no PCR band on an agarose gel.

Sequence and statistical analyses

Raw FASTQ files were processed in mothur v1.35.1 [34] and indicated in the MiSeq SOP (www.mothur.org/wiki/MiSeq_SOP). Default settings were used to minimize sequencing errors [35]. Clean sequences were aligned to the SILVA123-based bacterial reference alignment at www.mothur.org. Chimeras were removed using uchime [36] and non-chimeric sequences were classified using a naïve Bayesian classifier [37]. Sequences were clustered into Operational Taxonomic Units (OTUs) at the 0.03 similarity threshold (species level). A consensus taxonomy was generated based on the classification of sequences clustered within an OTU. OTU sequence representatives and taxonomy were then converted to a BIOM file for subsequent analyses and all OTU singletons (n = 1) were eliminated. We normalized our samples using the negative binomial distribution as recommended by McMurdie and Holmes [38] and implemented in the Bioconductor package DESeq2 [39]. This approach simultaneously accounts for library size differences and biological variability. Microbial normalized counts generated this way are referred to as taxon abundances throughout the text. Trees for phylogenetic diversity calculations were constructed using FastTree and midpoint rooting [40].

Taxonomic alpha-diversity was estimated using Shannon, Fisher and ACE indices, while phylogenetic alpha-diversity was calculated by the Faith’s phylogenetic diversity index [41]. Beta-diversity was estimated using phylogenetic Unifrac (unweighted and weighted), Bray-Curtis and Jaccard distances. Dissimilarity between samples was explored using principal coordinates analysis (PCoA). Linear mixed-effects (LME) models analysis, as implemented in the lmer4 R package [42], was applied to both alpha-diversity indices and taxa (genera and phyla) abundances (response) while accounting for non-independence of subjects (random effect) and LRI (predictor). Time was modeled as “number of days since the collection of the first sample” (days). Additionally, we were also interested in looking at microbial variation across meteorological seasons (quarters) because microbial epidemics and patients usually change their behavior through seasons and because previous studies (e.g., [11, 43, 44]) showed significant associations between microbial variation and season. Hence, we also included meteorological seasons or quarters in our models as either a numerical or a categorical variable. We performed multiple rounds of analysis including all of these factors and the co-variables in Table 1 and S1 Table (age, gender, feeding route, ventilator use, oxygen requirement, tracheostomy change frequency, prophylactic antibiotics and daily inhaled steroids). We also tested the interaction between LRI and Time and between LRI and Quarters. We also compared models with random intercepts and random slopes and the order of our factors. Initial LME models were compared using the function lmerTest, which performs automatic backward elimination of factors. ANOVA type II and III (if interactions were included in the model) tests were also carried out for hypothesis testing. Model assumptions in final LME models were validated using residual vs fit plots and a normal probability plots. Beta-diversity Unifrac indices were compared using permutational multivariate analysis of variance (adonis) as implemented in the vegan R package [45]. Adonis models were compared using the Akaike Index Criterion [46]. We treated our random factors as fixed factors and put them first in the model. Significance was determined through 10,000 permutations. Our preliminary analyses showed that random slopes, Time (alone or interacting with LRI), Quarters (alone or interacting with LRI) and all co-variables in Table 1, except daily inhaled steroids, did not have a significant impact on any representation of microbial diversity or taxon abundance. Hence, our final (most parsimonious) LME and adonis models and analyses included one predictor (LRI) and one co-variable (daily inhaled steroids). Bonferroni or Benjamini-Hochberg FDR multiple test correction methods were applied. All analyses were performed in mothur, QIIME [47], R [48] and RStudio [49].

Table 1

Clinical characteristics of the studied cohort.

NLRI = patients with no lower respiratory infections. YLRI = patients with lower respiratory infections.

Variables	NLRI (N = 20)	YLRI (N = 20)	P-value
Years of Age, median (range)	8.5 (<1–30)	12.5 (<1–34)	0.342
Daily Inhaled Steroids	6	12	0.068
Feeding Route			0.210
G-tube	12	7
GJ-tube	4	2
Oral	2	2
Oral & G-tube	2	9
Gender Male	15	14	0.999
Oxygen Requirement	3	8	0.090
Prophylactic Antibiotics	5	2	0.408
Multiple Tracheostomy Changes	9	6	0.514
Ventilator Use			0.113
Continuous	14	7
Night/Nap	4	11
No	2	2

Clinical characteristics were compared between NLRI and YLRI patients using Fisher’s exact, or Mann-Whitney tests, as appropriate.

Results

Forty patients were enrolled in this study of whom 20 had at least one clinically-evident LRI (YLRI) and 20 had no LRI (NLRI). Clinical characteristics for the study cohort are presented in Table 1 and S1 Table. When “healthy” or at their clinical baseline, we collected 62 tracheal samples from YLRI participants and 65 tracheal samples from NLRI participants. There were 33 (20.6%) tracheal samples missing because families missed sample collection times and because two participants died from their LRIs during the study. All 127 tracheal aspirates were analyzed via MiSeq sequencing of 16S rRNA V4 amplicons. A total of 2,743,049 sequences ranging from 1,514 to 75,719 sequences per sample (mean = 21,598.8; median = 20,727) were obtained after quality control analyses and OTU filtering. From these data, we identified a total of 950 OTUs and 10–154 OTUs (mean = 66.2) per sample (S2 Table; OTU taxa).

The taxonomic composition of the tracheal microbiome

The tracheal microbiome across all 40 patients (127 samples) included sequences that corresponded to the following 18 genera: Streptococcus (16.5%), Neisseria (11%), Haemophilus (8.7%), Moraxella (8.1%), Pseudomonas (7.8%), Corynebacterium (6.6%), Staphylococcus (4.4%), Acinetobacter (3.3%), Prevotella (3%) and Stenotrophomonas (3%). All the other detected genera accounted for <3% of the total sequences. Each of the 127 tracheal microbiomes contained 6 to 18 (mean = 14.8 genera) of these bacterial genera. This assortment of genera reflects membership from both the oral cavity (e.g., Prevotella, Streptococcus and Neisseria) [27, 28, 50] and upper airways (e.g., Streptococcus, Moraxella, Haemophilus, Corynebacterium, and Staphylococcus) [28, 50, 51]. This is not surprising given the anatomic location of the trachea (i.e., between the upper airway and the lung).

Tracheal microbiomes in YLRI patients differ from those in NLRI patients

Microbiotas from YLRI patients showed greater alpha-diversity than microbiotas from NLRI patients (Fig 1A and S1 Fig). These differences resulted significant (0.017≤P≤0.003) for three of the four indices compared in our LME analyses when accounting for steroid use (Table 2). Alpha-diversity indices also varied between YLRI and NLRI patients across quarters, particularly in Q1 (fall) and Q2 (winter) (Fig 1B); however, these differences did not result significant (P>0.05) in our LRI*Quarter LME model and analyses.

Fig 1

Box plots of phylogenetic alpha-diversity of microbiotas from patients with (YLRI) and without (NLRI) lower respiratory infections (LRI) (A) and of microbiotas from YLRI and NLRI patients across meteorological seasons (B).

Table 2

Mean alpha-diversity indices and mean relative proportions of dominant phyla and genera (>3%) in decreasing order of abundance for ALL samples (NLRI+YLRI), patients with no lower respiratory infections (NLRI) and patients with lower respiratory infections (YLRI).

Linear mixed-effects (LME) models results are shown for alpha-diversity indices and taxa proportions, while permutational multivariate analysis of variance (adonis) results are shown for beta-diversity indices. Significance of LME models analyses was estimated using ANOVA type II or III with Satterthwaite approximation. For each test we report the relevant F statistic (F), degrees of freedom (DF) and significance (P(>F)). Significant associations are indicated in bold.

Taxon	ALL	NLRI	YLRI	F	DF	P(>F)
Alpha-diversity
ACE	75.1	64.1	87.1	6.61	39	0.017
Fisher	8.7	7.1	10.4	6.22	40	0.014
PD	6.7	5.9	7.5	9.74	40	0.003
Shannon	2.1	2.0	2.1	0.20	40	0.655
Beta-diversity
Unifrac-unw				2.75	2	0.023
Unifrac-w				7.85	2	0.006
Bray-Curtis				3.40	2	0.011
Jaccard				2.46	2	0.018
Phyla
Proteobacteria	48.8	47.1	50.5	1.28	41	0.265
Firmicutes	26.5	28.7	24.1	1.22	41	0.277
Bacteroidetes	12.5	10.6	14.5	1.11	38	0.298
Actinobacteria	8.3	9.1	7.5	2.02	42	0.163
Fusobacteria	3.2	3.5	2.9	0.65	45	0.424
Genus
Streptococcus	16.5	16.7	16.2	0.69	36	0.411
Neisseria	11.0	9.1	13.6	2.13	37	0.153
Haemophilus	8.7	5.9	11.5	4.13	37	0.009
Moraxella	8.1	9.1	7.0	1.85	36	0.171
Pseudomonas	7.8	9.0	6.5	3.70	40	0.025
Corynebacterium	6.6	7.4	5.6	3.46	41	0.034
Staphylococcus	4.4	6.0	2.7	1.72	34	0.199
Acinetobacter	3.3	4.9	1.5	3.21	37	0.030
Prevotella	3.0	2.6	3.3	0.11	35	0.746
Stenotrophomonas	3.0	2.6	3.5	0.07	37	0.797

PCoA did not reveal clear dissimilarities in beta-diversity between YLRI and NLRI or any other variables in Table 1, since samples were not clearly depicted in discrete groups (see S2 Fig for an example). However, our adonis analyses detected significant differences (0.023≤P≤0.006) in beta-diversity between YLRI and NLRI for all of the four distances when accounting for steroid use. No significant differences (P>0.05) were observed across quarters (LRI*Quarter) in our adonis model.

Phyla abundances varied little between microbiotas from YLRI and NLRI patients (Fig 2A), but genera abundances varied more (Fig 2B). Our LME analyses showed significant associations (0.034≤P≤0.009; Table 2) with LRI for the following four bacterial genera: Haemophilus, Pseudomonas, Corynebacterium and Acinetobacter. Haemophilus was more abundant in YLRI patients than in NLRI patients, while the other three genera showed the inverse relationship. As before, our LME analyses did not detect significant (P>0.05) differences in taxon abundances between YLRI and NLRI samples across meteorological seasons (quarters).

Fig 2

Alluvial plots of mean relative proportions of most abundant (≥3%) phyla and genera in microbiomes from patients with (YLRI) and without (NLRI) lower respiratory infections (LRI) across meteorological seasons.

Discussion

In this study we investigated the composition and temporal dynamics of microbial communities inhabiting the trachea of tracheostomised children and young adults. Our direct tracheal sampling avoids contamination during sampling with microbes from the upper respiratory airways that might be picked up when trying to access the lower airways. We detected changes in the composition and structure of the tracheal microbiotas according to whether or not individuals developed LRI. To our knowledge this is the first study to examine the tracheal microbiome in patients with a tracheostomy. Our results provide a direct assessment of the tracheal microbiome diversity and its temporal dynamics in tracheostomised patients with and without LRIs.

Our diversity analyses show that microbiotas collected at times of health from YLRI patients had significantly higher intra-sample (alpha-diversity) diversity than the microbiotas of patients who did not develop a LRI (NLRI). All YLRI patients acquired their LRIs in the first three quarters (Q1 to Q3), while no LRI occurred during Q4 (summer). If “healthy” microbiomes were fully restored after LRI, one would expect that microbial communities in YLRI and NLRI patients would have similar levels of diversity; our results, however, seem to suggest otherwise. We suspect YLRI healthy samples are comprised of mixed-microbiotas including OTUs prevalent during “infective and healthy times”, while NLRI samples are comprised of OTUs only prevalent during health. In other words, we postulate that LRIs leave a detectable microbial signature.

The mean relative proportions of Haemophilus, Pseudomonas, Corynebacterium, and Acinetobacter varied significantly between microbiotas from YLRI and NLRI patients (Table 2). There are conflicting data about the association between Haemophilus and the frequency and severity of respiratory infections. Indeed, Haemophilus was found by Kloepfer et al. [52] to neither be associated with increased severity of rhinovirus respiratory infections, nor by Carlsson et al. [53] to be associated with the duration of wheezing in young children. By contrast, Teo et al. [43] found that infants with a Haemophilus-dominant nasopharyngeal microbiota had both a higher incidence and higher severity of respiratory infections. Similarly, our group found that infants hospitalized with bronchiolitis who had a Haemophilus-dominant nasopharyngeal microbiota at the time of hospitalization had higher rates of intensive care use and longer length of stay compared with infants who had Moraxella-dominant microbiota [7]. Hence, although there are discrepant findings about Haemophilus’ association with LRI, the weight of the current evidence suggests that the higher abundance of Haemophilus found in the trachea of YLRI patients is indicative of a LRI risk-microbiota [54].

The genus Pseudomonas is of particular importance because it includes several opportunistic human pathogens of clinical relevance. Our metataxonomic (from [55]) analyses detected 1–5 Pseudomonas OTUs (mean relative proportion <8%) of which OTU00005 (S2 Table) accounted for 96.7% of the 156,932 reads and was present in 90% of the patients (in 50% of them at a ≥5% mean relative proportion). BLAST searches identified this OTU as P. aeruginosa and gave it the highest E value (1e-127) in the top 1,000 references. All other Pseudomonas species had marginally lower E values (>1e-127). Although this result suggests that P. aeruginosa is the predominant Pseudomonas species, it must be interpreted with caution since 16S rRNA data cannot distinguish between different strains of P. aeruginosa. Moreover, the presence of P. aeruginosa should be validated by metagenomic approaches such as shotgun sequencing. Culture-based and microbiome analyses of bacterial communities suggest that P. aeruginosa is a common colonizer of the lower and upper respiratory airways in intubated patients [28, 56, 57]. Although P. aeruginosa may also be a devastating pathogen, especially given its frequent multidrug resistance [58]. In our study, a higher abundance of Pseudomonas was associated with a lower risk of developing LRI; in fact, none of the patients showed symptoms of nosocomial infections. Future work is needed to understand how P. aeruginosa transitions from a presumed asymptomatic colonizer of the lower airways to a harmful pathogen and how microbe interactions relate to disease outcomes.

Corynebacterium may be part of a LRI resistance-microbiota [54]. In the present study, individuals with NLRI had significantly higher abundance of Corynebacterium than those with YLRI. This observation is consistent with findings from Biesbroek et al. [26] who found that infants colonized with Corynebacterium (in addition to Moraxella) had a more stable microbiota over the first two years of life and fewer parent-reported respiratory infections. Teo et al. [43] also found that Corynebacterium dominant microbiota in the infant nasopharynx was associated with fewer respiratory infections. Furthermore, Corynebacterium dominance in the nasopharynx has also been associated with a reduced incidence of otitis media [59].

The genus Acinetobacter includes some opportunistic pathogenic species (e.g., A. baumannii) that are becoming increasingly recognized as important in nosocomial infections [60]. Moreover, Acinetobacter species are of increasing concern in association with ventilator associated pneumonia (VAP) [61]. In fact, a previous microbiome study by our group found representatives of this genus in the lung microbiome from mechanically ventilated patients with suspected pneumonia [56]. In our study, however, a higher abundance of Acinetobacter was associated with a lower risk of developing LRI; concordantly, none of the patients showed symptoms of nosocomial infections. Future work is needed to understand what and how Acinetobacter species become pathogenic and their interactions with other microbes during disease.

All these results combined suggest associations between community membership of specific microbial genera and health outcomes, but cause and/or effect has not been established. However, similar findings from multiple studies have helped identify a starting point for functional and/or interventional investigation that will begin unraveling the possible contributions of these bacteria to microbial network dynamics and health outcomes.

The composition (alpha-diversity) and structure (beta-diversity) of the tracheal microbiotas in tracheostomised patients with and without LRI did not change significantly over the course of one year (four meteorological seasons). Similarly, we did not find significant variation in taxon abundance between NLRI and YLRI patients across seasons. Previous studies of the nasopharyngeal microbiomes of healthy [44] and asthmatic [43] infants have revealed seasonal changes in diversity and mean relative proportions of two of the aforementioned genera (Haemophilus and Corynebacterium). A previous study by our group of the dynamics of nasopharyngeal microbiotas in asthmatic children [11] did not find differences in diversity among seasons, but found significant differences in Haemophilus mean relative proportions. Bogaert et al. [44] grouped fall and winter microbial samples and compared them against spring samples (no summer samples were collected), while Teo et al. [43] compared winter-fall to summer-spring samples, although no statistical justification for these groupings was provided. For the sake of comparison, we created two new variables in our study to match the seasonal groupings in those two published studies, while maintaining all the other factors in the model the same. Our LME and adonis analyses detected significant differences (0.04

Our study of the composition and temporal dynamics of microbial communities inhabiting the trachea of tracheostomised children and young adults has one main limitation. We collected samples when patients were healthy since our goal was to determine whether microbiotas from individuals who develop a LRI are different from those who do not develop a LRI. It would be also interesting to compare how the microbiomes change during heath and infection in the same individuals.

Conclusions

We directly characterized the diversity and temporal dynamics of the tracheal microbiota in patients with and without LRIs by bypassing the nose and mouth via a tracheostomy. We demonstrated that the composition and structure of tracheal microbiotas and normalized abundances of Haemophilus, Pseudomonas, Corynebacterium and Acinetobacter differ significantly according to whether individuals developed LRIs. We also showed that tracheal microbiotas’ diversity does not change significantly over meteorological seasons. This result disagrees with previous evidence suggesting seasonal variation in airway microbiotas. Further study is needed to address the interaction between microbes and LRI during times of health and disease.

Clinical and demographic characteristics of the cohort analyzed in this study.

(XLSX)

Click here for additional data file.

Number of cleaned sequences per OTU and sample and corresponding taxonomic identification up to the genus level.

(XLSX)

Click here for additional data file.

Box plots of Shannon, ACE, Fisher and phylogenetic alpha-diversity of microbiotas from patients with (YLRI) and without (NLRI) lower respiratory infections (LRI).

(EPS)

Click here for additional data file.

3D Principal coordinates analyses of unweighted and weighted Unifrac, Bray-Curtis and Jaccard distances between microbiotas from patients with (YLRI) and without (NLRI) lower respiratory infections (LRI).

(EPS)

Click here for additional data file.