Computational Prediction of RNA–RNA Interactions between Small RNA Tracks from Betacoronavirus Nonstructural Protein 3 and Neurotrophin Genes during Infection of an Epithelial Lung Cancer Cell Line: Potential Role of Novel Small Regulatory RNA

Whether RNA–RNA interactions of cytoplasmic RNA viruses, such as Betacoronavirus, might end in the biogenesis of putative virus-derived small RNAs as miRNA-like molecules has been controversial. Even more, whether RNA–RNA interactions of wild animal viruses may act as virus-derived small RNAs is unknown. Here, we address these issues in four ways. First, we use conserved RNA structures undergoing negative selection in the genomes of SARS-CoV, MERS-CoV, and SARS-CoV-2 circulating in different bat species, intermediate animals, and human hosts. Second, a systematic literature review was conducted to identify Betacoronavirus-targeting hsa-miRNAs involved in lung cell infection. Third, we employed sophisticated long-range RNA–RNA interactions to refine the seed sequence homology of hsa-miRNAs with conserved RNA structures. Fourth, we used high-throughput RNA sequencing of a Betacoronavirus-infected epithelial lung cancer cell line (Calu-3) to validate the results. We proposed nine potential virus-derived small RNAs: two vsRNAs in SARS-CoV (Bats: SB-vsRNA-ORF1a-3p; SB-vsRNA-S-5p), one vsRNA in MERS-CoV (Bats: MB-vsRNA-ORF1b-3p), and six vsRNAs in SARS-CoV-2 (Bats: S2B-vsRNA-ORF1a-5p; intermediate animals: S2I-vsRNA-ORF1a-5p; and humans: S2H-vsRNA-ORF1a-5p, S2H-vsRNA-ORF1a-3p, S2H-vsRNA-ORF1b-3p, S2H-vsRNA-ORF3a-3p), mainly encoded by nonstructural protein 3. Notably, Betacoronavirus-derived small RNAs targeted 74 differentially expressed genes in infected human cells, of which 55 upregulate the molecular mechanisms underlying acute respiratory distress syndrome (ARDS), and the 19 downregulated genes might be implicated in neurotrophin signaling impairment. These results reveal a novel small RNA-based regulatory mechanism involved in neuropathogenesis that must be further studied to validate its therapeutic use.


Introduction
Over the last two decades, there have emerged highly pathogenic and deadly Betacoronaviruses (Beta-CoVs), including severe acute respiratory syndrome coronavirus (SARS-CoV) and Middle East respiratory syndrome coronavirus (MERS-CoV), which can cause severe human disease but have only been observed in a one-time event and limited outbreaks [1,2]. More recently, a novel Beta-CoV named SARS-CoV-2 has been identified as the source of coronavirus disease 2019 (COVID- 19), which continues to threaten the lives of hundreds of millions of people [3]. These positive-sense, single-stranded RNA (+ ssRNA) viruses belong to the Beta-CoV genus and carry one of the largest RNA genomes (~30 kilobases, kb) capped at the 5 ends and poly-A tail among all RNA virus families [4][5][6]. Upon cell entry, study [39]. Briefly, the prediction of structured regions across the whole viral genome was scanned in windows of length 120 nucleotides sliding by 40 nucleotides with RNAz (v2.1) [40]. Conserved RNA structures with p > 0.98, and z < −3 were retrieved and then tested to estimate their selection pressures using SSS-test (v1.0) [11].
RNA-seq expression profiling was obtained from the Gene Expression Omnibus (GEO) database (NCBI Reference series: GSE148729) of the epithelial lung cancer cell line  in mock and infected with SARS-CoV (GER-FRA/2003) and SARS-CoV-2 (USA-WA1/2020). For MERS-CoV, a complete transcriptional profile (RNA-seq) of circRNA, miRNA, and mRNA expression was analyzed in epithelial lung cancer cells (Calu-3) for mock and infected with MERS-CoV (HCoV-EMC/2012), downloaded from the GEO database (NCBI Reference series: GSE139516).

