SARS-CoV-2 encoded microRNAs are involved in the process of virus infection and host immune response

The outbreak of COVID-19 caused by SARS-CoV-2 is spreading worldwide, with the pathogenesis mostly unclear. Both virus and host-derived microRNA (miRNA) play essential roles in the pathology of virus infection. This study aims to uncover the mechanism for SARS-CoV-2 pathogenicity from the perspective of miRNA. We scanned the SARS-CoV-2 genome for putative miRNA genes and miRNA targets and conducted in vivo experiments to validate the virus-encoded miRNAs and their regulatory role on the putative targets. One of such virus-encoded miRNAs, MR147-3p, was overexpressed that resulted in significantly decreased transcript levels of all of the predicted targets in human,i.e., EXOC7, RAD9A, and TFE3 in the virus-infected cells. The analysis showed that the immune response and cytoskeleton organization are two of the most notable biological processes regulated by the infection-modulated miRNAs. Additionally, the genomic mutation of SARS-CoV-2 contributed to the changed miRNA repository and targets, suggesting a possible role of miRNAs in the attenuated phenotype of SARS-CoV-2 during its evolution. This study provided a comprehensive view of the miRNA-involved regulatory system during SARS-CoV-2 infection, indicating possible antiviral therapeutics against SARS-CoV-2 through intervening miRNA regulation.


Introduction
The current pandemic of COVID-19, caused by a novel virus, SARS-CoV-2, has led to over 24 000 000 infections and 830 000 fatalities in over 100 countries since its emergence in late 2019 (data collected on △ These authors contributed equally to this work. August 30, 2020). SARS-CoV-2 is an enveloped, positive-sense, single-stranded RNA betacoronavirus. Human coronaviruses cause respiratory tract infections that can be mild, such as some cases of the common cold. However, over the past two decades, highly pathogenic human coronaviruses have emerged, including SARS-CoV in 2002 and 2003, which infected 8000 people worldwide and with a fatality rate of approximately 10% and MERS-CoV in 2012 with 2500 confirmed cases and a death rate of 37% [1] . Compared to MERS and SARS, SARS-CoV-2 appears to spread more efficiently, although with relatively lower mortality of 2% to 3% [1][2] .
The SARS-CoV-2 gains entry to cells by engaging the angiotensin-converting enzyme 2 (ACE2) as receptor with the spike (S) protein, and then using the host serine protease TMPRSS2 for S priming [3] . Pneumonia is the common symptom of SARS-CoV-2 infection [4] ; however, damages to the liver [5] and spleen [6] , as well as gastrointestinal symptoms [7] were also reported. The SARS-CoV-2 infection causes increased secretion of cytokines, and the leading death cause of SARS-CoV-2 is ARDS, for which one of the main mechanisms is the cytokine storm, the deadly uncontrolled systemic inflammatory response resulting from the release of large amounts of proinflammatory cytokines and chemokines by immune cells [4] . Nevertheless, the mechanism of virus invasion and release, and the causes of these exuberant inflammatory responses in SARS-CoV-2 infection remain mostly unknown.
MicroRNAs (miRNAs) are widely found in plants, animals, and some viruses and involved in a variety of biological processes. In most cases, miRNAs interact with the 3 ′ untranslated region (3 ′-UTR) of target mRNAs to induce mRNA degradation and translational repression. However, except for the role as a repressor, increasing evidence has shown that miRNAs could activate gene expression by interacting with other regions, including the 5 ′-UTR and enhancers [8][9] . As crucial regulators of gene expression, miRNA-targeted therapeutics have been proposed for the treatment of cancers, virus infection, and other diseases. For example, antimiR targeted at has-miR-122 is in trials for treating hepatitis [10] .
It has shown that both DNA and RNA viral genomes can encode miRNAs [11][12] . The virus-derived miRNAs can be expressed in host cells and participate in the lifecycle and cellular consequences of infection. To date, over 300 viral encoded miRNAs have been described, and many of them were involved in the host immune response, cell cycle, and vesicle trafficking [13] . Viral miRNAs can target a large number of host genes involved in regulating cell proliferation, apoptosis, and host immunity [11] . The miR-HA-3p encoded by H5N1 accentuates the production of antiviral cytokines, resulting in a high level of cytokine production and high mortality [14] . Also, viral miRNAs target the virus genes, which also help the virus to remain hidden from the host immune response. For example, miR-BART2 encoded by EBV lie antisense to the EBV DNA polymerase gene BALF5 and inhibit DNA polymerase expression by inducing cleavage within its 3′-UTR [15] .
Furthermore, viruses have evolved mechanisms to degrade, boost, or hijack cellular miRNAs to benefit the viral life cycle [16] . The non-translated region of North American eastern equine encephalitis virus genome base pairs with the miR-142-3p, resulting in a miR-142-3p mediated innate immune suppression [17] . Moreover, recent studies indicate the virus-encoded miRNAs, host-encoded miRNAs, and miRNA targets together form a novel regulatory system between the virus and the host, which contributes to the outcome of infection [18][19] . Hence, it is crucial to comprehensively understand the role of the regulatory system in SARS-CoV-2 infection. Here, we showed that the virus-encoded miRNAs might participate in the process of virus infection and host immunity, and the attenuated phenotype of mutant SARS-CoV-2 might partly be attributed by the mutation-introduced alteration of the miRNA repository and targets.

