Expression Profiles of Circular RNAs in Human Papillary Thyroid Carcinoma Based on RNA Deep Sequencing.

BACKGROUND: Papillary thyroid carcinoma (PTC) is the most prevalent type of thyroid cancer. Herein, we purposed to explore the expression patterns of circRNAs in PTC with the overarching goal of improving early diagnosis rates for individuals with PTC.
METHODS: We used RNA deep sequencing to determine the expression patterns of circRNAs in PTC. Besides, RT-qPCR was employed to confirm circRNAs. The diagnostic potential of the circRNAs was explored by constructing ROC curves. GO along with KEGG pathway analyses were utilized to elucidate the potential biological roles of differentially expressed circRNAs. Moreover, we predicted cross talks among circRNAs, miRNAs, and mRNAs, followed by establishment of a ceRNA network.
RESULTS: Deep sequencing of four PTC pairs and neighboring nontumor tissues identified 16569 circRNAs, of which, 301 were upregulated and 419 were downregulated. The RT-qPCR data demonstrated that the expression of chr5: 38481299-38530666-, chr2: 159932176-159945082-, chr10: 179994-249088+, chr3: 121378716-121381532+, and chr1: 237423092-237445522+ was downregulated, while the expression of chr4: 25665378-25667298+, chr5: 161330883-161336769-, chr1: 12578718-12579412-, chr7: 116695750-116700284+, and chr7: 116699071-116700284+ was upregulated. The stability test exhibited that circRNAs were more tolerant to temperature, RNase R, and time. On the other hand, ROC curves illustrated that chr4: 25665378-25667298+, chr1: 12578718-12579412-, chr7: 116699071-116700284+, chr7: 116695750-116700284+, chr5: 161330883-161336769-, and chr10: 179994-249088+ were effective as diagnostic indicators. However, a logistic regression model combining the six indicators achieved a better combined prediction index, with 97.7% sensitivity and 95.3% specificity. Moreover, GO along with KEGG pathway analyses illustrated that differentially expressed circRNAs were linked to tumorigenesis. Furthermore, bioinformatics analyses established a promising ceRNAs network among mRNAs, circRNAs, and miRNAs.
CONCLUSION: Herein, we demonstrated that several circRNAs are promising PTC diagnostic biomarkers. Further study on the functions and mechanisms of these circRNAs may contribute to the understanding of PTC.

Introduction

The incidence rate of thyroid cancer is increasing annually, with the major histopathologic type still being papillary thyroid carcinoma (PTC). Most PTCs have an excellent prognosis, with 10-year survival rates that can reach 95%.1 However, survival is mainly affected by extrathyroid infiltration, lymph node metastasis, and distant metastasis due to advanced disease, which are the main causes of death.2 Currently, the proportion of microcarcinoma detected using ultrasound is rising. Follow-up studies have shown that some thyroid cancer nodules progress slowly and can be observed with dynamic follow-up.3 However, the tumorigenesis of PTC has not yet been fully elucidated, alongside other microcarcinomas which also lead to severe symptoms.4 At the molecular biology level, these microcarcinomas may still have similar characteristics with the large ones. Therefore, this calls for the identification of specific indicators, and design of specific diagnostic and therapeutic regimens that target these features.

Previous studies have identified several genetic indicators associated with thyroid cancer diagnosis and prognosis including BRAF and TERT mutations.5,6 However, it is possible that there are many other ways through which tumor development is regulated at the RNA level in addition to DNA mutations.7 Non coding RNAs belong to an important family of RNAs, and recent investigations have shown that miRNAs and lncRNAs have many modulatory functions in cancer development.

Recently, circRNAs (circular RNAs) have been rediscovered as a novel class of RNA, which initially were recognized as “scrambled” exons in the early 1990s.8 The special circular structure of circRNAs, closed circular sequences, is more stable than linear RNAs, which have 5ʹ and 3ʹ ends.9 Moreover, longer circRNAs may also locally appear as double stranded RNAs in order to increase their stability.10 This stable structure implies a more stable regulatory role in tumorigenesis. Current functional studies on circRNAs are more focused on adsorbing miRNAs as sponges which influence expression of mRNAs.11,12 circRNAs resistance to exonucleases makes them relatively stable in cells, and thus they may constitute ideal biomarkers for cancer diagnosis and can even act as liquid biopsy biomarkers for human diseases.13 Researchers should identify novel early diagnostic biomarkers and elucidate the responsible molecular mechanism of PTC with the overarching goal of improving early detection rates and reducing mortality in advanced PTC patients.

