A disease named Covid-19 by the World Health Organization and caused by the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has been recognized as responsible for the pneumonia outbreak that started in December, 2019 in Wuhan City, Hubei, China 1 and spread in February to Milan, Lombardy, Italy 2 becoming pandemic. As of April 2020, the virus infected >2'000'000 people in >200 countries.
SARS-CoV-2 is a positive-sense single-stranded RNA virus that shares similarities with other betacoronavirus such as severe acute respiratory syndrome coronavirus (SARS-CoV) and Middle East respiratory syndrome coronavirus (MERS-CoV) 3 . Bats have been identified as the primary host for SARS-CoV and SARS-CoV-2 4,5 but the intermediate host linking SARS-CoV-2 to humans is still unknown, although a recent report indicates that pangolins could be involved 6 .
Coronaviruses use species-specific proteins to mediate the entry in the host cell and the spike S protein activates the infection in human respiratory epithelial cells in SARS-CoV, MERS-CoV and SARS-CoV-2 7 . Spike S is assembled as a trimer and contains around 1,300 amino acids within each unit 8, 9 . The receptor binding domain (RBD) of Spike S, which contains around 300 amino acids, mediates the binding with angiotensin-converting enzyme, (ACE2) attacking respiratory cells. Another region upstream of the RBD, present in MERS-CoV but not in SARS-CoV, is involved in the adhesion to sialic acid and could play a key role in regulating viral infection 7, 10 .
At present, very few molecular details are available on SARS-CoV-2 and its interactions with the human host, which are mediated by specific RNA elements 11 . To study the RNA structural content, we used CROSS 12 that was previously developed to investigate large transcripts such as the human immunodeficiency virus HIV-1 13 . CROSS predicts the structural profile of RNA molecules (single-
and double-stranded state) at single-nucleotide resolution using sequence information only. Here, we performed sequence and structural alignments among 62 SARS-CoV-2 strains and identified the conservation of specific elements in the spike S region, which provides clues on the evolution of domains involved in the binding to ACE2 and sialic acid.
As highly structured regions of RNA molecules have strong propensity to form stable contacts with proteins 14 and promote assembly of specific complexes 15, 16 , SARS-CoV-2 domains enriched in double-stranded content are expected to establish interactions within host cells that are important to replicate the virus 17 . To investigate the interactions of SARS-CoV-2 RNA with human proteins, we employed catRAPID 18, 19 . catRAPID 20 estimates the binding potential of a specific protein for an RNA molecule through van der Waals, hydrogen bonding and secondary structure propensities allowing identification of interaction partners with high confidence 21 . The unbiased analysis of more than 100000 protein interactions with SARS-CoV-2 RNA reveals that the 5' of SARS-CoV-2 has strong propensity to bind to human proteins involved in viral infection and especially reported to be associated with HIV infection. A comparison between SARS-CoV and HIV reveals indeed similarities 22 , but the relationship between SARS-CoV-2 and HIV is still unexplored. Interestingly, HIV and SARS-CoV-2, but not SARS-CoV nor MERS-CoV, have a furin-cleavage site occurring in the spike S protein, which could explain the spread velocity of SARS-CoV-2 compared to SARS-CoV and MERS-CoV 23, 24 . Yet, many processes related to SARS-CoV-2 replication are unknown and our study aims to suggest relevant protein interactions for further investigation.
We hope that our large-scale calculations of structural properties and binding partners of SARS-CoV-2 will be useful to identify the mechanisms of virus replication within the human host.
Structured elements within RNA molecules attract proteins 14 and reveal regions important for interactions with the host 25 . Indeed, each gene expressed from SARS-CoV-2 is preceded by conserved transcription-regulating sequences that act as signal for the transcription complex during the synthesis of the RNA minus strand to promote a strand transfer to the leader region to resume the synthesis. This process is named discontinuous extension of the minus strand and is a variant of similarity-assisted template switching that operates during viral RNA recombination 17 .
To analyze SARS-CoV-2 structure (reference Wuhan strain MN908947.3), we employed CROSS 12 to predict the double-and single-stranded content of RNA genomes such as HIV-1 13 . We found the highest density of double-stranded regions in the 5' (nucleotides 1-253), membrane M protein (nucleotides 26523-27191), spike S protein (nucleotides 23000-24000), and nucleocapsid N protein (nucleotides 2874-29533; Fig. 1 ) 26 . The lowest density of double-stranded regions were observed at nucleotides 6000-6250 and 20000-21500 and correspond to the regions between the nonstructural proteins nsp14 and nsp15 and the upstream region of the spike surface protein S (Fig. 1 ) 26 . In addition to the maximum corresponding to nucleotides 23000-24000, the structural content of spike S protein shows minima at around nucleotides 20500 and 24500 (Fig. 1) .
We used the Vienna method 27 to further investigate the RNA secondary structure of specific regions identified with CROSS 13 . Employing a 100 nucleotide window centered around CROSS maxima and minima, we found good match between CROSS scores and Vienna free energies ( Fig.   1 ). Strong agreement is also observed between CROSS and Vienna positional entropy, indicating that regions with the highest structural content have also the lowest structural diversity.
Our analysis suggests the presence of structural elements in SARS-CoV-2 that have evolved to interact with specific human proteins 11 . Our observation is based on the assumption that structured regions have an intrinsic propensity to recruit proteins 14 , which is supported by the fact that structured transcripts act as scaffolds for protein assembly 15, 16 .
We employed CROSSalign 13 to study the structural conservation of SARS-CoV-2 in different strains (Materials and Methods).
In our analysis, we compared the Wuhan strain MN908947.3 with around 2800 other coronaviruses (data from NCBI) with human ( Fig. 2) or other hosts (Supp. Fig. 1) . When comparing SARS-CoV-2 with human coronaviruses (1387 strains, including SARS-CoV and MERS-CoV), we found that the most conserved region falls inside the spike S genomic locus (Fig. 2) . More precisely, the conserved region is between nucleotides 23000 -24000 and exhibits an intricate and stable To better investigate the sequence conservation of SARS-CoV-2, we compared 62 strains isolated from different countries during the pandemic (including China, USA, Japan, Taiwan, India, Brazil, Sweden, and Australia; data from NCBI and in VIPR www.viprbrc.org; Materials and Methods).
Our analysis aims to determine the relationship between structural content and sequence conservation.
Using Clustal W for multiple sequence alignments 28 , we observed general conservation of the coding regions (Fig. 3A) . The 5' and 3' show high variability due to experimental procedures of the sequencing and are discarded in this analysis 29 . One highly conserved region is between nucleotides 23000 -24000 in the spike S genomic locus, while sequences up-and downstream are variable (red bars in Fig. 3A) . We then used CROSSalign 13 to compare the structural content (Materials and Methods). High variability of structure is observed for both the 5' and 3' and for nucleotides between 21000 -22000 as well as 24000 -25000, associated with the S region (red bars in Fig. 3A) . The rest of the regions are significantly conserved at a structural level (p-value < 0.0001; Fisher's test).
We then compared protein sequences coded by the spike S genomic locus (NCBI reference QHD43416) and found that both sequence (Fig. 3A) and structure (Fig. 2) of nucleotides 23000 -24000 are highly conserved. The region corresponds to amino acids 330-500 that contact the host receptor angiotensin-converting enzyme 2 (ACE2) 30 promoting infection and provoking lung injury 24, 31 . By contrast, the region upstream of the binding site receptor ACE2 and located in correspondence to the minimum of the structural profile at around nucleotides 22500-23000 ( Fig. 1) 
is highly variable 32 , as indicated by T-coffee multiple sequence alignments 32 (Fig. 3A) . This part of the spike S region corresponds to amino acids 243-302 that in MERS-CoV binds to sialic acids regulating infection through cell-cell membrane fusion ( Fig. 3B ; see related manuscript by E.
Milanetti et al.) 10, 33, 34 .
Our analysis suggests that the structural region between nucleotides 23000 and 24000 of Spike S region is conserved among coronaviruses (Fig. 2) and that the binding site for ACE2 has poor variation in human SARS-CoV-2 strains (Fig. 3B) . By contrast, the region upstream, which has propensity to bind sialic acids 10, 33, 34 , showed poor structural content and high variability (Fig. 3B) .
In order to obtain insights on how the virus replicates in human cells, we predicted SARS-CoV-2 interactions with the whole RNA-binding human proteome. Following a protocol to study structural conservation in viruses 13 , we first divided the Wuhan sequence in 30 fragments of 1000 nucleotides each moving from the 5' to 3' and then calculated the protein-RNA interactions of each fragment with catRAPID omics (3340 canonical and putative RNA-binding proteins, or RBPs, for a total 102000 interactions) 18 . Proteins such as Polypyrimidine tract-binding protein 1 PTBP1 (Uniprot P26599) showed the highest interaction propensity (or Z-score; Materials and Methods) at the 5' while others such as Heterogeneous nuclear ribonucleoprotein Q HNRNPQ (O60506) showed the highest interaction propensity at the 3', in agreement with previous studies on coronaviruses ( Fig.   4A ) 35 .
For each fragment, we predicted the most significant interactions by filtering according to the Z score. We used three different thresholds in ascending order of stringency: Z ³ 1.50, 1.75 and 2 respectively and we removed from the list the proteins that were predicted to interact promiscuously with more than one fragment. Fragment 1 corresponds to the 5' and is the most contacted by RBPs (around 120 with Z³2 high-confidence interactions; Fig. 4B ), which is in agreement with the observation that highly structured regions attract a large number of proteins 14 . Indeed, the 5' contains multiple stem loop structures that control RNA replication and transcription 36, 37 . By contrast, the 3' and fragment 23 (Spike S), which are still structured but to a lesser extent, attract fewer proteins (10 and 5, respectively), while fragment 20 (between orf1ab and Spike S) that is predicted to be unstructured, does not have binding partners.
The interactome of each fragment was then analysed using cleverGO, a tool for Gene Ontology (GO) enrichment analysis 38 . Proteins interacting with fragments 1, 2 and 29 were associated with annotations related to viral processes ( Fig. 4C ; Supp. Table 1 ). Considering the three thresholds applied (Materials and Methods), we found 22 viral proteins for fragment 1, 2 proteins for fragment 2 and 11 proteins for fragment 29 (Fig. 4D) .
Among the high-confidence interactors of fragment 1, we discovered RBPs involved in positive regulation of viral processes and viral genome replication, such as double-stranded RNA-specific editase 1 ADARB1 (Uniprot P78563 39 ), 2-5-oligoadenylate synthase 2 OAS2 (P29728) and 2-5Adependent ribonuclease RNASEL (Q05823). Interestingly, 2-5-oligoadenylate synthase 2 OAS2 has been reported to be upregulated in human alveolar adenocarcinoma (A549) cells infected with SARS-CoV-2 (log fold change of 4.2; p-value of 10 -9 and q-value of 10 -6 ) 40 . While double-stranded RNA-specific adenosine deaminase ADAR (P55265) is absent in our library due to its length that does not meet catRAPID omics requirements 18 , the omiXcore extension of the algorithm specifically developed for large molecules 41 attributes the same binding propensity to both ADARB1 and ADAR, thus indicating that the interactions might occur (Materials and Methods).
Moreover, experimental works indicate that the family of ADAR deaminases is active in bronchoalveolar lavage fluids derived from SARS-CoV-2 patients 42 and is upregulated in A549 cells infected with SARS-CoV-2 (log fold change of 0.58; p-value of 10 -8 and q-value of 10 -5 ) 40 .
We also identified proteins related to the establishment of integrated proviral latency, including Xray repair cross-complementing protein 5 XRCC5 (P13010) and X-ray repair cross-complementing protein 6 XRCC6 (P12956; Fig. 4B and 4E). In accordance with our calculations, comparison of A549 cells responses to SARS-CoV-2 and respiratory syncytial virus, indicates upregulation of XRRC6 in SARS-CoV-2 (log fold-change of 0.92; p-value of 0.006 and q-value of 0.23) 40 .
Nucleolin NCL (P19338), a protein known to be involved in coronavirus processing, was also predicted to bind tightly to the 5' (Supp. Table 1) 43 .
Importantly, we found proteins related to defence response to viruses, such as ATP-dependent RNA helicase DDX1 (Q92499), that are involved in negative regulation of viral genome replication.
Some DNA-binding proteins such as Cyclin-T1 CCNT1 (O60563), Zinc finger protein 175 ZNF175 (Q9Y473) and Prospero homeobox protein 1 PROX1 (Q92786) were included because they could have potential RNA-binding ability (Fig. 4E) 44 . As for fragment 2, we found two canonical RBPs: More than simple scaffold elements, Gag proteins are versatile elements that bind to viral and host proteins as they traffic to the cell membrane (Supp. Table 1 ) 45 .
Analysis of functional annotations carried out with GeneMania 46 revealed that proteins interacting with the 5' of SARS-CoV-2 RNA are associated with regulatory pathways involving NOTCH2, MYC and MAX that have been previously connected to viral infection processes ( Fig. 4E) 47, 48 .
Interestingly, some proteins, including DDX1, CCNT1 and ZNF175 for fragment 1 and TRIM32 for fragment 2, have been shown to be necessary for HIV functions and replication inside the cell. 
Recently, Gordon et al. reported a list of human proteins binding to Open Reading Frames (ORFs) translated from SARS-CoV-2 58 . Identified through affinity purification followed by mass spectrometry quantification, 332 proteins from HEK-293T cells interact with viral ORF peptides.
By selecting 274 proteins binding at the 5' with Z score ³1.5 (Supp . Table 1) , of which 140 are exclusively interacting with fragment 1 (Fig. 4B) , we found that 8 are also reported in the list by Gordon et al. 58 , which indicates significant enrichment (representation factor of 2.5; p-value of 0.02; hypergeometric test with human proteome in background). The fact that our list of protein-
RNA binding partners contains elements identified also in the protein-protein network analysis is not surprising, as ribonucleoprotein complexes evolve together 14 and their components sustain each other through different types of interactions 16 .
We note that out of 332 interactions, 60 are RBPs (as reported in Uniprot 39 ), which represents a considerable fraction (i.e., 20%), considering that there are around 1500 RBPs in the human proteome (i.e., 6%) and fully justified by the fact that they involve association with viral RNA.
Comparing the RBPs present in Gordon et al. 58 Table 2) . Interestingly, SRP72, LARP7 and LARP4B proteins assemble in stress granules [61] [62] [63] that are the targeted by RNA viruses 53 . We speculate that sequestration of these elements is orchestrated by a viral program aiming to recruit host genes 49 protein kinase A Radixin RDX (P35241; in addition to those mentioned above; Supp. Table 2 ).
In the list of 274 proteins binding to the 5' (fragment 1) with Z score ≥1.5, we found 10 hits associated with HIV (Supp. Table 3 (Q06787) and RNA polymerase-associated protein RTF1 homologue (Q92541; Supp. Table 3 ). By contrast, no significant enrichments were found for other viruses such as for instance Ebola. 65 , among other targets. In addition, HVB-related targets are Nuclear receptor subfamily 5 group A member 2 NR5A2 (CHEMBL3544), Interferon-induced, double-stranded RNA-activated protein kinase EIF2AK2 (CHEMBL5785) and SRSF protein kinase 1 SRPK1 (CHEMBL4375). We hope that this list can be the starting point for further pharmaceutical studies.
A number of proteins identified in our catRAPID calculations have been previously reported to coalesce in large ribonucleoprotein assemblies such as stress granules. Among these proteins, we found double-stranded RNA-activated protein kinase EIF2AK2 (P19525), Nucleolin NCL (P19338), ATP-dependent RNA helicase DDX1 (Q92499), Cyclin-T1 CCNT1 (O60563), signal recognition particle subunit SRP72 (O76094), LARP7 (Q4G0J3) and La-related protein 4B Heterogeneous nuclear ribonucleoprotein Q HNRNPQ (O60506) 63 . To further investigate the propensity of these proteins to phase separate in stress granules, we used the catGRANULE algorithm (Materials and Methods) 66 . We found that the 274 proteins binding to the 5' (fragment 1) with Z score ≥1.5 are highly prone to accumulate in stress-granules (274 proteins with the lowest Z score are used in the comparison; p-value<0.0001; Kolmogorov-Smirnoff; Fig. 4G ; Supp. Table   4 ).
This finding is particularly relevant because RNA viruses are known to antagonize stress granules formation 53 . Indeed, the role of stress granules and processing bodies in translation suppression and RNA decay have impact on virus replication 67 .
spreading. Using advanced computational approaches, we investigated the structural content of SARS-CoV-2 RNA and predicted human proteins that bind to it.
We employed CROSS 13, 68 to compare the structural properties of 2800 coronaviruses and identified elements conserved in SARS-CoV-2 strains. The regions containing the highest amount of structure are the 5' as well as glycoproteins spike S and membrane M.
We found that the spike S protein domain encompassing amino acids 330-500 is highly conserved across SARS-CoV-2 strains. This result suggests that spike S must have evolved to specifically interact with its host partner ACE2 30 and mutations increasing the binding affinity are highly infrequent. As the nucleic acids encoding for this region are enriched in double-stranded content, we speculate that the structure might attract host regulatory elements, thus further constraining the variability. The fact that the Spike S region is highly conserved among all the analysed SARS-CoV-2 strains suggests that a specific drug can be designed against it to prevent interactions within the host.
By contrast, the highly variable region at amino acids 243-302 in spike S protein corresponds to the binding site of sialic acids in MERS-CoV (see manuscript by E. Milanetti et al.) 7, 10, 34 and could play a role in infection 33 . The fact that the binding region changes in the different strains might indicate a variety of binding affinities for sialic acids, which could provide clues on the specific responses in the human population. Interestingly, the sialic acid binding site is absent in SARS-CoV but present in MERS-CoV, which represents an important difference between the diseases.
Both our sequence and structural analyses of spike S protein indicate high conservation among coronaviruses and suggest that human engineering of SARS-CoV-2 is highly unlikely.
Using catRAPID 18, 19 we computed >100000 protein interactions with SARS-CoV-2 and found previously reported interactions such as Polypyrimidine tract-binding protein 1 PTBP1, Heterogeneous nuclear ribonucleoprotein Q HNRNPQ and Nucleolin NCL 43 . In addition, we discovered that the highly structured region at the 5' has the largest number of protein partners including ATP-dependent RNA helicase DDX1 that was previously reported to be essential for HIV-1 and coronavirus IBV replication 54,55 , and the double-stranded RNA-specific editases ADAR and ADARB1 that catalyse the hydrolytic deamination of adenosine to inosine 49 . Other predicted interactions are XRCC5 and XRCC6 members of the HDP-RNP complex associating with ATP-dependent RNA helicase DHX9 69 as well as and 2-5A-dependent ribonuclease RNASEL and 2-5-oligoadenylate synthase 2 OAS2 that control viral RNA degradation 70, 71 . Interestingly, DDX1, XRCC6 and OAS2 are upregulated in human alveolar adenocarcinoma cells infected with SARS-CoV-2 40 . In agreement with our predictions, recent experimental work indicates that the family of ADAR deaminases is active in bronchoalveolar lavage fluids derived from SARS-CoV-2 patients 42 .
A significant overlap exists with the list of protein interactions reported by Gordon et al. 58 , and among the candidate partners we identified AKAP8L, involved as a DEAD/H-box RNA helicase binding protein involved in HIV infection 59 . In general, proteins associated with retroviral replication are expected to play different roles in SARS-CoV-2. As SARS-CoV-2 massively represses host gene expression 49 , we hypothesize that the virus hijacks host pathways by recruiting transcriptional and post-transcriptional elements interacting with polymerase II genes and splicing factors such as for instance A-kinase anchor protein 8-like AKAP8L and La-related protein 7 LARP7 that is upregulated in human alveolar adenocarcinoma cells infected with SARS-CoV-2 40 .
The link to proteins previously studied in the context of HIV and other viruses, if further confirmed fro, is particularly relevant for the repurposing of existing drugs 65 .
The idea that SARS-CoV-2 sequesters different elements of the transcriptional machinery is particularly intriguing and is supported by the fact that a large number of proteins identified in our screening are found in stress granules 63 . Indeed, stress granules protect the host innate immunity and are hijacked by viruses to favour their own replication 67 . Moreover, as coronaviruses transcription uses discontinuous RNA synthesis that involves high-frequency recombination 43 , it is possible that pieces of the viruses resulting from a mechanism called defective interfering RNAs 72 could act as scaffold to attract host proteins 14, 15 .
We predicted the secondary structure of transcripts using CROSS (Computational Recognition of Secondary Structure 13, 68 . CROSS was developed to perform high-throughput RNA profiling. The algorithm predicts the structural profile (single-and double-stranded state) at single-nucleotide resolution using sequence information only and without sequence length restrictions (scores > 0 indicate double stranded regions). We used the Vienna method 27 to further investigate the RNA secondary structure of minima and maxima identified with CROSS 13 .
We used CROSSalign 13,68 an algorithm based on Dynamic Time Warping (DTW), to check and evaluate the structural conservation between different viral genomes 13 . CROSSalign was previously employed to study the structural conservation of ~5000 HIV genomes. SARS-CoV-2 fragments (1000 nt, not overlapping) were searched inside other complete genomes using the OBE (open begin and end) module, in order to search a small profile inside a larger one. The lower the structural distance, the higher the structural similarities (with a minimum of 0 for almost identical secondary structure profiles). The significance is assessed as in the original publication 13 .
The FASTA sequences of the complete genomes of SARS-CoV-2 were downloaded from Virus Pathogen Resource (VIPR; www.viprbrc.org), for a total of 62 strains. Regarding the overall coronaviruses, the sequences were downloaded from NCBI selecting only complete genomes, for a total of 2862 genomes. The reference Wuhan sequence with available annotation (EPI_ISL_402119) was downloaded from Global Initiative on Sharing All Influenza Data. (GISAID https://www.gisaid.org/).
Interactions between each fragment of target sequence and the human proteome were predicted using catRAPID omics 18, 19 , an algorithm that estimates the binding propensity of protein-RNA pairs by combining secondary structure, hydrogen bonding and van der Waals contributions. As 
We used Clustal W 28 for 62 SARS-CoV-2 strains alignments and T-Coffee 32 for spike S proteins alignments. The variability in the spike S region was measured by computing Shannon entropy on translated RNA sequences. The Shannon entropy is computed as follows:
S(a) = -Sum_i P(a,i) log P(a,i)
Where a correspond to the amino acid at the position i and P(a,i) is the frequency of a certain amino-acid a at position i of the sequence. Low entropy indicates poorly variability: if P(a,x) = 1 for one a and 0 for the rest, then S(x) =0. By contrast, if the frequencies of all amino acids are equally distributed, the entropy reaches its maximum possible value.
catGRANULE 66 was employed to identify proteins assembling into biological condensates. Scores > 0 indicate that a protein is prone to phase separate. Structural disorder, nucleic acid binding propensity and amino acid patterns such as arginine-glycine and phenylalanine-glycine are key features combined in this computational approach 66 . Fig. 1 