Genome-Wide Analysis of circular RNAs and validation of hsa_circ_0006719 as a potential novel diagnostic biomarker in congenital scoliosis patients.

Congenital scoliosis (CS) is a form of spinal curvature resulting from anomalous development of vertebrae. Recent studies demonstrated that circRNAs could serve as potential biomarkers of disease diagnosis. Genome-wide circRNAs expression in seven CS patients and three healthy controls was initially detected. Bioinformatics analysis was conducted to explore the potential pathological pathway of CS. Quantitative PCR (qPCR) was performed to validate the selected circRNAs in the replication cohort with 32 CS patients and 30 healthy controls. Logistic regression controlling for gender was conducted to compare the expression difference. Receiver operating characteristic (ROC) curve analysis was performed to evaluate the diagnostic value. Twenty-two differentially expressed circRNAs were filtered from genome-wide circRNA sequencing. Seven circRNAs were validated by qPCR. Only hsa_circ_0006719 was confirmed to have a higher expression level in the CS group than the healthy control group (P = 0.036). Receiver operating characteristic curve also suggested that hsa_circ_0006719 had significant diagnostic value for CS (AUC = 0.739, P = 0.001). We described the first study of circRNAs in CS and validated hsa_circ_0006719 as a potential novel diagnostic biomarker of CS.

BACKGROUND

Congenital scoliosis (CS) is a form of spinal malformation resulting from abnormal axial bone development in embryo. The incidence of CS is approximately 0.5‐1 per 1000 newborn infants. , , CS patients always have multiple system deformity, such as missing ribs, split cord malformation and pulmonary dysfunction, which physically and psychologically affect their lives and daily activities. Although the studies to explore the pathogenesis of CS have been conducted for many years, the aetiology is still elusive. Previous studies have proved that genetic factors could lead to the development of CS, such as rearrangement of chromosome 16p11.2, mutation of TBX6 or FBN1. There is also evidence that environmental factors could contribute to the development of CS. Li et al found that vitamin A deficiency in pregnancy may induce CS in rats. Some researchers illustrated that hypoxia and high altitude are associated with a higher risk of CS. , However, there were few types of research focusing on the linkage of these aetiological factors. Thus, it is urgent and necessary to identify novel molecular markers of the mechanism, even for diagnosis and therapeutics.

Circular RNAs (circRNAs) are a class of non‐coding RNAs that manifested with stable circular RNA structure. Recent studies have reported that circRNAs play a pivotal role in different kinds of diseases. CircRNAs could act as efficient microRNA sponges, which function as post‐transcriptional regulators. Some researchers also found that some circRNAs could work as protein sponges, which can influence gene regulation. Other studies also found in cellular responses to environmental stress, some circRNAs can be translated. , Because of the unusually stable circular structure, circRNAs were also reported to be promising diagnostic biomarkers in different diseases. , There are very few studies focusing on circRNAs and CS. Our previous study has identified significantly different expressed circRNAs in vitamin A deficiency‐induced CS rats. However, the function and characteristics of circRNAs in CS patients are still unclear.

In this study, we enrolled seven patients and three healthy controls to identify and compare the expression of genome‐wide circRNAs. After filtering, annotation and validation, we compared the differentially expressed circRNAs aiming to find the potential diagnostic and therapeutic biomarkers of CS.

MATERIALS AND METHODS

Patients and materials

Seven patients diagnosed as CS and three healthy controls were recruited from Peking Union Medical College Hospital (PUMCH). Thirty‐two CS patients and 30 healthy controls were enrolled as a replication cohort. The inclusion and exclusion criteria of CS were as follows.

Inclusion criteria

Patients diagnosed as CS.

Age of onset under 18 years old.

Having complete imaging data, including X‐ray, three‐dimensional imaging of the spine CT or spinal MRI.

Exclusion criteria

Other types of scoliosis including adolescent idiopathic scoliosis, neuromuscular scoliosis, scoliosis secondary to skeletal dysplasia or connective tissue abnormalities.

Incomplete imaging data.

Having a chronic disease that influenced skeletal development.

All patients were diagnosed as CS by at least two experienced orthopaedic surgeons. The clinical information of the study subjects was summarized in Table 1. Written informed consent was obtained from all the participants or their guardians. All controls were healthy without any spinal deformity or other diseases. The Ethical Review Board of Peking Union Medical College Hospital approved this study (Protocol Number: JS‐1901).

TABLE 1

Basic characteristics of participants in the study

Characteristics	Discovery cohort		Replication cohort
	CS group (n = 7)	Control group (n = 3)	CS group (n = 32)	Control group (n = 30)
