This article is mentioned in:

Introduction

Osteosarcoma (OS) is one of the most common primary bone malignancies; in 2003–2007, for most countries, OS represented 20–40% of all bone cancers) (1). OS is characterized by its aggressiveness and high metastatic potential, coupled with rapid rates of progression and chemoresistance (2). Due to developments in novel resection surgical techniques and multidrug adjuvant chemotherapeutic methods, the 5-year survival rates have improved to 70–80% by the mid-1980s (3). However, the 5-year survival rate of patients showing chemoresistance is markedly lower at <20%; patients with OS may demonstrate superior outcomes with additional therapies, including small molecule targeted agents (such as endothelin-1) (4), but these strategies are not making breakthroughs in clinical trials because of the complex biology of OS (5,6). Therefore, novel therapeutic approaches for OS, especially those with complex gene regulatory networks, are needed.

Over the past decades, new classes of non-coding RNAs (ncRNAs) have been discovered, including circular RNAs (circRNAs), microRNAs (miRNAs) and long non-coding RNAs (lncRNAs) (3–5). During this time period, different regulatory functions associated with ncRNAs have been revealed. In particular, ncRNAs have been reported to effectively regulate essential protein effectors of cellular function that are important for the development of OS (6–9). Knockout of MALAT1 can reduce the expression of RhoA, cell nuclear antigen and vascular growth factor (angiomotin), thereby inhibiting the proliferation and invasion of human OS cells and inhibiting their metastasis (10). Amongst the various known ncRNA families, circRNAs are conserved types of special RNAs with covalent closed-loop properties, rendering them highly stable and more resistant to degradation by endonucleases (11). CircRNAs are widely found in various types of tissue and organs such as circRNA ZKSCAN1 in non-small cell lung cancer, circRNA CAMSAP1 in colorectal cancer, circRNA FBXW4 in trophoblast cell and circRNA MAN1A2 in esophageal squamous cell carcinoma (12). Circular RNAs may function similarly to regulate the activity of other miRNAs. circular RNAs may bind and sequester RNA-binding proteins or even base pair with RNAs besides microRNAs, resulting in the formation of large RNA-protein complexes (13). Previous studies have reported that circRNAs participate widely in disease development such as type2 diabetes, hypertension, spondyloarthropathies, osteoarthritis and serve important roles in cancer, such as thyroid and cervical cancer (14,15). In addition, circRNAs have been reported contribute to the onset and development of OS by controlling proliferation, invasive and metastatic ability, in addition to apoptosis, tumor cell metabolism and drug resistance (16–19). Therefore, understanding the function of circRNAs is of importance for elucidating the molecular mechanism underlying the development of OS whilst overcoming its chemoresistance potential.

In the present study, the differential expression of hsa_circ_0064636 in OS was assessed and its regulatory mechanism was evaluated using bioinformatic methods to determine its expression profile and prognostic significance. The aim was to provide a novel experimental basis and direction for future studies on molecular markers of OS.

Materials and methods

Materials and equipment

Human normal osteoblast (hFOB1.19) and OS (HOS, SJSA-1 and MG63) cell lines were purchased from Shenzhen Haodi Huatuo Biotechnology Co., Ltd. TRIzol reagent was purchased from Invitrogen; Thermo Fisher Scientific, Inc. PrimeScript™ RT reagent Kit with gDNA Eraser and TB Green® Fast qPCR Mix were purchased from Takara Bio, Inc. A LightCycler® 96 real-time quantitative PCR instrument was purchased from Roche Diagnostics.

Cell culture

hFOB1.19 cells were grown in DMEM/F12 supplemented with 10% FBS with the addition of geneticin (0.3 mg/ml; all Thermo Fisher Scientific, Inc.), penicillin (100 U/ml) and streptomycin (100 U/ml). By contrast, the OS cell lines were grown in DMEM/F12 medium supplemented with 10% FBS, streptomycin (100 U/ml) and penicillin (100 U/ml). All cells were incubated at 37°C in a humidified atmosphere with 5% CO2. Medium was replaced with fresh complete medium every 2–3 days for further incubation and when the cell density reached >90%, before the cells were passaged at a ratio of 1:3.

RNA extraction and reverse transcription-quantitative (RT-q)PCR

