Alteration of circulating ACE2-network related microRNAs in patients with COVID-19

Angiotensin converting enzyme 2 (ACE2) serves as the primary receptor for the SARS-CoV-2 virus and has implications for the functioning of the cardiovascular system. Based on our previously published bioinformatic analysis, in this study we aimed to analyze the diagnostic and predictive utility of miRNAs (miR-10b-5p, miR-124-3p, miR-200b-3p, miR-26b-5p, miR-302c-5p) identified as top regulators of ACE2 network with potential to affect cardiomyocytes and cardiovascular system in patients with COVID-19. The expression of miRNAs was determined through qRT-PCR in a cohort of 79 hospitalized COVID-19 patients as well as 32 healthy volunteers. Blood samples and clinical data of COVID-19 patients were collected at admission, 7-days and 21-days after admission. We also performed SHAP analysis of clinical data and miRNAs target predictions and advanced enrichment analyses. Low expression of miR-200b-3p at the seventh day of admission is indicative of predictive value in determining the length of hospital stay and/or the likelihood of mortality, as shown in ROC curve analysis with an AUC of 0.730 and a p-value of 0.002. MiR-26b-5p expression levels in COVID-19 patients were lower at the baseline, 7 and 21-days of admission compared to the healthy controls (P < 0.0001). Similarly, miR-10b-5p expression levels were lower at the baseline and 21-days post admission (P = 0.001). The opposite situation was observed in miR-124-3p and miR-302c-5p. Enrichment analysis showed influence of analyzed miRNAs on IL-2 signaling pathway and multiple cardiovascular diseases through COVID-19-related targets. Moreover, the COVID-19-related genes regulated by miR-200b-3p were linked to T cell protein tyrosine phosphatase and the HIF-1 transcriptional activity in hypoxia. Analysis focused on COVID-19 associated genes showed that all analyzed miRNAs are strongly affecting disease pathways related to CVDs which could be explained by their strong interaction with the ACE2 network.

Table 1.Participant's characteristics.The data is presented in two ways, encompassing both the mean ± standard deviation (SD) and the median with the interquartile range, which provide insights into the characteristics of the data distribution.The P-value is calculated by employing a One-way ANOVA test for the purpose of comparing three groups.Abbreviations used include: AF for Atrial fibrillation, CAD for coronary artery disease, MI for myocardial infarction, BMI for body mass index, COPD for chronic obstructive pulmonary disease, CRP for C-reactive protein, CVD for cardiovascular disease, DM for diabetes mellitus, PCT for procalcitonin, SD for standard deviation, and WBC for white blood cells.Significant values are in bold.www.nature.com/scientificreports/post-admission when compared to those of healthy controls (p < 0.0001 at all-time points).Similarly, miR-10b-5p displayed reduced expression levels in COVID-19 patients compared to healthy individuals at both the baseline and 21 days after admission (p = 0.001 at both time points).Moreover, the expression of this miRNA demonstrated a notable increase 7 days after admission in comparison to the baseline, for which p-value was 0.003.Additionally, there was an observed downward trend in the expression levels of miR-200b-3p among COVID-19 patients.Specifically, it was significantly lower at 7 and 21 days after admission when compared to the baseline expression (p < 0.0001 and p = 0.003, respectively) and when contrasted with healthy individuals (p < 0.001 at both time points).The opposite situation was observed in miR-124-3p and miR-302c-5p concentrations.MiR-124-3p expression was significantly higher at the baseline, 7-days and 21-days post-admission when compared to the control group (P < 0.001, P < 0.001 and P = 0.017 respectively).In miR-302c-5p the significant differences were shown in the decrease of its expression in patients 21-days post-admission when compared to the control group (P < 0.0001) and patients 7-days after admission (P = 0.010).We found that MiR-200b-3p expression was significantly and negatively correlated with the hospital length of stay and/or death.Patients who were hospitalized over 21 days and/or died had significantly lower expression levels of miR-200b-3p when compared to those hospitalized under 21 days (P = 0.002).The ROC curve analysis demonstrated that reduced delta miR-200b-3p expression has predictive value in evaluating both the length of hospital stay and the likelihood of mortality, with an AUC of 0.730 and a p-value of 0.002 (Fig. 3b).

