circRNA-TBC1D4, circRNA-NAALAD2 and circRNA-TGFBR3: Selected Key circRNAs in Neuroblastoma and Their Associations with Clinical Features

Introduction

Neuroblastoma (NB) is the most common extracranial malignant solid tumor in children and accounts for about 15% of cancer-associated mortality in childhood.¹ The outcome remains poor due to distant metastasis at initial diagnosis.² The prognosis for NB patients in advanced stage has not improved. Thus, exploiting effective drug targets and better biomarkers for advanced NB is of great necessity.

The pathogenic mechanism of NB remains unclear. A recent study showed that non-coding RNAs (ncRNAs) act as key regulators of physiological processes in the occurrence and development of diseases.³ Particularly relevant in cancer, ncRNAs have been identified as tumor suppressors and oncogenic drivers in diverse cancer types. NcRNAs are functional RNA species that include microRNAs (miRNAs), a class of short ncRNA (∼21–25 nucleotides), long ncRNAs (lncRNAs) consisting of more than 200 nucleotides, and circular RNAs (circRNAs) characterized by covalently closed-loop structures. They regulate gene expression post-transcriptionally and participate in diverse physiological programs such as apoptosis, autophagy, proliferation and angiogenesis.^4,5 We identified hsa-miR-21 (miR-21) as a target gene of tumorigenesis in NB⁶ and found that LINC01296 was the most upregulated lncRNA by lncRNA microarray in advanced-stage NB tumors which was associated with poor prognosis.⁷

Recently, circRNAs have been reported to play a critical role in the pathogenesis of several cancers, such as ovarian cancer, hepatoblastoma and bladder cancer.^8–10 The aim of this study was to investigate and validate potential therapeutic targets of circRNAs in NB.

Methods

Patients and Specimens

Tumor samples of NB and adjacent normal adrenal medulla were collected from 20 patients undergoing surgery from September 2014 to January 2019 in Children’s Hospital of Fudan University (Shanghai, China). This study was approved by the Institute Research Ethics Committee of the hospital and written informed consent was obtained from every patient. All patients were pathologically diagnosed with primary NB without preoperative chemotherapy and the samples were immediately stored at −80°C until use. Among the samples, five pairs were randomly selected for high-throughput RNA sequencing (RNA-seq) to identify differentially expressed circRNAs. The clinical information for five NB samples used is described in Table 1. We selected several of the highly upregulated and downregulated circRNAs, as well as miR-21-related circRNAs, for further analysis, and all 20 paired specimens were used for further expression validation by real-time quantitative reverse transcription polymerase chain reaction (real-time qRT-PCR). The clinical features of patients, including age, INSS stage, risk group, histologic classification, MYCN status, concentration of lactate dehydrogenase (LDH), neuron-specific enolase (NSE) and serum ferritin (SF), were also collected for further statistical analysis.

Table 1 Clinical Information for Five Paired Neuroblastoma Samples Used for circRNA Sequencing

circRNA Sequencing