Total cellular RNA was extracted using TRIzol and cDNA was synthesized using the miRNA First Strand cDNA Synthesis Kit (Sangon Biotech Co., Ltd.) at 25°C for 5 min, 42°C for 30 min, 85°C for 5 min. qPCR was performed using the SYBR green PCR mix (TB Green® Fast qPCR Mix) and divergent primers from hsa_circ_0064636 were used with reaction conditions as follows: Pre-denaturation at 95°C for 10 min, followed by denaturation at 95°C for 15 sec and annealing at 60°C for 30 sec for 45 cycles. The expression levels of GAPDH (circRNA) and U6 (miRNA) were used as internal controls. The relative expression of the screened genes was calculated using the 2−ΔΔCq method (20). The primers used in the present study are presented in Table I; primer sequence for U6 is not available.

Table I. Primer sequences used for reverse transcription-quantitative PCR.

Table I.

Primer sequences used for reverse transcription-quantitative PCR.

Gene	Sequence (5′-3′)
Hsa_circ_0064636	F: GCTTCCCCTGTCTCCACATA
	R: ATGTCCAAAGGGTTTCAGCA
GAPDH	F: CCACTCCTCCACCTTTGAC
	R: ACCCTGTTGCTGTAGCCA
Ubiquitination	F: TCCAGAGAACCTGCTACCCT
factor E4A	R: AGTTACATCTTCAAAATGGG
	CTCC
Voltage dependent	F: GGAAGGCAGAAGATGGCTGT
anion channel 1	R: GTCCACGTGCAAGCTGATCT

[i] F, forward; R, reverse.

Construction of the circRNA/miRNA/mRNA network

miRNAs interacting with hsa_circ_0064636 were predicted using the circRNA target miRNA prediction tools RNAhybrid (21), TargetScan (22) and miRanda (23) before being subjected to intersection analysis to screen out miRNAs that were commonly identified by the three prediction tools. For miRNA targeting, gene prediction tools, including RNAhybrid, TargetScan, miRanda, miRWalk (24), miRMap (25) and miRNAMap (26) were used to predict target genes for the intersecting miRNAs. These were compared with the screened miRNAs, before seven (six database prediction ensembles and a variance analysis) intersects from the differentially expressed miRNAs were finally identified. The selected circRNA-miRNAs and miRNA-mRNAs were then used to construct and visualize the regulatory network of circRNA/miRNA/mRNA using Cytoscape 3.8.0 (cytoscape.org/). Venn diagrams of predicted miRNAs or target genes were plotted using the R (version 4.2.3) package ‘Venn’ (version 1.11).

Analysis of the differential expression of miRNAs

The GSE65071 (9 normal samples from healthy individuals, 14 OS samples) and GSE16088 (9 normal samples, 14 OS samples) expression datasets and RNA-seq ere downloaded from the Gene Expression Omnibus (GEO; http://www.ncbi.nlm.nih.gov/geo/) database and the ‘affy’ (to read gene expression microarray data from affymetix) package (version 3.18) was used. The ‘Limma’ (version 3.54.2).package was used to identify differentially expressed genes between normal and tumor samples with log2[fold change (FC)]>1 and adjusted P<0.05 set as the significance criteria. The differentially expressed genes were plotted using the ‘ggplot2’ (version 3.4.3) package and a volcano plot was produced.

Survival analysis

PROGgeneV2 (progtools.net/gene/filter.php) (27) was used along with aforementioned GEO data to investigate the prognostic significance of genes. The target genes obtained from the aforementioned intersection) were input into a PROGgeneV2 before the dataset for OS was selected and the overall survival map was constructed based on the median expression of a given gene for classification into high and low expression groups for plotting Kaplan-Meier curves with PROGgeneV2. ‘PROGeneV2’ was used for hypothesis testing using the ‘Suvival’ (version 3.5.3) package in R (version 4.2.3). Statistical analysis was performed using the log-rank test and the threshold for a meaningful survival prognosis was set as P<0.05. Based on the survival analysis, a new circRNA/miRNA/mRNA visual regulatory sub-network was obtained using Cytoscape (cytoscape.org/) 3.8.0.

Statistical analysis

GraphPad Prism 8 (Dotmatics) software was used to statistically analyze the experimental data. The test data were normally distributed. Data are presented as the means ± SD of three independent experiments. P<0.05 was considered to indicate a statistically significant difference. Statistical comparison between groups was evaluated using one-way ANOVAwith the Dunnett's post hoc test.