Mean age (y)	13.57	26.67	14.49	23.67
Gender M:F	1:6	3:0	17:15	10:20
Main Cobb angle (°)	88.86	NA	56.26	NA

RNA extraction and Genome‐wide circRNAs sequencing

Total RNA was extracted from the peripheral blood using a Qiagen PAXgene Blood miRNA Kit (QIAGEN, Eastwin Scientific, Inc, Beijing, China) according to the manufacturer's instructions. The quality and integrity of the extracted RNA were evaluated using an RNA 6000 Nano Lab Chip Kit (Agilent Technologies, CA, USA) and an Agilent 2100 Bioanalyzer (Agilent Technologies, CA, USA). The ribosomal RNAs were removed using a Ribo‐Zero Gold rRNA Removal Kit (Human/Mouse/Rat) (Epicenter Company, Madison, WI, USA). The linear RNAs were digested by RNase R to separate the circRNAs. The RNAs were then reverse‐transcribed to create the final cDNA library using the mRNA‐Seq sample preparations kit (Illumina, San Diego, CA, USA). Sequencing was carried out by HiSeq X‐Ten with 150‐bp pair‐end reads mode.

Bioinformatics analysis and annotation of circRNAs

The low‐quality reads were removed, and the clean reads were aligned to the reference human genome (hg19). Based on the tool of find_circ, find_circ_enhance was used to identify circRNAs with default parameters (circRNA read count ≥2 and unique alignment reads ≥2). Filtered circRNAs were used for further annotation. After comparing to circBase, novel circRNA and known circRNA were separated. Transcripts per million clean tags (TPM) were applied to quantify the expression of circRNAs. Differentially expressed circRNAs were estimated by DESeq2. The expression level was considered significantly different when |log2 fold change| > 1 and P value < 0.05. Functional enrichment analysis of the differentially expressed circRNA‐related genes at Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) was performed using clusterProfiler (v3.6.0) Bioconductor package.

Validation of selected circRNAs

Quantitative PCR (qPCR) was performed to validate the circRNAs identified by the bioinformatics analysis. Divergent primers were designed using circPimer (v1.2). First‐strand cDNA was synthesized using SuperScript III First‐Strand (Thermo Fisher Scientific, Beijing, China) according to the manufacturer's instructions. qPCR condition was as follows: an initial 1‐min denaturation at 95°C (Ramp Rate 4.4°C/s) followed by 40 cycles of 15 seconds at 95°C, 15 seconds at the appropriate annealing temperature (depending on the divergent primer set used) and 30 seconds extension at 72°C, then 5 seconds at 95°C, 60 seconds at 60 and 95°C (Ramp rate 0.11°C/s) for melt curve analysis, with a final cooling at 50°C for 30 seconds. For candidate circRNAs, the PCR products were separated from the electrophoresis agarose gel. Sanger sequencing was conducted to validate the back‐spliced site of each circRNAs. Each qPCR of each circRNA was replicated three times using β‐actin as the internal control for normalization of expression. The ΔΔCt method was performed to calculate the relative expression of circRNAs.

Statistical analysis

Data of qPCR were analysed using SPSS software (16.0 version, SPSS Inc, Chicago, IL, USA). The expression differences of qPCR were evaluated using logistic regression controlling for gender. Receiver operating characteristic (ROC) curve was performed to estimate the diagnostic application of circRNA.

RESULTS

Genome‐wide circRNA expression profiles

There were 551,278 circRNAs identified by circRNA sequencing in seven CS patients and three healthy controls. After filtering, there were 126 907 candidate circRNAs left for further analysis (Figure 1A). In order to identify the source of those circRNAs, breakpoint type was also annotated (Figure 1B). Most of them were derived from exonic or intronic circRNAs. Differentially expressed circRNAs were estimated by the fold changes (FC) of circRNA expression (Figure 1C). In total, 394 circRNAs with significantly different expression level were analysed by the volcano plot (Figure 1D, Table S1).

FIGURE 1

Characteristics of genome‐wide circRNA expression profiles of the CS group and healthy control group. A, The total number of filtered and candidate circRNAs screened. B, The source of circRNAs based on the breakpoint. C, The scatterplot of the differentially expressed circRNAs. The circRNAs above the top red line and below the bottom red line indicate more than a 1.0‐fold change between the two groups. D, The volcano plot of circRNA expression profile. The vertical lines correspond to 1.0‐fold change up‐regulation and down‐regulation and the horizontal line represents P = 0.05

