Expression Profile Screening and Bioinformatics Analysis of Circrna, Lncrna and Mrna in Acute Myeloid Leukemia Drug-Resistance Cells

Objective
Acute myeloid leukemia (AML) is a highly heterogeneous hematological malignancy, drug resistance and relapse are key factors in the failure of leukemia treatment. More and more studies have shown that circRNA and LncRNA play an important role in the development of tumors, but it is remains unclear in the mechanism of AML resistance.


Materials and Methods
AML resistant cell line HL-60/ADM (Adriamycin, ADM) was constructed, and circRNA, LncRNA and mRNA expression profiles were screened follow by high-throughput sequencing, then bioinformatics analysis was carried out, circRNA-miRNA ceRNA network was constructed, and confirmed using qRT-PCR analysis.


Results
A total of 1824 circRNA, 2414 LncRNA, and 5346 mRNA were screened as differentially expressed genes, GO and KEGG mainly involved that Protein domain specific binding, transforming growth factor-β (TGF-β) receptor and cellular metabolism. mTOR signaling pathway, MAPK signaling pathway, RAP1 signaling pathway and Akt signaling pathway were closely related to multidrug.


Conclusion
Our study provides a systematic outlook on the potential function of ncRNA in the molecular mechanisms of AML resistant cells. Hsa-circ-0000978 and hsa-circ-0000483 might be served as potential prognostic biomarkers and therapeutic targets of AML resistance.


Introduction
Acute myeloid leukemia (AML) is a highly heterogeneous hematological malignancy. Although its treatment has made significant progress, the prognosis is still unsatisfactory. Recurrence and drug resistance are the main factors [1]. At present, there are many studies on the molecular mechanisms of AML resistance [2,3]. However, with the development of bioinformatics, the epigenetic mechanism in the pathogenesis of AML still remains unclear.
Among human transcripts, about 10%-20% are proteinencoding RNA, and the remaining 80%-90% are noncoding RNAs (ncRNAs) [4,5]. Long noncoding RNAs (LncRNAs) are a class of noncoding RNAs that regulate gene expression at the transcriptional or posttranscriptional level [6]. LncRNA plays an important regulatory role in the drug resistance process. Li et al. [7] reported that the LncRNA HOTTIP can promote the development of pancreatic cancer and regulate gemcitabine resistance by regulating HOXA13, while HOTTIP regulates cisplatin resistance in osteosarcoma cells by activating the Wnt/β-catenin pathway [8]. In addition, Qu et al. [9] found that in sunitinib-resistant renal cell carcinoma, when FOX0 and AKT expression decreased, LncRNA increased, knocking out LncRNA and then reversing drug resistance. Endogenous competition between mir-34 and mir-449 promotes the expression of AXL and c-MET in sunitinib-resistant renal cell carcinoma to regulate the drug resistance process, confirming that LncRNA can be used as a target to repair drug resistance. Circular RNAs (circRNAs) are novel noncoding RNAs characterized by a covalently closed structure with nonrandom spiking and RNase degradation resistance [10,11]. CircRNA is present in the cytoplasm and is extremely abundant and highly conserved and stable in the blood [12]. CircRNAs are increasingly found in various diseases and show cell or tissue specificity [13,14,15].
At present, the molecular mechanism of LncRNA and circRNA in resistant AML cells remains unclear. In this study, highthroughput sequencing of HL-60 and HL-60/ADM (adriamycin, ADM) was performed utilizing Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG). A circRNA-miRNA ceRNA network was constructed to provide new therapeutic targets and the theoretical basis for treatment of drug resistance in AML.

Cell Culture and Drug Resistance Induction
HL-60 cells were incubated in basic RPMI 1640 Medium containing 10% FBS, 100 µg/mL streptomycin, and 100 U/mL penicillin at 37 °C and 5% CO 2 under saturated humidity conditions after recovery. The liquid was changed once every 2 days, and 10 6 cells/mL were amplified at a 1:3 ratio. ADM induction was performed in HL-60 cells by combining the concentration gradient increasing method and the impact method (high-dose intermittent induction) as referenced in the literature [16]. The initial induction concentration was 0.1 µg/mL, shock induction was performed for 1 h, and culturing was continued until HL-60 cells grew and proliferated normally at 1 µg/mL ADM. It took 8 months to successfully induce ADM resistance in HL-60 cells.

