The majority of hepatocellular carcinoma cases are caused by infection with hepatitis B (HBV) or C (HCV) viruses. CTNNB1 is the most mutated oncogene in HBV- and HCV-associated tumors. CTNNB1 mutations can lead to β-catenin accumulation, resulting in tumor progression. Small interfering RNAs (siRNAs) can be used to silence CTNNB1 mRNA. After prediction and evaluation, four siRNAs were found to have the highest silencing potential. All four siRNAs had an acceptable GC content, no palindromic sequences, no off-targets, were thermostable, and had accessible target sites. Molecular docking of the siRNAs to Argonaute 2 demonstrated favorable docking scores within the binding pocket for three siRNAs. Molecular dynamics simulations and binding energy calculations demonstrated that the siRNAs steadily remained in the binding pocket. In this study, three siRNAs were successfully designed to silence oncogenic CTNNB1 mRNA as a therapeutic strategy against hepatocellular carcinoma and warrant further in vitro and in vivo validation.

Argonaute 2 (Ago2), molecular docking, molecular dynamics, small interfering RNA (siRNA), RNA interference (RNAi)

Hepatocellular carcinoma (HCC) is the most common type of primary liver cancer and accounts for approximately 90% of all cases (Llovet et al. 2021). HCC progression involves a series of genetic, epigenetic, and molecular alterations in hepatocytes that are mainly associated with chronic inflammation in the liver (Panneerselvam et al. 2023). Risk factors linked to the development of HCC include chronic hepatitis B virus (HBV) and hepatitis C virus (HCV) infections, which represent the most predominant etiology (Farazi and DePinho 2006; Akinyemiju et al. 2017). In particular, HCV accounts for 29.5% of HCC cases, compared to 8.2% of HCC cases caused by HBV (Shiels et al. 2019).

Somatic mutations in various oncogenes and tumor suppressor genes are considered drivers of HCC development. Some of these mutations may be caused by the integration of HBV DNA into the host genome (Yeh et al. 2023) or by dysregulation of important cellular processes by HCV-associated proteins (Irshad et al. 2017). CTNNB1 is the most frequently mutated oncogene in HCC, with mutations occurring in more than 10% of tumors. CTNNB1 encodes β-catenin, an important intracellular signal transducer in the Wnt signaling pathway. In the absence of Wnt, β-catenin is phosphorylated at the N-terminal domain and subsequently degraded. However, mutations in serine and threonine residues in exon 3 of CTNNB1 impair the phosphorylation site of β-catenin, preventing its degradation. This results in the aberrant activation of the Wnt pathway, which drives HCC progression (Tornesello et al. 2013) (Fig. 1).

Figure 1. 

Graphical representation of how HBV/HCV-induced CTNNB1 mutations drive HCC progression. HBV/HCV infection can introduce somatic mutations in exon 3 of CTNNB1, impairing the phosphorylation site and subsequently preventing the degradation of β-catenin, allowing aberrant activation of transcription factors, and driving HCC progression. hepatitis B virus; HCV, hepatitis C virus; TCF/LEF, T-cell factor/lymphoid enhancer factor; HCC, hepatocellular carcinoma.

RNA interference (RNAi) is a naturally occurring biological process in which short, double-stranded RNAs, such as microRNAs (miRNAs) and small interfering RNAs (siRNAs), recognize and degrade mRNA via post-transcriptional gene silencing. In the siRNA RNAi pathway, double-stranded RNA (dsRNA) is recognized by Dicer and cleaved into a siRNA duplex (Katoh and Suzuki 2007). The siRNA duplex is incorporated into the endonuclease Argonaute 2 (Ago2), a component of the RNA-induced silencing complex (RISC), and is unwound by an ATP-dependent helicase. One siRNA strand, known as the passenger strand, is discarded and degraded, while the other strand, known as the guide strand, recognizes and binds to a complementary mRNA sequence, allowing it to be cleaved by RISC and degraded (Leung and Whittaker 2005) (Fig. 2).

Figure 2. 

Graphical representation of the RNAi pathway for siRNAs. The siRNA guide strand is loaded onto Ago2, while the passenger strand is discarded. Ago2 associates with other proteins to form RISC and facilitates binding of the guide strand to a complementary mRNA sequence. The endonuclease activity of Ago2 cleaves the mRNA, thereby silencing its expression. Ago2, Argonaute 2.