Currently, RNA deep sequencing is progressively replacing microarrays because it has a greater dynamic range relative to microarrays.14 Accumulating research evidence has opined that circRNAs are usually aberrantly expressed in human cancers, and participate in oncogenesis via multiple approaches.15 Herein, we used RNA deep sequencing (RNA-seq) to identify the circRNAs in PTC and verified differentially expressed circRNAs in PTC samples using RT-qPCR. Furthermore, we employed bioinformatics analyses to predict the potential role of circRNAs in PTC. Our data will provide relevant reference for future studies.

Materials and Methods

Patient Tissue Samples and Cell Lines

Here, 99 fresh PTC and neighboring non-tumor tissues were collected from patients treated with thyroid surgery at the First Hospital of China Medical University between 2018 and 2019. All tissue samples were washed with saline and immediately frozen in liquid nitrogen until later use. PTC diagnosis was confirmed by two pathologists postoperative. The TNM (tumor‐node‐metastasis) staging system (8th ed.) was employed to classify tumors. High-throughput sequencing was performed on first four pairs of tissue samples (Table 1). The remaining 95 patients were randomly divided into the discovery cohort and the validation cohort according to the ratio of 1:2. (Table 2)

Table 1

Clinic Pathological Features of Patients Whose Samples Were Obtained for Sequencing

Patient No.	Sex	Age	TNM (AJCC, 8th, ed., 2017)
1	Female	41	T4aN1bM0, I
2	Male	49	T3bN1aM0, I
3	Male	54	T3bN1aM0, I
4	Female	53	T2N1aM0, I

Table 2

Clinic Pathological Features of Patients Whose Samples Were Obtained for qRT-PCR

Characteristics	Discovery Cohort (n)	Validation Cohort (n)	P-value
Total	32	63
Sex
Male	12	22	0.804
Female	20	41
Age
<55	24	42	0.404
≥55	8	21
Extra thyroidal extension
Yes	6	13	0.828
No	26	50
Primary tumor
≤1cm	11	23	0.837
>1cm	21	40
Lymph node metastasis
Yes	24	39	0.201
No	8	24
TNM stage
I–II	28	52	0.531
III–IV	4	11

The TPC1 cell line was generously provided by Professor Meiping Shen (Department of General Surgery, The First Affiliated Hospital of Nanjing Medical University, Nanjing, Jiangsu).

RNA Extraction and Sequencing

The TriZol reagent (Cat No. 15596018, Life Technologies, Carlsbad, CA, USA) was employed to isolate total RNA from the tissues as described by the manufacturer. RNA quality was checked through determination of 260/280 OD values on a NanoDrop ND-2000 instrument (Thermo Fisher Scientific, Waltham, MA, USA). Total RNAs of four paired PTC and neighboring non-malignant tissues were treated with VAHTS Total RNA-seq (H/M/R) Library Prep Kit for Illumina as described in the manufacturer provided manual. Subsequently, sequencing of the products was run on the Illumina HiSeqTM 2500 platform by Gene Denovo Biotechnology Co. (Guangzhou, China). The processed clean reads were then aligned to the human genome (version: hg38_GRCh38).

RT-qPCR Validation

The remaining 95 sample pairs were employed for confirmation via RT-qPCR. Concisely, the TriZol reagent (Cat No. 15596018, Life Technologies, Carlsbad, CA, USA) was employed to isolate total RNA from the tissues as described by the manufacturer. After that, cDNA was regenerated from the total RNA by using PrimeScript RT Master Mix (Cat No. RR036A, Perfect Real Time, TaKaRa) with random primers. Thereafter, SYBR Premix Ex Taq II (Cat No. RR820A, Tli RHaseH Plus, TaKaRa) was employed to prepare a qPCR reaction mixture and the amplification was run on the LightCycler 480 system (Roche, Basel, Switzerland). We used the following RT-qPCR protocol: denaturation at 95°C for 30 seconds, followed by 50 cycles of: 95°C for 5 seconds and 60°C for 30 seconds. GAPDH served as the normalization standard. The ΔCt method was employed to assess relative RNA expression. Primers for RT-qPCR analysis were designed to target the circRNA specificity junction site, and the specificity of primers was tested using CircPrimer v1.2.0.5.16 The sequences of all primers utilized in RT-qPCR assays as Table 3. The diagnostic potential of the circRNAs for PTC patients was explored through the construction of ROC (Receiver operating characteristic) curves.

