Association of exacerbation phenotype with the sputum microbiome in chronic obstructive pulmonary disease patients during the clinically stable state

Abstract

Background

Chronic obstructive pulmonary disease (COPD) is a progressive, life-threatening lung disease with increasing prevalence and incidence worldwide. Increasing evidence suggests that lung microbiomes might play a physiological role in acute exacerbations of COPD. The objective of this study was to characterize the association of the microbiota and exacerbation risk or airflow limitation in stable COPD patients.

Methods

The sputum microbiota from 78 COPD outpatients during periods of clinical stability was investigated using 16S rRNA V3-V4 amplicon sequencing. The microbiome profiles were compared between patients with different risks of exacerbation, i.e., the low risk exacerbator (LRE) or high risk exacerbator (HRE) groups, and with different airflow limitation severity, i.e., mild to moderate (FEV1 ≥ 50; PFT I) or severe to very severe (FEV1 < 50; PFT II).

Results

The bacterial diversity (Chao1 and observed OTUs) was significantly decreased in the HRE group compared to that in the LRE group. The top 3 dominant phyla in sputum were Firmicutes, Actinobacteria, and Proteobacteria, which were similar in the HRE and LRE groups. At the genus level, compared to that in the LRE group (41.24%), the proportion of Streptococcus was slightly decreased in the HRE group (28.68%) (p = 0.007). However, the bacterial diversity and the proportion of dominant bacteria at the phylum and genus levels were similar between the PFT I and PFT II groups. Furthermore, the relative abundances of Gemella morbillorum, Prevotella histicola, and Streptococcus gordonii were decreased in the HRE group compared to those in the LRE group according to linear discriminant analysis effect size (LEfSe). Microbiome network analysis suggested altered bacterial cooperative regulation in different exacerbation phenotypes. The proportions of Proteobacteria and Neisseria were negatively correlated with the FEV1/FVC value. According to functional prediction of sputum bacterial communities through Phylogenetic Investigation of Communities by Reconstruction of Unobserved States (PICRUSt) analysis, genes involved in lipopolysaccharide biosynthesis and energy metabolism were enriched in the HRE group.

Conclusion

The present study revealed that the sputum microbiome changed in COPD patients with different risks of exacerbation. Additionally, the bacterial cooperative networks were altered in the HRE patients and may contribute to disease exacerbation. Our results provide evidence that sputum microbiome community dysbiosis is associated with different COPD phenotypes, and we hope that by understanding the lung microbiome, a potentially modifiable clinical factor, further targets for improved COPD therapies during the clinically stable state may be elucidated.

Background

Chronic obstructive pulmonary disease (COPD) is a common respiratory disease characterized by airflow limitation, lung function impairment, and airway inflammation. The incidence of COPD has increased worldwide, and it will become the third most prevalent cause of death by 2030 [1]. The Epidemiology and Impact of COPD (EPIC) Asia survey concluded that there was a high prevalence of COPD and a substantial socioeconomic burden of the disease in nine Asia–Pacific regions [2]. Airway obstruction, which is confirmed by spirometry, leads to air trapping and shortness of breath in response to physical exertion, and the poorly irreversible airway obstruction that characterizes COPD is progressive. In 2007, the grading of COPD severity proposed by the Global Initiative for Chronic Obstructive Lung Disease (GOLD) was based on forced expiratory volume in 1 s (FEV1) only. This classification could not adequately predict clinical outcomes [3, 4]. Exacerbations of COPD are important events in the disease course, and in particular, mortality in the year after an exacerbation requiring hospital admission is estimated to be as high as 21% [5]. The Evaluation of COPD Longitudinally to Identify Predictive Surrogate End-points (ECLIPSE) study suggested that individuals with two or more exacerbations in a given year represent a distinct frequent exacerbation phenotype [6]. However, exacerbations can occur across all stages of airflow limitation measured by FEV1, which emphasizes the need to identify other predictors of high exacerbation risk [6]. The GOLD 2011 revision presented an ABCD classification that combined respiratory symptoms, risks of exacerbations, and airflow limitations as indicated by FEV1 [7]. In 2017, The Global Initiative for Chronic Obstructive Lung Disease (GOLD) guide to COPD diagnosis and management used a threshold of two or more acute exacerbations in the previous year or at least one hospital admission related to an acute exacerbation to identify individuals at high risk of future events (groups C and D) and separated spirometric grades from the “ABCD” groupings [8].

