Interplay between SARS‐CoV‐2‐derived miRNAs, immune system, vitamin D pathway and respiratory system

Abstract The new coronavirus pandemic started in China in 2019. The intensity of the disease can range from mild to severe, leading to death in many cases. Despite extensive research in this area, the exact molecular nature of virus is not fully recognized; however, according to pieces of evidence, one of the mechanisms of virus pathogenesis is through the function of viral miRNAs. So, we hypothesized that SARS‐CoV‐2 pathogenesis may be due to targeting important genes in the host with its miRNAs, which involved in the respiratory system, immune pathways and vitamin D pathways, thus possibly contributing to disease progression and virus survival. Potential miRNA precursors and mature miRNA were predicted and confirmed based on the virus genome. The next step was to predict and identify their target genes and perform functional enrichment analysis to recognize the biological processes connected with these genes in the three pathways mentioned above through several comprehensive databases. Finally, cis‐acting regulatory elements in 5′ regulatory regions were analysed, and the analysis of available RNAseq data determined the expression level of genes. We revealed that thirty‐nine mature miRNAs could theoretically derive from the SARS‐CoV‐2 genome. Functional enrichment analysis elucidated three highlighted pathways involved in SARS‐CoV‐2 pathogenesis: vitamin D, immune system and respiratory system. Our finding highlighted genes' involvement in three crucial molecular pathways and may help develop new therapeutic targets related to SARS‐CoV‐2.

SARS-CoV-2 (also known as coronavirus disease 2019 ) is a positive-stranded enveloped RNA virus. It has a +ssRNA genome with nearly 30 kb of size (ranging from 26 000 to 37 000 bases) in length. Genomic structure of SARS-CoV-2 comprises of a 5′-leader-untranslated region(UTR)-replicase-S (Spike)-E (Envelope)-M (Membrane)-N (Nucleocapsid)-3′UTRpoly (A) tail. 2 Each part of SARS-CoV-2 has a unique role in viral pathogenesis, and isolation of SARS-CoV-2 genomic sequence has demonstrated 88% of identity with two bat-derived SARS-like coronaviruses. 3 Immune system is the first and currently the only defence against the virus. Once, the virus enters the body and is detected by the angiotensin-converting enzyme 2(ACE2) factor, the virus enters the host cell and inserts its RNA genome into cytoplasm of that cell.
There, it multiplies and a large number of new viral particles germinate out of the host cell and move to other cells. In the first stage, the immune system responds by releasing a large number of inflammatory cytokines, such as interleukin 8(IL-8), eventually activating T lymphocytes and neutrophils and then uses the acquired immune system alongside innate immune system. 4 This results in activation of a vast range of signalling pathways, such as Janus kinase (JAK)/ signal transducer and activator of transcription (STAT) and nuclear factor kappa B (NF/κB). 5 In this case, if pathogenicity of the virus continues and more inflammatory cytokines are produced, the immune system activity can lead to serious tissue damage and even tissue death. Therefore, the immune system is one of the most challenging cases in the field of pathogenicity of this virus, which needs to be studied extensively. 6 Another issue in pathogenesis of the virus is the influence of vitamin D on its prevention. Many studies have demonstrated that the decreased levels of vitamin D are associated with the increased disease severity and mortality, and vitamin D supplements are widely recommended during the pandemic period. [7][8][9] The symptoms of SARS-CoV-2 are different from mild to severe.
Investigations have shown that the patients with COVID-19 infection have low lymphocyte counts, unusual respiratory results and elevated levels of pro-inflammatory cytokines in plasma. 10,11 MicroRNAs play a role in important biological processes, such as cellular metabolism, cell division, death, cell movement, intracellular signalling, immunity, apoptosis and oncogenesis. 12 They can interact with different regions of target mRNA including 5-UTR, gene promoter and coding sequences. Binding with 3-UTR causes posttranscriptional silencing rather than pairing with 5'-UTR or coding sequence, while miRNA interaction with promoter regions has been proved to induce transcription. 13 Viruses can also produce miRNAs in their genomes that can target host genes. To date, more than 250 new viral miRNAs have been discovered, making it possible to detect function and origin of virus-encoded miRNAs. However, a large proportion of these miRNAs belong to DNA virus, and only 30 mature miRNAs to RNA viruses have been identified. 14 Canonical biogenesis of miRNA is a multi-step process involving transcription of viral miRNA gene by RNA pol II and formation of a structure called as pri-miRNA, its breakdown by Drosha endonuclease and DiGeorge syndrome critical region 8 (DGCR8) complex in the form of a hairpin structure called as pre-miRNA, its extraction from the nucleus and maturation through Dicer endonuclease, eventually joining the RNAinduced silencing complex (RISC) and suppressing expression of host genes. 15 Cell penetration by SARS-CoV-2 is mediated by membrane fusion process, in which a +ssRNA-sense genome is delivered into the cytoplasm, serving as a template for virus replication process utilizing virally encoded RNA-dependent RNA polymerase (RdRp), and ultimately, protein synthesis in fact, miRNA biogenesis process in cytoplasmic viruses is performed by non-canonical miRNA biogenesis pathways in Drosha-or Dicer-independent manner. 16,17 It has recently been shown that one of the main mechanisms of virus pathogenesis is targeting key human genes and suppressing their expression through viral miRNAs. 18 Therefore, study of viral microRNAs allows us to become more familiar with mechanism of coronavirus pathogenesis. In fact, viral microRNAs play a role in completing virus life cycle. In this way, they try to maintain their survival and proliferation of their genomes by targeting genes of the immune system.  [20][21][22][23][24][25] Using bioinformatics approach, herein, novel SARS-CoV-2 encoded miRNAs were identified. Our findings demonstrated that SARS-CoV-2 miRNAs probably have putative role in virus pathogenesis and affect the host immune system and various physiological processes to take advantage of the prolonged refuge in host cell.
It was found that SARS-CoV-2 miRNA-targeted human host genes are involved in viral pathogenesis, such as respiratory system, cellular and immune pathways, and vitamin D pathways. Finding target genes of these miRNAs and their pathways could provide new insights into SARS-CoV-2 infection, pathogenesis and treatment design also; enrichment analysis of each of the host's target genes provides more information about the SARS-CoV-2 ( Figure 1).
As a result, the pre-miRNAs' sequence and positions were obtained. 27