Results

Analysis of differentially expressed miRNAs and mRNAs in OS

The expression of hsa_circ_0064636 was found to be significantly higher in the three human OS cell lines compared with that in the hFOB1.19 osteoblasts (Fig. 1A). A total of 114 miRNAs were shown to be upregulated, whilst 117 were downregulated in OS (compared with normal group) based on the miRNA (GSE65071) expression profile data. Differentially expressed miRNAs were visualized using a volcano plot (Fig. 1B). Differential analysis of the mRNA expression (GSE16088) spectrum data identified 716 mRNAs that were downregulated and 1,924 mRNAs that were upregulated(OS versus norma). Differentially expressed mRNAs were visualized using a volcano plot (Fig. 1C).

Figure 1. Hsa_circ_006463 expression validation and prediction of differentially expressed miRNAs. (A) Expression of hsa_circ_006463 in osteosarcoma cell lines. **P<0.01, ***P<0.001. Volcano plots of (B) Differential miRNA expression, with the positions of miR-326 and miR-503-5p labeled and (C) Differential mRNA expression with the positions of VDAC1 and U8E4A labeled. (D) Venn diagram of hsa_circ_0064636 target miRNAs predicted using the TargetScan, miRanda and RNAhybrid databases, with the two most commonly miRNAs labeled on the right side of the figure. DEG, differentially expressed gene; miR, microRNA.

Hsa_circ_006463 targets the predicted miRNA

To assess the regulatory mechanism of hsa_circ_006463 in OS, miRNAs targeted by hsa_circ_006463 were predicted using three databases. TargetScan, miRanda and RNAhybrid predicted 498, 965 and 226 target miRNAs, respectively. Differential expression analysis yielded 88 target (downregulated) miRNAs. Intersecting miRNAs targeted by hsa_circ_006463 according to all four different databases used were found to be miR-326 and miR-503-5p (Fig. 1D).

Prediction and analysis of miR-326 and miR-503-5p targets

The regulatory action of hsa-miR-326 and miR-503-5p in OS was then assessed based on the possible binding of target mRNAs. miR-326 target genes were predicted using the miRWalk, miRanda, miRMap, miRNAMap, RNAhybrid and TargetScan databases, with 5,071, 3,364, 6,616, 6,716, 16,373 and 4,830 target genes identified, respectively. Intersection with 1,924 upregulated differential genes identified 31 target genes using Venn diagram analysis (Fig. 2A). Target genes of miR-503-5p were also analyzed using the miRWalk, miRanda, miRMap, miRNAMap, RNAhybrid and TargetScan, with 2,034, 844, 6,692, 1,496, 12,416 and 2,196 target genes identified, respectively. These were subjected to intersection analysis with 1,924 differentially upregulated genes and 31 target genes were identified (Fig. 2B).

Figure 2. Target gene prediction and differential gene expression of miRNAs. Venn diagram showing the intersection of data from six databases and differentially expressed genes(OS compared to normal group) to predict binding target genes for (A) miR-326 and (B) miR-503-5p. The names of the target mRNAs are shown on the right. DEG, differentially expressed gene; UBE4A, ubiquitination factor E4A; VDAC1, voltage dependent anion channel 1; miR, microRNA.

Construction of the circRNA/miRNA/mRNA network

In total, six databases were used to predict binding targets of miR-326 and miR-503-5p, yielded 31 and 6 target gene relationships, respectively. These were used to construct relationship networks (Fig. 3). Hsa_circ_0064636 was therefore predicted to target miR-326 and miR-503-5p. The potential target genes of miR-326 included TNC, HNRNPA2B1, ATXN1, DPYSL2, RGS3, ubiquitination factor E4A (UBE4A), MRC2, PSMC1, FOXO3, SRRM1, SAMD4A, MYO1D, ADAM17, PDIA4, VGLL4, FYN, NFATC3, VCL, SPAG9, MED13L, RAI14, MYOF, NDEL1, TUBB, ALPL, USP11, CEP350, LSM1, GOLPH3, C9orf3 and CUX1 (31 target genes). By contrast, miR-503-5p target genes were PSMD7, INSR, PTPN12, TCF7L2, voltage dependent anion channel 1 (VDAC1) and TPM1 (6 target genes).