Recent data have suggested that COPD is a complex and heterogeneous disease resulting from a number of different pathological processes, including infections [9]. Bacterial pathogens are commonly identified in the respiratory tracts of patients both in the stable state and during acute exacerbations, with significant changes in the prevalence of airway bacteria occurring during acute exacerbations of COPD [10,11,12]. An increased bacterial load has been associated with a decline in lung function [13] and increased rates of exacerbations in COPD patients [14], suggesting an important role of bacteria in the pathogenesis of COPD. Advances in next-generation sequencing platforms for 16S rRNA gene sequencing have provided opportunities to study the lung microbiomes in COPD patients, and the results have suggested that changes in the lung microbiota may be associated with enhanced airway inflammation and disease progression [15]. In patients with COPD, two large studies (AERIS study and COPDMAP study) have recently investigated the value of respiratory microbiome research to understand the association of microbiome changes and COPD exacerbations [16, 17]. Interestingly, both studies also identified alteration in the taxonomic composition of the lung microbiome related to the frequency of exacerbation. Pragman et al. reported that even during periods of clinical stability, the frequent exacerbation phenotype is associated with decreased alpha diversity, beta diversity clustering, and changes in taxonomic abundance [18].

Most patients with COPD are stable outpatients, and an important clinical challenge is to provide COPD outpatients appropriate education and prescribe appropriate therapy to prevent exacerbations. The current parameters for evaluating the severity of COPD patients are airflow limitation, symptoms, and exacerbation risks. Therefore, we undertook the present study to determine if the lung microbiome is associated with COPD clinical assessment parameters. Our hypothesis is that the lung microbiota alteration is important factor of exacerbation risk or airflow limitation in COPD outpatients during periods of clinical stability. We hope that by understanding the lung microbiome, a potentially modifiable clinical factor, further targets for improved COPD therapies may be elucidated.

Methods

Study subjects and study design

Seventy-eight subjects with COPD were enrolled in this study. The study was approved by the Institutional Review Board of the Tri-Service General Hospital Taiwan, and all subjects were enrolled from April 2015 to April 2016 and provided written informed consent. The patients were included if they were aged ≧ 40 years and had been diagnosed with COPD. We categorized each patient in our cohort according to GOLD 2017 classifications. Patients had no exacerbations or infections for at least 30 days prior to sample collection, and their sociodemographic and clinical data were recorded. For full inclusion and exclusion criteria see Additional file 1: Table S1. Symptoms were quantified with the modified Medical Research Council (mMRC) scale. The severity of exacerbations was classified as mild in the case of self-management with short-acting bronchodilators only; moderate if the patient was not hospitalized but received a prescription of systemic corticosteroids, antibiotics or both; and severe if the patient was hospitalized. The 78 COPD patients were classified by exacerbation risk as low-risk exacerbators (LREs, < 2 moderate exacerbations and no severe exacerbations per year, n = 60) and high-risk exacerbators (HREs, ≥ 2 moderate or severe exacerbations or ≥ 1 hospitalizations for COPD exacerbation, n = 18). Furthermore, according to pulmonary function testing, the 78 COPD patients were grouped as mild-to-moderate airflow limitation with a value of forced expiratory volume in the post-bronchodilator (post-BD) second FEV1% predicted ≥ 50 (PFT I, n = 43) and as severe-to-very severe airflow limitation with a value of post-BD FEV1% predicted < 50 (PFT II, n = 35).

Sputum sampling and processing

Sputum was collected from patients by induction during stable visits after pulmonary function testing. All participants reported that they had no special dietary habits, had no known periodontal disease, had not taken systemic antibiotics in 4 weeks, and had not used antiseptic mouthwash before sample collection. Sputum was induced according to a previous protocol, with slight modifications [19]. For detailed methods for sputum sampling in this study, see Additional file 1.