The use of computational tools to predict and evaluate potential siRNA candidates allows for the efficient screening of hundreds to thousands of short RNA sequences at once. The selection of siRNAs with the highest silencing potential was successfully validated experimentally in previous studies. As such, siRNAs that were computationally designed against the MERS-CoV-orf1ab gene using siRNA prediction algorithms and evaluation of GC content, off-targets, target site accessibility, thermodynamic stability, and secondary structure formation were validated and were found to inhibit viral replication in Vero cells (Sohrab et al. 2021). Another study by Sohrab et al. (2022) designed siRNAs against SARS-CoV-2-S-RBD using similar prediction and evaluation methods and molecular docking of the predicted siRNAs to their target sequences. The predicted siRNAs demonstrated antiviral efficiency without cytotoxicity in Vero E6 cells. Similar methods were used to predict and evaluate siRNAs against a highly conserved region of the MERS-CoV membrane (M), with the top candidates demonstrating antiviral efficacy without cytotoxicity in Vero cells (El-Sayed et al. 2023).

The regulation of β-catenin expression through siRNA-mediated silencing of CTNNB1 may serve as a potential therapeutic strategy against HCC. To rapidly screen and evaluate all possible sequences, this study aimed to design siRNAs targeting the coding sequence (CDS) of oncogenic CTNNB1 using computational methods. In particular, this study aimed to predict and evaluate potential siRNAs targeting CTNNB1 using known algorithms and model the behavior of the top siRNA candidates through molecular docking with Ago2 and molecular dynamics simulation.

The use of computational tools to efficiently predict and evaluate potential siRNAs has been previously employed in targeting various genes, whether within the host or the pathogen, for the treatment of several diseases. This study offers a novel, more stringent, and broadly targeting screening method for predicting and evaluating siRNAs that are highly likely to have high silencing efficacy in experimental setups. The consideration of several transcript variants of CTNNB1 instead of just the most common ones guarantees that the siRNAs will be able to treat a larger cohort of patients. Furthermore, extensive molecular docking and molecular dynamics analyses of the siRNA candidates with Ago2 and comparison to an experimentally derived control offer more extensive insights than prior siRNA design studies.

To facilitate knockdown of β-catenin, siRNAs were designed around the coding region of CTNNB1. The CDS of CTNNB1 was chosen as the region to design siRNAs as it is retained in the mature mRNA and is the region that is translated into β-catenin (Tornesello et al. 2013). Although the 3’UTR and 5’UTR are also retained after the maturation process and could be used for siRNA design, these regions were not conserved across the retrieved CTNNB1 transcript variants and could not have been used to design broadly active siRNAs.

Different siRNA prediction servers implement different algorithms that consider various parameters that affect the siRNA efficacy. A total of 45 siRNAs were initially predicted using siDirect 2.1. siDirect 2.1 uses the Ui-Tei, Reynolds, and Amarzguioui algorithms (Amarzguioui and Prydz 2004; Reynolds et al. 2004; Ui-Tei 2004). These algorithms are based on observations made regarding siRNA efficiency in the presence or absence of certain nucleotides at certain positions on both the guide and passenger strands (Table 1). While these algorithms are widely used in siRNA prediction studies, they do not consider the interaction of siRNAs with other elements in the RNAi pathway (Yasmin et al. 2022). Therefore, in this study, a siRNA candidate meeting the conditions of at least one algorithm was considered sufficient, as other siRNA prediction algorithms were employed. Out of the 45 initially predicted siRNAs, 27 were cross-validated by OligoWalk, siPRED, and siExplorer. OligoWalk and siPRED predict efficient siRNA candidates based on the mechanics of the interaction between the siRNA and the target mRNA (Lu and Mathews 2008; Pan et al. 2011), while the siExplorer algorithm predicts siRNAs based on the finding that specific nucleotides at every third position (3n + 1) greatly enhance their binding affinity to transactivation response element RNA-binding protein (TRBP), a partner protein of Dicer (Katoh and Suzuki 2007). Using this highly selective method of siRNA prediction, candidate siRNAs are guaranteed to generate favorable results in downstream evaluations.

The 27 predicted siRNAs underwent evaluation for GC content, presence of palindromic sequences, presence of off-targets, secondary structure formation, target site accessibility, and thermostability. The GC content and presence of palindromic sequences were checked using OligoCalc. A high GC content has been suggested to inhibit the unwinding of the siRNA duplex, preventing RISC loading (Liu et al. 2012). A GC range of 31.6–57.9%, is considered the optimal range for siRNAs, as suggested by Amarzguioui and Prydz (2004). Palindromic sequences can prohibit the proper formation of the mRNA-siRNA duplex and must be avoided (Fakhr et al. 2016). The presence of off-targets was checked using BLAST and cross-validated using GESS. siRNAs and other small RNAs can tolerate mismatches at the target site and may induce silencing even with imperfect complementarity (Yasmin et al. 2022).

