MicroRNA mapping of bronchial aspirate for molecular phenotyping and prognostication in patients on mechanical ventilation

Abstract

The application of microRNA (miRNA) profiling in respiratory biospecimens, particularly bronchial aspirate (BAS), remains underexplored. Here, we aimed to validate and refine miRNA quantification in BAS samples to establish its suitability for molecular phenotyping. This was a multicenter study including 288 COVID-19 patients on invasive mechanical ventilation. Respiratory biospecimens included BAS, tracheal aspirate, and bronchoalveolar lavage fluid samples. A predesigned miRNA panel was evaluated using RT-qPCR. Biomarker evaluation and functional assessment were subsequently conducted. An initial technical validation phase corroborated the reproducibility of miRNA profiling in BAS samples. Comparative analyses of miRNA expression profiles across respiratory samples revealed distinct miRNA patterns among biospecimens. In the biomarker analysis, two miRNA ratios, miR-34c-5p/miR-34a-5p and miR-34c-5p/miR-125b-5p, were inversely associated with intensive care unit (ICU) survival (hazard ratio [HR]: 0.18 and 0.17, respectively) during the discovery phase. Risk and survival analyses in the test phase confirmed the reproducibility of the miR-34c-5p/miR-34a-5p ratio (hazard ratio [HR] = 0.17). Functional analyses revealed the utility of miRNA profiling in BAS for identifying pathogenic pathways and developing therapeutic targets. Overall, these findings position miRNA profiling in BAS samples as a valuable approach for biomarker discovery, identification of pathophysiological mechanisms, and development of targeted pulmonary therapies.

Keywords: COVID-19; MT: Non-coding RNAs; SARS-CoV-2; biomarkers; bronchial aspirate; bronchoalveolar lavage fluid; invasive mechanical ventilation; microRNAs; noncoding RNA; tracheal aspirate.

Graphical abstract

Figure 1

MicroRNA profiling of respiratory samples (A) Bland-Altman plot comparing the microRNA expression levels of independent aliquots from each patient in the technical validation sample set (n = 57). Each point represents a patient. The bold and dashed lines represent the mean difference and ±1.96 SD, respectively. In this step, the expression levels were normalized using the external reference miRNA cel-miR-39-3p to avoid technical variability: ΔCq = Cq_microRNA−Cq_{cel-miR-39-3p}. (B) Dual-axis chart including the microRNA expression levels of bronchial aspirate (BAS) samples (n = 57), tracheal aspirate (TA) samples (n = 42), and bronchoalveolar lavage fluid (BALF) samples (n = 96). The bars represent the percentage of expression of each microRNA in the different respiratory specimens (expression rate axis), and the points represent the median of the quantification cycle (Cq) of each microRNA (Cq axis). Expression rate of each miRNA was calculated as the percentage of samples with detectable expression (Cq < 35). The orange line indicates the Cq threshold of 35, above which miRNAs were considered undetectable.

Figure 2

MicroRNA profile of bronchial aspirate samples from survivors and nonsurvivors admitted to the ICU in the discovery cohort (n = 61) (A) Volcano plot showing the differences in the miRNA ratios between survivors and nonsurvivors in the discovery cohort (n = 61) after adjusting for age, sex, and batch. Each point represents a detected microRNA. The light green dots represent the significantly expressed microRNA ratios with a p value <0.05. Dark green dots represent the microRNA ratios with an FDR <0.10. The dashed lines are set at FDR = 0.10 and FC = 0.67 and 1.5. (B) Boxplot including bronchial aspirate (BAS) levels of candidate microRNA ratios between study groups. Each point represents a patient. (C) Boxplot showing the individual microRNAs that compose each differentially expressed microRNA ratio between survivors and nonsurvivors. Each point represents a patient. The gray lines connect both microRNAs for the same patient. (D) GAM modeling for risk mortality and the levels of candidate microRNA ratios. Effective degrees of freedom (edf) and p values are displayed. The odds ratio represents the risk change per 1 SD in continuous predictors. (E) Volcano plot showing differential microRNA expression in a human postmortem lung biopsy RNA-seq dataset (GEO: GSE235130) between COVID-19 patients and non-COVID-19 patients. The differential expression criteria were set at FC ≥ 1.5 and FDR <0.05. Individual microRNAs of candidate microRNA ratios are labeled. (F) Boxplot of the candidate microRNA ratio in BAS samples measured in tracheal aspirate (TA) samples. (G) Boxplot of the candidate microRNA ratios in BAS samples measured in bronchoalveolar lavage fluid (BALF) samples. Ratios are represented as ΔCq values, calculated as the difference between the Cq of the numerator miRNA and the Cq of the denominator miRNA in each ratio (ΔCq = Cq_numerator−Cq_denominator).

Figure 3

Biomarker potential of candidate microRNA ratios of bronchial aspirate samples to discriminate survivors and nonsurvivors admitted to the ICU (A) Kaplan-Meier estimations of candidate microRNA ratios in the discovery cohort (n = 61). The cutoff values, hazard ratios (HRs), and p values are displayed. (B) Kaplan-Meier estimations for the candidate microRNA ratio in the test population (n = 32). The cutoff, hazard ratio (HR), and p value are displayed. (C) Discrimination value for the miR-34c-5p/miR-34a-5p ratio and the main clinical predictors in the whole population (n = 93). The data are presented as the C-index. Clinical predictors are divided according to two temporal points: intensive care unit (ICU) admission and invasive mechanical ventilation (IMV) initiation.

Figure 4

Functional analysis of miR-34c-5p and miR-34b-5p targets (A) Volcano plot representing the differential expression of miR-34b-5p and miR-34c-5p targets between SARS-CoV-2-infected and noninfected samples in an external dataset (GEO: GSE210223) of human bronchial epithelial (HBE) cell cultures. The dashed lines are set at FDR = 0.05 and FC = 0.67 and 1.5. Each dot represents a target. The green dots represent the differentially expressed targets. The common targets of the GEO: GSE210223 and GEO: GSE147507 datasets are labeled. The analysis was based on the intersection (I = 2) of the total number of target genes of miR-34c-5p and miR-34b-5p. (B) Volcano plot representing the differential expression of miR-34b-5p and miR-34c-5p targets between SARS-CoV-2-infected and noninfected samples in an external dataset (GEO: GSE147507) of human postmortem lung biopsies. The dashed lines are set at FDR = 0.05 and FC = 0.67 and 1.5. Each dot represents a target. The green dots represent the differentially expressed targets. The common targets of the GEO: GSE210223 and GEO: GSE147507 datasets are labeled. The analysis was based on the intersection (I = 2) of the total number of target genes of miR-34c-5p and miR-34b-5p. (C) STRING protein-protein interaction network. The analysis included the seven upregulated transcripts common in the GEO: GSE210223 and GEO: GSE147507 external databases. Edges indicate both physical and functional associations (interaction score cutoff is set at 0.90). The US Food and Drug Administration (FDA)-approved drugs for the selected targets are displayed. (D) Functional analysis of the top 25 terms ranked by the FDR from the gene ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), and Reactome databases. The Rich factor is the ratio of differentially expressed gene numbers annotated in this pathway term to all gene numbers annotated in this pathway term. (E) Enrichment analysis of lung cell types on the basis of single-cell RNA-seq data from the GTEx Project database. Each column shows a cell type, and each row shows a gene. The point size indicates the number of cells where the gene was detected, and the color represents the expression level.

Figure 5

Study design and main findings