key: cord-334891-4jgtxg07
authors: Choudhury, Abhigyan; Chandra Das, Nabarun; Patra, Ritwik; Bhattacharya, Manojit; Mukherjee, Suprabhat
title: In silico analyses on the comparative sensing of SARS-CoV-2 mRNA by intracellular TLRs of human
date: 2020-11-11
journal: bioRxiv
DOI: 10.1101/2020.11.11.377713
sha: 
doc_id: 334891
cord_uid: 4jgtxg07

The worldwide outbreak of COVID-19 pandemic caused by SARS-CoV-2 leads to loss of mankind and global economic stability. The continuous spreading of the disease and its pathogenesis takes millions of lives of peoples and the unavailability of appropriate therapeutic strategy makes it much more severe. Toll-like receptors (TLRs) are the crucial mediators and regulators of host immunity. The role of several TLRs in immunomodulation of host by SARS-CoV-2 is recently demonstrated. However, the functionality of human intracellular TLRs including TLR3,7,8 and 9 is still being untested for sensing of viral RNA. This study is hoped to rationalize the comparative binding and sensing of SARS-CoV-2 mRNA towards the intracellular TLRs, considering the solvent-based force-fields operational in the cytosolic aqueous microenvironment that predominantly drive these reactions. Our in-silico study on the binding of all mRNAs with the intracellular TLRs shown that the mRNA of NSP10, S2, and E proteins of SARS-CoV-2 are potent enough to bind with TLR3, TLR9, and TLR7 and trigger downstream cascade reactions, and may be used as an option for validation of therapeutic option and immunomodulation against COVID-19.

The worldwide outbreak of Coronavirus disease or COVID-19 pandemic caused by the 2 Severe Acute Respiratory Syndrome Coronavirus-2 (SARS-CoV-2), leads to the infection of 3 about 0.33% of the total world's population causing the death of 1.24 million people by the 4 first week of November 2020. 1 The continuous expansion of contagion and pathogenesis, it is 5 taking the shape of another dark age in the history of mankind, not only to health crises but 6 also in bankrupting the global socio-economic status. The SARS-CoV-2 belongs to the β-7 Coronavirus genus of a 2B group of the Coronaviridae family and is considered the third most 8 virulent type, leading to the highest fatality rate in humans followed by the SARS-CoV and 9 MERS-CoV. 2 It transmits from person to person mainly via close physical contact and by 10 respiratory aerosols that produces during coughing, sneezing, and even talking within 11 proximity, although some recent studies indicate transmission through fecal matters and 12 fomite-borne contaminations. 3-5 13 The proteome of the virus consists of structural proteins, existing in three formsthe 14 Spike 'S' protein, the Envelope 'E' protein and the Membrane 'M' protein. Along with these, 15 16 forms of non-structural proteins (NSPs) combining together to generate different catalytic 16 models. 6 The genome, however, is much simpler and comprises of a 29,903 base long, positive-17 sense single-stranded RNA molecule. This (+)ssRNA genome makes it further feasible to be 18 detectible by the intracellular Toll-like receptors (TLRs) which have a high affinity towards 19 nucleic acid-related pathogen-associated molecular patterns (PAMPs). 7 Several empirical 20 studies have presented various propositions relating the intracellular TLR 3, 7, 8, and 9 to the 21 cytokine storm produced by the virus that majorly owes it the lethality. 8 Cytokine storm is 22 apparently the incessant extreme activation of cytokine production leading to prolonged and 23 consistent inflammatory response, which becomes almost continuous due to positive feedback 24 4 loops in the TLR signaling pathways. The studies are quite successful in relating the role of 1 TLRs in recognizing oligonucleotide PAMPs and triggering the cytokine storm. 9 2