Figure 3. CircRNA/miRNA/mRNA competing endogenous RNA regulatory network starting with hsa_circ_0064636. The red circle represents circRNA, blue triangles represent miRNAs and yellow squares represent mRNAs. CircRNA, circular RNA; miR, microRNA.

Survival analysis

Survival curves for overall survival were plotted for 31 of the 37 target genes, as shown in Table II. For survival analysis of the target genes, only two differences were found significant, namely UBE4A and VDAC1 (Fig. 4). This implies that patients with osteosarcoma have poorer prognostic survival when UBE4A and VDAC1 are highly expressed. According to Fig. 2, UBE4A was a target gene of miR-326 whereas VDAC1 was a target gene of miR-503-5p.

Figure 4. Survival curves of target genes. Survival curves were plotted for (A) VDAC1 and (B) UBE4A based on Kaplan-Meier analysis using the ‘PROGeneV2’ online platform. HR, hazard ratio; UBE4A, ubiquitination factor E4A; VDAC1, voltage dependent anion channel 1.

Table II. Target gene survival analysis.

Table II.

Target gene survival analysis.

Gene	miRNA	Hazard ratio	LCI (95%)	UCI (95%)	P-value
TNC	hsa-miR-326	0.80	0.59	1.09	0.16087
ATXN1	hsa-miR-326	1.06	0.45	2.51	0.89614
DPYSL2	hsa-miR-326	0.72	0.31	1.66	0.44617
RGS3	hsa-miR-326	0.59	0.21	1.61	0.30120
UBE4A	hsa-miR-326	1.60	1.00	2.57	0.04894
MRC2	hsa-miR-326	0.81	0.38	1.73	0.59453
PSMC1	hsa-miR-326	0.38	0.10	1.40	0.14712
SRRM1	hsa-miR-326	1.11	0.35	3.51	0.85323
SAMD4A	hsa-miR-326	0.75	0.25	2.28	0.61097
MYO1D	hsa-miR-326	0.76	0.4	1.45	0.40852
ADAM17	hsa-miR-326	0.4	0.06	2.86	0.36201
PDIA4	hsa-miR-326	0.81	0.38	1.74	0.58828
VGLL4	hsa-miR-326	1.03	0.46	2.27	0.95088
FYN	hsa-miR-326	1.66	0.7	3.95	0.25158
NFATC3	hsa-miR-326	4.94	0.76	32.03	0.09380
VCL	hsa-miR-326	0.60	0.29	1.24	0.16888
SPAG9	hsa-miR-326	0.45	0.18	1.17	0.10087
RAI14	hsa-miR-326	0.33	0.16	0.68	0.00261
NDEL1	hsa-miR-326	0.81	0.35	1.9	0.62756
TUBB	hsa-miR-326	1.59	0.92	2.78	0.09927
ALPL	hsa-miR-326	1.07	0.83	1.39	0.59196
USP11	hsa-miR-326	0.5	0.18	1.38	0.18324
CEP350	hsa-miR-326	0.88	0.35	2.21	0.78379
LSM1	hsa-miR-326	1.58	0.78	3.22	0.20742
GOLPH3	hsa-miR-326	1.54	0.94	2.5	0.08360
PSMD7	hsa-miR-503-5p	0.70	0.32	1.54	0.37880
INSR	hsa-miR-503-5p	0.13	0.01	1.79	0.12555
PTPN12	hsa-miR-503-5p	1.62	0.89	2.94	0.11362
TCF7L2	hsa-miR-503-5p	3.99	0.79	20.11	0.09396
VDAC1	hsa-miR-503-5p	4.09	1.26	13.29	0.01919
TPM1	hsa-miR-503-5p	0.73	0.47	1.16	0.185767
HNRNPA2B1	hsa-miR-326	-	-	-	-
FOXO3	hsa-miR-326	-	-	-	-
MYOF	hsa-miR-326	-	-	-	-
C9orf3	hsa-miR-326	-	-	-	-
CUX1	hsa-miR-326	-	-	-	-

[i] UBE4A, ubiquitination factor E4A; VDAC1, voltage dependent anion channel 1; miR or miRNA, microRNA; LCI, lower confidence interval; UCI, upper confidence interval; -, data not available.

CircRNA/miRNA/mRNA network has prognostic significance in OS