Differential Expression Analysis
Gene-level differential expression analysis was carried out using the edgeR package (v3. 36) in R [41] to determine the set of differentially expressed genes (DEGs) as a response to the viral infection. All samples were prefiltered to discard genes with a read count less than 2, and then for each group, genes with a minimum requirement of 20 counts per million (CPM) across libraries were kept. Count data were then normalized using trimmed means of M values (TMM) to adjust for sequencing library size difference [42]. For each virus, three comparisons were performed: (i) mock cells between the final and initial time (MM), (ii) infected cells upon post-infection against onset (II), and (iii) a comparison of infected cells with the respective mock cells (IM). We employed a quasi-likelihood F test to assess the significance of group differences [42]. Significant DEGs were defined with |log 2 (fold change, FC)| > 1.5 and false discovery rate (FDR)-adjusted p-value < 0.05 using Benjamini-Hochberg's procedure for multiple comparison adjustment.

hsa-miRNAs Targeting Beta-CoVs in Respiratory Epithelial Cells
Given that the upper respiratory tract probably represents the onset site for Beta-CoV infection [43], a systematic literature review was conducted to identify hsa-miRNAs involved in lung cell infection for SARS-CoV, MERS-CoV, and SARS-CoV-2. In this systematic review, databases such as Pubmed, Lilacs, EMBASE, Scopus, Web of Science Core Collection (WosCC), and EBSCO were searched using the terms "microRNA" and its synonyms (microRNA OR hsa-miRNA OR miR OR miRNA OR small noncoding RNA OR small ncRNA), "lung epithelial infection" (lung epithelial infection OR pulmonary epithelial infection OR respiratory epithelial infection) and in titles/abstracts as "SARS-CoV, MERS-CoV, and SARS-CoV-2", separately. Research articles with either robust computational or experimental validation were considered for the inclusion of hsa-miRNAs until July 2022.

RNA-RNA Interactions
For each virus, the putative 3 UTR regions of genes targeted by hsa-miRNAs were predicted using TargetScan [44], miRDB [45], and miRTarBase [46], implemented in the miRWalk 3.0 [47] online tool with a binding probability ≥ 1.0. These binding pairs were cross-validated with intercepted DEGs between comparisons II and IM. For a rigorous target conservation analysis, hsa-miRNAs sequences were retrieved from miRBase [48], while 3 UTRs sequences of their respective targets as DEGs from the UTRdb database [49]. Thermodynamic interaction was undertaken for each hsa-miRNAs:3 UTR pairing, preventing certain hsa-miRNAs from binding to the 3 UTR of other DEGs. We conducted these hybridization duplexes using RNAhybrid [50], miRanda [51], and mirTarP [52]. RNAhybrid queries were considered with a strict binding to seed region (nucleotides 2-8), a minimum free energy (MFE) ∆G cutoff ≥ −18 kcal/mol, and a maximum bulge and internal loop length of 1. The parameters for miRanda included a hybridization alignment score ≥140 and an MFE ∆G cutoff ≤ −18 kcal/mol, and for mirTarP, a consecutive base match of 7 and an MFE ∆G cutoff ≤ −18 kcal/mol. Finally, we focused only on hsa-miRNAs overlapped by the three methods with the ggVennDiagram package (v2.1) [53] implemented in R.

Prediction of Putative vsRNAs
Upon the determination of potential hsa-miRNA candidates, we sought the most energetically favorable hybridization sites targeting conserved RNA structures with the negative selection at each host to identify putative vsRNAs derived from genomes of SARS-CoV, MERS-CoV, and SARS-CoV-2. The interactions between hsa-miRNA:viralRNA were predicted using RNAhybrid, miRanda, and mirTarP with similar parameters as specified in the hsa-miRNAs:3 UTR pairing. The hsa-miRNA binding site to viralRNA was considered as a possible candidate vsRNA.

