Circulating microRNA profiling is altered in the acute respiratory distress syndrome related to SARS-CoV-2 infection

One of the hallmarks of SARS-CoV-2 infection is an induced immune dysregulation, in some cases resulting in cytokine storm syndrome and acute respiratory distress syndrome (ARDS). Several physiological parameters are altered as a result of infection and cytokine storm. Among them, microRNAs (miRNAs) might reflect this poor condition since they play a significant role in immune cellular performance including inflammatory responses. Circulating miRNAs in patients who underwent ARDS and needed mechanical ventilation (MV+; n = 15) were analyzed by next generation sequencing in comparison with patients who had COVID-19 poor symptoms but without intensive care unit requirement (MV−; n = 13). A comprehensive in silico analysis by integration with public gene expression dataset and pathway enrichment was performed. Whole miRNA sequencing identified 170 differentially expressed miRNAs between patient groups. After the validation step by qPCR in an independent sample set (MV+  = 10 vs. MV− = 10), the miR-369-3p was found significantly decreased in MV+ patients (Fold change − 2.7). After integrating with gene expression results from COVID-19 patients, the most significant GO enriched pathways were acute inflammatory response, regulation of transmembrane receptor protein Ser/Thr, fat cell differentiation, and regulation of biomineralization and ossification. In conclusion, miR-369-3p was altered in patients with mechanical ventilation requirement in comparison with COVID-19 patients without this requirement. This miRNA is involved in inflammatory response which it can be considered as a prognosis factor for ARDS in COVID-19 patients.

www.nature.com/scientificreports/ cytolytic activity, likely due to genetic factors and other acquired conditions, this may lead to the inability of NK and cytolytic CD8 + T-cells to lyse infected and activated antigen-presenting cells. This may result in prolonged and exaggerated interactions between innate and adaptive immune cells with secretion of many pro-inflammatory cytokines, including Tumor Necrosis Factor (TNF), interferon-γ, interleukin (IL)-1, IL-6, IL-18, and IL-33 in an uncontrolled manner causing a cytokine storm, ARDS, and multiorgan failure 3 .
One of the hallmarks of SARS-CoV-2 infection is an induced immune dysregulation. Severe lymphocytopenia is a very early sign of the disease, preceding lung involvement, which tends to normalize as the patient improves 4 . In parallel, monocytes and macrophages are increased, contributing to raised levels of proinflammatory cytokines such as IL-6, IL-1, TNF-α, and IL-8, which in some patients turn out to be a cytokine storm. The great majority of the inflammatory cells infiltrating the lungs are monocytes and macrophages 5 .
There are many potential factors regulating this immune response. Among them, microRNAs (miRNAs) are small non-coding RNA molecules that are partially complementary to messenger RNA and downregulate gene expression in a variety of manners 6 . Interestingly, miRNA expression is modified by different deleterious processes such as viral infection 7 and several immune system alterations 8 .
In this context, patients who require invasive mechanical ventilation display a general deterioration involving organic failure. A number of physiological parameters are altered as a result of infection and cytokine storm. Among them, miRNA levels might reflect this poor condition since they play a significant role in immune cellular performance including inflammatory responses [9][10][11] . Few studies attempted to identify microRNAs involved in the ARDS related to COVID-19. Among them, a study performed by Tang et al. 11 found several miRNAs specific for severe COVID-19 which may serve as potential biomarkers and therapeutic targets such as miR-15b-5p. Identification of dysregulated miRNAs during ARDS may help to understand molecular mechanisms and immunity pathways involved in regulation of excessive inflammatory responses and in the pathogenesis of many human inflammatory diseases to develop potentially useful treatments.