CCK-8 Assay and Cell IC 50 Values
HL-60 cells were collected and centrifuged at 1000 rpm/min for 5 min, and then they were resuspended in RPMI 1640 Medium and counted. Furthermore, 10 4 cells of cell suspensions of 100 µL were placed in 96-well plates at 37 °C in a 5% CO 2 incubator for culturing for 24 h. ADM was added with differences in concentrations of 10 µL and cells were incubated for 24 h, and then 10 µL of CCK-8 solution was added to each well. Culture plates were further incubated in the incubator for 4 h. OD values were measured and data were collected to calculate IC 50 values, or ADM concentrations required for 50% inhibition in vitro.

Western Blot Detection of the Expression of Drug-Resistant Protein
RIPA buffer was added with phenylmethanesulfonyl fluoride to collect the cells showing logarithmic growth, and proteins were extracted from the cells. Protein concentration was determined using the BCA Protein Assay Kit. Protein samples containing sample buffer were denatured for 5 min in boiling water. SDS-PAGE electrophoresis was performed with 25 µg of sample in each hole with the addition of 5 µL of prestained protein marker. When the bromophenol blue dye ran off the gel layer, the electrophoresis was terminated, and further experiments were performed on a 35 mA transmembrane overnight. The One-Step Western Kit HRP was used according to the manufacturer's instructions. Rabbit P-gp antibody, antibody pretreatment solution, and dilution buffer solution were added, Li M, et al: Expression Profile of Resistant Leukemia Cells mixed, and poured onto the membrane. Immobilon Western Chemiluminescent HRP Substrate was used for color reaction.

RNA Library Construction and Sequencing
High-throughput sequencing service was provided by CloudSeq Biotech (Shanghai, China). Transcriptome high-throughput sequencing and subsequent bioinformatics analyses were also performed by CloudSeq Biotech (Shanghai, China). Briefly, total RNA was used to remove the rRNAs using the Ribo-Zero rRNA Removal Kit (Illumina, USA) according to the manufacturer's instructions. RNA libraries were constructed using rRNAdepleted RNAs with the TruSeq Stranded Total RNA Library Prep Kit (Illumina, USA) according to the manufacturer's instructions. Libraries were checked for quality, and they were quantified using the Bioanalyzer 2100 system (Agilent Technologies, USA). Furthermore, 10 pM libraries were denatured as single-stranded DNA molecules, captured on Illumina flow cells, amplified in situ as clusters, and finally sequenced for 150 cycles on an Illumina HiSeq Sequencer according to the manufacturer's instructions.

Bioinformatics Analysis
For circRNA, high-quality reads were aligned to the reference genome/transcriptome with STAR software (v2.5.1b) and circRNAs were detected and identified with DCC software (v0.4.4). edgeR software (v3.16.5) was used to normalize the data and perform analysis of differentially expressed circRNAs. GO and KEGG analyses were performed for the differentially expressed circRNA-associated genes.
For LncRNA and mRNA, high-quality reads were aligned to the human reference genome (UCSC hg19) with HISAT2 software (v2.0.4). Then, guided by the Ensembl gtf gene annotation file, the Cuffdiff program (v2.2.1, part of Cufflinks software) was used to get the FPKM (fragments per kilobase of exon model per million reads mapped) for the expression profiles of LncRNA and mRNA. Accordingly, fold change and p-values were calculated based on FPKM, and differentially expressed LncRNAs and mRNAs were identified. LncRNA target genes were predicted by locations in relation to nearby genes, and GO and pathway analyses were performed on these target genes.

Validation of Differentially Expressed circRNAs
Total RNA was extracted by TRIzol (Invitrogen Life Technologies, Shanghai, China) to synthesize cDNA via reverse transcription. Quantitative real-time PCR was performed on the ViiA 7 Real-Time PCR System (Applied Biosystems) using the qPCR SYBR Green Master Mix (CloudSeq). Primer design was as follows: chr13: The target RNA and internal parameters of each sample were subjected to real-time PCR, which was repeated three times. The data were analyzed by the 2 -ΔΔC T method.