Validation of vsRNAs and Prediction of Their Potential Targets
To verify whether there is a significant similarity between predicted vsRNAs and host miRNAs, a search of miRNAs deposited in miRbase 22 [48] based on the miRNA seed sequence was performed. Then, to better understand what processes may affect this subset of Beta-CoV vsRNAs, assuming that they are acting as host-gene regulators, human genes potentially targeted were identified using two independent tools including Diana (software MR-microT) [54], and miRDB (Custom prediction) [45]. Only predicted targets with a score ≥ 70 were considered in both cases. In addition, a reliable set of targets was obtained by overlapping predicted genes with the two algorithms. Finally, this list of targets was cross-validated with the repertoire of intercepted DEGs, and those overlapping were retrieved.

Functional Enrichment Analysis
The functionality of these DEGs was analyzed by an over-representation analysis (ORA) of gene ontology (GO) terms using Protein Analysis Through Evolutionary Relationships (PANTHER) [55]. GO terms enriched in biological process (BP) and Reactome pathways were considered using a Fisher's exact tests with an FDR-adjusted p-value < 0.05. In addition, the Gene Cards database (https://www.genecards.org/, accessed on 16 March 2023) was used to reveal the human lung epithelial tissue, in which DEGs targeted by vsRNAs are expressed, and also to retrieve information on the pathologies that might be linked to these genes.

Over 80% of Conserved RNA Structures in Genomes of SARS-CoV, MERS-CoV, and SARS-CoV-2 Are under Negative Selection
This study employs conserved RNA structures that showed negative selective pressures for genomes of the three Beta-CoVs (s ≤ 2.99) in a wide variety of hosts, including bats, intermediate animals, and humans. We retrieved 505 conserved loci or regions containing RNA structures acting under negative selection (Table 1). Among these loci, 133 (26%) were detected in SARS-CoV-2, followed by SARS-CoV and MERS-CoV with 159 (31%) and 213 (42%), respectively. We found that 81% of 619 loci under different selection pressures were conserved RNA structures with negative signals.  1 Genome coverage percentage was calculated by multiplying the total number of nucleotides of all predicted loci by 100 and then dividing the viral genome length of a given host, shown in Figure S1.

Identification of DEGs in Beta-CoV-Induced Epithelial Lung Cancer
To validate downstream analyses, we conducted a differential expression analysis using RNA-seq data obtained from Calu  As the II and IM comparisons revealed the most significant number of DEGs, those in common were used for the validation of RNA-RNA interactions, target prediction, and functional enrichment. Accordingly, overlapping DEGs were identified, where 466 (331 up-and 135 downregulated) and 593 (374 up-and 219 downregulated) were accounted for SARS-CoV (Figure 2A,B), and SARS-CoV-2 ( Figure 2E,F), respectively. Conversely, as MERS-CoV obtained the most considerable amount of DEGs, only 1423 (727 up-and 696 downregulated) were retrieved ( Figure 2C,D)). The complete list of DEGs for the three viruses is available in Table S1. As the II and IM comparisons revealed the most significant number of DEGs, those in common were used for the validation of RNA-RNA interactions, target prediction, and functional enrichment. Accordingly, overlapping DEGs were identified, where 466 (331 upand 135 downregulated) and 593 (374 up-and 219 downregulated) were accounted for SARS-CoV (Figure 2A,B), and SARS-CoV-2 ( Figure 2E,F), respectively. Conversely, as MERS-CoV obtained the most considerable amount of DEGs, only 1423 (727 up-and 696 downregulated) were retrieved ( Figure 2C,D)). The complete list of DEGs for the three viruses is available in Table S1.