Patients.
A case-control study was conducted with COVID-19 patients admitted in the Hospital del Mar, Barcelona, Spain. All participants (n = 48) were hospitalized in COVID wards due to severe COVID-19 symptoms and a subset were admitted to Intense Care Unit (ICU) according to criteria listed in supplemental material. Of them, 28 patients were included in the discovery sample and analysed by next generation sequencing (NGS): Fifteen patients required mechanical ventilatory suport (MV+) and they were admitted to ICU whereas the other thirteen patients remained in COVID ward without mechanical ventilation (MV−). Co-morbidities and cumulative risk factors of each individual patients was reported in Supplementary Table 1. Twenty patients were included in the validation step by qPCR: 10 MV+ and 10 MV−.
All hospitalized patients received the hospital's standard-of-care therapy at the time of the study (first wave of COVID-19 pandemic in Spain), consisting in hydroxychloroquine 400 mg/24 h first day and 200 mg/24 h 4 days with azithromycin 500 mg/24 h 3 days, plus ceftriaxone 1 or 2 g/24 h 7 days when there was bacterial superinfection. Patients with severe or critical conditions of pulmonary affectation or suspected cytokine storm were additionally treated with dexamethasone bolus (20 mg/day × 4 days) according to hospital protocol. Demographic and anthropometric variables, full medical history and blood samples were obtained for all participants at time of hospital admission. The blood draw was done systematically in fasting between 8 and 10am. A tube BD Vacutainer® Plus Serum Tubes / CAT Tubes that are coated with silicone and micronized silica particles to accelerate clotting were used. We followed the manufacturer instructions for sample preparation and processing in all the samples. In addition, the same personnel obtained the samples and processed them to minimize interindividual variability.
Statistical analysis of clinical variables was done using SPSS23. Welch Two Sample t-test and Fisher's Exact Test were used to make pairwise comparisons between groups of quantitative and qualitative variables respectively. The sum risk factors, which encompasses all comorbidities described in Table 1, was analysed by Wilcoxon rank sum test. Results with p < 0.05 were considered statistically significant.
The study protocol was approved by the ethics committee of Parc de Salut Mar (CEIm 2020/9232/I) and it was carried out in accordance with the Declaration of Helsinki and the Good Clinical Practice guidelines of the International Conference on Harmonization. All patients were verbally informed about the treatments, by formally obtaining their consent, and its acceptance was recorded in the electronic medical record of the Hospital.