Participant's characteristics
The COVID-19 patient group was stratified into two subgroups using ROC curve analysis, considering changes in miR-200b-3p expression (delta), specifically, distinguishing between low and high values.This stratification was based on factors such as extended hospitalization or the occurrence of mortality (Fig. 3a).
The threshold of -0.522 was set to define a low delta miR-200b-3p expression level, encompassing 44% of the population.The area under the curve (95% CI) was 0.730 (0.61-0.86), with a positive predictive value (PPV) of 60% and a negative predictive value (NPV) of 79%.The sensitivity and specificity for this threshold were 69% and 72%, respectively (Supplemental Table 1).

Results of machine learning based analysis
In order to identify features strongly associated with COVID19 death and long hospitalisation we performed SHAP analysis which calculated the contribution of each feature to the prediction indicating their importance.The SHAP values we acquired enabled us to understand how having a particular value for a specific feature influences the prediction, as opposed to what the prediction would be with a baseline value for that feature.This method helps with the interpretability problem and can be used to improve the classification models (Fig. 4A  and B).From the Laboratory section miR-200b-3p on day 7th (low expression), miR-302c-5p on day 7th (low and moderate expression), and CRP on day 7th (high level), neutrophils on day 7th (high and average), D-Dimer on day 0 (very low and very high) appeared in the top important features associated with death + 21daygstay outcome according to SHAP analysis (Fig. 4B).Actually, considering all the executions, miRNA's miR-200b-3p, miR-302c-5p and miR10b-5p were always present among the top features, which reinforces their importances in predicting the death + 21 day stay outcome.Moreover, while most miR's have high values associated with not death, miR10b-5p on day 0 and miR-10 on day 7th specifically tend to present high values associated with death and long term hospitalization.Additional hierarchical analysis of SHAP values allowed us to identify associations between most important features.We observed redundancy between levels of D-dimer (day 7th) and miR-26-5p (day 0).In summary this machine learning based analysis identified miR-200b-3p levels on the www.nature.com/scientificreports/7th day, miR-302c-5p levels on the 7th day, CRP levels on the 7th day, neutrophil levels on the day of admission, and D-Dimer levels on the day of admission as the most reliable predictors of extended hospitalization in COVID-19 patients.

Target predictions of analyzed miRNAs in context of clinical data
To clarify the role of analyzed miRNAs in COVID19 hospitalization we performed target predictions in the context of several lists of genes associated with important features identified in clinical data, and processes recognized as playing a role in COVID19 disease, especially related to thrombosis.Prioritization based on the number of regulated targets within each of the analyzed gene lists highlighted miR-124-3p and miR-26b-5p as having an especially high number of targets of our interest.All analyzed miRNAs, except for miR-124-3p were identified as targeting ACE2 directly, besides other components of the ACE2 network.Three miRNAs, miR-124-3p, miR-26b-5p and miR-302c-5p were found targeting D-dimer related genes.The summary results of this analysis are shown in Fig. 4C.A complete list of genes from each list targeted by analyzed miRNAs is available in Supplementary File 2.

