Analyzing the Interactions of mRNAs and ncRNAs to Predict Competing Endogenous RNA Networks in Osteosarcoma Chemo-Resistance

Please cite this article in press as: Zhu et al., Analyzing the Interactions of mRNAs and ncRNAs to Predict Competing Endogenous RNA Networks in Osteosarcoma Chemo-Resistance, Molecular Therapy (2019), https://doi.org/10.1016/j.ymthe.2019.01.001

Original Article

Analyzing the Interactions of mRNAs and ncRNAs to Predict Competing Endogenous RNA Networks in Osteosarcoma Chemo-Resistance Kun-Peng Zhu,1,2,3,4 Chun-Lin Zhang,1,2,3,4 Xiao-Long Ma,1,2 Jian-Ping Hu,1,2 Tao Cai,1,2 and Lei Zhang1,2 1Department

of Orthopedics, Shanghai Tenth People’s Hospital, Tongji University, School of Medicine, Shanghai 200072, PR China; 2Institute of Bone Tumor, Tongji

University, School of Medicine, Shanghai 200072, PR China

Chemo-resistance is a huge obstacle encountered in the osteosarcoma (OS) treatment. Protein-coding mRNAs, as well as non-coding RNAs (ncRNAs), including long ncRNA (lncRNA), circular RNA (circRNA), and microRNA (miRNA), have been demonstrated to play an essential role in the regulation of cancer biology. However, the comprehensive expression profile and competing endogenous RNA (ceRNA) regulatory network between mRNAs and ncRNAs in the OS chemo-resistance still remain unclear. In the current study, we developed whole-transcriptome sequencing (RNA sequencing [RNAseq]) in the three paired multi-drug chemo-resistant and chemo-sensitive OS cell lines to comprehensively identify differentially expressed lncRNAs, circRNAs, miRNAs, and mRNAs. Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were performed for mRNAs with significantly different expression. Then the ceRNA networks combining lncRNAs, circRNAs, miRNAs, and mRNAs were predicted and constructed on the basis of the authoritative miRanda and TargetScan databases combined with the widely accepted vital drug resistance-related genes and signal transduction pathways. In addition, two constructed ceRNA regulatory pathways, lncRNAMEG3/hsa-miR-200b-3p/ AKT2 and hsa_circ_0001258/hsa-miR-744-3p/GSTM2, were randomly selected and validated by real-time qPCR, RNA immunoprecipitation (RIP), RNA pull-down assay, and dual luciferase reporter gene system. Taken together, our findings may provide new evidence for the underlying mechanism of OS chemo-resistance and uncover some novel targets for reversing it.

chemo-resistance in OS needed to be further illuminated to provide the clinician with some novel treatment targets and tools.1,4 With the rapid development of high-throughput sequencing technology, many non-coding RNAs (ncRNAs), previously regarded as junk molecules, have been identified and demonstrated to play a vital role in the cellular multiple physiological and pathological processes.5,6 Actually, less than 2% of the total genome encodes protein-coding genes, suggesting that ncRNAs represent most of the human transcriptome.7 On basis of the size and shapes, these ncRNAs can be subdivided into small ncRNAs (<200 nt long), such as microRNAs (miRNAs), and long ncRNAs (lncRNAs), with a length >200 bp, as well as circular RNAs (circRNAs) with a covalently closed circular structure.8 Increasing evidence has shown that these ncRNAs could participate in the regulation of tumor progression, including cell proliferation, migration, metastasis, and chemo-resistance, through regulating gene expression at the levels of pre-transcription, transcription, and post-transcription.9–11 Actually, some ncRNAs have been reported to be involved in the regulation of OS multi-drug chemo-resistance. For example, Pu et al.12 reported that miR-34a-5p promotes the multi-drug chemo-resistance of OS via repression of the AGTR1 gene. Xu et al.13 showed that miR-34c suppresses OS metastasis and chemoresistance by targeting Notch1 and LEF1. Besides, Wang et al.14 found that lncRNA CTA sensitizes OS cells to doxorubicin (DXR) through inhibition of autophagy. Zhu et al.15 showed that lncRNA FENDRR sensitizes DXR resistance of OS cells through downregulating ABCB1 and ABCC1. However, there is seldom a report about circRNA in the OS chemo-resistance. Recently, competing endogenous RNA (ceRNA) networks have been reported as an important mechanism to explain a post-transcriptional

INTRODUCTION Osteosarcoma (OS) is the most common type of primary malignant bone tumor among children and adolescents, which is characterized by high incidence in the metaphysis, such as the distal femur and proximal tibia.1 During the few decades, the 5-year survival rate has increased to 70%–80% by the wide resection surgery combined with multi-drug adjuvant chemotherapy. However, the 5-year survival rate of patients who experience chemo-resistance markedly decreased to less than 20%, owing to the poor sensitivity to secondline chemotherapy.2,3 The mechanism underlying the multi-drug

Received 13 September 2018; accepted 2 January 2019; https://doi.org/10.1016/j.ymthe.2019.01.001. 3

Senior author

4

These authors contributed equally to this work.

Correspondence: Chun-Lin Zhang, MD, Department of Orthopedics, Shanghai Tenth People’s Hospital, Tongji University, School of Medicine, 301, Yan-chang Middle Road, Shanghai 200072, PR China. E-mail: [email protected]

Molecular Therapy Vol. 27 No 3 March 2019 ª 2019 The American Society of Gene and Cell Therapy.

1

Please cite this article in press as: Zhu et al., Analyzing the Interactions of mRNAs and ncRNAs to Predict Competing Endogenous RNA Networks in Osteosarcoma Chemo-Resistance, Molecular Therapy (2019), https://doi.org/10.1016/j.ymthe.2019.01.001

Molecular Therapy

layer of gene translation regulation.16 These ceRNAs include various types of RNAs, such as lncRNAs, circRNAs, mRNAs, and pseudogenes, which could compete for the same miRNA response elements (MREs) to mutually regulate.17 The potential regulatory pathways including lncRNA-miRNA-mRNA and circRNA-miRNA-mRNA could be constructed on the basis of their shared bridge miRNAs. Actually, emerging evidence has revealed that ceRNA networks play a crucial role in tumor development, including chemo-resistance.18,19 For example, An et al.20 reported that lncRNA NEAT1/miR-194/ ZEB1 contributes to paclitaxel resistance of ovarian cancer cells. Gao et al.21 found that lncRNA CASC2/miR-183/Sprouty2 modulates the sensitivity of prostate cancer cells to docetaxel. Besides, Gao et al.22 discovered a regulatory role of the circ_0006528-miR-7-5p-Raf1 axis in DXR-resistant breast cancer. However, the comprehensive expression profile of lncRNA, circRNA, miRNA, and mRNA, as well as the underlying lncRNA or circRNA-miRNA-mRNA ceRNA regulatory networks in the OS chemo-resistance, has not been reported to date. In the current study, we developed whole-transcriptome sequencing (RNA sequencing [RNA-seq]) with Illumina HiSeq2000 in the three paired multi-drug chemo-resistant and chemo-sensitive OS cell lines (MG63/DXR versus MG63, KH-OS/DXR versus KH-OS, U2-OS/ DXR versus U2-OS) to comprehensively identify differentially expressed lncRNAs, circRNAs, miRNAs, and mRNAs. Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were performed for mRNAs regulated by the lncRNA-miRNA or circRNA-miRNA interactions with significantly different expression in OS chemo-resistance. Then the ceRNA networks of lncRNAs, circRNAs, miRNAs, and mRNAs were predicted and constructed on the basis of the authoritative miRanda and TargetScan databases combined with the widely accepted vital drug resistance-related genes and signal transduction pathways. Finally, two constructed ceRNA regulatory pathways were selected and validated by real-time qPCR, RNA immunoprecipitation (RIP), RNA pull-down assay, and dual luciferase reporter gene system. Our findings may provide new evidence for underlying mechanisms of ncRNAs and related ceRNAs networks in OS multi-drug chemoresistance and uncover some novel targets for reversing it.