Based on the original circRNA/miRNA/mRNA competing endogenous RNA interaction network, though RAI14 was significant in the survival analysis, RAI14 is highly expressed in osteosarcoma compared with normal group, but survival analysis of RAI14 low expression group has better survival prognosis). Thus only two target genes with prognostic significance for survival (UBE4A and VDAC1) were retained before a new sub-network with prognostic implications was created (Fig. 5A).

Figure 5. CircRNA/miRNA/mRNA competing endogenous RNA regulatory network based on the target genes with prognostic significance. (A) Red circles represent circRNAs, Blue triangles represent miRNAs and yellow squares represent mRNAs. Expression of (B) hsa-miR-326, (C) hsa-miR-503-5p, (D) UBE4A and (E) VDAC1 in osteosarcoma and normal osteoblast cell lines. *P<0.05, **P<0.01, ***P<0.001 vs. hFOB1.19. (F) Schematic representation of target genes and miRNA binding sites. CircRNA, circular RNA; UBE4A, ubiquitination factor E4A; VDAC1, voltage dependent anion channel 1; miRNA, microRNA.

Validation of differential miRNA and mRNA expression in OS cell lines

Differential expression of hsa_circ_0064636 was next assessed using RT-qPCR in human (HOS, SJSA-1, MG63) OS cell lines and the human osteoblast (hFOB1.19) cell line. The expression of hsa_circ_0064636 was found to be significantly higher in the human OS cell line compared with that in the osteoblast cell line (Fig. 1A). The expression levels of hsa-miR-326 (Fig. 5B) and hsa-miR-503 (Fig. 5C) were significantly lower in OS cell lines compared with those in the osteoblast cell line. UBE4A (Fig. 5D) and VDAC1 (Fig. 5E) expression were found to be significantly increased in the OS cell lines compared with that in the human osteoblast cell line. A schematic representation of the target gene and corresponding miRNA binding site are presented in Fig. 5F.

Discussion

OS is one of the most common malignancies in the bone, Treatment of OS remains a major challenge. Existing treatment options for OS include as surgery and chemotherapy, but the prognosis of patients remains unsatisfactory (28). Although four genes [BRCA1, mutS homolog 2, CCND1 (cyclin D1) and ITGA5 (integrin subunit alpha 5)] have been documented to associate with the development and progression of OS, the focus has been only on protein-coding genes or miRNAs (29). In addition, molecular mechanisms underlying the development and progression of OS remain poorly understood. A number of ncRNAs have been reported to be important in the development of OS (30–33). The lncRNA TMPO-AS1 was previously found to decreases miR-199a-5p/WNT7B(elevated) axis, thereby functioning as a ceRNA that promotes tumorigenesis in OS (34). In addition, miR-1236-3p overexpression inhibits cell proliferation and induces apoptosis by targeting Krueppel-like factor 8 (35).

Unlike other ncRNAs, including lncRNAs and miRNAs, circRNAs are highly conserved and stable in mammalian cells. These properties render them potentially ideal biomarkers and therapeutic targets (36). Previous studies have reported the important role of circRNAs in different processes, including tumorigenesis, development and metastasis, of tumors such as gastric (37), colorectal (38), lung (39) and cervical (40) cancer. Circ-transcriptional Adaptor 2A was previously found to be differentially expressed between OS cell lines (HOS, 143B,U2OS,SJSA-1,MG63) and corresponding non-cancer cell lines (HEK-293 and hFOB1.19), with higher expression in OS (41). In the present study, the significantly upregulated expression of hsa_circ_0064636 in OS cell lines compared with normal tissue cell lines was found before a regulatory network of its miRNA targets and downstream mRNA targets was constructed.

In the present study, miR-326 and miR-503-5p were predicted to be the target miRNAs of hsa_circ_0064636 using multiple databases, yielding circRNA/miRNA interactions. miR-326 and miR-503-5 were identified from the GSE65071 dataset and were downregulated in OS samples compared with normal samples, before their expression in OS samples was verified. Here, hsa_circ_0064636 was found to be significantly upregulated in OS, miR-326 to promote cervical cancer progression through up-regulation of ELK1 (42). The overexpression of miR-326 can lead to the inhibition of proliferation and invasion, whilst inducing apoptosis and autophagy in cervical cancer cells (42). In addition, miR-326 expression was previously found to be significantly downregulated in prostate cancer (43), which associated with the prognoses of these patients. By contrast, miR-503-5p has been reported to inhibit tumorigenesis, angiogenesis and lymphangiogenesis in colon cancer through the direct inhbition of VEGF-A expression (44). In hepatocellular carcinoma, miR-503-5p was documented to regulate epithelial-to-mesenchymal transition, metastasis and prognosis through the inhibition of WEE1 (45). miR-326 and miR-503-5p are key for the suppression of numerous tumors such as prostatic carcinoma and colon cancer (44,46), where their expression is low in certain cancers such as miR-326 is low expressed in lung adenocarcinoma and miR-503-5p is low expressed in colon cancer. Therefore, hsa_circ_0064636 may promote OS development by suppressing the expression of miR-326 and miR-503-5p. To the best of our knowledge, studies on miR-326 and miR-503-5p in OS are lacking.

