Pyrosequencing Unveils Cystic Fibrosis Lung Microbiome Differences Associated with a Severe Lung Function Decline

Abstract

Chronic airway infection is a hallmark feature of cystic fibrosis (CF) disease. In the present study, sputum samples from CF patients were collected and characterized by 16S rRNA gene-targeted approach, to assess how lung microbiota composition changes following a severe decline in lung function. In particular, we compared the airway microbiota of two groups of patients with CF, i.e. patients with a substantial decline in their lung function (SD) and patients with a stable lung function (S). The two groups showed a different bacterial composition, with SD patients reporting a more heterogeneous community than the S ones. Pseudomonas was the dominant genus in both S and SD patients followed by Staphylococcus and Prevotella. Other than the classical CF pathogens and the most commonly identified non-classical genera in CF, we found the presence of the unusual anaerobic genus Sneathia. Moreover, the oligotyping analysis revealed the presence of other minor genera described in CF, highlighting the polymicrobial nature of CF infection. Finally, the analysis of correlation and anti-correlation networks showed the presence of antagonism and ecological independence between members of Pseudomonas genus and the rest of CF airways microbiota, with S patients showing a more interconnected community in S patients than in SD ones. This population structure suggests a higher resilience of S microbiota with respect to SD, which in turn may hinder the potential adverse impact of aggressive pathogens (e.g. Pseudomonas). In conclusion, our findings shed a new light on CF airway microbiota ecology, improving current knowledge about its composition and polymicrobial interactions in patients with CF.

About the Authors:

Giovanni Bacci

Affiliation: Department of Biology, University of Florence, Florence, Italy

Patrizia Paganin

Affiliation: Department for Sustainability of Production and Territorial Systems, Biotechnologies and Agro-Industry Division, ENEA Casaccia Research Center, Rome, Italy

Loredana Lopez

Affiliation: Department of Energy Technologies, Bioenergy, Biorefinery and Green Chemistry Division, ENEA Trisaia Research Center, Rotondella (MT), Italy

Chiara Vanni

Affiliation: Department of Biology, University of Florence, Florence, Italy

Claudia Dalmastri

Affiliation: Department for Sustainability of Production and Territorial Systems, Biotechnologies and Agro-Industry Division, ENEA Casaccia Research Center, Rome, Italy

Cristina Cantale

Affiliation: Department for Sustainability of Production and Territorial Systems, Biotechnologies and Agro-Industry Division, ENEA Casaccia Research Center, Rome, Italy

Loretta Daddiego

Affiliation: Department of Energy Technologies, Bioenergy, Biorefinery and Green Chemistry Division, ENEA Trisaia Research Center, Rotondella (MT), Italy

Gaetano Perrotta

Affiliation: Department of Energy Technologies, Bioenergy, Biorefinery and Green Chemistry Division, ENEA Trisaia Research Center, Rotondella (MT), Italy

Daniela Dolce

Affiliation: Department of Pediatrics, Cystic Fibrosis Center, Meyer Hospital, Florence, Italy

Patrizia Morelli

Affiliation: Department of Pediatrics, Cystic Fibrosis Center, G. Gaslini Institute, Genoa, Italy

Vanessa Tuccio

Affiliation: Cystic Fibrosis Microbiology and Cystic Fibrosis Center, Children's Hospital and Research Institute Bambino Gesù, Rome, Italy

Alessandra De Alessandri

Affiliation: Department of Pediatrics, Cystic Fibrosis Center, G. Gaslini Institute, Genoa, Italy

Ersilia Vita Fiscarelli

Affiliation: Cystic Fibrosis Microbiology and Cystic Fibrosis Center, Children's Hospital and Research Institute Bambino Gesù, Rome, Italy

Giovanni Taccetti

Affiliation: Department of Pediatrics, Cystic Fibrosis Center, Meyer Hospital, Florence, Italy

Vincenzina Lucidi

Affiliation: Cystic Fibrosis Microbiology and Cystic Fibrosis Center, Children's Hospital and Research Institute Bambino Gesù, Rome, Italy

Annamaria Bevivino

* E-mail: [email protected]

Affiliation: Department of Biology, University of Florence, Florence, Italy

ORCID http://orcid.org/0000-0003-4277-7048

Alessio Mengoni

Affiliation: Department of Biology, University of Florence, Florence, Italy

Introduction

Despite recent improvements in cystic fibrosis (CF) treatments, some individuals with CF still experience a rapid progression of lung disease, which usually leads to irreversible morbidity and mortality. Bacterial airway infections in CF patients are currently monitored through routine microbiology procedures and cures are mainly based on culturable microbial pathogens colonizing CF airways [1]. Unfortunately, the selection and administration of antimicrobial therapies based on the in vitro susceptibility of classic CF pathogens are not necessarily connected with clinical outcomes [2]. In the last years, there has been an increasing awareness that airway colonization is not achieved by a single strain or species but from a complex mixture of microorganisms [3]. Indeed, CF airway microbiota has been shown to be significantly more complex than initially considered [4–7]. Recent works using culture-independent identification methods have revealed that CF patients harbor a vast array of bacterial species, not previously identified and suspected to be involved in CF progression [8–10].

Our understanding of this microbial “black box” has evolved thanks to the development of next-generation sequencing technologies and bioinformatics tools capable of providing new insights into airway microbial communities [11]. Recently, it was proposed that microbial lung community might be considered as a unique, distinct pathogenic entity, whose impact on the host may be greater than the combined effects of its individual component species alone [12]. Even if these findings have drastically altered our understanding of CF lung disease [13,14], their translation into clinical improvements remains a substantial obstacle for enhancing the quality of care [15]. Despite that, a detailed perspective focused on the polymicrobial nature of CF infections would permit to unravel bacterial dynamics providing biomarkers of disease progression, as well as novel bacterial targets for antibiotic treatment [16,15].