RESULTS Model Identification of Multi-Drug-Resistant OS Cell Lines

Chemo-resistance was identified by comparing the half maximal inhibitory concentration (IC50) value of chemo-resistant OS cell lines with that of the chemo-sensitive cell lines. IC50 values of paired MG63/DXR (KH-OS/DXR or U2-OS/DXR) and MG63 (KH-OS or U2-OS) when exposed to DXR, cisplatin (DDP), or methotrexate (MTX) for 48 h were shown in Table S1. Obviously, all the three chemo-resistant OS cell lines were more resistant to the three commonly used chemotherapy drugs than the paired paternal chemo-sensitive cell lines, laying a solid foundation for further study. Expression Profiles of lncRNAs, circRNAs, miRNAs, and mRNAs

The whole-transcriptome sequencing data (lncRNA, circRNA, miRNA, mRNA) were acquired from the two groups by the Illumina

2

Molecular Therapy Vol. 27 No 3 March 2019

Hiseq2500 platform. We analyzed differently expressed (DE) ncRNAs (lncRNAs, circRNAs, miRNAs) and mRNAs with the cutoff fold changes R2 or %0.5 along with p < 0.05 and false discovery rate (FDR) <0.05 in the three paired chemo-resistant and chemo-sensitive OS cell lines. The analysis strategy and procedure of the current study is shown in Figure 1. DE ncRNAs and mRNAs in the samples were shown using Volcano plot (Figure 2A) and heatmap (Figure 2B). The results showed that there were 339 DE lncRNAs (150 upregulation and 189 downregulation), 80 DE circRNAs (57 upregulation and 23 downregulation), 21 DE miRNAs (9 upregulation and 12 downregulation), and 2,606 mRNAs (1,607 upregulation and 999 downregulation) in the chemo-resistant group compared with the chemo-sensitive group, respectively. As was displayed in Table 1, the most upregulated lncRNA was lnc-DNMT1-2:1 with 85-fold change, and the most downregulated lncRNA was lnc-RPLP0-4:1 with 55-fold change. Besides, the most upregulated circRNA, miRNA, and mRNA were hsa_circ_0004674 (17-fold change), hsa-miR-337-3p (12-fold change), and C6orf106 (124-fold change), respectively. The most downregulated circRNA, miRNA, and mRNA were hsa_circ_ 0068458 (12-fold change), hsa-miR-203a-3p (10-fold change), and AKR1C2 (197-fold change), respectively. The top 10 upregulated and downregulated lncRNAs, circRNAs, and mRNAs were illustrated in Tables S4–S6. The top 5 upregulated and downregulated miRNAs were listed in Table S7. GO and KEGG Pathway Analysis in These DE ncRNAs and mRNAs

According to the analysis strategy and procedure previously described, we constructed the lncRNA-miRNA-mRNA networks, including lncRNA, miRNA, and mRNA, and the circRNA-miRNAmRNA networks, including circRNA, miRNA, and mRNA. Considering there were so many relationship pairs, we performed the GO and KEGG pathway analysis contrary to the mRNAs involved in the two complex interaction networks, respectively. As was shown in Figure 3A, GO analysis showed that these DE mRNAs in the lncRNA-miRNA-mRNA networks are mainly related to intramembranous ossification (biological_process), collagen type XIV trimer (cellular_component), and keratin filament binding (molecular_ function). KEGG pathway analysis revealed that there were 20 pathways enriched in these DE mRNAs involved in the lncRNA-miRNAmRNA networks. Of them, chemical carcinogenesis, regulation of actin cytoskeleton, and phosphatidylinositol 3-kinase (PI3K)-Akt signaling pathway were the most enriched. In addition, GO analysis demonstrated that these DE mRNAs in the circRNA-miRNA-mRNA networks are mainly related to bone trabecula formation (biological_process), collagen type XIV trimer (cellular_component), and peptide antigen-transporting ATPase activity (molecular_function). KEGG pathway analysis revealed that there were 61 pathways enriched in these DE mRNAs involved in the circRNA-miRNA-mRNA networks. Of them, the AGE-RAGE (advanced glycation end products-receptor for advanced glycation

Please cite this article in press as: Zhu et al., Analyzing the Interactions of mRNAs and ncRNAs to Predict Competing Endogenous RNA Networks in Osteosarcoma Chemo-Resistance, Molecular Therapy (2019), https://doi.org/10.1016/j.ymthe.2019.01.001

www.moleculartherapy.org

Figure 1. Flow Chart of This Study Design DE, differently expressed; MRE, miRNA response element.

end products) signaling pathway in diabetic complications, osteoclast differentiation, and cGMP-PKG (cyclic guanosine monophosphate-protein kinase G) signaling pathway were the most enriched (Figure 3B).

and circRNA. On the basis of the analysis strategy previously described, we predicted target mRNAs of these miRNAs based on miRanda and compared these target mRNAs with the differentially expressed mRNAs that were identified in our RNA-seq results.

ceRNA Networks Were Constructed on the Basis of the Screened mRNAs and Bioinformatics Prediction

To fully explore the key genes related to drug resistance and make the networks more concise and effective, we further filtered the key DE mRNAs closely related to the drug resistance based on the literature and the results of KEGG pathway analysis in the previous step. Then

As is known to us, the ceRNA regulatory mechanism is very important between the mRNA and ncRNAs, including miRNA, lncRNA,

Molecular Therapy Vol. 27 No 3 March 2019

3

Please cite this article in press as: Zhu et al., Analyzing the Interactions of mRNAs and ncRNAs to Predict Competing Endogenous RNA Networks in Osteosarcoma Chemo-Resistance, Molecular Therapy (2019), https://doi.org/10.1016/j.ymthe.2019.01.001

Molecular Therapy

Figure 2. Expression Profiles of lncRNAs, circRNAs, miRNAs, and mRNAs DE ncRNAs and mRNAs in the three paired chemo-sensitive and chemo-resistant OS cell lines were shown using volcano plot (A) and heatmap (B).

we selected 128 mRNAs in the lncRNA-miRNA-mRNA network, which were involved in the 12 pathways, such as chemical carcinogenesis, metabolism of xenobiotics by cytochrome P450, PI3K-Akt signaling pathway, and so on. Besides, the selected 66 mRNAs in the circRNA-miRNA-mRNA network were involved in the 10 pathways, such as ABC transporters, cGMP-PKG signaling pathway, drug metabolism-cytochrome P450, and so on. Then we constructed the lncRNA/cirRNA-miRNA-mRNA networks based on the selected mRNAs. 2,138 lncRNA-miRNA-mRNA pathways were constructed including 222 lncRNAs, 10 miRNAs, and 33 mRNAs. Furthermore, 3,318 cirRNA-miRNA-mRNA pathways were established, including 111 circRNAs, 13 miRNAs, and 145 mRNAs. Of them, there were 333 lncRNA/cirRNA-miRNA-mRNA pathways that shared the same miRNAs and mRNAs simultaneously, as well as 222 lncRNA/cirRNA-miRNA-mRNA pathways that just shared the same mRNAs. These ceRNA regulatory relationships may be more important in the process of chemo-resistance considering the complexity among the ncRNAs and mRNAs. Owing to the huge and complicated ceRNA networks, we further constructed and illustrated the ceRNA networks on the basis of 11 key DE

4

Molecular Therapy Vol. 27 No 3 March 2019