hsa-miRNAs Associated with Lung Physiopathology in Beta-CoVs
For the prediction of putative RNA-RNA interactions, our strategy first consisted of searching research articles published until July 2022 that focused on describing the role, function, and/or association of hsa-miRNAs in the lung physiopathology of Beta-CoVs. Based on the literature review, we found a total of 256 hsa-miRNAs, including for SARS-CoV (36 hsa-miRNAs in 4 articles from 2002), MERS-CoV (70 hsa-miRNAs in 5 articles from 2012), and SARS-CoV-2 (150 hsa-miRNAs in 22 articles from 2019) ( Table 2). Table  S2 describes these hsa-miRNAs for each virus.

hsa-miRNAs Associated with Lung Physiopathology in Beta-CoVs
For the prediction of putative RNA-RNA interactions, our strategy first consisted of searching research articles published until July 2022 that focused on describing the role, function, and/or association of hsa-miRNAs in the lung physiopathology of Beta-CoVs. Based on the literature review, we found a total of 256 hsa-miRNAs, including for SARS-CoV  Table 2). Table S2 describes these hsa-miRNAs for each virus. For each virus, hsa-miRNAs defined through the literature review were queried into miRWalk 3.0 to predict their putative 3 UTR regions of the target genes. Then, the hsa-miRNA:3 UTR pairings were cross-validated with previously identified DEGs, in which SARS-CoV obtained 52 pairings, MERS-CoV 210, and SARS-CoV-2 76, thus MERS-CoV had the most predicted pairings. To further decrease the impact of false favorable rates, the predicted interaction of hsa-miRNA:3 UTR was corroborated using scores from highly cited prediction tools such as RNAhybrid, miRanda, and mirTarP. These three tools matched 20 pairings in SARS-CoV, of which 15 (75%) showed to be of the let-7 family ( Figure 3A). In a lower proportion, from 31 and 21 pairings of MERS-CoV and SARS-CoV-2, 19 (61%) and 5 (24%) concerned the same family of hsa-miRNAs, respectively ( Figure 3B,C) (Table S3). In addition, it is worth mentioning that hsa-miRNA:3 UTR interactions showed a high binding homology to the seed sequence. For instance, for SARS-CoV, the predicted target sites were those where the MFE ranged from −19.4 to −30.7 kcal/mol, MERS-CoV between −18.0 and −30.6 kcal/mol, and SARS-CoV-2 between −18.2 and −32.2 kcal/mol (Table S3).

The let-7 Family of hsa-miRNAs Is the Most Frequently Predicted in hsa-miRNA:3′UTR Interactions
For each virus, hsa-miRNAs defined through the literature review were queried into miRWalk 3.0 to predict their putative 3′UTR regions of the target genes. Then, the hsa-miRNA:3′UTR pairings were cross-validated with previously identified DEGs, in which SARS-CoV obtained 52 pairings, MERS-CoV 210, and SARS-CoV-2 76, thus MERS-CoV had the most predicted pairings. To further decrease the impact of false favorable rates, the predicted interaction of hsa-miRNA:3′UTR was corroborated using scores from highly cited prediction tools such as RNAhybrid, miRanda, and mirTarP. These three tools matched 20 pairings in SARS-CoV, of which 15 (75%) showed to be of the let-7 family ( Figure 3A). In a lower proportion, from 31 and 21 pairings of MERS-CoV and SARS-CoV-2, 19 (61%) and 5 (24%) concerned the same family of hsa-miRNAs, respectively ( Figure  3B,C) (Table S3). In addition, it is worth mentioning that hsa-miRNA:3′UTR interactions showed a high binding homology to the seed sequence. For instance, for SARS-CoV, the predicted target sites were those where the MFE ranged from −19.4 to −30.7 kcal/mol, MERS-CoV between −18.0 and −30.6 kcal/mol, and SARS-CoV-2 between −18.2 and −32.2 kcal/mol (Table S3). From the final hsa-miRNA:3 UTR interactions, hsa-miRNAs were retrieved to identify which binding sites are shared with the loci seed sequence showing negative selection in individual virus hosts using RNAhybrid, miRanda, and mirTarP ( Figure S3). Interestingly, the analysis matched a repertoire of 25 hsa-miRNA:virus RNA pairings, where 20 of them were found in SARS-CoV-2, whereas MERS-CoV and SARS-CoV obtained 3 and 2 pairings, respectively (Table S4). However, there were hsa-miRNA:virus RNA interactions that matched identically in MERS-CoV and SARS-CoV-2, resulting in one and six pairings for each virus (Table S5). Finally, nine hybridization duplexes ranging from 18 to 23 nucleotides were analyzed.

