Designing potential siRNA molecule for the nucleocapsid(N) gene silencing of different SARS-CoV-2 strains of Bangladesh: Computational approach.

SARS-CoV-2 is a single-stranded RNA (+) virus first identified in China and then became an ongoing global outbreak. In most cases, it is fatal in humans due to respiratory malfunction. Extensive researches are going to find an effective therapeutic technique for the treatment of SARS-CoV-2 infected individuals. In this study, we attempted to design a siRNA molecule to silence the most suitable nucleocapsid(N) gene of SARS-CoV-2, which play a major role during viral pathogenesis, replication, encapsidation and RNA packaging. At first, 270 complete N gene sequences of different strains in Bangladesh of these viruses were retrieved from the NCBI database. Different computational methods were used to design siRNA molecules. A siRNA molecule was built against these strains using the SiDirect 2.0 server. Using Mfold and the OligoCalc server, the siRNA molecule was tested for its secondary structure and GC material. The Clustal Omega tool was employed to evaluate any off-target harmony of the planned siRNA molecule. Herein, we proposed a duplex siRNA molecule that does not fit any off-target sequences for the gene silencing of SARS-CoV-2. To treat SARS-CoV-2 infections, currently, any effective therapy is not available. Our engineered siRNA molecule could give an alternative therapeutic approach against various sequenced SARS-CoV-2 strains in Bangladesh.

Introduction

A deadly bat-borne Coronaviridae family virus SARS-CoV-2 was identified first in mid-December of 2019 in Wuhan, China (Chan et al., 2020). Since then, it has become an ongoing outbreak worldwide. To this date, it has affected 213 countries and territories across the globe (Centers-of-Disease-Control-and-Prevention, 2020). In this positive sense, the single-stranded RNA virus has already infected about thirteen million people worldwide till July 2020 (Anon., 2021). Patients suffering from SARS-CoV-2 infection are usually found with complications in the respiratory tract, diarrhea, high fever, thrombocytopenia, lymphopenia, and increased C-reactive protein and lactate dehydrogenase levels within the first 3–6 days of viral exposure (Chan et al., 2020; Zhang et al., 2020). With a genome size of ∼30000 bases, one of the largest nonsegmented genomes among all RNA viruses, this virus has already compiled plentiful changes in its genome and procures more (Forni et al., 2017; Woo et al., 2009). Due to the sequencing technologies advancement, various databases preserves the available whole genome sequence of these deadly virus strains now. Researchers are looking for effective molecular therapy against the SARS-CoV-2 using these resources. Bangladeshi scientists also sequenced different SARS-CoV-2 strains, where 273 sequences were available in the NCBI public database until December 2020 (Saha et al., 2020; Moniruzzaman et al., 2020). A comparative analysis computed the origin of the various strains of this virus in Bangladesh among publicly available SARS-CoV-2 genome sequences from 27 countries. Their prediction about Bangladeshi strains based on phylogenetic analysis reported that the pathogen appeared from multiple countries (Shishir et al., 2021). This study was on 64 SARS-CoV-2 isolated genomic sequences from Bangladesh with the reference (NC_045512.2 first Chinese SARS-CoV-2 sequence) and 433 sequences of the whole world available from several countries. The sequences of different Bangladeshi strains were showed three clusters, among them 43 of the 64 sequences shared a comparable node with Germany, which carried a common ancestor with the United Kingdom. Two other clusters of Bangladesh had 4 and 5 sequences, respectively, and in both cases, they shared the same node with the sequence of SARS-CoV-2 reported from India and Saudi Arabia.

On the other hand, 12 sequences that did not belong to any clusters were found to be similar with sequences from Europe, including United Kingdom, Germany, France, Italy, and Russia, while one of these sequences was close to the sequence reported from the USA (Shishir et al., 2021; Hasan et al., 2020). The underlying link of Bangladeshi SARS-CoV-2 isolates with part of a haplotype observed high in Europe (Shishir et al., 2021; Mahmud et al., 2020). Though the virus is a mutated particle over time, the SARS-CoV-2 sequence shares 79.6 % similarity to SARS-CoV (Zhou et al., 2020; Marra et al., 2003), with the fewer mutations and higher (90 %) homology showing N gene, which is more stable and conserved over time, but the S gene has 76 % similarity (Marra et al., 2003; Drosten et al., 2003; Grifoni et al., 2020; Holmes and Enjuanes, 2003; Rota et al., 2003; Zhu et al., 2005). Nucleocapsid gene silencing Small Interfering RNA (siRNA) molecules potentiality has already been reported (Ge et al., 2003). Double-stranded 20–25 base pairs containing RNA generally work through the RNA interference (RNAi) pathway to silence a specific gene expression, where it causes mRNA breakage into the following transcription (Agrawal et al., 2003). These mRNA degrading small molecules play an important role in the post-transcriptional gene silencing (PTGS) (Hamilton and Baulcombe, 1999).