Table 3

PCR Primer Sequences for circRNAs

circRNA	Primer F (5ʹ-3ʹ)	Primer R (5ʹ-3ʹ)
chr5: 38481299–38530666-	AAGGCTAATACCTACTTGACTGACTGC	CCACTGGAAATTTGAAGCAGTCCT
chr2: 159932176–159945082-	GAAACTGGTGCTCAAGATGAAACCAT	GGGATCCTCAAGGCAACACAATG
chr10: 179994–249088+	CAGACCAAGAAGAAGCAGTGGCTAA	TGGCAATCTGCTTCTGGTTCCGTATA
chr3: 121378716–121381532+	AGTCATGGTATTGCCAAGTGACCA	GGTCTTCCAGGACCCCAGC
chr1: 237423092–237445522+	CGAGATCCCAGCATGAAGAATCAC	TCTTATGTGGCTTCCACTCCACCTTA
chr4: 25665378–25667298+	GATATCCAGCACAGAAATGATGCTTCC	AGACCAGCCTCAGGTGAGC
chr5: 161330883–161336769-	GGATGTCAACAAGATGGATACACAACT	AGAGAACTGTGGAAGTTCAATTTTCG
chr1: 12578718–12579412-	CATGAGAGTCAGGTGGGTGACAT	AGGGCATCATCATCACTGTCCAT
chr7: 116695750–116700284+	TGCTTTAATAGGGGTGGTGATGAAGAG	GCACGAGGATGCCAGGT
chr7: 116699071–116700284+	AGGATAAACCTCTCATAATGAAGGCCC	CTCTTTACACTCCCCATTGCTCCT

Gene Ontology and Kyoto Encyclopedia of Genes and Genomes Pathway Analysis

Paired-end reads were harvested from an Illumina HiSeqTM2500 sequencer and Bowtie2 v2.4.217 was employed to map the reads to the ribosome RNA (rRNA) data resource. Next, high-quality reads were aligned to the reference genome using Hisat2 v2.2.1.18 We then collected the unmapped reads for circRNA identification using Find_circ v1,19 followed by blasting against circBase20 for annotation. The circRNAs that could not be annotated were defined as novel circRNAs. Moreover, circRNAs that exhibited fold changes ≥ 2 along with adjusted P< 0.05 were considered to be remarkably expressed differentially.

GO along with KEGG pathway analyses were then employed to identify the source genes for the differentially expressed circRNAs.

circRNA-miRNA-mRNA Network Construction

We predicted the miRNAs that might target cirRNAs. For circRNAs that have been annotated in circBase, we used StarBase (v2.0) to determine the target relationship with miRNAs. On the other hand, we used three softwares; Mireap, Miranda (v3.3a), as well as TargetScan (Version:7.0) for novel circRNAs. miRTarBase v8.021 was employed to predict circRNAs/mRNAs targeted by the miRNAs sponge and the resulting correlation of circRNAs-miRNAs-mRNAs was visualized using Cytoscape v3.7.2.22

Statistical Analyses

Statistical analyses were implemented in the GraphPad Prism 6 software (La Jolla, CA) and SPSS 19.0 (IBM, Chicago, IL, USA). P<0.05 signified statistical significance. The Student’s t-test was employed to explore the differences between groups where appropriate.

Compliance with Ethical Standards

We obtained written informed consent for the experimental use of surgical tissue samples from all participants, and the study was approved by the ethics committee of the First Hospital of China Medical University, Shenyang, China.

Results

Overview of circRNA Profiles in PTC