The binding of Spike protein with the human ACE2 receptor triggers the pathogenesis 3 of the SARS-CoV-2, leading to the activation of TLRs to activate the proliferation and 4 production of pro-inflammatory cytokines causing cytokine storm, those results in 5 inflammations. 10, 11 From previous studies, it has been found that the spike protein shows 6 binding efficiency with the extracellular domains of TLRs including TLR1, TLR4 and, TLR6, 7 with the strongest affinity with TLR4. 12 Furthermore, the development of several in-silico 8 multi-epitope-based peptide vaccine candidates against the SARS-CoV-2 has shown to be 9 effectively binding with TLR3, TLR4 and, TLR5 to regulate the TLR signaling pathways 10 activation and proliferation. 13 It has been found that targeting human TLRs, in an order to either 11 block the binding of SARS-CoV-2 or by inhibiting the TLR activation and proliferation that 12 induce the production of pro-inflammatory cytokines using certain TLR agonists, might be 13 used as an effective therapeutic strategy against coronavirus disease. 11 Additionally, the 14 continuous replication and generation of new virus within the host system indicates that there 15 must be binding efficiency and potency with the intracellular TLRs including TLR3, 7, 8, and 16 9 in order to induce the severity and pathogenesis of the virus and further introduction to the 17 new host. TLR3 is well known for its sensing capabilities of viral PAMPs and exists as a 18 monomer that is attached to the membrane of endosomes, thereby it detects and binds to certain 19 motifs in the invading viral RNA. 14 TLR7 expressed majorly in the cerebral cortex, lung, 20 bronchus, breast, kidney, rectum, and smooth muscle tissues, and functions by adhering to the 21 endosome and thereby binding to the viral RNAs with high guanosine and uridine content 15 to 22 initiates a MyD88 signal transduction resulting in activation of NF-κB, mitogen-activated 23 protein kinase (MAPK) cascades, as well as IRF-7 and IRF-5 activation via IL-1 receptor-24 associated kinases (IRAK)-1/2/4 and TNF receptor-associated factor-3/6. 16 The signaling 25 5 finally induces the production of pro-inflammatory cytokines including IL-1β, IL-6, IL-12, 1 TNF-α, and IFN-α. 17 TLR8 is expressed predominantly in the lung and the peripheral blood 2 leucocytes and plays a major role in recognizing GU-rich viral ssRNAs including those of 3 SARS-CoV-2. 18 TLR9 is predominantly expressed in the spleen, lymph node, bone marrow, 4 and peripheral blood leukocytes and it recognizes CpG motifs in viral DNAs. However, several 5 studies suggest its role in sensing ssRNA fragments generated by the SARS-CoV-2 genome. 19 6 However, a lack of precise knowledge regarding the nature of oligonucleotides and 7 their ligating affinity towards the TLRs still pertains to exist. This in context is abstaining the 8 medical research community from developing certain therapeutic interventions that would have 9 been vitally important in this hour of severity. Our study is hoped to rationalize the picture and 10 provide clues regarding the interaction range of the oligonucleotides towards the intracellular 11 TLRs, considering the solvent-based force-fields operational in the cytosolic aqueous 12 microenvironment that predominantly drive these reactions. 13 15 Current literatures suggest the predominant expression of ten proteins of SARS-CoV-2. 20 Four 16 of them are structural proteinsthe spike protein subunits S1 and S2, the envelope protein (E), were obtained from the RCSB PDB database. However, the structures of TLR7 and TLR9 were 2 obtained by performing homology modelling in the SWISS-MODEL server from ExPASy 3 (https://swissmodel.expasy.org/) by using protein sequences obtained from GenBank with 4 Accession No. AAZ99026 and AAZ95518.1 respectively. All the crystal surface structures of 5 TLRs and the whole genome of SARS-CoV-2 were then prepared to visualize ( Figure 1 ).

To perform further studies, PyMOL is used for removing water molecules, aboriginal 7 heteroatoms as well as any xenobiotic ligands whenever present, and adding polar hydrogens 8 and Kollman charges to the structures. 10 The retrieved RNA sequences were subjected to the imRNA tool developed by IIIT-Delhi 11 (https://webs.iiitd.edu.in/raghava/imrna/). This tool is based on Motif-EmeRging and with 12 Classes-Identification (MERCI) program and scans through the RNA sequences for motifs 13 that are potent to have immunomodulatory properties. 21 

The HDOCK server (http://hdock.phys.hust.edu.cn/) presents a novel algorithm which is a 1 hybrid of template-dependent along with template-independent ab initio free docking. 2 Moreover, it is one of the advance programme which support protein docking against 3 DNA/RNA molecules. 24 Oligonucleotides have great sizes and demand heavier computing 4 resources for the rendering of molecular models, which is seldomly feasible for any 5 supercomputing server to provide. Thus, HDOCK was a program of choice. The retrieved PDB 6 structures along with the RNA sequences were used for the purpose, however, the species- 16 The so retrieved docked structures were then fed into the iMOD webserver from ChaconLab 17 (http://imods.chaconlab.org). This program has a user-friendly GUI and is a well-recognized 18 tool for performing normal mode analysis (NMA) and simulating various modal trajectories of 19 protein dynamics. NMA in dihedral coordinates naturally mimics the combined functional 20 motions of protein molecules modelled as a set of atoms connected by harmonic springs. 26 As 21 output the server delivers affine modelled arrow, vector field and a modal animation to signifies 22 the motions. Moreover, study also provide more detail profile about mobility using B-factor 23 and deformability plots whereas eigenvalue helps to measure the relative modal stiffness of the 3 In-silico molecular docking between receptor and ligand is not sufficient to conclude the nature In this work, Chemistry at Harvard Macromolecular Mechanics (CHARMM) force field was 10 accessed initially to generate necessary input files for MD simulation. 29-31 Solution builder 11 approach under CHARMM-GUI web tool was first add water box (TIP3 216) and then 12 neutralizing atoms to solvate the system. In order to remove bad contacts and generate more 13 specific outcomes periodic boundary conditions (PBC) were analysed and minimization steps 14 were performed sequentially. Server provided outcomes were then applied on GROMACS v. 15 2020.1 to achieve the equilibrium (NVT-constant volume and temperature) of the system. 16 Finally, module gmx_mdrun was accessed to run the MD simulation of the system and xmgrace 21 Protein with conformational changes from its native form suggests it is in bounded state with 22 any ligand. 32 Thus, to analyse the changes in TLRs, firstly we extract the PDBs from 23 trajectories and then visualized using PyMOL. In order to assesse the quality of modelled structure of TLR7 and TLR9, we access 3 Structure Assesement tool under SWISS-MODEL server and found out both the structures are 4 significant as more than 92% residues from both structures reside in Ramachandran plot 5 favoured region signifying a stable stereochemical structure ( Figure S1 ). 15 Our study focuses principally on the interactions of the intracellular TLRs with the mRNA 16 fragments, and molecular docking studies play a pivotal role in this study by precisely Table S1 . However, only the four topmost scoring docked complexes were selected for 1 further experimentation and analysis. 

3

In accord to our docking studies, TLR3 binds with a much proficient binding pose with 4 the RNA fragment encoding NSP10, with a high docking score of -404.77 which is higher than 5 any other docking operation involving TLR3 like that of the RNA fragment of papain-like 6 protease which binds with a score of -346.51 or that of the main protease which has a score of 7 -333.76 (Table S1 ). While in search of insights of molecular docking PLIP study found the 8 involvement of 251Arg, 252Asn, 277Thr residues in hydrogen bonding, 326Tyr, 359His in π- 9 Cation interactions and Salt bridges from 32 and 319His ( Figure 2A , Table 1 ) 11 In this study, the docking operations determine a good binding pocket for E-protein Figure 2B and Table 1 . 18 Herein, the most probable binding pose of NSP8 RNA fragment and TLR8 is built with 19 a definite score of -416.84, which is much higher than those poses comprising of the RNAs of Table 1 . 3 Among the docked complexes with TLR9 possessing the S2 subunit mRNA fragment 4 stood at the highest position with a score of -440.33. The molecular insights of this ligation are 5 well supported by the PLIP analysis as finding involves 224Tyr, 247Arg, 248Val, 292Lys 6 residues in hydrogen bonding, 311Arg in π-Cation interactions and Salt bridges from 207Lys, 7 287Glu, and 292Lys ( Figure 3B , Table 1 ). 8 Figure S3I , S3J, S4I and S4J, where red colour represents individual and green to 21 cumulative variance. In the covariance matrix, motions are categorized in different modes and 22 visualized in Figure S3K and S3M with red, blue and white colour that signify variations 23 among correlated, anti-correlated and uncorrelated motions respectively. In the end, the last 24 outcome of NMA study, elastic network map graph is represented by Figure S4L and S4N   1 where grey dots indicate the stiffness of the motions and springs direct the atomic connections 2 of the complex. 4 Four most active complexes, those were primarily strained out according to their significant Herein, RMSD plot in Figure 4A reflects the structural stability of the backbone of different 10 receptor TLRs of the four complexes and suggests complex TLR8-NSP8 is less stable among 11 others. Where Figure 4B and 4C define the compactness and solvent accessible area of that 12 complexes as Rg plot and SASA plot, respectively. At last the graphical plot of hydrogen bond 13 clearly describes the insights of the complexes and supports the postulation reflected on 14 RMSD, Rg and SASA plots ( Figure 4D ). 15 16 Apart from RMSD, Rg, SASA and hydrogen bond studies, we also analyse conformational Figure S2E ).

14 Further, the molecular docking studies showed a similar trend by reflecting highest docking 15 score to the complexes selected through RPISeq tool and in turn support the prediction of best 16 docked structures (Table S1 ). After getting confirmation on best predicted structure through 17 HDOCK we forwarded to next step of simulation to analyze the conformational stability and 18 flexibility of that selected four structures viz. TLR8-NSP8 mRNA complex ( Figure 3A ),

TLR3-NSP10 mRNA complex as shown in Figure 2A , TLR7-E mRNA complex in Figure 2B 20 and TLR9-S2 mRNA complex ( Figure 3B) . Initially, the NMA study reveals all the compounds 21 have better stability as it reflects quite similar and significant eigenvalue scores, those are 22 3.154899x10 -5 for TLR8-NSP8 mRNA complex as shown in Figure S4G , 1.27056x10 -5 TLR3-23 NSP10 mRNA complex ( Figure S3G) , 7.042567x10 -5 TLR7-E mRNA complex as in Figure   24 S3H , and 9.131746 x10 -5 TLR9-S2 mRNA complex ( Figure S4H ). Later, MD simulation 25 studies figure out postulation of NMA studies is not fully accepted as TLR8-NSP8 mRNA 1 complex found as unstable throughout the simulation process in the RMSD plot Figure 4A . 2 The PLIP analysis Table 1 and hydrogen bond plot Figure 4D We do acknowledge the efforts of all the doctors, health workers, social workers, scientists, 22 and researchers currently working endlessly against coronavirus disease worldwide. We have 23 kept selected studies as reference due to the limitation in space, but we appreciate all the uncited 24 articles dedicated to the research on coronavirus. RP thanks the Department of Higher π-Cation interactions, while dotted yellow lines denote Salt bridges. It is to be noted that for 12 the sake of perceptual clarity some residues have been intentionally omitted from the diagram, 13 however the same have been furnished in table 1. 

World Health Organization. COVID-19 Pandemic update. Situation report 10 website

The continuing 2019-nCoV epidemic threat of novel 12 coronaviruses to global health-The latest 2019 novel coronavirus outbreak in Wuhan

Coughs and Sneezes: Their Role in Transmission of Respiratory Viral 15

Including SARS-CoV-2

Potential fecal transmission of SARS-CoV-2: Current evidence and 18 implications for public health

Transmission of 21 SARS and MERS coronaviruses and influenza virus in healthcare settings: the possible role of 22

Toll-like 11 receptors and COVID-19: a two-faced story with an exciting ending

COVID-19: towards understanding of pathogenesis

Targeting human TLRs to combat COVID-19: A solution?

In silico studies on the comparative characterization of the 18 interactions of SARS-CoV-2 spike glycoprotein with ACE-2 receptor homologs and human 19

Development of epitope-based peptide vaccine 21 against novel coronavirus 2019 (SARS-COV-2): Immunoinformatics approach

The Toll for Trafficking: Toll-Like Receptor 7 Delivery to 1 the Endosome

Is a Dual Receptor for Guanosine and Single-Stranded RNA

Small anti-viral compounds activate immune cells 7 via the TLR7 MyD88-dependent signaling pathway

TLR2 and TLR4 mediated host immune responses 10 in major infectious diseases: a review. The Brazilian Journal of Infectious Diseases

Recognition of GU-rich polyadenylation regulatory 13 elements by human CstF-64 protein

Coronavirus infections and immune responses

The Proteins of Severe Acute Respiratory Syndrome Coronavirus-2 18 (SARS CoV-2 or n-COV19), the Cause of COVID-19

Identifying discriminative classification-based motifs 21 in biological sequences

Learning Framework for 1 Robust and Accurate Prediction of ncRNA-Protein Interactions Using Evolutionary 2 Information

Predicting RNA-Protein Interactions Using Only 5

Sequence Information

HDOCK: a web server for protein-protein 8 and protein-DNA/RNA docking based on a hybrid strategy

PLIP: fully automated 11 protein-ligand interaction profiler

iMODS: internal coordinates 14 normal mode analysis server

Designing of a novel multi-epitope peptide based 17 vaccine against Brugia malayi: An in silico approach

Molecular Modeling of Nucleic Acid Structure: 20 Setup and Analysis

Protein conformational switches: from nature to design

Computational Prediction of the 13

Immunomodulatory Potential of RNA Sequences

In silico study the inhibition of 16 angiotensin converting enzyme 2 receptor of COVID-19 by Ammoides verticillata components 17 harvested from Western Algeria

Molecular Dynamics Simulation for All