Statistical Analysis
All of the experimental data are presented as mean ± standard deviation (SD), and the t-test was used for comparisons between the two groups. Values of p<0.05 were considered to be significantly different. GraphPad Prism 5 software was used for statistical analysis.

Construction of HL-60/ADM Drug-Resistant Cell Lines
Adriamycin resistance was induced in HL-60 cells via a combination of the concentration gradient method and the impact method (high-dose intermittent induction method) according to the literature [16], and it took 8 months to obtain HL-60/ADM resistant cell lines. IC 50 detection is shown in Figure 1A. The expression of the drug resistance protein P-gp in HL-60/ADM resistant cells was significantly higher than that in HL-60 cells, with significance at p<0.01 ( Figures 1B and 1C).

Expression of circRNA, LncRNA, and mRNA in Resistant Cells
To analyze the gene expression of resistant AML cells, we performed high-throughput sequencing of circRNA, LncRNA, and mRNA, and we screened differentially expressed genes. Hierarchical cluster analysis showed that the expression patterns  Figure 2A; LncRNA, Figure 2D; mRNA, Figure  2G). Scatter plots were used to evaluate circRNA between drug-resistant and drug-sensitive cells. LncRNA and mRNA signal values were normalized to log2 values for visualization of expression differences (Figures 2B, 2E, and 2H, respectively). Volcanic maps were constructed based on fold change (FC≥1.2) and p-value (<0.05), and volcano maps of the differentially expressed genes between these two different conditions are provided in Figures 2C, 2F, and 2I. The general characteristics of RNA include RNA type, length, and localization distribution, as shown in Figures 2J-2L. The results showed that resistant AML cells differentially expressed 1824 circRNAs, 2414 LncRNAs, and 5346 mRNAs.

Functional Analysis of Differentially Expressed circRNA, LncRNA, and mRNA
To explore the underlying genomics mechanisms involved in the developmental disorders of AML tumorigenesis, GO and KEGG pathway enrichment analyses of differentially expressed genes were used to evaluate candidate RNA functions. The GO terms with the highest enrichment scores for upregulated circRNA targeting were ribonucleoside triphosphate catabolic process and purine ribonucleoside triphosphate catabolic process (BP), intracellular (CC), and adenyl ribonucleotide binding (MF); for LncRNA, the GO terms involved macromolecule modification (BP), cornified envelope (CC), and protein homodimerization activity (MF); and anatomical structure morphogenesis (BP), cytoplasm (CC), and protein binding (MF) belonged to the GO analysis of mRNA. Furthermore, KEGG pathway analysis was performed to predict potential module functions. The KEGG analysis results were as follows: for circRNA, the B cell receptor signaling pathway (hsa04662), T cell receptor signaling pathway (hsa04660), MAPK signaling pathway (hsa04010), and mTOR signaling pathway (hsa04150); for LncRNA, signaling pathways regulating pluripotency of stem cells (hsa04550); and for mRNA, the Wnt signaling pathway (hsa04310), Rap1 signaling pathway (hsa04015), p53 signaling pathway (hsa04115), and VEGF signaling pathway (hsa04370). These are closely related to cancer progression and are significantly enriched in AML (Figure 3).

Construction of circRNA-miRNA ceRNA Network
To fully understand the underlying mechanisms of circRNA and AML development, based on differentially expressed circRNA data, we used a database to predict target miRNAs interacting with circRNA, and Cytoscape was used to construct a circRNA-targeted miRNA gene network map (Figure 4). For a particular miRNA, circRNA has many targets, and the network map illustrates the first five predicted miRNA targets that differentially express circRNA.

