Potential treatment for chronic myeloid leukemia using microRNA: in silico comparison between plants and human microRNAs in targeting BCR-ABL1 gene

Chronic myeloid leukemia (CML) is a myeloproliferative disorder characterized by the expression of the BCR-ABL1 fusion gene. Tyrosine kinase inhibitors (TKI) are used to treat CML, but mutations in the tyrosine kinase domain contribute to CML chemo-resistance. Therefore, finding alternative molecular-targeted therapy is important for the comprehensive treatment of CML. MicroRNAs (miRNA) are small non-coding regulatory RNAs which suppress the expression of their target genes by binding to the 3′ untranslated region (3′UTR) of the target mRNA. Hypothetically, the miRNA-mRNA interaction would suppress BCR-ABL1 expression and consequently reduce and inhibit CML cell proliferation. Thus, our objective was to determine the target interaction of human and plant miRNAs targeting the 3′UTR region of BCR-ABL1 in terms of miRNA binding conformity, protein interaction network, and pathways using in silico analysis. The 3′UTR sequence of BCR-ABL1 is obtained from Ensembl Genome Browser while the binding conformity was determined using the PsRNATarget Analysis Server, RNA22, Target Rank Server, and DIANA TOOLS. Protein-protein interaction network and pathway analysis are determined using STRING, Cytoscape, and KEGG pathway analysis. Five plants and five human miRNAs show strong binding conformity with 3′UTR of BCR-ABL1. The strongest binding conformity was shown by Oryza sativa’s Osa-miR1858a and osa-miR1858b with −24.4 kcal/mol folding energy and a p value of 0.0077. Meanwhile, in human miRNA, the hsa-miR-891a-3p shows the highest miTG score of 0.99 with −12 kcal/mol folding energy and a p value of 0.037. Apart from ABL1, osa-miR1858a/osa-miR1858b and hsa-miR891a-3p also target other 720 and 645 genes, respectively. The interaction network of Osa-miR1858a/osa-miR1858b and hsa-miR891a-3p identifies nineteen and twelve ABL1’s immediate neighboring proteins, respectively. The pathways analysis focuses on the RAS, MAPK, CML, and hematopoietic cell lineage pathway. Both plant and human miRNAs tested in this study could be a potential therapeutic prospect in CML treatment, but thermodynamically, osa-miR1858a/osa-miR1858b binding to ABL1 is more favorable. However, it is important to carry out more research in vitro and in vivo and clinical studies to assess its efficacy as a targeted therapy for CML.


Background
Chronic myeloid leukemia (CML) is a type of cancer that affects blood and bone marrow cells. It is a clonal disorder of multipotential hematopoietic stem cells that excessively induces proliferation of immature and mature myeloid cells [1]. CML patients are characterized by the accumulation of Philadelphia chromosome-positive (Ph+) myeloid cells. Ph chromosome is the result of a reciprocal translocation between Abelson (ABL) protooncogene on chromosome 9 and the breakpoint cluster area (BCR) on chromosome 22, t(9;22) (q34;q11) [2,3]. This translocation generates the BCR-ABL1 fusion gene, which is the key cause of CML by the development of highly active tyrosine kinase [4]. As the frontline treatment of CML, TKIs are used against the kinase activity of the BCR-ABL protein. However, because TKI resistance in some CML patients are due to tyrosine kinase domain (TKD) mutations, alternative treatments targeting other than TKD are needed [5][6][7]. Hence, miRNA is proposed as an alternative treatment for CML. This is due to the association between miRNA expression and the onset of TKI resistance [4]. MicroRNAs are small non-coding regulatory RNAs which suppress the expression of their target genes. It reduces the stability of messenger RNA (mRNA) and thus inhibits translation [8]. miRNA comprised 20 to 24 nucleotides that are involved in gene regulation and cell signaling. In cancer cells, miRNAs can modulate multiple dysregulated genes or signaling pathways. This is because gene regulation is controlled by the incomplete complementary binding between the seed sequence of miRNA with the 3′ untranslated region (3′UTR) of mRNA. This enables miRNA to simultaneously control the expression of many target genes [9]. Certain miRNAs like hsa-miR-21 and hsa-miR-17-92 clusters are elevated in cancer cells. Meanwhile, other cancer-associated miRNA families such as the hsa-let-7 and hsa-miR-34 with hsa-miR-34 are downregulated [10,11]. These outcomes suggested that miRNAs play a very important role in promoting or preventing cancer. miRNAs can be manipulated as cancer therapeutic agents using miRNA restoration agents known as miRNA mimics. These mimics are used to compensate for the loss of some miRNA [12,13]. Meanwhile, another approach is by using the miRNAinhibiting agents, the antagomiR which functions in blocking the aberrantly overexpressed miRNA. Despite the controversy of cross-kingdom inhibition of cancer cells by plant miRNA, there is a report which indicated that synthetic mimic of plant miRNAs of ath-miR159a (from Arabidopsis thaliana), gma-miR159a-3p, and gma-miR159e-3p (from Glycine max) was able to inhibit the proliferation of breast cancer cells [14]. Meanwhile, another study confirmed the ability of mol-miR168 (from Moringa oleifera) to downregulate the protein expression of SIRT1, when its mimic was transfected into the human hepatoma cell line G2 (HEPG2) cells [15]. Therefore, to understand the RNA silencing mechanism, this study utilized in silico analysis in predicting the miRNA-mRNA interactions. The current study aimed to predict miRNAs that targeted the 3′UTR of the BCR-ABL1 gene and compare their binding affinity, target genes annotation, pathways involved, and proteinprotein interaction in relation to CML pathogenesis.