DNA extraction

Total bacterial genomic DNA was isolated from the sputum samples using the QIAamp DNA Microbiome Kit (Qiagen, USA). Briefly, 250 µl of AHL buffer was added to 500 µl of sample for host cell lysis, followed by digestion of the host nucleic acids with 1.25 μl of benzonase and 10 μl of proteinase K. The host DNA was separated by centrifugation, after which 100 μl of ATL buffer was added to the bacterial cells in a pathogen lysis tube L and the sample was vortexed using a TissueLyser LT for 10 min at 30 Hz. The bacterial DNA was washed, eluted using nuclease-free water, and stored at − 80 °C. The concentrations and qualities of the purified DNA were determined with a Qubit high-sensitivity dsDNA assay (Life Technologies).

Sputum microbiota profiling by 16S rRNA gene sequencing

A 16S rRNA gene amplicon library targeting the 16S rRNA V3-V4 region was constructed as in a previous report [20]. Illumina adaptor overhang nucleotide sequences were added to the gene-specific forward and reverse primers. Two-round PCRs were performed, and the final amplicon libraries were approximately 630 bp in length. The multiplex amplified libraries were pooled equally and sequenced on a MiSeq system with 2 × 300 paired-end v3 sequencing reagents (Illumina, USA). For detailed methods for 16S rRNA gene sequencing, see Additional file 1.

Bioinformatics analysis

The sequencing reads were processed, and the taxonomic classification was performed using FLASH (version 1.2.11) [21]. Low-quality reads were filtered [22], and only sequence tags with lengths > 400 bp were retained for subsequent analysis. The operational taxonomic units (OTUs) were clustered at 97% sequence similarity using USEARCH (version 9.2.64) [23] against the Greengenes 16S rRNA gene database (13_8 release), and final taxonomic assignments were performed using the RDP classifier [24]. Furthermore, our results were validated by another bioinformatic pipeline using the DADA2 package [25] for modelling and amplicon error correction, which was followed by quality filtering, dereplication, denoising, merging and chimaera removal. A naïve Bayes classifier [24] was trained using the most recent available version of the Silva (version 132) sequences for taxonomic assignments. A bivariate correlation analysis of the 15 most abundant genera using Spearman’s correlation coefficient was performed in R. We then constructed a co-occurrence network of the predominant sputum microbiota with different COPD disease severities. The network was generated using Cytoscape (version 3.7.0) and visualized using a circular layout [26]. Potential biomarkers were determined using linear discriminate analysis effect size (LEfSe) [27]. Microbial functionality profiles were predicted using Phylogenetic Investigation of Communities by Reconstruction of Unobserved States (PICRUSt) to generate the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways [28]. For detailed methods for bioinformatics analysis, see Additional file 1.

Statistical analyses

Box and whiskers plots (10–90 percentile) of alpha diversity indices and taxonomic abundances comparing two groups were plotted using GraphPad Prism 6 (GraphPad Software, Inc., La Jolla, CA, USA). Principal component analysis plots for LRE vs HRE, PFT I vs PFT II, and smoker vs non-smoker were prepared using unweighted UniFrac distance. The Mann–Whitney Wilcoxon test followed by Bonferroni correction was used to test for significant differences in alpha diversity or taxonomic levels between groups. Bonferroni-adjusted p values were calculated as an alpha error of 0.05 divided by the numbers of parameters in each table. Spearman's correlation was used to determine the association of the dominant phyla or genera with FEV1/FVC, and p values were adjusted by using Bonferroni correction for multiple tests. Statistical analyses were performed by using R or SPSS (SPSS Inc., Chicago, IL, USA).

Results

COPD patient demographics and sputum microbiota profiling

