Design of Potential RNAi (miRNA and siRNA) Molecules for Middle East Respiratory Syndrome Coronavirus (MERS-CoV) Gene Silencing by Computational Method.

The Middle East respiratory syndrome coronavirus (MERS-CoV) is a virus that manifests itself in viral infection with fever, cough, shortness of breath, renal failure and severe acute pneumonia, which often result in a fatal outcome. MERS-CoV has been shown to spread between people who are in close contact. Transmission from infected patients to healthcare personnel has also been observed and is irredeemable with present technology. Genetic studies on MERS-CoV have shown that ORF1ab encodes replicase polyproteins and play a foremost role in viral infection. Therefore, ORF1ab replicase polyprotein may be used as a suitable target for disease control. Viral activity can be controlled by RNA interference (RNAi) technology, a leading method for post transcriptional gene silencing in a sequence-specific manner. However, there is a genetic inconsistency in different viral isolates; it is a great challenge to design potential RNAi (miRNA and siRNA) molecules which can silence the respective target genes rather than any other viral gene simultaneously. In the current study, four effective miRNA and five siRNA molecules for silencing of nine different strains of MERS-CoV were rationally designed and corroborated using computational methods, which might lead to knockdown the activity of virus. siRNA and miRNA molecules were predicted against ORF1ab gene of different strains of MERS-CoV as effective candidate using computational methods. Thus, this method may provide an insight for the chemical synthesis of antiviral RNA molecule for the treatment of MERS-CoV, at genomic level.

Introduction

