Sputum Exosomal microRNAs Profiling Reveals Critical Pathways Modulated By Pseudomonas aeruginosa Colonization In Bronchiectasis

Introduction

Chronic airway infection and inflammation are implicated in the pathogenesis and progression of bronchiectasis.¹ The dilated bronchi provide the niche for the expansion of opportunistic pathogens that elicit airway inflammation.² Compared with other potentially pathogenic bacteria, Pseudomonas aeruginosa (PA) conferred significantly greater adverse impacts on bronchiectasis,³ including the more rapid lung function decline and frequent exacerbations.⁴ In fact, colonization of PA correlated with a substantial morbidity and mortality.⁵

Micro-ribonucleic acids (miRNA), the 18–25nt non-coding ribonucleic acids, have been implicated in the airway inflammatory responses.⁶ Up-regulated expression of miR-301b correlated with the heightened inflammation associated with PA infection.⁷ Exosomes are extracellular vesicles laden mainly with non-coding RNAs such as miRNAs,⁸ which are transferred to the target cells for modulating the host immunity.⁹ For instance, exosomal miR-155 enhanced whereas miR-146a suppressed inflammatory gene expression,¹⁰ suggesting that exosomal miRNA might be capable of modulating airway inflammation.

Currently, the mechanisms underlying the adverse impacts associated with PA colonization in bronchiectasis are unclear. Here, we sought to identify the differentially expressed miRNAs (DEMs) in sputum exosomes that were associated with PA colonization. This might help clarify whether and which miRNAs could be the biomarkers associated with PA colonization, and through which pathways PA modulates chronic airway inflammation in bronchiectasis.

Methods

Participants

Between November 2016 and March 2018, we consecutively recruited 98 bronchiectasis patients and 17 healthy subjects. Bronchiectasis was diagnosed with chest high-resolution computed tomography in the out-patient clinics of The First Affiliated Hospital of Guangzhou Medical University. Eligible patients remained clinically stable (respiratory symptoms not significantly exceeding normal daily variations) for at least 4 weeks. Healthy subjects had normal chest X-ray findings, and had no history of chronic respiratory diseases. None of the participants had malignancy, acute upper airway infection or antibiotic use within 4 weeks before recruitment. Exacerbations were defined based on the criteria stated in international guidelines.¹¹ The study protocol was approved by The Ethics Committee of The First Affiliated Hospital of Guangzhou Medical University (Medical Ethics Year 2012, the 33rd). Written informed consent was obtained.

Study Design

Patients underwent history inquiry, spirometry (COSMED Inc., Italy) and collection of spontaneous sputum when clinically stable. Patients were followed-up at 3-to-6-month intervals to determine the bacterial colonization status, and record exacerbations and hospitalizations. We excluded patients followed-up for less than 6 months. Modified Reiff score was used to evaluate the radiologic severity of bronchiectasis.^12,13 Etiologies, Bronchiectasis Severity Index (BSI), and E-FACED score were assessed^5,14,15 at the initial stable visit. We validated the DEMs with quantitative polymerase chain reaction (qPCR) in 34 bronchiectasis patients who had sufficient paired sputum.

Procedures

Patients thoroughly gargled their mouth with distilled water, followed by forceful coughing. Sputum was expectorated in a sterile plastic pot. We selected the most purulent portion and performed cytology assessment under the light microscope. Sputum that met quality control criteria was split for bacterial culture (within 2 hrs), sputum sol preparation for inflammatory biomarker [interleukin (IL)-8, IL-1β, and tumor necrosis factor-α (TNF-α)] assays, and storage in −80°C freezers for subsequent extraction of the exosomes and miRNA detection. Sputum exosomes were identified with the H-7650 transmission electron microscope (HITACHI, Japan), and the DEMs associated with PA colonization were validated with qPCR by using ABI 7500 Real-Time System (Thermo Scientific, USA).

Definition Of PA Colonization

PA colonization was defined as sputum culture positive of PA for at least two occasions (at least 3 months apart) within 1 year, when clinically stable during the longitudinal follow-up. Bronchiectasis patients were stratified into PA and non-PA group (unmet the criteria of PA colonization).

Sputum miRNA Expression Profiling