High-throughput RNA sequencing was performed in four pairs of PTC and neighboring nontumor tissues. In total, 16569 circRNAs were identified (Figure 1A), with most of them coming from exons (Figure 1B). For consistency, we named these circRNAs according to the location of their source chromosome.

Figure 1

(A) Distribution of identified circRNAs on the chromosome; (B) Type of circRNAs; (C) Hierarchical cluster analysis of differentially expressed circRNAs; (D) Volcano plot of differential expressed circRNAs in PTC.

On the basis of the screening criteria (adjusted P values ≤ 0.05 along with fold change ≥ 2), we screened out 301 upregulated and 419 downregulated circRNAs. Figure 1C shows the heat-map of all differentially expressed circRNAs, while Table 4 summarizes the genetic information of the top 20 differentially expressed circRNAs. Moreover, the volcano plot exhibits the statistical significance of differentially expressed circRNAs between PTC and neighboring nontumor tissues (Figure 1D).

Table 4

Top 20 Up-Regulated and Down-Regulated Expressed circRNAs

CircRNA ID	Log2(FC)	P-value	Regulation	Length	Type	Source Gene	Circbase ID
chr4: 25665378–25667298+	18.2542369	4.63E-09	Up	1267	Exon_intron	SLC34A2	Novel
chr5: 161330883–161336769-	6.438140432	4.38E-14	Up	536	Exons	GABRB2	hsa_circ_0128482
chr1: 12578718–12579412-	5.087388418	2.79E-09	Up	359	Exons	DHRS3	hsa_circ_0010023
chr7: 116695750–116700284+	5.078815027	0.0000162	Up	1264	Exons	MET	hsa_circ_0002699
chr7: 116699071–116700284+	4.854229053	5.08E-12	Up	1214	One_exon	MET	hsa_circ_0082002
chr4: 119594170–119607297-	4.354257892	0.0000222	Up	3032	Exon_intron	PDE5A	hsa_circ_0125184
chrX: 50595698–50608184-	4.279529272	0.00000361	Up	2522	Exons	SHROOM4	Novel
chr6: 48008634–48009114-	4.002030905	0.0000107	Up	481	One_exon	PTCHD4	hsa_circ_0076710
chr8: 18799295–18804898-	3.130162035	0.00000338	Up	643	Exons	PSD3	hsa_circ_0004458
chr15: 94356137–94385525+	2.885963915	5.97E-11	Up	783	Exons	MCTP2	Novel
chr8: 18765449–18804898-	2.790965789	0.000000218	Up	908	Exons	PSD3	hsa_circ_0002111
chr15: 85664563–85693451+	2.500143534	0.00000539	Up	1450	Exons	AKAP13	Novel
chr15: 94356137–94402019+	2.431752776	0.000000321	Up	1080	Exons	MCTP2	hsa_circ_000609
chr17: 20204333–20205912+	2.346145666	0.0000154	Up	1580	One_exon	SPECC1	hsa_circ_000013
chr1: 232460887–232471528-	2.230642825	0.000115831	Up	1235	Exons	SIPA1L2	Novel
chr18: 8076455–8143779+	2.040147551	0.0000291	Up	934	Exons	PTPRM	hsa_circ_0007144
chr15: 94384022–94402019+	2.031631614	0.000161259	Up	503	Exons	MCTP2	hsa_circ_0037010
chr5: 73058433–73077493+	1.943054875	0.000118498	Up	594	Exons	FCHO2	hsa_circ_0003571
chr21: 15014344–15043574-	1.666373052	0.000129326	Up	203	Exons	NRIP1	hsa_circ_0004771
chr5: 124701014–124701269-	1.489592426	0.00000313	Up	256	One_exon	ZNF608	hsa_circ_000005
chr8: 91351810–91394039+	−4.869558703	5.19E-08	Down	874	Exons	SLC26A7	Novel
chr8: 91348299–91394039+	−4.975266108	0.00000211	Down	4385	Exons	SLC26A7	Novel
chr8: 109412265–109427156+	−4.744834547	0.00000525	Down	915	Exons	PKHD1L1	Novel
chr8: 109409865–109420690+	−4.937439877	0.000017	Down	726	Exons	PKHD1L1	Novel
chr7: 152194440–152224179-	−4.54916372	0.0000619	Down	1349	Exons	KMT2C	Novel
chr7: 137571175–137587210-	−4.53873155	0.0000107	Down	636	Exons	DGKI	Novel
chr6: 154199668–154214276-	−4.546343115	1.85E-08	Down	518	Exons	IPCEF1	hsa_circ_0078336
chr5: 38510464–38530666-	−4.120223175	0.00000208	Down	1010	Exons	LIFR	Novel
chr5: 38481299–38530666-	−18.05460834	0.00000014	Down	3752	Exons	LIFR	hsa_circ_0129039
chr3: 89209860–89210520+	−4.562945635	0.0000851	Down	661	One_exon	EPHA3	hsa_circ_0066596
chr3: 121378716–121381532+	−16.77573593	0.0000165	Down	411	Exons	STXBP5L	hsa_circ_0121603
chr2: 159967539–159987358-	−5.607787748	1.34E-09	Down	1070	Exons	PLA2R1	Novel
chr2: 159932176–159945082-	−17.06211207	0.00000249	Down	10027	Exons	PLA2R1	Novel
chr18: 33626862–33683568+	−4.581834619	0.0000598	Down	995	Exons	ASXL3	Novel
chr12: 98773042–98801125-	−4.39889169	0.000139626	Down	471	Exons	ANKS1B	Novel
chr11: 30535994–30580494-	−5.181240634	0.000226299	Down	431	Exons	MPPED2	hsa_circ_0021553
chr11: 30495296–30580494-	−4.429826268	0.0000192	Down	755	Exons	MPPED2	hsa_circ_0021550
chr10: 4830675–4847299+	−4.451248834	0.000028	Down	1400	Exons	AKR1E2	hsa_circ_0017523
chr10: 179994–249088+	−17.05639739	0.00000379	Down	1975	Exons	ZMYND11	hsa_circ_0093332
chr1: 237423092–237445522+	−16.39703526	0.000199672	Down	444	Exons	RYR2	hsa_circ_0112640