miRNA targets enrichment analysis
Enrichment analysis is a computational method allowing us to identify which biological terms affected are by analysed miRNAs, based on all genes they are targeting.We screened multiple databases provided by EnrichR tool, in order to identify top signalling pathways (Bioplanet_2019), and phenotypes (GWAS_Catalog_2019), affected by ACE2-related miRNAs.Pathway enrichment analysis showed distinct regulation between analysed miRNAs (Supplemental Fig. 1A).Among top shared pathways was BDNF signaling including miR-10b-5p, miR-124-3p, miR-200b-3p.Interleukin-2 signaling pathway, TGF-beta regulation of extracellular matrix and Pathways in cancer were overrepresented for targets of all analysed miRNAs.miR-200b-3p showed regulation of targets associated with developmental pathways, including signaling events mediated by HGFR, Insulin signaling, BDNF pathway and EGF/ EGFR signaling.
GWAS enrichment analysis for targets of analysed miRNAs, showed multiple significantly enriched terms related to CVD but only for miR-200b-3p and miR-302c-5p (Supplemental Fig. 1B).miR-200b-3p showed association with CVD-related sodium levels, systolic blood pressure and atrial fibrillation.While miR-302c-5p targeted genes associated with systolic blood pressure and QRS duration.
In order to evaluate which signalling pathways and disease-related pathways can be directly affected during COVID19 by our miRNAs of interest, we first subset COVID19 related targets of those miRNAs.Those targets were mentioned in at least 2 publications according to the DisGeNet database https:// www.disge net.org/ covid/ disea ses/ summa ry/ (837 genes in total).Further on selected targets we performed enrichment analysis using Bioplanet 2019 database and Elsevier pathway collection.Signalling pathway analysis pointed out again IL-2 signalling, but didn't include BDNF-related signaling (Fig. 5A).Other especially interesting affected pathways included Complement and coagulation cascades for miR-10b-5p; oncostatin M for miR-124-3p; signaling events mediated by T cell protein tyrosine phosphatase (TC-PTP) and HIF-1 transcriptional activity in hypoxia for miR-200b-3p; immune system for miR-26b-5p; and mTOR signaling for miR-302c-5p.Disease related ontological analysis showed that all analysed miRNAs are strongly targeting terms associated with CVDs (Fig. 5B).The most significantly enriched were Proteins Involved in Myocardial Ischemia (miR-200b-3p, miR-124-3p, miR-10b-5p), Proteins Involved in Myocarditis (miR-26b-5p), and Proteins Involved in Arterial Hypertension (miR-302c-5p).Additionally, for miR-200b-3p we observed enrichment of Proteins Involved in Prostate Cancer.