After thawing, sputum was mixed with phosphate buffer solution. Exosomes were isolated with exoEasy Maxi Kit, and the total RNA was extracted with exoRNeasy Serum/Plasma Maxi Kit (Qiagen, Germany). Agilent Bioanalyzer 2100 system (Agilent Technologies, USA) was used to assess the purity and integrity of RNA, followed by sequencing with Hiseq 2500 platform (Illumina, USA). After denoising the raw reads, we calculated the false discovery rate (FDR), P value and fold-change for high-quality reads. MiRNA expression was corrected by transcripts per million.^16,17 Mature human miRNAs were annotated with miRBase 22.0. See Figure S1 regarding the flowchart of sequencing schemes. See further details in Online Supplement.

Statistical Analysis

Statistical analyses were conducted with R software version 3.5.0 and Graphpad Prism 5.0 (Graphpad Inc., USA). Kolmogorov–Smirnov test was applied to determine the normality of continuous variables, which were presented as mean ± standard deviation or median (interquartile range) as appropriate. Count (percentage) was demonstrated for categorical data.

Principal component analysis (PCA) was applied for quality control and assessment of the global gene expression patterns. Identification of the DEMs was based on the FDR (cut-off: 0.10).^18–20 The DEMs were initially screened via comparison between all bronchiectasis patients and healthy subjects, followed by the comparison between PA and non-PA group in bronchiectasis patients. Venn diagram was applied to demonstrate the unique or overlap DEMs. Multiple-group comparisons were made using analysis-of-variance or non-parametric analysis and adjusted with Benjamini–Hochberg correction. The diagnostic value of miRNA in differentiating PA group from non-PA group was analyzed with the C-statistic and 95% confidence interval (95% CI). To select the minimum set of DEMs that best discriminated PA colonization, machine learning was applied to perform least absolute shrinkage and selection operator linear regression analysis (LASSO)²¹ to avoid overfitting. Multiple linear regression was performed to analyze the correlation between DEMs and PA colonization. We finally included the miRNA for validation by performing qPCR assays, provided that the C-statistic was greater than 0.6, and that the expression level correlated significantly with PA colonization in the above-mentioned multiple regression models. The biological function and pathways of DEMs were analyzed with Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes Pathway Analysis (KEGG). Significant enrichment was defined as the unadjusted P value <0.01 and P_adj <0.05 (corrected with the Benjamini–Hochberg algorithm).

Results

Subject Recruitment And Clinical Characteristics

Of the 115 participants who underwent screening, 82 bronchiectasis patients and 9 healthy controls were included in the final analysis. Twenty-four participants were excluded because of incomplete bacterial culture data or low-quality sequencing results (Figure 1).

Figure 1 Study participants recruitment flow chart. We screened for 115 study participants, 22 of whom yielded sputa with a poor quality and were therefore excluded for sequencing. Two patients whose data for sequencing with low quality were also excluded from further analysis. Hence, the final cohort consisted of 91 study participants (82 with stable bronchiectasis and 9 healthy subjects).

Table 1 shows the demographic and clinical characteristics of the study participants. Most patients had post-infectious and idiopathic bronchiectasis. Thirty-four bronchiectasis patients (41.5%) had PA colonization, who had worse lung function and greater disease severity than those without. The demographic and clinical characteristics of PA and non-PA group are shown in Table S1, and those in RT-qPCR cohort are shown in Table S2.

Table 1 Clinical Characteristics Of Patients With Clinically Stable Bronchiectasis And Healthy Subjects

Main Sequencing Findings

We first applied the transmission electron microscope to confirm the successful extraction of extracellular vesicles (50–200 nm, Figure S2). On average, 10 million raw reads per sample were sequenced. We identified 1254 miRNAs in PA group, 1403 miRNAs in non-PA group, and 1132 miRNAs in healthy subjects. Next, we applied PCA and dendrogram to demonstrate the global miRNA expression patterns (Figure S3). There was some separation between PA group and healthy subjects, but less so between PA and non-PA group. The flowchart of the analysis is shown in Figure S4.

DEMs Between PA And Non-PA Group