The present study performed miRNA prediction using miRWalk, miRanda, miRMap, miRNAMap, RNAhybrid and TargetScan. Significantly differentially expressed genes in OS were used to screen for potential target genes of miR-326 and miR-503-5p to obtain a circRNA/miRNA/mRNA ceRNA interaction network. Subsequently, mRNAs with prognostic significance were also screened to obtain a circRNA/miRNA/mRNA regulatory sub-network. UBE4A was identified to be a potential direct target of miR-326 whereas VDAC1 was found to be a potential direct target of miR-503-5p. A significant difference in prognosis between UBE4A and VDAC1 was confirmed in a survival analysis of OS patients.. The results of survival analysis showed that the prognosis of patients was significantly worse in the high UBE4A and VDAC1 expression group compared with that in the low expression group.

UBE4A belongs to the U-box ubiquitin ligase class of enzymes. The encoded protein is involved in polyubiquitin chain assembly and regulates chromosome condensation and polyubiquitination segregation through securing (47). Auto-antibodies against the encoded protein (recombinant C-terminal UBE4A Protein) are potential markers for scleroderma and Crohn's disease (47,48). UBE4A was found to be significantly upregulated in ovarian plasmatic cystic carcinoma compared with the adjacent normal tissues. The VDAC1 channel forms a major component of the outer mitochondrial membrane (49). VDAC1 is expressed in all compartments, including mitochondria, plasma membrane. This protein regulates major metabolic and energetic functions in the cell, including Ca2+ homeostasis, oxidative stress and mitochondria-mediated apoptosis (50,51). VDAC1 may be associated with the destruction of nerve cells, where its overexpression triggers cell death (52). Furthermore, VDAC1 was found to be a tumor promoter in cervical cancer, where its knockdown in cervical cancer cell line S12 increased the rate of apoptosis in a manner that was partially reversed by the overexpression of VDAC1 (53).

In the present study, RT-qPCR analysis demonstrated that hsa_circ_00063636, VDAC1 and UBE4A were highly expressed in OS cell lines, whereas miR-326 and miR-503-5p expression was significantly lower in OS cell lines. These results were consistent with those of the bioinformatics analysis. A limitation of the present study is that the regulatory network was not validated experimentally and further investigations are required such as experimentally verifying miRNA binding to mRNAs. In the present study, bioinformatic methods are used to predict potential regulatory networks, with emphasis on discovering new targets and potential networks. However, further verification of the regulation between miRNA and mRNA is required.

In conclusion, a sub-regulatory network of hsa_circ_0064636-miR-326/miR-503-5p-UBE4A/VDAC1 was identified with significance for survival prognosis. RT-qPCR experiments demonstrated that hsa_circ_0064636, UBE4A and VDAC1 were significantly and differentially overexpressed in OS cell lines, whilst miR-326 and miR-503-5p were downregulated. It was hypothesized that hsa_circ_0064636 may be involved in the development of OS by acting as a sponge, thereby suppressing miR-326 and miR-503-5p to facilitate the upregulation of VDAC1 and UBE4A.

Acknowledgements

Not applicable.

Funding

The present study received financial support from the National Natural Science Foundation of China (grant no. 81960400).

Availability of data and materials

The data generated in the present study may be requested from the corresponding author.

Authors' contributions

GHY and NCH confirm the authenticity of all the raw data. GHY and NCH designed the study and drafted the manuscript. CTC, JWC and HJH interpreted data. JWC and HJH revised the manuscript. All authors have read and approved the final version of the manuscript.

Ethics approval and consent to participate

Not applicable.

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Glossary

Abbreviations

Abbreviations:

References