SnoRNAs and miRNAs Networks Underlying COVID-19 Disease Severity

There is a lack of predictive markers for early and rapid identification of disease progression in COVID-19 patients. Our study aims at identifying microRNAs (miRNAs)/small nucleolar RNAs (snoRNAs) as potential biomarkers of COVID-19 severity. Using differential expression analysis of microarray data (n = 29), we identified hsa-miR-1246, ACA40, hsa-miR-4532, hsa-miR-145-5p, and ACA18 as the top five differentially expressed transcripts in severe versus asymptomatic, and ACA40, hsa-miR-3609, ENSG00000212378 (SNORD78), hsa-miR-1231, hsa-miR-885-3p as the most significant five in severe versus mild cases. Moreover, we found that white blood cell (WBC) count, absolute neutrophil count (ANC), neutrophil (%), lymphocyte (%), red blood cell (RBC) count, hemoglobin, hematocrit, D-Dimer, and albumin are significantly correlated with the identified differentially expressed miRNAs and snoRNAs. We report a unique miRNA and snoRNA profile that is associated with a higher risk of severity in a cohort of SARS-CoV-2 infected patients. Altogether, we present a differential expression analysis of COVID-19-associated microRNA (miRNA)/small nucleolar RNA (snoRNA) signature, highlighting their importance in SARS-CoV-2 infection.


Introduction
At the beginning of 2020, the World Health Organization (WHO) declared coronavirus disease-2019 (COVID-19) as a global pandemic. The causative organism, severe acute Vaccines 2021, 9, 1056 2 of 16 respiratory syndrome coronavirus 2 (SARS-CoV-2), exhibits a wide spectrum of clinical manifestations in disease-ridden patients. Differences in the severity of COVID-19 ranges from asymptomatic infections and mild cases to the severe form, leading to acute respiratory distress syndrome (ARDS) and multiorgan failure with poor survival [1]. Moreover, the mortality rate is influenced by aging, viral strain, pre-existing comorbidities, and the degree of immunocompromise. Indeed, the health and socioeconomic implications of the COVID-19 pandemic are enormous, and thus warrants the search for new interventions and treatment measures. Recent research has suggested that a unique non-coding RNA signature can aid in identifying the likelihood of developing specific disease outcomes [2]. Alterations of miRNA levels in the blood have been described in multiple inflammatory and infectious diseases, including SARS-related coronaviruses [3][4][5][6][7][8][9]. MiRNAs are endogenous small non-coding RNAs, around 22 nucleotides, that bind specific messenger RNAs (mRNAs) through complementary base-pairing [10]. Hence, miRNAs can regulate various cellular processes, including proliferation, apoptosis, and differentiation, by binding to the 3 UTR of target mRNAs inducing their degradation, thus serving a fundamental role in post-transcriptional repression [11,12]. In this context, a single miRNA can target several genes, and multiple miRNAs may regulate a single gene. Hence, identification of miRNAs as well as characterization of miRNA-mRNA interactions in SARS-CoV-2 infection is important to understand their role in disease pathogenesis, progression, and severity [13][14][15]. Accumulating evidence is also implicating snoRNAs in numerous physiological and pathological processes, including their interactions with some RNA viruses [16,17]. SnoRNAs have been shown to play a role in ribosomal RNAs (rRNAs) modification and maturation. These include pseudouridylation, 2 -O-methylation, polyadenylation, alternative splicing, and formation of protein complexes with fibrillarin that is a component of several ribonucleoproteins [18]. There are no reports yet with a comprehensive characterization of snoRNAs in SARS-CoV-2 infection. Here, we report for the first time the snoRNA signature, and similar to previous studies the miRNA signature in the peripheral blood of severe COVID-19 cases (n = 9), as compared to mild (n = 10) and asymptomatic (n = 10) patients. This study analyzed a total of 29 COVID-19 patients that were matched for age and comorbidities using Affymetrix GeneChip miRNA 4.0 array and were recruited between July and October 2020 in Qatar.