After quality control of sequencing data, we identified 190 (129 up-regulated and 61 down-regulated) DEMs between all bronchiectasis patients and healthy subjects, and 22 DEMs (FDR<0.10) between PA and non-PA group in bronchiectasis patients. Of these, 10 (5 up-regulated and 5 down-regulated) DEMs in PA group discriminated from both non-PA group and healthy subjects (Figure 2A–C). Figure 2D demonstrates the expression levels of these 10 DEMs in PA and non-PA group. We noted significantly higher levels of miR-628-5p in PA group than in non-PA group (fold-change: 1.72, FDR: 0.03). The expression levels of the overlap DEMs are shown in Table S3.

Figure 2 miRNA expression differentiating PA colonization status. (A) 190 DEMs were differentially expressed between PA group and healthy subjects, and 22 DEMs between PA group and non-PA group. There were 10 overlapped DEMs for the comparison of PA and non-PA group. Five of the DEMs were down-regulated, whereas five were up-regulated. (B) The volcano plot showing the log2(fold-change) and log2FDR of the 190 differentially expressed miRNAs between PA group and healthy subjects. The vertical line was set at the Log2(fold-change) of +1 and −1, respectively. (C) The volcano plot showing the log2(fold-change) and log2FDR of these 22 differentially expressed miRNAs between bronchiectasis patients with PA colonization and those without (non-PA colonization). The vertical line was set at the Log2(fold-change) of +1 and −1, respectively. (D) A heat map demonstrating 10 DEMs (5 up- and 5 down-regulated) after correction for multiple testing. Each column represents an individual patient sample, and each row represents an miRNA. The color code indicates the expression level, with red indicating higher expression levels and blue indicating lower expression levels. (E, F) Receiver operation characteristic curve demonstrating the diagnostic value of the 10 DEMs and the LASSO regression model to differentiate bronchiectasis patients with PA colonization from those without PA colonization. Abbreviations: PA, P. aeruginosa; DEM, differential expressed miRNA; FDR, false-discovery rate; AUC, area under the curve; Down, down-regulated (lower expression levels) in PA group compared with non-PA; Up, up-regulated (higher expression levels) in PA group compared with non-PA; LASSO, least absolute shrinkage and selection operator linear regression analysis.

The Diagnostic Value Of DEMs

Next, we analyzed how well the 10 DEMs discriminated against the PA colonization status (Figure 2E and F). Overall, most miRNAs had the modest diagnostic value in discriminating the PA group from non-PA group. miR-628-5p (C-statistic: 0.71, 95% CI: 0.60, 0.83), miR-223-3p (C-statistic: 0.71, 95% CI: 0.59, 0.83), miR-23a-3p (C-statistic: 0.71, 95% CI: 0.59, 0.83) yielded the highest diagnostic value. Other miRNAs yielded the C-statistics ranging from 0.40 to 0.69. To improve the discriminative capacity, we also performed the LASSO regression model in which we included 74 bronchiectasis patients (>90%) in the training set and eight in the validation set. The accuracy of the LASSO model (see details in Online Supplement) to discriminate PA colonization in the validation set was 87.5%. The LASSO model yielded higher diagnostic value (C-statistic: 0.87, 95% CI: 0.79, 0.95) than any single miRNA in discriminating PA colonization.

DEMs Associated With PA Colonization

We performed multiple linear regression analysis to further validate the correlation between PA colonization and miRNA expression. miR-223-3p, miR-23a-3p, miR-628-5p, miR-92b-5p and miR-548f-3p expression significantly and independently correlated with PA colonization (all P<0.05, Table 2), but the remaining five DEMs did not (all P>0.05, Table S4). We validated the former five DEMs with qPCR assay among 34 bronchiectasis patients who had sufficient paired sputum samples. The correlation coefficient of the readings in duplicate wells was 0.943. miR-92b-5p was significantly down-regulated in PA group (fold-change: −2.94, P<0.05), while miR-223-3p trended towards an increased expression in PA group (fold-change: 2.06, P=0.063) compared with non-PA group (Figure S5, Table S5). In light of the discriminative capacity, the significant correlation with PA colonization, we have included miR-223-3p in pathway analysis despite that it narrowly missed the statistical significance in qPCR assays.

Table 2 Multiple Linear Regression Between The DEMs Associated With PA Colonization And Other Clinical Parameters

miR-92b-5p And miR-223-3p Correlated With The Clinical Parameters In Bronchiectasis