| Prediction of potential pre-miRNAs and mature miRNAs
One of the most critical steps in predicting mature miRNA is selecting bioinformatics tools that could differentiate between two miRNA precursors, pre-miRNA and pri-miRNA. pri-miRNA, the precursor of pre-miRNAs, is structurally characterized by a terminal loop, a stem of approximately three helical turns that flanked by two basal unpaired sequences. 28 Sequence-Structure Motif Base (http://www.regul atory rna.org/webse rver/SSMB/pre-miRNA/ index.html), a pre-miRNA prediction web server, was used to differentiate between pri-miRNA and pre-miRNAs and predict mature miRNAs between a complex of the conserved stem-loop structures. 29 It utilizes effective PriMir and Mirident programs, accurately approved to predict pre-miRNA and RNAfold software to plot its F I G U R E 1 Flow chart of the methodology used in the study secondary structure. This resource consists of a set of Perl, Python and PHP programs.
Besides, miREval 2.0, a web server available at http://mimir na.cente nary.org.au/mirev al/, was applied to further validate resulting pri-miRNAs. It uses support vector machines (SVMs) as a useful machine learning method for predicting accurate miRNAs. 30 The least bases for the stem-loop were set to 22

| Obtaining SARS-CoV-2 miRNA target genes
The miRDB database (http://mirdb.org/) was implemented to identify the targeted genes by the predicted miRNAs with score 80 as a cut-off. 32 Human genes that are related to SARS-CoV-2 infection were extracted from the gene ontology database (http://geneo ntolo gy.org/). 33 Subsequently, a list of hub genes targeted by these miR-NAs was identified followed by running the overlap analysis between miRNA target genes and corresponding human genes in SARS-CoV-2 infection using the Venn diagram (http://bioin forma tics.psb.ugent. be/webto ols/Venn/).
A principal hypothesis regarding miRNA/target functional interaction is a conserved base pairing (Watson-Creek match) between miRNA and target at miRNA positions 2-7 or 8 (5´ upstream of miRNA). 34 Recognizing and confirming interplays between miRNAs and their downstream targets is essential in describing regulative capacity of miRNAs in complex systems controlling biological processes. 35 RNAhybrid program (https://bibis erv.cebit ec.uni bielefeld.de/rnahybrid/), which can determine the possible binding sites of multiple miRNAs in large target RNAs based on the minimum free energy of hybrid structure, was performed for reliable prediction of the most reliable energy hybridization of v-miRNAs and 3´UTR of their target mRNAs. 36 The higher the degree of complementarity between the seed region sequence in miRNA and the 3 UTR in the mRNA target gene, the more reliable the prediction, so according to the other previous studies, the minimum six-base seed match length was considered. For this evaluation, input data, that is 3 UTR of host candidate target genes (extracted from the miRDB database) and mature miRNA, were uploaded.

| Analysis of gene, function and pathway enrichment
For gaining mechanistic insight into our gene list, gene ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) (KEGG; https://www.genome.jp/kegg/) and over-representation analyses were conducted on the list of target genes using the R software cluster profiler package. 37 Significant groups were determined according to Benjamini-Hochberg correction and the adjusted cut-off levels of <0.05. In addition, Cytoscape tool (version 3.8.0) was used to create and visualize a miRNA-target gene network.

| Expression analysis
Gene expression omnibus (GEO), an organized database for functional genomic studies, is a free source of the results regarding microarray expression gene stored at NCBI (https://www.ncbi.nlm.nih.

| Cis-acting regulatory elements in 5′ regulatory regions
Investigation of cis-acting regulatory elements (CREs) can provide useful information on transcriptional regulation of the genes in-

| Identification of SARS-CoV-2 miRNA target genes
Herein, a bidirectional analysis was used in order to achieve genes targeted with these miRNAs. In the first step, 5,584 genes targeted by the miRNAs were found with an 80% of cut-off and they were  Tables 1 and 2). According to the results, all 32 genes were common between three considered pathways ( Figure 2). The miRNA-target gene network was also visualized using Cytoscape tool (Figure 3).

| Gene annotation and pathway analysis
For investigating the relevance of the target genes, these genes were grouped into GO terms. A total of 278 GO categories were significantly enriched including 189 terms of biological processes (BP), 20 molecular function (MF) terms and 69 terms of cellular components (CC). The most important GO terms are shown in Figure 4. The most highly significant GO categories within BP terms were identified as protein targeting, heat response and establishment of protein localization to organelle ( Figure 4A). In addition, regulation of cell cycle G2/M phase transition term was also overrepresented (supplementary data file 3). The significantly enriched MF terms were guanosine diphosphate (GDP) binding and glucosyltransferase activity ( Figure 4B). In the category of CC, organelle inner membrane, mitochondrial inner membrane and endoplasmic reticulum lumen were the top GO terms ( Figure 4C). The pathway analysis also revealed that target genes were enriched only in protein processing in endoplasmic reticulum and RNA transport pathways.

| Identification of promoter motifs
The 1000 bp of 5' regulatory region of target genes was scanned to identify the conserved motifs and consensus CREs. The top 10 significant motifs were discovered by MEME Suite, and motif-based sequence analysis tools, with widths ranging from 15 to 50aa in the target genes' promoter. GoMo tool predicted several interesting biological processes by further analysing the identified motifs using MEME tool ( Table 3). The flow chart of results is provided in Figure 5. For investigating this process, viral miRNAs and human target genes were studied and predicted to determine whether the SARS-CoV-2 miRNAs could play a significant role in the virus's pathogenicity.

| D ISCUSS I ON
Using bioinformatics tools, each algorithm has a special rate of false-positive and false-negative results. 45 Therefore, exploiting more than one algorithm is necessary to make reliable predictions    Results of a recent in vitro study confirmed that Ebola virus (EBOV)-derived pre-miRNAs are dependently processed by cellular miRNA processing machinery into subsequent mature miRNAs that like host miRNAs can directly silence their target mRNA. 47 Islam et al predicted several novel miRNAs produced by Zika virus.
They showed that one of the most important Zika virus pathogeneses might be caused by viral microRNAs. These miRNAs target the genes associated with cellular immunity and neurogenic functions. 47 Interestingly, results of two experimental studies, in which tickborne encephalitis virus (TBEV) and Sindbis virus were artificially modified by inserting an exogenous miRNA hairpin, confirmed that cytoplasmic RNA virus could express functional miRNA motifs, for example it was discovered that KUN-miR-1, produced by West Nile RNA virus, targets GATA binding protein 4(GATA4) and leads to virus replication. 48,49 In another study, it was shown that Dengue virus miRNA 'DENV-vsRNA-5' targets non-structural protein 1(NS1) and has a role in autoregulation of the virus. 50 Since transcription factors (TFs) play a role in immune response,

| CON CLUS ION
In this study, we investigated viral miRNAs as a principal factor in the pathogenicity of SARS-CoV-2. MicroRNAs can act on pathogenesis of the disease by targeting and reducing the expression of key genes involved in the immune system against viruses and the pathway of the respiratory system and vitamin D. We have predicted several new microRNAs produced by SARS-CoV-2 using bioinformatics tools, and as expected, the target genes of these miRNAs play an important role in lung tissue cells, vitamin D and inflammatory processes.

ACK N OWLED G EM ENTS
This work has been supported by grant (990268) from the office of the Vice Chancellor for Research, Hormozgan University of Medical Sciences, Bandar Abbas, Iran.

CO N FLI C T O F I NTE R E S T
The authors declare that they have no conflict of interests. F I G U R E 4 Gene Ontology (GO) enrichment analysis conducted using clusterProfiler R package for 32 target genes related to the respiratory system, cellular and immune pathways, and vitamin D pathways. All GO terms were identified with an adjusted enrichment P value of less than .05. A) A bubble plot of the enriched Biological Process (BP), B) a bubble plot of the enriched Molecular Function (MF), and C) a bubble plot of the enriched Cellular Component (CC). The size and colour of the nodes represent the gene number and the p. adjust of each GO terms, respectively F I G U R E 5 The flow chart of results used in the study