Verification of circRNA Expression

The ten most differentially expressed circRNAs were selected to verify the sequencing quality in the discovery data set using qRT-PCR. Data illustrated that the expression of chr5: 38481299–38530666-, chr2: 159932176–159945082-, chr10: 179994–249088+, chr3: 121378716–121381532+, and chr1: 237423092–237445522+ was downregulated, while the expression of chr4: 25665378–25667298+, chr5: 161330883–161336769-, chr1: 12578718–12579412-, chr7: 116695750–116700284+, and chr7: 116699071–116700284+ was upregulated (Figure 2A–C).

Figure 2

Ten circRNAs were selected to verify the sequencing quality by qRT-PCR. (A) Comparison of log2 fold changes in circRNAs between RNA-Seq and qRT-PCR results; (B and C) The relative expression levels of circRNAs in 32 PTC tissues and paired normal thyroid tissues.

Diagnostic Significance of RT-qPCR Verified circRNAs

To further explore the diagnostic capacity of circRNAs, we chose chr7: 116699071–116700284+ which had the highest expression in PTC. After quantification, we extracted fresh RNA from TPC1 cell line and divided it into six samples. Next, three samples were inoculated with RNase R according to the manufacturer’s instructions, and placed in an enclosed environment without RNA enzyme pollution at 4°C, 25°C, and 37°C, respectively. After 1 h, 6 h, and 12 h, the same volume of sample RNA was taken for reverse transcription with random primers. The expression levels of circRNA and source gene mRNA were then determined using qRT-PCR. Results indicated that circRNAs were more stable than mRNAs under the influence of temperature, RNase R, and time (Figure 3).

Figure 3

The stability of circRNAs (A) and its source mRNA (B) in cell line TPC1.

We then developed ROC curves for RT-qPCR verified circRNAs using discovery data set (n = 32, Figure 4A) and verification data set (n = 63, Figure 4B). Our results indicated that the ten circRNAs all showed potential diagnostic capabilities for diagnosing PTC, with consistent areas under the curve (AUC) in both the verification and discovery data sets (Table 5).

Figure 4

ROC curves qRT-PCR verified circRNAs: (A) Discovery cohort (n = 32); (B) Validation cohort (n = 63); (C) The ROC of Combined Prediction Index.

