Circular RNA regulatory network reveals cell–cell crosstalk in acute myeloid leukemia extramedullary infiltration

Background Acute myeloid leukemia can develop as myoblasts infiltrate into organs and tissues anywhere other than the bone marrow, which called extramedullary infiltration (EMI), indicating a poor prognosis. Circular RNAs (circRNAs) are a novel class of non-coding RNAs that feature covalently closed continuous loops, suggesting their potential as micro RNA (miRNA) “sponges” that can participate in biological processes and pathogenesis. However, investigations on circRNAs in EMI were conducted rarely. In this study, the overall alterations of circRNAs and their regulatory network between EMI and non-EMI AML were delineated. Methods CircRNA and whole genome microarrays derived from EMI and non-EMI AML bone marrow mononuclear cells were carried out. Functional analysis was performed via Gene Ontology and KEGG test methods. The speculated functional roles of circRNAs were based on mRNAs and predicted miRNAs that played intermediate roles. Integrated bioinformatic analysis was conducted to further characterize the circRNA/miRNA/mRNA regulatory network and identify the functions of distinct circRNAs. The Cancer Genome Atlas (TCGA) data were acquired to evaluate the poor prognosis of distinct target genes of circRNAs. Reverse transcription-quantitative polymerase chain reaction was conducted to identify the expression of has_circRNA_0004520. Connectivity map (CMap) analysis was further performed to predict potential therapeutic agents for EMI. Results 253 circRNAs and 663 genes were upregulated and 259 circRNAs and 838 genes were downregulated in EMI compared to non-EMI AML samples. GO pathways were enriched in progress including cell adhesion (GO:0030155; GO:0007155), migration (GO:0016477; GO:0030334), signal transduction (GO:0009966; GO:0007165) and cell–cell communication. Overlapping circRNAs envolved in pathways related to regulate cell–cell crosstalk, 17 circRNAs were chosen based on their putative roles. 7 target genes of 17 circRNAs (LRRK1, PLXNB2, OLFML2A, LYPD5, APOL3, ZNF511, and ASB2) indicated a poor prognosis, while overexpression of PAPLN and NRXN3 indicated a better one based on data from TCGA. LY-294002, trichostatin A and SB-202190 were identified as therapeutic candidates for EMI by the CMap analysis. Conclusion Taken together, this study reveals the overall alterations of circRNA and mRNA involved in EMI and suggests potential circRNAs may act as biomarkers and targets for early diagnosis and treatment of EMI. Electronic supplementary material The online version of this article (10.1186/s12967-018-1726-x) contains supplementary material, which is available to authorized users.


Background
Acute myeloid leukemia (AML) is one of the most common hematologic malignancies that aggressive myeloid blasts derived from bone marrow (BM) have disseminated into the peripheral blood (PB), leading to the accumulation of myeloblasts in both the BM and PB [1]. Remarkably, BM-derived blasts may also represent extramedullary infiltration (EMI) in areas other than the BM, including cutis injury, myeloid sarcoma, lymphoid infiltration and central nervous system (CNS) infiltration [2]. EMI may be present at the initial diagnosis of AML, during relapse, or even when patients achieve a complete remission in the BM. EMI predicts a poor outcome of AML due to the destruction of important organs (e.g., in the CNS) and indicates that refractory and relapse leukemia are likely to develop [3]. EMI is often misdiagnosed as inflammation, lymphoma or melanoma [4].
Current studies on the molecular and cellular manifestations of EMI center on cytogenetic backgrounds [1,5], gene mutations [6], cluster of differentiation (CD) markers and significant non-coding RNAs [7]. Early diagnosis before the appearance of clinical manifestations and interventional treatment may delay disease progression. Furthermore, clarifying the alterations of biological function underlying EMI may provide support to overcome the disease.
Circular RNAs (circRNAs) were considered outcomes of a "splicing error" for decades [8]. However, with the development of high throughput analysis, the functions of circRNAs are being identified. Their specific covalently closed loop structure facilitates their resistance to RNase R, which leads to a half-life of approximately 24 h, compared to the 4-6 h half-life of mRNA at 4 °C [9]. Additionally, circRNAs exhibit tissue and temporal expression specificity and participate in multiple biological processes and human diseases. Emerging evidence reveals that circRNAs can act as miRNA "sponges" [10], interact with RNA binding proteins [11], and regulate parental gene expression, and some can even be translated [12]. Deregulated circRNAs exhibit close relationships with hematological malignancies, especially in AML, where circRNAs are derived from the NMP1 gene [13], and fusion-circRNA (f-circRNA) is derived from the fusion gene MLL-AF9 [14]. However, neither circRNA analysis nor microarray analysis of EMI in AML have been previously conducted.
This study aims to provide an overview of circRNA profiles in EMI to identify biomarkers to aid early diagnosis and treatment. A circRNA/miRNA/mRNA regulatory network was constructed and subjected to a pathway analysis, which revealed that cell-cell interactions are instrumental in EMI; a connectivity map (CMap) analysis was performed based on mRNA dysregulated in the whole regulatory network that indicates three potential therapeutic agents (LY-294002, trichostatin A and SB-202190) for EMI; these results improve our understanding of the cellular progression and pathogenesis of EMI.