mRNAs including IKBKB, GSTM2, FGFR4, FGFR1, FGF5, FGF7, FGF10, CSF1, AKT2, MAP4K4, and UGT1A7, which at most occurred in all of the acquired pathways. Then 1,269 lncRNAmiRNA-mRNA pathways were constructed, including 61 lncRNAs, 101 miRNAs, and 11 mRNAs (Figure 4A). 90 circRNA-miRNAmRNA pathways were constructed, including 9 circRNAs, 37 miRNAs, and 10 mRNAs (Figure 4B). Of them, there were 3 circRNAs and 11 lncRNAs that shared the same miRNAs and mRNAs, and 8 circRNAs and 132 lncRNAs that just shared the same mRNAs, not the same miRNAs. All of the circRNA/lncRNAmiRNA-mRNA relationship pairs were shown in Tables S8 and S9. The lncRNA-miRNA-mRNA Pathway Was Verified by the RIP Assay and Dual Luciferase Reporter Gene System

To demonstrate the feasibility and value of the constructed ceRNA networks, we chose the lncRNA MEG3/hsa-miR-200b-3p/AKT2 pathway to validate. For the pathway, the RNA-seq results showed that MEG3 was upregulated with 9-fold change, and AKT2 was upregulated with 13-fold change. We first examined the expression level of MEG3, hsa-miR-200b-3p, and AKT in the chemo-sensitive and chemo-resistant OS cell lines and tissues. The results showed that the mRNA levels of lncRNA MEG3 and AKT2 gene were significantly

Please cite this article in press as: Zhu et al., Analyzing the Interactions of mRNAs and ncRNAs to Predict Competing Endogenous RNA Networks in Osteosarcoma Chemo-Resistance, Molecular Therapy (2019), https://doi.org/10.1016/j.ymthe.2019.01.001

www.moleculartherapy.org

Table 1. Statistical Analysis of All of the Differently Expressed ncRNAs and mRNAs DE RNAs

Total No.

No. of Upregulated

No. of Downregulated

The Most Upregulated (Fold Change)

The Most Downregulated (Fold Change)

lncRNA

339

150

189

lnc-DNMT1-2:1 (85)

lnc-RPLP0-4:1 (55)

circRNA

80

57

23

hsa_circ_0004674 (17)

hsa_circ_0068458 (12)

miRNA

21

9

12

hsa-miR-337-3p (12)

hsa-miR-203a-3p (10)

mRNA

2,606

1,607

999

C6orf106 (124)

AKR1C2 (197)

increased in the chemo-resistant OS cell lines (Figures 5A and 5C) and tissues (Figures 5D and 5F) compared with the controlled. However, hsa-miR-200b-3p expression level was decreased in the chemoresistant OS cell lines (Figure 5B) and tissues (Figure 5E), which was in contrast with the expression of MEG3 and AKT2 (Figure 5G). To further investigate the role of lncRNA MEG3 in OS chemo-resistance, small interfering RNAs (siRNAs) were constructed and transfected into the MG63/DXR (U2OS/DXR) cell lines. qPCR was conducted to examine the transfection effect, and si-MEG3-1 (termed as si-MEG3) with highest inhibition efficacy was chosen for further assays (Figure 5H). The transfection effects of miR-200b mimics and inhibitor were also examined in the MG63/DXR cells (Figure 5I). In the MG63/DXR cell line, AKT2 mRNA and protein expression was significantly decreased when downregulating the lncRNA MEG3 by siRNA, which could be rescued by inhibiting the hsa-miR-200b-3p expression (Figure 5J). Besides, cell viability and IC50 value of OS cells (MG63/DXR) to DXR was regulated by MEG3 and miR-200b-3p, although affecting the AKT2 expression (Figure 5K). Bioinformatics prediction tools further demonstrated that hsa-miR200b-3p targeted 30 UTRs of lncRNA MEG3 with complementary binding site (Figure 5L). Luciferase reporter plasmids containing wild-type (WT) MEG3 and mutant-type (MUT) MEG3 were constructed. Obviously, co-transfection of the luciferase reporter plasmid containing the WT-MEG3 with miR-200b-3p mimics into MG63/ DXR cells resulted in a decreased reporter activity (Figure 5M). Subsequently, RIP assay and RNA pull-down were carried out to determine whether MEG3 can sponge miR-200b-3p. As was shown in Figure 5, lncRNA MEG3 and miR-200b-3p were more abundant in Ago2 pellet compared with in IgG pellet in MG63/DXR cells (Figure 5N). In addition, RNA pull-down assay using biotin-labeled miR-200b-3p (miR-200b-3p-bio) probe increased the expression of lncRNA MEG3 compared with control (NC-bio) or miR-200b-3p-mut-bio probes (Figure 5O). The above results revealed that hsa-miR-200b3p targeted lncRNA MEG3 30 UTR with molecular binding. Besides, bioinformatics prediction tools also illustrated that miR-200b-3p could target 30 UTR of AKT2 mRNA with complementary binding site (Figure 5P). Luciferase reporter gene assay also validated the molecular interaction within miR-200b-3p and AKT2 (Figure 5Q). Taken together, the above results identified the lncRNA MEG3/hsamiR-200b-3p/AKT2 pathway and demonstrated its role in the OS chemo-resistance.

The circRNA-miRNA-mRNA Pathway Was Randomly Chosen and Verified by Dual Luciferase Reporter Gene System

Similarly, the predicted hsa_circ_0001258/hsa-miR-744-3p/GSTM2 pathway was randomly chosen to test the feasibility and value of the current network. We first tested the expression of hsa_ circ_0001258, hsa-miR-744-3p, and GSTM2 in the chemo-resistant and chemo-sensitive OS cell lines and tissues, respectively, and found that hsa_circ_0001258, as well as GSTM2, was downregulated, whereas hsa-miR-744-3p was upregulated in the chemo-resistant OS cell lines and tissues relative to the controlled groups, which is consistent with the results of RNA-seq (Figures 6A–6G). Then, MG63/DXR (or KH-OS/DXR) cells were transfected with the hsa_circ_0001258 expression vector or empty vector, and the transfection effect was further verified by qPCR (Figure 6H). The effect of mimics or inhibitor on the expression of miR-744-3p was also examined in the MG63/DXR cells by qPCR (Figure 6I). We found that, in MG63/DXR cells transfected with hsa_circ_0001258-vector, hsa-miR-744-3p expression was significantly decreased, whereas GSTM2 mRNA expression was increased compared with empty-vector transfection. In MG63/DXR cells transfected with hsa-miR-7443p inhibitor, GSTM2 mRNA and hsa_circ_0001258 expression levels were significantly increased compared with control transfection (Figure 6J). Besides, cell viability and IC50 value of OS cells (MG63/DXR) to DXR was regulated by hsa_circ_0001258 and miR-744-3p, although affecting the GSTM2 expression (Figure 6K). In follow-up assay, bioinformatics prediction tools further demonstrated that hsa_circ_0001258 holds eight binding sites with hsamiR-744-3p (Figure 6L), and hsa-miR-744-3p could target the 30 UTR of GSTM2 mRNA with complementary binding sites (Figure 6M). As was previously described, luciferase reporter gene assay showed that hsa_circ_0001258 could bind with hsa-miR-744-3p, and hsa-miR-744-3p could bind with the 30 UTR of AKT2 mRNA (Figures 6N and 6O). Then we constructed the hsa_circ_0001258/ hsa-miR-744-3p/GSTM2 pathway in the OS chemo-resistance.

DISCUSSION Chemo-resistance is a huge obstacle encountered in OS treatment.2 During the past few decades, regardless of mRNAs, increasing ncRNAs, such as miRNA, lncRNA, circRNA, have been demonstrated to play an essential role in the drug resistance of tumor.23,24 Actually, some pivotal mRNAs and ncRNAs have been reported to be involved in OS chemo-resistance, including RFC related to MTX

Molecular Therapy Vol. 27 No 3 March 2019

5