We next explored the correlation between miR-92b-5p and miR-223-3p and the clinical parameters in bronchiectasis (corrected with Benjamini–Hochberg algorithm, Figure 3). The expression level of miR-223-3p correlated negatively with that of miR-92b-5p (r=−0.65 in PA group and −0.61 in non-PA group, both P_adj<0.05). The expression level of miR-223-3p positively correlated with that of sputum IL-8 (r=0.61, P_adj<0.05), and the expression level of miR-92b-5p negatively correlated with that of sputum IL-1β (r=−0.63, P_adj<0.05) and the BSI (r=−0.58, P_adj<0.05) in PA group. Moreover, the expression level of miR-92b-5p negatively correlated with that of sputum IL-1β (r=−0.49, P_adj<0.05). Neither the expression level of miR-223-3p nor miR-92b-5p correlated with the number of bronchiectatic lobes, lung function or the E-FACED scores. We noted a less prominent correlation between the expression level of miR-223-3p and miR-92b-5p and that of sputum IL-1β in non-PA group than in PA group. Nonetheless, the BSI correlated with the expression level of miR-223-3p (r=0.18, P<0.001) but not miR-92b-5p (Table S6). However, there was no marked difference in the expression levels of miR-223-3p and miR-92b-5p between severe and mild-to-moderate bronchiectasis (both P>0.05, Figure S6).

Figure 3 miR-92b-5p and miR-223-3p correlated with the clinical parameters of bronchiectasis. (A) Correlation matrix in PA group. (B) Correlation matrix in non-PA group. In both matrices, the lower left quadrant of the matrix represents the crude P values, and the upper right quadrant of the matrix presents the P values corrected with the Benjamini–Hochberg algorithm. The intersection within the circle represents the adjusted P value >0.05. Blue circle indicate positive correlation, while red circles denote negative correlation. Greater color intensity indicates a greater magnitude of correlation (correlation coefficient). Abbreviations: IL-1β, interleukin-1β; IL-8, interleukin-8; TNF-α, tumor necrosis factor-α; lobe, The number of affected lobes; FEV1%, forced expiratory volume in 1 s; BSI: bronchiectasis severity index; E-FACED, exacerbation-FACED score.

Biological Pathways Of miR-92b-5p And miR-223-3p

Finally, we analyzed the canonical pathways of miR-223-3p and miR-92b-5p by using clusterProfiler package,²² based on the 826 target mRNAs of these two DEMs that were predicted by the DIANA-microT web server version 5.0 and validated with the DIANA-TarBase version 8.0 database.^23,24

Gene Ontology Analysis

We identified the biological processes associated with inflammation, such as leukocyte differentiation, cellular response to oxidative stress. See further details in Figure S7.

KEGG And Genomic Analysis

KEGG analysis revealed 12 canonical pathways, five of which were not further analyzed because they were related to cancer pathways. Table 3 shows the most significant pathways of the annotated mRNAs. The interaction network of mRNAs was classified into four domains: inflammatory processes, hormonal responses, lung aging and stem cell (in which the nuclear factor-κB [NF-κB] signaling pathway, phosphatidylinositol signaling system were the top enriched pathways in the inflammatory processes). See Figure S8 regarding all enriched signaling pathways.

Table 3 The Significant Biological Pathways Of The DEMs Associated With PA Colonization In Bronchiectasis Patients

Discussion

PA is a common opportunistic pathogen which correlates with poor prognosis in bronchiectasis.²⁵ The known mechanisms included airway epithelial injury^26,27 and toxin secretions (such as cyanide and pyocyanin).^28,29 However, little is known regarding other candidate biological pathways underlying the adverse impacts of PA in bronchiectasis. We have detected sputum exosomal DEMs and analyzed their ability to discriminate PA colonization status and the correlation with clinical parameters in bronchiectasis. miR-223-3p and miR-92b-5p, which correlated significantly with sputum inflammatory biomarkers, were implicated in nuclear factor-κB, phosphatidylinositol and longevity regulation pathways. However, neither miR-223-3p nor miR-92b-5p correlated significantly with bronchiectasis severity (the BSI).