Study Design and Data Collection
This research is in accordance with the Reporting of Observational Studies in Epidemiology (STROBE) recommendations, the Code of Ethics of the World Medical Association. This study was granted ethical approval from the Medical Research Center at Hamad Medical Corporation (MRC-05-084, Immunological and immune-genetic investigations in COVID-19 patients with varying disease severity, 06/21/2020). All study participants gave written informed consent where possible, and deferred consent was obtained for ICU cases. From 173 recruited patients in this prospective cohort study, 29 male age-matched patients were included. All patients were previously diagnosed with COVID-19 using TaqPath COVID-19 Combo Kit (Thermo Fisher Scientific, Waltham, MA, USA), or Cobas SARS-CoV-2 Test (Roche Diagnostics, Rotkreuz, Switzerland), with a CT value < 30. Additional criterion for selection was age between 35 and 75 years. Participants were grouped into severe, mild, and asymptomatic. Classifying severe cases was based on the requirement of high-flow oxygen support and ICU admission (n = 9). Whereas mild patients were identified based on symptoms, and positive radiographic findings with pulmonary involvement (n = 10). Patients with no clinical presentation were labelled as asymptomatic cases (n = 10). In the severe group, three patients died with respiratory failure listed as the primary cause of death. Blood samples were collected in the PAXgene Blood RNA Tubes (PreAnalytiX) at the time of diagnosis, prior to isolation, or hospitalization. Routine laboratory tests that were performed include complete blood cell counts, electrolytes, glucose, albumin, total protein, C-reactive protein, procalcitonin, IL-6, D-dimers, ferritin, urea, and liver enzymes.

RNA Isolation and Quality Control
Peripheral venous blood (2.5 mL) in PaX gene tubes was inverted 8-10 times to ensure complete mixing with the lysis reagent. Tubes were then stored upright at room temperature for a minimum of 2 h, before transferring them to a freezer at −80 • C until RNA isolation. The RNA (both total and miRNA) was isolated with a blood miRNA kit from Qiagen (PreAnalytiX GmBH Hombrechtikon, Switzerland) following the manufacturer's instructions. The concentrations and purity of the RNA samples were evaluated spectrophotometrically (Nanodrop ND-1000, Thermo, Wilmington, DE, USA). The RNA isolation process was validated by analyzing the integrity of RNA with the RNA 6000 Nano Chip Kit (Agilent, CA, USA); the presence of the small RNA fraction was confirmed by the Agilent Small RNA Kit (Agilent).

Microarray and Data Analysis
Affymetrix GeneChip miRNA 4.0 array was used following the manufacturer's instructions for miRNA expression analysis. For each sample, 250 ng of RNA was labelled using the FlashTag™ Biotin RNA Labeling Kit (Genisphere, Hatfield, PA, USA). Following this, the labelled RNA was quantified, fractionated, and hybridized to the miRNA microarray with continuous agitation at 60 rpm for 16 h at 48 • C on a GeneChip Hybridization Oven 640. The miRNA microarray chips were then washed and stained using the GeneChip Fluidics Station 450 (Affymetrix, Santa Clara, CA, USA). Finally, the miRNA microarray chips were scanned using an Affymetrix GCS 3000 scanner (Affymetrix, Santa Clara, CA, USA), and the signal values were evaluated using the Affymetrix ® GeneChip™ Command Console software.
Raw data were extracted using the Affymetrix data extraction protocol in the Affymetrix GeneChip ® Command Console ® (AGCC) Software. All statistical analyses including raw data import, annotation, and quality control were conducted with R statistical software (version 4.0.4, R Foundation, Boston, MA, USA). R Bioconductor packages "oligo" [19], "oligoCalsses", and "BioBase" were used for the pre-processing of microarray datasets [20]. Raw expression data in CEL files were read and parsed into ExpressionFeatureSet using the "read.celfiles" function, followed by sample and probe annotations with "Annotated-DataFrame" function. Next background correction, normalization, and log 2 transformation using the robust multiarray average (RMA) method was performed. The differentially expressed transcripts were screened out via the "limma" package [21]. The analysis was performed using the design matrix "~severity" as a factor, where it was either "asymptomatic", "mild" or "severe". A comparative analysis between different severity groups was carried out using thresholds of |log 2 FC| > 1.5 and FDR < 0.05. Hierarchical clustering of differentially expressed miRNAs was carried out using the "pheatmap" R package. The prediction of miRNA-target genes was performed using the "multiMiR" R package filtering with validated databases, including miRecords, miRTarBase, and TarBase, and the list of 26 differentially expressed miRNAs in severe vs asymptomatic and severe vs mild comparisons [22]. The package "miRBaseConverter" was used to retrieve the miRNA sequence based on the accession number [23]. The database miRTarBase was used with a threshold of 2 for the minimum number of miRNA-target interactions, 1 for the adjusted p-value (FDR), and with filter by evidence categories set to strong [24]. The miRNA-target genes network was built using the "MIENTURNET" tool [25], while the final network was visualized using Cytoscape software [26]. Functional enrichment analysis of target genes list was performed using the "Enrichr" platform including Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis. Spearman rank correlation tests were calculated to assess the correlation between blood parameters and the levels of DE miRNAs and snoRNAs employing "ggstatsplot" R package. Other helper packages in plotting and formatting include "dplyr", "kableExtra", "gtsummary", "flextable", "officer", "grid", "VennDiagram", "pheatmap", "gplots", "ggplot2", "stringr", "reshape2", "corrplot", "tidyr", "tidyverse". Age was normally distributed and hence mean ± standard deviation was reported. Since all clinical variables had a non-normal distribution, comparisons between different groups of severity were performed using the Kruskal-Wallis rank sum test for continuous variables, and Fisher's Exact Test for categorical. Clinical values were reported as medians and interquartile range [IQR]. Two-tailed p-values were calculated, and p-value < 0.05 was considered statistically significant. The R code that was used to generate the figures and analysis is available as a supplementary file.