In-vitro studies have verified that siRNA can significantly repress gene expression in mammalian cells (Lee et al., 2002; Kapadia et al., 2003; Lee et al., 2003). In addition, effective in vivo silencing of the endogenous gene and transgene expression is already revealed (McCaffrey et al., 2002; Giladi et al., 2003; Sorensen et al., 2003; Tompkins et al., 2004a). The capability of small interfering RNA is already proved in a in vitro study through the inhibition of SARS-CoV N gene expression in cultured cells and mouse muscles (Zhao et al., 2005). More than 80 % blocking activity was shown utilizing chemically synthesized siRNA duplexes targeting genomic RNA of SARS-CoV (Zheng et al., 2004). Like other viruses, SARS-CoV-2 nucleocapsid(N) protein is a multifunctional protein and plays a crucial role in encapsidation, viral RNA transcription, and replication (Chang et al., 2016).

For this reason, the study of the nucleocapsid protein or N gene has become much popular as a diagnostic and therapeutic target (Liu et al., 2020; Eshaghi et al., 2005). However, some studies show a mutation in the N gene, which is not the focus of this study. Moreover, our main objective is to identify a common therapeutic target for different SARS-CoV-2 strains of Bangladesh. Targeting this sequence may cause SARS-CoV-2 viral inhibition like SARS-CoV by impeding the N gene's translation and thus preventing the RNA transcription and replication (Kapadia et al., 2003; Wu et al., 2005).

Materials and methods

The complete workflow is summarized into the methodology (Fig. 1 )

Fig. 1

The potential siRNA molecule prediction methodology.

Viral strain selection