Target site accessibility was assessed using Sfold. Sfold returns various scores that assess the functionality of the siRNA throughout the entire RNAi pathway, from loading onto RISC to binding and cleavage of the target mRNA. A siRNA that meets the thresholds for all scores is highly likely to be functional. The target accessibility score is based on the antisense (guide) siRNA binding energy, which is the weighted sum of the RNA/RNA stacking energies of the mRNA-siRNA duplex. The more negative the antisense siRNA binding energy, the higher the score. According to the target accessibility rule, the antisense siRNA binding energy must be less than or equal to -10 kcal/mol. The AIS values of the candidate siRNAs indicate that the mRNA-siRNA duplexes were predicted to have favorable stacking interactions.

To silence the intended target site, the correct siRNA strand must first be loaded onto RISC. The DSSE is the difference in stability between the 5’ end of the guide strand and the 5’ end of the passenger strand. The DSSE must be greater than 0 kcal/mol (i.e., the 5’ end of the guide strand is less stable than the 5’ end of the passenger strand), as RISC forms a complex with the strand with the less stable 5’ end (Elhefnawi and Mysar 2011). The DSSE values of the candidate siRNAs confirm that the desired guide strands will be preferably loaded into RISC. The duplex feature score is based on the Reynolds algorithm (Reynolds et al. 2004) and adds or deducts points based on whether certain criteria are met, with the highest possible score being 10 and the lowest possible score being -2. siRNAs with scores greater than 6 are considered acceptable.

After ensuring that the guide strand is loaded onto RISC, its ability to bind to the target mRNA so that Ago2 can be correctly positioned to perform cleavage activity is assessed. The AIS is the average of internal stability values for positions 9–14 of the guide strand and should be greater than or equal to -8.6 kcal/mol.

For the remaining four (4) siRNA candidates, the secondary structures of the siRNA guide strands, and the mRNA-siRNA duplexes were predicted. A positive free energy of folding for siRNA guide strands indicates higher binding efficiency as the molecule is less prone to secondary structure formation (Chowdhury et al. 2021). In contrast, a negative binding energy for mRNA-siRNA duplexes indicates a stronger interaction between the mRNA and siRNA (Yasmin et al. 2022). Tm (Conc) is the temperature at which the concentration of the double-stranded mRNA-siRNA duplex is half its maximum value, while Tm (Cp) is the temperature at which the collective heat capacity reaches its maximum value. Higher Tm values indicate the effectiveness of the siRNA, with a higher melting temperature denoting a more effective siRNA (Yasmin et al. 2022).

The tertiary structures of the siRNAs and Ago2 were modeled in preparation for molecular docking. The siRNAs were modeled by mutating the bases of the template miRNA so as to retain the conformation of the small RNA when it is within the binding pocket of Ago2. The tertiary structure of Ago2 was predicted and refined. GalaxyRefine returns various scores that describe the refined models with respect to the initial model. The global distance test-high accuracy (GDT-HA) and root mean square deviation (RMSD) are measures of similarity between the positions of the alpha carbons of two protein structures with the same amino acid sequence. GDT-HA calculates the proportion of alpha carbons located within a specified proximity to one another (Zemla 2003), while RMSD measures the displacement of alpha carbons from their initial positions (Sargsyan et al. 2017). The MolProbity score is a combination of the clash score, poor rotamers score, and Rama-favored scores. The clash score represents the number of unusually close atom pairs in a protein model that may result in steric clashes; the poor rotamers score is the number of amino acids with poor conformations; and the Rama-favored score is the percentage of residues with favorable backbone ϕ and ψ angles (Yasmin et al. 2022). A lower MolProbity score indicates a higher-quality refined structure; therefore, Model 1 was considered the best-refined structure of Ago2.

ProSA-web assesses the overall and local model qualities of the protein structure. The overall model quality plots the z-score of the query structure against that of the experimentally derived protein structures arranged by size. A good-quality structure is expected to have a z-score similar to that of protein structures of similar sizes. The local model quality plots the knowledge-based energy of the protein structure as a function of amino acid sequence position, where negative values indicate good quality regions (Wiederstein and Sippl 2007). ERRAT compares the distribution of non-bonded atomic interactions in the query structure to distributions in a database of known, high-resolution structures; a structure with a score of 95% or greater is considered to be high-resolution (Yasmin et al. 2022). PROCHECK plots each residue in the query structure in the Ramachandran plot based on their backbone ϕ and ψ angles and categorizes them based on their favorability. Protein structures with more than 90% of residues in the most favored regions are considered to be of good quality (Sumitha et al. 2020).

Out of the four (4) candidate siRNAs, three (3) docked to Ago2 in the expected orientation. These three siRNAs had common interacting residues within the MID and PAZ domains. Y529 in particular is highly conserved and serves as an insertion site for the phosphorylated 5’-end of the guide strand. This interaction is critical for the binding and cleavage activity of the guide strand. Meanwhile, the non-specific anchoring of the hydroxylated 3’ overhang of the guide strand to the PAZ domain is crucial for preventing its degradation (Wu et al. 2020).