Results
First, we identified differentially expressed miRNAs and snoRNAs in severe patients as compared to mild and asymptomatic cases, a |log 2 (fold change)| > 1.5 and an FDR < 0.05 were used as cutoffs. We found 12 differentially expressed miRNAs and snoRNAs in severe versus asymptomatic, and 31 differentially regulated miRNAs and snoRNAs in severe versus mild cases. Of the 4 miRNAs and snoRNAs that were unique in severe versus asymptomatic, there were 3 upregulated and 1 downregulated transcripts, whereas, of the 23 uniquely differentially expressed miRNAs and snoRNAs in severe versus mild, 7 were upregulated and 16 downregulated. Interestingly, we also observed 3 miRNAs and 5 snoRNAs that were common in both severe versus mild, and severe versus asymptomatic comparisons ( Figure 1A). Altogether, differential expression (DE) analysis showed a signature of 27 unique miRNAs and snoRNAs that was associated with COVID-19 severity ( Figure 1B,C). Notably, the most highly differentially expressed miRNAs and snoRNAs in severe versus asymptomatic are hsa-miR-1246, ACA40, hsa-miR-4532, hsa-miR-145-5p, and ACA18. Furthermore, the top 5 differentially expressed miRNAs and snoRNAs in severe versus mild are ACA40, hsa-miR-3609, ENSG00000212378 (SNORD78), hsa-miR-1231, hsa-miR-885-3p (Tables S1 and S2). Additionally, we report for the first time 9 differentially expressed snoRNAs, namely ACA18 (SNORA18), ACA20 (SNORA20), ACA40 (SNORA40), ACA57 (SNORA57), ENSG00000212378 (SNORD78), U17b, U44, U78, and U79 (Table S3). Our novel findings suggest that SARS-CoV-2 infection benefits from gene expression programs regulated by snoRNAs, including C/D box snoRNAs (SNORDs) and H/ACA box snoRNAs. Interestingly, no differentially expressed miRNAs or snoRNAs were found in mild versus asymptomatic comparison with the applied cut-offs. These data suggest that both miRNAs and snoRNAs dysregulation is associated with a severe form of COVID-19.
Next, we identified putative targets for the identified differentially expressed miRNAs and snoRNAs in a multi-ethnic cohort of 29 COVID-19 patients residing in Qatar and presenting with varying disease severity. Then we used miRTarBase package validated databases (mirecords, mirtarbase, and tarbase) to identify putative pathways that may be implicated in severe cases of COVID-19. To this end, we used the miRTarBase database, which is curated with >50,000 miRNA-target interactions and validated experimentally by reporter assay, Western blot, microarray, and next-generation sequencing. In total, we found 2584 unique validated target genes for the differentially expressed miRNAs. The resulting target genes list was used in functional enrichment analyses to identify top deregulated pathways in severe COVID-19 patients. Remarkably, key enriched KEGG pathways for the 2584 miRNA-target genes included mainly cancer-related pathways such as pathways in cancer, hepatocellular carcinoma, chronic myeloid leukemia, microRNAs in cancer, breast cancer, gastric cancer, pancreatic cancer, colorectal cancer, proteoglycans in cancer, cellular senescence ( Figure 2B). Furthermore, other hallmark gene sets in the Molecular Signatures Database with important biological implications at the top of the list include TNF-alpha Signaling via NF-kB, UV Response Dn, IL-2/STAT5 Signaling, Hypoxia, Inflammatory Response, Apoptosis, PI3K/AKT/mTOR Signaling, G2-M Checkpoint (Table S4) [27]. Interestingly, miRNA-target gene network analysis (MTGN) elucidated 7 nodes of differentially expressed miRNAs namely hsa-miR-145-5p, hsa-miR-199a-5p, hsa-miR-98-5p, hsa-miR-139-5p, hsa-let-7i-5p, hsa-miR-1246, and hsa-miR-572 (Figure 2A,C). We demonstrate 4 regulated genes targeted by more than 6 differentially expressed miRNAs (Table S5).    The prospect of using blood miRNAs and snoRNAs as biomarkers can be instrumental in identifying patients with a higher risk of severity and mortality. In this context, we calculated a correlation matrix to examine individual associations between the identified common ( Figure 3) and unique ( Figure S1) differentially expressed miRNAs and snoRNAs in severe cases, and routine laboratory tests such as complete blood count (CBC), glucose, electrolytes, liver, and kidney function parameters and inflammatory markers. Several of these clinical markers have been linked to prognosis in COVID-19 patients. We found a significant correlation between several differentially expressed miRNAs and snoRNAs and certain hematological and serological parameters, including WBC count ( Figure 4A), Vaccines 2021, 9, 1056 9 of 16 lymphocyte (%) ( Figure 4B), ANC, hematocrit, and albumin ( Figure S2). Table 1 represents a summary of patient's demographics and clinical parameters used in the correlation analysis. These observations can be validated in larger cohorts to help in the early detection of patients at higher risk of severity, and thus in selection of dynamic treatment regimens.   orange), and severe (n = 9, blue) cases. The comparison between different severity groups was done using the non-parametric Kruskal-Wallis one-way ANOVA test, the median for asymptomatic, mild, and severe groups is reported as a measure of centrality in each scatter plot. In each respective comparison, the p-adjusted values were calculated using Holm's sequential Bonferroni procedure for multiple hypothesis tests at an alpha level of 0.05. (B) A correlation matrix of the 8 miRNAs and snoRNAs with associated clinical markers. Each cell contains a correlation coefficient between the possible pairs of variables, namely Spearman's rho statistic calculated at a significance level of 0.05. Color scale ranges from blue (r = −1) to white (r = 0) to red (r = 1).