The selection of SARS-CoV-2 strains and their associated information, including their genus, family, host, transmission pattern, disease pathogenicity, genome, proteome, and the available therapeutic agents against them, was identified ViralZone (http://viralzone.expasy.org/) of the Ex-PASy Bioinformatics Resource Portal.

Sequence retrieval and evolution analysis

Two hundred and seventy-three sequences of different SARS-CoV-2 (BDG) strains were collected from the viral gene bank database available at NCBI (http://www.ncbi.nlm.nih.gov/). Among them, 270 complete cds of nucleocapsid(N) protein gene was selected, and the other 3 incomplete sequences were disqualified for the next experiment. Multiple sequence alignment of the retrieved sequences was done by Clustal Omega (Sievers et al., 2011) (https://www.ebi.ac.uk/Tools/msa/clustalo/), and phylogeny was analyzed by Phylogeny analysis (http://www.phylogeny.fr/simple_phylogeny.cgi) tool.

Target prediction and allowable siRNA molecule devising

A web server siDirect 2.0 (http://siDirect2.RNAi.jp/) (Naito et al., 2009) was used for efficient and target-specific siRNA design for mammalian RNAi. It utilized Ui-Tei, Amarzguioui, and Reynolds rules combined (Ui‐Tei et al., 2004; Reynolds et al., 2004; Amarzguioui and Prydz, 2004), and as a parameter, melting temperature (Tm) below 21.5 °C for siRNA duplex was also used. Besides, these other parameters were taken on the concept of algorithms given in Table 1 .

Table 1

Algorithms or rules for the rational design of siRNA molecules.

No.	Rules Name	Rules
1.	Ui-Tei rules	1 A/U at the 5′ terminus of the sense strand 2 G/C at the 5′ terminus of the antisense strand 3 At least 4 A/U residues in the 5′ terminal 7bp of the sense strand 4 No GC stretch longer than 9nt
2.	Amarzguioui rules	1 Duplex End A/U differential > 0. 2 Strong binding of 5′ sense strand 3 No U at position 1. 4 Presence of A at position 6. 5 Weak binding of 3′ sense strand.
3.	Reynolds rules	1 GC content 30 %–52 % (1 point) 2 Occurrence of 3 or more 3 A/U base pair at position 15−19 of the sense strand (Each A/U base pair in this region earns 1 point) 4 Low internal stability at the target site (melting temperature Tm>-20°c) (1 point) 5 Presence of A at position 19 of the sense strand (1 point) 6 Presence of A at position 3 of the sense strand (1 point) 7 Presence of U at position 10 of the sense strand (1 point) 8 Absence of G at position 13 of the sense strand (1 point) 9 Threshold for efficient siRNAs score> = 6

Off target harmony investigation

Blast tool (Johnson et al., 2008) (http://www.ncbi.nlm.nih.gov/blast) was used to identify off-target similarity with any sequence on whole Gene bank datasets other than the target sequence by applying expected threshold value 10 as a parameter.

GC content count and secondary structure prophecy

Oligonucleotide Properties Calculator for GC content calculation of predicted siRNA, OligoCalc (Kibbe, 2007) (http://basic.northwestern.edu/biotools/OligoCalc.html) tool was used while Mfold server (Zuker, 2003) (http://www.mfold.-rna.albany.edu/) was used for secondary structure prediction aimed to compute the free energy of folding.

Calculation of RNA-RNA interaction through thermodynamics

To study the thermodynamics of interaction between predicted siRNA and target gene RNAcofold program (Gruber et al., 2008) (http://rna.tbi.univie.ac.at/cgi-bin/RNAcofold.cgi) was used. The hybridization energy and base-pairing form of two RNA sequences was calculated by it. It functions as an extension of McCaskill's partition function algorithm to compute base-pairing probabilities, realistic interaction energies, and equilibrium concentrations of duplex structures.

Verification of the considered siRNA

Finally, the siRNAPred (Kumar et al., 2009) (http://imtech.res.in/raghava) server from Imtech was used to validate the predicted siRNA further. Here we used efficacy prediction for 21 mers. The predicted siRNA was screened against the Main 21 dataset using a binary pattern. siRNAPred score greater than 0.9 predicts very high efficacy, a score ranging 0.8−0.9 indicates high efficacy, and a score ranging 0.7−0.8 predicts moderate efficacy.

Results and discussion

At first, N gene cds of different SARS-CoV-2 strains in Bangladesh were picked to design a siRNA molecule. The N gene sequence of the SARS-CoV-2 Bangladeshi strain (MT476385.1) was first sequenced by the child health research foundation (CHRF). The 270 complete sequences were collected from NCBI. The accession numbers and Phylogenetic tree of these complete cds were shown by date released (Supplementary file 1 and Supplementary file 2). It was stated from the phylogenetic analysis tool that all the sequences had a common predecessor and some significant harmony reserved during evolution (Fig. 2 ).

Fig. 2

Phylogenetic analysis (treedyn) of the nucleocapsid(N) gene of different SARS-CoV-2 (BDG) strains.

Extremely fruitful small interfering RNA with maximal target specificity from the retrieved sequences computed by the siDirect web-based online software system. The proposed consensus target for siRNA with location is shown in Table 2 . The Clustal Omega server did multiple sequence alignment (MSA) among the different SARS-CoV-2 cds of the N gene (Supplementary file 3). Web server siDirect predicts siRNA by calculating the Tm of the seed target duplex using the nearest neighbor model and the thermodynamic parameters in Table 3

Table 2

Selected siRNA and their location with the consensus target of the N gene cds.

Total Accession	Target	Location of the target within the gene	siRNA target sequence within the gene
270	Consensus (270/270)	314−336	GTCCAAGATGGTATTTCTACTAC
	Consensus (270/270)	727−749	GGCCAAACTGTCACTAAGAAATC
	Consensus (269/270)	789−811	TGCCACTAAAGCATACAATGTAA
	Consensus (270/270)	892−914	TACAAACATTGGCCGCAAATTGC

Table 3

Top four RNA oligo sequences against the targeting 21 nt of target sequences.

Target sequence 21 nt target + 2 nt overhang	RNA oligo sequences 21 nt guide (5′→3′) 21 nt passenger (5′→3′)	Functional siRNA selection: Ui-Tei Reynolds Amarzguioui
GTCCAAGATGGTATTTCTACTAC	AGUAGAAAUACCAUCUUGGAC	U R A
	CCAAGAUGGUAUUUCUACUAC
GGCCAAACTGTCACTAAGAAATC	UUUCUUAGUGACAGUUUGGCC	U R A
	CCAAACUGUCACUAAGAAAUC
TGCCACTAAAGCATACAATGTAA	ACAUUGUAUGCUUUAGUGGCA	U R A
	CCACUAAAGCAUACAAUGUAA
TACAAACATTGGCCGCAAATTGC	AAUUUGCGGCCAAUGUUUGUA	U R A
	CAAACAUUGGCCGCAAAUUGC

The formula for calculating the Tm is:Where ΔH (kcal/ mol) is the sum of the nearest neighbor enthalpy change, A is the helix initiation constant (-10.8), ΔS is the sum of the nearest neighbor entropy change.37 R is the gas constant (1.987 cal/deg/mol), and CT is the total molecular concentration of the strand (100 μM). [Na+] was fixed at 100 mM.

Our targeted siRNA's GC content was calculated 38 %, 40 %, 36 %, 43 % by the OligoCalc calculator, an oligonucleotide features Counter. Mfold server predicts RNA secondary structure through widely used algorithms based on a minimal free energy state (Zuker, 1989). The RNAcofold server from Vienna RNA web services calculated the hybridization energy and base-pairing pattern of the RNA sequences (Table 4 ). The siRNAPred server assessed the 21mer siRNA through the Support vector machine-based methods with high accuracy. The Main21 dataset of the siRNAPred server consists of 2182 siRNAs (21mer) derived from homogeneous experimental conditions. The binary model was chosen for the justification. These results in the top calculation score of each 21mer siRNA were 0.946, 0.861, 0.986, 0.793. This score lies within the range of very high efficacy, high efficacy, and moderate prediction score for siRNAPred server (≥0.9, very high efficacy; 0.8−0.9, high efficacy; 0.7−0.8, moderate efficacy). Potential treatment can be provided through siRNA therapeutic approaches for the diseases hereditary genetic defects, viral infectious diseases, immune disorders, and cancers caused by a particular gene set (Aagaard and Rossi, 2007). The siRNA to silence specific genes against several viruses such as hepatitis C virus (Gruber et al., 2008), HIV-1 infection (Kumar et al., 2009), and herpes simplex virus 2 infections is already approved as a potential (Palliser et al., 2006).

Table 4

The proposed siRNA molecule with GC%, the free energy of binding with target and validity.

Target	Predicted duplex siRNA candidate at 37 °C	GC%	Free energy of binding (Kcal/mol)	Validity(Binary)
Consensus target	AGUAGAAAUACCAUCUUGGAC	38 %	−31.50 Kcal/mol	0.946
	CCAAGAUGGUAUUUCUACUAC
	UUUCUUAGUGACAGUUUGGCC	40 %	−34.54 Kcal/mol	0.861
	CCAAACUGUCACUAAGAAAUC
	ACAUUGUAUGCUUUAGUGGCA	36 %	−30.74 Kcal/mol	0.986
	CCACUAAAGCAUACAAUGUAA
	AAUUUGCGGCCAAUGUUUGUA	43 %	−31.61 Kcal/mol	0.793
	CAAACAUUGGCCGCAAAUUGC

Designing a fruitful siRNA molecule to target a particular gene has a number of challenges. Finding an operative delivery manner is the main challenge; further challenges are off-target silencing, the creation of the immune reaction, and finally, the siRNA stability (Gavrilov and Saltzman, 2012). Thus, scientists have improved some models to prophesy the efficacy for a deliberate siRNA fragment before in vivo investigation– named Ui-Tei rules, Amarzguioui rules, and Reynolds rule (Taxman et al., 2006). In this study, we also follow all of these rules to designing a siRNA molecule. The proposed siRNA is covered the threshold score 6 of Reynolds rules, which indicates a proficient siRNA. The GC substance was also found to be within the supported range of 30–52 % of the proposed siRNA (Chan et al., 2009). Another important parameter for siRNA efficiency is the prediction of thermodynamics of RNA-RNA interaction which was predicted by the RNAcofold server. Here the two RNA sequences are concatenated, and an ampersand specifies the point of concatenation. The two heterodyne sequences A and B binding free energy ΔGAB can be calculated by the equation (Bernhart et al., 2006),

Now the binding free energy for the interaction into the selected siRNA with its consensus target is

The binding free energy (ΔGAB) of the above four siRNA meets the qualified range.

Off-target silencing ability of the proposed siRNA was checked using the NCBI blast program and there is no effects were identified. But in vivo experimentation is mandatory for the measurement of the immune response activity which is the most important challenge for siRNA designing. The MSA result showed a single code mutation in a single strain for 3rd siRNA the target locality in the N gene cds identified by siDirect server but the other three 1st, 2nd, and 4th were completely conserved within 270 different sequences (Supplementary file 2). Also, 1st siRNA showed very high efficacy, which is greater than 2nd and 4th. Its predicted GC content percentage (38 %) and free binding energy (-31.50Kcal/mol) were in the excellent range (Oany et al., 2015). On top of all the analysis, it can be decided that the proposed 1st siRNA molecule meets all the desired criteria to be considered a potential siRNA and might play an important role to combat SARS-CoV-2. These results support the possible utilization of the anticipated novel siRNA hostile to the 270 strains of Bangladesh considered in our study. Scientists are trying to show outstanding collective effort due to characterize SARS-CoV-2 and recognize effective treatment as well as therapeutic alternatives for this rapidly proliferative, highly infective and potentially death causing virus (Iacob and Iacob, 2020). Multiple tentative drugs having viral defeated potentiality have been introduced since the beginning of this virulent disease. Nevertheless, effective molecular agents is not identified currently for the treatment of the increased number of people because of their lack of an appropriate medical response (Iacob and Iacob, 2020; Chowdhury et al., 2021), especially in Bangladeshi people. Although this study is not the first molecular therapeutic approach for treating SARS-CoV-2, our predicted siRNA can be the most suitable molecular therapeutic approach after passing experimental validation.

The last years, many studies have been revealed that explain various mechanism of in vivo uses of siRNA in systemically or locally. The most significant numerous articles focus on using unchanged siRNAs administration through intravenously by hydrodynamic transfection. Intravenous and Intranasal unmodified siRNAs delivery into lung targeting nucleoprotein of Influenza virus and SARS corona virus gives relief from Influenza virus infections (Tompkins et al., 2004b) and SARS corona virus fever (Li et al., 2005) in vivo, which represents strong protection against these lethal viruses.

In recent times, application of polyethyleneimine has been extended towards the complexation and delivery of RNA molecules, especially siRNAs (Urban-Klein et al., 2005). While chemically unchanged RNA molecules are very unstable with having high degrading possibility, but the enzymatic or non-enzymatic degradation is completely protected by PEI complexation. Before or after influenza virus infection, PEI promoted siRNAs' delivery into the lungs through IV route decrease virus production (Ge et al., 2004). In vitro activity also proved for siRNA nanoplexes which prevents siRNAs from serum devastation (Schiffelers et al., 2004). Polyethyleneimine and siRNA complexes are effectively delivered into particular cells in vitro, specifically low molecular weight PEI exhibit high siRNA protection and delivery efficacies (Werth et al., 2006). Lyophilization preserves physical stability and biological activity of the PEI/RNA (siRNA or ribozyme) complexes under some conditions (Werth et al., 2006; Brus et al., 2004).

Conclusions

Molecular therapy is replacing conventional therapeutic approaches as a promising alternative way because of its specificity and successive nature. Potentiality of siRNAs have already afford to effectual cure for different diseases. In this study, we have anticipated that a siRNA molecule to apply as therapeutic agent to combat SARS-CoV-2. Functionality and stability of the predicted siRNA molecule is also defined using various algorithms with their suitable parameters. Molecular therapies are precise to DNA sequences that are sometimes nonfunctional against different strains. Our designed siRNA is expected to overcome this issue as it has targeted a conserved sequence found in different SARS-CoV-2 strains, which Bangladeshi scientists sequence. While proper in vitro and in vivo support is still mandatory, even so we expect that this siRNA molecule will afford a practical treatment method hostile to the targeted Bangladeshi SARS-CoV-2 strains.

Funding

No funding for this study.

Author statement

All authors are equally contributed and approved it for Publication.

Declaration of Competing Interest

None.