Please cite this article in press as: Zhu et al., Analyzing the Interactions of mRNAs and ncRNAs to Predict Competing Endogenous RNA Networks in Osteosarcoma Chemo-Resistance, Molecular Therapy (2019), https://doi.org/10.1016/j.ymthe.2019.01.001

Molecular Therapy

Figure 3. Bioinformatics Analysis in the ceRNA Network GO and KEGG pathway analysis for the mRNAs regulated by the lncRNA-miRNA network (A) and circRNA-miRNA network (B).

resistance,25 miR-33a related to DDP resistance,26 and lncRNA FOXC2-AS1 related to DXR resistance.27 As is known to us, most ncRNAs function through regulating the expression of some vital mRNAs and the ceRNA regulatory networks among these so many coding and ncRNAs were very complex and intertwined. Actually, there have been several reports about the comprehensive expression profile and related ceRNA networks in the carcinogenesis or differentiation. Huang et al.28 reported the comprehensive analysis of differentially expressed profiles of lncRNAs and circRNAs with associated co-expression and ceRNA networks in bladder carcinoma. Yuan et al.29 reported the interactions of mRNAs, miRNAs, lncRNAs, and circRNAs to predict ceRNA networks in glioblastoma. Xue et al.30 constructed an esophageal cancerspecific ceRNA network based on miRNA, lncRNA, and mRNA expression data. Besides, Dou et al.31 showed the changing expression profiles of lncRNAs, mRNAs, circRNAs, and miRNAs during osteoclastogenesis. Gu et al.32 identified and integrated analysis of differentially expressed lncRNAs and circRNAs to reveal the potential ceRNA networks during periodontal ligament stem cell osteogenic differentiation. However, there is seldom a report about the comprehensive expression profile and analysis of mRNAs and these ncRNAs involved in the multi-drug chemo-resistance of OS. In the current study, we simultaneously identified the lncRNA, circRNA, miRNA,

6

Molecular Therapy Vol. 27 No 3 March 2019

and mRNA expression profile in the three paired multi-drug chemo-resistant and chemo-sensitive OS cell lines by the whole-transcriptome RNA-seq. We found that there were 2,718 lncRNAs, 80 circRNAs, 35 miRNAs, and 3,389 mRNAs DE with 2-fold change and p < 0.05. The subsequent GO and KEGG pathway analysis showed that these DE molecules mainly related to the chemical carcinogenesis, osteoclast differentiation, PI3K-AKT pathway, and cGMPPKG signaling pathway. Salmena et al.33 first proposed the “ceRNA hypothesis” describing how mRNAs, transcribed pseudogenes, and lncRNAs “talk” to each other using MREs as letters of a new language. Many lncRNAs and circRNAs have been reported to function as ceRNAs to sponge miRNAs to further regulate the expression of mRNAs and their related proteins in the regulation of carcinogenesis and tumor progression, including chemo-resistance. Han et al.34 found that lncRNA CRNDE promotes colorectal cancer cell proliferation and chemoresistance via miR-181a-5p-mediated regulation of Wnt/b-catenin signaling. Pan et al.35 reported that lncRNA UCA1 promotes DDP or gemcitabine resistance through CREB modulating miR-196a-5p in bladder cancer cells. Besides, Han et al.36 found that lncRNA LUCAT1 modulates MTX resistance in OS via the miR-200c/ ABCB1 axis. Gao et al.22 reported that the circ_0006528-miR-7-5pRaf1 axis may play a role in the DXR resistance of breast cancer.

Please cite this article in press as: Zhu et al., Analyzing the Interactions of mRNAs and ncRNAs to Predict Competing Endogenous RNA Networks in Osteosarcoma Chemo-Resistance, Molecular Therapy (2019), https://doi.org/10.1016/j.ymthe.2019.01.001

www.moleculartherapy.org

Figure 4. Constructed ceRNA Networks Based on the Selected Vital Genes Related to Drug Resistance lncRNA-miRNA-mRNA network (A) and circRNA-miRNA-mRNA network (B) were shown, respectively. The green nodes represented mRNAs, blue nodes lncRNAs or circRNAs, red frames miRNAs, and the edges represented the competing interactions among them.

In the present study, we provided a systematic perspective on the potential function of ncRNAs in the multi-drug chemo-resistance of OS. The lncRNA- or circRNA-miRNA-mRNA regulatory networks in the multi-drug chemo-resistance of OS were comprehensively integrated and predicted on the basis of the screening and bioinformatics analysis. The selected mRNAs in the networks conformed to the ceRNA rules and were also included in the well-known drug-resistance-related genes or signaling pathways. Then we constructed the ceRNA regulatory networks including 388 lncRNAs, 28 circRNAs, 15 miRNAs, and 11 mRNAs. Of them, there are 533 lncRNA-miRNA-mRNA and another 388 circRNA-miRNAmRNA pathways. Moreover, there are 20 lncRNAs and 10 circRNAs with shared miRNA and mRNA, 30 lncRNAs and 20 circRNAs with shared miRNAs not mRNA, and 35 lncRNAs and 22 circRNAs with shared mRNAs not miRNAs. These intertwined ceRNA networks may suggest a more important role in the regulation. In addition, to further demonstrate the credibility of our results, two predicted pathways, lncRNAMEG3/hsa-miR-200b-3p/AKT2 and hsa_circ_ 0001258/hsa-miR-744-3p/GSTM2, were selected and verified by the PCR, RIP, RNA pull-down, and dual luciferase reporter gene system. Actually, there have been several reports about lncRNA MEG3, miR200b-3p, and AKT2 in OS progression or cancer chemo-resistance. Zhang et al.37 previously reported that lncRNA MEG3 was upregulated in the OS cells and tissues, promoted OS cell growth and metas-

tasis in vitro through sponging miR-127, and may be a potential therapeutic target for OS. Besides, Lv et al.38 found that disruption of the c-Myc-miR-200b-3p-PRDX2 regulatory loop enhances tumor metastasis and chemotherapeutic resistance in colorectal cancer. In addition, Zhu et al.39 showed that elevated expression of AKT2 correlates with disease severity and poor prognosis in human OS, and Zhang et al.40 reported that knockdown of AKT2 sensitizes OS cells to apoptosis induced by DDP treatment. In the current study, we found that lncRNA MEG3 and AKT2 expression was simultaneously upregulated, whereas miR-200b-3p expression was downregulated in the chemo-resistant OS cell lines and tissues compared with the controlled. Function assay revealed that lncRNA MEG3 could promote the OS DXR resistance through upregulating the expression of AKT2, which could be inhibited by overexpressing the miR200b-3p expression. Mechanism investigation demonstrated that lncRNA MEG3 as a ceRNA upregulated the AKT2 expression by sponging miR-200b-3p. Then we constructed the lncRNA MEG3/ hsa-miR-200b-3p/AKT2 pathway and identified its role in OS chemo-resistance. However, there has been no report about hsa_circ_0001258 or hsa-miR-744-3p related to OS or cancer chemo-resistance. We first reported the expression, function, and regulatory mechanism of hsa_circ_0001258 and hsa-miR-744-3p in OS chemo-resistance. We found that hsa_circ_0001258, as well as GSTM2, was downregulated, whereas hsa-miR-744-3p was upregulated in the

Molecular Therapy Vol. 27 No 3 March 2019

7

Please cite this article in press as: Zhu et al., Analyzing the Interactions of mRNAs and ncRNAs to Predict Competing Endogenous RNA Networks in Osteosarcoma Chemo-Resistance, Molecular Therapy (2019), https://doi.org/10.1016/j.ymthe.2019.01.001

Molecular Therapy