Discovery stage
Next generation sequencing (NGS) of serum microRNAs. miRNA isolation, library construction and sequencing. RNA was extracted from serum samples using miRNeasy Serum/Plasma kit (QIAGEN). MiRNA quantity and purity were determined on Qubit miRNA assay and RNA integrity was assessed using Agilent 2100 Bioanalyzer (Agilent Technologies). Small RNA libraries were prepared by processing 16 13 and posterior counting of unique UMIs with UMI tools version 0.5.4 (count, default method "directional"). The differential gene expression analysis (DEG) was assessed with voom + limma in the limma package version 3.44.1 14 and using R version 4.0.0. Genes having less than 3 counts in at least 2 samples were excluded from the analysis. Raw library size differences between samples were treated with the weighted "trimmed mean method" TMM 15 implemented in the edgeR package 16 . The normalized counts were used in order to make unsupervised analysis, PCA, clusters and heatmaps. For the differential expression (DE) analysis, read counts were converted to log2-counts-per-million (logCPM) and the mean-variance relationship was modelled with precision weights using voom approach in limma package using a linear model design adjusting for the sum of risk factors. The AUC (Area Under The Curve)-ROC (Receiver Operating Characteristics) curve for miR-369-3p was obtained using SPSS 22 and STATA software.
Validation stage qPCR analysis of candidate miRNAs. All experiments were conducted by QIAGEN Genomic Services.
Sample preparation. An aliquot of 200 μl of serum was transferred to a FluidX tube and 60 μl of Buffer RPL containing 1 μg carrier-RNA and RNA spike-in template mixture were added to the sample and mixed for 1 min and incubated for 7 min at room temperature, followed by addition of 20μL Buffer RPP. Total RNA was extracted from the samples using miRNeasy Serum/Plasma Advanced Kit; high-throughput bead based protocol v.1 (Hilden, Germany) in an automated 96 well format. The purified total RNA was eluted in a final volume of 50 μl. miRNA real-time qPCR 2 μl RNA was reverse transcribed in 10 μl reactions using the miRCURY LNA RT Kit (QIAGEN). www.nature.com/scientificreports/ cDNA was diluted 50× and assayed in 10 μl PCR reactions according to the protocol for miRCURY LNA miRNA PCR; each miRNA was assayed once by qPCR on the miRNA Ready-to-Use PCR, Custom panel using miRCURY LNA SYBR Green master mix. Negative controls excluding template from the reverse transcription reaction was performed and profiled like the samples. The amplification was performed in a LightCycler® 480 Real-Time PCR System (Roche) in 384 well plates. The amplification curves were analyzed using the Roche LC software, both for determination of Cq (by the 2nd derivative method) and for melting curve analysis.
Hemolysis was tested in all samples using miR-451 (is expressed in red blood cells) and miR-23a-3p (relatively stable in serum and plasma and not affected by hemolysis). dCq(23a-451) was lower than 7 in all samples showing minimal signs of hemolysis.
Data analysis. The amplification efficiency was calculated using algorithms similar to the LinReg software. All assays were inspected for distinct melting curves and the Tm was checked to be within known specifications for the assay. Furthermore, assays must be detected with 5 Cq less than the negative control, and with Cq < 37 to be included in the data analysis. All data was normalized with hsa-miR-502-3p and hsa-miR-23a-3p using the average of the assays detected in all samples (n = 20 samples).
Integration with public gene expression dataset. Overmyer et al. 17 performed RNA-Seq from leukocytes from 102 patients with COVID-19 (51 patients needed for admission into ICU and 51 hospitalized in COVID-19 floor). We downloaded the count matrix from GEO (GSE157103) and filtered genes having less than 10 counts in at least 50 samples. Read counts were converted to log2-counts-per-million (logCPM) and the mean-variance relationship was modelled with precision weights using voom approach in limma package adjusting by sex and age (as done in the original publication). The contrast performed for this analysis was ICU. requirement vs. no ICU.
Genes were considered differentially expressed (DE) if the adjusted p value < 0.05 and |FC| > 1.5. To find which DE genes could potentially be regulated by the miR-369-3p, we assumed that down-regulation in miR-369-3p will lead to up-regulation in genes. Hence, we retrieved targets for miR-369-3p from miRTarBase (v.8), miRDB (v.6) and TargetScan (v.7.2) and intersected these targets with DE genes (up-regulated in ICU patients).
Finally, we performed enrichment analysis of the list of intersected genes implemented in clusterProfiler 18 R package version 3.18.0 using the molecular signatures database collection c5.bp 19,20 : Gene sets derived from the Biological Process Gene Ontology (GO), version 7.2. Gene sets were considered significantly enriched if the adjusted p value < 0.05 (after multitesting correction 21 ).

