In-Depth Bioinformatic Analyses of Nidovirales Including Human SARS-CoV-2, SARS-CoV, MERS-CoV Viruses Suggest Important Roles of Non-canonical Nucleic Acid Structures in Their Lifecycles

Non-canonical nucleic acid structures play important roles in the regulation of molecular processes. Considering the importance of the ongoing coronavirus crisis, we decided to evaluate genomes of all coronaviruses sequenced to date (stated more broadly, the order Nidovirales) to determine if they contain non-canonical nucleic acid structures. We discovered much evidence of putative G-quadruplex sites and even much more of inverted repeats (IRs) loci, which in fact are ubiquitous along the whole genomic sequence and indicate a possible mechanism for genomic RNA packaging. The most notable enrichment of IRs was found inside 5′UTR for IRs of size 12+ nucleotides, and the most notable enrichment of putative quadruplex sites (PQSs) was located before 3′UTR, inside 5′UTR, and before mRNA. This indicates crucial regulatory roles for both IRs and PQSs. Moreover, we found multiple G-quadruplex binding motifs in human proteins having potential for binding of SARS-CoV-2 RNA. Non-canonical nucleic acids structures in Nidovirales and in novel SARS-CoV-2 are therefore promising druggable structures that can be targeted and utilized in the future.

The small size of RNA virus genomes is in principle linked to their limited ability for RNA synthesis, which is directly connected to a replication complex containing RNA-dependent RNA polymerase (RdRp) without reparation mechanisms (Gorbalenya et al., 2006). These RNA viruses can generate one mutation per genome per replication round (Drake and Holland, 1999). This combination of features means that RNA viruses are able to adapt to new environmental conditions, but they are limited in expanding their genomes because they must keep their mutation load low so that their survival is possible (Eigen and Schuster, 1977;Lauring et al., 2013;Carrasco-Hernandez et al., 2017). Nidovirales bind to their host cell receptors on the cell surface, after which fusion of the viral and cellular membranes is mediated by one of the surface glycoproteins. This mechanism, which progresses in the cytoplasm or endosomal membrane, releases the nucleocapsid into the host cell's cytoplasm. After genome "transportation, " translation of two replicase open reading frames (ORFs) is initiated by host ribosomes. This results in large polyprotein precursors that undergo autoproteolysis to eventually produce a membrane-bound replicase/transcriptase complex. This complex initiates synthesis of the genome RNA and controls the synthesis of structural and some other proteins. New virus particles are assembled by association of the new genomes with the cytoplasmic nucleocapsid protein and subsequent envelopment of the nucleocapsid structure. Subsequently, the viral envelope proteins are inserted into intracellular membranes and targeted to the site of virus assembly (most often membranes between the endoplasmic reticulum and Golgi complex) and then they meet up with the nucleocapsid and trigger the budding of virus particles into the lumen of the membrane compartment. The newly formed virions then leave the cell by following the exocytic pathway toward the plasma membrane (Lai and Cavanagh, 1997;Gorbalenya, 2001;Snijder et al., 2003;Ziebuhr, 2004;Gorbalenya et al., 2006).
Today, SARS-CoV-2 is being studied intensively by scientists all over the world due to the ongoing 2020 coronavirus pandemic (Cohen, 2020). The origin of this virus is unknown, but recently a two-hit scenario was proposed wherein SARS-CoV-2 ancestors in bats first acquired genetic characteristics of SARS-CoV by incorporating a SARS-like receptor-binding domain through recombination prior to 2009; subsequently, the lineage that led to SARS-CoV-2 accumulated further, unique changes specifically within this domain (Patino-Galindo et al., 2020). As true of SARS-CoV, cell entry by SARS-CoV-2 depends upon angiotensin-converting enzyme 2 (ACE2) and transmembrane protease, serine 2 (TMPRSS2) for viral spike protein priming (Hoffmann et al., 2020). The genome of SARS-CoV-2 was sequenced and annotated in early January 2020 . A recent study revealed that the transcriptome of SARS-CoV-2 is highly complex due to numerous canonical and non-canonical recombination events. Moreover, it was found that SARS-CoV-2 produces transcripts encoding unknown ORFs and at least 41 potential RNA modification sites with an AAGAA motif were discovered in its RNA (Kim et al., 2020).
G-quadruplex binding proteins (QBPs) play crucial roles in many signaling pathways, including such biologically highly relevant activities as cell division, dysregulations of which lead to cancer development (Wu et al., 2008;Brázda et al., 2014). QBPs have been found to be involved in various viral infection pathways. An interesting example is HIV-1 nucleocapsid protein NCp7. It has been described how this protein helps to resolve an otherwise very stable G-quadruplex structure in viral RNA that stalls reverse transcription (Butovskaya et al., 2019). Various Gquadruplex-forming aptamers are used as drugs against many different viral proteins, suggesting a prominent role for QBPmediated regulation. Among many other viruses, in particular Hepatitis C virus, HIV-1, and SARS-CoV are targets of these Gquadruplex-forming aptamers (Platella et al., 2017). Quadruplex binding domain has been found in non-structural protein 3 (Nsp3) (Lei et al., 2018). This so-called SARS-UNIQUE-DOMAIN (SUD), and especially its M subdomain, was observed to be essential for SARS replication in host cells. Deletion or even substitution mutations in key RNA-interacting amino acids were shown to result in viral inability to replicate within host cells (Kusov et al., 2015). Moreover, this subdomain was found also in MERS and several other coronaviruses (Kusov et al., 2015).
G-quadruplexes are secondary nucleic acid structures formed in guanine-rich strands (Burge et al., 2006;Vorlíčková et al., 2012;Kolesnikova and Curtis, 2019). These have been detected in various genomes, but most extensively they have been described in human genomes (Chambers et al., 2015;Bedrat et al., 2016;Hänsel-Hertsch et al., 2016). They are present also in viruses (Lavezzo et al., 2018;Frasson et al., 2019). G-quadruplexes probably play an important role in regulating replication in most viral nucleic acids (Lavezzo et al., 2018), and these structures have been suggested as targets for antiviral therapy (Métifiot et al., 2014;Richter, 2018, 2020). Along with cruciforms and hairpins, which can be formed in nucleic acids by inverted repeats (IRs), G-quadruplexes are significant genome elements playing specific biological roles. They are involved, for example, in the regulation of DNA replication and transcription (Bagga et al., 1990;Yu, 2009;Brázda et al., 2011). It has been demonstrated that IRs are important for various processes in viruses, including their genome organization (Li and Li, 2010;Ishimaru et al., 2013;Xie et al., 2017;Bridges et al., 2019). Another interesting RNA motif that has been used as a drug target and was found in SARS-CoV targeted by 1,4-diazepame is the "slippery sequence" followed by a pseudoknot (Plant et al., 2005). This structure, common among all coronaviruses, works based on ribosomal −1 frameshifting that switches on viral fusion proteins (Plant et al., 2005).
In all prokaryotic and eukaryotic cells, as well as viruses, there have been found sequence motifs such as IR sequences forming cruciforms and hairpins or G-rich sequences that form G-quadruplexes (Brázda et al., 2011;Lavezzo et al., 2018;Bartas et al., 2019). In the present research, we conducted a systematic and comprehensive bioinformatic study searching for the occurrence of IRs and putative quadruplex sites (PQSs) within the genomes of all known Nidovirales. The aim was to find one or more potential druggable RNA targets to address the present COVID-19 threat (Figure 1).