The docked complexes of the four siRNAs underwent 100-ns molecular dynamics simulations. The results of the MD simulations showed that s3 and s22 behaved in a manner more similar to the reference miRNA, while s19 and s28 displayed more erratic fluctuations. The relatively higher RMSD values of the s28-Ago2 complex were most likely caused by the random movements of the nondescript loops of Ago2 interacting with the mobile water molecules. Meanwhile, the high RMSD and RMSF values of the s19-Ago2 complex are likely attributed to the incorrect binding pose of s19, where the 5’ end docked against the PAZ domain and the 3’ end docked against the MID domain. The interaction between the mismatched Ago2 domains and the ends of the siRNA might have generated repulsive forces that weakened the binding affinity and stability of the Ago2-s19 complex compared to the other siRNA complexes. Furthermore, the constant movement of s19 may have hindered its ability to establish a consistent number of hydrogen bond interactions with the protein. Despite its stronger affinity compared to the reference miRNA, as shown by the binding energy calculations, the incorrect orientation of s19 against Ago2 renders it unsuitable as a viable siRNA to silence CTNNB1.

It is interesting to note that the majority of the residues with the strongest binding energy contribution are arginine, a positively charged amino acid. Arginine residues in the PAZ, MID, and PIWI domains are known to significantly interact with siRNAs (Bhandare and Ramaswamy 2016a). These residues may have contributed to enhancing the electrostatic forces from the charge-dipole interaction against the polar siRNA molecules, thereby leading to stable and strong complex binding. An additional arginine R761 interacts with the three correct siRNA molecules and s0, as well as Y529, which, as mentioned earlier, is a highly conserved insertion site for the guide strand. Aside from the flipped conformation, the complete absence of an interaction between s19 and Y529 likely contributed to the stark difference in results.

These results show that in silico methods can be employed for the rational design of siRNAs targeting oncogenic CTNNB1. The computational design of therapeutic siRNAs has previously been employed to target wild-type or mutant human genes to combat various diseases, such as mutant tardbp against amyotrophic lateral sclerosis (Bhandare and Ramaswamy 2016b), CLEC5A against dengue and Japanese encephalitis virus (Yasmin et al. 2022), and oncogenic KRAS against cancer (Ramalingam and Arumugam 2023). Knockdown of CTNNB1 using siRNAs has previously been shown to decrease β-catenin activity in hepatoma cell lines (Zeng et al. 2007), although clinical trials have yet to be conducted. In silico screening allows us to rapidly characterize thousands of nucleotide sequences from a single gene and identify those with the highest silencing potential. Likewise, molecular docking and molecular dynamics studies with Ago2 provide a visualization of how siRNAs insert and interact with the binding pocket, as well as how this interaction changes throughout a minuscule timeframe under simulated conditions in the human body. The use of in silico methods in developing therapeutics is a rapidly growing field, and as technology evolves, computational results become more accurate, precise, and informative and greatly supplement experimental validation.

In this study, siRNAs were predicted from the CDS of CTNNB1 and evaluated for their silencing capability. The top four siRNA candidates underwent molecular docking and molecular dynamics simulation. s3, s22, and s28 demonstrated stability similar to the reference miRNA (s0), while s19 exhibited relatively weaker stability due to its incorrect binding orientation. Further analysis using MMPBSA and MMGBSA revealed that all four siRNA molecules exhibited stronger binding affinity towards the Ago2 protein compared to the miRNA, underscoring their potential for silencing CTNNB1. Although s19 showed promising binding energy scores, its mismatched binding pose indicates its unsuitability as a siRNA molecule, prompting the need for additional investigation. Additionally, we identified five arginine residues (R68, R179, R710, R792, and R795) with the highest binding energy contribution to siRNA binding. These specific residues hold potential as crucial targets for the optimization of siRNA design. Considering all evaluations, three (3) siRNAs were successfully designed to silence oncogenic CTNNB1, as demonstrated by their favorable properties and stable complex formation with Ago2. Although the results warrant further in vitro and in vivo validation, the current findings provide valuable insights and a strong foundation for further research.

The authors would like to thank the Philippine Council for Agriculture, Aquatic, and Natural Resources Research and Development (PCAARRD), which funded the project under the Virology and Vaccine Research and Development Program. The authors would also like to thank the Department of Science and Technology—S&T Fellows Program, the Department of Science and Technology—Science Education Institute (DOST-SEI) Career Incentive Program, the Department of Science and Technology—Industrial Technology Development Institute (DOST-ITDI), and the Department of Science and Technology—Advanced Science and Technology Institute (DOST-ASTI) Computing and Archiving Research Environment (COARE) for their support of this research.