Discussion
At present, chemotherapy is still one of the main treatments for leukemia. However, multidrug resistance and treatment are key factors in the failure of leukemia treatment. Several factors are involved in the mechanism of leukemia resistance, including ABC transporter-mediated multidrug resistance [17], DNA repair abnormalities [18], variations in the bone marrow microenvironment [19], and abnormal expression of noncoding RNAs including circRNA, miRNA, and LncRNA [20]. Undeniably, ncRNAs have opened up new prospects for AML diagnosis, prognosis, and treatment. Indeed, the expression of specific ncRNAs such as circRNAs and LncRNAs could assist clinicians in classifying subtypes, evaluating prognosis, and predicting the response to drug treatment in AML. Garzon et al. [21] evaluated the associations of LncRNA expression with clinical characteristics, gene mutations, and outcomes and constructed an LncRNA score including 48 LncRNAs for independently predicting outcome prognosis, confirming that LncRNAs can assist in predicting clinical outcomes in older patients with CN-AML. Moreover, Li and Sun [22] reported that SNHG5 overexpression was frequently observed in AML patients with advanced FAB classification and unfavorable cytogenetics. Furthermore, a higher SNHG5 expression level was also associated with shorter overall survival. However, comprehensive analyses of the profiles of differentially expressed circRNAs, LncRNAs, and mRNAs in resistant AML cells have not been studied. Thus, we explored the expression profiles and predicted the potential functions of circRNAs, LncRNAs, and mRNAs in resistant AML cells by utilizing RNA high-throughput sequencing and bioinformatics analysis.
The numbers of differentially expressed genes of circRNAs, LncRNAs, and mRNA in resistant AML cells are 1824, 2414, and 5346, respectively. GO and KEGG pathway analyses of differentially expressed LncRNAs mainly revealed protein domain specific binding and protein dimerization activity. This provides a basis for important contributions to the development and resistance of leukemia. Although the understanding of the nature and function of circRNAs is still limited, it is undeniable that circRNA has always been a research hotspot in the field of ncRNAs, which particularly regulate miRNA-targeted gene expression as ceRNA molecules [23,24]. GO and KEGG pathway analyses predicted that these differentially expressed circRNA functions were related to tumor development, drug-resistant regulation, and metabolism-related pathways.
GO analysis mainly revealed involvement with PI3K activity, transforming growth factor-β (TGF-β) receptor, and cellular metabolism. Sui et al. [25] reported that activation of the PI3K/Akt/NF-κB pathway promotes P-gp expression, and the inhibition of this pathway reverses P-gp-mediated multidrug resistance. Zhou et al. explored the effect of the PI3K-specific inhibitor ZSTK474 on K562/A02 cells and their results showed that ZSTK474 reversed the resistance of K562/A02 cells to ADM and imatinib by downregulating P-gp expression; accordingly, the target of ZSTK474 for CML treatment is PI3K [26]. It can be  seen that PI3K activity plays an important role in the regulation of drug resistance in leukemia. Results of KEGG analysis indicated that the upregulated circRNA in drug-resistant cells was mainly related to the mTOR signaling pathway, MAPK signaling pathway, and Akt signaling pathway. The PI3K/Akt signaling pathway maintains a close relationship between tumor cell multidrug resistance and P-gp. Studies have shown that activation of the PI3K/Akt signaling pathway increases drug efflux via the ATPbinding cassette (ABC) transporter [27], while the blocking of PI3K/Akt signaling pathways leads to downregulation of P-gp and MRP1 expression, restoring sensitivity to chemotherapeutic drugs [28].
According to the ceRNA theory, a circRNA-miRNA regulation network is present in cases of resistant AML. Moreover, differentially expressed circRNAs-miRNA interactions were predicted and the potential molecular mechanisms were further explored. Among these predicted potential target miRNAs, hsa-miR-24-2-59 is reported to be upregulated during hematopoietic cell terminal differentiation, suppressing MYC expression [29]. Hsa-miR-181b-5p may play a prominent role in pituitary adenoma as an effective biomarker and therapeutic target [30]. However, there are some limitations to this study, such as a small sample size and the in vitro research being conducted only on HL-60 cells. The next step of this work will be to verify the expression of LncRNAs and circRNAs in AML patients and to study the mechanisms of LncRNAs and circRNAs in the development of resistant AML.

Conclusion
Even though only resistance to adriamycin was assessed in this study, the results suggest that the expression changes of circRNA/LncRNA regulate the cell resistance of AML, providing a new theoretical basis for the further understanding of multidrug resistance mechanisms and targeted therapies in AML.

Ethics Committee Approval:
No human studies were involved in this work.

Informed Consent:
This study did not involve animal and human ethics.