Genome Source
Linear genomes of 109 viruses from the order Nidovirales were downloaded from the genome database of the National Center for Biotechnology Information (NCBI) (Federhen, 2011). Full names, phylogenetic groups, exact NCBI accession numbers, and further information are summarized in Supplementary Material 1.

Nidovirales Phylogenetic Tree Construction
Exact taxonomic identifiers of all analyzed Nidovirales species (obtained from Taxonomy Browser via NCBI Taxonomy Database Drake and Holland, 1999 were downloaded to phyloT: a tree generator (http://phylot.biobyte.de) and a phylogenetic tree was constructed using the function "Visualize in iTOL" in the Interactive Tree of Life environment (Letunic and Bork, 2019).

Analyses of PQSs Occurrence in Nidovirales Sequences
Nidovirales sequences were analyzed using our G4Hunter Web Tool , and selected sequences were verified also by QGRS mapper (Kikin et al., 2006) and G4screener (Garant et al., 2018). G4Hunter web is a more recent tool compared to Quadparser and QGRS mapper, and this algorithm allows quantitative analyses whereby G4 propensity is calculated depending on G richness and G skewness. Moreover, a new implementation of G4Hunter allows for performing batch analyses. In general, the earlier algorithms did not consider atypical quadruplex-forming structures, as have been described in various outstanding articles (Guedin et al., 2010;Kocman and Plavec, 2017;Lightfoot et al., 2019). The G4Hunter Web software is capable to read the NCBI identifier of the sequences uploaded in a.csv file. The parameters for G4Hunter were set to 25 as window size and G4Hunter score above 1.2. The results for each analyzed sequence contained information about the size of the genome and number of putative PQSs. All the results were merged into a single Microsoft Excel file where statistical analysis was then made. We also downloaded features tables of each virus from the NCBI database. These tables contain information about known features in the genome of each species. We searched the occurrence of G-quadruplex-forming sequences before, inside, and after the specific features of each genome using a predefined feature neighborhood ±100 nucleotides (nt). Data were then plotted in Excel, where the statistical analysis was made. Complete analyses of PQS occurrence in Nidovirales, including a list of those PQSs found, are provided in Supplementary Material 2.

Analyses of IRs Occurrence in Nidovirales Sequences
software was modified to read NCBI identifiers of sequences from text files. The size of IRs was set to 6-30 nt, size of spacers to 0-10 nt, and maximally one mismatch was allowed. A separate list of IRs in each of the 109 sequences and an overall report were exported. The overall report contained the lengths of the analyzed sequences, total number of IRs found, numbers of IRs grouped by size of IR (6-30 nt), and sum of IRs longer than 8, 10, and 12 nt. The software also counted frequency of IRs in each sequence. Frequencies of IRs were normalized per 1,000 nt. Features tables of 109 Nidovirales genomes were downloaded from the NCBI database and grouped by their names as stated in the feature table file. Analyses of IRs occurrence inside and around (before and after) these features was performed. The search for IRs took place in predefined feature neighborhoods ±100 nt around and inside feature boundaries. We calculated the numbers of all IRs and of those longer than 8, 10, and 12 nt in regions before, inside, and after the features. The categorization of an IR according to its overlap with a feature or feature neighborhood is demonstrated by the example shown in Supplementary Material 3. Complete analyses of IRs occurrence in Nidovirales are provided in Supplementary Material 4.
FIGURE 2 | Phylogenetic relationships in Nidovirales. SARS-CoV-2 is not yet recognized by iToL, but, based upon current evidence, it probably will be placed close to the Bat Hp-betacoronavirus/Zheijang2013 and Severe acute respiratory syndrome-related coronavirus (SARS-nCoV) branches (Lu et al., 2020).

RNA Fold Predictions
In order to be able to display higher structures of the coronavirus genome, we used Galaxy's free-online webserver (Afgan et al., 2018) and its RNA fold tool (Lorenz et al., 2011). This tool allows quick calculation of minimum free energy of secondary structures. We left the default parameters (Temperature 37 • C, Unpaired bases to participate in all dangling ends, Naview layout). We used SARS-CoV-2 genomic RNA sequence (NC_045512.2) in FASTA format as the input format. The output data were then displayed using the RNA plot tool (Lorenz et al., 2011). We again left the default parameters (Naview layout, Output format Postscript.ps). We then downloaded the displayed secondary structures in high-resolution format. The raw data are provided in Supplementary Material 5.

Prediction of Human RNA-Binding Proteins Sites in SARS-CoV-2 RNA
The human SARS-CoV-2 RNA sequence was downloaded from NCBI (accession NC_045512.2) in FASTA format and inserted into the RBPmap (Mapping Binding Sites of RNA-binding proteins) web-based tool (Paz et al., 2014). The database of 114 human experimentally validated motifs was used. Both the "High stringency" and "Apply conservation filter" options were used. The output was further filtered in Excel to keep only those hits below p = 1.10 −6 . The complete results are provided in Supplementary Material 6.

Statistical Analyses
A cluster dendrogram of PQS characteristics was constructed in the program R, version 3.6.3, using the pvclust package to further reveal and graphically depict similarities between particular Nidovirales species (Figure 4). The following values were used as input data: frequency of PQS per 1,000 nt with threshold G4Hunter score of 1.2; frequency of PQS per 1,000 nt with threshold G4Hunter score of 1.4; and % PQS in genome (coverage) (Supplementary Material 7). A cluster dendrogram of IRs (Figure 6) was constructed from these values: frequency of IRs per 1,000 nt; frequency of IRs per 1,000 nt of length 8+; frequency of IRs per 1,000 nt of length 10+; and frequency of IRs per 1,000 nt of length 12+. The following parameters were used for both analyses: cluster method "ward.D2, " distance "euclidean, " number of bootstrap resampling was 10,000. Statistically significant clusters (based on AU values above 95, equivalent to p < 0.05) are highlighted by rectangles marked with broken red lines. R code is provided in Supplementary Material 7.
To determine whether SARS-CoV-2 significantly differs in frequencies of PQSs and IRs compared to randomly shuffled sequences (N = 10) of the same length and nucleotide composition, we performed a two-sided Wilcoxon signed-rank test (with p-value cutoff at 0.05). This test was run for plus and minus strands separately.

Phylogenetic Relationships in Nidovirales
According to the current state of knowledge, the order Nidovirales can be divided into six distinct families (Liò and Goldman, 2004;Hanada et al., 2005;Chen et al., 2020). There are just two unclassified species among all those species recorded within NCBI Genome (Figure 2). Coronaviridae is the largest family and consists of 56 species. Among these are 12 species able Seq (total number of sequences), Median (median length of sequences), Short (shortest sequence), Long (longest sequence), PQS (total number of predicted PQSs), Mean f (mean frequency of predicted PQS per 1,000 nt), Min f (lowest frequency of predicted PQS per 1,000 nt), Max f (highest frequency of predicted PQS per 1,000 nt), GC% (mean GC content). Note that 132 PQSs were in Nidovirales with Undefined host.
to infect humans, including SARS-CoV, MERS-CoV, and SARS-CoV-2. Second largest is the Arteriviridae family with 22 species, third largest is the Mesnidovirineae family with 13 species, fourth is the Tobaniviridae family with 11 species, including 1 species able to infect humans. Fifth largest is the Ronidovirineae family with 5 species, and sixth is the Mononiviridae family with only 1 species.

Variation in PQS Frequency in Nidovirales
We analyzed the occurrence of PQSs using G4Hunter in all 109 known genomes of Nidovirales. , and there were no PQSs above the 1.6 G4Hunter score threshold. In general, a higher G4Hunter score means a higher probability of G-quadruplexes forming inside the PQS (Bedrat et al., 2016). Genomic sequence sizes, total PQS counts, and PQS frequencies characteristics are summarized in Table 1.
The highest PQS frequencies were found in Beihai Nido-like virus 1, which is an intracellular parasite of sea snails in genus Turritella. Beihai Nido-like virus 1 has in total 184 PQSs in its genomic sequence 20,278 nt long and PQS frequency of 9.07 PQS per 1,000 nt. On the other hand, no PQSs were found in the genomic sequences of 16 other Nidovirales species. The mean PQS value for all Nidovirales was 0.48 PQS per 1,000 nt. By viral families, the highest mean PQS frequency per 1,000 nt was in Arteriviridae (1.18), followed by in Mesnidovirineae (1.03), much lower in Tobaniviridae (0.43) and in Ronidovirineae (0.39), and the lowest was in Coronaviridae (0.10). When grouped and analyzed by host organisms, new information became apparent. The highest mean PQS frequency was in Nidovirales infecting Invertebrates (0.80), then in Nidovirales infecting Vertebrates including Humans (0.38), and the last in Nidovirales infecting Humans (0.06). In pathogenic human coronaviruses, the total PQSs counts were as follow: 4 PQS in SARS-CoV, 0 PQS in MERS-CoV, and 1 PQS in SARS-CoV-2. To test whether the PQS frequency per 1,000 nt in SARS-CoV-2 is significantly different than in a randomly shuffled sequence, we generated 10 scrambled sequences with length and nucleotide content the same as in SARS-CoV-2. In that test, the mean PQS frequency per 1,000 nt was 0.191 and standard deviation was 0.101. We performed Wilcoxon signed-rank test with continuity correction and the resulting p-value was 0.008 Thus, the PQS frequency per 1,000 nt in SARS-CoV-2 (0.033) is significantly lower than expected (only about one-sixth).
All PQSs found in ranges of G4Hunter score intervals and precomputed PQS frequencies per 1,000 nt are summarized in Table 2. The relationships between observed PQS frequency per 1,000 nucleotides and GC content in all analyzed Nidovirales sequences are depicted in Figure 3. The Beihai Nido-like virus has both the highest GC content and highest PQS frequency. Within Nidovirales with Humans as host, however, the highest PQS frequency was found in Breda virus, which does not have the highest GC content in the dataset.
Most of the PQSs have G4Hunter scores between 1.2 and 1.4. Only a few sequences have G4Hunter scores above 1.4. In comparison with other Nidovirales, the members of the Coronaviridae and especially pathogenic human Coronaviridae have the lowest PQS frequency in the dataset.
A cluster dendrogram based on the PQS characteristics (Figure 4) further revealed interesting information. Beihai Nido-like virus 1 was confirmed as an outlier from the rest of the Nidovirales. All viral species belonging to the largest Coronaviridae family are located together in the middle of the cluster dendrogram. The relatively abundant family Arteriviridae clustered together mainly on the right side of the cluster dendrogram. Based FIGURE 3 | Relationship between observed PQS frequency per 1,000 nucleotides and GC content in all analyzed Nidovirales sequences. In each G4Hunter score interval miniplot, frequencies were normalized according to the highest observed PQS frequency. Organisms with maximum frequency per 1,000 nt >50% are described and highlighted in color. on dendrogram distances, we highlighted two main branches-one consisting mainly of Coronaviridae (light blue background) and the second of Arteriviridae (light orange background). Figure 5 shows the differences in PQS frequency by annotated loci. We downloaded features for every virus genome and analyzed the presence of PQS in each annotated sequence and in its proximity (before and after). The most notable enrichment of PQSs was located inside 5 ′ UTR and before and inside 3 ′ UTR. The lowest PQS frequencies were found after 3 ′ UTR, before 5 ′ UTR, and after miscellaneous features.

Variation in Frequency of Inverted Repeats in Nidovirales
To find short IRs in Nidovirales genomic sequences, we utilized the core of the Palindrome analyzer (Brázda et al., 2016). We used values for Palindrome analyzer returning inverted repeats capable to form stable hairpins (size of IRs was set to 6-30 nt, size of spacer to 0-10 nt, maximally one mismatch). Genomic sequence sizes, total IR counts, and IR frequencies characteristics are summarized in Table 3. In total, 93,369 IRs were found and the mean IRs frequency was 33.90 per 1,000 nt. The maximal IR frequency was found in Planarian secretory cell nidovirus (56.17), FIGURE 5 | Differences in PQS frequency by annotated locus. The figure shows the PQS frequencies between annotations from the NCBI database. We analyzed frequencies of all PQSs inside, before, and after the annotations. Seq (total number of sequences), Median (median length of sequences), Short (shortest sequence), Long (longest sequence), GC% (mean GC content), IRs (total number of predicted IRs), Mean f (mean frequency of predicted IRs per 1,000 nt), Min f (lowest frequency of predicted IRs per 1,000 nt), and Max f (highest frequency of predicted PQS per 1,000 nt) are shown for all 109 Nidovirales sequences, families, host groups, and for SARS-CoV-2 alone.
the species with the largest genome among Nidovirales. It has been suggested that Planarian secretory cell nidovirus diverged early from multi-ORF Nidovirales and acquired additional genes, including those typical of large DNA viruses or hosts (RNAse T2, Ankyrin, and Fibronectin type II), which might modulate virushost interactions (Saberi et al., 2018). The lowest IR frequency was noticed in Beihai Nido-like virus 1 (19.23). Noteworthy is that this is the species of Nidovirales with the highest GC content and PQS frequency. A cluster dendrogram based on the IR characteristics (Figure 6) further revealed significant differences in IR presence within Nidovirales genomes. Based on dendrogram distances, we highlighted two main branches-one consisting mainly of Coronaviridae (light blue background) and the second of Arteriviridae (light orange background). It is noteworthy that SARS-CoV and MERS-CoV are clustered adjacent to one another (green asterisk), but the novel coronavirus SARS-CoV-2 is located relatively far away from them (red asterisk).
Differences in IR frequency by annotated locus are depicted in Figure 7. By families, the highest mean IRs frequency was found in Coronaviridae (37.08) and lowest in Ronidovirineae (26.37). Novel human coronavirus SARS-CoV-2 has relatively high IR frequency per 1,000 nt in comparison with other Nidovirales. To test if the frequency of IRs per 1,000 nt in SARS-CoV-2 is significantly greater than in randomly shuffled sequences, we generated 10 scrambled sequences with length and nucleotide content the same as in SARS-CoV-2 (for both positive and negative strands). For both positive and negative strands, the results were very similar (mean frequencies of IRs per 1,000 nt were 34.43 ± 0.54 and 34.47 FIGURE 7 | Differences in IR frequency by annotated locus. The chart compares IR frequencies per 1,000 nt between "gene" annotation and other annotated locations from the NCBI database that were found in genomes more than 10 times. We analyzed frequencies of all IRs (all) and of IRs with lengths 8 nt and longer (8+), 10 nt and longer (10+,) and 12 nt and longer (12+) within annotated locations (inside) as well as before and after annotated locations.
± 0.61, respectively). Wilcoxon signed-rank test with continuity correction showed the IR frequency per 1,000 nt in SARS-CoV-2 (40.23) to be significantly higher than expected (pvalue 0.003). When we inspected differences of IR frequencies according to hosts, we found the highest mean IR frequency per 1,000 nt to be in Nidovirales infecting Humans (38.13), then in Vertebrates (34.25), and the lowest mean frequency is in Invertebrates (31.77).
A summary of all IRs found in ranges of different IR sizes and precomputed IR frequencies per 1,000 nt is shown in Table 4. Although generally the frequency of IR presence decreases with the IR length, there are notable differences between groups and also between viruses with different hosts. The most as well as longest IRs occur in the Coronaviridae group and in viruses having humans as a host. IRs 12 bp long and longer are very rare in the Ronidovirineae group. The relationship between IRs length and stability of resulting secondary structure is not simple. While some authors believe that longer IRs are more stable, others suggest that there is an energy optimum defined by arm and spacer length (Sinden et al., 1991;Brázda et al., 2016;Georgakopoulos-Soares et al., 2018).
Differences in IR frequency according to annotated loci are shown in Figure 5. The most notable enrichment of IRs was found inside 5 ′ UTR for IRs of size 12+ nt, and this is the most frequently occurring location of 12+ IRs in Nidovirales genomes. Noteworthy is that there are no 12+ IRs around 5 ′ UTR loci, but there is an abundance of IRs 10+ nt long in these locations. The 5 ′ UTR are abundant for 12+ IR, but there are no 12+ IRs within 3 ′ UTR. This points to functional relevance of these IRs in viral genomes.

Prediction of Human SARS-CoV-2 RNA Folding
It is very likely that coronavirus RNA in the viral capsid is compactly folded into a spiral-like structure due to nucleocapsid phosphoprotein N, as was proposed earlier for SARS-CoV (Chang et al., 2014). But what structure does the RNA takes after the capsid content is released into the cytoplasm of the host cell? We made a global prediction of RNA folding for SARS-CoV-2 and in a negative control (random RNA sequences with the same GC content and length). It was clearly visible that there arose a much more compact RNA structure in SARS-CoV-2 (Figure 8). The RNA folding algorithm shows that inverted repeats forming hairpins have the main importance for its folding. Due to the relative lack of PQSs with high G4Hunter scores in the SARS-CoV-2 positive strand, we assume that G-quadruplexes are not important for its folding. For the 10 longest IRs obtained from Palindrome analyzer, we prepared local predictions of their 2D structures using the Mfold webserver (Zuker, 2003) (Supplementary Material 8).

Comparison of SUD Domain M Region of Nsp3 in Three Pathogenic Human Coronaviruses
Consisting of nearly 2,000 amino acids, Nsp3 is the largest multi-domain protein in Nidovirales. Nevertheless, Nsp3's role is still largely unknown. It is believed that the protein plays various roles in coronavirus infection. Nsp3 interacts with other viral Nsps as well as RNA to form the replication/transcription complex. It also acts on posttranslational modifications of host proteins to antagonize the host innate immune response. On the other hand, Nsp3 is itself modified in host cells by Nglycosylation and can interact with host proteins to support virus survival (Lei et al., 2018). Both SARS-CoVs (2003 and have PQSs in their genomes and have a retained M region of the SUD domain in protein Nsp3 that is reported to be critical for interacting with G-quadruplexes, and particularly all G-quadruplex interacting residues, as proposed by Kusov et al. (2015), are 100% conserved (Figure 9). MERS has no PQS in its RNA sequence and, strikingly, it has a deletion in this region of the SUD domain suggesting parallel evolution simplifications. Shown are data for all 109 representative Nidovirales sequences, for families, host groups, and for SARS-CoV-2 alone. Frequency was computed using total number of IRs in each category divided by total length of all analyzed sequences and multiplied by 1,000.
FIGURE 8 | RNA fold (Lorenz et al., 2011) prediction for SARS-CoV-2 RNA molecule (Left) and random sequence negative control (Right). This figure shows a high level of SARS-CoV-2 genome folding via complementarity of particular RNA regions and forming of hairpins and/or cruciforms. RNA fold prediction was carried out using default parameters via Galaxy webserver (Afgan et al., 2018), which enables queries of lengths longer than 10,000 nucleotides.
FIGURE 9 | Comparison of SUD domain M region in three pathogenic human coronaviruses. Nsp3 proteins of three pathogenic human viruses (SARS-CoV-2, SARS-CoV, and MERS-CoV) were aligned and part of the critical G-quadruplex binding M region of the SUD domain (525-577 aa according to the SARS-CoV reference sequence) was visualized. As predicted according to Kusov et al. (2015), the G-quadruplex-interacting residues K565, K568, and E571 are highlighted in green (upper consensus panel).

Prediction of Human RNA-Binding Proteins Sites in SARS-CoV-2 RNA
Using the RBPmap tool, we predicted human RNA-binding proteins sites in SARS-CoV-2 RNA. While holding to very stringent thresholds, three highly promising candidate human RNA-binding proteins were predicted. SRSF7 was predicted to bind viral RNA in nucleotide position 26,199 (NC_045512.2), which is in fact exactly the binding motif for this protein (ACGACG). Protein HNRNPA1 was predicted to bind viral RNA in nucleotide position 23,090-23 097 (exact binding motif GUAGUAGU). The last protein found was TRA2A in nucleotide position 3,056-3,065 and position 3,074-3,083, in both cases with one mismatch (the motifs found were GAAGAAGAAG and the experimentally validated motif for TRA2A is GAAGAGGAAG). All three of these proteins share multiple RGG-rich novel interesting quadruplex interaction (NIQI) motifs, which are common to most human G-quadruplexbinding proteins (Brázda et al., 2018;Huang et al., 2018;Masuzawa and Oyoshi, 2020). All find individual motif occurrence (FIMO) hits below p = 1.10 −4 are enclosed in Supplementary Material 9. The proteins SRSF7, HNRNPA1, and TRA2A are involved in mRNA splicing via the spliceosome pathway (Szklarczyk et al., 2015).

DISCUSSION
Local DNA structures have been shown to play important roles in basic biological processes, including replication and transcription (Brázda et al., 2011;Travers and Muskhelishvili, 2015;Surovtsev and Jacobs-Wagner, 2018). PQS analyses of human viruses have clearly demonstrated that these sequences are conserved also in the viral genomes (Lavezzo et al., 2018) and could be targets for antiviral therapies (Wang et al., 2015;Krafčíková et al., 2017;Ruggiero and Richter, 2018). For example, the conserved PQS sequence present in the L gene of Zaire ebolavirus and related to its replication is inhibited by interaction with G-quadruplex ligand TMPyP4, and this finding has led to suggestions that Gquadruplex RNA stabilization could constitute a new strategy against Ebola virus disease (Wang et al., 2016). A similar strategy has been proposed against Hepatitis C virus belonging to the Flaviviridae family of positive ssRNA viruses. Stabilization of G4s by Phen-DC3 ligand has been shown to lead to inhibited replication of this virus (Jaubert et al., 2018). Viruses with singlestranded genomes are often non-symmetrically distributed with guanines and cytosines. For example, Zika viral genome (Flaviviridae family) has a G-rich positive-sense genome and a C-rich negative-sense strand. Fleming et al. found more than 60 PQSs in the Zika virus genome on the positive strand but no PQSs on the negative strand. This observation identifies a large asymmetry with respect to PQS content between the two strands. The strand asymmetry for PQS sites likely results from the high guanine content relative to cytosine content in the positive-sense strands (Fleming et al., 2016). There are many PQS searching algorithms (reviewed by Puig Lombardi and Londoño-Vallejo, 2020). We used the G4Hunter algorithm.
Because it had originally been written to analyze DNA (Bedrat et al., 2016), the G4Hunter algorithm analyzed both strands simultaneously . Analyses of both strands are therefore presented in our manuscript. The results with positive G4Hunter scores show PQSs in the plus strand and the results with negative G4Hunter scores show PQSs in the minus strand. In our study, we analyzed all accessible genomes of Nidovirales, thus including also the contemporary pandemic SARS-CoV-2 genome. Our analyses found only one PQS in the genomic sequence of SARS-CoV-2, that having a G4Hunter score of −1.24 and the following sequences: 5 ′ -CCCCAAAAUCAGC GAAAUGCACCCC-3 ′ for a positive-strand intermediate and 5 ′ -GGGGUGCAUUUCGCUGAUUUUGGGG-3 ′ for a negativestrand intermediate (where G-quadruplex could theoretically arise). This sequence was not identified by other prediction algorithms that use searching based upon regex sequence or that are pattern-based (PQSfinder or QGRS mapper), and it has been given only a low score by the new machine-learning G4screener algorithm (Garant et al., 2018). Nevertheless, none of these algorithms consider possible substitutions of guanine by adenine in quadruplex tetrads, as has been described by Kocman and Plavec (2017), or other untypical quadruplexes as reviewed by Lightfoot et al. (2019) and stable quadruplexes with long loops that have been described (Guedin et al., 2010). Moreover, there is another GGG track in the proximity of our predicted sequence. On the other hand, analysis of the SARS-CoV-2 genome by QGRS algorithms showed 25 hits on the positive and 12 hits on the negative RNA strand (Supplementary Material 10). These hits are almost exclusively with two G-repeats only and have relatively low scores (maximum 19 for the QGRS algorithm and 1.111 for the G4Hunter algorithm). The best PQS suggested by the G4Hunter algorithm is located in the position 28,289-28,313 within the nucleocapsid phosphoprotein coding sequence. It is noteworthy that all viral RNAs are produced through negativestrand intermediates, which are only about 1% as abundant as their positive-sense counterparts (Fehr and Perlman, 2015). It would be of great value to know whether TMPyP4 or other G-quadruplex stabilizing compounds can inhibit replication processes of SARS-CoV-2. We have aligned three pathogenic human coronaviruses (SARS-CoV, SARS-CoV-2, and MERS-CoV) to see if there are differences in the SUD domain, which was earlier proposed to be G-quadruplex binding (Kusov et al., 2015). Comparison of three key amino acid residues involved in G-quadruplex binding revealed that the SUD domain of MERS-CoV lacks these residues. This correlates with the fact that no G-quadruplex was predicted that is within its genome. It has been demonstrated, however, that the RGG domain can play roles in various nucleic acid and protein interactions (Thandapani et al., 2013), so this correlation can be G4independent for SARS-CoV-2. RNA hairpins, which are formed by IRs, are basic structural elements of RNA and play crucial roles in gene expression and intermolecular recognition. Conserved palindromic RNA structures have been found in many viral genomes, including HIV-1, and play a crucial role in their replication (Liu et al., 2018). We have found an abundance of IRs inside 5 ′ UTR in Nidovirales genomes. In general, 5 ′ UTR is an important locus for regulation of viral replication and gene expression. It has been demonstrated that stem integrity of phylogenetically conserved stem-loop structure located in 5 ′ UTR of the PRRSV virus from the Arteriviridae family is crucial for viral replication and subgenomic mRNA synthesis. Similar secondary structures have been proposed for several viruses from the Arteriviridae and Coronaviridae families (Lu et al., 2011). Our analyses of annotated features in Figure 5 therefore support this report. The discrete locations of specific IRs in viral genomes could therefore be additional targets for their regulation. Significant differences of PQS and IR frequencies among various ssRNA viruses in Nidovirales groups show that their genome organization and regulation are not identical but that for some Nidovirales the presence of the G-quadruplexes most probably does not play an essential role in their biological regulation. Moreover, our analyses suggest that G-quadruplexes have been evolutionally eliminated in some genomes of Nidovirales. This is quite surprising, considering that G-quadruplexes have been found in all evolutionary groups, including such CG low-level organisms as Saccharomyces cerevisiae Brázda et al., 2019;Singh and Lakhanpaul, 2019;Gage and Merrick, 2020). In fact, it could be an evolutionary advantage not to present G-quadruplex in viral genome because a number of cellular proteins interact with G-quadruplexes (Brázda et al., 2014;Mishra et al., 2016) and therefore the presence of G-quadruplex in viral genome could serve as a structure recognized by the innate immune system (Unterholzner et al., 2010;Hároníková et al., 2016;Voter et al., 2020). Their RNAs are therefore not recognized as alien nucleic acids and are processed by the cellular machinery. On the other hand, presence of IRs in Nidovirales genomes constitutes an inseparable part of their genomes and allows their correct folding and structure-specific regulation of their functions (Lorenz et al., 2011;Dutkiewicz et al., 2016).
Targeting viral proteins is usually effective only against specific viral strains and fails even for closely related viral species. It appears that targeting host proteins should be able to provide a response toward a wider spectrum of viruses inasmuch as different viruses exploit common cellular pathways. Many cellular RNA-binding proteins (RBPs) containing wellestablished RNA-binding domains (RBDs) are known to be critical for infection by different viruses. Recently, 472 RBPs were reviewed for their linkage to viruses (Garcia-Moreno et al., 2018). It has been demonstrated that G-quadruplex formation in HIV-1 viral genome stalls RNA polymerase, thereby limiting viral replication in host cells. HIV-1 nucleocapsid protein NCp7 helps to resolve G-quadruplex formation and therefore enables virus to spread. Stabilization of this quadruplex has been targeted by several experimental compounds. This treatment slowed or inhibited viral growth (Butovskaya et al., 2019).
Virus reproduction is dependent on the cellular transcription machinery, and therefore the interaction of the cellular proteins with viral RNA could be another target for antiviral therapy (Roberts et al., 2009). Our analysis predicted several QBPs to be capable to bind SARS-CoV-2 RNA. One of these, HNRNPA1 protein, is responsible for nuclear-cytoplasm shuttling (Garcia-Moreno et al., 2018). Like some other RNA-binding proteins, HNRNPA1 forms so-called membraneless organelles. These organelles are assemblies of proteins along with RNA or DNA that condense in specific cellular loci. The organelles can undergo transition from liquid-like droplets to amyloid fibrils, and mutations in so-called low complexity domains of these proteins lead to formation of amyloid aggregations that are found in many neurodegenerative diseases (Gui et al., 2019). HNRNPA1 is involved in many different cellular processes and has been targeted in various diseases. One such example is the varicella zoster virus that causes chicken pox. Moreover, this response to varicella zoster virus has been connected with autoimmune disease that complicates multiple sclerosis (Kattimani and Veerappa, 2018). HNRNPA1 is known to bind and resolve G-quadruplex formed in TRA2B promoter and to promote its transcription. Dysregulation of this binding leads to progression of colon cancer. This interaction has been targeted by the well-known G-quadruplex stabilizer pyridostatin, which led to decreased transcription from the TRA2B promoter (Nishikawa et al., 2019). Furthermore, HNRNPA1 is co-expressed with bromodomain and extraterminal domain protein BRD4 in human tumor samples. It has been shown experimentally that the well-described, naturally occurring polyphenolic flavonoid quercetin inhibited this protein and thereby led to better susceptibility to treatment in cancer patients (Pham et al., 2019). Both SRSF7 and TRA2A proteins, which are predicted to interact with SARS-CoV-2 RNA, play roles in alternative splicing (Ghosh et al., 2016). SRSF7 is a serine and arginine-rich splicing factor and is part of the spliceosome. Its expression has been connected to several types of cancer and it has been shown that its knockdown induced p21 expression and thus reduced cancer development (Saijo et al., 2016). Little is known about TRA2A protein. It was first identified in insects together with its paralog TRA2B, which has been studied to a greater extent . It has been found that TRA2A can promote paclitaxel therapy and promote cancer progression in triple-negative breast cancers (Liu et al., 2017). A role of TRA2A protein in regulation of HIV1 virus replication has been described. Both TRA2A and TRA2B bind to a specific HIV1 sequence and regulate its replication within the cell through alternative splicing of viral RNA (Erkelenz et al., 2013).
Scientists around the world are united in their efforts to find an effective therapy against coronavirus disease . Among the most promising candidates are remdesivir and chloroquine. Remdesivir is an adenosine analog that incorporates into nascent viral RNA chains and results in premature termination. Preliminary data have shown that remdesivir effectively inhibited virus infection in a human liver cancer cell line . Chloroquine is a potential broadspectrum antiviral drug and it already has been widely used as a low-cost and safe anti-malarial and autoimmune disease drug for more than 70 years. Application of chloroquine causes elevation of endosomal pH and also interferes with terminal glycosylation of the cellular receptor ACE2. This probably has a negative influence on virus-receptor binding and abrogates the infection (Vincent et al., 2005). Drug repurposing seems like a very good strategy for quickly and at low cost finding a new therapy for new human diseases (Oprea and Mestres, 2012). To date, there are few substances with G-quadruplex stabilizing features. One example is topotecan (also known by its brand name Hycamtin R ), which frequently is used for treating ovary cancer. If everything else fails, it seems to be an option. On the other hand, topotecan cause unpleasant side effects, such as nausea, vomiting, and diarrhea (Topotecan -Chemotherapy Drugs -Chemocare, n.d.). Many studies have described strong G-quadruplex stabilization effects (Li et al., 2018;Satpathi et al., 2018), which might be one possible mode of action. Moreover, it has been proven that higher structures of nucleic acids, and G-quadruplex especially, might be stabilized by use of various natural substances. Berbamine is one such substance and is a component of traditional Chinese medicine. It frequently is used for treating chronic myeloid leukemia or melanoma and has strong binding affinity to G-quadruplex structures (Tan et al., 2014). Viral nucleic acids and their loci with G-quadruplexforming potential are in all cases very specific and are promising molecular targets for treating serious diseases. Evidence of Gquadruplex formation as a potential target for therapy was proven for Hepatitis A, flu virus, and HIV-1 (Métifiot et al., 2014). Because all the aforementioned cases concern RNA viruses, the viral G-quadruplexes are localized in the cytoplasm of the host cell. Our comparative analyses of IRs and PQSs in Nidovirales show that, in contrast to IRs, which are presented abundantly in all Nidovirales genomes, the sequences able to form G-quadruplex structures are very unequally distributed and are very rare especially in Nidovirales species capable to infect humans. This suggests intentional suppression during evolution in order to simplify viral RNA replication. Finding the proper stabilizers of viral higher RNA structures might be crucial for inhibiting or stopping viral RNA replication in order to gain time for the immune system to deal successfully with an infection.

DATA AVAILABILITY STATEMENT
All datasets presented in this study are included in the article/Supplementary Material.

AUTHOR CONTRIBUTIONS
MB and VB contributed to the conceptualization, formal analysis, and methodology. NB, AC, AV, and TS helped with the data curation. PP was involved in with funding acquisition and project administration. MB, VB, AV, and JČ carried out the investigation. MB, NB, OP, VB, and EJ helped with the resources. VB, KM, and PP supervised the study. VB and EJ validated the study. MB, NB, and AC worked on the visualization. MB, VB, NB, OP, and JČ wrote the original draft. VB, JČ, KM, and PP reviewed and edited the manuscript. All authors contributed to the article and approved the submitted version.

FUNDING
This work was supported by the Ministry of Education, Youth and Sports of the Czech Republic in the National Feasibility Program I (LO1208 TEWEP); by the EU structural funding Operational Programme Research and Development for innovation, project No. CZ.1.05/2.1.00/19.0388; and by the projects SGS/01/PrF/2020 and SGS/07/PrF/2020 financed by University of Ostrava.