Samples
Twelve samples from 4 matched AML patients with and without EMI and healthy volunteers were used in our study. The study was approved by the Ethics Committee of the First Affiliated Hospital of Harbin Medical University. Informed consent was obtained in accordance with the Declaration of Helsinki. Patients with EMI found at diagnosis or at relapse were included, and the clinical features included leukemic cutis, imaging evidence of sarcoma, lymphadenectasis, or manifestations of other types of organ destruction, together with morphological characteristics of myoblast infiltration. The diagnosis was established using the WHO diagnostic criteria [15]. Samples of AML patients without EMI were collected during diagnostic biopsy. The detailed patient information is described in Table 1.

Total RNA extraction
The total RNA was isolated from the samples of both groups using TRIzol (Invitrogen, CA, USA) and a miRNeasy mini kit (QIAGEN, CA, USA) according to the manufacturer's instructions. The RNA quality and quantity were accurately measured using a NanoDrop spectrophotometer (ND-1000, Nanodrop Technologies, DE, USA). The RNA integrity was assessed by gel electrophoresis.

Whole genome microarray
The Human Whole Genome Oligo Microarrays (onecolor) (Agilent, CA, USA), which covers 44,000 probes detecting approximately 27,958 genes obtained from NCBI and other authoritative data sources. Total RNA was then subjected to reverse transcription and converted to double-strand cDNA (ds-cDNA) using an Invitrogen SuperScript ds-cDNA synthesis kit with oligo dT primers. The ds-cDNA was then purified and labeled with Cy3 using the NimbleGen One-Color DNA labeling kit in accordance with the NimbleGen Gene Expression Analysis Procedure (NimbleGen Systems, Inc., Madison, WI, USA). Thus, the Cy3 labeled-ds-cDNA was hybridized with the microarrays, which were then washed with the Wash Buffer kit (NimbleGen Systems, Inc., Madison, WI, USA).

CircRNA microarray
CircRNAs from the 12 samples were detected by the Arraystar Human circRNA Array V2 (8 × 15 K, Arraystar, Rockville, MD, USA) including 15,000 probes detecting 13,617 circRNAs ever obtained from circbase and other databases. Total RNA from each sample was quantified using the NanoDrop ND-1000 (Thermo Scientific, DE, USA). The sample preparation and microarray hybridization were performed based on the Arraystar's standard protocols. Briefly, total RNAs were digested with Rnase R (Epicentre, WI, USA) to remove linear RNAs and enrich circular RNAs. Then, the enriched circular RNAs were amplified and transcribed into fluorescent cRNA utilizing a random priming method (Arraystar Super RNA Labeling Kit; Arraystar, Rockville, MD, USA). The labeled cRNAs were hybridized onto the Arraystar Human cir-cRNA Array V2. After having washed the slides, the arrays were scanned by the Agilent Scanner G2505C.