Bat
Intermediate Human

ORF1a Is a Crucial Viral Gene in Modulating the Host Transcriptome
First, we asked whether these nine potential vsRNAs had significant similarities with miRNAs from any organism. To this end, the miRBase database was used to search based on the miRNA seed sequence. Surprisingly, none of our predicted vsRNAs showed miRNAs associated with other organisms, which is a promising result (Table S5). In light of this, the following purpose was to understand what processes may affect this subset of Beta-CoV vsRNAs, assuming that they act as host gene regulators. Diana (MR-microT software) [54] and miRDB (Custom prediction) [45] tools with a score ≥ 70 were employed to obtain a reliable set of candidate human genes potentially targeted by these vsRNAs. A total of 1998 human genes were predicted as possible vsRNA targets for the three viruses. Despite SARS-CoV having two vsRNA candidates, it registered the most significant number of targets compared to the other viruses, resulting in 1076 ranging from 102 (SB-vsRNA-S-5p) to 974 (SB-vsRNA-ORF1a-3p), 538 being the average number of targets across the vsRNAs predicted ( Figure S5A). Instead, SARS-CoV-2 with five potential vsRNAs scored 879 targets, between 82 (S2H-vsRNA-ORF1a-3p) and 393 (S2H-vsRNA-ORF1b-3p) averaging 175.8 ( Figure S5C). As MERS-CoV only has one vsRNA, 43 targets were predicted ( Figure S5B). The full details of each vsRNA targeting human genes can be consulted depending on the virus: SARS-CoV (Table S6), MERS-CoV (Table S7), and SARS-CoV-2 (Table S8).

vsRNAs Are Probably Shutting down Genes Associated with Neurotrophin Signaling Impairment upon Infection
An ORA was performed using GO terms with the PANTHER online tool to better understand the functional alterations induced by up-and downregulated DEGs following lung epithelium infection with Beta-CoVs. Figure 6A shows the DEGs regulated by SARS-

