The study included 525 subjects referred to the sleep unit due to suspected OSA (ClinicalTrials.gov identifier: NCT03513926). The recruitment was conducted at University Hospital Arnau de Vilanova-Santa María of Lleida and Hospital San Pedro Alcántara of Cáceres in Spain. Patients aged older than 18 years were recruited for the present study. The exclusion criteria were the presence of a previously diagnosed sleep disorder, a history of CPAP treatment, inability to complete the questionnaires and any medical, social or geographic factors that could compromise the eligibility of the subject, i.e., pregnancy, substance abuse, cancer, renal insufficiency, severe chronic obstructive pulmonary disease, chronic depression and other very limiting chronic diseases or life expectancy less than 1 year. General physical and anthropometric parameters were documented, and information on sociodemographic characteristics, medical history, medication use and lifestyle habits was collected by trained clinicians.

All eligible participants underwent a comprehensive polysomnographic (PSG) sleep study (Sleepware G3, Philips, Amsterdam, Netherlands). All procedures were conducted according to national guidelines and regulations governing clinical practice.27 Trained personnel who adhered to standard criteria28 analyzed the results of the sleep studies. Apnea was characterized as an interruption or reduction in oronasal airflow by ≥90% for at least 10s, while hypopnea was defined as a 30–90% reduction in oronasal airflow for at least 10s, associated with either oxygen desaturation by at least 3% or an arousal on the electroencephalogram. The apnea–hypopnea index (AHI), which is indicative of OSA severity, was calculated based on the average number of apnea and hypopnea events per hour of sleep. Following the PSG study, subjects were categorized according to the International Consensus Document on Obstructive Sleep Apnea29 into non-OSA (AHI<15events/h) and OSA (AHI≥15events/h) groups.

All subjects provided written informed consent for study participation prior to the collection of blood samples. This study was revised and approved by the Clinical Research Ethics Committee of the University Hospital Arnau de Vilanova-Santa María of Lleida (no. 1153/1411). The study was performed according to the Declaration of Helsinki.

The study design is presented in Fig. 1. The experimental approach was composed of three clearly defined phases: (i) Discovery phase: from the initial population of 525 participants with suspected OSA, 27 subjects with confirmed OSA were randomly selected. Subsequently, 26 controls without OSA from the same study population were selected using nearest neighbor propensity score matching by age, sex and body mass index (BMI). miRNA sequencing was conducted on these 53 subjects. The number of samples aligns with the recommendations for transcriptomic profiling studies.30 (ii) Technical validation phase: to confirm the reliability of the miRNA sequencing findings, a rigorous assessment of miRNA candidates was conducted using RT-qPCR, which is acknowledged as the gold standard for circulating miRNA quantification. This technical validation was performed on the same set of samples employed in the discovery phase (n=53). (iii) Comprehensive evaluation: miRNA candidates that successfully passed the technical validation phase were further evaluated in the entire study population, excluding samples utilized in the previous steps (n=472).

Fig. 1.

