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 lifethreatening 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][9][10][11][12][13] and influence both immune [14][15][16][17][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][22][23][24][25][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. 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][31][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 Table 1. Clinical characteristics of the studied cohort. NLRI = patients with no lower respiratory infections. YLRI = patients with lower respiratory infections. 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

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; OTU taxa).

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.
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).

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

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).
https://doi.org/10.1371/journal.pone.0182520.g001 The temporal dynamics of the tracheal microbiome in tracheostomised patients 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.

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 mixedeffects (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. 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<P<0.01) in alpha-and beta-diversity, respectively, between NLRI and YLRI samples across seasonal groups (LRI Ã Seasonal Group), but not between Seasonal Groups alone (e.g., winter-fall and summer-spring). We find this outcome interesting and worth exploring, but since at this point, we do not have statistical support to group individual meteorological seasons, we advise the reader to interpret this result with caution.
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.
Supporting information S1 Table. Clinical and demographic characteristics of the cohort analyzed in this study.