To characterize the lung microbiome constituents that differentiate the low risk vs high risk exacerbation groups or mild-to-moderate vs severe-to-very severe airway limitation groups under stable COPD conditions, sputum specimens were collected from 78 COPD patients during stable visits. The characteristics of the patients are summarized in Table 1, and the detailed information regarding inclusion and exclusion criteria are available in Additional file 1: Table S1. The patients’ age range was 40–93 years, and none of the patients had been on antimicrobial therapy for the last 4 weeks prior to sample collection. Ages and smoking histories were similar in COPD patients at different stages. Among these 78 patients, 31 patients (39.7%) with COPD received a long-acting bronchodilator (LAB) and 45 patients (57.6%) received both long-acting β2 agonists (LABA) and inhaled corticosteroid (ICS) treatment. The 78 COPD patients were divided into two subgroups, including LRE (n = 60) and HRE (n = 18). Furthermore, the 78 COPD patients were grouped as mild-to-moderate airflow limitation with a value of FEV1% predicted ≥ 50 (PFT I, n = 43) and as severe-to-very severe airflow limitation with a value of FEV1% predicted < 50 (PFT II, n = 35).

Table 1 Clinical characteristics of the study population (N = 78) in this study

DNA isolated from the sputum specimens was subjected to two-round PCR amplification and library construction of the 16S V3-V4 region and was sequenced on an Illumina MiSeq system. The average number of raw reads was 228,741 for all COPD patients (Additional file 1: Table S2). After selecting the qualified reads, the average number of paired quality-filtered reads was 133,386 for all COPD patients (Additional file 1: Table S2). The rarefaction curve showed that the sequencing depth per sample was enough to represent most of the community diversity and reached a saturated plateau phase (data not shown).

We first compared the bacterial compositions and diversity of the sputum microbiomes in different COPD phenotypic subgroups. Regardless of exacerbation phenotype, the unweighted principal component analysis (PCA) plot revealed a similar sputum microbiota between the LRE and HRE groups (Fig. 1a). Species richness according to Chao1 index (p = 0.002) and observed OTUs (p = 0.002) was lower in the HRE group, indicating that the sputum microbiota in these patients was characterized by a lower diversity than that in the LRE patients (Fig. 1b). The species evenness according to Shannon index was similar in the LRE and HRE groups (p = 0.652) (Fig. 1b). We then analysed the microbiome diversity in patients with different airway limitation severity. The unweighted PCA plot (Fig. 1c) and alpha diversity according to Chao1, observed OTUs, and Shannon indices (Fig. 1d) were similar between the PFT I vs PFT II groups. The bacterial compositions and diversity of the sputum microbiomes in smokers and non-smokers were similar (Additional file 1: Figure S1). These results indicated that the sputum microbiome was altered and that these compositional changes may have been associated with exacerbation risk.

Fig. 1

Analysis of the sputum microbiome and diversity in stable COPD. The unweighted PCA plot (a) and Chao1 index, observed OTUs, and Shannon index (b) were calculated in the LRE vs HRE groups. The unweighted PCA plot (c) and the diversity index (d) were presented in the PFT I vs PFT II groups. The box and whiskers plots show the median, 10th and 90th percentile in each group. Bonferroni-adjusted p-values < 0.05/6 = 0.0083 indicate significance. LRE low-risk exacerbator, HRE high-risk exacerbator

Most abundant bacterial taxa changes in stable COPD patients with high risk of exacerbation

Figure 2 shows the dominant taxa at the phylum and genus levels. The sputum microbiome dataset from our cohort revealed a total of five phyla, which accounted for 99% of all the bacteria (Fig. 2a; Table 2). The top five dominant phyla were Firmicutes, Actinobacteria, Proteobacteria, Bacteroidetes, and Fusobacteria. Firmicutes was the most dominant phylum in stable COPD patients, with a relative abundance of 43.5–54.6%, with Actinobacteria being the second most dominant phylum, with a relative abundance of 16.6–16.7% (Table 2). The proportions of the dominant phyla were similar in the LRE vs HRE or PFT I vs PFT II groups (Fig. 2b, c; Table 2). The relative abundance of Proteobacteria was slightly increased but not significantly different in the HRE group compared with that in the LRE group (p = 0.031) (Fig. 2b).

Fig. 2