Study workflow. Figure created with BioRender.com (https://app.biorender.com/), license number: TH26UP0M73. Abbreviations: CPAP: continuous positive airway pressure; OSA: obstructive sleep apnea.

(0.71MB).

Total RNA was extracted from 200μL of frozen plasma utilizing the miRNeasy Serum/Plasma Advanced Kit (Qiagen, Hilden, Germany). To establish a normalization strategy, synthetic Caenorhabditis elegans miR-39-3p (cel-miR-39-3p) (Qiagen) was introduced as an external reference miRNA at a concentration of 1.6×108copies/μL. Additionally, 1μg of MS2 carrier RNA (Roche Diagnostics, Mannheim, Germany) was added to the mixture to increase the yield of extracellular miRNAs. Quality control for RNA isolation was assured by the inclusion of the RNA Spike-In Kit (synthetic UniSp2, UniSp4, and UniSp5) (Qiagen). All reagents were spiked into samples during RNA isolation after incubation with the denaturing solution. The isolated RNA was then eluted in 20μL of RNAse-free water and stored at −80°C.

For miRNA quantification, the miRCURY LNA Universal RT microRNA PCR System (Qiagen) protocol was followed. Reverse transcription (RT) to synthesize complementary DNA (cDNA) was performed using a miRCURY LNA RT Kit (Qiagen) in a total reaction volume of 10μL. The spike-in UniSp6 (Qiagen) was added to monitor the RT reaction. The RT conditions included incubation at 42°C for 60min, inactivation at 95°C for 5min and immediate cooling at 4°C. Subsequently, the cDNA was stored at −20°C. We used the miRCURY LNA miRNA Custom Panels (384-well plates) (Qiagen), which contained the selected miRNAs and spike-in primers (Qiagen). A QuantStudio™ 7 Flex Real-Time PCR System (Applied Biosystems, Waltham, MA, USA) was used for qPCR in a 10μL reaction volume. The RT-qPCR conditions included an initial incubation at 95°C for 2min, followed by 40 cycles of 95°C for 10s and 56°C for 1min and a final melting curve analysis. Synthetic UniSp3 served as an interplate calibrator and qPCR control. Amplification curves were assessed by melting curve analysis using QuantStudio Software v1.3 (Thermo Fisher Scientific, MA, USA), ensuring the presence of single products and the absence of primer dimers.

The quantification cycle (Cq) was determined as the fractional cycle number at which the fluorescence surpassed a defined threshold. Cq values of spike-in RNA templates were initially scrutinized to monitor the homogenous efficiencies of RNA extraction procedures, the robustness of RT and the absence of PCR inhibitors. The ΔCq ratio (miR-23a-3p–miR-451a), as proposed by Blondal et al.,31 was used to exclude hemolyzed samples. Cq values exceeding 35 cycles were considered undetectable and censored at the minimum level observed for each miRNA. miRNAs detected below the limit of detection in more than 80% of the samples were categorized as nonexpressed. Samples that did not pass the quality control test were excluded from subsequent statistical analysis. Relative quantification was accomplished using the 2−Cq method (ΔCq=CqmiRNA−Cqcel-miR-39-3p). The expression levels were log-transformed for subsequent statistical analysis.

We retained the predicted miRNA–target interactions utilizing the web-based TargetScan (Release 8.0, https://www.targetscan.org/vert_80/).32 The search settings were set to all predicted targets, regardless of conservation. To elucidate the biological relevance and pathways associated with the identified miRNA targets, we conducted a functional enrichment analysis using the R package ClusterProfiler 4.6.2.33ClusterProfiler is released through the Bioconductor project and can be accessed via https://bioconductor.org/packages/clusterProfiler/. The Kyoto Encyclopedia of Genes and Genomes (KEGG), Gene Ontology (GO) and Reactome were used as reference databases. To further explore the expression patterns of selected miRNA targets, tissue patterns from the Genotype-Tissue Expression (GTEx) Portal (https://www.gtexportal.org/home/) project v.8 were used. To explore potential therapeutic interventions, we investigated the repositioning potential of upregulated miRNA targets using the Drug–Gene Interaction database (DGIdb v.5.0.5) (https://dgidb.org).34

All the statistical analyses were performed using R version 4.2.2 (R Foundation for Statistical Computing). Descriptive statistics were employed to explore the characteristics of the study population. Continuous variables were compared between groups using the Mann–Whitney U test, whereas categorical variables were compared using Fisher's exact test. The data are presented as frequencies (percentages) for categorical variables and as medians (25th and 75th percentiles) for continuous variables. Spearman correlation was utilized to estimate the correlation between continuous variables. A correlation coefficient (rho) greater than 0.3 was considered as biologically relevant.35 The p-value threshold defining statistical significance was set at <0.05.

Three distinct ML-based feature selection methods were implemented simultaneously: Random forest with Boruta,36 variable selection using random forests (VSURF)37 and sparse partial least squares (sPLS).38 To ensure the robustness and consistency of our selection process, we repeated each method 50 times. In each iteration, miRNAs selected in at least one run were considered important variables. Finally, miRNAs consistently identified across all three methods were designated as candidates for further investigation through bioinformatics analysis, as previously described.39

For the external datasets, we obtained the preprocessed gene expression data and normalized them using quantile normalization. The raw data were processed and normalized using methods implemented in the Limma R package40 for the Illumina and Agilent microarrays and the oligo R package41 for the Affymetrix chips. Subsequently, we applied a base 2 logarithmic scale transformation. Following data normalization, we averaged the expression data from multiple probes associated with the same gene. An extensive outlier analysis was performed to assess the plausibility of extreme values. After thorough evaluation, outliers were excluded from further analysis. Differential expression analysis was performed with processed data between groups of each microarray dataset using the Limma R package. For the longitudinal analysis, we used a paired-samples design to account for interindividual differences with respect to the pretreatment state. For each comparison, we utilized the Limma moderated t-test (implemented with the “lmFit” and “eBayes” functions) to calculate differential expression for each gene.

miRNA sequencing was conducted in specific matched subgroups of subjects with and without confirmed OSA (n=53). The main characteristics of the study groups are detailed in Table 1. miRNA candidates were chosen based on the following criteria: over the 90th decile of expression (to ensure their detection in the following phases), ≥1.3-fold differential expression and p-value≤0.2 (to identify a broader range of potential biomarkers). A total of 38 miRNAs fulfilled these criteria, as depicted in Fig. 2A.

Fig. 2.

Identification of microRNAs related to obstructive sleep apnea. (A) Volcano plot showing the p-value versus the fold change for the microRNAs analyzed in the discovery phase. microRNAs over the 90th decile of expression are depicted in dark gray. The green dots indicate microRNA candidates. (B) Violin plots including microRNAs that fulfilled the selection criteria in the technical validation phase. Plots depict log102−ΔCq as the microRNA expression unit. The gray color indicates the nonobstructive sleep apnea group, and the red color indicates the obstructive sleep apnea group. The plot presents the median (25th and 75th percentiles) estimator of the density (as density curves) and individual values (black dots). Fold changes and p-values are displayed. (C) Correlogram showing spearman correlation coefficients between polysomnographic parameters and the microRNAs that fulfilled the selection criteria in the technical validation phase. The color scale illustrates the degree of correlation and ranges, indicating the negative to positive correlations. Abbreviations: AHI: apnea–hypopnea index; OSA: obstructive sleep apnea; SaO2: oxygen saturation; TSat90: time during the night with SaO2 below 90%.

(0.95MB).

To validate the miRNA sequencing results, the selected candidates were further examined through RT-qPCR using the same sample set (n=53). Eight samples were excluded from further analysis due to failing the quality control test, which was attributed to the presence of hemolysis and/or low quality of the spike-ins. miR-6877-5p was expressed at low levels (Cq≥35 in more than 80% of the samples) and was therefore excluded from further analysis. Among the remaining 37 miRNA candidates, those with a ≥1.3-fold difference between study groups, with a p-value≤0.01 and that correlated with the AHI (rho>0.3) were retained (Fig. 2B, C and Supplemental Table S1). Ultimately, 11 miRNAs—let-7d-5p, miR-15a-5p, miR-15b-5p, miR-18a-5p, miR-20b-5p, miR-30e-5p, miR-107, miR-199a-5p, miR-425-5p, miR-484 and miR-625-5p—were identified as candidates for the subsequent feature selection process. With the exception of the AHI and the apnea index, weak correlations were observed between the miRNA panel and the parameters obtained from PSG (Fig. 2C).

Following the identification of miRNA candidates associated with OSA, quantification was performed using RT-qPCR in the entire study population. A total of 9.5% of the samples did not pass quality control due to hemolysis or variability in spike-ins and were discarded from further analysis. The detailed characteristics of the study population are presented in Table 2. As anticipated, patients with OSA were older, predominantly men, had a higher BMI and had a more unfavorable profile of comorbidities and disease-related sleep parameters.

Due to the intricate dynamics of miRNA networks, the circulating levels of miRNAs were entered into a supervised feature selection protocol that integrates three distinct ML-based methods (Fig. 3A). The outcome of each method resulted in a unique feature subset: (i) VSURF selected eight candidates: let-7d-5p, miR-15a-5p, miR-18a-5p, miR-30e-5p, miR-107, miR-199a-5p, miR-425-5p and miR-484; (ii) sPLS selected six candidates: let-7d-5p, miR-15a-5p, miR-15b-5p, miR-20b-5p, miR-107 and miR-625-5p; and (iii) Boruta selected all candidates. Notably, let-7d-5p, miR-15a-5p and miR-107 demonstrated robust and consistent selection across the three methods.

Fig. 3.

Identification of microRNA-mediated molecular mechanisms associated with obstructive sleep apnea. (A) Identification of microRNAs associated with obstructive sleep apnea using a consensus of three supervised machine learning feature selection algorithms (Boruta, VSURF and sPLS) with the intersected microRNAs among the methods selected as candidates. (B–D) Functional enrichment analysis of the predicted downstream target genes using KEGG (B), GO (C) and Reactome (D) in the R package ClusterProfileR 4.0. The plots present the rich factors of the pathways with the most significant differences, considering their p-values. The intensity of the colors of the bars denotes the p-values. The rich factor consists of a ratio of target genes annotated in the molecular processes to all genes annotated in the processes.

(0.94MB).

Target prediction for let-7d-5p, miR-15a-5p and miR-107 utilizing TargetScan v8.0 revealed a total of 3436 transcripts as potential targets (1993 targets for let-7d-5p, 1605 targets for miR-15a-5p and 101 targets for miR-107), with 257 targets shared by at least two of the selected miRNAs. To elucidate the underlying molecular pathways and biological processes associated with these identified target transcripts, an exhaustive enrichment analysis was conducted employing the KEGG, GO and Reactome databases. The analyses identified 9 KEGG pathways, 220 GO processes and 15 Reactome pathways (Fig. 3B–D and Supplemental Tables S2–S4). The predominant significant pathways and processes were closely linked to cellular events such as cell death, cell differentiation, extracellular remodeling, autophagy and metabolism.

To assess the robustness of our findings, we systematically validated the results by analyzing the expression of the target transcripts in external datasets obtained from GEO. This process was divided into two steps. Initially, we conducted a comparative examination of miRNA targets between individuals with and without OSA. Subsequently, we analyzed the influence of CPAP therapy on the identified miRNA targets.

First, three distinct external datasets, including OSA patients and non-OSA subjects, were evaluated (GSE226379, GSE135917 and GSE75097). Among the targetomes of let-7d-5p, miR-15a-5p and miR-107, one gene, the transcription factor Dp-2 (TFDP2), which is a target of let-7d-5p and miR-15a-5p, exhibited significant differential expression in the three comparisons (fold change≥1.5 and p-value≤0.05) (Fig. 4A–C). To further elucidate the expression profile of the miRNA targets, we retrieved data from the GTEx project, which revealed that TFDP2 is expressed in most human tissues (Fig. 4D).

Fig. 4.

External datasets. (A–C) Volcano plots showing the p-value versus the fold change for the predicted targetome in three external RNA sequencing datasets obtained from the GEO database from patients with obstructive sleep apnea. Significantly differentially expressed genes are presented in green (fold change>1.5 and p-value≤0.05). Intersected genes among the datasets are labeled. (A) Plasma, GSE226379; (B) subcutaneous fat, GSE135917; and (C) peripheral blood mononuclear cells (PBMCs), GSE75097. (D) Tissue enrichment analysis using GTEx (https://www.gtexportal.org/home/). Each column denotes a tissue type, and the row represents the intersected gene among the external datasets. The size and color of the points reflect the expression level of the gene, represented as normalized transcripts per million (nTPM) values. The GTEx Project was supported by the Common Fund of the Office of the Director of the National Institutes of Health and by the NCI, NHGRI, NHLBI, NIDA, NIMH, and NINDS. (E, F) Volcano plots showing the p-value versus the fold change for the predicted targetome in two external RNA sequencing datasets, including obstructive sleep apnea patients allocated to CPAP treatment, obtained from the GEO database. Intersected genes among the A–C datasets are labeled. (E) PBMCs, GSE133601; and (F) peripheral blood leukocytes, GSE49800.

(0.88MB).

To explore the impact of the current therapeutic strategy, the gene expression levels of TFDP2 were evaluated in two independent external datasets (GSE133601 and GSE49800), including samples from patients with OSA collected before and after CPAP therapy. No alterations in the transcript levels were observed following therapy (fold change≥1.5 and p-value≤0.05) (Fig. 4E and F). Next, we evaluated the repositioning potential of the miRNA target using DGIdb. No FDA-approved drugs were identified.

The current investigation has several strengths, including the adoption of a “real clinical context”, multicentric design, comprehensive clinical assessments, the analysis of the plasma miRNome and the employment of diverse feature selection algorithms. However, it is imperative to acknowledge the inherent limitations of this study. Although the discovery phase included a subset of 53 subjects, which may appear limited, this sample size is relatively large compared to typical RNA sequencing studies. Nevertheless, we acknowledge that this limitation may affect the robustness of our findings. To further investigate the relevance of our results, we utilized different external datasets. However, validation in additional cohorts with more diverse demographic and clinical characteristics and an appropriate sample size to ensure adequate statistical power is essential. A significant limitation of our study is the lack of mechanistic investigations to confirm the biological roles of the identified miRNAs and their target genes in the pathogenesis of OSA. While we have established associations between specific miRNAs and OSA using various datasets, experimental studies are vital to elucidate the causal relationships and underlying signaling pathways influenced by these miRNAs. Finally, the criteria for selecting candidate miRNAs in the discovery and technical validation phases were intentionally permissive. This approach was chosen for two reasons. First, the inherent biology of miRNAs allows them to regulate biological responses even with minor fluctuations in their levels,26 which justifies the use of lower fold-change thresholds. Second, the primary goal of these phases was to “feed” the feature selection analysis rather than to identify individual miRNAs specifically associated with OSA.