Computational identification of virus-encoded miRNAs
The complete genome sequence of SARS-CoV-2 (accession number: MN908947) was obtained from the National Center for Biotechnology Information (NCBI). The genome sequence of SARS-CoV, MERS-CoV-2, and HCoV-NL63 were obtained from NCBI with accession numbers of NC_004718.3, NC_005831.2, and NC_019843.3, respectively. The VMir software package (v2.2) [20] was used to identify potential pre-miRNAs within the virus genome. The parameters of window and step sizes were set at 500 and 1 nucleotide (nt), respectively. The minimum free energy (MFE) of the predicted hairpin structures was predicted with RNAFold algorithm from the Vienna package [21] . Then HuntMi software was used to further distinguish miRNA hairpins from pseudo hairpins with a virus-based model [22] . Maturebayes web server [23] was used to identify the mature miRNAs within the pre-miRNAs based on sequence and secondary structure information of their miRNA precursors.

Prediction of miRNA targets and functional annotation
Human 3′-and 5′-UTR, and the enhancer sequence for five tissues were downloaded from the UTRdb database [24] and EnhancerAtlas database [25] , respectively. Human miRNA sequence was downloaded from the TargetScan database. The targets of virus-encoded miRNAs and human miRNAs were predicted with TargetFinder [26] and miRanda [27] . The miRNA targets on human genes were extracted from the TargetScanHuman 7.2 database [28] , with context++ score <-0.4 and percentile >0.99 as previously applied [29] . Then UTR and enhancer associated genes were extracted from annotations in UTRdb database [24] and EnhancerAtlas database [25] . Genes were functionally annotated using Gene Ontology (GO) terms using R packages, and the function enrichment tests were conducted with GOstats [30] .

Virus variation processing
The genomic variations were downloaded from China National Center for Bioinformation; data were collected until Mar 19, 2020. The mutations were obtained based on sequence alignment against the reference genome MN908947.3, as we used for miRNA scanning.

SARS-CoV-2 infection of Vero E6 cell line
The African green monkey kidney Vero E6 cell line was purchased from the Cell Resources Center of Shanghai Institute of Life Science, Chinese Academy of Sciences (China) and cultured in DMEM medium (Cat. No. 12430-054, Gibco, USA) containing 10% fetal bovine serum (FBS; Gibco) at 37 °C with 5% CO 2 atmosphere.
BetaCoV/JS03/human/2020 (EPI_ISL_411953), a SARS-CoV-2 virus strain, was isolated from nasopharyngeal swab of a 40-year-old female confirmed with COVID-19 by reverse transcriptase PCR (RT-PCR) in December 2019. The virus was propagated in Vero E6 cells, and the viral titer was determined by the 50% tissue culture infective dose based on microscopic observation of cytopathic effects. All the in vitro SARS-CoV-2 infection experiments were performed in a bio-safety level-3 laboratory of Jiangsu Provincial Center for Diseases Control and Prevention, Jiangsu, China.
Vero E6 cells were seeded into 12-well plates with a density of 4×10 4 cells/well for incubation in DMEM medium supplemented with 10% FBS for 16 hours in an incubator with 5% CO 2 at 37 °C for cells to reach 80% confluent. Six wells were infected with the SARS-CoV-2 virus with multiplicity of infection of 0.05 at 37 °C with a 5% CO 2 atmosphere for 48 hours. Then the infected cells of each well were harvested for RNA extraction and quantitation. RNA was extracted from cells using the RNeasy Plus Mini Kit (Cat. No. 74134, Qiagen, Germany) according to the manufacturer's instructions.

Statistical analysis
Data are presented as the mean±SD. P-values less than 0.05 indicated statistical significance. Differences between two samples were analyzed by Student's t test or Wilcoxon signed-rank test.