Relative abundance of the most prevalent bacteria at the phylum and genus levels in stable COPD patients. The bacterial taxonomic distributions are shown at the phylum and genus levels in the sputum microbiomes of COPD patients. The stacked bar represents differentially expressed bacteria at the phylum level (a) in different COPD phenotypic subgroups. The box and whisker plots (10–90 percentile) show the relative abundance of Proteobacteria in the LRE vs HRE groups (b) or in the PFT I and PFT II groups (c). The stacked bar represents differentially expressed bacteria at the genus level (d) in the different COPD phenotypic subgroups. The box and whisker plots show the relative abundance of Streptococcus in the LRE vs HRE groups (e) or in the PFT I and PFT II groups (f). The p value was analyzed using the Mann–Whitney Wilcoxon test. Bonferroni-adjusted p-values < 0.05/15 (5 phyla and 10 genera) = 0.0033 indicate significance. LRE low-risk exacerbator, HRE high-risk exacerbator

Table 2 Taxonomic identification at the phylum and genus levels in COPD

The top 10 dominant genera in sputum were Streptococcus, Rothia, Haemophilus, Neisseria, Veillonella, Granulicatella, Porphyromonas, Leptotrichia, Actinomyces, and Capnocytophaga (Fig. 2d; Table 2). Among the top 10 genera, the relative abundance of Streptococcus was slightly decreased but not significantly different in the HRE group (28.68%) compared to that in the LRE group (41.24%) (p = 0.007) (Fig. 2e; Table 2). The relative abundance of Streptococcus was similar between the PFT I (39.94%) and PFT II groups (34.83%) (p = 0.226) (Fig. 2f). We also classified the participants into older-aged adults (60 years and above, n = 67). Among the older-aged COPD patients, the proportions of the dominant phyla and genera were similar in the LRE (n = 51) vs HRE groups (n = 16) (Additional file 1: Table S3). A trend of slightly increased Proteobacteria and decreased Streptococcus abundances was also observed in older-aged COPD patients with different risks of exacerbation. However, the taxonomic distribution was similar in smoker (n = 61) vs non-smoker (n = 17) COPD patients (Additional file 1: Table S4).

Differential taxa in the microbiome of stable COPD patients

A LEfSe analysis was performed to identify the differences in taxonomic distributions associated with different COPD phenotypic subgroups. The cladogram plotted based on LEfSe analysis shows five taxonomic levels, with the phyla levels and genera levels plotted in the innermost ring and outermost ring, respectively (Additional file 1: Figure S2). The amount of Pseudomonadales, which is an order of the phylum Proteobacteria, was enriched and had the highest linear discriminant analysis (LDA) score in HRE subjects. Otherwise, the amounts of Bacilli and Lactobacillales, which is a subclass of phylum Firmicutes, were enriched with in LRE subjects (Additional file 1: Figure S2).

To identify the differential bacterial taxa at the species level, we performed LEfSe analysis based on the known OTUs at the species level. Many bacteria at the species level were enriched in different subgroups of COPD; therefore, we selected a LDA score higher than 3 or lower than -3 to represent the most significantly enriched species in each group (Fig. 3). The relative abundances of Gemella morbillorum (G. morbillorum), Prevotella histicola (P. histicola), and Streptococcus gordonii (S. gordonii) were significantly decreased in HRE compared to LRE subjects (Fig. 3a). Furthermore, the relative abundances of Gemella morbillorum (G. morbillorum), Veillonella atypica (V. atypica), and Corynebacterium durum (C. durum) were significantly decreased in patients with severe-to-very severe airflow limitation compared with patients with mild-to-moderate airflow limitation (Fig. 3b). In our cohort, the percentages of patients who were treated with a long-acting bronchodilator (LAB) or long-acting β2 agonists (LABA) and an inhaled corticosteroid (ICS) were approximately 39.7% and 57.7%, respectively. We also compared the microbiomes of COPD patients who received different medications. The Chao1 index and observed OTUs were similar between these two medications (Additional file 1: Figure S3). The proportions of Fretibacterium fastidiosum and Oribacterium sinus were enriched in LAB and LABA plus ICS, respectively (Fig. 3c).