Discussion
Our previous bioinformatic analysis revealed the ability of miR-200b-3p to regulate ACE2 interaction network and differentially expressed genes in heart and Human-induced pluripotent stem cell-derived cardiomyocytes (hiPSC-CMs) infected with SARS-COV-2 1 .Another in silico analysis has shown miR-200b targeting host cell receptor ACE2 and was predicted to be highly expressed in respiratory epithelial cells, which indicates its role in COVID-19 infection 17 .Moreover, lower miR-200 expression is correlated with a lower survival rate in severe and critical clinical conditions 18 .Corresponding with that, in our study group, miR-200b levels were significantly lower in the COVID-19 patient group.Interestingly, the miR-200b-3p transfection alleviated the pulmonary fibrosis resulting in regeneration and slowdown of senescence through restoration of damaged alveolar cells 19 .Our enrichment analysis aligned with those results showing that an important target of this miRNA, are HGFR, EGFR and BDNF pathways associated with fibrosis.Analysis focused only on COVID-19-related targets identified signalling events mediated by TC-PTP and HIF-1 transcriptional activity in hypoxia as regulated by miR-200b-3p.Protein tyrosine phosphatases PTP-1B and TC-PTP are known to play nonredundant roles in macrophage development and IFN-γ signaling and JAK/STAT pathway 20 .Also HIF-1α is a key activator for both SARS-CoV-2 infection and inflammatory response and plays a role in aggravation of COVID19 21 .Those results strongly support the importance of miR-200b-3p as a biomarker and an important player in COVID19 progression.
In previous studies, the role of miR-26b was evaluated only by in silico analysis in the context of SARS-CoV2 infection.All three bioinformatic analyses showed miR-26b as a promising monitoring biomarker in COVID-19 and predicted the ACE2 as a strong miR-26b target 22 .Virus S-protein, ACE along with histone deacetylase (HDAC) pathway which are known to have a core role in SARS-CoV2 cellular entry and pathogenicity were mentioned to modulate miR-26b-5p expression levels.Further evaluation of top pathways targeted by miR-26b-5p occurred to be also highly linked and modulated by HDAC.This indicates a particular cross-talk between miR-26b-5p and HDAC that can be used as a therapeutic strategy 23 .Additionally, our previous publication emphasized miR-26b-5p potential role in COVID-19 associated heart failure (HF) progression 1 .In this current study for the first time, we confirmed the bioinformatic analysis in the samples from COVID-19 patients which showed a significant downregulation trend over time.Utilization of Machine learning related SHAP tool and hierarchical clustering allowed us to identify association of hsa-miR-26b-5p level on the first day of hospitalization with the later level of D-dimer (7th day of hospitalization).Target prediction analysis identified that this miRNA targets F3 which according to GWAS studies is associated with D-dimer levels 24 .D-dimer value on admission is a well-recognized and accurate biomarker for predicting mortality in patients with COVID-19 25 .In our study enrichment analysis of COVID-19 related targets of miR-26b pointed out strong regulation of proteins involved in myocarditis.Which is also related to SARS-CoV-2 infection 26 .Literature data shows multiple associations of miR-26b with MI and other CVDs, but so far there are no direct links with myocarditis 27,28 .
MiR-124-3p influences a variety of processes including neural functional recovery, adipogenic differentiation and wound healing.In the central nervous system (CNS) miR-124-3p was shown to promote functional repair in animal models of neurodegenerative diseases 29 .However, the role of miR-124 in SARS-CoV2 infection is highly limited.In the human study, the upregulation of miR-1-3p/miR-124-3p ratio was observed among patients who died due to SARS-CoV2 infection, but considering the number of days under invasive mechanical ventilation, the correlation was inverted 30 .Importantly, the latest study aimed to analyze the difference of miR-124 expression in the 5 different groups of COVID-19 patients and showed that the expression of miR-124 in the critically ill group was the highest compared to all the other groups 31 .Our results are in line with this observation as we have found significantly higher miR-124 expressions in the patients in three different time points compared to controls.Our enrichment analysis pointed out myocardial ischemia for its COVID-19 related targets.Literature data are consistent with this finding showing that miR-124 targets various pathways including SIRT1 and STAT3 which promotes ischemia-reperfusion induced cardiomyocyte apoptosis 32 .Additionally, enrichment analysis pointed out oncostatin M as the 1st rank signalling pathway for miR-124-3p.Oncostatin M is a member of the IL-6 family and stands as a mediator in various CVDs such as myocardial infarction or myocarditis and regulates cardiomyocyte remodeling 33 .Additionally for miR-124-3p, miR-302c-5p and miR-26b-3p we observed a high number of targets associated with iron overload, suggesting inhibition of this process.This finding is interesting taking into account that iron level increases risk of sepsis and severe COVID-19 34 .
Similarly, to the results obtained from our study, the miR-10b expression levels were found significantly lower in COVID-19 patients compared to controls.Additionally, the miR-10b expression significantly negatively correlated with the IL-2 and IL-8 serum levels as well as with the age of patients, erythrocyte sedimentation rate and CRP level, suggesting its role in viral infection-induced enhanced cytokine storm 35 .Our results confirmed this relationship showing that miR-10b-5p targets genes significantly associated with IL-2 signalling.An even stronger association was observed for IL-1 regulation of the extracellular matrix.This result points out its potential role of miR-10b-5p the fibrotic process and lung repair 36 .Our enrichment analysis showed that miR-10b-5p strongly regulates COVID-19 targets associated with Myocardial ischemia.Literature data shows that overexpression of miR-10b-5p in the murine model of MI significantly reduced MI size, improved cardiac function, and inhibited apoptosis 37 .miR-10b-5p was identified by us as a key regulator of ACE2 network which aligns with COVID-19 cardiovascular effects.
MiR-302c-5p role is mainly investigated in viral infections, nervous system diseases and cancers 38,39 .The miR-302 has not yet been described specifically in the context of the SARS-CoV2 infection.Lower expression of miR-302a and increased levels of NF-κB direct target, interferon regulatory factor-5 (IRF-5) was found in patients with influenza A infection compared to controls.IRF-5 promotes viral replication, the opposite effect was shown by miR-302a treatment, which inhibits viral replication 40 .This can be the explanation of our observation as in our study we found that miR-302c expression decreased over time in patients with COVID-19.Moreover, our Vol:.(1234567890 www.nature.com/scientificreports/previous bioinformatic analysis of the ACE2 interaction network showed that miR-302c-5p modulated a high number of virus-linked proteins 1 .Therefore, we suggest that miR-302c-5p might be a prognostic biomarker for COVID-19 patients.Additionally, members of miRNA-302 family members are considered potential biomarkers for the diagnosis of acute heart failure 41 .This aligns with our results regarding the CVDs related ontological terms associated with targets of miR-302c-5p, especially Arterial hypertension. The enrichment analysis performed in this study focused on COVID-19 related targets of analyzed miRNAs showed that all of them strongly regulated CVDs-related phenotypes.This result is consistent with our in silico findings regarding the role of ACE-2 interaction network in Cardiovascular outcomes in COVID19.That strong regulation of the ACE-2 network can have an effect in developing, or even aggravating, the CVDs phenotypes 42,43 .But also we observed that analyzed miRNAs regulated COVID-19 related genes associated with other phenotypes like Chronic Obstructive Pulmonary Disease (COPD) (miR-10b-5p), Glomerulonephritis (miR-124-3p), Prostate Cancer (miR-200b-3p), Amyotrophic Lateral Sclerosis (miR-26b-5p), and Endometriosis (miR-302c-5p) which could be important factors for consideration of post-COVID19 complications.
We identified the following top targeted pathways affected by analyzed miRNAs: IL-2 signalling pathways, signaling events mediated by T cell protein tyrosine phosphatase (TC-PTP) and HIF-1 transcriptional activity in hypoxia which would require further investigation.Regarding IL-2, certain studies propose that there is an elevation in the IL-2 cytokine family levels during the infection, leading to a heightened inflammatory response and the potential occurrence of a cytokine storm 44 .Elevated cytokine levels worsen prognosis in immunocompromised patients, while in those without underlying conditions, IL-2 family upregulation correlates with disease severity.Conflicting studies exist, with some showing no significant cytokine changes.Healthy individuals postvaccination exhibit IL-2 upregulation, but immunocompromised patients require booster doses for efficacy.The IL-2 cytokine family holds potential as immunotherapy for COVID-19 thus it seem necessary to validate those findings on protein level.
TC-PTP is an enzyme involved in the regulation of immune responses and cellular signaling, especially in JAK-STAT pathway, important in COVID19 progression 45 .Crucial in this process ORF3a induces renal tubule injuries by activating NF-κB and STAT3 signaling.It promotes STAT3 activation by interacting with the ubiquitin E3 ligase TRIM59.This interaction leads to the dissociation of the phosphatase TC-PTP from STAT3, inhibiting STAT3 dephosphorylation and resulting in prolonged STAT3 activation 46 .Investigation of the role of TC-PTP signalling seems important for better understanding the COVID progression.
As mentioned before, HIF-1α is a key activator for both SARS-CoV-2 infection and inflammatory response and plays a role in aggravation of COVID19 21 .There are significant roles of SARS-CoV-2 ORF3a and HIF-1α in virus infection and pro-inflammatory responses.In one study RNA sequencing reveals dysregulation in HIF-1α signaling, immune response, and metabolism pathways in COVID-19 patients 21 .Clinical analyses highlight increased HIF-1α production, inflammatory responses, and mortality rates in elderly patients.HIF-1α TC-PTP and HIF-1α ORF3a related induces mitochondrial damage and Mito-ROS production, promoting HIF-1α expression, facilitating SARS-CoV-2 infection and cytokine production.Importantly, HIF-1α broadly promotes infection in various viruses.Overall, ORF3a-induced HIF-1α aggravates viral infection and inflammatory responses during SARS-CoV-2 infection, establishing HIF-1α's pivotal role in promoting SARS-CoV-2 infection and inducing pro-inflammatory responses in COVID-19 21,47 .Thus it seem important to evaluate the better the role of this factor in COVID19 progression.
In this study we identified and validated miRNAs which could serve as novel, predictive biomarkers and their targets of the COVID-19 long term hospitalisation, and can be used for early stratification of patients and prediction of severity of infection development in an individual.We found that low delta miR-200b-3p expression presents predictive utility in assessment of the hospital length of stay and/or death.Analysis focused on COVID-19 associated genes showed that all analysed miRNAs are strongly affecting multiple signalling and disease pathways related to CVDs which could be explained by their strong interaction with the ACE2 network.Identifying novel biomarkers and implementing machine learning based tools enabling feature identification and bioinformatics may improve clinical data interpretation and prediction of the outcome in patients with COVID-19.

Study group
Seventy-nine patients diagnosed (positive nasopharyngeal swab PCR test) with COVID-19 were enrolled to the study during the third pandemic wave in Poland (January-May 2021).Patients were admitted to the Military Institute of Medicine in Warsaw 2 .Thirty-two age and sex and cardiovascular risk factors matched SARS-CoV-2 infection-free participants were sampled as a control group in out-patients clinic.Blood samples of COVID-19 patients were collected at three different time points, whereas control group patients were sampled only once 2 .
Medical records were obtained and demographic, clinical and laboratory data were presented in Table 1.Whole blood samples were collected by using Tempus™ Blood RNA Tube (Applied Biosystems) and were kept in -80 ˚C until the day of experiments in Medical University of Warsaw.No freeze-thaw cycles were performed during the experiments.

RNA preparation, detection, and quantification of miRNAs by Quantitative PCR
We isolated total RNA using Tempus™ Spin RNA Isolation Kit (Invitrogen).Next, we performed a reverse transcription reaction using the TaqMan Advanced miRNA cDNA synthesis kit (ABI, California, USA) following the guidelines provided by the manufacturer.Afterwards, miRNA expressions were detected by quantitative polymerase chain reaction (qPCR) using TaqMan miRNA Assay kits (ABI, California, USA, catalog number A25576, assay ID: 478418_mir; 480901_mir; 478494_mir; 477963_mir; 478800_mir) by using CFX384 Touch Real-Time PCR Detection System (BioRad Inc.Hercules, California, USA).Cel-miR-39 was spiked-in during miRNA extraction phase and was used as an external normalizer.All reactions were performed in triplicates, and mean values were used in statistical analysis as described before 48,49 .MiRNA expressions were expressed as 2-ΔΔCT 50 , results were log-10 transformed for statistical analysis.

Statistical analysis
Categorical variables were presented as a number and percentage, continuous variables were expressed as mean ± standard deviation (SD) or median and interquartile range (IQR) based on the distribution of the data.They were selected based on the normality of distribution evaluated by Shapiro-Wilk test.Depending on the normality of the distribution, the Student's t-test or Mann-Whitney test was used for unpaired samples and the Wilcoxon test for paired samples.One-way ANOVA or Kruskal-Wallis tests were used to compare more than two groups, accordingly.To assess the predictive value of miR-200b-3p changing over time (delta miR-200b-5p expression) for increased hospital length of stay or death in follow-up as a composite endpoint we used receiver operating characteristic (ROC) analysis.Low delta miR-200b-3p expression (< -0.522 expression), age (years), male sex, BMI > 30, hypertension, diabetes, smoking, and coronary artery disease were included in the multivariate logistic regression analysis model.All tests were two-sided with a significance level of p < 0.05.Calculations were performed using SPSS version 22.0 (IBM Corporation, Chicago, USA).Graphs were improved by Adobe Illustrator 24.0.2.

SHAP analysis of the clinical data
We applied SHAP analysis to the cleaned clinical dataset considering the 'Demographics' and 'Laboratory Results' Section.For this analysis, we considered as output of interest (target) the column 'death + 21day_stay' , which is a binary feature referring to patients who died or not after 21 days after the hospital admission.Then, we ran the XGBoost classifier ten times to check robustness, which led to similar results with low variation in the top 10 features.

Gene lists selection
We downloaded the following lists of genes related to our processes and phenotypes of interest.From the Gene Ontology (GO) database we obtained gene lists associated with coagulation (17 GO terms) platelet activity (26 GO terms), and inflammatory response (22 GO terms) using the biomartR R package 51 (Supplementary file 2, spreadsheet GO terms).From the Malacards database we obtained gene lists related to thrombosis (https:// www.malac ards.org/; ID:THR024).From the Disgenet database we obtained multiple gene lists: related to thrombosis (C0040053), Diabetes Mellitus Insulin-Dependent (genes with at least 2 NCBI publications; C0011854), COVID19 (genes with at least 2 NCBI publications; https:// www.disge net.org/ covid/ genes/ summa ry/), White blood.cellabnormality (C0152009), C reactive protein.measurement(C0201657), Iron Overload (C0282193), Iron deficiency (C0240066), Lymphocyte Count measurement (C0200635) related to clinical data.List of ACE2 related genes was obtained from previous publication regarding ACE2 network in COVID19 1 .D-dimer related gene lists were obtained from literature data 52 .All gene lists are available in Supplementary file 2 (spreadsheet gene lists combined).All gene symbols were unified using the NCBI annotation file.

miRNA target prediction and ranking
We used on all steps of bioinformatic analyses we used "wizbionet" R-package https:// github.com/ wizbi onet/ wizbi onet/ blob/ master/ doc/ vigne tte_ wizbi onet.md 53 .To identify targets of analyzed miRNAs we used the multiMiR R package 54 .We searched the top 20% hits among all conserved and non-conserved target sites in 14 target prediction databases.To evaluate how many targets from each gene list of interest were regulated by analysed miRNAs we used wizbionet R package 53 .

Enrichment analysis
Enrichment analysis was done using the API for EnrichR database by implementing Hypergeometric test with Benjamini and Hochberg correction, while the reference was the human genome.For all statistical analyses, the significance cutoff was set to adjusted p-value ≤ 0.05.We performed enrichment analyses using the following databases: ARCHS4_Tissues, BioPlanet_2019, Elsevier_Pathway_collections, and GWAS_Catalog_2019.Data visualization of the overlap between top five most significant pathways for each analysed miR was performed in R using ggplot2 and ggrepel libraries.

Figure 1 .
Figure 1.Top predicted miRNA modulators of the ACE2 network in COVID-19 and workflow leading to their identification.Statistics for top miRNAs regulating the highest number of genes within ACE2 related-networks based on previous in silico predictions 1 .Red squares show if miRNA was present among the top miRNAs in a analyzed dataset.MiRNAs validated in present study using qRT-PCR have blue labels.

Figure 3 .
Figure 3. Delta of miR-200b-3p (7 days after admission compared to the day of admission) (a) box-plots comparing increased hospital length of stay and/or mortality, and (b) ROC curves for predicting hospital length of stay and/or death.Abbreviations used: AUC for Area under the ROC Curve, COVID-19 for coronavirus disease 2019, and miR for microRNA.

Figure 4 .Figure 5 .
Figure 4. Results of SHAP method applied by using XGBoost classifier to the panel (A) demographic and panel (B) laboratory part of COVID19 data.Results show the contribution locally per sample of each feature to the prediction and the most important features associated with the outcome death + 21 day stay.The features with highest importance for a given outcome are localized on the top of the graph.We highlight the features BMI (from very high and very low spectrum), age (high) and diabetes (positive) from the Demographics section are associated with death + 21 day stay outcome (Fig.3a).Panel (C) Numbers of targets of interest regulated by analysed miRNAs.Analysis was performed separately for each miRNA and then combined.Genes from specific gene lists were divided into four clusters based on their amount using the OneR algorithm.

Table 2 .
A multivariable logistic regression model was employed to predict extended hospitalization and/or mortality, incorporating the delta of miR-200b-3p alongside clinical variables.Significant values are in bold.OR odd ratio, CI confidence interval, miR microRNA, BMI body mass index, CVD cardiovascular disease.