Figure 5. lncRNAMEG3/hsa-miR-200b-3p/AKT2 Pathway Was Randomly Chosen to Verify Feasibility of the Constructed lncRNA-miRNA-mRNA Network (A–F) The mRNA expression levels of MEG3 (A and D), miR-200b-3p (B and E), and AKT2 (C and F) were examined by real-time qPCR in OS cells and tissues. (G) The correlation between MEG3 and miR-200b-3p, as well as miR-200b-3p and AKT2, in OS tissues. (H) Real-time qPCR analysis of the effect on knockdown of MEG3 expression by siRNA in the MG63/DXR and KH-OS/DXR cell lines. si-MEG3-1 was chosen for the siRNA used in the study because of the highest knockdown efficiency. (I) Real-time qPCR analysis of the effect on mimics or inhibitor of miR-200b-3p in the MG63/DXR cells. (J) AKT2 expression was regulated by the MEG3 and miR-200b-3p. (K) Cell viability and IC50 value of OS cells to doxorubicin were regulated by MEG3 and miR-200b-3p, although affecting the AKT2 expression. (L) Potential binding sites of lncRNA MEG3 and miR-200b-3p. (M) Luciferase activity assay showed that MEG3 could combine with miR-200b-3p. (N) RIP assays showed that MEG3 and miR-200b-3p both could bind to the Ago2 protein. (O) RNA pull-down assay showed that using biotin-labeled miR-200b-3p (miR-200b-3p-bio) probe increased the expression of lncRNA MEG3 compared with control (NC-bio) or miR-200b-3p-Mut-Bio probes. (P) Potential binding sites of miR-200b-3p and 30 UTR of AKT2 mRNA. (Q) Luciferase activity assay showed that miR-200b-3p could combine with 30 UTR of AKT2 mRNA.

8

Molecular Therapy Vol. 27 No 3 March 2019

Please cite this article in press as: Zhu et al., Analyzing the Interactions of mRNAs and ncRNAs to Predict Competing Endogenous RNA Networks in Osteosarcoma Chemo-Resistance, Molecular Therapy (2019), https://doi.org/10.1016/j.ymthe.2019.01.001

www.moleculartherapy.org

Figure 6. The hsa_circ_0001258/hsa-miR-744-3p/GSTM2 Pathway Was Randomly Chosen to Verify Feasibility of the Constructed circRNA-miRNA-mRNA Network (A–F) The mRNA expression levels of hsa_circ_0001258 (A and D), miR-744-3p (B and E), and GSTM2 (C and F) were examined by the real-time qPCR in OS cells and tissues. (G) The correlation between hsa_circ_0001258 and miR-744-3p, as well as miR-744-3p and GSTM2, in OS tissues. (H) Real-time qPCR analysis of the effect on overexpression of hsa_circ_0001258 by vector transfection in the MG63/DXR and KH-OS/DXR cell lines. (I) Real-time qPCR analysis of the effect on mimics or inhibitor of miR-744-3p in the MG63/DXR cells. (J) GSTM2 expression was regulated by hsa_circ_0001258 and miR-744-3p. (K) Cell viability and IC50 value of OS cells to doxorubicin was regulated by hsa_circ_0001258 and miR-744-3p, although affecting the GSTM2 expression. (L) Potential binding sites of hsa_circ_0001258 and miR-744-3p. (M) Luciferase activity assay showed that hsa_circ_0001258 could combine with miR-744-3p. (N) Potential binding sites of miR-744-3p and 30 UTR of GSTM2 mRNA. (O) Luciferase activity assay showed that miR-744-3p could combine with the 30 UTR of GSTM2 mRNA.

Molecular Therapy Vol. 27 No 3 March 2019

9

Please cite this article in press as: Zhu et al., Analyzing the Interactions of mRNAs and ncRNAs to Predict Competing Endogenous RNA Networks in Osteosarcoma Chemo-Resistance, Molecular Therapy (2019), https://doi.org/10.1016/j.ymthe.2019.01.001

Molecular Therapy

Table 2. Clinical Parameters of Osteosarcoma Patients Enrolled in This Study Chemo-Sensitive Group

Chemo-Resistant Group

p Value

48

32

–

Male

30

26

Female

18

6

Age, y (range)

17.2 ± 0.8 (6–32)

18.4 ± 0.6 (8–30)

Proximal of femur

8

4

Distal of femur

20

16

Proximal of tibia

14

11

Other

6

1

Follow-up time, months (range)

28.0 ± 1.2 (6–102)

27.4 ± 1.4 (3–96)

Number Gender

0.08 0.1

Location

0.07

To the best of our knowledge, this is the first report examining the changing expression of lncRNAs, circRNAs miRNAs, and mRNAs in the OS chemo-resistance. In our study, we acquired the comprehensive expression profile detected by the whole-transcriptome RNA-seq in three paired multi-drug chemo-resistant and chemo-sensitive OS cell lines. GO and KEGG pathway analyses allowed us to annotate the potential functions of differentially expressed mRNAs regulated by the lncRNA- or circRNA-miRNA. ceRNA regulatory networks were constructed for both lncRNA-miRNA-mRNA and circRNA-miRNA-mRNA interactions, which contributes to our understanding of the mechanism of OS chemo-resistance and finding some novel targets reversing it.

MATERIALS AND METHODS Cell Lines

Three human OS cell lines (MG63, KHOS, U2OS) were purchased from American Type Culture Collection and cultured in DMEM supplemented with 10% FBS (GIBCO, Gran Island, NY, USA), 100 U/mL penicillin, and 100 mg/mL streptomycin (Invitrogen) at 37 C in a humidified CO2 (5%) atmosphere. The paired three DXR-resistant OS cell lines (MG63/DXR, KHOS/DXR, U2OS/DXR) were kindly donated by Dr. Oda Yoshio, Dr. E.S. Gonos, and Dr. Z.F. Duan, respectively, which have been described in our previous study.44

Molecular Therapy Vol. 27 No 3 March 2019

A total of 80 primary OS patients who received the same chemotherapy regimen (MTX, DXR, DDP, and ifosfamide [IFO] for two cycles) before surgery and underwent complete resection surgery at Shanghai Tenth People’s Hospital between 2006 and 2016 were included in this study. The study was approved by the Ethics Committee of Shanghai Tenth People’s Hospital, and written informed consent was obtained from all of the patients. All patients’ slides were reviewed to confirm the diagnosis by the same pathologist. All of the resected specimens were placed immediately into liquid nitrogen and stored at 80 C. According to the Huvos scoring system,13 the patients were classified as chemo-resistant and chemo-sensitive groups. The clinical parameters of OS patients in this study are shown in Table 2. CCK-8 Assay

0.55

chemo-resistant OS cell lines and tissues relative to the controlled groups. Further study showed that hsa_circ_0001258 could serve as a ceRNA to sponge the hsa-miR-744-3p to regulate the expression of GSTM2. Actually, glutathione S-transferase mu 2 was encoded by GSTM2, which could detoxify electrophilic compounds, including carcinogens, therapeutic drugs, environmental toxins, and products of oxidative stress, by conjugation with glutathione.41–43 Then we found that hsa_circ_0001258 suppressed the DXR resistance of OS cells through upregulating the GSTM2 expression via sponging hsa-miR-744-3p.

10

Human Tissue Samples

After 48 h of transfection in 96-well plates, freshly prepared medium containing several final concentrations of DXR (0, 2, 4, 8, 16, 32, and 64 mg/mL) was added to the wells, with three replicate wells for each concentration. After incubating for another 48 h, cell viability was measured using Cell Counting Kit-8 (CCK-8; Dojindo, Japan), according to the manufacturer’s instructions. A microplate reader (Biotek, Winooski, VT, USA) was used to quantify the absorbance of each well at the wavelength of 450 nm. All reactions were repeated in triplicate. RNA Extraction and Quality Control