Fig. 3

Linear discriminant analysis (LDA) effect size (LEfSe) revealed differentially abundant bacterial taxa in stable COPD. The plots from LEfSe analysis indicated differentially enriched bacteria in the LRE vs HRE groups (a), PFT I vs PFT II groups (b), and different medical treatments (c). LRE low-risk exacerbator, HRE high-risk exacerbator

Correlation between lung function and the microbiome

To explore the potential bacterial co-existence and co-exclusion relationships, we performed an interaction network analysis. We first selected the top 15 most abundant genera in each group, and the specific network was built and estimated based on the relative abundances of bacterial genera using SparCC correlation coefficients. Each node represents a genus of bacteria, and the red and blue lines represent positive and negative correlations, respectively. All plotted nodes of the networks with significant coefficients are shown in Fig. 4. In total, 12 and 9 nodes were constructed in the LRE (Fig. 4a) and HRE groups (Fig. 4b), respectively. In HRE subjects, Moraxella was included in a closed negatively correlated network containing Streptococcus, Haemophilus, Moraxella, Capnocytophaga, Lactobacillus, and Porphyromonas. (Fig. 4b). Furthermore, Actinomyces showed a negative correlation with Moraxella in HRE subjects (Fig. 4b). The heat map of Spearman correlation coefficients between the top 15 genera is shown in the LRE and HRE groups (Fig. 4c, d).

Fig. 4

Sputum microbiome networks in COPD patients in different phenotypic subgroups. The networks of the top 15 genera were built using Spearman correlation coefficients in COPD patients. In total, 12 and 9 nodes were constructed in the LRE (a) and HRE groups (b), respectively. The nodes represent bacterial genera. The red and blue lines represent positive and negative correlations, respectively. Heatmap of the Spearman correlation matrix from the LRE (c) or HRE groups (d) with the top 15 genera. In the figure, the larger the circle, the higher is the correlation coefficient. Red dots represent positive correlations, and blue dots represent negative correlations. The correlation values ranged from − 1.00 (blue) to 1.00 (red). LRE low-risk exacerbator, HRE high-risk exacerbator

To explore the potential relationship between different bacterial taxa and lung function in COPD patients, we performed Spearman’s correlation analyses using the 16S rRNA gene sequence dataset. Among the most abundant phyla and genera, the relative abundances of Proteobacteria and Neisseria were negatively correlated with FEV1/FVC (Fig. 5a, b). In a comparison of the Chao1 and Shannon indices with the FEV1/FVC value, there was no significant association with bacterial diversity and lung function (Fig. 5c).

Fig. 5

Correlation between the sputum microbiome and lung function in stable COPD. a Spearman correlations between the relative abundances of bacteria and the values of post-BD FEV1/FVC were calculated and plotted. b Proteobacteria and Neisseria show negative correlations with the values of post-BD FEV1/FVC. c The association of the Chao1 and Shannon indices and the values of post-BD FEV1/FVC. Bonferroni-adjusted p-values < 0.05/27 (5 phyla, 20 genera, 2 diversity indices) = 0.0019 indicate significance

Functional analysis of the microbiome in stable COPD patients by PICRUSt analysis

To explore the predicted functional capacity of the microbiome involved in COPD, we performed a PICRUSt analysis to predict the sputum microbiome functions in stable COPD patients. The results predicted a number of KEGG pathways that were slightly enriched or depleted associated with COPD exacerbation risk. The levels of metabolism, such as glycan biosynthesis and metabolism, lipopolysaccharide biosynthesis, sulphur metabolism, and biotin metabolism, were positively associated with high exacerbation risk of COPD, as their abundances were higher in the HRE group than in the LRE group (Fig. 6a, b). In contrast, parameters related to the phosphotransferase system, fructose and mannose metabolism and galactose metabolism were increased in the LRE group (Fig. 6c, d). These results may suggest that an altered sputum microbiome affect the nutrient availability, sugar metabolism, or growth conditions.

Fig. 6