We enriched the differentially expressed circRNA‐related genes and then subjected these genes to GO annotation and KEGG pathway analysis. It indicated that the circRNA‐related genes were enriched in several pathways. The top outcomes of GO enrichment were organelle organization in biological process (BP), intracellular organelle part in cellular component (CC) and transcription coactivator activity in molecular function (MF) (Figure 2A). In the KEGG pathway analysis, three pathways with the most significant association were the ubiquitin‐mediated proteolysis signalling pathway, endocytosis signalling pathway and oocyte meiosis signalling pathway (Figure 2B), indicating that these pathways could be associated with CS.

FIGURE 2

GO and KEGG analysis of the differentially expressed circRNA‐related genes in the CS group and healthy control group. A, The enriched results of biological process (BP), cellular component (CC) and molecular function (MF) by GO analysis. B, Pathways enriched in the KEGG analysis

Combing with the P value and the read count of each circRNAs, we selected 22 circRNAs with the most stable and significant different expression for validation (Table 2). Divergent primers were designed (Table S2), and qPCR was performed, and the back‐spliced sites of seven circRNAs were confirmed by Sanger sequencing including hsa_circ_0006856, hsa_circ_0006719, hsa_circ_0006208, hsa_circ_0002785, hsa_circ_0002692, hsa_circ_0002372 and hsa_circ_0000225 (Figure 3). Then, we compared the expression level with ΔΔCt method (Table S3). All of those seven circRNAs had higher expression levels in the CS group than the healthy control (HC) group (Table S4). After compared with the expression level of the discovery cohort, only hsa_circ_0006719 had the same expression tendency in the replication cohort. After logistic regression controlling for gender, we found that the expression level was statistically significant different between the CS group and the HC group (95% CI: 0.59‐0.98, P = 0.036; Table 3).

TABLE 2

The 22 selected circRNAs with the most stable and significant different expression

CircRNA	Control	Case	log2(fold change)	P
hsa_circ_0006719	77.29	276.76	1.84	4.25E−05
hsa_circ_0037173	2.17	11.78	2.44	0.01
chr3_105389076_105404310_‐	1.01	10.41	3.37	0.01
hsa_circ_0001243	1.24	9.48	2.94	0.02
hsa_circ_0041267	1.22	8.58	2.82	0.02
chr19_50840790_50865349_‐	10.55	30.27	1.52	0.02
hsa_circ_0006208	55.92	5.43	−3.36	5.41E−06
hsa_circ_0002692	120.66	22.85	−2.4	3.54E−05
hsa_circ_0000225	98.56	35.15	−1.49	0.01
chr5_72311452_72333042_+	15.36	2.69	−2.51	0.01
hsa_circ_0023233	12.17	1.34	−3.18	0.01
chr14_59730158_59758024_+	18.61	3.58	−2.38	0.01
hsa_circ_0002785	15.16	2.24	−2.76	0.01
hsa_circ_0007128	14.49	1.83	−2.99	0.01
chr9_3630937_3651867_+	27.01	5.89	−2.2	0.01
hsa_circ_0002372	13.75	2.25	−2.61	0.02
hsa_circ_0006856	69.08	27.04	−1.35	0.02
chr11_65202523_65211534_+	32.54	9.9	−1.72	0.02
chr9_3647337_3651867_+	87.83	35.67	−1.3	0.02
hsa_circ_0002280	12.89	2.25	−2.52	0.03
chr21_40578033_40584633_‐	44.01	17.36	−1.34	0.03
chr6_37250657_37284982_+	12.18	2.27	−2.43	0.04

FIGURE 3

Sanger sequencing of hsa_circ_0006856 (A), hsa_circ_0006719 (B), hsa_circ_0006208 (C), hsa_circ_0002785 (D), hsa_circ_0002692 (E), hsa_circ_0002372 (F) and hsa_circ_0000225 (G). The arrows represented the back‐spliced sites

TABLE 3

Logistic regression predicting likelihood of CS based on gender and hsa_circ_0006719

	B	SE	Wald	df	P	Odds Ratio	95% CI for Odds Ratio
							Lower	Upper
Gender	−0.95	0.59	2.603	1	0.107	0.39	0.12	1.23
has_circ_0006719	−0.28	0.12	4.40	1	0.036	0.76	0.59	0.98
Constant	1.04	0.45	5.37	1	0.020	2.83

Gender is for males compared to females.

ROC curve analysis of hsa_circ_0006719 in CS patients