Total RNA was isolated from cells and tissue using TRIzol reagent (Life Technologies, CA, USA) according to the manufacturer’s protocol. The quantity and quality of total RNA samples were measured using NanoDrop ND-1000 (Wilmington, DE, USA). RNA integrity was evaluated using the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). The samples with RNA Integrity Number (RIN) R7 were subjected to the subsequent analysis. Library Preparation for Small RNA-Seq

A total of 3 mg of total RNA per sample was used as input material for the small RNA library. Sequencing libraries were generated using NEBNextR Multiplex Small RNA Library Prep Set for IlluminaR (NEB, USA) following the manufacturer’s recommendations, and index codes were added to attribute sequences to each sample. Briefly, the NEB 30 SR adaptor was directly and specifically ligated to the 30 end of miRNA. After the 30 ligation reaction, the SR RT primer hybridized to the excess of the 30 SR adaptor (that remained free after the 30 ligation reaction) and transformed the single-stranded DNA adaptor into a double-stranded DNA molecule. This step is important to prevent adaptor-dimer formation. In addition, dsDNAs are not substrates for ligation mediated by T4 RNA Ligase 1, and thus do not ligate to the 50 SR adaptor in the subsequent ligation step. Next, the 50 end adaptor was ligated to 50 ends of miRNAs. First-strand cDNA was synthesized using M-MuLV Reverse Transcriptase (RNase H). PCR amplification was performed using LongAmp Taq 2X Master Mix, SR Primer for Illumina, and index (X) primer. PCR products were purified on an 8% polyacrylamide gel (100 V,

Please cite this article in press as: Zhu et al., Analyzing the Interactions of mRNAs and ncRNAs to Predict Competing Endogenous RNA Networks in Osteosarcoma Chemo-Resistance, Molecular Therapy (2019), https://doi.org/10.1016/j.ymthe.2019.01.001

www.moleculartherapy.org

80 min). DNA fragments corresponding to 140–160 bp (the length of small noncoding RNA plus the 3 and 5 adaptors) were recovered and dissolved in 8 mL of elution buffer. At last, library quality was assessed on the Agilent Bioanalyzer 2100 system using DNA High Sensitivity Chips.

coefficient) and the predicted pairs of lncRNA(circRNA)-mRNA based on ceRNA score principle were then considered as the true ceRNAs. Based on established co-expression data, deregulated ceRNA networks and circRNA, lncRNA, or miRNA interactions of interest were mapped using Cytoscape software (V. 3.2.1).

Construction of cDNA Libraries and High-Throughput Sequencing

Real-Time qPCR

Total RNA was extracted from two groups using TRIzol (Invitrogen) according to the manufacturer’s protocols. Strand-specific cDNA libraries were constructed following a previously described protocol and were sequenced using an Illumina HiSeq2000 sequencer (OEbiotech, Shanghai, China) by performing a paired-end run with a 100-bp read length. The raw reads were processed by removing the adaptor reads and low-quality tags. All subsequent analyses were performed using clean reads. Clean paired-end reads were aligned to the reference genome using TopHat (version 2.0.6). The transcriptome of each sample was constructed using Cufflinks (version 2.0). The reads per kilobase of model per million base pairs sequenced (RPKM) was used to quantify the expression levels of a gene or lncRNA (circRNA), and transcript per million (TPM) was used to determine the expression levels of miRNA. The difference between different groups was determined by DEGseq software. GO and Pathway Analysis