vsRNAs Are Probably Shutting Down Genes Associated with Neurotrophin Signaling Impairment upon Infection
An ORA was performed using GO terms with the PANTHER online tool to better understand the functional alterations induced by up-and downregulated DEGs following lung epithelium infection with Beta-CoVs. Figure 6A (Table S9).  (Table S9).
Unlike SARS-CoV and MERS-CoV, most of the SARS-CoV-2 vsRNAs were predicted in human virus genomes, allowing a closer approximation of how infection may alter gene expression in human epithelial cells ( Figure 6C). For instance, upregulated DEGs were enriched mainly by BP terms such as the regulation of natural killer cell proliferation (GO:0032819), interleukin-6-mediated signaling pathway (GO:0070102), regulation of tumor necrosis factor production (GO:0032760), and regulation of cytokine production involved in immune response (GO:0002720). DEGs engaged in these processes were highly  (Table S9).
Unlike SARS-CoV and MERS-CoV, most of the SARS-CoV-2 vsRNAs were predicted in human virus genomes, allowing a closer approximation of how infection may alter gene expression in human epithelial cells ( Figure 6C). For instance, upregulated DEGs were enriched mainly by BP terms such as the regulation of natural killer cell proliferation (GO:0032819), interleukin-6-mediated signaling pathway (GO:0070102), regulation of tumor necrosis factor production (GO:0032760), and regulation of cytokine production involved in immune response (GO:0002720). DEGs engaged in these processes were highly overexpressed, namely, the cytokine family (CXCL8 Although these BP terms are already widely known in human SARS-CoV-2 infection [56][57][58], the most remarkable findings concern downregulated DEGs, where ORF1a vsRNAs targeted 78%. Similar to SARS-CoV, these DEGs were involved in the regulation of presynaptic membrane potential (GO:0099505), which together with GRIK2  (Table S9).

Discussion
Since a picture is emerging in which various DNA viruses have been documented to encode v-miRNAs and undergo transcription and biogenesis similar to host miRNAs [59], whether RNA-RNA interactions of RNA viruses might function as putative vsRNAs has been a matter of controversy for several years. Nevertheless, several reports have demonstrated the functional identification v-miRNAs encoded by nuclear RNA viruses like retroviruses, which may benefit from host miRNA biogenesis machinery [21][22][23]60], as they have a DNA replication intermediate in their replication cycle and maintain long-term persistent infections [20,61]. Cytoplasmic RNA viruses do not confer this evolutionary advantage, as they may not access the nuclear enzyme Drosha needed for v-miRNA processing [62]. However, a recent study has shown that the hepatitis A virus may encode two novel functional v-miRNAs, hav-miR-1-5p and hav-miR-2-5p, which might hijack the host miRNA pathway [63]. Under these premises, we used RNA-RNA interactions of negative selection RNA structures to predict potential vsRNAs proposed as miRNA-like molecules, which might not fall into the more widely accepted category of miRNAs in the field [38].
In the wake of the SARS-CoV-2 pandemic, the discovery of vsRNA candidates in different Beta-CoVs has increasingly taken precedence. Unlike multiple reports that have shown computationally the presence of putative Beta-CoV-derived v-miRNAs (Table 2), our study employed RNA-RNA bindings of conserved RNA structures undergoing negative selection in the genomes of SARS-CoV, MERS-CoV, and SARS-CoV-2 circulating in different bat species, intermediate animals, and human hosts across the globe. Additionally, most studies predict v-miRNAs using the complete list of hsa-miRNAs from miRbase v22 [48] to be aligned against Beta-CoV genomes [64][65][66][67][68][69], and even in the Kyoto Encyclopedia of Genes and Genomes (KEGG) map pathways database [70], increasing the false positive rate of v-miRNAs considerably. Instead, we focused only on Beta-CoVs targeting hsa-microRNAs previously predicted upon complete validation by computational or experimental reports, coupled with sophisticated long-range RNA-RNA interactions from six highly cited prediction tools such as TargetScan [44], miRDB [45], and miRTarBase [46], integrated into miRWalk 3.0 [47] and corroborated with RNAhybrid [50], miRanda [51], and mirTarP [52], to refine the seed sequence homology of hsa-miRNAs with RNA structures under negative selection. More importantly, our results always are cross-validated with DEGs induced by human lung cells infected with Beta-CoVs for each RNA-RNA interaction.
Considering our RNA-RNA strategy, nine potential vsRNAs conserved in negatively selected RNA structures of Beta-CoVs are proposed. Among them, two vsRNAs are encoded by SARS-CoV circulating in bats (SB-vsRNA-ORF1a-3p; SB-vsRNA-S-5p), one vsRNA for MERS-CoV circulating in bats (MB-vsRNA-ORF1b-3p), and six vsRNAs for SARS-CoV-2 distributing in bats (S2B-vsRNA-ORF1a-5p), intermediate animals (S2I-vsRNA-ORF1a-5p), and humans (S2H-vsRNA-ORF1a-5p; S2H-vsRNA-ORF1a-3p; S2H-vsRNA-ORF1b-3p; S2H-vsRNA-ORF3a-3p). Remarkably, we report for the first time that two SARS-CoV-2 vsRNAs from intermediate animals (S2I-vsRNA-ORF1a-5p) and humans (S2H-vsRNA-ORF1a-5p) at positions 4153 to 4189 (ORF1a) match in both sequence and structure with an identical MFE, indicating the RNA structure is strongly under negative selection during viral evolution [9,11,71]. Efforts addressed to identify vsRNAs in the different Beta-CoV hosts start from scratch. Most are dedicated to human-associated SARS-CoV-2, given the context of the current clinical crisis. Nevertheless, a recent study has been reported for SARS-CoV, which falls into the intermediate host category. Morales and coworkers identified three vsRNAs in lung derived from SARS-CoV-infected mice using deep sequencing. This study demonstrates that vsRNA molecules derived from N protein and nonstructural protein 3 (nsp3) of ORF1a contribute to SARS-CoV lung pathogenesis, which is reduced by the antagonization of these vsRNAs [72]. There are no findings of possible vsRNA candidates concerning MERS-CoV, and SARS-CoV-2 is highlighted from infected human cells. For instance, a Northern blot assay-based study demonstrated that SARS-CoV-2 could encode miRNA-like vsRNA, namely vmiR-5p, from its ORF7a associated with the pathogenesis of the virus [73]. In addition, the deep sequencing of infected Calu-3 and Vero E6 cells showed that the SARS-CoV-2 N protein might encode vsRNAs v-miRNA-N-28612, v-miRNA-N-29094, and v-miRNA-N-29443 with a positive association with viral load in COVID-19 patients. More recently, Singh and coworkers identified a CoV2-miR-O7a as vsRNA in the ORF7a using small sequence RNAs from SARS-CoV-2-infected human A549-ACE2 cells, which have the potential to interfere with host transcripts involved in IFN signaling [74]. These studies were the stepping stones that paved the way for further exploring that SARS-CoV-2, a cytoplasmic RNA virus, might generate vsRNAs during infection.
Interestingly, it has been suggested that many RNA viruses contain "hotspots" that serve as sites for vsRNA production [62], whereby sequences surrounding hotspots might be expected to adopt stable hairpin-like structures predicted by RNAz within the conserved RNA structures [75,76]. This pattern was recently observed by comparing virally derived small RNAs in different SARS-CoV-2-infected samples that are prone to generate vsRNAs [77]. In our case, the vast majority of our Beta-CoV vsRNAs encoded in ORF1a at genome positions ranging from 4153 to 5544 seem to show possible "hotspots", indicating that the nsp3 of ORF1a is deeply involved in vsRNA production [78,79]. Intriguingly, Fu and coworkers found a vsRNA in the 3 arm of a hairpin within nsp3, namely miR-nsp3-3p, which appears to indicate the risk of occurrence of critical illness among COVID-19 patients [80]. Indeed, a computational study discovered the expression of eight putative vsRNAs in SARS-CoV-2-infected human lung cancer cells using deep learning. They also found vsRNAs in viral ORF1a that may regulate the host transcriptome upon infection [34]. These results indicate that the nsp3 of ORF1a appears to be a prominent hotspot to produce specific vsRNAs from bat virus transmission to humans. However, further functional studies are required to validate this result, especially for Omicron SARS-CoV-2 variants [81][82][83].
Regarding the function of potential vsRNAs, it is tempting to speculate that microRNAlike molecules might modulate host gene expression upon infection. A limitation of this study relies on RNA-seq data from the human lung cancer cell line (Calu-3) infected with SARS-CoV, MERS-CoV, and SARS-CoV-2, which assumes that the predicted vsRNAs in viruses isolated from bats and different animal species may, in part, also infect the human upper respiratory tract and promote their post-transcriptional regulation. Remarkably, among 74 DEGs, we found 55 up-and 19 downregulated DEGs upon infection according to the II and IM comparisons between the time points assessed. This may suggest that our vsRNA candidates would be produced either at different stages of the viral replication cycle or depending on the viral load in the cell, triggering probably post-transcriptional gene silencing under particular conditions [84]. Regarding upregulated DEGs, it would be expected that vsRNAs encoded by bat-associated SARS-CoV increased the expression levels of KLF4 and PPARGC1A in human lung epithelium, which are involved in the gene transcription pathway. Likewise, the unique vsRNA predicted in MERS-CoV bats appears to induce the differentiation of human T-helper 17 cells via SMAD7 and RC3H1.
In contrast, as most SARS-CoV-2 vsRNAs were detected in human viruses, we observed an over-expression of critical genes associated with the acute respiratory distress syndrome (ARDS) seen in COVID-19 patients, such as cytokines (CXCL8 and CXCL11) [85], interleukins (IL-6) [86], and proinflammatory genes (JAK2, STAT2, and MXD1) [87]. Nevertheless, the premise of this study is that 83% and 73% of downregulated DEGs are targeted by vsRNAs encoded by the nsp3 of ORF1a, produced in bat SARS-CoV and human SARS-CoV-2 genomes, respectively. Most importantly, these vsRNAs appear to be shutting down genes implicated in the mechanisms regulating synaptic function, such as GRIK2, L1CAM, and NEFL, which are largely responsible for the neurotrophin signaling impairment of SARS-CoV-2 [88,89]. Although the role of SARS-CoV-2 in neurotrophins is an emerging trend, it has been reported that patients positive for SARS-CoV-2 have increased levels of brain-derived neurotrophic factor (BDNF) [90,91], a pivotal regulator of synaptic neuroplasticity, where its dysregulation might mediate in the pathogenesis of diverse neurological diseases such as Alzheimer's disease (AD) [92]. To the best of our knowledge, we report the first potential RNA-RNA interactions capable of acting as vsRNAs linked to SARS-CoV-2 neuropathologies. Although recent publications show experimental evidence that there exists an interphase of RNA-RNA interactions on A549, Vero E6, and A549hACE2 cells infected with SARS-CoV-2, their findings cannot draw conclusions in other cell types [93,94]. Due to the limitation of our survey, we claim that further experimental studies are warranted to confirm the legitimate identity of these putative vsRNAs detected in Calu-3 cells infected with SARS-CoV-2 to understand their function and to implement strategies to be used as potential biomarkers to interrupt the course of SARS-CoV-2-induced neurological manifestations.

