A bioinformatic approach of targeting SARS-CoV-2 replication by silencing a conserved alternative reserve of the orf8 gene using host miRNAs.

The causative agent of the COVID-19 pandemic, the SARS-CoV-2 virus has yielded multiple relevant mutations, many of which have branched into major variants. The Omicron variant has a huge similarity with the original viral strain (first COVID-19 strain from Wuhan). Among different genes, the highly variable orf8 gene is responsible for crucial host interactions and has undergone multiple mutations and indels. The sequence of the orf8 gene of the Omicron variant is, however, identical with the gene sequence of the wild type. orf8 modulates the host immunity making it easier for the virus to conceal itself and remain undetected. Variants seem to be deleting this gene without affecting the viral replication. While analyzing, we came across the conserved orf7a gene in the viral genome which exhibits a partial sequence homology as well as functional similarity with the SARS-CoV-2 orf8. Hence, we have proposed here in our hypothesis that, orf7a might be an alternative reserve of orf8 present in the virus which was compensating for the lost gene. A computational approach was adopted where we screened various miRNAs targeted against the orf8 gene. These miRNAs were then docked onto the orf8 mRNA sequences. The same set of miRNAs was then used to check for their binding affinity with the orf7a reference mRNA. Results showed that miRNAs targeting the orf8 had favorable shape complementarity and successfully docked with the orf7a gene as well. These findings provide a basis for developing new therapeutic approaches where both orf8 and orf7a can be targeted simultaneously.

Introduction

The causative agent of the 2019 global COVID pandemic that has engendered severe public health consequences [1,2] and claimed over 5.3 million lives as of December 10, 2021 [3], is the seventh coronavirus that has infected the human race [4]. This virus, known as the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), was first reported to have infected humans in Wuhan, Hubei province, China [1,2]. Ever since the first reported case in December 2019, the virus has yielded relevant mutations [[5], [6], [7], [8], [9], [10], [11], [12]] resulting in multiple variants of interest and variants of concern [13]. The evolutionary and functional significance of many of these variations are still under research, however, some viral components responsible for carrying out crucial interactions with the host cell, are much more susceptible to mutations [14]. One such component, highly susceptible to deletions and nucleotide substitutions [15], is the non-structural protein orf8, coded by one of the most variable accessory gene [15,16]⁠. Although dispensable for viral replication [17]⁠, this protein play a role in interactions with the infected host cell by modulating vesicular trafficking, viral particle packaging as well as by modifying the host innate immunity [18].

The orf8 gene of SARS-CoV-2 comprises of 366 nucleotides spanning from position 27,894 to 28,259 of the viral genome. These nucleotides encode for the 121 amino acid long orf8 protein (NCBI reference sequence NC_045512.2; Gene ID 43740577; Protein ID YP_009724396.1) [19] which takes up the role of a central organizer of the virus-host hybrid network [20]. Evolution of the SARS-CoV-2 genome, since its initial human-to-human transmission has caused divergence of the virus, globally, into at least three major phylogenetic types. One of these diverged groups emerged after single point mutations at positions 28077 and 28144 of the orf8 gene, causing V62L and L84S substitutions respectively [[21], [22], [23], [24]]. In addition, a H112Q substitution has also been reported [25]⁠. Apart from single point alterations, the orf8 gene has undergone multiple deletions resulting in frame-shift and substitution mutations [19]⁠. Among those deletions is a 382-nucleotide deletion between positions 27848 and 28229 of the genome (Δ382) that was first detected in clinical specimens from Singapore and Taiwan [26]⁠. The emergence of this mutation was traced back to Wuhan at the beginning of the pandemic in mid-December 2019 [26,27]⁠. The outcome of this mutation is the deletion of 40 nucleotides from the end of orf7b gene sequence and 336 nucleotides from that of orf8 region including the removal of a regulatory element from the gene. This ultimately leads to an inhibition of the orf8 expression [27]⁠ without however affecting the replicative fitness of the virus [26,28]. Rather, it was observed that, the Δ382 orf8 variant caused milder disease symptoms with a later onset, a lower probability of developing hypoxia and improved clinical outcome when compared to the wild type strain of the virus. Since orf8 plays a role in host innate immunity evasion, the reason behind a better disease outcome of the Δ382 variant could be a more enhanced host immune response [28]⁠. Other similar deletions in the orf8 gene were later reported in samples from Bangladesh (Δ345), Australia (Δ138), and Spain (Δ62) [27,28]⁠. Among recently emerging alterations, the B.1.1.529 variant named ‘Omicron’, has been designated by the World Health Organization (WHO) as a Variant of Concern (VOC) [29] which is still under investigation. The complete molecular functioning and functional significance of orf8 is still under research, however, studies reported so far support the hypothesis that the truncated orf8 gene has proven to be advantageous in the adaptation of SARS-CoV-2 to humans [30]⁠.