Table 5

Area Under the Curve (AUC) of qRT-PCR Verified circRNAs in the Discovery and Validation Cohort

	Discovery Cohort			Validation Cohort
circRNA	AUC	95% CI	P	AUC	Se(%)	Sp(%)	Youden	95% CI	P
chr4: 25665378–25667298+	0.893	0.787–0.999	<0.0001	0.889	92.5	80	72.5	0.812–0.967	<0.0001
chr1: 12578718–12579412-	0.901	0.792–1.000	<0.0001	0.861	82.5	80	62.5	0.774–0.945	<0.0001
chr7: 116699071–116700284+	0.966	0.901–1.000	<0.0001	0.941	95	92.5	87.5	0.876–1.000	<0.0001
chr7: 116695750–116700284+	0.931	0.835–1.000	<0.0001	0.919	97.5	85	82.5	0.849–0.989	<0.0001
chr5: 161330883–161336769-	0.932	0.837–1.000	<0.0001	0.924	90	90	80	0.852–0.994	<0.0001
chr5: 38481299–38530666-	0.741	0.579–0.904	0.007	0.788	62.5	92.5	55	0.683–0.894	<0.0001
chr10: 179994–249088+	0.961	0.912–1.000	<0.0001	0.856	82.5	82.5	65	0.731–0.932	<0.0001
chr3: 121378716–121381532+	0.755	0.599–0.911	0.005	0.775	67.5	85	52.5	0.668–0.882	<0.0001
chr1: 237423092–237445522+	0.669	0.502–0.836	0.061	0.706	52.5	82.5	35	0.593–0.821	0.001
chr2: 159932176–159945082-	0.658	0.486–0.829	0.081	0.663	47.5	92.5	40	0.513–0.764	0.012

Furthermore, a combined prediction index was established based on a logistic regression model which combined the expression levels of the most efficient circRNAs (AUC > 0.8, sensitivity > 80%, and specificity > 80%). The coefficients in the prediction equation were as follows: Combined Prediction Index =−0.321× (chr4: 25665378–25667298+)+0.23×(chr1: 12578718–12579412-) −2.818× (chr7: 116699071–116700284+)-1.078× (chr7: 116695750–116700284+)+0.689× (chr5: 161330883–161336769-)−0.723× (chr10: 179994–249088+)+15.627.

The results showed that the performance of AUC, sensitivity, and specificity was greater when the Combined Prediction Index was used as a diagnostic marker in PTC. Moreover, the Combined Prediction Index showed the best diagnosis efficiency compared to the single markers (Figure 4C).

GO Enrichment Analysis and Pathway Enrichment Analysis

GO annotation reveals the source genes of the differentially expressed circRNAs. The most dramatically enriched GO terms with regard to the biological process, cellular component, and molecular function categories were modulation of small GTPase mediated signal transduction (GO: 0051056, gene count = 27, P = 5.75E-07), organelle (GO: 0043226, gene count = 392, P = 1.14E-04), and adenyl nucleotide-binding (GO: 0030554, gene count = 73, P = 2.72E-06), respectively (Figure 5).

Figure 5

GO terms analysis for differential expressed host gene of circRNAs.

On the other hand, KEGG pathway analysis demonstrated that the differentially expressed circRNA genes were primarily associated with axon guidance, lysine degradation, bacterial invasion of epithelial cells, adherens junction, sphingolipid signaling cascade, choline metabolism in cancer, cellular senescence, thyroid hormone signaling pathway, tight junction, and autoimmune thyroid disease (Figure 6).

Figure 6

KEGG pathways analysis for differential expressed host gene of circRNAs.

In addition, miRNAs sponge has been widely studied in order to determine the function of circRNAs. The top five miRNAs were selected for each circRNAs, the top ten mRNAs were selected for each miRNAs. We constructed the predicted circRNA-miRNA-mRNA network for the top three upregulated and downregulated circRNAs in PTC (Figure 7).

Figure 7

The potential competing endogenous RNA (ceRNA) relationships in PTC. Yellow represents circRNAs, green represents miRNAs, and blue represents mRNAs.

Discussion