Conclusions
We implemented a novel computational approach to predict nine potential vsRNAs from negatively selected RNA structures of SARS-CoV, MERS-CoV, and SARS-CoV-2 circulating in bats, intermediate animals, and humans, coupled with a complete review of Beta-CoVs targeting hsa-miRNAs and sophisticated long-range RNA-RNA interactions as well as validation using RNA-seq data. We found that most of our vsRNA candidates are encoded by the nsp3 of ORF1a, which transcriptionally dysregulated 74 DEGs in human epithelial cells (Calu-3). Among them, 55 upregulate molecular mechanisms underlying ARDS, whereas 19 downregulated DEGs might be implicated in neurotrophin signaling impairment. However, further experimental studies are needed to consider our predicted vsRNAs as potential biomarkers against Beta-CoV pathologies.
Supplementary Materials: The following supporting information can be downloaded at https: //www.mdpi.com/article/10.3390/v15081647/s1, Figure S1: Flowchart showing approach to data collection, curation process, and range of sequence lengths for the three Beta-CoVs used in the study by [39], Figure S2: Volcano plots showing the differential expression profiling of mock cultures of Calu-3 cells infected with Beta-CoVs between the final and initial time (MM), Figure S3: Potential common miRNA:virusRNA pairings predicted by RNAhybrid, miRanda, and mirTarP in Beta-CoVs, Figure S4: Comparison of vsRNAs encoded by Beta-CoVs across hosts of bats, intermediate animals, and humans, Figure S5: Number of predicted targets in the human genome for each vsRNA candidate in the three Beta-CoVs, Table S1: DEGs of Beta-CoVs using a |log 2 FC| > 1.5 and FDR < 0.05, Table S2: Full list of hsa-miRNAs targeting Beta-CoVs highly supported by either computational or experimental studies [35,[66][67][68][69], Table S3: Number of hsa-miRNA:3 UTR interactions resulting