Coronaviruses are enveloped single-stranded positive-sense RNA viruses with genomes of 25–32 kb, and the group includes the biggest known genomes among the RNA viruses [1, 2], which infect and cause disease in a wide variety of species, including bats, birds, cats, dogs, pigs, mice, horses, whales and humans. Five types of coronaviruses are found to be circulating the human population [3]. A second coronavirus known as the Middle East respiratory syndrome coronavirus (MERS-CoV) has been recognized as the causal organism of a highly lethal pneumonia in patients in the Middle East and in travellers from this region about 10 years after the severe acute respiratory syndrome became epidemic [4]. The virus MERS-CoV belongs to the genus betacoronavirus [5] which is also known as novel coronavirus (nCoV). After the first report of causing human infection was published in September 2012, the World Health Organization (WHO) International Health Regulations Emergency Committee denied the fact that MERS-CoV did not meet the criteria of becoming “public health emergency of global concern” in July 2013. Nevertheless it was termed “serious and great concern” and was also reported as a threat to the entire world. In June, 2012, for the first time virus was identified in a patient who died in Saudi Arabia after developing acute respiratory distress syndrome (ARDS) and multiple organ dysfunctional syndrome (MODS) [6]. In same year, the first MERS disease broke out has been reported [7]. Although the name suggests that it is limited to the Middle East, the death scenario has been also reported in France, UK, Italy and Tunisia [8]. The symptoms include fever, cough and shortness of breath. Although MERS-CoV patients have primarily had respiratory disease, a number of secondary complications have also been informed, including acute renal failure, multi-organ failure, acute respiratory distress syndrome (ARDS) and consumptive coagulopathy. In addition, many patients have also reported gastrointestinal symptoms, including diarrhea [7]. As mentioned earlier, this virus has the capacity to transmit from human to human. It was also shown that MERS-CoV was capable of infecting bats and porcine cell lines in addition to human cells. These possessions would indicate a low barrier for transmission between hosts [9-11]. Recently, dipeptydyl peptidase 4 (DPP4) of primary human bronchial epithelial cells and huh-7 cells are identified as the receptor of MERS-CoV virus entry. After entering into the cell, the virus specially modulates the innate immune response, antigen presentation and mitogen-activated protein kinase. The genome structure of MERS-CoV is similar to other coronaviruses, with the 5′ two-thirds of the genome encoding the nonstructural proteins (NSPs) required for viral genome replication, the remaining 3′ third of the genome encoding the structural genes that make up the virion (spike, envelope, membrane and nucleocapsid proteins) and four accessory genes intermingled within the structural gene region. The viral genome expression starts with the translation of two large replicative polyproteins, pp1a and pp1b, which are encoded by the viral replicase gene that compromises ORF 1a and 1b. This replicase complex not only mediates genome replication but also causes transcription of nested set of subgenomic mRNA responsible for structural proteins [12, 13]. After transcribing pp1a, −1 ribosomal frame shifting is occurred to transcribe pp1b which is subsequently cleaved into 15 or 16 nonstructural proteins [12]. The downstream region of ORF1b is characterized by containing a variable number of smaller genes, always including those encoding the spike (S), envelope (E), membrane (M) and nucleocapsid (N) structural proteins. These genes are translated from subgenomic (sg) mRNAs that form a 5′ and 3′ co-terminal nested set with the viral genome. Subgenomic mRNAs are composed of a common 5′ leader sequence that is identical to the genomic 5′ region and a variable part of the 3′ quarter of the genome, with different sg mRNAs making different ORFs available for translation. The complement of the leader and “body” segments of the sg mRNAs are assumed to be joined during discontinuous negative-strand RNA synthesis. This step produces the templates for sg mRNA synthesis and is directed by a base-pairing interaction between conserved transcription regulatory sequences (TRSs) [14-16]. Such TRSs are found at the 3′ end of the leader sequence (leader TRS) and at different positions upstream of genes in the genomic 3′ proximal domain (body TRSs). The synthesis of subgenome-length negative-stranded RNAs is directed by the complement of a body TRS at the 3′ end of a nascent minus-strand base pairing with the leader TRS, with the extent of sequence complementarity being an important determinant of the level at which a given sg mRNA is produced [13]. In this study, two types of gene silencing techniques were implemented computationally, siRNA and miRNA. siRNA also called silencing RNA is generally a double-stranded RNA and contains 21 base pairs in length. The central role of siRNA is posttranscriptional gene silencing by interfering with the protein, and their basic mechanism is to bind and promote the degradation of messenger RNA (mRNA) at specific sequences [17]. An effective siRNA has to contain some basic features like A/U at the 5′ end of the antisense strand, G/C at the 5′ end of the sense strand, at least 5 A/U residues in the 5′ terminal end or one third of the antisense strand and absence of any GC stretch of more than 9 nucleotide in length. If designed siRNA fails to fulfill above criteria, then gene silencing is not possible by that particular siRNA [18]. Each siRNA strand contains a 5′ phosphate group and a 3′ hydroxyl (OH) group, produced from dsRNA or hairpin looped RNA. By using nanoparticles like polyethyleneimine and other approaches, siRNA can be introduced into the cell for specific knockdown of a gene of interest, after entering into the cell siRNA cleaved by an enzyme called dicer (RNase III-like enzyme). Then the siRNA incorporated with multi-subunit protein complex called RNAi-induced silencing complex (RISC). RISC seeks out the correct target mRNA, and then siRNA unwinds and is believed. After this process, the antisense strand directs degradation of the complementary strand of appropriate mRNA using different kinds of endo and exogenous enzymes [19]. siRNA may also be used for therapeutic purpose like vaccine, and their efficacy can be ensured to more than 75 % to knockout or kill any targeted gene or protein [20]. Mature miRNA regulates diverse cellular processes in plant and animal species through the assembly of miRNA-induced RISC complex which binds to the complementary targets in the cellular mRNA and subsequently cleaves or transcriptionally represses the mRNA depending on the complementarity. In addition, recent studies suggest that miRNA also play a role in the regulation of viral infection [21]. Therefore, in the present study an effort has been made to identify potential siRNA and miRNA molecules for silencing of ORF1ab coding mRNA or gene in MERS-CoV using computational approach.

Materials and Methods

Sequence Retrieval and Analysis