Until now, many efforts have been made to expand our understanding of the microbial ecology of the CF airways and new insight into the impact of antibiotic treatment, patient age increasing, and periodic pulmonary exacerbation on CF microbiology have been obtained [8,14,17,18]. Furthermore, a correlation between decreasing diversity of the airways bacterial community and advancing CF lung disease has been found [8] suggesting that reduced diversity may play a role in CF lung disease progression [19]. However, it is still not well clear how the airway microbiota composition changes following a severe decline in lung function, as evidenced by a clinically important drop in FEV₁ (i.e., with a decrease of 5% or more from baseline values). CF airway infection is viewed as an ecosystem where a climax community dominates during relatively stable periods [3]. This ecosystem is dominated by stable populations that are well-adapted to their particular biological niche. A sharp decline in lung function can alter this stable community, possibly shifting it to an alternative state. Our preliminary study, performed by culture-based and culture-independent methods, has provided first insights into relationships between the airway microbiota structure and severe lung function decline [20]. In particular, we identified significant changes in bacterial community diversity when stable (S) and substantial decliner (SD) patients were compared and among patients with different lung disease status (normal/mild, moderate and severe) suggesting that patients with a moderate lung disease (FEV₁ ranging from 40% to 69%) have experienced changes in their airway bacterial communities [20]. However, the lower discrimination power of molecular methods previously used did not permit to clearly address the CF airway microbiota associated with a severe lung function decline. So, in this study, we employed a deep-sequencing approach based on the 16S rRNA gene-targeted analysis to provide a more-in-depth investigation of the airway microbiota of CF patients with stable lung function (S) versus that detected in patients with a sharp decline in lung function (SD). We considered a subset of previously investigated patients (a total of 52 patients) [20] focusing our attention on S and SD patients only, without considering the differences among the different pulmonary physiopathological status (normal/mild, moderate and severe). New insights into bacterial community features such as biological diversity, taxa composition and interaction among taxa were given.

Materials and Methods

Ethics Statement

Protocols for the collection and use of sputum samples from cystic fibrosis patients, and for the procedure of written informed consent were approved by the Ethics Committee of Bambino Gesù Children's Hospital (Rome, Italy), Cystic Fibrosis Center, Meyer Children's Hospital (Florence, Italy) and Giannina Gaslini Children's Hospital (Genoa University, Genoa, Italy) (Prot. N. 681 CM of November 2, 2012; Prot. N. 85 of February 27, 2014; Prot. N. FCC 2012 Partner 4-IGG of September 18, 2012), as previously stated [20]. Informed written consent was obtained from all subjects aged 18 years and over and from parents of all subjects under 18 years of age before enrollment in the study. All sputum specimens were produced voluntarily. The study protocol was in agreement with the Guidelines of the European Convention on Human Rights and Biomedicine for Research in Children and to those of the Ethics Committees of Bambino Gesù, Meyer and Giannina Gaslini Hospitals. All measures were taken to ensure patient data protection and confidentiality.

Clinical characteristics of patients

Fifty-two patients with CF (aged 8–59 years) who regularly attended the three Italian CF Centers were included in this study. Patients were recruited between September 2012 and April 2013 as follows: Bambino Gesù Children's Hospital (Rome, Italy) (n = 29), Giannina Gaslini Children's Hospital (Genoa, Italy) (n = 14) and Meyer Children's Hospital (Florence, Italy) (n = 9). Patients older than six years of age, who had been diagnosed with CF according to the published Guidelines [21], were treated according to current standards of care [22]. Patients were seen in the CF Centres at least four times per year on a regular basis [1]. At each visit, clinical data collection and microbiological status (colonizing germs) were performed [1]. All patients were clinically stable, without any pulmonary exacerbation (as defined by a cluster of symptoms and signs as previously indicated [23,24]) and antibiotic therapy (i.v. or oral) in the previous four weeks before specimen collection. A subset of a previously reported [20] cohort of patients was selected. These patients were classified in our previous work [20] into two groups according to their annual rate of FEV₁ decline by measuring the difference between the best FEV₁% registered within the previous year and the best FEV₁% registered two-years before specimen collection. In total, 29 subjects with a rate decline lower than 1.5% were flagged as “stables” (S), whereas 23 subjects with a rate decline higher than 5%, were flagged as “substantial decliners” (SD). Each FEV₁ value was the average of three repeated measurements obtained 2–3 minutes from each other. FEV₁ values were measured according to the American Thoracic Society (ATS) and the European Respiratory Society (ERS) standards [25]. The overall description of the patient dataset is reported in Table 1.

[Figure omitted. See PDF.]

Table 1. Demographic and clinical characteristics of patients enrolled in the study.

https://doi.org/10.1371/journal.pone.0156807.t001

Sample processing and DNA extraction

Spontaneously expectorated sputum (SES) was used in this study, as it represents by far the most widely used sample in productive patients [26]. Samples processing was performed as previously described [20]. About 400 μl aliquots of frozen sputum were subjected to genomic DNA extraction using the commercially available Kit QIAamp DNA Mini Kit. Sample aliquots were spun at 10,000×g to pellet cellular material. After removal of the supernatant, cell pellets were resuspended in 180 μl of the appropriate enzyme solution (20 mg/ml lysozyme; 20 mM Tris·HCl, pH 8.0; 2 mM EDTA; 1.2% Triton), incubated for 30 min at 37°C and then processed according to the manufacturer’s protocol. Quantity and purity of extracted DNA were checked by NanoDrop (NanoDrop Technologies, USA), PicoGreen fluorescent assay (Life Technologies), and gel electrophoresis.

DNA amplification and sequencing