Differentially enriched functions in COPD patients by PICRUSt analysis. Comparison of the relative abundance of the PICRUSt-generated functional profile of the sputum microbiome in COPD patients. Upregulated KEGG pathways associated with HRE subjects are plotted as a heatmap (a) and bar graph (b). Bonferroni-adjusted p-values < 0.05/7 = 0.0071 indicate significance. Upregulated KEGG pathways associated with LRE subjects are plotted as a heatmap (c) and bar graph (d). Bonferroni-adjusted p-values < 0.05/11 = 0.0045 indicate significance. GOLD A, B, C, and D were identified via the GOLD 2017 classification. LRE low-risk exacerbator, HRE high-risk exacerbator

Discussion

In 2017, the GOLD announced another major revision of the COPD guidelines. The recommendations for medications used in each patient group were also updated [8]. In this study, we extensively investigated the sputum microbiota in stable COPD patients with a high-throughput 16S rRNA gene sequencing analysis. Overall, the microbiome diversity decreased in HRE subjects, and the composition of the sputum microbiome changed between the two different exacerbation risk subgroups. We further explored the significant taxa and the predicted functional analyses associated with different COPD phenotypic subgroups. Collectively, our findings suggest notable airway microbiome changes in stable COPD patients and highlight its cooperative network and functional capacity.

Our results demonstrated that the bacterial richness was reduced in HRE COPD patients with even during periods of clinical stability. Mayhew et al. reported that the bacterial diversity index was reduced in very severe COPD patients compared with moderate COPD patients in sputum samples obtained in both stable and exacerbated states [16]. Diao et al. also reported that the OTU richness in throat swab samples from COPD patients was lower than that observed in samples obtained from healthy controls [29]. These results are consistent with our findings in this study. In contrast, Pragman et al. reported an increase in the microbial diversity index with the development of COPD in bronchoalveolar lavage fluid (BALF) samples in a cohort that included 14 moderate and 8 severe COPD patients [30]. These discrepancies may be due to the different geographical areas and the different sampling methods used.

In our study, the most dominant phylum in sputum samples was Firmicutes (approximately 43.5–54.6%), which was consistent with other reports that analysed sputum samples and identified the proportion of Firmicutes as approximately 40–50% [16, 31]. Our finding contrasts with the results of Garcia-Nunez et al. who also analysed sputum samples during stable COPD and observed that Proteobacteria (44%) and Firmicutes (16%) were the first and second most abundant phyla, respectively, in 17 moderate to advanced COPD patients [32]. The most dominant phylum in bronchial wash samples from stable COPD patients was Firmicutes [33]. Otherwise, Bacteroides (approximately 40–60%) was the most abundant phylum in lung and throat swab samples [29, 34]. There are differences in the microbiome at different locations within the respiratory tract [35]. The microbiota of the induced sputum might represent the upper airway microbiota, which would explain the difference in the microbiome compositions from the BALF or lung samples. Pragman et al. profiled the microbiome in oral, bronchial, and lung tissue samples from individual patients and observed that oral bacteria are true members of the early-stage COPD lung microbiota and exhibit ecological drift [36].

In the present study, the top 7 dominant genera present in COPD patient sputum samples were Streptococcus, Rothia, Haemophilus, Neisseria, Veillonella, Granulicatella, and Porphyromonas, genera to which some pathogenic bacteria belong. Streptococcus was the most common genus in the oral, bronchial, and lung tissue samples in COPD patients [36]. Of the 5 dominant genera identified from sputum samples by Mayhew et al. our study identified 3 common taxa (Streptococcus, Haemophilus, and Veillonella) among our 7 most common taxa. In the report by Tangedal et al. the 7 most dominant genera in induced sputum samples were identical to those observed in our results [37]. Furthermore, the relative abundance of Streptococcus was decreased in the stable HRE COPD patients compared with that in the LRE patients in our cohort. In contrast, the relative abundance of Streptococcus was not altered between the groups with different lung function levels. Streptococcus pneumoniae infection is reported to be important pathological bacteria in COPD in association with disease exacerbation or airway limitation [38]. However, the 16S rRNA sequencing platform was unable to discriminate between the various species in the genus of Streptococcus because the sequences of 16S rRNA gene in this genus are relatively similar.