Identification of SARS-CoV-2 encoded miRNAs
A total of 808 potential pre-miRNAs in the SARS-CoV-2 genome were detected with VMir Analyzer [20] , an ab initio prediction program explicitly designed to identify pre-miRNAs in viral genomes. Then the candidate pre-miRNA sequences were further filtered using HuntMi [22] with the virus model, and the sequences with MFE larger than −20 kcal/mol were removed. A total of 45 candidate virus pre-miRNAs were finally obtained (Fig. 1A, The innermost cycle demonstrated the distribution of the candidate pre-miRNAs across the SARS-CoV-2 genome with red and green representing pre-miRNAs in the direct direction and reserve direction, respectively. The middle layer showed the mutation in the genome with red denoting SNPs, and blue denoting deletions. Outside of it is the frequency of the variations. C: The hairpin structure of the example pre-miRNAs putatively encoded by SARS-CoV-2. D: RT-qPCR analysis of the candidate SARS-CoV-2-encoded miRNAs expression levels in Vero E6 cells infected with SARS-CoV-2 (n=3, data are presented as mean±SD). Two-sided Student's t-test, * P<0.05, *** P<0.001. MFE: minimum free energy. Table 2, available online), with 30 in the direct orientation and 15 in the reverse orientation (Fig. 1B). The average length of the candidate pre-miRNA sequence was 78 nt with an average MFE of −28.1 kcal/mol. Finally, the sequence and position of mature miRNA within the candidate pre-miRNA were identified with MatureBayes, resulting in 90 putative mature miRNAs.

Virus miRNAs regulate host gene expression
Forty candidate miRNAs were predicted to bind at the 3 ′-UTR of 73 genes (Supplementary Table 3, available online). Gene ontology analysis revealed that the target genes were enriched in the apoptosisrelated functions, such as Notch binding, singlestranded DNA endodeoxyribonuclease activity, and deoxyribonuclease activity ( Fig. 2A). Apoptosis is a powerful mechanism to curtail viral spread [31] . However, viruses have also evolved sophisticated molecular strategies to subvert host cell apoptotic defenses. A couple of studies have investigated the mechanism by which viruses modulate apoptosis signalling [13,32] . CHAC1 and RAD9A were targeted by the virus-encoded miRNAs, MD2-5p and MR147-3p, respectively (Fig. 2B). CHAC1 is a proapoptotic enzyme and a regulator of Notch signalling. RAD9A is a cell cycle checkpoint protein and regulates cell death and promotes apoptosis [33] . FOXO3, which was predicted to be targeted by MR359-5p, is required for antioxidant responses and autophagy, and altered expression of FOXO3 is observed in hepatitis C infection and fatty liver [34] . The suppressive role of the candidate SARS-CoV-2-encode miRNAs on these genes suggested their possible role in reducing the host cell apoptosis to subvert host defense.
Eleven candidate virus miRNAs were predicted to bind at the 5′-UTR of 13 target genes (Supplementary Table 3). For example, MR385-3p was predicted to bind at the 5′-UTR of TGFBR3, a gene widely expressed on cells of the innate and adaptive immune system and promote Th1 differentiation and regulation of regulatory T-cell activation and survival [35] . MD2-5p was targeted at ENOX1, which encodes a protein that presents in the cytoskeletal compartment of endothelial cells and suppression of it in human or mouse endothelial cells inhibits migration and the ability to form tubule-like structures in matrigel [36] . miRNAs can also target the enhancer regions to activate their expression [9] . Targets of the candidate virus-encoded miRNAs on the enhancer sequence of five tissues (lung, gut, spleen, liver, and heart) were searched. Target prediction revealed that the lung had got the most significant number of genes targeted by the candidate miRNAs on the enhancer, followed by spleen and gut (Supplementary Fig. 1, Supplementary  Table 4, available online). Functional annotation of these targeted genes revealed enrichment in distinct functions in these tissues (Supplementary Fig. 1B-E). Significant enrichment in cell skeleton-related and immune response-related functions was observed in the lung (Supplementary Fig. 1B) and spleen (Supplementary Fig. 1C), respectively.
The four validated miRNAs above targeted the UTR sequences of six genes in the Vero E6 cell (Supplementary Table 5, available online). We then verified the expression of them in the virus-infected cells (two genes that encode putative proteins were not tested, i.e., LOC103242070 and LOC103217999). MR359-5p targeted at 3′-UTR of the same genes (i.e., FOXO3, GHR, and TIFA) in the green monkey genome as in human. Two of these three genes were significantly downregulated in the infected cells (Fig. 2C). Besides, GPR1 (G-protein-coupled receptor 1), which was targeted by MR359-5p at the 5 ′-UTR, was highly upregulated in the infected cells (Fig. 2C). G-protein-coupled receptors (GPCRs) constitute the largest group of cell-surface proteins that participate in a wide variety of physiological functions. GPCRs are involved in facilitating viral propagation and thereby contributing to viral pathogenesis [37][38] , and GPCR antagonists have broad antiviral activity [39] .
To further validate the regulatory role of these miRNAs on human genes, the construct of pGCMV/EGFP/MR147-3p was transfected into the pulmonary epithelial BEAS-2B cells. The increased MR147-3p levels relative to empty vector were observed (Supplementary Fig. 2, available online). Real-time qRT-PCR showed that overexpression of MR147-3p resulted in significantly decreased transcript levels of all of the predicted targets in human, i.e., EXOC7, RAD9A, and TFE3 (Fig. 2D). EXOC7 is a component of the exocyst complex, a protein complex that guides membrane addition and polarized exocytosis [40] . TFE3 is a transcription factor that promotes the expression of genes downstream of transforming growth factor beta signalling and also plays a crucial role in the regulation of lipid and glucose metabolism [41] .

Virus genome hijacks host miRNAs to regulate the immune response
Several studies have reported that the virus genome can hijack host miRNAs to modulate host biological processes [38] . A total of 28 human miRNAs were predicted to target the SARS-Cov-2 genome (Supplementary Table 6, available online). Since these miRNAs were hijacked by the virus genome, we supposed that the genes originally targeted by these miRNAs could be affected. There were over 800 human genes predicted to be regulated by these miRNAs, and a significant enrichment at the immune system process was observed ( Fig. 3A and B). Among the hijacked miRNAs, miR-146 is one of the critical modulators of the immune response. Studies in mice have shown that miR-146b −/− mice displayed enlarged spleens and increased myeloid cell populations in both spleen and bone marrow, and spontaneously developed intestinal inflammation [42] . Another miRNA, miR-939, was proposed to regulate proinflammatory   genes. Increasing miR-939 levels should restore homeostasis by decreasing inflammatory protein synthesis [43] . Our result posed a possibility that the hijack of human miRNAs by the SARS-Cov-2 genome might contribute to the abnormal immunity activation in patients with COVID-19.

Virus genomic region targeted by virus and host miRNAs
We then explored the targets of both virus and host cell-derived miRNAs on the SARS-CoV-2 genome. There were 27 targets of the putative virus miRNAs at the virus genome (Fig. 4A). Forty-three targets were predicted for these miRNAs, and most of them were bound to the ORFlab region (Fig. 4B), the longest ORF in the SARS-CoV-2 genome. There were two miRNAs target at the 5 ′-UTR of the virus genome, and two at the S gene. Twenty-eight human miRNAs were predicted to have 30 targets on the SASRS-CoV-2 genome (Fig. 4B). Most of the human miRNAs were bound to the ORF1ab, followed by seven miRNAs targeted at the N genes of the virus, and two at the 5′-UTR. The hsa-miR-4661-3p was predicted to target at the genomic region 25296 to 25320 within S genes (21563 to 25384), which is the potential 3′-UTR of the S gene transcript (Fig. 4C), suggesting a possible repressor role on the S gene expression. Since the transcriptome landscape of the SARS-Cov-2 is mostly unknown, further research in characterizing its genomic and transcriptomic features may enlighten the role of miRNAs in regulating the virus gene expression and infection.

Mutations in virus genome affect miRNA function
Mutations accumulated since the SARS-CoV-2 emerged, and studies have implicated that the variations on the virus genome may lead to an attenuated phenotype of SARS-CoV-2 [44] . Therefore, we explored whether miRNAs encoded by either virus or human contributed to the fitness of SARS-CoV-2 to human. A total of 810 mutations on the SARS-CoV-2 genome were collected from China National Center for Bioinformation, including 646 SNPs, 19 indels, and 145 deletions (Fig. 1B). There were 33 candidate pre-miRNAs overlaps the mutant region, including 30 pre-miRNAs overlapped with SNPs and four located in the deletion region. Compared with the wild-type pre-miRNAs, the MFE of the mutant pre-miRNAs were significantly increased (paired Wilcox test Pvalue <0.001, Fig. 5A), indicating an overall decreased  stability of the mutant-genome encoded pre-miRNAs. Next, we investigated the alteration of host miRNA targets on the mutant virus genome. Twenty-two human miRNAs can bind to the mutant genome, with 16 overlapped with those binding to the wild type genome (Fig. 5B). The genes targeted by the newly recruited miRNAs were annotated to the function of epithelial cell differentiation and actin filament fragmentation (Fig. 5C). In contrast, the genes targeted by the miRNAs specifically binding to the wild type genome were annotated to the immune response-related progress (Fig. 5D). The two miRNAs discussed, hsa-miR-939-5p and hsa-miR-146b-3p, which were involved in maintaining homeostasis by decreasing inflammation, failed to be hijacked by the mutant virus genome. Taken together, we speculate that the mutant type of virus might relieve the high level of immunity response, even the cytokine storm. Nevertheless, it enhances the invasion of the virus by intensifying the regulation of the epithelial differential and actin filament fragmentation, through attracting different batches of host miRNAs.

miRNA view of the difference between SARS-CoV-2 with other human coronaviruses
Compared with MERS-CoV and SARS-CoV, SARS-CoV-2 appears to have high transmission rates, but low mortality [1][2] . To explore whether and how miRNAs implicated in this difference between them, we predicted the miRNAs encoded by MERS-CoV, SARS-CoV, as well as HCoV-NL63, a coronavirus that causes mild symptoms in human and uses the same cell acceptor with SARS-CoV and SARS-CoV-2 (ACE2). The largest putative virus-encoded miRNA repository (212 miRNAs) was detected from the MERS-CoV genome ( Table 1), followed by SARS-CoV with 168 miRNAs and SARS-CoV-2 with 90 miRNAs. HCoV-NL63 encodes the least number of miRNAs (64 miRNAs). Additionally, the number of targets on the virus genome bound by the host miRNAs follows the same trend ( Table 1). Targets and functional analysis also revealed a dramatic difference between these virus (Supplementary Fig. 3  and 4, available online). These results implicated that differences in miRNA repositories and miRNA-target interaction might be involved in the pathological difference among the four coronaviruses.

Discussion
In this study, we revealed a systematic view of the possible role of the candidate virus-encoded miRNAs as well as the host-derived miRNAs in the SARS-CoV-2 infection (Fig. 6). The analysis revealed that the immune response was the function most affected by the virus-infection introduced miRNA change in  Fig. 6 A graphic summary of the findings in this study. The SARS-CoV-2 encoded miRNAs are able to target the regulatory regions (3′-UTR, 5′-UTR or enhancer) of host genes, and the host miRNAs might be hijacked by the virus genome. The involved miRNAs could target genes that participate in variety of biological process that induced by virus infection. repository and targets, which might contribute to the immune evasion of the virus as well as the abnormal immunity activation in the host. The SARS-CoV-2encoded miRNAs were also implicated in the cytoskeleton dynamics facilitating the virus envision, trafficking within the cell, and release. The virusencoded miRNAs were also predicted to be able to repress the expression of genes involved apoptosis by binding at the 3′-UTR of host genes. We validated the virus-encoded miRNAs and their targets in the Vero E6 and human cell lines. Additionally, the genomic mutation-introduced alterations on the coding ability of virus miRNAs and the targets of both virus and host miRNAs might contribute to the attenuated phenotype of SARS-CoV-2 during its evolution.
Two virus miRNAs were predicted to target at potent drug targets, which were MR147-3p that bind to the enhancer of TMPRSS2 and MD40-5p that bind to the enhancer of NUP214 (Supplementary Table 4, available online). TMPRSS2 enhances SARS-CoV-2 infection, together with ACE2 [3] . NUP214 is predicted as a drug target of Verdinexor, a selective inhibitor of nuclear export for virus RNP under clinical trial. Animal experiments have proved the efficiency of Verdinexor in reducing lung viral loads and proinflammatory cytokine expression in influenza A virus-infected mice [45] . We have quantified the expression of TMPRSS2 in the MR147-3p transfected cells, and elevated average expression was observed, but the variation between different replicates was large (data do not show). Since the mechanisms of miRNAs mediated RNA activation are largely unknown, for example, how the miRNAs being imported back into the nucleus [9] , more studies are needed to explore the regulation of gene expression of miRNAs by binding enhancer.
In this work, we did exhaustive in silico analysis of the SARS-CoV-2 genome from the perspective of miRNA and aided by simple validations for the virus miRNAs and their targets. Though several interesting phenomena have been observed, more detailed and sophisticated experiments are needed to validate them. With further experiments to validate the role of candidate miRNAs in vivo, there is a reasonable prospect to develop antiviral therapeutics against SARS-CoV-2 through targeting the candidate miRNAs.