DNA samples were subjected to 16S rRNA gene amplification. Primers and barcodes were designed according to the Human Microbiome Project Consortium (HMP) (http://www.hmpdacc.org/tools_protocols/tools_protocols.php). The V3, V4, and V5 hypervariable regions of the 16S rRNA gene were amplified using primers 357F (5′-CCTACGGGAGGCAGCAG-3′) and 926R (5′-CCGTCAATTCMTTTRAGT-3′) modified with the addition of the 454 FLX-titanium adaptors B and A, respectively. Unique 7-nucleotide barcode sequences (MID 1–18) were added to the A adaptor as reported in S1 Table. PCR amplification was performed on 2 μl of DNA template in a total volume of 25 μl containing 1× AccuPrime Buffer II (Life Technologies), 10 μM of 357F fusion-primer, 10 μM of 926F fusion-primer and 0.03 U/μl AccuPrime High-Fidelity Taq DNA polymerase (Life Technologies). PCR reactions were heated at 95°C for 2 min followed by 30 cycles of 95°C for 20 seconds, 55°C for 30 seconds, and 72°C for 5 minutes. All reactions were prepared in a sterile PCR hood. Three independent PCR reactions were performed for each DNA-template. Negative control reactions were also performed. Each replicate reaction was examined by electrophoresis on 1.5% agarose gel. PCR amplicons were cleaned using the Agencourt AMPure XP Beads according to the manufacturer’s specifications (Beckman Coulter). To ensure removal of primers and any non-specific amplicons, purified amplicon libraries were analyzed using the Agilent Bioanalyzer 2100 employing the Agilent DNA 1000 Kit. The Quant-iT PicoGreen dsDNA fluorescent assay Kit (Life Technologies) was employed to establish the concentration of the purified amplicons. Three independent purified and quantified PCR reactions were pooled in equimolar proportion foe each sample. Pools were quantified using KAPA Library Quantification Kits (KAPA Biosystems) to determine the number of amplifiable molecules in the libraries. Emulsion PCR, emulsion breaking, and amplicon pyrosequencing were performed applying the 454 GS FLX+ chemistry following supplier protocols (454 Life Sciences Roche Corporation). The GS FLX+ software 2.9 version was used for sequencing and the pipeline 3 for long amplicon was used to data processing (454 Life Sciences Roche Corporation). Sequence files have been uploaded to the NCBI Sequence Read Archive (SRA) with the accession number: PRJNA290694.

Sequence processing and data generation

Amplicon sequences were demultiplexed using Mothur [27], discarding those with mismatches in the barcode, primers, linkers or spacers. A first quality check step was performed using StreamingTrim [28] to remove ambiguous base calls (Phred’s quality score < 25) and convert sequences into FASTQ format. Filtered sequences were processed following the UPARSE pipeline included in the USEARCH package (version 7.0.1090) [29]. In particular, sequences were trimmed to a fixed length of 250 base-pairs (bp) through the “fastq_trunclen” command. Chimeras were removed from pooled sequences with the “uchime_denovo” algorithm, and Operational Taxonomic Units (OTUs) were generated with an identity cutoff of 97% (which is approximately close to the taxonomic level of species according to [30]. SINA standalone classifier was used for taxonomic assignment of representative sequences in combination with the most recent version of SILVA non-redundant database available (SSU Ref NR 99, release 115). Genera were categorized as aerobic/anaerobic and pathogenic/non-pathogenic as reported in [10,31]. OTU oligotyping [32] was performed following the procedures reported in http://merenlab.org/2013/11/04/oligotyping-best-practices/ with a minimum substantive abundance criterion (M) of 50. We considered an OTU fully resolved when all its olygotypes had a purity index of at least 90%. For each oligotype, the first 5 BLAST best hits [33] corresponding to a defined taxonomy with 100% identity and 100% coverage were reported.

Statistical analyses

To analyze FEV₁ decline between patient groups, mixed-effect models with random intercept and slope were used as data contained repeated measures of the same variable (FEV₁ index) for each unit (patient) in the dataset. Shapiro-Wilk test was used to determine to determine if data were normally distributed. Thus different linear mixed-effect models were compared through the likelihood ratio test choosing the one which best fitted our data; for additional information about fitted models see S1 Appendix and S1 Fig. Mixed-effect models were fitted and compared in R version 3.2.5 using the lme4 [34] package version 1.1.

Bacterial diversity was inspected through Shannon, Chao 1, Evenness and Richness indices (the latter expressed as the count of OTUs observed in each sample) on both rarefied and non-rarefied counts. Since differences between rarefied and non-rarefied values were really low (less than 0.01%, data not shown), non-rarefied values were reported. All indices were computed with the R [35] package vegan [36]. The percentage of coverage was calculated by Good's estimator [37] using the formula: [1—(n/N)] × 100, where n is the number of sequences found once in a sample (singletons), and N is the total number of sequences in that sample. The Evenness index was calculated using the formula E = S/log(R), where S is the Shannon diversity index and R is the number of OTUs in the sample (the Richness). Bacterial communities from different groups of patients were compared using Multivariate Analysis of Variance (MANOVA) with Pillai-Bartlett statistic. Single OTUs and oligotypes were compared through the Analysis of Variance (ANOVA) in combination with the Tukey's honest significant difference (HSD) test. Canonical Correlation Analysis (CCA) was performed using the log transformed OTUs abundance. All statistical analyses and graphical representations were carried out using the R software.

Network construction

Co-occurrence networks were constructed following previously reported methods [38,39] Undirected graphs reporting the correlation between OTUs were produced. Each node corresponds to a different OTU, whereas each edge represents a significant (p < 0.05) Spearman’s correlation. The diameter of the nodes was directly proportional to the OTU occurrence in the dataset whereas the width of the edges was directly proportional to the Spearman’s correlation index. Produced graphs were arranged using the Fruchterman Reingold algorithm [40]. All statistical analyses were carried out in the R environment using the igraph package for network generation [41].

Results

Airways bacterial community structure in S and SD patients

Sputum DNA from 52 specimens, obtained during routine medical care, in accordance with the ethical Guidelines, was analyzed by pyrosequencing of the V3–V5 hypervariable region of the 16S rRNA gene. A total of 201903 reads were generated from these samples with an average length of 541 bp. After demultiplexing, quality control, chimera detection, and length filtering, 158107 high-quality sequences of 250 bp each were collected with an average of 3040 sequences per sample (ranging from 880 to 6381). The number of sequences and OTUs detected for each patient has been reported in S2 Table.

To investigate the depth of our 16SrRNA amplicon libraries, rarefaction curves were computed on OTU assignments. Most rarefaction curves approached asymptote (S2 Fig), suggesting that our sampling efforts were able to represent bacterial community assemblages of all patients. Moreover, Good's estimator reported a satisfactory coverage (higher than 99% for all 52 samples, S2 Table). Forty-four OTUs (based on 97% sequence similarity) were detected with each sample ranging from 10 to 35 OTUs (S2 Table) (median = 25). Twenty-two genera were detected with a high intra- and inter- groups variance (Fig 1). All the genera recovered so far from the respiratory tracts of patients with CF were detected, with the only exception of the anaerobic genus Sneathia. The analysis of variance of diversity indices (Shannon index, Chao1 and Evenness) based on OTU counts, showed that different conditions did not alter the bacterial diversity of the lung microbiota which remained almost constant in each group (S3 Fig).

[Figure omitted. See PDF.]

Fig 1. Relative abundance of sputum microbiota among patients.

Relative OTUs abundance was computed dividing the number of 16S sequences assigned to each OTU by the total number of sequences obtained for each sample. Boxes denote the interquartile range (IQR) between the 25^th and the 75^th percentile (first and third quartiles), whereas the inner line represents the median. Whiskers represent the lowest and highest values within 1.5 times IQR from the first and third quartiles Outliers were reported using white circles. CP: common CF pathogens; NC: non-classical but commonly identified genera in CF; OG: other genera in CF; NDG: not yet described genera in CF; S: stable patients; SD: severe decliner patients; [un] unclassified; [unc] uncultured.

https://doi.org/10.1371/journal.pone.0156807.g001

OTU 5 was the most abundant (35052 sequences), and it was the only one classified as Pseudomonas. Moreover, Pseudomonas was the most abundant genus in 33% of patients, followed by Staphylococcus (OTU 2, 12990 reads in total and the most abundant genus in 19% of patients). Ten OTUs were classified as Prevotella, which was found in 96% of the subjects with a total of 16353 sequences (ranging from 0.02% to 74.49% of the total sequences in each sample). In one sample, Prevotella was the dominant genus followed by Staphylococcus (8.5% of sequences). Two OTUs could not be identified at the genus level, and they were reported using their family attribution with the unclassified flag “un” (Carnobacteriaceae, 5143 sequences in 43 samples and Alcaligenaceae, 741 sequences in 4 samples). Only one OTU was classified as an “unculturable organism belonging to the Ruminococcaceae family” (18 sequences in 6 samples). Interestingly, when considering genus abundance, Pseudomonas showed a significant negative correlation with microbial diversity in both groups of patients (Shannon index: r < -0.4, p < 0.05; Chao 1 index: r < -0.4, p < 0.05; Richness index: r < -0.5, p < 0.01, Spearman's rank correlation coefficient for both S and SD group). The abundance of aerobic and anaerobic genera did not show significant differences between the two groups of patients. Similarly, the abundance of the most common CF pathogens did not seem to change in S and SD groups (Wilcoxon signed-rank test: p > 0.05 for all contrasts).

The relationship between global bacterial community composition and patient groups was then inspected using CCA and two-way MANOVA. The CCA analysis (which explained the 15.7% of the overall inertia, CCA1: 1.6%; CA1: 14.1%), indicated a slight overlap between patient groups (Fig 2). However, a significant correlation was found between the ordination analysis and patients’ conditions (Environmental fitting based on 1 000 permutations: p < 0.01).

[Figure omitted. See PDF.]

Fig 2. Canonical Correlation analysis based on the log-transformed abundance of 44 OTUs.

Patient conditions were used as constraining variable. Two first components (CCA1 and CA1) were plotted accounting for 15.6% of overall inertia of the data set. Individuals (represented by points) were clustered, and centroids were computed for each group. Ellipses were drawn using the standard deviation of points belonging to the same cluster with a confidence limit of 95%.

https://doi.org/10.1371/journal.pone.0156807.g002

The presence of a microbiota signature of patients groups was also supported by two-way MANOVA (Pillais’ Trace = 0.99, p < 0.05; S3 Table). In particular, two OTUs (OTU 2, classified as Staphylococcus and OTU 36, classified as Prevotella) showed a different distribution between groups of patients (ANOVA: p < 0.05; Tukey HSD: p < 0.05; Fig 3a).

[Figure omitted. See PDF.]

Fig 3. Differences in microbial community distribution between stable and severe decliner groups.

Bacterial taxonomic profiling revealed differences between stable and severe decliner communities. a) Standardized abundances of two OTUs (OTU 2 and 36) showing a different distribution between patient groups. b) Differences between the two most represented oligotypes identified for OTU 2. Standardize abundances were calculated as: [x - mean(x)]/sd(x), where "sd" is the standard deviation and “mean” is the mean value. As a result, all OTUs have equal means and standard deviations (0 and 1, respectively) but different ranges. Reported p values were obtained with a Wilcoxon signed-rank test.

https://doi.org/10.1371/journal.pone.0156807.g003

OTUs were further characterized by oligotype analysis aiming to detect the presence of different species/strains within the same genus. For OTU 2 (classified as Staphylococcus), by using five entropy locations (68^th, 73^th, 80^th, 175^th and 178^th position), we identified three oligotypes: ACTTC (8 218 reads in 36 subjects), ATTCC (1278 reads in 33 subjects) and ATTCT (104 reads in one subject) (Table 2).

[Figure omitted. See PDF.]

Table 2. BLAST analysis of OTU 2 and OTU 36 oligotypes.

https://doi.org/10.1371/journal.pone.0156807.t002

The first oligotype (the most abundant) was the only one to show a different distribution between S and SD patients (Wilcoxon signed-rank test: p < 0.05), as reported in Fig 3b. BLAST analysis reported different taxonomic annotations for the three oligotypes, with the ACTTC oligotype matching Staphylococcus aureus sequences in the first three BLAST best hits (Table 2). OTU 36 was totally resolved with three entropy locations (12^th, 196^th and 246^th positions) accounting for a single olygotype. Its taxonomic attribution was confirmed by BLAST analysis, which reported three different members of Prevotella genus (Table 2). In addition to OTU 2 and OTU 36, we performed oligotype analysis on the OTU 5. This OTU was the only one assigned to the Pseudomonas genus, which has been correlated negatively with bacterial diversity as reported above. Oligotype analysis produced ten different oligotypes most of which reporting P. aeruginosa as their first BLAST hit (for additional details see S2 Appendix).

Network analysis

Aiming to investigate if patients’ groups may reflect ecological niche processes at the airway level that drive coexistence within microbiota, co-occurrence relationships were analyzed. Four networks were produced reporting positive and negative correlations between OTUs detected in S and SD patients (Fig 4).

[Figure omitted. See PDF.]

Fig 4. OTU Networks based on correlation analysis.

Networks report positive and negative correlations between OTUs found in the airway microbiota of S and SD patients. Node color corresponds to taxonomic assignments, whereas node size reflects the log-transformed abundance. Edge thickness is proportional to the modulus of Spearman's rank correlation coefficient. Average degree values are reported in the upper-right corner of each network.

https://doi.org/10.1371/journal.pone.0156807.g004

The positive correlation could indicate cooperative interactions or the presence of common biological functions or ecological niche between taxa. Negative correlations could be indicative of either competitive interactions or non-overlapping ecological niches or processes between taxa [42]. Spearman's rank correlation coefficient obtained in our dataset ranged from 0.37 to 0.81 and from -0.42 to -0.70 in positive and negative networks, respectively. Considering global network properties, positive correlation networks showed a higher number of edges than the negative ones. Interestingly, the networks of SD patients showed lower average degree values in both positive and negative correlation on those of S patients. Finally, Pseudomonas showed no positive correlations with any other member of the microbial community of both S and SD patients but represented the node with the highest number of links in the negative correlation networks (with 18 and 11 edges for S and SD group, respectively).

Discussion

During the last five years, the approach to the study of CF airway infections has changed. The analysis of single pathogens has been overcome by the analysis of lung bacterial community as a whole. Indeed, several studies have shown that bacterial community in CF lung is very heterogeneous and that its variability may be influenced by multiple factors, such as patient age, exacerbations status, antibiotic treatment and pulmonary function decline [5,8,18,43]. Given the importance of lung function in CF patients’ health, it is by extension important to understand the complexity of CF microbiota especially in patients displaying fluctuating levels of lung function, identifying ecological factors associated with higher/lower pulmonary function decline. As we continue to enhance our knowledge of the airway microbiota, these factors and other unanswered questions will guide future research efforts directed towards understanding the complex interactions governing the host–microorganism relationship.

In a previous work [20] we provided first insights into the CF airway microbiota of S and SD patients. By combining culture-based methods and Terminal Restriction Fragment Length Polymorphism (T-RFLP) analysis of amplified 16S rRNA gene from sputum with clinically relevant information, we found that the bacterial community in SD patients had clear differences in alpha diversity parameters. However, due to the limitation of T-RFLP in taxa identification, we did not detail which bacterial taxa mainly contribute to determining such differences between S and SD patients’ airways bacterial communities, nor we did provide a clearer picture regarding the dynamic interactions between bacterial communities. The approach used in this study allowed us to assess exhaustively the bacterial composition in CF samples showing the presence of 22 distinct bacterial genera, with a massive presence of Pseudomonas and Staphylococcus representatives. In particular, Staphylococcus aureus was more abundant in S patients than in SD ones. Interestingly, we were able to identify several members of Prevotella genus, which is considered as an emerging CF pathogens [43,44]. Moreover, thanks to the culture-free approach here adopted we were able to detect an unclassified member of Alcaligenaceae (a family in the order Burkholderiales of the β-Proteobacteria [45]) and Carnobacteriaceae families along with the unusual anaerobic genus Sneathia. Representatives of the first family have been isolated from various clinical samples, including respiratory secretions of CF patients [46]. Among them, Achromobacter xylosoxidans represents one of the most important emerging CF pathogens [47]. The recovery of organisms that are not reported by using the routine clinical protocols such as members of the Carnobacteriaceae were previously reported by expanded culture-dependent profiling of CF airway microbiology [48]. Conversely, to the best of our knowledge, members of Sneathia have never been associated with CF airways microflora. However, reports of Sneathia occurrence in lung transplant are present [49], suggesting that this unusual anaerobic genus could indeed colonize lung tissue. It is worth mentioning that, due to the resolution limit of 16S rRNA based analyses, further taxonomic characterization was not possible for members of these families within our subjects. In general, results obtained allowed us to confirm the value of culture-based methods. However, the most commonly identified non-classical genera in CF, as well as other minor genera described in CF have been detected by pyrosequencing analysis, highlighting the polymicrobial nature of CF infection [7,50], but their relation with lung decline deserves further attention and investigation.

The dynamic balance in the airways microbiota ecosystem depends on interactions among bacterial species as well as between bacteria and the host environment [51]. Until now, microbe-microbe interactions in CF have been evaluated under in vitro conditions, focused mainly on the traditional pathogen P. aeruginosa, and have revealed the interplay among the classical CF pathogens [50]. Applying established ecological principles to the study of CF airway infection has produced important biological and pathophysiological insights into CF airway infection [3]. Network analysis of taxon co-occurrence patterns gives new insight into the structure of complex microbial communities [39] and offers one of the best apparatus available today to tackle the tangled human microbiota [51]. Positive correlation networks revealed widely associated bacterial communities, showing no connections with Pseudomonas. A significant negative correlation was found between bacterial diversity and the presence of members of Pseudomonas genus. Furthermore, negative correlation networks showed a massive number of negative correlations between Pseudomonas and other members of lung microbiota, suggesting a strong antagonistic potential of Pseudomonas over the airway microbiota and the other CF pathogens [52]. Interestingly, the network degree values (i.e. the number of interaction each taxon has with the others) is higher in S than in SD group, suggesting that the S airway microbiota could be more resilient to environmental changes (as for instance those due to the spread of a novel opportunistic pathogen, such as Pseudomonas). Indeed, a direct link between degree of the correlation network and community resilience has been shown in other systems [53–55]. We can consequently speculate that airway bacterial community, as a whole, can confer a benefit to S patients with respect to the perturbation caused or associated with the spread of an opportunistic pathogen (e.g. Pseudomonas). The interactions within the microbiota network can influence the healthy or disease status of the human body. The loss of any interconnections between bacteria belonging to Pseudomonas genus and the OTUs co-occurrence offers a piece of concrete evidence to support the role of this important pathogen in CF disease. However, diversity measurements are not sufficient to predict patient status in a cross-sectional study [44]. Indeed, longitudinal studies, following a cohort of patients over time, may help to elucidate the microbial factors that can contribute to their status changes.

In conclusion, we provided an in-depth description of bacterial community composition based on the analysis of 16S rRNA gene sequence at different resolution levels from genera to olygotype characterization. Our findings revealed a different structure and composition of airway microbiota in CF airways of S and SD patients. The severe decline in lung function experienced by CF patients was associated with a lower complexity of microbial community, which in turn could have paved the way for the proliferation of pathogenic bacteria (such as members of Staphylococcus genus), whereas the stable lung function condition favours a relatively complex bacterial community. The correlation network analysis showed a decreased number of positive correlations in SD patients microbiota with respect with S ones. In such analysis Pseudomonas emerged as the most negatively connected genus, suggesting the presence of a high number of competitive interactions between Pseudomonas and the other taxa. A high presence of Pseudomonas would consequently reduce the number of other taxa inhabiting the airways bacterial community, and then the overall community diversity. Cross-sectional analysis is, of course, unable to show how community dynamics change, as a result of immune response or treatment. This is why a metagenomic longitudinal analysis is mandatory as it might lead to the identification of predictors of clinical change. The possibility to analyze the meta-community dynamics and to identify signatures for S and SD patients can give us a set of tools to unlock the potential of microbiome-based personalized medicine in major disease areas including CF.

Supporting Information

[Figure omitted. See PDF.]

S1 Appendix. Supplementary Methods.

https://doi.org/10.1371/journal.pone.0156807.s001

(DOCX)

S2 Appendix. Supplementary Results.

https://doi.org/10.1371/journal.pone.0156807.s002

(DOCX)

S1 Fig. FEV₁ values for each patient in the two years before the enrollment.

Enrollment date is reported next to the subject ID whereas the sampling year is reported in each x-axis. Blue lines represent predicted values based on mixed-effect models whereas red crosses were used to mark the best values in each sampling year. The difference between the best FEV₁ values registered one year and two years before the enrollment was reported with dashed lines.

https://doi.org/10.1371/journal.pone.0156807.s003

(PDF)

S2 Fig. Rarefaction curves based on the number of OTUs found in each sample.

https://doi.org/10.1371/journal.pone.0156807.s004

(PDF)

S3 Fig. Diversity indices in the two groups of patience observed.

Each p-values obtained through the analysis of variance (ANOVA) was reported along with the name of the index analyzed. Boxes denote the interquartile range (IQR) between the 25^th and the 75^th percentile (first and third quartiles), whereas the inner line represents the median. For a more detailed description about boxes see Fig 1.

https://doi.org/10.1371/journal.pone.0156807.s005

(PDF)

S1 Table. Primers and barcodes used for 16S rRNA gene amplifications and sequencing.

https://doi.org/10.1371/journal.pone.0156807.s006

(DOCX)

S2 Table. Number of sequences, OTUs, diversity indices and Good’s coverage estimator for all patients enrolled in the study.

https://doi.org/10.1371/journal.pone.0156807.s007

(DOCX)

S3 Table. Multivariate analysis of variance (MANOVA) on OTUs dataset.

https://doi.org/10.1371/journal.pone.0156807.s008

(DOCX)

Acknowledgments

We thank the Italian Cystic Fibrosis Research Foundation (FCC) for its support and administrative tasks, and Gabriella Ricciotti and Silvia Campana for technical assistance.

Author Contributions

Conceived and designed the experiments: AB VL GT EVF AM. Performed the experiments: GB PP LL CV LD. Analyzed the data: GB AM AB. Contributed reagents/materials/analysis tools: CD CC LL LD GP VT PM DD. Wrote the paper: GB AM AB. Provided comments and recommendations that improved the manuscript: AM GT CD EVF VL AB. Supervised research: AB VL EVF GT ADA AM.

Citation: Bacci G, Paganin P, Lopez L, Vanni C, Dalmastri C, Cantale C, et al. (2016) Pyrosequencing Unveils Cystic Fibrosis Lung Microbiome Differences Associated with a Severe Lung Function Decline. PLoS ONE 11(6): e0156807. https://doi.org/10.1371/journal.pone.0156807

References

1. Smyth AR, Bell SC, Bojcin S, Bryon M, Duff A, Flume P, et al. European cystic fibrosis society standards of care: Best practice guidelines. J Cyst Fibros. 2014;13.

2. Sibley CD, Parkins MD, Rabin HR, Surette MG. The relevance of the polymicrobial nature of airway infection in the acute and chronic management of patients with cystic fibrosis. Curr Opin Investig Drugs. 2009;10: 787–794. pmid:19649923

3. Conrad D, Haynes M, Salamon P, Rainey PB, Youle M, Rohwer F. Cystic fibrosis therapy: a community ecology perspective. Am J Respir Cell Mol Biol. 2013;48: 150–6. pmid:23103995

4. Lipuma JJ. The changing microbial epidemiology in cystic fibrosis. Clin Microbiol Rev. 2010;23: 299–323. pmid:20375354

5. Lynch SV, Bruce KD. The cystic fibrosis airway microbiome. Cold Spring Harb Perspect Med. 2013;3: a009738. pmid:23457293

6. Rogers GB, Shaw D, Marsh RL, Carroll MP, Serisier DJ, Bruce KD. Respiratory microbiota: addressing clinical questions, informing clinical practice. Thorax. 2014; 1–8.

7. Huang YJ, LiPuma JJ. The Microbiome in cystic fibrosis. Clin Chest Med. 2015;37: 59–67. pmid:26857768

8. Zhao J, Schloss PD, Kalikin LM, Carmody LA, Foster BK, Petrosino JF, et al. Decade-long bacterial community dynamics in cystic fibrosis airways. Proc Natl Acad Sci U S A. 2012;109: 5809–14. pmid:22451929

9. Parkins MD, Floto RA. Emerging bacterial pathogens and changing concepts of bacterial pathogenesis in cystic fibrosis. J Cyst Fibros. 2015;14: 293–304. pmid:25881770

10. Bittar F, Rolain J-M. Detection and accurate identification of new or emerging bacteria in cystic fibrosis patients. Clin Microbiol Infect. 2010;16: 809–20. pmid:20880410

11. Delmont TO, Malandain C, Prestat E, Larose C, Monier J-M, Simonet P, et al. Metagenomic mining for microbiologists. ISME J. International Society for Microbial Ecology; 2011;5: 1837–43.

12. Rogers GB, Carroll MP, Hoffman LR, Walker AW, Fine DA, Bruce KD. Comparing the microbiota of the cystic fibrosis lung and human gut. Gut Microbes. 2010; 85–93. pmid:21326915

13. Carmody LA, Zhao J, Schloss PD, Petrosino JF, Murray S, Young VB, et al. Changes in Cystic Fibrosis Airway Microbiota at Pulmonary Exacerbation. Ann Am Thorac Soc. American Thoracic Society; 2013;10: 179–187.

14. Zemanick ET, Harris JK, Wagner BD, Robertson CE, Sagel SD, Stevens MJ, et al. Inflammation and airway microbiota during cystic fibrosis pulmonary exacerbations. PLoS One. 2013;8.

15. Filkins LM, O’Toole GA. Cystic fibrosis lung infections: polymicrobial, complex, and hard to treat. PLoS Pathog. 2015;11: e1005258. pmid:26719892

16. Sibley CD, Parkins MD, Rabin HR, Duan K, Norgaard JC, Surette MG. A polymicrobial perspective of pulmonary infections exposes an enigmatic pathogen in cystic fibrosis patients. Proc Natl Acad Sci U S A. 2008;105: 15070–5. pmid:18812504

17. Cox MJ, Allgaier M, Taylor B, Baek MS, Huang YJ, Daly RA, et al. Airway microbiota and pathogen abundance in age-stratified cystic fibrosis patients. PLoS One. 2010;5: e11044. pmid:20585638

18. Zhao J, Murray S, Lipuma JJ. Modeling the impact of antibiotic exposure on human microbiota. Sci Rep. 2014;4: 4345. pmid:24614401

19. VanDevanter DR, LiPuma JJ. Microbial diversity in the cystic fibrosis airways: where is thy sting? Future Microbiol. 2012;7: 801–803. pmid:22827300

20. Paganin P, Fiscarelli EV, Tuccio V, Chiancianesi M, Bacci G, Morelli P, et al. Changes in cystic fibrosis airway microbial community associated with a severe decline in lung function. PLoS One. 2015;10: e0124348. pmid:25898134

21. Farrell PM, Rosenstein BJ, White TB, Accurso FJ, Castellani C, Cutting GR, et al. Guidelines for diagnosis of cystic fibrosis in newborns through older adults: Cystic Fibrosis Foundation consensus report. J Pediatr. 2008;153: S4–S14. pmid:18639722

22. Kerem E, Conway S, Elborn S, Heijerman H, Committee C, others. Standards of care for patients with cystic fibrosis: a European consensus. J Cyst Fibros. 2005;4: 7–26.

23. Fuchs HJ, Borowitz DS, Christiansen DH, Morris EM, Nash ML, Ramsey BW, et al. Effect of aerosolized recombinant human DNase on exacerbations of respiratory symptoms and on pulmonary function in patients with cystic fibrosis. The Pulmozyme Study Group. N Engl J Med. 1994;331: 637–42. pmid:7503821

24. Ramsey BW, Pepe MS, Quan JM, Otto KL, Montgomery AB, Williams-Warren J, et al. Intermittent administration of inhaled tobramycin in patients with cystic fibrosis. Cystic Fibrosis Inhaled Tobramycin Study Group. N Engl J Med. 1999;340: 23–30. pmid:9878641

25. Flume PA, O’Sullivan BP, Robinson K a, Goss CH, Mogayzel PJ Jr, Willey-Courand DB, et al. Cystic fibrosis pulmonary guidelines: chronic medications for maintenance of lung health. Am J Respir Crit Care Med. 2007;176: 957–969. pmid:17761616

26. Rogers GB, Skelton S, Serisier DJ, van der Gast CJ, Bruce KD. Determining cystic fibrosis-affected lung microbiology: comparison of spontaneous and serially induced sputum samples by use of terminal restriction fragment length polymorphism profiling. J Clin Microbiol. 2010;48: 78–86. pmid:19906901

27. Schloss PD, Westcott SL, Ryabin T, Hall JR, Hartmann M, Hollister EB, et al. Introducing mothur: open-source, platform-independent, community-supported software for describing and comparing microbial communities. Appl Environ Microbiol. 2009;75: 7537–7541. pmid:19801464

28. Bacci G, Bazzicalupo M, Benedetti A, Mengoni A. StreamingTrim 1.0: A Java software for dynamic trimming of 16S rRNA sequence data from metagenetic studies. Mol Ecol Resour. 2014;14: 426–434. pmid:24128146

29. Edgar RC. UPARSE: highly accurate OTU sequences from microbial amplicon reads. Nat Methods. 2013;10: 996–8. pmid:23955772

30. Konstantinidis KT, Tiedje JM. Prokaryotic taxonomy and phylogeny in the genomic era: advancements and challenges ahead. Current Opinion in Microbiology. 2007;10: 504–509. pmid:17923431

31. Rogers GB, Cuthbertson L, Hoffman LR, Wing PAC, Pope C, Hooftman DAP, et al. Reducing bias in bacterial community analysis of lower respiratory infections. ISME J. International Society for Microbial Ecology; 2013;7: 697–706.

32. Eren AM, Borisy GG, Huse SM, Mark Welch JL. Oligotyping analysis of the human oral microbiome. Proc Natl Acad Sci U S A. 2014;111: E2875–84. pmid:24965363

33. Camacho C, Coulouris G, Avagyan V, Ma N, Papadopoulos J, Bealer K, et al. BLAST+: architecture and applications. BMC Bioinformatics. 2009;10: 421. pmid:20003500

34. Bates D, Mächler M, Bolker B, Walker S. Fitting linear mixed-effects models using lme4. J Stat Softw. 2015;67: 1–48.

35. R Development Core Team. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. ISBN 3-900051-07-0, URL http://www.R-project.org/. R Found Stat Comput Vienna, Austria. 2012;

36. Oksanen J, Blanchet F, Kindt R, Legendre P, Minchin P, O’Hara R, et al. Vegan: community ecology package. R package version 2.0–10. R package version. 2013. p. 10.4135/9781412971874.n145

37. Good IJ. The population frequencies of species and the estimation of population parameters. Biometrika. 1953;40: 237–264.

38. Williams RJ, Howe A, Hofmockel KS. Demonstrating microbial co-occurrence pattern analyses within and between ecosystems. Front Microbiol. 2014;5: 358. pmid:25101065

39. Barberán A, Bates ST, Casamayor EO, Fierer N. Using network analysis to explore co-occurrence patterns in soil microbial communities. ISME J. 2012;6: 343–51. pmid:21900968

40. Fruchterman TMJ, Reingold EM. Graph drawing by force-directed placement. Softw Pract Exp. 1991;21: 1129–1164.

41. Csardi G, Nepusz T. The igraph software package for complex network research. InterJournal. 2006; Complex Sy: 1695.

42. Faust K, Raes J. Microbial interactions: from networks to models. Nat Rev Microbiol. 2012;10: 538–50. pmid:22796884

43. Mahenthiralingam E. Emerging cystic fibrosis pathogens and the microbiome. Paediatr Respir Rev. 2014;15: 13–15.

44. Surette MG. The cystic fibrosis lung microbiome. Ann Am Thorac Soc. 2014;11: 61–65.

45. Rosenberg E, DeLong EF, Lory S, Stackebrandt E, Thompson F, editors. The Prokaryotes. Berlin, Heidelberg: Springer Berlin Heidelberg; 2014. https://doi.org/10.1007/978-3-642-30197-1

46. Coenye T, Vanlaere E, Samyn E, Falsen E, Larsson P, Vandamme P. Advenella incenata gen. nov., sp. nov., a novel member of the Alcaligenaceae, isolated from various clinical samples. Int J Syst Evol Microbiol. 2005;55: 251–256. pmid:15653883

47. Jakobsen TH, Hansen MA, Jensen PØ, Hansen L, Riber L, Cockburn A, et al. Complete genome sequence of the cystic fibrosis pathogen Achromobacter xylosoxidans NH44784-1996 complies with important pathogenic phenotypes. PLoS One. 2013;8.

48. Sibley CD, Grinwis ME, Field TR, Eshaghurshan CS, Faria MM, Dowd SE, et al. Culture enriched molecular profiling of the cystic fibrosis airway microbiome. PLoS One. 2011;6: e22702. pmid:21829484

49. Bittinger K, Charlson ES, Loy E, Shirley DJ, Haas AR, Laughlin A, et al. Improved characterization of medically relevant fungi in the human respiratory tract using next-generation sequencing. Genome Biol. 2014;15: 487. pmid:25344286

50. Magalhães AP, Azevedo NF, Pereira MO, Lopes SP. The cystic fibrosis microbiome in an ecological perspective and its impact in antibiotic therapy. Appl Microbiol Biotechnol. 2016;100: 1163–81. pmid:26637419

51. Sam Ma Z, Guan Q, Ye C, Zhang C, Foster JA., Forney LJ. Network analysis suggests a potentially ‘evil’ alliance of opportunistic pathogens inhibited by a cooperative network in human milk bacterial communities. Sci Rep. 2015;5: 8275. pmid:25651890

52. Smith DJ, Badrick AC, Zakrzewski M, Krause L, Bell SC, Anderson GJ, et al. Pyrosequencing reveals transient cystic fibrosis lung microbiome changes with intravenous antibiotics. Eur Respir J. 2014;44: 922–30. pmid:25034564

53. Shade A, Peter H, Allison SD, Baho DL, Berga M, Bürgmann H, et al. Fundamentals of microbial community resistance and resilience. Front Microbiol. 2012;3: 417. pmid:23267351

54. Wang Y, Zhang R, Zheng Q, Deng Y, Van Nostrand JD, Zhou J, et al. Bacterioplankton community resilience to ocean acidification: evidence from microbial network analysis. ICES J Mar Sci. 2016;73: 865–875.

55. Duran-Pinedo AE, Paster B, Teles R, Frias-Lopez J. Correlation network analysis applied to complex biofilm communities. PLoS One. 2011;6: e28438. pmid:22163302