Total RNA was extracted using TRIzol reagent (Life Technologies, Carlsbad, CA, USA). After using the ribo-zero kit to digest ribosomal RNA and RNase R to remove linear RNA, the enriched circRNAs were reverse-transcribed into cDNA with random hexamer primers, and then used for library construction and RNA-seq via the Illumina sequencer (Oebiotech, Shanghai, China). The raw data were processed by Trimmomatic software (http://www.usadellab.org/cms/index.php?page=trimmomatic)¹¹ to produce high quality clean data which were then aligned to circBase¹² to predict the circRNA sequences by CIRI software (https://sourceforge.net/projects/ciri/).¹³ Spliced reads per millions, which is the ratio of number of back-spliced junction reads to number of total reads, was utilized to quantitate the expression level of each circRNA in each sample. Correlation tests, principal component analysis (PCA) and hierarchical clustering analysis were used to show the degree of similarity of circRNA expression patterns among different specimens.

Bioinformatics Analysis

Junction read counts of each circRNA in every sample were standardized to base mean to estimate the expression level via DESeq software (in R, enter citation [“DESeq”;  http://www.R-project.org]).^14,15 Two criteria were applied to estimate the discrepancy of circRNA expression level between any two groups of samples: FoldChange, the multiple of change of the specific circRNA expression level in tumor compared with normal tissue, and P value, generated via negative binomial distribution test. The screening criteria set for up- and downregulated circRNAs were fold change >2.0 and a P value ≤ 0.05. Hierarchical clustering analysis was performed to display the correlation of dysregulated circRNA expression profiles among different samples. Gene Ontology (GO) analysis (http://geneontology.org/) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis (http://www.genome.jp/kegg/) were performed for functional annotation of the host genes of differentially expressed circRNAs.

Validation of Dysregulated circRNAs by Real-Time qRT-PCR

Total RNA was extracted using TRIzol reagent and RNA concentration was measured by NanoDrop photospectrometer (Thermo Scientific, Waltham, MA, USA). Total RNA (1000 ng) from each sample was added into the reverse transcription system to synthesize cDNA using the PrimeScript RT reagent Kit (TaKaRa, Beijing, China). Specific primers spanning the back-splice junction site of circRNAs were designed and then synthesized (Table S1). The qRT-PCR analyses were performed using a LightCycler^® 480 System (Roche, Rotkreuz, CH) with TB GREEN Premix Ex Taq II (TaKaRa). Relative RNA expression level was calculated by the 2^−ΔΔCt method with GAPDH serving as internal control.

Cell Culture and circRNA-TBC1D4 Overexpression

SH-SY5Y cells were purchased from the Cell Bank of the Chinese Academy of Science (Shanghai, China). The neuroblastoma cell line SH-SY5Y was cultured in Dulbecco’s Modified Eagle’s Medium/Nutrient Mixture F-12 (DMEM/F12) (Meilunbio, Dalian, China) supplemented with 10% fetal bovine serum (FBS) (Gibco) in a humidified incubator with 5% CO₂ at 37 °C. Lentiviruses LV18-circRNA-TBC1D4 and LV18NC were purchased from GenePharma (Shanghai, China). SH-SY5Y was infected with LV18-circRNA-TBC1D4, according to the manufacturer’s instructions, to overexpress circRNA-TBC1D4, and the cells were denoted SH-SY5Y-circRNA-TBC1D4. SH-SY5Y was also infected with LV18NC to serve as a control, and the cells were denoted SH-SY5Y-NC.

Cell Proliferation Assays

Cell proliferation assays were performed with Cell Counting Kit-8 (CCK-8) (Dojindo Laboratories, Kumamoto, Japan). Approximately 5×10³ cells with 100 μL DMEM/F12 were seeded in 96-well plates. At 0, 24, 48, 72, 96 and 120 h, 10 μL CCK-8 reagent was added to each well and incubated at 37 °C for 2 h. The optical density (OD) at 450 nm was read. Every 24 h, five replicates were examined for SH-SY5Y-circRNA-TBC1D4 and SH-SY5Y-NC.

Colony Formation Assays

Approximately 2×10³ SH-SY5Y-circRNA-TBC1D4 and SH-SY5Y-NC cells were seeded into six-well plates and incubated for approximately 2 weeks. Then, the colonies were fixed with 4% paraformaldehyde for 30 min and stained with 1% crystal violet for 15 min at room temperature. The cell colonies were photographed.

Cell Migration Assays

Cell migration assays were performed with Transwell permeable supports with 6.5-mm insert, 24-well plate and 8.0-μm polycarbonate membrane (Corning, NY, USA). Approximately 1×10⁴ SH-SY5Y-circRNA-TBC1D4 and SH-SY5Y-NC cells in 100 μL DMEM/F12 without FBS were placed into the upper well, and 600 μL DMEM/F12 with 10% FBS was placed into the lower well. After incubation for 36 h, the cells on the inserts were fixed with 4% paraformaldehyde for 30 min, stained with 1% crystal violet for 15 min and then photographed.

Statistical Analysis

Statistical analysis was performed by SPSS 22.0 (IBM, Armonk, NY, USA). Graphical depiction of data was generated by GraphPad Prism 8 (GraphPad Software, Inc., La Jolla, CA, USA). Comparison of circRNA expression levels between NB and adjacent normal fetal adrenal medulla was analyzed by paired Student t-test. Correlations between dysregulated FoldChange of specific circRNAs and other clinical data were performed by Spearman correlation analysis. A P value <0.05 was considered statistically significant.

Results

Expression Profiles of circRNAs in NB (Reference PRJNA721263 in SRA Database)

A total of 57,857 circRNAs were identified in five randomly selected pairs of samples by RNA-seq. The number of identified circRNAs in each sample is shown in Figure S1A. The average length was 3315.65 bp, range 38–99,918 bp (detailed length distribution is shown in Figure S1B). GC content frequency distribution and chromosome/scaffold distribution were also analyzed (Figure S1C and D). The most common circRNA form analyzed by FEELnc software¹⁶ was sense (direction), genetic (type) and exonic (location), which accounted for 90.7% of total circRNA (Figure S1E). Correlation tests, PCA, as well as hierarchical clustering analysis, all suggested that different expression patterns of circRNAs could be distinguished in the NB group and adjacent normal fetal adrenal medulla control group, but expression patterns were very similar in the same group (Figure S1F, G and H).

Bioinformatics Analysis of Differentially Expressed circRNAs

With the screening criteria for up- and downregulated circRNAs set as FoldChange >2.0 and P < 0.05, 4704 differentially expressed circRNAs were identified between NB and adjacent normal fetal adrenal medulla, including 2462 upregulated and 2242 downregulated circRNAs. Volcano plots visualized the significant differentially expressed circRNAs between the two groups (Figure 1A). The unsupervised hierarchical clustering analysis of differentially expressed circRNAs showed that the expression patterns of dysregulated circRNAs were distinguishable between the NB group and adjacent normal fetal adrenal medulla group (Figure 1B). Functional annotation of the host genes of differentially expressed circRNAs was conducted through GO analysis and KEGG pathway analysis. The most significantly enriched GO terms in biological process, cellular component and molecular function were mRNA pseudouridine synthesis, plasma membrane and acetyl-CoA carboxylase activity, respectively (the top 10 significantly enriched GO terms in the above three categories are shown in Figure 1C). KEGG pathway analysis showed that the most related pathways were fluid sheer stress and atherosclerosis (Figure 1D).

Figure 1 (A) Volcano plot showing dysregulated circRNAs between neuroblastoma specimens and adjacent normal fetal adrenal medulla. The X-axis represents log2FoldChange, and the Y-axis represents -log10P value. The horizontal gray line indicates a P value of 0.05. The two vertical gray lines indicate a fold increase or decrease of 2 (left, log2FoldChange = −1; right, log2FoldChange = 1). The red and gray points represent significantly up-regulated (n = 2462) and down-regulated (n = 2242) circRNAs in neuroblastoma, respectively. (B) Hierarchical clustering of the dysregulated circRNAs in neuroblastoma specimens and adjacent normal fetal adrenal medulla. Each red or blue bar represents higher or lower expression of circRNA in each specimen. The circRNA expression patterns were distinguishable between the neuroblastoma group and the adjacent normal fetal adrenal medulla group, but similar within the same group. (C) Functional annotation of the host genes of differentially expressed circRNAs, conducted on the basis of gene ontology (GO) analyses. The X-axis represents the GO terms, and the Y-axis represents -log10Pvalue. The green, blue and red bars represent the top ten significantly enriched GO terms in biological processes, cellular components and molecular functions. (D) Bulb map of KEGG pathway analysis for the host genes of differentially expressed circRNAs. The X-axis represents enrichment score, and the Y-axis represents the top 20 enriched pathways. The size of each bulb represents the number of enriched dysregulated circRNAs, and the P values are represented by a color scale.

Validation of Selected Differentially Expressed circRNA by Real-Time qRT-PCR

Among the 4704 differentially expressed circRNAs, 10 highly upregulated circRNAs and five highly downregulated circRNAs, which are associated with tumorigenesis and progression, were selected for expression validation by real-time qRT-PCR. Eight of 10 upregulated circRNAs and two of five downregulated circRNAs were finally verified in 20 paired samples (Figure 2A and B). Our previous study had shown that miR-21 plays a significant role in the regulation of proliferation and apoptosis in NB,^6,17 so we selected another seven circRNAs for analysis, whose predicted target miRNA was miR-21. Real-time qRT-PCR confirmed that three miR-21-related circRNAs, including circRNA-TBC1D4 (average FoldChange = 18.7; P = 0.0002), circRNA-NAALAD2 (average FoldChange = 5.6; P = 0.0032) and circRNA-TGFBR3 (average FoldChange = 23.3; P = 0.0079), were significantly downregulated in NB samples (Figure 2C). The binding sites are shown in Figure 3. Detailed information about the 22 selected circRNAs are described in Table S1.

Figure 2 Real-time qRT–PCR comparison of the relative expression (normalized to GAPDH) of 22 selected circRNAs in 20 paired neuroblastoma specimens and adjacent normal fetal adrenal medulla. Student t-test P values are shown in the figure. (A) Among ten selected highly up-regulated circRNAs from circRNA-seq, eight were confirmed to be significantly up-regulated in neuroblastoma, with P < 0.05. (B) Among five selected highly down-regulated circRNAs from circRNA-seq, two were confirmed to be significantly down-regulated in neuroblastoma, with P < 0.05. (C) Among seven selected miR-21 related circRNAs from circRNA-seq, three were confirmed to be significantly down-regulated in neuroblastoma, with P < 0.05.

Figure 3 Binding sites between miR-21 and (A) circRNA-TBC1D4, (B) circRNA-NAALAD2 and (C) circRNA-TGFBR3.

Correlation Between Significantly Downregulated miR-21-Related circRNAs and Clinical Features

To determine whether circRNA-TBC1D4, circRNA-NAALAD2 and circRNA-TGFBR3 were associated with clinical characteristics of NB, Spearman correlation analysis was performed to calculate the correlation coefficients between downregulated FoldChange of circRNA-TBC1D4, circRNA-NAALAD2, circRNA-TGFBR3 and clinical features in 20 patients (Table 2). The results showed that circRNA-TBC1D4 was associated with MYCN status (P = 0.0488; R = 0.6344) and the concentration of LDH (P = 0.0109; R = 0.7590), circRNA-NAALAD2 was associated with the concentration of LDH (P = 0.0214; R = 0.7101), and circRNA-TGFBR3 was associated with the histology category (P = 0.0001; R = 0.9692). In addition, we also found that downregulated circRNA-TBC1D4, circRNA-NAALAD2, circRNA-TGFBR3 were associated with each other, indicating that the three circRNAs may be downregulated simultaneously in NB (Table 2).

Table 2 Correlations of Fold Decreases Among circRNA-TBC1D4, circRNA-NAALAD2, circRNA-TGFBR3 and Clinical Features in 20 Patients

Overexpression of circRNA-TBC1D4 Promotes CRC Cell Migration

In order to investigate the biological functions of circRNA-TBC1D4 in neuroblastoma, we overexpressed circRNATBC1D4 in the SH-SY5Y neuroblastoma cell lines, because circRNA-TBC1D4 was downregulated in neuroblastoma samples. LV18-circRNA-TBC1D4 infection significantly increased circRNA-TBC1D4 expression in SH-SY5Y cells (Figure 4A). Cell proliferation assays using CCK8 showed that overexpression of circRNA-TBC1D4 did not decrease the proliferation of SH-SY5Y cells (Figure 4B). Colony formation assays showed that overexpression of circNRA-TBC1D4 did not decrease the colony-forming ability of SH-SY5Y cells (Figure 4C). However, cell migration assays showed that overexpression of circNRA-TBC1D4 significantly decreased the migration ability of SH-SY5Y cells (Figure 4D).

Figure 4 Overexpression of circRNA-TBC1D4 inhibited SH-SY5Y cell migration, but not proliferation and colony-formation. (A) SH-SY5Y cells were infected with lentivirus LV18-circTBC1D4 and LV18NC. qRT-PCR was conducted to detect circRNA-TBC1D4 levels. The bar graph showed that circRNA-TBC1D4 level of SH-SY5Y infected with LV18-circTBC1D4 was significantly higher than SH-SY5Y infected with LV18NC. (B) CCK-8 assay of SH-SY5Y with circRNA-TBC1D4 overexpression and NC on the indicated days. The growth curve showed that overexpression of circRNA-TBC1D4 did not inhibit the proliferation of SH-SY5Y cells. There is no significant difference between SH-SY5Y with circRNA-TBC1D4 overexpression and NC. (C) Colony formation assays of SH-SY5Y with circRNA-TBC1D4 overexpression and NC. The left panel visually showed that overexpression of circNRA-TBC1D4 did not inhibit the colony-forming ability of SH-SY5Y cells. And the bar graph showed that there was no significant difference between SH-SY5Y with circRNA-TBC1D4 overexpression and NC. (D) Transwell migration assay of SH-SY5Y with circRNA-TBC1D4 overexpression and NC. The left panel visually showed that overexpression of circNRA-TBC1D4 significantly inhibited the migration ability of SH-SY5Y cells. The bar graph showed that there was significant difference between SH-SY5Y with circRNA-TBC1D4 overexpression and NC. Data are shown as means ± s.d. of at least three independent experiments. ****P < 0.0001; ns, P > 0.05.

Discussion

The circRNAs belong to the family of ncRNAs implicated in cancer progression and tumorigenesis. Emerging studies have shown that circRNAs are ubiquitous in eukaryotes and play essential roles in modulating gene expression and signaling pathways by sponging miRNA or other molecules including proteins that participate in carcinogenesis.¹⁸ Some studies have found that circRNA can be used as a diagnostic or prognostic indicator. For instance, Peng et al reported that hsa_circ_0081001 may act as a diagnostic biomarker in osteosarcoma because it is upregulated in tissues and serums and is correlated the poor prognosis.¹⁹ Chen et al suggested that hsa_circ_0000190 may have better sensitivity and specificity in the diagnosis of gastric cancer compared with commonly used biomarkers such as CEA and CA19-9²⁰ Other studies have found that circRNA can be a potential therapeutic target which is associated with cancer pathogenesis. Liu et al found that circ0001361 plays an oncogenic role in bladder cancer invasion and metastasis by targeting the miR-491-5p/MMP9 axis, and this circRNA may be a potential novel target for bladder cancer therapy.²¹ Hsa_circ_0001944 is another tumor-promoting circRNA in bladder cancer that functions as an endogenous miRNA competitor to regulate PROK2 expression by binding to miR-548.⁸ Chen et al suggested that oncogenic circAGO2 drives cancer progression by facilitating HuR-repressed functions of AGO2-miRNA complexes. Hsa_circ_000984 affects colorectal cancer cell invasion, migration and growth by competing with cell cycle-associated proteins for sponging of miR-106b.²²

Despite the above studies, the expression and function of circRNAs in NB remain unclear. In the present study, we screened circRNA expression profiles in five paired NB samples by RNA sequencing. The results indicated there were 2462 upregulated and 2242 downregulated circRNAs. Our previous study focused on the oncogenic gene miR-21 in NB and found that miR-21 may regulate the potential targets PTEN/PDCD4 to promote apoptosis.^6,17 In addition, miR-21 functions as a significant common carcinogen in a variety of cancers, such as gliomas,²³ colon cancer,²⁴ cutaneous malignant melanoma²⁵ and retinoblastoma.²⁶ Three miR-21-related circRNAs, circRNA-TBC1D4, circRNA-NAALAD2 and circRNA-TGFBR3, were confirmed to be overexpressed in adjacent adrenal tissue compared to tumor tissue and this was verified by qRT-PCR. Furthermore, we correlated these circRNAs with clinical features, which indicated that circRNA-TBC1D4 was associated with MYCN status and value of LDH, while circRNA-NAALAD2 was associated with value of LDH and circRNA-TGFBR3 was associated with the Favorable Histology category. These correlations showed that patients with NB with downregulation of circRNA-TBC1D4, circRNA-NAALAD2 and circRNA-TGFBR3 expression might have poorer clinical outcomes. Meanwhile, circRNA-TBC1D4, circRNA-NAALAD2 and circRNA-TGFBR3 were correlated with each other, indicating that the three circRNAs may synergistically downregulate the expression of miR-21 in NB. We additionally overexpressed circRNA-TBC1D4 in the NB cell line SH-SY5Y. Cell proliferation assays, colony formation assays and cell migration assays showed that overexpression of circRNA-TBC1D4 promotes NB cell migration, but not proliferation and colony-formation in vitro. These results suggested that circRNA-TBC1D4, circRNA-NAALAD2 and circRNA-TGFBR3 are often downregulated simultaneously in NB, and may serve as biomarkers indicating an unfavorable NB type.

Conclusion

In conclusion, we performed high-throughput sequencing of NB. Through clinical correlation analysis and bioinformatics analysis, we determined that circRNA-TBC1D4, circRNA-NAALAD2 and circRNA-TGFBR3 were cancer suppressor genes, which act by sponging miR-21. We intend in future studies to further investigate the effects of these genes on the biological behavior of NB in vivo and in vitro.