Nine genes of ORF1ab of nine different strains of MERS-CoV were collected from the gene bank database available at National Centre for Biotechnological Information (http://www.ncbi.nlm.nih.gov/). This gene has the crucial role for further down regulation of other necessary genomic and structural replicator components. The accession numbers of these nine genes are JX869059.2, KF192507.1, KC667074.1, KC164505.2, KC776174.1, KF186567.1, KF186566.1, KF186565.1 and KF186564.1 for different strains of MERS-CoV.

Open Reading Frame (ORF) Screening and Multiple Sequence Alignment

With a view to screening out the coding regions of these nine different gene sequences of MERS-CoV, a web-based interactive tool for predicting and analyzing ORFs from nucleotide sequences was used known as ORF Finder (http://www.ncbi.nlm.nih.gov/projects/gorf/). Afterward, using CLASTALW (http://www.genome.jp/tools/clustalw/), individually all the coding regions for nine different strains were analyzed. The analysis showed 99.99 % similarities among the coding regions.

Potential Designing of siRNA and miRNA

To design a potential siRNA, siDirect 2.0 was used [22]. While using this tool, some rules were marked as parameters including Ui-Tei, Amarzguioui, Renold rules as well as melting temperature which must be below 21.5 °C to allow potential siRNA duplex. Apart from these, other factors were also kept in mind on the concept of algorithms (Table 1). In similar fashion, for the purpose of designing a potential miRNA, another web-based tool called Integrated DNA Technology (IDT), (http://sg.idtdna.com/Scitools/Applications/RNAi/RNAi.aspx) was utilized. Following this, the RNAis were submitted to RISC Binder for the prediction of functional guide strand of miRNA [23].

Table 1

Algorithms or rules for rational design of siRNA molecules

Ui–Tei rules	Amarzguioui rules	Reynolds rules
A/U at the 5′ terminus of the sense strand	Duplex end A/U differential >0 strong binding of 5′ sense strand	Each rule is assigned a score which is summed up to a total duplex score to improve the efficiency of siRNA
G/C at the 5′ terminus of the antisense strand	No U at position 1. Presence of A at position 6
At least 4 A/U residues in the 5′ terminal 7 bp of sense strand	Weak binding of 3′ sense strand. No G at position 19
No GC stretch longer than 9 nt

Similarity Search and Target Alignment

For checking any off-target sequence resemblance in other nontargeted organism’s genome, BLAST tool (http://www.ncbi.nlm.nih.gov/blast) was used against whole Genebank database by applying expected thresholds value 10 and BLOSUM 62 matrix as parameter.

GC Calculation and Secondary Structure Prediction

DNA/RNA GC content calculator (http://www.endmemo.com/bio/gc.php) was used so as to calculate the GC content of the predicted siRNA. Consequently, for prediction of secondary structure and for computation of the free energy of folding for both siRNA and miRNA, mfold server [24] was used.

Calculation of RNA–RNA Interaction Through Thermodynamics

To study the thermodynamics of interaction between predicted siRNA and target gene, RNAcofold program [25] was used. It calculates the hybridization energy and base-pairing form of two RNA sequences. It functions as extension of McCaskill’s partition function algorithm to compute probabilities of base pairing, realistic communication energies and equilibrium concentrations of duplex structures. Similar interaction study was performed for miRNA with RNA hybrid program [26].

Calculation of Heat Capacity and Concentration Plot of Duplex

For both cases (miRNA and siRNA), the heat capacity plot and concentration plot have been calculated. The ensemble heat capacity is plotted as a function of temperature, with the melting temperature indicated. The detailed heat capacity plot also shows the contributions of each species to the ensemble heat capacity, whereas concentration plot—(Conc.), the point at which the concentration of double-stranded molecules of one-half of its maximum value defines the melting temperature (Conc.). Both of these were performed by DNAMelt Web server [27].

Prediction of Secondary Structure

The mfold Web server contained UNAfold program [24] which was also used to predict the secondary structure formed between mRNA and RNAi guide strands.

Results and Discussion

ORF1ab gene is the first line open reading frame genes in its whole genome length in MERS-CoV which covers two-third of length. In addition to its gene, its translational products also play active part in different functional aspects such as viral genome regulating transcriptase and protease and also encode replicase machinery for own replication of virus [12]. So to speak, it can be considered as an immediate early gene and silencing of this gene will result in termination of its essential function. In accordance with the present study, nine ORF1ab gene sequences from nine different strains of MERS-CoV were obtained from gene bank of NCBI. All sequences were subjected in ORF Finder to determine the cds of genes. ClastalW was brought into service in order to find out similarities among the cds. Results illustrate that they were more or less similar (99.99 %). In next step of the study, siDirect, a web-based online software system for computing highly effective siRNA, was used. siDirect 2.0, an updated version of siDirect, provides functional and off-target minimized siRNA design. All cds examined individually for siRNA against themselves through siDirect 2.0. As Ui–Tei, Amarzguioui and Reynolds rules (U, R, A rules) govern siRNA sequence preference and in accordance with which highly effective siRNA can be readily designed. 469 siRNA suggested by siDirect 2.0 which were tested for whether they follow URA rules. 93 siRNA found to follow URA rules. Similarly for miRNA against same genes, all sequences were posed to IDT where 45 RNAi found to meet the standards of IDT. As it is important to ensure whether RNAi binds RISC, results presented in Table 3 prove the successful incorporation of guide strand with RISC. In this study, RISC Binder tool was used to confirm the incorporation of miRNA guide strand with RISC. Binding ability of RNAi with RISC proves their capability of working as miRNA [23]. Only four miRNA from forty-five persists after the selection. ClastalW results of four miRNA revealed that three miRNA holds identical sequence based on which two groups were maiden. One group contains identical three miRNAs and the remaining one subjected to another group. For further validation of reduced off target, all siRNAs and miRNAs were subjected to BLASTN. All the siRNAs and miRNAs hold GC content within 36–48 % determined by GC content calculator. GC content of siRNA and miRNA is a determinant of the perfection of their business as GC-rich target sites were avoided because of their potential to be involved in strong bonding in their secondary structure while designing siRNA and miRNA [28]. Lower GC content leads to stronger inhibitory effects, and 30–57 % of GC content is considered sufficient for the execution of RNAi’s action [29]. Because function of RNAi is determined by molecular structure, there is a significant investment in computational methods for predicting RNA secondary structure [30] and these computational methods for modeling RNA secondary structure have proven to be valuable in many cases in which crystal structures are not available [31]. Minimum free energy (MFE) is considered as benchmark of RNA’s molecular structural accuracy [32]. Free energy minimization is a long-established paradigm in computational structural biology that is based on the assumption that, at equilibrium, the solution to the underlying molecular folding problem is unique and that the molecule folds into the lowest energy state [31]. To measure the stability of structure of the guide strand, the minimum free energy of folding of was computed using mfold and those who possess least minimum free energy was tabulated (Tables 2, 3, 4). Therefore, negative free energy of folding may have lower accessibility for target site. Reliable prediction of RNA–RNA binding energies is crucial, e.g., for the understanding on RNAi, microRNA–mRNA binding and antisense interactions. The thermodynamics of such RNA–RNA interactions can be understood as the sum of two energy contributions: (1) the energy necessary to “open” the binding site and (2) the energy gained from hybridization [33]. Hence the accurate computational treatment of RNA–RNA binding therefore lies at the heart of target prediction algorithms. The RNAcofold Web server was used to figure the hybridization energy and base-pairing pattern of two RNA sequences. RNAcofold, therefore, guides the structural prediction process by a combination of thermodynamic and kinetic considerations. In addition to GAB, the free energy of the heterodimer of sequence A and sequence B can be calculated using the equation Gbinding = GAB − GA − GB [25]. So, for the interaction of a target mRNA and its predicted siRNA, free energy of binding is Gbinding = GAB − GA − GB (Table 2). Free energy of binding for miRNA was determined through RNA hybrid. siRNA and miRNA having highest energy of binding with target sequence are shown in Tables 2, 3 and 4. The DINAMelt Web server predicts entire equilibrium melting profiles as a function of temperature with melting temperature for a hybridized pair of nucleic acids. The heat capacity as a function of temperature is given in Tables 2, 3 and 4. (Conc.) is the temperature at which the temperature of nucleic acid reaches half of its maximum value [27]. It is known that the secondary structure can indeed cause inhibition of RNAi [34-36]. Also the secondary structure of RNAi can be involved in the RISC cleavage. That is why predicting secondary structures in order to guide selection of sites for RNAi targeting is more important [28]. Secondary structures of RNAi provide a convenient and computationally tractable estimation not only to structures but also to the thermodynamics of RNA–RNA interaction. UNAfold predicts the most stable secondary structure of an oligonucleotide by minimizing folding free energy. The folded structures are predicted at the specified temperature which is 37 °C. The secondary structures of siRNA and miRNA are shown in Figs. 1a–e, 2a–c and 3. After analysis of all siRNA and miRNA, finally five siRNA and four miRNA pass all barriers of effectiveness and those are worthy to fight against MERS-CoV. RNA interference (RNAi) is a novel regulatory and powerful silencing approach for molecular therapeutics [37]. Though both siRNA and miRNA have novel function for gene silencing, siRNA is comparatively better than miRNA due to its sequence specificity with targets and easier insertion techniques (via nano-particle, lipid based transfection, etc.) into the cell. Imbued by the mechanism and suitability of RNAi, scientists and pharmaceuticals are working to use it as next-generation therapeutics. To treat a wide variety of disease, RNAi technology is the riding lamp to researchers. Current discoveries on miRNAs’ involvement in cancer suggest that RNA inhibition can be used as an advanced mean for treating cancer patients [38]. miRNA-based therapeutics are also in consideration for the treatment of cardiovascular disorders [39].

Table 3

Effective miRNA molecule with GC%, free energy of folding and free energy of binding with target (Group 1)

Target no.	Location of target within mRNA	miRNA target within mRNA	Predicted miRNA duplex candidate at 37 °C	GC% (%)	Free energy of folding with target	RB score	mfe (kcal/ mol)	T_mC_p (°C)	T_mConc. (°C)
01	18,483–18,509	TGCTACTAATCACGATCGTTATTGCTC	CUACUAAUCACGAUCGUUAUUGCTC	40	0.9	+	−18.9	67.0	66.5
			GAGCAAUAACGAUCGUGAUUAGUAGCA
02	4717–4743	AAGTGGAAGATGGTTGTGTGTGATAAG	GUGGAAGAUGGUUGUGUGUGAUAAG	44	0.3	+	−7.8	54.3	57.7
			CUUAUCACACACAACCAUCUUCCACUU
03	3652–3678	AAGTGCTATAAGGCTATGAATGCATAT	GUGCUAUAAGGCUAUGAAUGCAUAT	36	1.0	+	−17.5	69.1	65.2
			AUAUGCAUUCAUAGCCUUAUAGCACUU

RB score RISC binder score

Table 2

Effective siRNA molecule with GC%, free energy of folding and free energy of binding with target

Target no.	Location of target within mRNA	siRNA target within mRNA	Predicted siRNA duplex siRNA candidate at 37 °C	GC% (%)	Free energy of folding with target	ΔGbinding (k/cal)	T_mC_p	T_mConc.
01	1169–1191	AGGGTGATAGCTGTTCTTCTAAT	UAGAAGAACAGCUAUCACCCU	42.85	0.5	−36.05	86.1	84.9
			GGUGAUAGCUGUUCUUCUAAU
02	5458–5480	GACCTATCTGCTTTCTATGTTAA	AACAUAGAAAGCAGAUAGGUC	38.09	0.8	−34.53	82.9	81.5
			CCUAUCUGCUUUCUAUGUUAA
03	9520–9542	GACACTTATGCAGCTCTTAGAAA	UCUAAGAGCUGCAUAAGUGUC	42.85	0.0	−35.31	84.8	83.7
			CACUUAUGCAGCUCUUAGAAA
04	14,284–14,306	TCCTCTAACGCTTTTCTTGATTT	AUCAAGAAAAGCGUUAGAGGA	38.09	1.0	−33.97	82.9	81.6
			CUCUAACGCUUUUCUUGAUUU
05	14,892–14,914	TTGCGTCATTGGTACTACAAAGT	UUUGUAGUACCAAUGACGCAA	38.09	0.8	−32.50	84.2	83.0
			GCGUCAUUGGUACUACAAAGU

Table 4

Effective miRNA molecule with GC%, free energy of folding and free energy of binding with target (Group 2)

Target no.	Location of target within mRNA	miRNA target within mRNA	Predicted miRNA duplex candidate at 37 °C	GC% (%)	Free energy of folding with target	RB score	mfe (kcal/mol)	T_mC_p (°C)	T_mConc. (°C)
04	6668–6694	ACCTAGGAATCTCTTCTGCTTGTGACG	CUAGGAAUCUCUUCUGCUUGUGACG	48	0.3	+	−11.7	71.5	73.0
			CGUCACAAGCAGAAGAGAUUCCUAGGU

RB score RISC binder score

Fig. 1

a–e Secondary structures of target–siRNA duplex

Fig. 2

a–c Secondary structures of target–miRNA duplex (Group 1)

Fig. 3

Secondary Structure of target–miRNA duplex (Group 2)

Conclusion

With RNAi technology, it is possible to design a number of siRNA and miRNA molecules for silencing of substantial genes in numerous biological systems. Their connections with target can also be calculated computationally. Hence, in our study five siRNA and four miRNA molecules were predicted against ORF1ab gene of different strains of MERS-CoV as effective aspirant using computational approaches. These molecules may be a key to a novel antiviral therapy against MERS-CoV. This research will open the milestone corridor in MERS therapy. The consequence of this study would also afford a foundation to the researchers and pharmaceutical industry to advance antiviral therapeutics at genomic level.