Previous studies demonstrated that sputum exosomes could be successfully extracted for biomarker assays,^30–34 which was supported by our findings. Exosomal miRNAs play vital roles in mediating infection and inflammation in many chronic airway inflammatory diseases. For instance, miRNA expression in airway epithelial cells differed significantly between patients with cystic fibrosis and healthy subjects, and between patients with and without PA infection.³⁵ Specifically, miR-223 was up-regulated along with increased IL-8 levels in the cells from patients with ΔF508 mutation³⁵ in those with PA infection. miR-223 reportedly modulated the physiologic function of cystic fibrosis transmembrane conductance regulator and airway IL-8 expression, which might be the target of PA colonization in bronchiectasis. In neutrophilic asthma, miR-223 was up-regulated and positively correlated with sputum IL-1β levels.³⁶ These observations hinted on how PA colonization might interact with the airway inflammation via miRNA modulation in bronchiectasis.

The inflammatory signaling pathways and the accelerated lung senescence might be associated with PA colonization in bronchiectasis. Upon PA infection, airway epithelial-derived exosomes³⁷ stimulated toll-like receptor 4/MyD88/NF-κB signaling pathways.⁷ Increased levels of IL-8 of epithelial cells have also been detected in cystic fibrosis patients with PA infection.³⁸ Congruent with these findings, our findings hinted on the potential link among PA colonization, altered miRNA expression, and NF-κB signaling pathway in bronchiectasis. Our findings suggested the significant correlation between miR-223-3p expression and sputum IL-8 levels, and between miR-92b-5p expression and sputum IL-1β levels. Phosphatidylinositol signaling pathway might also play a role in modulating inflammation associated with PA colonization in bronchiectasis. In fact, PA infection reportedly triggered the formation of airway epithelial cell protrusions leading to activation of phosphatidylinositol and NF-κB signaling pathways.³⁹ Meanwhile, the toxin pyocyanin released by PA could mediate oxidative stress of airway epithelial cells by inducing the expression of antioxidative gene via the epidermal growth factor receptor-PI3K signaling pathway,⁴⁰ which implied that phosphatidylinositol signaling pathway might play a role in bronchiectasis patients with PA colonization.

Apart from the inflammation-related pathways, longevity regulation pathway (ageing) might have also affected the prognosis of bronchiectasis associated with PA colonization. Ageing of the immune system predisposed to infection of opportunistic pathogens such as PA.⁴¹ In turn, PA infection accelerated cell senescence via NF-κB signaling pathway, which further hampered the eradication of PA.⁴² Collectively, our findings provide some hints on the complex network underlying the pathogenesis of bronchiectasis associated with PA colonization.

Our findings have clinical implications. Identification of sputum exosomal miRNA helps explore the important biological pathways underlying the pathogenesis of bronchiectasis, thus miRNAs and their pathways might serve as the complementary tool to explore for the candidate targets for future interventions of PA colonization. Apart from canonical inflammatory pathways, our results have implied the modulation of longevity regulation pathway which opens a new avenue for exploring how PA colonization interacts with the airway epithelium. The significant correlation between sputum inflammatory biomarkers and miR-92b-5p and miR-223-3p provided further evidence on the unresolved inflammation in the PA-colonized microenvironment.

There are notable limitations of our study. First, we compared miRNA expression in spontaneous sputum in bronchiectasis and in induced sputum in healthy subjects. Because all bronchiectasis patients could spontaneously produce sputum, we did not collect induced sputum in bronchiectasis patients. Salivary contamination has been minimized based on our experimental collection procedures. Second, the differential expression of miR-223-3p between PA and non-PA group did not reach to statistical significance (P=0.063) according to the qPCR assays given the insufficient sputum volume in the remaining patients. However, based on the discriminative capacity and the significant correlation with PA colonization and the fold-change of differential expression in both sequencing and qPCR assays, our findings of miR-223-3p expression patterns remained valid. We did not investigate the interactions between miRNA and mRNA. Besides, the association among airway neutrophilic inflammation, miRNA expression and PA colonization in bronchiectasis merits further investigation. Despite the efforts of identifying patients with PA colonization during longitudinal follow-up, we might have still underestimated the risk of misclassification (false-negative cases). Finally, all biological pathways were explored via bioinformatic analysis, further confirmatory experiments are warranted.

In summary, we have identified the differentially expressed sputum exosomal miRNA (miR-223-3p and miR-92b-5p) in bronchiectasis patients with PA colonization. The activation of inflammation and aging regulation pathways might provide clues for the intervention of PA colonization in bronchiectasis.