To estimate the diagnostic value of hsa_circ_0006719 as candidate biomarkers of CS, ROC curve analysis was performed. The area under the curve (AUC) was 0.739 (95% CI: 0.611‐0.866, P = 0.001) (Figure 4). Therefore, hsa_circ_0006719 may be a potential diagnostic biomarker for CS.

FIGURE 4

ROC curve analysis of hsa_circ_0006719 in CS patients. The AUC value is 0.739 with P value of 0.001

DISCUSSION

In this study, we performed genome‐wide circRNA sequencing and compared the circRNA expression profiles of peripheral blood from CS patients and healthy controls. Gene Ontology and KEGG pathway analysis were conducted to explore the potential function of the differentially expressed circRNAs. Previous studies have revealed that circRNAs are outputs of protein‐coding genes in eukaryotes. , They could function as post‐transcriptional regulators, including working as microRNA sponge and protein sponge, even some of them could translate into proteins. Thus, it is feasible to predict the functions of circRNAs via their related genes. It is postulated that GO and KEGG pathways are extensively used in predicting gene function and enrichment analysis. , According to GO annotation, we hypothesized that circRNAs might contribute to the development of CS by biological process (BP), intracellular organelle part in cellular component (CC) and transcription coactivator activity in molecular function (MF). Kyoto Encyclopedia of Genes and Genomes pathway analysis also indicated that differentially expressed circRNA‐related genes were enriched in the ubiquitin‐mediated proteolysis signalling pathway, endocytosis signalling pathway and other several pathways. Previous studies indicated that the ubiquitin‐mediated proteolysis signalling pathway is involved in osteogenic differentiation and may play essential roles in disorders manifested with scoliosis. , The endocytosis signalling pathway has also reported to play an pivotal role in the development of scoliosis. , We hypothesized that these differentially expressed circRNAs could regulate ubiquitin‐mediated proteolysis and endocytosis signalling pathway, which could influence the development of CS.

To explore specific circRNAs involved in CS, we compared the expression level in the replication cohort. We selected 22 circRNAs with the most stable and significant different expression for the validation. Seven circRNAs were confirmed, and only hsa_circ_0006719 had higher expression level in the CS group than the HC group with significant difference. Moreover, we conducted ROC curve analysis to evaluate the value of hsa_circ_0006719 for differentiating CS and HC. The results suggested that hsa_circ_0006719 was a potential novel diagnostic biomarker of CS. hsa_circ_0006719 originates from the exon 3 and intron 3 of gene VKORC1 with the spliced sequence length of 568 bases. VKORC1 locates in 16p11.2 encoding vitamin K epoxide reductase (VKOR) complex subunit 1. VKORC1 plays an important role in vitamin K metabolism. It is a catalytic subunit of the VKOR complex which reduces vitamin K 2,3‐epoxide to active vitamin K. , Vitamin K is proved to be essential in bone development, including affecting the function of osteoblasts and osteoclasts. The supplementation of vitamin K has an effect of osteoporosis prevention. Several studies also delineate that mutations of VKORC1 are associated with bone mineral density and osteoporosis. , Thus, we hypothesized that hsa_circ_0006719 could contribute to the aetiology of CS by regulating the expression of VKORC1. Further functional studies are needed to identify this hypothesis.

There are several limitations in this study. First, since all the CS patients were under 18 years old, the age of the healthy controls was not matched with them both in the discovery cohort and in the replication cohort. The results could be biased by the unmatched age. Second, the circRNAs always act as post‐transcriptional regulators in target tissues. In our study, the circRNAs in plasma and tissues were not detected, and further studies are needed to explore the underlying mechanism. Third, the sample size of the study was relatively small. Further studies with larger sample size are necessary to validate our findings.

CONCLUSION

In conclusion, we described the first study of circRNAs in CS patients. After genome‐wide circRNAs analysis and validation in a replication cohort, we found that the differentially expressed circRNA‐related genes were enriched in several pathways, and we also found that hsa_circ_0006719 had a higher expression level in CS and could act as a potential novel diagnostic biomarker of CS.

CONFLICT OF INTEREST

The authors declare no conflict of interest.

AUTHOR CONTRIBUTION

GL and JS performed the research, analysed and interpreted the data. GL and JS drafted the manuscript. CC, YJ, ZL, HT, YL and TR helped sample collection. GL and CC performed helped analysis and interpretation of the data. HT provided technique support. JS offered professional discussions and instructions. CC, YL and TR helped bioinformatic analyses. JS conceived and designed the study, revised the manuscript and provided final approval of the manuscript.

Table S1

Click here for additional data file.

Table S2

Click here for additional data file.

Table S3

Click here for additional data file.

Table S4

Click here for additional data file.