Microarray data analysis
Agilent feature extraction software (version 11.0.1.1) was used to analyze acquired both genome expression and circRNA array images. Quantile normalization and subsequent data processing were performed using the R software limma package. Differentially expressed genes and circRNAs with statistical significance between two groups were identified through Volcano Plot filtering. Differentially expressed circRNAs between two samples were identified through Fold Change filtering.
Hierarchical Clustering was performed to show the distinguishable circRNAs expression pattern among samples.

GO and KEGG enriched pathway analysis
To explore the underlying pathways and biological processes in EMI, we conducted GO and KEGG pathway analyses based on all differentially expressed mRNAs using the Database for Annotation Visualization and Integrated Discovery (DAVID version 6.8; https ://david .ncifc rf.gov/) [18] and the latest KEGG database (http:// www.genom e.jp/kegg/) [19].

Target miRNA and mRNA prediction and regulatory network construction
MiRNAs targeted by circRNAs were predicted using TagetScan [20] and miRanda [21] (Threshold: targetscan ≤ 0.05, sitecoverage ≥ 3). miRanda was used for MRE identification and evaluation, and TargetScan was used to supplement the MRE evaluation. Putative miRNAs were listed based on competitive binding ability, and the top 5 miRNAs were selected for further mRNA predictions. mRNA was predicted by TargetScan, microT and Tarbase on the DIANA-miRPath v.3 platform [22] (Threshold: TargetScan: − 0.4, microT:0.8), and only the target mRNAs presented more than 2 times among the 3 databases were considered target genes of the given miRNAs. The targeted mRNAs were then overlapped with our whole genome microarray data.

Validation the expression of hsa_circRNA_0004520 by reverse transcription-quantitative polymerase chain reaction
RT-PCR was conducted on ViiA 7 Real-time PCR System (Applied Biosystems, Foster City, CA, USA) based on the above 4 pairs of EMI and non-EMI AML samples. Hsa_circRNA_0004520 was amplified using the specific primers shown in Table 2. The PCR reaction system with 10 μl volume consisted of 2 μl cDNA, 0.5 μl of primers (10 μM), 5 μl of 2× Master Mix and 2 μl of H 2 O. The amplification conditions of qRT-PCR reaction was carried out at 95 °Cfor 10 min, followed by 40 cycles at 95 °C for 10 s, 60 °C for 60 s to collect fluorescence, and finally followed by the melting program at 95 °C for 10 s, 60 °C for 60 s,95 °C for 15 s, and 60 °C to 99 °C slowly.

TCGA data accession and prognostic analysis
To further validate that EMI-related circRNAs and targeted genes have the potential to predict a poor prognosis in AML patients, clinical and RNA-seq data were acquired from TCGA database via LinkedOmics (http:// www.linke domic s.org/) [24]. Thirty-five genes targeted by 17 candidate circRNAs were further analyzed according to the protocol, and the parameters were set as follows: Cancer Type: acute myeloid leukemia (AML); RNA-seq based on HiSeq RNA conducted on 01/28/2016 was selected; the dataset for an ethnicity of not Hispanic or Latino was chosen; the target dataset was a clinical dataset conducted on 01/28/2016. Statistical analysis was performed by a non-parametric test, with a value of p ≤ 0.05 considered statistically significant.

CMap analysis
Differentially expressed genes involved in the circRNA/ miRNA/mRNA regulatory network were listed into upand down-regulated tags and then were uploaded into the CMap web tool (https ://porta ls.broad insti tute.org).
Matches between the signatures of interest and chemicals from CMap were assessed through a connectivity score from −1 to 1 that a positive score indicates the stimulative effect of compound on the RNAs; while a negative score implicates a inhibiting effect of a compound [25].

CircRNA and differential gene expression in EMI and non-EMI AML BM samples
We examined 4 matched samples each of AML patients with or without EMI and healthy adults using cir-cRNA (Arraystar Human circRNA Array V2, 8 × 15 K, Arraystar, Rockville, MD, USA) and whole genome microarrays (Agilent, CA, USA) at the Department of Hematology of The First Affiliated Hospital of Harbin Medical University. CircRNAs and mRNAs with fold changes ≥ 2.0 and p-values ≤ 0.05 were considered significantly differentially expressed. Two hundred fifty-three circRNAs and 663 genes were upregulated, while 259 circRNAs and 838 genes were downregulated in AML patients with EMI compared to those without EMI (Additional file 1: Tables S1 and Additional file 2: Table S2). Of the differentially expressed circRNAs detected, 84.77% were covered in the exon of the genome (Fig. 1a). Hierarchical clustering was performed between of the AML with and without EMI groups to demonstrate the circRNA expression patterns (Fig. 1b, c). Furthermore, clustering overlapping was conducted to identify the circRNAs likely responsible for EMI. There were 353 upregulated circRNAs in the EMI versus non-EMI AML groups, which further overlapped with circRNAs upregulated in the AML with EMI versus the healthy groups. Finally, 82 upregulated and 27 downregulated circRNAs were considered to have closer relationships with EMI ( Fig. 2a-d).

Pathway enrichment of differentially expressed genes in EMI based on GO and KEGG analysis
To clarify the biological processes related to EMI, Gene Ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomics (KEGG) pathway enrichment were performed on all 663 upregulated and 838 downregulated genes aberrantly expressed in the EMI group.
The GO analysis results revealed that these upregulated genes participated in vital biological processes including cell adhesion (GO:0030155; GO:0007155), migration (GO:0016477; GO:0030334), signal transduction (GO:0009966; GO:0007165) and cell-cell communication (GO:0010646; 0007154) ( Fig. 3a, b, Additional file 3: Table S3 and Additional file 4:  Table S4), which is consistent with current reports on the relationship between EMI and the microenvironment [26]. KEGG enriched pathway analysis indicated that several significative pathways might be aberrantly activated in response to EMI, including the RAP1 signaling pathway (hsa04015), Cell adhesion molecules (hsa04514) and other immunology related pathways (Additional file 5: Table S5 and Additional file 6: Table S6).

Construction and visualization of the circRNA/miRNA/ mRNA regulatory network
Given the potential regulatory roles of circRNAs on miRNA "sponges", two algorithms that focused on miRNA response elements (MREs) between circRNA and miRNA (miRanda and TargetScan) were used to predict miRNAs related to EMI. To further associate circRNAs with biological processes, we delineated the regulation of circRNAs and expressed genes indirectly, i.e., the circRNA/miRNA/mRNA network. We included Finally, we overlapped the predicted target mRNAs with our whole genome microarray profiles. A regulatory network involving circRNAs/miRNAs/genes was constructed using Cytoscape software (version 3.6.0; http:// www.cytos cape.org) Additional file 7: Figure S1), and the prediction of functional roles of circRNA was based on the downstream genes.

Seventeen circRNAs were associated with vital processes in EMI
To explore the candidate circRNAs that may play vital roles in EMI of AML, 82 upregulated and 27 downregulated circRNAs were selected for further analysis after overlapping. We then eliminated candidate circRNAs based on their associated biological processes. To determine which circRNAs were involved in the regulation of cell-cell interactions, genes related to these processes were retained. Thirty circRNAs were responsible for cell migration, 51 were responsible for cell-cell communication, 53 were responsible for signal transduction, and 27 were responsible for cell adhesion. CircRNAs involved in the processes mentioned above were further overlapped to identify the "kernel" circRNA. After excluding "irrelevant" genes, miRNAs and circRNAs, 17 circRNAs likely to play important roles in EMI remained, and the network of circRNAs related to EMI was diagrammed ( Table 3, Fig. 4). And the structure of 9 circRNAs (8 not available) were described in Fig. 5 based on the data from Cancer Specific CircRNA (CSCD, http://gb.whu.edu.cn/ CSCD/) [27].

Validation of the expression of hsa_circRNA_0004520
RT-PCR was conducted in 4 pairs of EMI and non-EMI AML BMMCs to validate the expression of key circR-NAs. Hsa_circRNA_0004520 was chosen for our interest on its target genes: PLXNB2 and VEGFA. The expression of hsa_circRNA_0004520 was upregulated by 4.38-fold (p = 0.144) (Fig. 6).

Validation of the prognostic value of 17 circRNAs that target genes involved in AML
To validate the prognostic value of circRNAs related to EMI, we obtained access to the Cancer Genome Atlas (TCGA) database to assess the relationship between the target genes of the 17 circRNAs and the prognosis in AML patients because EMI is associated with a poorer prognosis than non-EMI AML. A non-parametric test was applied to clinical and RNA-seq data, and a value of p ≤ 0.05 was considered a statistically significant difference. AML with overexpression of 7 genes (LRRK1, PLXNB2, OLFML2A, LYPD5, APOL3, ZNF511, and ASB2) indicated a poor prognosis, while overexpression of PAPLN and NRXN3 indicated a better prognosis (Fig. 7). The data obtained from TCGA suggested the potential roles of circRNA and differentially expressed genes in the  genes (a, b). The Y-axis shows the names of process that were significantly enriched (a, b). The size of each node represents the number of enriched differential genes, and the p-values are indicated by color changes from red to purple (b). Downregulated genes were not associated with specific EMI-related processes Page 8 of 15 Lv et al. J Transl Med (2018) 16:361 AML prognosis and further indicated that they might play crucial roles in EMI. However, this finding should be confirmed in further investigations.

Identification of three bioactive compounds for the treatment of EMI based on CMap analysis
All differentially expressed mRNA in the circRNA/ miRNA/mRNA regulatory network were loaded into the CMap web tool. Following the signature query, three compounds with potential anti-tumor effect, LY-294002, trichostatin A (TSA) and SB-202190, i.e. with the highest negative enrichment score were determined as the potential therapeutic agents for EMI ( Table 4). The 3D chemical structures of the three compounds obtained from NCBI are presented in Fig. 8.

Discussion
In this study, we delineated a circRNA/miRNA/gene regulatory network related to EMI using microarray and bioinformatic tools for the first time. The altered genes are downstream in the regulatory network, and the functions were predicted via bioinformatics methods. MiRNAs were used to link mRNAs and circRNAs to predict their function in EMI pathogenesis. CMap analysis revealed LY-294002, trichostatin A (TSA) and SB-202190 to be treatment candidates for EMI based on dysregulated mRNAs. We explored potential biological processes and pathogenesis and predicted conceivable roles for cir-cRNA throughout the EMI network.
EMI is relatively common in AML, with myeloid sarcoma appearing in nearly 1.4-9% of AML patients [28], and leukemic cutis present in approximately 15% of patients. It has been estimated that EMI occurs in 18-25% of childhood AML cases [29]. EMI exhibits a close relationship with relapse and refractory AML. Strikingly, after patients have undergone allo-hematopoietic stem cell transplantation (HSCT), EMI may still appear or even worsen [30].
Previous studies have investigated cytogenetic factors in EMI, such as alterations in chromosome 8 alteration, FLT3-ITD and NPM1 mutations. As the application of high throughput technology has increased, the roles of circRNAs are being clarified. Based on the circRNA/ miRNA/mRNA regulatory network, the functions of cir-cRNA in EMI were assessed, which primarily involved cell adhesion, migration, signal transduction and cellcell communication. Previous reports on EMI have also considered cell-cell communication to play an important role in EMI, with leukemic cells subverting normal interactions and converting a normal niche into a neoplastic niche [31]. circRNAs responsible for multiple processes were thought to play an even more important role. After overlapping of four processes, 17 circRNAs in our study were identified. Due to the interesting targets (PLXNB2 and VEGFA, which may take part in angiogenesis pathways) of has_circRNA-0004520, RT-PCR was performed to validate its expression, and the function needs to be identified in our further investigation.  Fig. 4 circRNAs involved in vital processes in EMI. a circRNAs related to cell-cell crosstalk were further overlapped to identify circRNAs that act as "kernels" in EMI. Seventeen circRNAs had the potential to regulate all 4 processes including cell migration, cell-cell communication, signal transduction, and cell adhesion. b The 17 "kernel" circRNAs and their target miRNAs and genes were finally delineated into a regulatory network. Round nodes represent circRNAs, V nodes represent differentially expressed genes, and red represents upregulated genes. Blue triangles represent predicted miRNAs. Lines indicate a regulatory relationship manifestation arises. Interestingly, 2 genes (PALPN and NRXN3) seemed to exhibit a protective effect against EMI thus prolong the overall survival, which needs to be investigated. We selected circRNA as a novel biomarker for early diagnosis and treatment, thus hinder the progression of EMI and further efforts to overcome this process. Interestingly, a CMap analysis based on mRNAs dysregulated in EMI was applied to explore practicable drugs for the treatment of EMI. Drugs on CMap were all available for human use and listed in the Food and Drug Administration. Data on CMap are based on genome-wide transcriptional expression profiles and accurately show the promoting and inhibiting effects on gene expression in transcriptional level. Three chemicals were identified to be options for EMI based on CMap results. LY-294002 is considered to be a conventional specific inhibitor of phosphatidylinositol 3-kinase (PI3K) [32]. PI3K/Akt pathways play important roles in multiple oncogenetic processes. LY-294002 has been applied for solid tumors for decades. While in hematological malignancies, there are not any PI3K inhibitors for clinical application. Recent studies on the resistance of chronic myeloid leukemia or to be a sensitizer for chemotherapy reveal that LY-294002 might work in refractory and relapsed leukemia [33,34], such as EMI. SB-202190 is inhibitor for p38 MAPK signaling pathway [35]. P38 MAPK pathways may even play opposite effect in diverse tumors. Compared to LY-294002, SB-202190 has more applications in AML. SB-202190 can diminish doxorubicin-induced resistance in K562 cells [36], but promote THP-1 and MV4:11 cells also [37], and the exact effect of SB-202190 in EMI needs to be clarified later. As an inhibitor of histone deacetylases (HDAC), trichostatin A (TSA) has been reported for a potent inducer of tumor cell growth arrest, and apoptosis in many diverse transformed cells and tumor-bearing animals by regulating the expression of tumor suppressor genes [38]. HDAC inhibitors are gradually in use for lymphoma and multiple myeloma in this decade [39,40]. The use of HDAC inhibitor in our EMI patients also appears a modest curative effect (data not shown). However, there are limitations in this study. First, EMI is a very complex disease. The mechanisms of different subtypes and organs may be different. A larger and more accurate group of samples is needed to clarify these differences. In addition, only parts of circRNAs function as miRNA "sponges". Some circRNAs may also regulate parent gene expression and participate in EMI. Additionally, classical post-transcriptional regulation of miRNAs allows them to bind to the 3′ UTR region of mRNAs and thus degrade them, while the targeted relationship still needs to be validated. Subsequent investigations will take this finding into account. Notably, our study focused on RNA levels, and further research on protein levels in cell lines and in vivo models is needed.

Conclusions
Our research is the first to use microarrays of circR-NAs and whole genome expression together with bioinformatic tools to identify a circRNA/miRNA/gene regulatory network of EMI in AML patients. Seventeen circRNAs may act as key regulators of cell-cell crosstalk in EMI, and most of their target genes predict a poor prognosis. LY-294002, trichostatin A (TSA) and SB-202190 are potential treatment candidates for EMI. The unique structure and function makes circRNA an