Some reports have indicated that bacterial diversity was positively associated with lung function [30, 32], whereas other researchers observed that the bacterial diversity was not correlated with lung function [29]. In our study, the Chao1 diversity index and the observed OTUs were not associated with lung function (FEV1/FVC) in the 78 stable COPD patients. Another interesting finding from our study was that the proportions of Proteobacteria and Neisseria were negatively correlated with FEV1/FVC in stable COPD patients. The proportion of Neisseria was slightly increased in the HRE group (3.42%) compared to that in the LRE group (3.07%) in our cohort. Neisseria are gram-negative bacteria that belong to the family Neisseriaceae and are present in mucosal surfaces in the upper respiratory and genitourinary tracts [39]. Commensal Neisseria can be opportunistic pathogens in humans, and some clinical cases of infections with Neisseria species such as N. bacilliformis in sputum samples from patients with bronchitis and with N. flavescens, N. lactamase, and N. mucosa in lung samples from patients with pneumonia have been reported [40, 41].

The proportions of Gemella morbillorum (G. morbillorum, LDA score in LRE = 3.64, p value = 0.033) and Prevotella histicola (P. histicola, LDA score in LRE = 3.54, p value = 0.016) were significantly decreased in HRE subjects in the present study. G. morbillorum and P. histicola are some of the normal flora of the mucous membranes, predominantly of the oropharynx, but can also be found in the upper respiratory and other sites. P. histicola may suppress the production of inflammatory cytokines, and P. histicola suppresses disease in the animal model of multiple sclerosis or arthritis [42, 43]. These results may suggest that the altered normal flora distribution may affect the colonization of pathogenic bacteria and enhanced inflammation, leading to airway obstruction and increased exacerbation risk.

Functional prediction showed that lipopolysaccharides (LPS) biosynthesis, which produces the main cell wall components of gram-negative bacteria, was slightly enriched in HRE subjects. LPS challenging may lead to air flow limitation, decreases in the level of FEV1 and enhanced pulmonary inflammation, suggesting that this occurs during exacerbations in COPD patients [44, 45]. Also, sulphur and biotin metabolism were enriched in HRE subjects [46, 47]. These pathways are essential for the survival of bacteria, including some pathogens. Furthermore, the levels of phosphotransferase system-related parameters, such as fructose and mannose metabolism and galactose metabolism were significantly reduced in COPD patients compared to the controls. The phosphotransferase system plays a pivotal role in the uptake of multiple sugars in bacteria. Glucose concentration in the airway might contribute to bacterial infections, and impaired glucose metabolism was also observed in COPD [48, 49]. These results may suggest that the normal flora balance is disrupted and that the energy and metabolic machinery of normal and/or pathogenic bacteria are altered.

An important limitation of our study was that no follow-up was conducted. Further studies with longitudinal sampling from each individual at both stable and exacerbation time points will be important for monitoring the microbiome dynamics, clinical phenotypes and treatment responses. Another limitation of our study was that the patients were enrolled at a single site; therefore, exploring other variables between different areas was not possible. A third limitation was that fewer female subjects were enrolled. We observed that among the PFT II group (n = 35), only about 40% of patients (n = 14) had severe airflow limitation (FEV1 < 50) and frequent exacerbations, although this finding could be due to the small sample size of this subgroup in our cohort during the one-year study. Interestingly, it was previously reported that “high-risk” COPD patients (GOLD groups C and D) are highly heterogeneous populations [50]. Moreover, environmental and occupational exposures may affect the exacerbation of COPD patients.

Conclusions

In conclusion, the present study revealed that the sputum microbiome changed in different exacerbation risk subgroups of COPD. Additionally, the bacterial cooperative networks were different in the different COPD phenotypic subgroups. An altered lung microbiome can have an important effect on the host immunity and initiates disease pathogenesis, promotes chronic inflammation, or merely serves as a marker of injury and inflammation. Understanding the mechanisms driving bacterial compositions and diseases will help us prevent or treat COPD.