Results
Sample description. The anthropometric features of the MV+ and MV− groups are shown in Table 1.
Statistical differences in high blood pressure and obesity were observed between both groups. In addition, an increased inflammatory component was detected in MV+ group according to IL-6 levels measured a baseline. Otherwise, there were not significant differences between groups at time of hospital admission in the different cell counts, which included: Leukocytes, Neutrophils, Lymphocytes, Eosinophils, Hemoglobin, Hematocrit, CD4 and CD8 cells.
Finally, the miR-369-3p was validated (Supplementary Table 3): it was found under-expressed in MV+ patients (Fold change − 2.7; p value = 0.0061). ROC curves and AUCs were used to assess the discriminative accuracy of the miR-369-3p (Fig. 1). The AUC (95% CI) for discriminating MV− vs. MV+ patients was 0.72 (0.53-0.91), standard error = 0.098, and p value = 0.045. Using a cut-off point of 2.75 in miR-369-3p levels according to NGS results, the sensitivity was 77% with a specificity of 53.3%. Pathway enrichment analysis. Pathway enriched results for the miR-369-3p displayed diverse significant signaling pathways, some of them related to immune system: T-Cell Receptor and Co-stimulatory Signaling (p value = 0.002), IL-4 Signaling Pathway (p value = 0.008), Notch signaling pathway (p value = 0.035), and TGF-beta signaling pathway (p value = 0.003). Figure 2 shows predicted targets for the miR-369-3p using Chair for Bioinformatics at the University of Saarland 27 . MirPath v.3 analysis, considering TarBase v7.0 that provides high-quality manually curated experimentally validated miRNA:gene interactions, the most significant pathway was circadian entrainment (p value < 0.001) but the Hedgehog signaling pathway was also significant (p value = 0.016). After integrating with gene expression results from GSE157103 (Supplementary Table 4), the most significant GO enriched pathways were fat cell differentiation, regulation of transmembrane receptor protein Ser/Thr, regulation of biomineralization and ossification, and acute inflammatory response (Fig. 3).