Predicting miRNA-mRNA target interaction in plants
The PsRNATarget and Schema V2 criteria were used in predicting the interaction. The 3′UTR of the BCR-ABL1 sequence was imported into the software, and miRNAs were selected based on expectation equal to or less than 5. Subsequently, RNA22 (https://cm.jefferson.edu/rna22/ Interactive/) was used to determine miRNA folding energy, p value, and heteroduplex structure by inserting miRNA sequences and single target sequence [19]. Next, miRNA sequences were inserted into Target Rank (http://genes.mit.edu/targetrank/) to obtain the list of target genes [20].

Predicting miRNA-mRNA target interaction in human
Screening of human miRNAs targeting 3′UTR of ABL1 was performed using DIANA TOOLS, with the microT-CDS [21] and miTG scores being set at more than 0.9. Target genes of each selected miRNA were also identified using microT-CDS with the threshold set at 0.7. Predicted miRNAs were later assessed based on their folding energy, p value, and predicted target site using the RNA22 software.

Protein-protein interaction (PPI) network
The interactions between genes could be analyzed quickly by constructing the protein interaction network. STRING which is an online website (http://string-db. org) was used to build the interactive networks of target genes. By restricting the condition to "human species," networks were constructed, and the "TSV" format file of the gene interaction relationship was downloaded, and then the file is imported into the Cytoscape software (http://www.cytoscape.org/) for further analysis and clustering [22]. The Cytoscape MCODE plug-in was employed to perform module analysis of the target network and protein clustering. The module selection criteria were as follows: degree cutoff ≥ 2, node score cutoff ≥ 0.2, K-score ≥ 2, and max depth = 100. The first neighbors of selected nodes used the "directed: outgoing" feature to further analyzed the neighboring protein of ABL1.

Gene Ontology and pathway enrichment analysis
Selected genes from Cytoscape networks were further analyzed in STRING using the Cytoscape plug-in for their pathway enrichment of the nodes. The roles of these genes in the metabolic pathway were elucidated based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway.

Plant miRNA-mRNA target interaction
Five plant miRNAs, namely cpa-miR8154, osa-miR1858a, osa-miR1858b, mtr-miR2654, and sof-miR159e, were predicted to show binding conformation with 3′UTR of BCR-ABL1. The minimum folding energy (MFE) levels needed for the binding of each miRNA to the 3′UTR of BCR-ABL1 were identified. The best MFE score was observed in osa-miR1858a and osa-miR1858b at −24.4 kcal/mol whereas the lowest MFE score was shown by sof-miR159e with −14.4 kcal/mol which is still considered as a strong and energetic binding affinity. The results are summarized in Table 1.

Target gene annotation
Osa-miR1858a and osa-miR1858b were found to have 720 target genes while hsa-miR891a-3p had 645 target genes apart from ABL1 (Table 1). All these genes were imported into STRING analysis for the protein-protein interaction and pathway analysis.

Protein interaction networks, key candidate targets, and pathways analysis
The results showed that the protein network of hsa-miR891a-3p produced 550 nodes and 1550 edges (Fig. 1a). The first neighbors of selected nodes of the ABL1 module were then chosen for further analysis. This module had 12 nodes and 16 edges (Fig. 1b). Pathway enrichment analysis revealed that these nodes were primarily enriched in pathways such as RAS and MAPK signaling pathways. In addition, the encoded proteins were involved mainly in multiple cancer pathways specific to various signaling pathways, including the ErbB signaling pathway (hsa04012), CML (hsa05220), and microRNAs in cancer pathways (hsa05206). The genes that were involved in the pathways were NTRK2, RIN1, NRAS, PPM1B, MAPT, and ABL1. The Gene Ontology of biological processes is emphasized on the cell population proliferation, regulation of apoptotic process, and MAPK cascade ( Table 2). Meanwhile, Osa-miR1858a and osa-miR1858b protein networks showed 588 nodes and 1925 edges (Fig. 2a) from 721 genes. The first neighbors of selected nodes of ABL1 displayed a module of 19 nodes and 39 edges (Fig. 2b). The pathway analysis exhibited an enriched pathway of hematopoietic cell lineage and CML. Furthermore, these nodes were also involved with pathways such as ErbB, RAS (hsa04014), Jak-Stat (hsa04630), and Wnt signaling pathways (hsa04310). The genes that were included in the pathways are CDL, STAT5B, BCR, SHC4, ABL1, RAC2, CCND2, PDGFB, XIAP, PAK2, CD5, and TFRC. Moreover, Gene Ontology of biological process focused on the cell population proliferation, hemopoiesis, regulation of apoptotic process, and regulation of MAPK cascade.

Discussion
CML is a myeloid clonal proliferation characterized by the translocation of the BCR-ABL1 gene in the Philadelphia chromosome. TKI is preferred as the frontline treatment to treat CML. However, in approximately 30% of CML patients, TKI resistance has become a challenging problem [23]. Therefore, finding the molecular targets for CML treatment is of primary importance to solve this problem. Thus, the current research aimed to determine and evaluate miRNAs from plants and humans that targeted the 3′ UTR of the BCR-ABL1 gene in CML in terms of folding energy, gene annotation, and gene enrichment analysis. It is known that miRNA expression is often altered in diseases, especially in cancers. Either acting as an oncogenic or a tumor-suppressive agent, miRNA had gained the researcher's interest to be used as an alternative therapeutic agent to enhance the current treatment for cancer. Furthermore, cross-kingdom studies of miRNAs between plant and animal kingdoms are intensely explored in recent days. Therefore, this research employed the in silico analysis to predict the specific plant and human miRNAs that can target 3′UTR of the human BCR-ABL1 gene in CML. This is because a cross-kingdom interaction is possible due to the presence of various conserved regions in the human and plant genome [14,24]. In addition, due to the highly conserved region in the composition of miR-NAs processing machinery, the biogenesis and working mechanism of miRNAs in plants showed a similar molecular and development roles as in humans [25]. Besides, it was found that sequences of diverse endogenous plant miRNAs were complementary to the components of the genomes and transcriptomes of humans and animals [26]. The present study predicted plant miRNAs by using the psRNATarget. psRNATarget is a platform that was developed to identify plant small RNA (sRNA) targets by  Table 2 The first neighbors of selected nodes (ABL1) of osa-miR1858a, osa-miR1858b, and hsa-miR891a-3p target genes, pathways, and Gene Ontology analyzing complementary matching between sRNA and target mRNA sequences. It employed a predefined grading schema that evaluated the target site accessibility between miRNA and mRNA [27]. The analysis was continued by exporting miRNA sequences into another platform including RNA22 and Target Rank for further validation. These platforms were specifically designed to predict the interaction between miRNA and mRNA in humans. Positive results from RNA22 and Target Rank indicated that the selected plant miRNAs were capable of binding with 3′UTR of BCR-ABL1 in humans. From the in silico analysis, the current study predicted five plants and five human miRNAs that have binding conformity with 3′UTR of the BCR-ABL1 gene. These miR-NAs were selected based on their folding energy, p value, psRNATarget score, and miTG score. For plant miRNAs, osa-mir1858a and osa-mir1858b showed the most promising results. They have the lowest folding energy of −24.4 kcal/mol and a significant p value of 0.007. Thus, the binding between osa-miR1858a and osa-miR1858b and 3′UTR of BCR-ABL1 was much favorable. Nevertheless, in human miRNAs, it was observed that hsa-miR891a-3p showed the highest miTG score of 0.99 with −12kcal/mol folding energy and a p value of 0.037. Therefore, binding of hsa-miR891a-3p towards 3′ UTR of BCR-ABL1 was less favorable compared to osa-miR1858a and osa-miR1858b towards 3′UTR of BCR-ABL1. A miRNA can have multiple targets, and one mRNA can be targeted by many miRNAs [28]. This is observed in Table 1, showing that hsa-miR891a-3p, osa-miR1858a, and osa-miR1858b have numerous target genes despite ABL1. miRNAs can have multiple targets mainly due to the partial complementary between miRNA and its target mRNA which is known as the seed pairing [29,30]. Because of this feature, in silico tools can predict thousands of miRNA-regulated genes. Recently, studies have found that hsa-miR-891a and hsa-miR-891b can be used as prognostic biomarker and therapeutic target in different types of cancer such as lung cancer [31], breast cancer [32], colorectal cancer [33], and gastrointestinal stromal tumor [34]. In concordance with the present study, hsa-miR891a-3p, osa-miR1858a, and osa-miR1858b were also predicted to have another 645 and 720 targets genes, respectively. Thus, it is important to analyze all the target genes of hsa-miR891a-3p, osa-miR1858a, and osa-miR1858b. Further analysis on the protein-protein interaction network and KEGG pathway showed that hsa-miR891a-3p, osa-miR1858a, and osa-miR1858b were all involved in the MAPK signaling pathway. The MAPK signaling pathway is very critical for cell existence, and its signaling cascade plays a central role in the pathogenesis of CML [35]. In addition, fundamental cellular processes such as cell proliferation and cell survival are regulated by this pathway [36,37]. Moreover, the RAS/MAPK signaling can also lead to the survival of positive BCR-ABL cells under imatinib pressure [38]. Previous findings suggested that apoptosis could be triggered by the inhibition of MAPK signaling in the K562 cell line, and the use of MAPK-specific inhibitors showed cooperative effects with imatinib in CML cells [39]. Furthermore, among the targeted genes of osa-miR1858a and osa-miR1858b, the STAT5 gene was one of the important targets. This is because STAT5 promoted imatinib resistance by significantly being involved in the activation of hematopoietic cells transformation by p210Bcr/Abl [40]. The activation of this gene would speed up cell division and inhibit DNA repair, which thus causes genomic instability such as the blast crisis of CML [40,41]. In addition, the current Gene Ontology analysis also showed that STAT5, CCND2, and PDGFB were involved in the JAK/STAT pathway. This finding is supported by a previous study which indicated that STAT5 could regulate the transcription factor and modulate cytokine receptors and growth factors [42]. Therefore, this pathway is crucial in mediating cytokine signaling by promoting the growth of leukemic cells [40]. According to the KEGG pathway enrichment analysis of osa-miR1858a and osa-miR1858b, the pathways are concentrated on the hematopoietic cell lineage and CML pathway. miRNAs are important regulators involved in many cell functions including the hematopoietic cell lineages. It has been reported that the miRNA expression is regulated during hematopoietic differentiation, suggesting miRNA involvement in the process [43]. Thirteen miRNAs with differential expression in peripheral blood cells have been documented in a previous study, and these hematopoieticspecific miRNAs are essential for blood cells to function properly. In addition, another study also observed that the differently expressed miRNAs allowed the cell-specific transgene expression in human hematopoietic cell lineages [44]. Apart from that, osa-miR1858a and osa-miR1858b additionally involved in the ErbB, Jak-Stat, and Wnt pathways. The Wnt pathway plays an important function during the growth of normal hematopoietic stem cells (HSC) and leukemic stem cells (LSC). Activation of the Wnt/beta-catenin pathway plays a pivotal part in maintaining the self-renewal potential in HSCs. In CML, the BCR-ABL1 oncogene stimulates the β-catenin within the blast stage thus encouraging the LSC reestablishment and repressing imatinib's activity [45]. In addition, another study demonstrated that the Wnt-3a pathway could promote the proliferation and survival of acute lymphoid leukemia cells [46]. Nevertheless, target gene annotation results observed that osa-miR1858a and osa-miR1858b targeted both the BCR and ABL1 genes which perhaps might help in further reduction of CML cell proliferation. Therefore, the determination of miRNA's target pathways and genes that are involved in imatinib resistance could be a crucial data that might perhaps be manipulated as alternative treatments for CML therapy.

Conclusion
In summary, Osa-miR1858a/Osa-miR1858b showed the strongest binding conformity to 3′UTR of BCR-ABL1 with −24.4 kcal/mol folding energy. Nevertheless, hsa-miR891a-3p was chosen based on an excellent miTG score of 0.99 with −12 kcal/mol folding energy. The interaction network of Osa-miR1858a/osa-miR1858b and hsa-miR891a-3p identified 19 and 12 ABL1's immediate neighboring protein, respectively. These proteins were involved in the RAS, MAPK, CML, and hematopoietic cell lineage pathway. miRNA application as a biomarker is being widely used. Therefore, the determination of plants and human miRNA determination using in silico analysis that could be used as the molecular targeted therapy for CML was our main objective in the current study. By using in silico analysis, this study successfully predicted miRNAs from plants and humans that were able to bind with 3′UTR of BCR-ABL1.
Theoretically, these miRNAs are capable of suppressing BCR-ABL1 expression, thus reducing CML cell proliferation. Therefore, it is predicted that these miRNAs could inhibit the BCR-ABL1 by binding to the 3′UTR and would be the best potential molecular targeted therapy for CML. In silico study is an informative and powerful tool in helping the researcher to form an estimation of miRNA-mRNA target interaction. For future study, in vitro work could be carried out using CML cells to validate the Osa-miR1858a/ osa-miR1858b and hsa-miR891a-3p capability in inhibiting CML cell proliferation by targeting 3′UTR of the BCR-ABL1 gene.