The orf8 of SARS-CoV-2 shows a significant divergence with less than 20% sequence similarity with the orf8 of SARS-CoV [31]. In fact, there have been evidences to showing that in the viral lineage, this protein is an evolutionary hot spot [27,32]⁠. Despite being one of the most variable regions of the viral genome [16]⁠, protein interaction mapping has pointed out orf8 as the most connected hub with 47 interactions with the host cell [18]⁠ (https://amp.pharm.mssm.edu/covid19/genesets/22). Hence, understanding all the functional attributes of orf8 of the SARS-CoV-2 might give us an edge in clarifying relevant significant virus-host interactions which can be utilized for developing therapeutic targets. As of now it is known that, in addition to its interactions with the host immune system [[33], [34], [35]]⁠, the non-structural accessory protein, orf8 of SARS-CoV-2, interacts with the endoplasmic reticulum (ER) of the host cell [36]⁠. Despite being significantly different from its counterpart in SARS-CoV, the orf8 of both the virus contains a signal sequence for ER import [36]⁠ and it is conjectured that instead of being retained, this protein is rather secreted out of the ER [37]⁠. SARS-CoV-2 orf8 variants are also responsible for inducing ER stress in the host cell by the activation of inositol-requiring enzymes 1 (IRE1) and activating transcription factor 6 (ATF6) pathways [34]⁠.

When it comes to defense against viral infections, the host innate immunity plays a crucial role [38]⁠. Several transcription factors like NF-κB and IRF-3 [39,40], start stimulating the expression of type-I interferons (IFN-α/β) [41], which bind to its respective receptors and activate the JAK/STAT pathway [39,40]. In addition, the interferons cause the nuclear translocation of IFN-responsive transcriptional factors. These factors ultimately manifest an antiviral state by causing the activation of a set of interferon-stimulated genes (ISGs) which contain IFN-stimulated response elements (ISREs) in their promoters [33,39,40]. A part of the immune evasion strategy adopted by SARS-CoV-2, is preventing the manifestation of this antiviral state by inhibiting the expression of type-I interferon [33,34], Recent studies have demonstrated that along with some other viral components of SARS-CoV-2, the orf8 protein acts as an inhibitor of ISGs [34], he IFN-β promoter, ISRE promoter and the NF-κB element [33]⁠, which is a clear indication of its involvement in the suppression of the type-I interferons.

A rather unique feature exhibited by the orf8 of SARS-Cov-2 but not by its counterpart in SARS-CoV [35], is the downregulation of MHC-1 expression [35,42]. The major histocompatibility complex class 1 (MHC-1) molecules are a group of cell surface proteins that play a crucial role in the presentation of antigenic peptides on the surface of infected host cells [43]. The MHC-1-peptide complex is further recognized by the T cell receptor on CD8+ T cells, which then initiate an immune response involving the release of multiple toxic substances (like perforins, FasL, and granzyme) and cytokines like IL-2, TNF-α, and interferon-γ [44,45]. Cell death, induced as a consequence of this response by the CTLs, leads to the elimination of the virus infected cells from the host system, resulting in a successful prohibition of the spread of viruses [45]. Hence, suppression of the MHC-1 molecules by the orf8 of SARS-CoV-2 completely conceals the virus in the host cell and protects it from identification by the host immune system [35]. This unique mechanism which was previously not detected in SARS-CoV [4,35]⁠, has given an evolutionary advantage to the COVID-19 virus, enabling it to survive better in the human host cells.

Acquiring different strategies for immune evasion is a well-known feature of viruses that are known to exploit the host system for its survival [46,47]. Previous studies suggest that in the case of SARS-CoV-2 the protein mediating such evasion by camouflaging the virus in the host cell is the orf8 protein [37,48,49], making it a potential target for antiviral therapies. Multiple evidence is suggestive of the fact that therapies, designed using host encoded micro-RNAs can serve as potential antiviral remedies since they are known to modulate viral pathogenesis and replication [[50], [51], [52], [53]]. miRNAs do so by interacting with the viral mRNA by binding to a complementary region, thereby silencing its expression [[54], [55], [56], [57]]. Even in the case of SARS-CoV-2, it is possible to subdue the viral replication and other virus-host interactions by using target specific miRNAs from the human genome [58]⁠. Given the functional attributes of orf8, silencing its expression using miRNA mediated gene regulation can potentially improve the efficiency of host immune surveillance [49,56].

While the reasons for targeting the orf8 region are quite evident, its versatility and high mutation rate in emerging strains [30,59] have made the identification of a specific target sequence quite challenging. Studies however have pointed out comparatively conserved regions of the SARS-CoV-2 genome including the orf7a gene [60,61], that is located before the orf8 and orf7b regions [27,62]. The Δ382 mutant of this virus exhibits a deletion in the ORF7 gene where only the orf7b segment has deletions with no changes in the orf7a segment [27]. This region is mostly conserved not only in different strains of SARS-CoV-2 but also in its previous homologs as determined from multiple phylogenetic analyses [63]. Functional studies on this gene have shown its involvement in the host immune system modulation via induction of a cytokine storm [64] and also by its role as an immune antagonist [[63], [65]]. Adding to its functional versatility, the orf7a protein also contains an immunoglobulin (Ig)-like fold [[63], [66]] which has been maintained in different lineages of SARS-CoV-2 [[63], [65]]. This Ig-like domain manifests an integrin binding site which further contributes to multiple virus-host immune interactions [[66], [67]]. What makes the orf7a gene more interesting is the fact that it might share a common origin with the orf8 gene of SARS-CoV-2, as pointed out by gene profile comparison studies [63,68,69]. Additionally, the orf8 protein also contains an Ig-like domain with a core fold similar to that of the orf7a protein [31]. It has been observed that emerging strains of the SARS-CoV-2 virus exist with deletion of the orf8 region, but conservation of the orf7a gene with which orf8 share a sequence homology [[27], [63]]. This peculiar course of evolution and functional similarity between orf7a and orf8 adopted by the SARS-CoV-2 virus hints towards a hidden viral strategy, targeting which might be beneficial for controlling the viral growth [63].

In this study, we have used computational methods to first check the binding affinity of specific miRNAs to different orf8 sequences. After a successful molecular docking, we proposed a hypothesis based on the sequence homology between the orf8 and orf7a genes. The relative evolutionary conservation of the orf7a gene and its functional similarity with both of the truncated and reference orf8 sequences of SARS-CoV-2, led us to hypothesize the following: the orf7a gene might be a similar alternative reserve for orf8 that was compensating for the deleted gene. Hence, there is a possibility that the same therapeutic strategy could target both these sequences simultaneously to impact the virus life cycle.

Materials and methods

Data mining

Viral genome (Accession number MT259249 and NC_045512.2) data was obtained from NCBI. Sequence of nine variants reported from different countries, of Clade G (GR2, G2, GH2, GK1, GH2) were downloaded along with sequence for the new B.1.1.529 variant named ‘Omicron’. (Table 1 ). The primary sequences of orf8 and orf7a were extracted from the whole genome sequences.

Table 1

List of reference and variants of SARS-CoV-2 used in this study.

Sl. No	Variant	Accession number	Country	Clade
1.	–	NC_045512.2	Wuhan	Reference
2.	Omicron	EPI_ISL_6959868	Omicron	GRA
3.	Beta	EPI_ISL_3299558	South Africa	GH
4.	Lambda	EPI_ISL_3298366	Peru	GR
5.	Kappa	EPI_ISL_3277130	India	G
6.	Eta	EPI_ISL_3278349	UK/Nigeria	G
7.	Iota	EPI_ISL_3298341	USA	GH
8.	Zeta	EPI_ISL_3266119	Brazil	GR
9.	Delta	EPI_ISL_3299048	India	G
10.	Gamma	EPI_ISL_3299430	Brazil	GR
11.	Alpha	EPI_ISL_3300249	UK	GR

Sequence alignment

Multiple sequence alignments were performed with the help of the CLUSTAL Omega version 1.2.4 among 11 sequences including the reference sequence of orf7a and orf8 region (reference and 9 variants) (Table 1). Clustal Omega is a multiple sequence alignment tool that generates alignments between three or more sequences using seeded guide trees and HMM profile-profile algorithms. Divergent sequences are aligned in biologically relevant multiple sequence alignments [70].

Screening of human miRNAs targeting orf8 region

Mature human miRNAs were screened using the sequences of orf8 region of reference sequences and variants. The “custom prediction” application of miRDB (http://mirdb.org/custom.html) [71] online server was used for the screening of miRNAs for 8 sequences (including the orf8 reference sequence and 7 variant sequences). Due to the high amount of nucleotide deletion in the variants EPI_ISL_3298366 and EPI_ISL_3299048, these two sequences were excessively short to be processed by the miRDB server. Hence, for these variants, the miRNAs with the highest target score were screened from the miRBase online server.

miRNA and mRNA preparation and duplex formation

The sequences of the orf8 reference gene, orf8 variants (Refer to Table 1 for sequence IDs), and the orf7a gene (NCBI Reference sequence: NC_045512) were obtained from the NCBI database. Screened miRNA sequences were obtained from the miRDB database (Refer to Table 2 for the list of miRNAs). The build and edit nucleic acid tool of Digital Studio v21.1.0.20298 was used to upload the sequence of the sense strand of both the mRNAs and miRNAs separately. The pdb file for the sense strand was saved for each case.

Table 2

Table showing list of miRNAs targeting reference orf8 region with respective target scores.

Accession number	Target Score	miRNA (hsa-miR)	Role
NC_045512.2	93	8073	Identification of novel serum microRNAs in sepsis patients [91]
	93	221-5p	Inhibits INF beta expression [97], inhibits porcine epidemic diarrhoea replication by targeting viral genome through NF-kβ pathway [95], suppresses osteogenic differentiation by smad 3 targetting [100], correlates with increased heart failure risk in coronary patients [99], miR 221 repression maintains the differentiation of cells [96], inhibits interleukin-6 expression in human colonic epithelial cells [101], Serum marker for lung cancer detection [98].
	84	145-3p	can directly inhibit S protein expression and SARS-CoV-2 replication [105], Involves in macrophage polarization via IL-16 [104], a Novel biomarker for an acute myocardial infarction [103], a tumor suppressor in lung adenoma carcinogenesis [102].
	82	758-5p	upregulation of miR-758 expression by HCV as a novel mechanism contributing to downregulation of TLR3 and TLR7 in patients with HCV infection [107]. Regulates cholesterol uptake by targeting CD36 3′ UTR [108], suppresses Glioblastoma proliferation [106].
	79	5047	Downregulates cervical cancer [110], upregulates HNSCC [109].
	79	570-5p	Upregulated in CSF in neurosyphilis patients [111].
	79	548ba	Expression in human ovarian granulosa cells [112].
	79	548ai	Differential expression in prostate cancer [113].
	79	548ag	Lower levels in stored platelets [89].
EPI_ISL_3277130	94	8073
	94	221-5p
	84	570-5p
	84	548ba
	84	548ai
	84	548 ag
EPI_ISL_3278349	65	1537-5p	activation of platelets [90]
	94	8073
	94	221-5p
	84	570-5p
	84	548ba
	84	548ai
	84	548 ag
EPI_ISL_3299558	81	3059-3p	Targeting gene MORC3. May required for the transcription of influenza A virus during infection [92]
	78	194-5p	negative regulation of interleukin-10 production
EPI_ISL_3298341EPI_ISL_3300249	95	1226-3p	Targets Vasohibin-1 gene
	91	4684-3p	unknown
EPI_ISL_3266119EPI_ISL_3299430	92	219b-3p	Targets PACRG (Parkin Coregulated). A protein coding gene positively regulate TNF signaling [93].
EPI_ISL_3299048EPI_ISL_3298366	60	888-5p	Targets RC3H1 that can modulate the activity of the IKK/NF-κB pathway [94].
	60	3926	unknown
	60	7849-3p	activation of platelets [90]

To form the mRNA-miRNA duplex strand, the nucleic acid docking tool HNADock (http://huanglab.phys.hust.edu.cn/hnadock/) was used where both the nucleic acid molecules were uploaded in the pdb file format. For the docking of the reference orf8 gene with its respective screened miRNAs, the seed region was specified (Table 3 ) in the receptor-binding site residue(s) section of the server. The job was submitted, and after completion, the best models were selected for the reference docking while the top models were chosen in case of variants. All models were hierarchically arranged according to docking scores where top models had the highest negative score. Pdb files of the selected models were saved for further analyses.

Table 3

Selected reference miRNAs of orf8 on selected variants.

miRNA (HSA-MIR)	Variant ID (EPI_ISL)	Sequence Details
		Seed sequence	miRNA interaction sequence on variant orf8 region
145-5p	3299430	AGGAAUCA	AUAACACUUCAAGGUAUUGGGAA
5047		AGCUGCA	UUAUGAGGCUAUGUACACA
8073		UGCCAGGA	AACACUUCAAGGUAUUGGGAA
145-5p	3300249	AGGAAUCA	AUAAAAAUACAUGUGUUGGUAGCG
5047		AGCUGCA	AUAUAAAAAUACAUGUG
8073		UGCCAGGA	CAACAUGUGACUGGACAAAUGCUG
145-5p	3299558	AGGAAUCA	AAGUGAAAUCAUAGGAUACAAGG
5047		AGCUGCA	CUUUUCAAGUGAAAUCAUAGGAUAC
8073		UGCCAGGA	ACCUAGAGUUUUUAGUGCAGUUG

Preparing 3D folded structure of the nucleotide duplex

The duplex models chosen from HNADock were viewed on Digital Studio. The miRNA molecule was specified as the ligand from the ‘Receptor-Ligand Interactions’ section of the software and the ligand interactions were viewed and the specific region of interaction between mRNA and miRNA was identified. This sequence was then uploaded on the UNAFold web server (http://www.unafold.org/mfold/applications/rna-folding-form.php) for RNA folding and formation of the secondary duplex structures. From the output files, the structure with the highest negative ΔG value was selected and its respective Vienna file was downloaded to get the structure details which included a dot-bracket form of the 2D structure.

The next step was to generate a 3D structure of the duplex for which the 3dRNA tertiary structure prediction method (http://biophy.hust.edu.cn/new/3dRNA/create) was used. The duplex sequence and the dot-bracket format of the 2D structure were uploaded to the server and the job was submitted to form the 3D models using the ‘3dRNA-Lib1’ as the Template Library Version. The obtained results showed the top five predictions where the model with the lowest score was selected for each case and its pdb file was saved for further analyses.

Protein preparation and molecular docking

For docking analyses, the human AGO-2 protein was used [72] and its crystal structure (PDB ID: 4F3T; resolution: 2.25 Å) was retrieved from the Protein DataBank (PDB) in pdb format, for the study. The protein complex was viewed on Digital Studio and was prepared by deleting the ligand groups containing hetero atoms (hetatm) and water molecules [72].

The prepared hAGO2 protein and the 3D mRNA-miRNA duplex structures obtained from 3dRNA were considered for the docking analyses for each case. For an illustration of the catalytic binding between the AGO-2 protein and the tertiary duplex structures, the online web server PatchDock [[73], [74]] was used for performing a receptor-ligand molecular docking. The prepared protein was uploaded as the receptor molecule in pdb format while the saved 3D duplex models were specified as the ligand. To perform a successful protein-small molecular docking, the clustering RMSD value was set to 1.5 Å before the final job submission. Following the PatchDock algorithm, the solutions on the server are ranked according to the Geometric shape complementarity score. This score provides information about the shape complementarity which is an important feature of the surface interface of the interacting biological molecules and is correlated with molecular associations like non-polar desolvation and van der Waals interactions [75]. The top-ranked model with the highest score was selected for interpretation and viewing using Digital Studio.

The pdb file of the final protein-duplex docking model was viewed on Digital studio where the mRNA-miRNA duplex was specified as the ligand to highlight the specific ligand interactions. The molecular complex was properly oriented and the final image was extracted by setting the pdb file to a publication-quality script on Digital Studio following which optimization was done using the Illustrator platform (https://www.adobe.com/in/products/illustrator/free-trial-download.html).

Results

The alignment of the reference and variant sequences resulted in nucleotide deletion and shortening of the orf8 region. The orf8 regions of Omicron (EPI_ISL_6959868), Kappa (EPI_ISL_3277130), and Eta (EPI_ISL_3278349) variants is identical to the reference, according to the similarity index. With 55.84% sequence similarity to the reference variant alpha (EPI_ISL_3300249) and Iota (EPI_ISL_3298341) showed insertions in two locations and deletions in six. Through nine insertions and six deletions, the beta (EPI_ISL_3299558) variant had a sequence similarity of 55.59%. Furthermore, deletion of nine and insertion of five sites, lambda (EPI_ISL_3298366) and delta (EPI_ISL_3299048) variants revealed 55.06% sequence similarity. Further, three insertions and six deletions in two Brazilian variants, zeta (EPI_ISL_3266119) and gamma (EPI_ISL_3299430), exhibited 51.38% sequence similarity (Fig. 1 ).

Fig. 1

A. Multiple sequence alignment view of orf8 regions using reference and variants. B. Similarity matrix of orf8 regions of reference and variants.

Screening of miRNA targeting the orf8 region

The reference SARS-CoV orf8 region is of 366 nucleotides long mRNA. Total 24 miRNA was predicted by miRDB with target score ranging from 51 to 93%. Among them miRNAs with target score more than 75 was considered for the further analysis. On the basis of the target score nine miRNAs were considered for the reference. Further, orf8 region of the variants were also screened for the potential miRNAs. The miRNA search of Kappa (κ) and Eta (η) variant resulted with 6 miRNAs same as reference. Further we have screened nine unique miRNAs for the variants (β – 2, ι and α – 2, ζ and γ-1, ε – 1, δ and λ-3) (Table 2).

Sequence similarity of orf7a and orf8 region

Multiple sequence alignment result revealed that the orf8 region has sequence similarities with the orf7a region. Delta and lambda variants showed the highest similarity (59.30%) with the conserved orf7a region among the studied variants. Recently emerging variant omicron showed 56.51% sequence similarity. Two Brazilian variants, zeta, and gamma exhibited 58.51% sequence similarity. Further, variant beta, iota and alpha showed 57.19%, 55.14%, and 54.21% sequence similarity respectively. Kappa and eta variants showed 51.58% sequence similarity with the orf7a region.

Molecular docking analysis

Duplex interaction between screened miRNAs and mRNAs of the reference orf8 gene

Given the therapeutic potential of the orf8 gene sequence, human encoded miRNAs that can target the orf8 gene were screened from the miRDB online database. Nine miRNAs were selected for the reference orf8 gene while nine others for the variant strains (Table 2). The duplex formation was done using all the miRNAs with the reference gene.

All the predicted secondary duplex structures which were saved had negative ΔG values which were indicative of a strong binding affinity between the strands at the molecular level (Refer to Table 4 for details). Following the tertiary structure formation, the top 5 predictions were provided, where, the model with the lowest score had the most favorable folded structure. The pdb files of the duplexes were generated and deciphered using the Digital Studio platform. Ligand interactions were studied and it was seen that, in the case of reference miRNAs binding with the reference orf8 gene, favorable binding occurred at the specified seeing region. The variant miRNAs also interacted with the reference orf8 sequence with a good binding affinity.

Table 4

Docking results using miRNAs on reference orf8 region.

miRNA (HSA-mir)	ΔG mRNA-miRNA 2D structure	Docking (Clustering RMSD value 1.5 Å)
		GSCaScore	Transformations
145-5p	−6.00 kcal/mol	23294	2.54	−0.43	−0.82	128.34	264.9	163.61
194-5p	−3.70 kcal/mol	20614	−2.47	0.42	3.1	59.1	23.88	−54.61
219b-3p	−4.50 kcal/mol	18706	−0.01	0.61	−1.38	75.07	73.63	69.9
221-5p	−14.80 kcal/mol	18710	−0.56	−0.61	−2.37	40.21	−67.06	37.98
548 ag	−5.80 kcal/mol	18680	−0.84	0.13	2.2	45.83	18.9	48.71
548ai	−6.50 kcal/mol	18720	0.35	0.05	−1.02	−9.21	156.2	113.64
548ba	−6.50 kcal/mol	22042	−1.42	−0.20	1.49	−31.79	−106.75	55.65
570-5p	−6.50 kcal/mol	18720	0.35	0.05	−1.02	−9.21	156.2	113.64
758-5p	−13.90 kcal/mol	18418	2.81	0.10	2.31	31.52	−337.27	104.27
888-5p	−6.80 kcal/mol	19950	−1.49	0.05	2.75	107.6	−18.13	38.72
1226-3p	−5.60 kcal/mol	20676	−3.02	−0.33	0.4	−67.59	−24.04	−1.47
1537-5p	−5.70 kcal/mol	19244	2.12	−0.34	0.89	−106.19	−66.7	87.67
3059-3p	−4.70 kcal/mol	19558	2.85	−0.04	−2.69	121.02	88.62	75.97
3926	−6.50 kcal/mol	21168	−0.76	−0.41	2.88	118.87	35.65	136.53
4684-3p	−6.40 kcal/mol	21140	1.08	−0.37	−0.67	−42.01	23.7	206
5047	−22.30 kcal/mol	16968	−0.60	0.45	1.91	38.50	10.34	−11.12
7849-3p	−4.80 kcal/mol	20078	−1.33	−0.83	−1.89	−105.51	91.53	403.79
8073	−16.60 kcal/mol	20802	−2.34	−0.84	2.77	58.70	1.91	173.97

GSC-Geometric shape Complementarity.

Study of the binding affinity of the orf8 mRNA-miRNA duplex with Ago2 protein

The tertiary structure of the silencing complex of human AGO2 protein containing a seed guide RNA was retrieved from the Protein Data Bank. Docking analysis on PatchDock gave multiple predictions which were arranged based on a geometric complementarity score. The model with the highest score had the best surface complementarity between the receptor (hAGO2 complex) and the ligand (mRNA-miRNA tertiary duplex). The Geometric Shape Complementarity (GSC) score along with three rotational and three translational transformation values for the top models have been provided in Table 4. The models were viewed and interpreted using Digital Studio (Fig. 2a, Fig. 2b, Fig. 2c a–c ). It was seen that all the 18 screened miRNA-mRNA (orf8) duplexes interacted with the catalytic orf8 protein as well as with the guide template within the protein complex indicating potential regulation of the gene.

Fig. 2a

Docking of orf8 region using miRNAs screened for both (A-I) and variants (J-R).

Fig. 2b

Docking of orf8 region using miRNAs screened for both (A-I) and variants (J-R).

Fig. 2c

Docking of orf8 region using miRNAs screened for both reference (A-I) and variants (J-R).

Duplex interaction between screened miRNAs and mRNA of the reference orf7a gene

Keeping in mind the functional and possible partial sequence homology between orf8 and the conserved orf7a regions of SARS-CoV-2 (Fig. 3 ), human encoded miRNAs targeting the orf8 gene (Table 2) were checked for favorable binding with the orf7a reference sequence (orf7a region of NC_045512.2). Out of 18 miRNAs, nine miRNAs were selected which generally targets the reference orf8 (orf8 region of NC_045512.2) gene while nine others for the variant strains.

Fig. 3

A. Multiple sequence alignment showing similarity between orf7a and orf8 regions. B. Similarity matrix showing percentage sequence similarity of orf7a with orf8 regions of reference and variants.

All the predicted secondary duplex structures which were saved had negative ΔG values. This indicates a strong binding affinity between the strands at the molecular level (Refer to Table 5 for details). Following the tertiary structure formation, the top 5 predictions were provided, where, the most preferable folding was selected based on the computed score. A lower score was indicative of a better tertiary folding. The pdb files of the duplexes were retrieved and using the Digital Studio platform, ligand interactions were studied. It was seen that all the miRNAs formed a successfully folded duplex with the orf7a mRNA strand.

Table 5

Docking results using miRNAs on reference orf7a region.

miRNA (HSA-mir)	ΔG mRNA-miRNA 2D structure	Docking (Clustering RMSD value 1.5 Å)
		GSCaScore	Transformations
145-5p	−2.4 kcal/mol	25820	2.06	0.41	0.52	16.55	−117.69	48.36
194-5p	−5.5 kcal/mol	20512	−0.87	−0.3	−1.64	120.31	92.41	69.67
219b-3p	−1.7 kcal/mol	20758	−1.76	0.2	1.48	9.55	−1.08	131.92
221-5p	−2.6 kcal/mol	19922	−0.2	0.93	2.97	27.29	6.16	47.51
548 ag	−3.3 kcal/mol	18566	2.01	−0.32	−0.95	−67.03	65.01	84.49
548ai	−2.2 kcal/mol	19728	1.25	−0.19	−0.66	−74.07	−43.10	62.96
548ba	−3.2 kcal/mol	19692	−1.36	−0.04	1.67	78.54	−68.70	203.92
570-5p	−2.2 kcal/mol	19728	1.25	−0.19	−0.66	−74.07	−43.10	62.96
758-5p	−3.6 kcal/mol	22634	0.90	−1.26	2.18	55.82	−614.51	1010.94
888-5p	−3.2 kcal/mol	19938	3.05	0.11	−0.35	−35.52	46.11	8.51
1226-3p	−3.4 kcal/mol	18570	−0.60	0.19	1.19	132.65	−180.78	29.6
1537-5p	−4.9 kcal/mol	21232	−0.83	−0.39	−2.24	114.63	19.12	138.83
3059-3p	−4.8 kcal/mol	19524	1.14	−0.91	2.74	78.23	−57.64	66.88
3926	−4.3 kcal/mol	25592	1.20	−0.82	1.85	82.87	−22.11	160.06
4684-3p	−7.1 kcal/mol	18706	0.30	0.74	2.01	33.06	56.52	−13.66
5047	−6.4 kcal/mol	19990	−1.00	−0.61	0.98	−11.44	17.73	141.17
7849-3p	−8.7 kcal/mol	19458	0.66	0.46	−0.52	−88.29	47.31	106.33
8073	−8.6 kcal/mol	20314	0.37	−1.01	−1.07	5.27	29.40	100.05

GSC-Geometric shape Complementarity.

Study of the binding affinity of the orf7a mRNA-miRNA duplex with ago2 protein

Prepared hAGO2 protein and miRNA-mRNA (orf7a) duplexes were uploaded for the final molecular docking. Docking analysis on PatchDock gave multiple predictions where the server algorithm arranged the results based on a geometric complementarity score. The model with the highest score had the best surface complementarity between the receptor (hAGO2 complex) and the ligand (mRNA-miRNA tertiary duplex). The Geometric Shape Complementarity (GSC) score along with three rotational and three translational transformation values for the top models have been provided in Table 5 .

The models were viewed and interpreted using Digital Studio (Fig. 4a, Fig. 4b, Fig. 4c ). It was seen that 15 out of 18 screened miRNA-mRNA (orf7a) duplexes interacted with the catalytic orf8 protein as well as with the guide template within the protein complex indicating potential regulation of the gene. Duplex formed with hsa-mir-221-5p, hsa-mir-548ag, and hsa-mir-888-5p had no interaction with the guide template of the hAGO2 complex. The highest duplex-protein complex interaction was observed in cases of hsa-mir-145-5p and hsa-mir-3059-3p (Fig. 4a, Fig. 4b, Fig. 4cC and K), moderate interactions were seen for hsa-mir-219b-3p, hsa-mir-1537-5p, hsa-mir-548ai, hsa-mir-548ba, hsa-mir-758-5p, hsa-mir-5047 and hsa-mir-8073 (Figure 4O, J, 4H, 4G, 4D, 4E, 4A) while the least interaction was exhibited by duplexes formed by hsa-mir-194-5p, hsa-mir-570-5p, hsa-mir-1226-3p, hsa-mir-3926, hsa-mir-4684-3p and hsa-mir-7849-3p (Figure 4L, F, 4 M, 4Q, 4 N, 4R). There was no interaction with the guide template in the hAGO2 protein strand in duplexes formed by hsa-mir-221-5p, hsa-mir-548ag, and hsa-mir-888-5p (Fig. 4a, Fig. 4b, Fig. 4ci).

Fig. 4a

Docking of orf7a region using miRNAs screened for both (A-I) and variants (J-R).

Fig. 4b

Docking of orf7a region using miRNAs screened for both (A-I) and variants (J-R).

Fig. 4c

Docking of orf7a region using miRNAs screened for both (A-I) and variants (J-R).

Discussion

Viruses and humans appear to be frequently involved in a dominating game, in which a virus infects host cells armed with a variety of methods to elude the host immunity [[76], [77]]. Host cells in return build a defence mechanism to keep the virus at bay [[78], [79]]. One such mechanism is miRNA silencing. The interplay between cellular and virus miRNA is a new and crucial player in the biological world [80]. A virus can generate its miRNAs to facilitate its life cycle in two ways; by targeting viral transcripts to control their lifecycle and targeting cellular transcripts to create a hospitable environment [[81], [82]]. Host cells, on the other hand, can thwart viral replication by releasing a miRNA that can either target viral transcripts directly or indirectly by altering host components implicated in viral infection and pathogenesis [[83], [84]]. Present study is focused on the finding of miRNAs that have the potential to target the viral orf8 along with the orf7a region due to their important roles in the immune system modulation and sequence homology with each other.

The multiple sequence alignment studies done in this analysis once again confirms the deletion of the orf8 region [[26], [27], [28]] in variants from different countries. Several genomic deletions have been observed from the starting of SARS pandemic [85]. During the middle and late phases of the COVID-19 pandemic, all SARS-CoVs had a 29-nt orf8 deletion, but at the end, entire or almost complete orf8 deletions were found [86]. The emergence of these deletions over time and their eventual prevalence led to the idea that orf8 was an evolutionary hotspot for SARS-CoV adaption to humans [86]. The sequence alignment study presented in this analysis showed deletion in orf8 ranging from 35 (β) to 288 (ε) in the selected variants (Table 1; Fig. 1). On the contrary, the Omicron, Kappa and Eta variants share the identical orf8 region in spite of being reported from different geographical locations, which is suggestive of their common ancestry from the wild strain of the virus. Through the multiple sequence alignment, we have also shown that isolates from the same geographical location like the Gamma and Zeta variants exhibit identical orf8 regions. Additionally, high sequence similarity is also observed between the Lambda (Peru) and delta (UK) variants, and the lota (USA) and alpha (UK) variants although they are isolated from different geographical regions (Fig. 1). These findings suggest that the evolution of orf8 region is not guided by zonal or geographical adaptations [87].

We studied the sequence homology of the orf8 region, obtained from 10 different SARS-CoV-2 variants, collected from different parts of the globe (Table 1; Fig. 1). From the multiple sequence alignment data, it was clear that the orf8 region is a hotspot for the virus evolution. However, given the role of this gene in host immune modulation [[27], [34], [35], [42], [63]], we wanted to see the binding efficacy of different miRNAs which targeted the orf8 region. The docking results obtained (Fig. 2a, Fig. 2b, Fig. 2c), showed favorable receptor-ligand interactions (Table 4). All the miRNAs interacted successfully with the AGO2 protein.

Following this docking, we wanted to check for the presence of the orf8 seed sequence in the different variants. For doing so we generated the miRNA-mRNA duplex using three randomly selected reference miRNAs and three variants with the least sequence similarity with the reference orf8 gene (Table 3). Following a thorough analysis of all the variant sequences, it was seen that the seed sequence was not present in any of the variants. However, the same could not be said for the Omicron strain. Since sequence alignment showed the presence of an identical orf8 sequence in the Omicron variant, a seed sequence identical to the reference was detected in its coding strand.

From previous studies, it was seen that the truncations and even silencing of the orf8 region of the SARS-CoV-2 did not affect the replicative fitness of the virus [[26], [27], [28]]. Given this finding, we hypothesized here that, there might be another similar alternative system concealed in the virus which was complementing for the lost orf8 in the deleted variants. While checking for the same, we came across another region of viral genome coding for the open reading frame 7a (orf7a). This non-structural protein of the virus also shows immune-modulating functions and helps the virus in camouflaging inside the host cell. The immune strategies adopted by this protein include but are not limited to induction of cytokine storms and interaction with host integrins via its Ig-like domain [[56], [63], [64], [67]]. What made this gene even more intriguing was its partial sequence homology with the reference orf8 gene and also the fact that the virus seems to be deleting the orf8 while maintaining the orf7a regions of its genome [63]. To further strengthen this hypothesis, we performed a multiple sequence alignment to check exactly how much similarity was being exhibited between the different strains used in this study. From the data obtained we concluded that there was 51.58–59.30% sequence similarity of different variant orf8 regions with the reference sequence of orf7a (Fig. 3). This was a clear indication of partial sequence similarity between the two genes and hence it was possible for the miRNAs screened for orf8 to exhibit a potential binding with the orf7a region.

For our next aim to check the mentioned potential binding, we performed similar docking studies using the reference orf7a sequence and all the 18 miRNAs. The same protocol that was used for orf8 docking was followed in this case. After the final docking on the hAGO-2 protein complex using PatchDock, the top models with the highest GSC scores were selected and analyzed (Fig. 4a, Fig. 4b, Fig. 4c; Table 5). On viewing the ligand interactions on Digital studio, it was seen that 15 out of the list of 18 miRNAs had interacted with the entire catalytic AGO2 complex including the guide template molecule. The three miRNAs where no interaction was seen between the nucleotide duplex and the guide template of hAGO2 were hsa-mir-221-5p, hsa-mir-548ag, and hsa-mir-888-5p (Fig. 4a, Fig. 4b, Fig. 4cB, I, 4P). The best interactions observed showed binding of the nucleotide duplex with up to 4 nucleotides of the template strand and multiple interactions with protein residues of hAGO2. Moderately efficient binding included the involvement of fewer protein residues or, interaction of the duplex molecule binding with up to 3 nucleotides of the guide template. The least efficacy was exhibited by the miRNAs where the nucleotide duplexes interacted with only up to 2 nucleotides of the guide strand present in the hAGO2 complex. Based on this it was concluded that the hsa-mir-145-5p and hsa-mir-3059-3p (Fig. 4a, Fig. 4b, Fig. 4cC and K) had the best binding potential to the orf7a sequence. hsa-mir-145-5p is also associated with other crucial molecules like TP53 which is responsible for cell-cycle arrest and apoptosis [88]. Since the miRNAs used in this study for orf7a molecular docking were the ones screened against the orf8 gene, here we present 15 miRNAs, as potential molecules which could target the orf7a gene as well.

If the orf7a gene acts as the conserved alternative reserve of orf8 in SARS-CoV-2, then targeting this region along with orf8 might severely hamper the viral life cycle. Additionally, if the same set of miRNAs can be used to effectively silence both the regions of the genome, then there is a high chance of improving the overall patient survival as the impaired host immune system will be restored to its full functionality. After further studies based on this computational analysis, if the proposed hypothesis is backed up with more concrete evidence and proven to be true, it might be possible to interfere in the viral replication process of SARS-CoV-2, by targeting the conserved orf7a region in the mutants which exhibit an inactivated orf8 gene (refer to Fig. 5 for hypothesis).

Fig. 5

Possible ways in which a single miRNA can be used to target both orf7a and orf8 region of variants and wild type for the regulation of SARS-COV2.

Limitation

In our study, we have used a bioinformatics approach to demonstrate the possibility of simultaneous targeting of two genes using one miRNA molecule. Our hypothesis was formulated on the basis of the assumption of partial sequence similarity and evolutionary conservation of genes. We have presented a compact analysis however, there are certain limitations as well. A lab-based experimental analysis was not carried out and hence the bioinformatics data has not been validated. Further research will be needed to arrive at conclusive results. Apart from that, there might be other even more conserved regions present in the viral genome which could act as potential targets. Those genes were not looked into and our study primarily focused only on two genes of the open reading frame of the SARS-CoV-2 genome (i.e orf7a and orf8). Additionally, the study could have been carried out using more variants that would have represented the global scenario to a better extent.

Conclusions

Mature miRNAs produced in humans have high therapeutic potential and can be screened against multiple foreign genes to provide immunity against different diseases. The same has also been found to be effective against certain portions of the COVID-19 virus genome where mature human miRNAs have been successfully used to silence viral genes so as to interfere with the life-cycle of the virus. However, given the rapid emergence of new mutant variants, it is crucial to develop a treatment that will be effective against all SARS-CoV-2 variants equally, irrespective of the mutation site. Here in this analysis, we have tried to approach this challenge using a novel strategy. We have performed computational molecular docking analysis to show how human miRNAs screened for a highly variable region of the genome (i.e orf8) can be used to target a similar but comparatively conserved gene of the virus (orf7a). This conserved sequence, not only has partial sequence homology with orf8 but also shares functional similarities with the same. Both these genes are involved in the host immune suppression thereby camouflaging the virus inside the infected cell. Hence, using one molecule for simultaneously silencing both these genes might effectively restore the efficacy of the host immune system. Since the SARS-CoV-2 virus seems to be maintaining its orf7a gene while gradually deleting the orf8 sequence, it might be possible that the emergence of orf8 was a result of gene duplication, with the gene being shortened as a result of evolution. If orf7a acts as the region which has been compensating for the lost orf8 gene function, then there is a possibility that targeting both these sequences might affect the virus cell cycle. With further experimental studies based on this analysis, concrete evidence on the same can be obtained following which a novel treatment for the pandemic can be designed.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.