Discussion
The aim of this study was to identify microRNAs altered in ARDS in the context of SARS-CoV-2 infection. Circulating miRNAs in patients who underwent ARDS and needed mechanical ventilation were analyzed in comparison with patients who had COVID-19 but without ICU requirement. First of all, NGS results identified 170 miRNAs differentially expressed between patient groups in a discovery sample. Twenty miRNAs with the best hit scores were chosen for validation by qPCR in an independent set of samples. Finally, the miR-369-3p was found significantly reduced in patients with ARDS. In silico analysis of this miRNA predicted a putative role in the regulation of immune pathways for example by targeting genes in the T-Cell Receptor and Co-stimulatory Table 2. The top ten upregulated and top ten downregulated miRNAs in MV+ patients compared to MV− patients in the discovery sample using NGS. miRNA name: Annotations related to the gene (miRNA) according to miRBase v22.1, genome version GRCh38. p value obtained with the moderated t-test. Chr chromosome, logFC log2 fold change, MV mechanical ventilation requirement. www.nature.com/scientificreports/ signaling, IL-4 signaling, Notch signaling, and TGF-beta signaling pathways. Considering only validated targets, it is interesting to highlight the Hedgehog signaling pathway that is involved in a wide variety of functions and in some diseases like pulmonary fibrosis 28 . Moreover, differentially expressed genes between ICU requirement vs. no ICU patients (GSE157103) that are potentially regulated by the miR-369-3p are enriched in pathways involved in the immune system like Serine/threonine kinase receptors which are mediators of TGF-beta family signals 29 , and in the acute inflammatory response. Previous literature validated the miR-369-3p in its role in the immune system. For example, miR-369-3p was found upregulated in severe and/or recurrent Herpes Simplex Virus infection 30 . In this study, miR-369-3p inhibition in NK-92 cells resulted in profound upregulation of 4 genes (APOBEC3G, MAP2K3, MAVS and TLR7) and downregulation of 36 genes taking part in antiviral response or associated with signaling pathways of Toll-like receptors (TLR), NOD-like receptors, the retinoic acid-inducible gene I (RIG-I)-like receptors (RLRs) and type I IFN-related response 30 .
Accordingly, miR-369 inhibition significantly increased cell apoptosis and inflammatory factor production triggered by hypoxia by targeting TRPV3 31 .
Another study demonstrated that the upregulation of miR-369-3p suppresses the LPS-induced inflammatory response, reducing C/EBP-β, TNFα and IL-6 production. Moreover, it was suggested that miR-369-3p plays a key role in negatively regulating the LPS-induced Dendritic cells responses mainly targeting iNOS. They found that the upregulation of miR-369-3p decreased the release of TNFα, IL-6, IL-12, IL-1α, IL-1β in response to LPS, and increased the production of anti-inflammatory cytokines such as IL-10 and IL-1RA 32 . Furthermore, a downregulation of miR-369-3p has been previously reported in Crohn's disease plasma compared with control plasma 33 . It is noteworthy that all these findings are in line with our results where miR-369-3p is significantly decreased in patients with poorer COVID-19 prognosis suggesting this microRNA as a potential central player in the inflammatory response. Accordingly, the MV+ group had significantly higher IL-6 making this hypothesis more plausible and showing a potential association between these facts. Interestingly, a recent study found that miR-369-3p expression was decreased in the lung tissues of mice with induced idiopathic pulmonary fibrosis. By silencing the lncRNA DLEU2, the idiopathic pulmonary fibrosis was suppressed through upregulating miR-369-3p expression and downregulating TRIM2 expression 34 .
On the other hand, the miR-369-3p has a target site in the SARS-CoV-2 genome 35 that could suggest an antiviral activity of this miRNA. Hence, a dual action of this miRNA, both on immune system and viral performance could be involved in the COVID-19 severity. Interestingly, the miR-369-3p does not target other SARS viruses even though they share approximately 78.7% of the genome 35 .
As a limitation, we cannot know if the under-expression of miR-369-3p in MV+ patients is cause or consequence of the condition, since we only analyzed the microRNA signature at a single time point, therefore we only have a vision of this particular moment. Noteworthy, samples were obtained at hospital admission, prior to ARDS and ICU requirement. Hence, our study points out that miR-369-3p is altered in severe COVID-19 patients prior developing ARDS and therefore it can play a role in this context. Thus, in an intra-hospital setting, miR-369-3p could be used as a potential prognostic marker in a time window between hospital admission and ICU requirement. Further studies should be addressed to elucidate whether miR-369-3p was downregulated prior to cytokine storm and therefore whether miR-369-3p can predict the transition from a severe to a critical phase of the disease. www.nature.com/scientificreports/ Moreover, since COVID-19 severity is highly associated with other comorbidities, as for example obesity and high blood pressure, we cannot rule out that miR-369-3p levels were also associated with a multiple comorbidity/ risk factors state, instead specifically associated with COVID-19-related ARDS. However, results from discovery sample were adjusted for the sum of risk factors in an attempt to minimize the effect of co-morbidities and a validation sample was chosen with minimum differences between groups, or at least, with no statistical differences in the most known COVID-19-related comorbidities.
These results and previous studies collectively contribute to a better knowledge of miR-369-3p and their target pathways which may be considered a potential target for new molecular therapeutic approaches 36,37 . Importantly, further studies should focus on pathway regulation by miRNAs rather than individual targets, in a cell or tissue context. Hence, new miRNA therapeutic tools could be designed in a tissue-directed approach.
It should be pointed out that although the other microRNAs found differentially expressed in the NGS step were not validated in the independent sample set, we cannot rule out that they could be involved in COVID-19 symptomatology, and patients MV+ had a particular circulating microRNA signature with more than one altered microRNA. Further analyses in a larger sample size should be addressed.
In this regard, a similar study (CIBERESUCICOVID) that included 84 participants identified a signature of three miRNAs (miR-148a-3p, miR-451a, and miR-486-5p) that distinguished between ICU patients and ward patients 38 . MiR-451a was also found to be altered in our discovery samples (Supplementary Table 2), suggesting a possible role of this miRNA in the severity of COVID-19. The different findings between both studies can be attributed to differences between the methodologies used to identify microRNAs (RT-PCR of candidate genes versus NGS of all circulating microRNAs).
In conclusion, miR-369-3p was found decreased in serum samples of COVID-19 patients with mechanical ventilation requirement compared to COVID-19 patients without this requirement. Along with in silico analyses  www.nature.com/scientificreports/