GO analysis is a functional analysis associating differentially expressed mRNAs with GO categories. The GO categories are derived from Gene Ontology (http://www.geneontology.org), which comprise three structured networks of defined terms that describe gene product attributes. Pathway analysis for differentially expressed mRNAs was provided, based on the latest KEGG (Kyoto Encyclopedia of Genes and Genomes; https://www.genome.jp/kegg) database, which allowed us to determine the biological pathway with the significantly enriched mRNAs involved. Overview of the Processes Used to Identify ceRNA Interaction Pairs

Based on the expression levels of mRNAs, miRNAs, and lncRNAs or circRNAs, Pearson’s correlation coefficient and p value were calculated for miRNA-target. Negatively correlated pairs with a p value <0.05 and Pearson’s correlation coefficient >0.8 were selected for further analysis. These were the predicted pairs for miRNAmRNA, miRNA-lncRNA(circRNA), and mRNA-lncRNA(circRNA). Subsequently, shared pairs from the predicted pairs from binding sites and the predicted pairs from the expression of mRNA, lncRNA (circRNA), and miRNA were used for further analysis. Based on the principle for ceRNA prediction, shared pairs of miRNA-mRNA and miRNA-lncRNA(circRNA) were used to predict ceRNA score according to the following formula: ceRNA_score = MRE_for_share_ miRNA/MRE_for_lncRNA(circRNA)_miRNA. The shared pairs from the predicted pairs of lncRNA(circRNA)-mRNA based on the expression of lncRNA(circRNA) and mRNA (Pearson’s correlation

Total RNA was isolated from OS cells and tissues using TRIzol reagent (Invitrogen). And then, the cDNA was synthesized from 200 ng of extracted total RNA using the PrimeScript RT reagent Kit (Takara Bio Company, Shiga, Japan) and amplified by real-time qPCR with an SYBR Green Kit (Takara Bio Company) on an ABIPRISM 7500 Sequence Detection System (Applied Biosystems, Foster City, CA, USA) with the housekeeping gene glyceraldehyde-3-phosphate dehydrogenase (GAPDH) as an internal control. The 2DDCt method was used to determine the relative quantification of gene expression levels. All of the primers were synthesized by Sangon (Shanghai, China), and the sequences of primers were shown in Table S2. Cell Transfection

The small interfering (si) RNA targeting against lncRNA MEG3 and corresponding controls (si-NC), as well as mimics and inhibitor of miR-200b-3p or miR-744-3p, were designed and synthesized by GenePharma (Shanghai, China). The sequence of hsa_circ_ 0001258 was acquired from the circBase (http://www.circbase.org/ cgi-bin/singlerecord.cgi?id=hsa_circ_0001258). MG-63 (KH-OS) and MG63/DXR (KH-OS/DXR) cells were transfected using Lipofectamine 2000 (Invitrogen) according to the manufacturer’s instruction. The designed sequences were shown in Table S3. Western Blot

Total proteins were extracted from OS cells and then transferred into 10% SDS-PAGE. Separation of protein was then transferred onto polyvinylidene difluoride (PVDF) membranes and blocked using 5% no-fat milk for 1 h. Then, PVDF membranes were incubated overnight with primary antibodies at 4 C (1:1,000 dilution; Santa Cruz), including anti-AKT2 and anti-GSTM2. Then, secondary antibody (goat anti-rabbit IgG, 1:500) was incubated for 1 h. Finally, protein expression was determined using a chemiluminescence detection kit (Amersham Pharmacia Biotech) correlated to GAPDH endogenous control. RNA Immunoprecipitation Assays

RIP experiments were performed according to a previous report by using a Magna RIP RNA-Binding Protein Immunoprecipitation Kit (Millipore, USA) according to the manufacturer’s instructions. The antibody used for RIP assays was Ago2 (Cell Signaling Technology, USA). When RIP was treated with RNase, lysates were incubated with RNase for 1 h in 37 C. The co-precipitated RNAs were detected by real-time PCR. RNA Pull-Down

RNA pull-down assay was conducted to determine the interaction between lncRNA MEG3 and miR-200b-3p. In brief, the miR-200b-3p or

Molecular Therapy Vol. 27 No 3 March 2019

11

Please cite this article in press as: Zhu et al., Analyzing the Interactions of mRNAs and ncRNAs to Predict Competing Endogenous RNA Networks in Osteosarcoma Chemo-Resistance, Molecular Therapy (2019), https://doi.org/10.1016/j.ymthe.2019.01.001

Molecular Therapy

miR-200b-3p-Mut probe was synthesized and biotinylated by GenePharma (Shanghai, China). RNA pull-down assay was carried out using the Magnetic RNA-Protein Pull-Down Kit (Thermo Fisher Scientific, Waltham, MA, USA) according to its specifications. First, 3 mg of miR-200b-3p that was labeled with streptavidin was diluted into 100 mL with structure buffer and then kept at 90 C for 5 min. Then, the streptavidin beads were washed with 60 mL of cell lysate and incubated at temperature for 1 h with agitation. After the pretreated RNA was incubated at 90 C for 1 h with agitation, the 60 mL of washed streptavidin beads was added for incubation at 90 C for 1 h with agitation. The beads were washed twice with 1 mL of high-concentration salt solution and subsequently washed twice with 1 mL of high-concentration salt solution. Finally, the RNA-binding protein complexes were washed and eluted for realtime qPCR analysis.

calculated to correct the p value on RNA-sequencing analysis. The expression level of each mRNA, miRNA, lncRNA, and circRNA was represented as fold change using the 2DDCt method on realtime qPCR analysis. p values <0.05 were considered statistically significant.

Luciferase Gene Reporter Assay

CONFLICTS OF INTEREST

MG63/DXR cells were plated in 96-well plates at a density of 5  104 cells/well and then transfected. The 30 UTR of lncRNA MEG3 cDNA and AKT2 were amplified using PCR amplification and then subcloned into the luciferase reporter vector pGL3-basic (Promega, Madison, WI, USA). Renilla luciferase acted as control reporter. miR200b-3p mimics or negative controls (100 nM) were transfected using Lipofectamine 2000 (Invitrogen, CA, USA) according to the manufacturer’s instructions. Cells were lysed and assayed for luciferase activity for 36 h. The values were normalized to the Renilla luciferase transfection control using the dual luciferase reporter gene assay kit (Promega) according to the manufacturer’s instructions. The same procedure was suitable to validate the relationship between hsa_ circ_0001258 and miR-744-3p or miR-744-3p and GSTM2. Analysis Strategy and Procedure of ceRNA Network Construction

1. All of the DE lncRNAs, circRNAs, and mRNAs (fold change > 2.0, p < 0.05, and FDR < 0.05), combined with all of the miRNAs (differently or non-DE). 2. Retain the regulatory relationship of lncRNA or circRNA (upregulation)-miRNAs (downregulation or invariably)-mRNA (upregulation) and lncRNA or circRNA (downregulation)-miRNAs (upregulation or invariably)-mRNA (downregulation). 3. There was a co-expression positive relationship between the differential lncRNA or circRNA and mRNA, and the correlation coefficient was >0.8. 4. lncRNA or circRNA and miRNA, as well as miRNA and mRNA, have complementary sequences, and there are at least three MREs between the lncRNA or circRNA and miRNA. 5. lncRNA or circRNA-miRNA-mRNA pathways with shared miRNAs and/or mRNAs were listed separately. Statistical Analyses

A statistical analysis was performed using Student’s t test to compare two variables of the RNA-sequencing data. The differences with fold change (FC) R2 or %0.5 and p < 0.05 were considered to be statistically significant on RNA-sequencing analysis. The FDR was

12

Molecular Therapy Vol. 27 No 3 March 2019

SUPPLEMENTAL INFORMATION Supplemental Information includes nine tables and can be found with this article online at https://doi.org/10.1016/j.ymthe.2019.01.001.

AUTHOR CONTRIBUTIONS K.-P.Z. and X.-L.M. carried out the sample collection. J.-P.H. and T.C. carried out the cell experiments. K.-P.Z., C.-L.Z., and L.Z. participated in the design of the study and performed the statistical analysis. K.-P.Z. drafted the manuscript, and C.-L.Z. helped to correct it.

The authors declare no competing interests.

ACKNOWLEDGMENTS This project was supported by grants from the National Natural Science Foundation of China (81572630 and 81872174), Fundamental Research Funds for the Central Universities (22120170216), the Shanghai Pujiang Program of Shanghai Science and Technology Commission (13PJD023), and Shanghai Jiaotong University Medical-Engineering Cross Research Fund (YG2012MS49).

REFERENCES 1. Saraf, A.J., Fenger, J.M., and Roberts, R.D. (2018). Osteosarcoma: accelerating progress makes for a hopeful future. Front. Oncol. 8, 4. 2. Chou, A.J., and Gorlick, R. (2006). Chemotherapy resistance in osteosarcoma: current challenges and future directions. Expert Rev. Anticancer Ther. 6, 1075–1085. 3. Kim, M., and Kim, D.J. (2018). GFRA1: a novel molecular target for the prevention of osteosarcoma chemoresistance. Int. J. Mol. Sci. 19, e1078. 4. Yan, G.N., Lv, Y.F., and Guo, Q.N. (2016). Advances in osteosarcoma stem cell research and opportunities for novel therapeutic targets. Cancer Lett. 370, 268–274. 5. Cech, T.R., and Steitz, J.A. (2014). The noncoding RNA revolution-trashing old rules to forge new ones. Cell 157, 77–94. 6. George, J., and Patel, T. (2015). Noncoding RNA as therapeutic targets for hepatocellular carcinoma. Semin. Liver Dis. 35, 63–74. 7. Wright, M.W., and Bruford, E.A. (2011). Naming ‘junk’: human non-protein coding RNA (ncRNA) gene nomenclature. Hum. Genomics 5, 90–98. 8. Hombach, S., and Kretz, M. (2016). Non-coding RNAs: classification, biology and functioning. Adv. Exp. Med. Biol. 937, 3–17. 9. Huang, T., Alvarez, A., Hu, B., and Cheng, S.Y. (2013). Noncoding RNAs in cancer and cancer stem cells. Chin. J. Cancer 32, 582–593. 10. Li, C.H., and Chen, Y. (2013). Targeting long non-coding RNAs in cancers: progress and prospects. Int. J. Biochem. Cell Biol. 45, 1895–1910. 11. Venkatesh, T., Suresh, P.S., and Tsutsumi, R. (2015). Non-coding RNAs: functions and applications in endocrine-related cancer. Mol. Cell. Endocrinol. 416, 88–96. 12. Pu, Y., Zhao, F., Li, Y., Cui, M., Wang, H., Meng, X., and Cai, S. (2017). The miR-34a5p promotes the multi-chemoresistance of osteosarcoma via repression of the AGTR1 gene. BMC Cancer 17, 45. 13. Xu, M., Jin, H., Xu, C.X., Bi, W.Z., and Wang, Y. (2014). MiR-34c inhibits osteosarcoma metastasis and chemoresistance. Med. Oncol. 31, 972.

Please cite this article in press as: Zhu et al., Analyzing the Interactions of mRNAs and ncRNAs to Predict Competing Endogenous RNA Networks in Osteosarcoma Chemo-Resistance, Molecular Therapy (2019), https://doi.org/10.1016/j.ymthe.2019.01.001

www.moleculartherapy.org

14. Wang, Z., Liu, Z., and Wu, S. (2017). Long non-coding RNA CTA sensitizes osteosarcoma cells to doxorubicin through inhibition of autophagy. Oncotarget 8, 31465–31477.

30. Xue, W.H., Fan, Z.R., Li, L.F., Lu, J.L., Ma, B.J., Kan, Q.C., and Zhao, J. (2018). Construction of an oesophageal cancer-specific ceRNA network based on miRNA, lncRNA, and mRNA expression data. World J. Gastroenterol. 24, 23–34.

15. Kun-Peng, Z., Xiao-Long, M., and Chun-Lin, Z. (2017). LncRNA FENDRR sensitizes doxorubicin-resistance of osteosarcoma cells through down-regulating ABCB1 and ABCC1. Oncotarget 8, 71881–71893.

31. Dou, C., Cao, Z., Yang, B., Ding, N., Hou, T., Luo, F., Kang, F., Li, J., Yang, X., Jiang, H., et al. (2016). Changing expression profiles of lncRNAs, mRNAs, circRNAs and miRNAs during osteoclastogenesis. Sci. Rep. 6, 21499.

16. Kartha, R.V., and Subramanian, S. (2014). Competing endogenous RNAs (ceRNAs): new entrants to the intricacies of gene regulation. Front. Genet. 5, 8.

32. Gu, X., Li, M., Jin, Y., Liu, D., and Wei, F. (2017). Identification and integrated analysis of differentially expressed lncRNAs and circRNAs reveal the potential ceRNA networks during PDLSC osteogenic differentiation. BMC Genet. 18, 100.

17. Sen, R., Ghosal, S., Das, S., Balti, S., and Chakrabarti, J. (2014). Competing endogenous RNA: the key to posttranscriptional regulation. Sci. World J. 2014, 896206. 18. Karreth, F.A., and Pandolfi, P.P. (2013). ceRNA cross-talk in cancer: when ce-bling rivalries go awry. Cancer Discov. 3, 1113–1121. 19. Qi, X., Zhang, D.H., Wu, N., Xiao, J.H., Wang, X., and Ma, W. (2015). ceRNA in cancer: possible functions and clinical implications. J. Med. Genet. 52, 710–718.

33. Salmena, L., Poliseno, L., Tay, Y., Kats, L., and Pandolfi, P.P. (2011). A ceRNA hypothesis: the Rosetta Stone of a hidden RNA language? Cell 146, 353–358. 34. Han, P., Li, J.W., Zhang, B.M., Lv, J.C., Li, Y.M., Gu, X.Y., Yu, Z.W., Jia, Y.H., Bai, X.F., Li, L., et al. (2017). The lncRNA CRNDE promotes colorectal cancer cell proliferation and chemoresistance via miR-181a-5p-mediated regulation of Wnt/b-catenin signaling. Mol. Cancer 16, 9.

20. An, J., Lv, W., and Zhang, Y. (2017). LncRNA NEAT1 contributes to paclitaxel resistance of ovarian cancer cells by regulating ZEB1 expression via miR-194. OncoTargets Ther. 10, 5377–5390.

35. Pan, J., Li, X., Wu, W., Xue, M., Hou, H., Zhai, W., and Chen, W. (2016). Long noncoding RNA UCA1 promotes cisplatin/gemcitabine resistance through CREB modulating miR-196a-5p in bladder cancer cells. Cancer Lett. 382, 64–76.

21. Gao, W., Lin, S., Cheng, C., Zhu, A., Hu, Y., Shi, Z., Zhang, X., and Hong, Z. (2018). Long non-coding RNA CASC2 regulates Sprouty2 via functioning as a competing endogenous RNA for miR-183 to modulate the sensitivity of prostate cancer cells to docetaxel. Arch. Biochem. Biophys. Published online January 31, 2018. https:// doi.org/10.1016/j.abb.2018.01.013.

36. Han, Z., and Shi, L. (2018). Long non-coding RNA LUCAT1 modulates methotrexate resistance in osteosarcoma via miR-200c/ABCB1 axis. Biochem. Biophys. Res. Commun. 495, 947–953.

22. Gao, D., Zhang, X., Liu, B., Meng, D., Fang, K., Guo, Z., and Li, L. (2017). Screening circular RNA related to chemotherapeutic resistance in breast cancer. Epigenomics 9, 1175–1188. 23. Fanale, D., Castiglia, M., Bazan, V., and Russo, A. (2016). Involvement of non-coding RNAs in chemo- and radioresistance of colorectal cancer. Adv. Exp. Med. Biol. 937, 207–228. 24. Zebisch, A., Hatzl, S., Pichler, M., Wölfler, A., and Sill, H. (2016). Therapeutic resistance in acute myeloid leukemia: the role of non-coding RNAs. Int. J. Mol. Sci. 17, e2080. 25. Wang, J.J., and Li, G.J. (2014). Relationship between RFC gene expression and intracellular drug concentration in methotrexate-resistant osteosarcoma cells. Genet. Mol. Res. 13, 5313–5321. 26. Zhou, Y., Huang, Z., Wu, S., Zang, X., Liu, M., and Shi, J. (2014). miR-33a is up-regulated in chemoresistant osteosarcoma and promotes osteosarcoma cell resistance to cisplatin by down-regulating TWIST. J. Exp. Clin. Cancer Res 33, 12.

37. Zhang, S.Z., Cai, L., and Li, B. (2017). MEG3 long non-coding RNA prevents cell growth and metastasis of osteosarcoma. Bratisl. Lek Listy 118, 632–636. 38. Lv, Z., Wei, J., You, W., Wang, R., Shang, J., Xiong, Y., Yang, H., Yang, X., and Fu, Z. (2017). Disruption of the c-Myc/miR-200b-3p/PRDX2 regulatory loop enhances tumor metastasis and chemotherapeutic resistance in colorectal cancer. J. Transl. Med. 15, 257. 39. Zhu, Y., Zhou, J., Ji, Y., and Yu, B. (2014). Elevated expression of AKT2 correlates with disease severity and poor prognosis in human osteosarcoma. Mol. Med. Rep. 10, 737–742. 40. Zhang, G., Li, M., Zhu, X., Bai, Y., and Yang, C. (2011). Knockdown of Akt sensitizes osteosarcoma cells to apoptosis induced by cisplatin treatment. Int. J. Mol. Sci. 12, 2994–3005. 41. Bräutigam, M., Teusch, N., Schenk, T., Sheikh, M., Aricioglu, R.Z., Borowski, S.H., Neudörfl, J.M., Baumann, U., Griesbeck, A.G., and Pietsch, M. (2015). Selective inhibitors of glutathione transferase P1 with trioxane structure as anticancer agents. ChemMedChem 10, 629–639.

27. Zhang, C.L., Zhu, K.P., and Ma, X.L. (2017). Antisense lncRNA FOXC2-AS1 promotes doxorubicin resistance in osteosarcoma by increasing the expression of FOXC2. Cancer Lett. 396, 66–75.

42. Cacciatore, I., Caccuri, A.M., Cocco, A., De Maria, F., Di Stefano, A., Luisi, G., Pinnen, F., Ricci, G., Sozio, P., and Turella, P. (2005). Potent isozyme-selective inhibition of human glutathione S-transferase A1-1 by a novel glutathione S-conjugate. Amino Acids 29, 255–261.

28. Huang, M., Zhong, Z., Lv, M., Shu, J., Tian, Q., and Chen, J. (2016). Comprehensive analysis of differentially expressed profiles of lncRNAs and circRNAs with associated co-expression and ceRNA networks in bladder carcinoma. Oncotarget 7, 47186– 47200.

43. Bhattacharjee, P., Paul, S., Banerjee, M., Patra, D., Banerjee, P., Ghoshal, N., Bandyopadhyay, A., and Giri, A.K. (2013). Functional compensation of glutathione S-transferase M1 (GSTM1) null by another GST superfamily member, GSTM2. Sci. Rep. 3, 2704.

29. Yuan, Y., Jiaoming, L., Xiang, W., Yanhui, L., Shu, J., Maling, G., and Qing, M. (2018). Analyzing the interactions of mRNAs, miRNAs, lncRNAs and circRNAs to predict competing endogenous RNA networks in glioblastoma. J. Neurooncol. 137, 493–502.

44. Kun-Peng, Z., Xiao-Long, M., and Chun-Lin, Z. (2018). Overexpressed circPVT1, a potential new circular RNA biomarker, contributes to doxorubicin and cisplatin resistance of osteosarcoma cells by regulating ABCB1. Int. J. Biol. Sci. 14, 321–330.

Molecular Therapy Vol. 27 No 3 March 2019

13