Discussion
Numerous differentially expressed miRNAs that are observed in this study have been highlighted in host-pathogen interactions [28]. For instance, the family of let-7e/miR-125a/miR-200 miRNAs have been reported to mediate ACE2 gene silencing [29]. Likewise, KEGG enrichment analysis of the differentially expressed miRNAs targets has shown pathways involved in cellular proliferation, invasion, and apoptosis such as PI3K/AKT/mTOR signaling. Whereas, inhibition of hsa-miR-1246 was demonstrated to reduce the cytotoxicity of the Ebola virus glycoprotein in vitro. With multiple viral miRNAs sharing similarities with the host miRNAs, in silico computational studies also uncovered several putative host miRNAs involved in controlling viral replication and limiting disease progression [14,[30][31][32]. These host miRNAs are predicted to bind viral sequences modulating cellular immune responses, metabolic pathways, and inhibiting host miRNA maturation, thus facilitating viral escape. Notably, we found 10 differentially expressed miRNAs targeting CDKN1A, a potent cyclin-dependent kinase inhibitor, also known as p21. Interestingly ivermectin inhibits p21-activated kinase 1 (PAK1), a serine/threonine kinase with oncogenic activity. Caly et al. reported that ivermectin decreased SARS-CoV-2 RNA viral load in vitro by 5000-fold with a single treatment [33]. Another target of multiple differentially expressed miRNAs is a suppressor of cytokine signaling 7 (SOCS7), which is known to exert an inhibitory effect on interferon responses, thereby facilitating viral replication. Moreover, insulin-like growth factor I receptor (IGF1R) was also regulated by 6 different differentially expressed miRNAs detected in severely ill patients. IGF1R has been reported as a critical regulator of transformation events, with high levels of expression in most malignancies, where it acts as an anti-apoptotic agent by enhancing cell survival [34]. The MTGN analysis suggests that SARS-CoV-2 infection is regulated by complicated miRNA regulatory networks, through multiple miRNAs targeting the same gene, and single miRNAs targeting multiple genes.
Several studies have reported the fundamental role of interferons in COVID-19 disease severity [35]. Similarly, the identified miRNA in this study ( Figure 2) suppresses interferon signaling via binding interferon target genes. Additionally, multiple viruses, including SARS-CoV-2 have been reported to enhance TGF-β signaling, which is known to induce fibrosis and suppress adaptive immunity [36]. Our data suggests a modulation of TGF-β signaling, via the surface receptors and canonical SMAD and MAPK pathways, regulating adaptive immune responses and tissue repair. These findings are in line with the relative lymphopenia reported in severe COVID-19 [1]. Captivatingly, siRNA studies targeting candidate snoRNAs provide an evidence of their functional roles in virus-host interactions against numerous viruses, while knock-down studies have demonstrated that RNA viruses require specific C/D box snoRNAs for optimal replication [17]. Thus, a validated COVID-19 miRNA and snoRNAs signature could be a useful tool to discriminate COVID-19 infections from other respiratory viral infections, identify asymptomatic infections, and develop proteins and metabolite-based tests. This report provides a systematic identification of miRNAs and snoRNAs profile in the blood of SARS-CoV-2 patients with different disease severity, which expands on the previous computational approaches in COVID-19-associated miRNAs and snoRNAs profiling. The patient samples were obtained at a single time point ranging between 2-5 days post COVID-19 diagnosis in Qatar. The patients' clinical history and lab investigations are well documented in the national healthcare database (CERNER), based on which none of the studied patients were previously diagnosed with any chronic illness, cancer, or immunological disorder, except 18 patients who had type 2 diabetes mellitus. Outcomes from our study can guide future studies involving longitudinal sampling and analysis of circulatory miRNAs/snoRNAs that play a role in disease resolution versus those involved in the development of long COVID. Moreover, future studies with a larger sample size are needed to delineate miRNAs and snoRNAs profiling in females, as well as in children. This study has extended our understanding of the miRNAs and snoRNAs regulatory mechanisms underlying the pathogenesis of SARS-CoV-2 infection, which can be explored further to produce promising therapies.  Figure S1: Spearman correlation matrix of clinical markers with uniquely DEMIs and snoRNAs in (A) severe versus mild and (B) severe versus asymptomatic comparisons. Each cell contains a correlation coefficient between the possible pairs of variables, namely Spearman's rho statistic calculated at a significance level of 0.05. Color scale ranges from blue (r = −1) to white (r = 0) to red (r = 1), Figure S2: Correlation between expression levels of common differentially expressed miRNAs and snoRNAs (in both severe versus asymptomatic, and severe versus mild comparisons) with (A) hematocrit, (B) absolute neutrophil count (ANC), and (C) albumin. Scatter plots representation with x-axis as Log 2 normalized transcript expression level and y-axis as the clinical variable measurements. Each dot represents a single miRNA/ snoRNAs transcript. Dots that are tagged with D letter represent a deceased patient measurement. For all statistical tests in the plots, the APA standard for statistical reporting is shown for Spearman correlation test including evidence in favor of null over alternative hypothesis, natural logarithm of Bayes Factor, p-value, Spearman's rank correlation coefficient, confidence intervals, and number of observations, Table S1: Differentially expressed miRNAs and snoRNAs in severe versus asymptomatic and severe versus mild COVID-19 patients. Probe ID and Name, Transcript ID and accession, Log 2 fold change, and FDR values are shown, Table S2: Mature sequence of annotated differentially expressed miRNAs that are available in miRbase. microRNA names, miRBase accession IDs, and microRNA mature sequences are shown. For hsa-mir-453, reads that map to the annotated mir-4532 locus (many with one mismatch) map exactly to annotated 28S rRNA sequences. The miRNA annotation is therefore likely to be false, and the miRNA was therefore removed from the database, Table S3: Sequences of differentially expressed small nucleolar RNAs in severe cases of COVID-19. Transcripts include snoRNA predicted using sequences from RFAM and miRbase. The following snoRNAs transcript IDs are not annotated in miRbase database, Table S4: Gene Set Enrichment Analysis (GSEA) of differentially expressed miRNAs-targets via Enricher platform and Molecular Signatures Database (MSigDB). Each hallmark gene set is an expressed signature derived by aggregating many MSigDB gene sets to characterize reported biological states or processes, Table S5: Candidate genes that connect with more than 6 miRNA-target interactions using miRTarBase. Informed Consent Statement: All study participants gave written informed consent where possible, and deferred consent was obtained for ICU cases.