Circular RNA is a class of endogenous RNA that widely exists in eukaryotes,23 and was once considered to be the product of wrong splicing. However, with the development of research, the characteristics and functions of circular RNA are becoming more and more clear.

CircRNAs may arise from exons or introns of pre‐mRNA after back‐splicing.24 Interestingly, they regulate gene expression and biological processes, serve as miRNA sponge, regulate transcription, are involved in rolling circle translation, generate pseudogenes, and affect alternative splicing.25 The progress of high-throughput sequencing technology in recent years has resulted in identification of several circRNAs. Accumulating evidence has opined that circRNAs play a core role in the onset of cancers and may be applied as novel biomarkers.26–29 It is worth noting that they are few studies on thyroid cancer, and thus it is essential to explore the expression of circRNA in PTC.

Herein, we chose ten differentially expressed circRNAs for RT-qPCR verification. The obtained results were consistent with the reliability of RNA-Seq. chr4: 25665378–25667298+ was the most up-regulated circRNA, while chr5: 38481299–38530666 was the most down-regulated circRNA in PTC. We then tested the stability of circRNA and its source gene mRNA in PTC cell line. Results indicated that circRNAs can tolerate the effects of RNase R treatment and temperature stimulation, and time extension demonstrated that the stability of circRNAs was stronger than that of linear RNAs. Next, we plotted ROC curves using the data verified by RT-qPCR as the discovery data set. ROC curves illustrated that the verified circRNAs are promising PTC diagnostic biomarkers. Moreover, we validated the ROC curves using 63 newly paired samples of cancer and neighboring nontumor tissues. The data obtained from the verification data set also showed that the circRNAs could stably predict PTC. To achieve better diagnostic ability, we constructed a Combined Prediction Index using a logistic regression model which combined the six most efficient circRNAs. Results indicated that the Combined Prediction Index achieved the highest AUC, sensitivity, and specificity compared to the original six indicators. With regard to verification of stability, we believe that the above-mentioned circRNAs can all be used as new biomarkers for PTC identification. However, combined use of multiple circRNAs has the highest diagnostic efficiency.

CircRNAs are derived from alternative splicing of RNA. Most of the circRNAs sequences are similar to their parent genes and their functions may also be similar. Moreover, a previous research documented that circRNAs can modulate the transcription of parental genes.30 To further explore the role of circRNAs, we identified differentially expressed circRNAs among groups and analyzed the function of source genes through GO along with KEGG pathway analyses. Our results illustrated that differentially expressed circRNAs were associated with tumorigenesis including adherens junction, cellular senescence, and tight junction. The most commonly reported function of circRNAs is their action as miRNA sponges.31 One study reported that circRNAs are involved in the pathogenesis of PTC through modulating miRNAs and their target mRNA.32 Herein, we predicted the possible miRNA targets of circRNAs and developed a circRNA-miRNA-mRNA network for determining the regulatory role of circRNAs in PTC. Our network indicated that some miRNAs have been identified in other tumors. For instance, hsa-miR-1299 regulates breast cancer,33 prostate cancer,34 hepatocellular carcinoma,35 and cholangiocarcinoma.36 In addition, Hsa-mir-22-3p regulates lung adenocarcinoma progression,37 while Wnt9A activation is related to repression of human colorectal cancer cell proliferation.38 Furthermore, FOXA2 suppresses the metastasis of hepatocellular carcinoma,39 and ZEB1 activation induces epithelial-mesenchymal transition.40 Therefore, future research should focus on these miRNAs and mRNAs identified in our network.

However, this study has some limitations. Firstly, we only selected four pairs of cancer tissues and nontumor tissues for RNA sequence analysis, which may lead to bias. Secondly, there is a long way to go before circRNAs can be used in the diagnosis of PTC.

In conclusion, this study has exhibited the circRNA expression profiles in PTC. To the best of our knowledge, this is the first study that has shown that chr4: 25665378–25667298+, chr1: 12578718–12579412-, chr7: 116699071–116700284+, chr7: 116695750–116700284+, chr5: 161330883–161336769-, and chr10: 179994–249088+ may be potential biomarkers for the diagnosis of PTC. Further studies on the functions and mechanisms of these circRNAs will contribute to our understanding of PTC tumorigenesis.