This article is mentioned in:

Introduction

Thyroid cancer (TC) is the commonest endocrine malignancy worldwide, accounting for 1–5% of all cancers in females and <2% of all cancers in males. TC incidence continues to increase globally and is neither confined to a particular region nor influenced by underlying rates of TC (1). Papillary thyroid carcinoma (PTC) is the most frequent type of differentiated thyroid cancer, accounting for 80% of all pathological types of TC (2,3). Studies on the incidence trends for TC subtypes and data from volumes 4–9 of Cancer Incidence Database on five continents (CI5) show that the incidence of PTC in both sexes has risen steadily in countries around the world over the past few decades (4–10). As indicated by a study of PTC in Manitoba, Canada, the age-standardized rate of PTC increased from 0.93/100,000 to 6.68/100,000 individuals per year over the period 1970–2010, showing a significant increase (11). At present, the main treatment for PTC is a combination of surgery, radioactive iodine and levothyroxine suppression, with a 50% 5-year survival rate (12–15). Although PTC is less malignant than other cancer types, the incidence and recurrence rates are still high (20–40%). The reasons for the apparent increase in PTC incidence are unclear and overdiagnosis and a small but actual increase in PTC incidence are indicated as the two main processes in some studies (16,17). Considering the expanding demands for targeted therapies, the role of biomarkers in PTC diagnosis and treatment remains to be elucidated. BRAF mutation is proposed to be a prognostic biomarker; however, controversy exists regarding the high positive ratio of BRAF mutation combined with the low ratio of PTC aggressive biological behaviors (18–20). Recently, non-coding RNAs and their processing machinery have been shown to be common hallmarks of cancer regulating tumorigenesis, progression and metastasis via various mechanisms (21,22). Micro (mi)RNA dysregulation has been discovered in numerous cancer types, including miRNA (miR)15a/16-1 cluster dysregulation in chronic lymphocytic leukemia, miR-127 in primary prostate cancer and bladder tumors and miR-29 family in lung cancer (23,24). Significant advances have been made regarding the role of miRNAs in diagnostics, monitoring and therapy. For example, the expression levels of miR-145 are inversely correlated with BCR-ABL rearrangement levels during diagnosis and treatment of chronic myeloid leukemia, which can be used to modify the clinical treatment for improved outcomes (25). As potential oncogenes or oncosuppressor genes, miRNAs can be used as anticancer therapies. For instance, miR-203 inhibits the spheroid formation and cancer stem cells marker expression by targeting SOCS3, indicating a novel miRNA-based clinical treatment for estrogen receptor-positive breast cancer (26–28). In contrast to miRNAs, circular RNAs (circRNAs) are a special type of non-coding RNA molecules that are produced from pre-mRNA by backsplicing mechanisms. Unlike traditional linear counterparts (containing 5′ and 3′ ends), circRNAs have a closed circular structure, are not affected by RNA exonucleases in eukaryotes and exhibit distinct advantages in their stability over canonical linear RNAs (29,30). Without free ends, circRNAs are expressed in a cell type-specific manner, display high stability and function in gene regulation in various pathologies. Recent studies have indicated that circRNAs can serve as miRNA sponges to regulate gene expression and sequester or recruit RNA/proteins; they even encode peptides or proteins (31,32).

Increasing evidence shows that circRNAs are stably expressed in exosomes, plasma and saliva, indicating their high capacity for detection and clinical utility in cancer (33,34). Evidence suggests that circRNA profiles are cancer specific and related to tumor occurrence and progression (35–37). The present study screened circRNA expression in PTC tumors and verified the dysregulated circRNAs in whole blood from patients with PTC to evaluate their potential value as PTC biomarkers. The possible functional associations of these circRNAs with PTC were also analyzed using bioinformatics.

Materials and methods

Ethics

The present study was approved by the Medical Ethics Committee of Harbin Medical University (approval no. IRB3011619). All participants involved in this study were informed of the research objectives and signed consent forms before the study.

Subjects

Patients who underwent thyroidectomy were recruited at The First, Second and Fourth Clinical Hospitals Affiliated with Harbin Medical University between January 2017 and July 2019. Patients with PTC were diagnosed and evaluated by neck ultrasonic imaging and pathological diagnosis. Histopathological analysis and diagnostic immunohistochemistry of all thyroid tissue specimens were independently performed by two licensed pathologists. The immunohistochemistry of NF-κB1 was used to confirm the diagnostic accuracy. Briefly, the morphological and pathological identification of PTC was performed by H&E staining at room temperature for 10 min. Paraformaldehyde-fixed (4°C overnight) and paraffin embedded thyroid tissue slides were dewaxed in xylene and rehydrated in ethanol. The endogenous peroxidase was blocked by 3% H2O2 at room temperature for 10 min. Heated slides in a microwave submersed in 1X citrate unmasking solution (Cell Signaling Technology, Inc.) until boiling was initiated, followed by 10 min at sub-boiling temperature (95–98°C). The slides were cooled on bench top for 30 min. The slides were then incubated with a 1:100 dilution of NF-κB1 primary antibody (cat. no. bs-1194R; Beijing Biosynthesis biotechnology Co., Ltd.) at 4°C overnight. After washing, the slides were incubated with goat anti-rabbit IgG (cat. no. PV-6001; OriGene Technologies, Inc.) for 20 min, followed by staining with DAB for 5 min at room temperature. The slides were counterstained with hematoxylin at room temperature for 20 sec, and then dehydrated with ethanol and xylene. The morphological changes were observed by a light microscope (magnification, ×40 or ×100; BM2000; Nanjing Jiangnan Yongxin Optics Co., Ltd.). Images were captured in three random fields. The number of cells was analyzed by ImageJ software (v1.8.0.112, National Institutes of Health). Ultimately, 57 cases of PTC and age- and sex-matched healthy controls were recruited for this study (Table I). Among these subjects, three pairs of PTC tumors and matching paracancerous tissues were collected for microarray analysis. In total, 57 pairs of whole-blood samples were obtained from all subjects and refrigerated at −80°C for further studies.

Table I. Correlation between hsa_circ_IPCEF1 expression and clinicopathological characteristics in 57 patients with papillary thyroid carcinoma.

Table I.

Correlation between hsa_circ_IPCEF1 expression and clinicopathological characteristics in 57 patients with papillary thyroid carcinoma.

Characteristics	Number of patients (%)	Relative expression	P-value
Age (years)
<45	26 (45.6)	0.75±0.16	0.1933
≥45	31 (54.4)	0.52±0.09
Sex
Male	17 (29.8)	0.76±0.20	0.1571
Female	40 (70.2)	0.45±0.11
Tumor size (cm)
<2	49 (86.0)	0.70±0.14	0.1246
≥2	8 (14.0)	0.17±0.08
TNM stage
I–II	42 (73.7)	0.77±0.20	0.1515
III–IV	15 (26.3)	0.28±0.11
Multifocality
Yes	30 (52.6)	0.61±0.09	0.4660
No	27 (47.4)	0.48±0.16
Extrathyroidal extension
Yes	12 (21.1)	0.69±0.26	0.8503
No	45 (78.9)	0.63±0.18
LNM
Yes	32 (56.1)	0.82±0.13	<0.001a
No	25 (43.9)	0.25±0.05

a P<0.001. circ, circular RNA; TNM, tumor-node-metastasis; LNM, lymph node metastasis.

RNA isolation and quality control

Total RNA was extracted from the tissue samples and whole blood with TRI Reagent BD, which was supplied by Sigma-Aldrich (Merck KGaA) in accordance with the manufacturer's protocol. RNA was precipitated with isopropanol, washed twice with 75% ethanol and resuspended in RNase/DNase-free water. The RNA concentration and purity were assessed with a Nano-200 system (Hangzhou Allsheng Instruments Co., Ltd.). RNA samples with OD260/OD280 ratios ≥1.8 and ≤2.1 and OD260/OD230 ratios >1.8 were used for further studies.

circRNA array

circRNA array analysis was performed by KangChen Bio-tech. Briefly, three pairs of PTC tumors and matching paracancerous tissues were used for Arraystar Human circRNA Array analysis. The linear RNAs were removed by RNase R (Epicenter Biotechnologies) to enrich circRNAs. The enriched circRNAs were amplified and transcribed into fluorescent complementary RNA (cRNA) by the random priming method [Quick Amp Labeling Kit, One-Color (Agilent p/n 5190-0442)], which was performed by KangChen Bio-tec. Subsequently, these labeled cRNAs were purified by an RNeasy Mini kit (Qiagen GmbH) and were hybridized to an Arraystar human circular RNA array V2 (8×15K, Arraystar, Inc.). The arrays were then scanned with an Agilent G2505C scanner (Agilent Technologies, Inc.) followed by washing. Array images and data analysis were performed using Agilent Feature Extraction software (version 11.0.1.1; Agilent Technologies, Inc.) and the R software package version 3.1.2 (R Foundation for Statistical Computing; http://www.R-project.org/) Differentially expressed circRNAs with P-values <0.05 between the two groups were visualized in a volcano plot and with hierarchical clustering. The relevant datasets have been submitted to National Center for Biotechnology Information Gene Expression Omnibus (GEO). The GEO accession no. is GSE173299 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE173299).

Reverse transcription-quantitative (RT-q) PCR validation

The top 20 most differentially expressed circRNAs from array analysis were further validated by RT-qPCR on 57 pairs of whole blood and three pairs of tissue samples from subjects. Total RNA (1 µg) was reverse-transcribed into cDNA using a First Strand cDNA Synthesis kit (Thermo Fisher Scientific, Inc.) in accordance with the manufacturer's protocol. The candidate circRNAs were amplified by Power SYBR Green PCR Master Mix (Applied Biosystems; Thermo Fisher Scientific, Inc.) with various primers. The primers for each circRNA are presented in Table II. RT-qPCR was performed on an ABI 7300 Real-Time PCR System (Applied Biosystems; Thermo Fisher Scientific, Inc.) in a 10 µl PCR mixture. The thermocycling conditions were as follows: 95°C for 10 min, 40 cycles of 95°C for 15 sec and 60°C for 60 sec. β-actin was used as an internal control. The relative expression of each group was analyzed using the 2−ΔΔCq method (38).

Table II. Primers fortop 20 dysregulated circRNAs.

Table II.

Primers fortop 20 dysregulated circRNAs.

circRNA	Regulation	Forward	Reverse
hsa_circ_0000277	Up	5′CAGTCTTCAAGGTGGGATCGTAA3′	5′TGGAAGGCTTGGATCAGTCAG3′
hsa_circ_0044556	Up	5′CTGGTCCTGATGGCAAAACTG3′	5′GGGGTCCTTGAACACCAACA3′
hsa_circ_0014234	Up	5′CCAGAGCTATGCTTTAGGTCTCA3′	5′AGTGGGAAGTGGGAGGTGTC3′
hsa_circ_0040773	Up	5′AAGTATTACCCCGTCTTTAAGCAG3′	5′TTCCAGACACGCCCATCAC3′
hsa_circ_0074530	Up	5′ATGAGCAGGCACTCCTTGGA3′	5′TCAGTGGCGGGTACACCTTC3′
hsa_circ_0074595	Up	5′GCCTATAAGGAGGACTATCACAAG3′	5′CTGCGGTGCGTGATGATA3′
hsa_circ_0001681	Up	5′AGAGGTGGCATCTGTGAACTGTC3′	5′GGGAAGGCGTATGTTCAAGGTA3′
hsa_circ_0049237	Up	5′CCATAGGCTCACAACACCACA 3′	5′CCCTGCGTGTCCACCTCTA3′
hsa_circ_0058230	Up	5′TGGATGGGGAGCCCTACAAG3′	5′CCAGGTGCGGGTGTACAGG3′
hsa_circ_0057691	Up	5′AGCAACCAAGTGCCAGGAGT3′	5′CGGGTGCATCTGTCACATAACT3′
hsa_circ_FGFR2	Down	5′GAATACGGGTCCATCAATCACA3′	5′ATCACGGCGGCATCTTTC3′
hsa_circ_0095448	Down	5′CCCCTGGAATAACATACAAACC3′	5′GTAGCTGCTCCCGTAAACTGAT3′
hsa_circ_0031968	Down	5′GTTTGGTGTCTCCCCGCTAT3′	5′GCCTTCTGCAACTGGAATCA3′
hsa_circ_IPCEF1	Down	5′AGATAAGCCTGCTGGATCAAAG3′	5′ACAGTGAAATCAGGCAGGTTG3′
hsa_circ_0021549	Down	5′TTCGGAGGTAACAGTGAAGGGA3′	5′AGGCATCTGGATACCATCTGTTCT3′
hsa_circ_0070098	Down	5′GCTCAGTGTCAGCCTCTATTTTG3′	5′GCATTGGTTGGCAGCTATTTG3′
hsa_circ_0079891	Down	5′GAATTGCTTTTGATGCTGAGTCTG3′	5′TTTCCATGAGTTTGGGGTAGG3′
hsa_circ_0001938	Down	5′CATCGTTCTTACAGTTCTGCACA3′	5′CAGCCACGAAGCCAAAGC3′
hsa_circ_0020396	Down	5′GCTTGATCGAAATCGTCCACA3′	5′GAAAGTTCATCCGCTCCTCTG3′
hsa_circ_0021550	Down	5′GATTCGGAGGTAACAGTGAAGG3′	5′TACAGAGCAAAAGATGGAAGCA3′
β-actin		5′GTGGCCGAGGACTTTGATTG3′	5′CCTGTAACAACGCATCTCATATT3′

[i] circRNA/circ, circular RNA.

Circularization of hsa_circ_IPCEF1 validation

To confirm the splicing junction of hsa_circ_IPCEF1, RT-PCR followed by sequencing was performed in whole blood from PTC subjects. Briefly, total RNA was extracted from whole blood from patients with PTC as described above. After reverse transcription, the cDNA was amplified using PrimeSTAR® Max DNA Polymerase (Takara Biotechnology Co., Ltd.). The thermocycling conditions were 35 cycles of 95°C for 10 sec, 55°C for 15 sec and 72°C for 10 sec. Using divergent primers (forward: 5′-GTTTGTCTGCTGCTGAAGATGAG-3′; reverse 5′-CCATCAGCTTTCTCTGCCTTTGTC-3′) to amplify the hsa_circ_IPCEF1 product containing a head-to-tail splicing junction site. The PCR product was purified by 1.5% agarose gel electrophoresis, visualized by ethidium bromide and identified by Sanger sequencing.

Competing endogenous RNA (ceRNA) network analysis

An interaction analysis of circRNA/miRNA was performed by KangChen Bio-tech based on TargetScan (39–41). Through merging the commonly targeted miRNAs, a circRNA-miRNA-mRNA interaction network of hsa_circ_IPCEF1 was constructed by Cytoscape 3.7.2 (https://cytoscape.org). Briefly, by merging the commonly targeted miRNAs, a ceRNA network was constructed via three conditions (42). First, the relative concentration of the ceRNAs and their microRNAs is clearly important; second, the effectiveness of a ceRNA depends on the number of microRNAs that it can ‘sponge’; and third, not all of the MREs on ceRNAs are equal. Therefore, only these ceRNA-pair relations that passed some filtering measures were accepted. In addition to a measure with the number of common microRNAs, a hypergeometric test was executed for each ceRNA pair separately, which was defined by four parameters: i) N is the total number of miRNAs used to predict targets; ii) K is the number of miRNAs that interact with the chosen gene of interest; iii) n is the number of miRNAs that interact with the candidate ceRNA of the chosen gene; and iv) is the common miRNA number between the two genes (43). The test calculates the P-value by using the following formula:

Functional group analysis

The molecular functional roles of the circRNA-target gene profiles were analyzed using Gene Ontology (http://www.geneontology.org), including biological process (BP), cellular component (CC) and molecular function (MF). The P-value produced by topGO (http://www.bioconductor.org/packages/release/bioc/html/topGO.html) denoted the significance of GO term enrichment in the circRNA-miRNA targeted genes. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis was then performed (https://www.genome.jp/kegg/). The P-value (EASE score, Fisher P-value, or hypergeometric P-value) denoted the significance of the pathway relevant to PTC. The P-value cutoff was 0.05. The ceRNAs of each differentially expressed circRNA were further annotated in terms of the diseases and pathways using the KEGG orthology-based annotation system (KOBAS), (version 3.0; http://kobas.cbi.pku.edu.cn/kobas3/?t=1). P<0.05 was considered to indicate a statistically significant difference.

Statistical analysis

All data are expressed as the means ± standard error of the mean. circRNAs with fold changes ≥2 and P-values ≤0.05 were selected as the significantly differentially expressed circRNAs. Receiver operating characteristic (ROC) curves and the area under the curve (AUC) value and 95% confidence intervals (CIs) were established using GraphPad Prism 8.0 (GraphPad Software, Inc.) to predict their potential association with PTC. Comparisons of data were acquired by one-way ANOVA followed by Bonferroni post hoc or Student's t-test. P<0.05 was considered to indicate a statistically significant difference.

Results

Characterization of human thyroid samples

A total of three pairs of PTC tumors and matching paracancerous tissues were collected for microarray. The PTC tumor showed the typical papillary architecture with branching and the enlarged irregular nuclei were oval-shaped and overlapping, showing ground-glass appearance with prominent nuclear grooves and pink cytoplasmic invaginations and intra-nuclear pseudoinclusion (Fig. 1A). NF-κB1-related cancer phenotype (increased expression of NF-κB1 in cell cytoplasm) were also observed in PTC tumors (Fig. 1B and C). Through pathological reports, the tumor tissues were confirmed to be PTC for circRNA chip analysis.

Figure 1. Immunohistochemistry of thyroid specimens. Hematoxylin-eosin stained images of (A) papillary thyroid cancer and cancer contralateral normal thyroid tissues, (B) NF-κB1 staining in papillary thyroid cancer and cancer contralateral normal thyroid tissues, (C) NF-κB1 expression in papillary thyroid cancer and cancer contralateral normal thyroid tissues. Results are described as ratios of NF-κB1 positive cells to all cells. ***P<0.001 (n=3 pairs of samples; ≤6–10 fields from each group were imaged and scored in a blinded manner). Comparisons of data were acquired by one-way ANOVA followed by Bonferroni post hoc. Magnification of ×40 and ×100 from left to right. Nor, PTC tumor-adjacent tissues; PTC, papillary thyroid cancer tissues.

circRNA expression profile in PTC

circRNA profiling was performed with Arraystar Human circRNA Array in three pairs of PTC tumor tissues and cancer contralateral normal thyroid tissues. Hierarchical cluster analysis revealed that the circRNA expression patterns were distinctive between PTC tumor and tumor-adjacent normal thyroid tissues (Fig. 2A). A cut-off value for differential expression fold change was set to 2.0 and 158 significantly dysregulated circRNAs (74 upregulated and 84 downregulated) were found in the PTC tumors (Fig. 2B-D, P<0.05 and false discovery rate <0.05).

Figure 2. Difference analysis of circRNA microarray between PTC tumors and the adjacent tissues. (A) Heatmap of dysregulated circRNAs in two groups. Red represents high expression, while green represents low expression. (B) The volcano plot shows the differential expression of circRNAs in two groups. Red points indicated circRNAs expressed more/less than 2-fold change in PTC patients (P<0.05). (C) The box plot shows the stable distributions between samples. The distributions were nearly the same after normalization. (D) Scatter plots were used to identify differentially expressed circRNAs in two groups. The circRNAs above the top green line and below the bottom green line indicated more than a 2.0-fold change of circRNAs between two compared samples. n=3 pairs of samples. Comparisons of data were acquired by the Student's t-test. circRNA, circular RNA; PTC, papillary thyroid cancer; Thy, PTC tumor tissues. Nor, PTC tumor-adjacent tissues.

Validation of the differential expression of circRNAs in PTC

From microarray results, thetop 20 most dysregulated circRNAs (10 upregulated and 10 downregulated, Table II) were selected for RT-qPCR validation on three pairs of PTC tumor tissues and cancer contralateral normal thyroid tissues as well as in whole blood from 57 pair of patients with PTC and healthy controls. The characteristics of the 57 patients with PTC are presented in Table I. The RT-qPCR results showed that 7 circRNA expression trends were consistent with the microarray assay, which included 3 upregulated circRNAs (hsa_circ_0000277, hsa_circ_0074530 and hsa_circ_0057691) and 4 downregulated circRNAs (hsa_circ_0020396, hsa_circ_0095448, hsa_circ_IPCEF1 and hsa_circ_0021549; Figs. 3 and 4).

Figure 3. Confirmation of the differential expression of circRNAs in PTC tissues by reverse transcription-quantitative PCR. The relative expressions of 12 specific circRNAs between PTC and adjacent tissues were validated. *P<0.05, **P<0.01, ***P<0.001, n=3 pairs of samples. Comparisons of data were acquired by the Student's t-test. circRNA, circular RNA; PTC, papillary thyroid cancer.

Figure 4. Confirmation of the differential expression of circRNAs in whole blood from PTC patients and healthy controls. (A) hsa_circ_0000277, (B) hsa_circ_0074530, (C) hsa_circ_0057691 were significantly upregulated, while (D) hsa_circ_0020396, (E) hsa_circ_0095448, (F) hsa_circ_IPCEF1 and (G) hsa_circ_0021549 were significantly downregulated in PTC patients. **P<0.01, ***P<0.001, n=57 pairs of subjects. Comparisons of data were acquired by one-way ANOVA followed by Bonferroni post hoc. circRNA/circ, circular RNA; PTC, papillary thyroid cancer; Ctl, healthy controls.

Functional annotation of the dysregulated circRNAs

By targeting miRNA response elements (MREs), circRNAs can regulate mRNA by acting as miRNA sponges. The predicted MREs for these seven validated circRNAs are presented in Table III. Based on TargetScan, the targets/miRNAs were analyzed to identify the potential targets of miRNAs (KangChen Bio-Tech). Based on the seven dysregulated circRNAs and their predicted MREs, the 256 predicted ceRNA genes are listed in Table SI.

Table III. Predicted MREs of seven validated circRNAs.

Table III.

Predicted MREs of seven validated circRNAs.

circRNA	FC (abs)	chrom	GeneSymbol	MRE1	MRE2	MRE3	MRE4	MRE5
hsa_circ_0057691	2.81	chr2	SATB2	hsa-miR-9-5p	hsa-miR-1264	hsa-miR-628-5p	hsa-miR-514a-5p	hsa-miR-30c-1-3p
hsa_circ_0074530	2.92	chr5	CD74	hsa-miR-4677-5p	hsa-miR-7161-3p	hsa-miR-515-3p	hsa-miR-1271-3p	hsa-miR-3663-5p
hsa_circ_0000277	2.55	chr11	PDE3B	hsa-miR-138-2-3p	hsa-miR-873-5p	hsa-miR-571	hsa-miR-205-3p	hsa-miR-421
hsa_circ_0020396	2.68	chr10	DOCK1	hsa-miR-103a-2-5p	hsa-miR-132-3p	hsa-miR-196a-5p	hsa-miR-138-2-3p	hsa-miR-136-5p
hsa_circ_0095448	3.02	chr11	SOX6	hsa-miR-1228-5p	hsa-miR-367-5p	hsa-miR-6501-5p	hsa-miR-3692-5p	hsa-miR-22-3p
hsa_circ_IPCEF1	2.69	chr6	IPCEF1	hsa-miR-3916	hsa-miR-3619-5p	hsa-miR-4314	hsa-miR-5589-5p	hsa-miR-6837-3p
hsa_circ_0021549	2.61	chr11	MPPED2	hsa-miR-574-5p	hsa-miR-452-5p	hsa-miR-145-3p	hsa-miR-29b-1-5p	hsa-miR-553

[i] MRE, miRNA response elements; miRNA, microRNA; circRNA/circ, circular RNA; FC, fold change; chrom, chromosome.

The molecular functional roles of the circRNA-target gene profiles were analyzed using GO and KEGG analyses. For the seven dysregulated circRNAs, the top 10 enriched GO terms are shown and ranked by Enrichment Score [-log10 (P-value)]. Through GO analysis of the different genes of dysregulated circRNAs, it is possible to identify genes that may be related to changes in gene function in PTC pathogenesis. As a result, the enriched MF, BP and CC terms were determined to be involved in ‘Olfactory receptor activity’, ‘G-protein coupled receptor signaling pathway’ and ‘Signal receptor activity in the plasma membrane’ (Fig. 5A). In KEGG analysis, the significantly enriched pathways were involved in the ‘cytosolic DNA-sensing pathway’, the ‘RIG-I-like receptor signaling pathway’ and the ‘Olfactory transduction’ pathway (Fig. 5B). These results indicate that these dysregulated circRNAs may cause PTC cells to undergo a series of molecular and signal changes.

Figure 5. GO and KEGG analysis of the target genes of seven validated dysregulated circRNAs. (A) GO includes three domains: Biological Process, Cellular Component and Molecular Function. The P-value produced by top GO indicates the significance of GO terms abundance in seven dysregulated circRNAs targeting genes. (B) The seven pathways of seven validated dysregulated circRNAs targeting genes were identified using KEGG analysis according to the P-value. GO, Gene Oncology; KEGG, Kyoto Encyclopedia of Genes and Genomes; circRNA, circular RNA.

ceRNA analyses of dysregulated circRNAs

Among these potential circular RNAs, hsa_circ_IPCEF1 was related to 207 ceRNA genes out of the 256 ceRNA genes, exhibiting the strongest role in gene regulation (Table SII; Fig. 6A). Therefore, 207 ceRNA genes related to hsa_circ_IPCEF1 were further annotated by KOBAS 3.0, which were enriched in a variety of pathways and diseases. The top 10 items of KEGG pathway and diseases are presented in Fig. 6B, including ‘MAPK signaling pathway’, ‘PI3K-AKT signaling pathway’, ‘TNF signaling pathway’, ‘Olfactory transduction’, ‘Salmonella infection’, ‘Relaxin signaling parhway’, ‘Acute myeloid leukemia’, ‘Fcepsilon RI signaling pathway’, ‘Shigellosis’, ‘Bacterial invasion of epithelial cells’ and ‘cancers of endocrine organs’, ‘cancer’, ‘endocrine, metabolic diseases’, respectively, indicating the role of hsa_circ_IPCEF1 in PTC regulation (Fig. 6B).

Figure 6. ceRNA analysis of hsa_circ_IPCEF1. (A) All 207 predicted ceRNA genes of hsa_circ_IPCEF1. (B) The top 10 enriched KEGG pathways and diseases of hsa_circ_IPCEF1 related ceRNAs. (C) The ceRNA network of hsa_circ_IPCEF1, which contains hsa_circ_IPCEF1 (represented by yellow node), 21 miRNAs (represented by blue nodes) and four target gene mRNAs related to thyroid cancer (represented by green nodes). ceRNA, competing endogenous RNA; KEGG, Kyoto Encyclopedia of Genes and Genomes; circ, circular RNA.

Based on a comprehensive search of other studies in human malignancies, it was found that four ceRNAs (CASR, CDC25B, NFκB1 and SHOC2) were cancer-related target genes, which were reported to be regulated by 21 miRNAs (hsa-miR-16-5p, hsa-miR-195-5p, hsa-miR-15a-5p, hsa-miR-6838-5p, hsa-miR-15b-5p, hsa-miR-497-5p, hsa-miR-424-5p, hsa-miR-4530, hsa-miR-1323, hsa-miR-3691-5p, hsa-miR-545-3p, hsa-miR-128-3p, hsa-miR-922, hsa-miR-3918, hsa-miR-3619-5p, hsa-miR-608, hsa-miR-1262, hsa-miR-761, hsa-miR-214-3p, hsa-miR-4651 and hsa-miR-370-3p, data not shown). Therefore, the circRNA-miRNA-mRNA module containing four mRNAs and 21 miRNAs was constructed in Cytoscape software to concentrate the ceRNA targets (Fig. 6C).

Among these miRNAs, hsa-miR-3619-5p, reported as a tumor inhibitor, was one of the top MREs of hsa_circ_IPCEF1, with two seed sequence matching regions (44). Bioinformatics screening suggested that CDC25B (NM_021873) and CASR (NM_000388) were potential direct target genes of hsa-miR-3619-5p, with complementary binding sites located in the 3′-UTRs of CDC25B and CASR mRNAs, respectively (Fig. 7). These bioinformatics analyses indicated that hsa_circ_IPCEF1/hsa-miR-3619-5p/target genes may be involved in PTC pathogenesis.

Figure 7. hsa_circ_IPCEF1, miRNAs and target genes interactions. (A) Predicted top five miRNAs for hsa_circ_IPCEF1 by ceRNA analysis. (B) Seed match regions of hsa_circ_IPCEF1, hsa_miR_3619-5p and target genes. miRNA, microRNA; ceRNA, competing endogenous RNA; circ, circular RNA.

Prediction value of dysregulated circRNAs

To assess the performance of these seven circRNAs in predicting their potential association with PTC, their ROC curves were calculated and the accuracy evaluated with the areas under the ROC curves (AUCs). The preferable AUC value for circRNAs were hsa_circ_0021549 (AUC: 0.8194; 95% CI: 0.7117 to 0.9270; P<0.0001) and hsa_circ_IPCEF1 (AUC: 0.801095% CI: 0.7108 to 0.8912; P<0.0001; Fig. 8A-G). Instead of other circRNAs, hsa_circ_IPCEF1 exhibited the strongest functional role in ceRNA analysis, therefore, the correlation between hsa_circ_IPCEF1 expression and clinicopathological characteristics was examined in 57 patients with PTC. Statistical analysis showed that the downregulation of hsa_circ_IPCEF1 expression was associated with tumor lymph node metastasis in patients with PTC, suggesting its clinical value for PTC diagnosis and prognosis (Table I). To further confirm the formation of hsa_circ_IPCEF1 in patients with PTC, the splicing junction of hsa_circ_IPCEF1 was validated using RT-PCR and sequencing methods. As shown in Fig. 8H, hsa_circ_IPCEF1 is generated from the IPCEF1 gene located on human chromosome 6 (154520801–154544377). The 245bp products of hsa_circ_IPCEF1 in whole blood from subjects were amplified using divergent primers (Fig. 8I). Sanger sequencing further confirmed that the head-to-tail junction sites (AG/GC) were consistent with the hsa_circ_IPCEF1 annotation (Fig. 8J).

Figure 8. Evaluation candidate circRNAs. (A-G) Receiver-operating characteristic curve analysis of seven circRNAs in 57 paired whole blood of PTC patients. n=57 pairs of subjects. Comparisons of data were acquired by one-way ANOVA followed by Bonferroni post hoc. Identification of hsa_circ_IPCEF1 circularization. (H) Annotation of hsa_circ_IPCEF1. (I) Reverse transcription PCR and (J) sequencing analysis of head-to-tail splicing junction of hsa_circ_IPCEF1 in whole blood from PTC and healthy subjects. circRNA/circ, circular RNA; PTC, papillary thyroid cancer.

Discussion

circRNAs, a special novel class of endogenous non-coding RNAs, have been documented to serve crucial roles in various diseases by regulating numerous cellular events (45). Although numerous studies have focused on circRNAs as cancer markers, the diagnostic ability and understanding of circRNAs in PTC are still unclear. Therefore, the present study aimed to profile the expression levels of all circRNAs in patients with PTC, determine their potential roles in PTC and explore new insights into the pathogenesis of TC.

From microarray analysis, an aberrant expression pattern of circRNAs in PTC tissues compared to benign thyroid tissues was found. In liquid biopsy, the usage of whole blood is the main source of body fluid and circRNAs in blood cells or whole blood have been proposed as biomarkers for human diseases due to their high stability, abundance and spatiotemporal specific expression (46). Therefore, three upregulated and four downregulated circRNAs were verified in whole blood from 57 pairs of patients with PTC compared to healthy subjects and validated by RT-qPCR.

PTC originates from thyroid follicular epithelial cells, accounting for 89.8% of differentiated thyroid carcinomas, while only 38% of thyroid cancers show clinical symptoms (47). However, the extensive implementation of cancer imaging studies in preoperative diagnostics has resulted in a growing incidence of low-risk PTCs, which has caused a continuing debate about the adequacy of aggressive surgical approaches and radioiodine therapy (48). The molecular predictors of preoperative diagnostics and accurate treatment for PTC are poorly defined. Although studies have focused on mRNA transcripts and non-coding RNAs as potential PTC markers, the expression profiles and roles of circRNAs in PTC remain to be elucidated (49–52). By contrast with recent studies in circRNA screening in PTC tissues, the present study validated seven dysregulated circRNAs (three upregulated and four downregulated) in whole blood from PTC subjects, which were found to be more suitable for biomarker detection in the future (53,54). To explore the mechanisms of these dysregulated circRNAs on the pathological processes of PTC, GO and KEGG pathway analyses were performed. KEGG pathway analysis demonstrated that the olfactory transduction pathway was the most significant pathway for enrichment, followed by the cytosolic DNA-sensing pathway and the RIG-I-like receptor signaling pathway. Moreover, according to the BP, CC and MF terms with substantial enrichment, the genes were mainly associated with ‘G-protein-coupled receptor activity’, ‘Signal receptor activity in the plasma membrane’ and ‘Olfactory receptor activity’. These results suggested that the deregulated circRNAs may participate in the pathological process of PTC. Recently, circRNAs have been well described as ceRNAs in human diseases, as has the possibility of using circRNAs as molecular markers for related diseases (55). Compared with miRNAs, circRNAs act on miRNA sponges to attenuate the effect of miRNAs on mRNA expression (42). For example, circCRIM1 promotes the expression of BTG2 by inhibiting the expression of miR-125b-5p in LUAC cells, thus affecting the growth of LUAC cells (37). The loss of circ-znf609 inhibits the proliferation and cell cycle transition and induces the apoptosis of NPC cells via modulation of the miR-188/elf2 axis (35). In bladder cancer, circRNA-cTFRC is upregulated and correlated with tumor grade and survival rate, acting as a miR-107 sponge to regulate TFRC expression, cancer cell invasion, proliferation, epithelial to mesenchymal transition and tumor growth (56). These findings suggest the key roles of abnormal circRNAs and target genes in the pathogenesis of disease, indicating a new perspective for the early diagnosis, treatment and prognosis of PTC.

Using a set of bioinformatics tools, the present study predicted MREs for the seven validated circRNAs and the potential targets of miRNAs. The present study found one promising downregulated circRNA (hsa_circ_IPCEF1) in PTC subjects, which was not indicated in any previously published literature but showed interactions with four cancer-related ceRNAs (CASR, CDC25B, NF-κB1 and SHOC2). Notably, NF-κB1 is a transcription regulator that serves important roles in numerous biological processes, such as survival and proliferation, inflammation and adaptive immune responses and has been widely implicated in the development and progression of cancer pathogenesis (57–60). In cancer cells, chronic activation and nuclear localization result in NF-κB1 accumulation and NF-κB activation, which stimulate cell survival by inducing anti-apoptotic gene expression (61). Recently, NF-κB1 has been proved to be upregulated in a variety of cancers, including subungual keratoacanthoma, urothelial carcinoma and breast cancer, suggesting that it is a biomarker for diagnosis and treatment (62–66). In the results of the present study, NF-κB1 was predicted to be a ceRNA of hsa_circ_IPCEF1 regulation and the downregulation of hsa_circ_IPCEF1 in PTC was accompanied by the upregulation of NF-κB1, suggesting that a combination of hsa_circ_IPCEF1 and NF-κB1 could be used as biomarkers for PTC diagnosis/prognosis. The ROC analyses of RT-qPCR results, as well as correlation analysis between hsa_circ_IPCEF1 and clinical parameters, further confirmed hsa_circ_IPCEF1 as one preferable candidate biomarker for PTC. Moreover, the circRNA-miRNA-mRNA module suggested that hsa-miR-3619-5p was a direct target for hsa_circ_IPCEF1 sponging. As a tumor suppressor, hsa-miR-3619-5p inhibits prostate cancer cell growth and hampers the proliferation of cutaneous squamous cell carcinoma and liver cancer cells, exerting tumor inhibitory effects on human malignancies (67–69). Bioinformatics analysis further confirmed the interactions between hsa-miR-3619-5p and CDC25B (an initiator of mitosis) and CASR (calcium-sensing receptor) mRNAs, which may be direct targets of hsa_circ_IPCEF1/hsa-miR-3619-5p axis regulation. By controlling the G2/M cell cycle phase transition and calcium metabolism, CDC25B and CASR serve important roles in the occurrence and metastasis of cancer, respectively. These results suggested that hsa_circ_IPCEF1 may participate in the pathological process of PTC through sponging hsa-miR-3619-5p. Clearly, the functions and mechanisms of the hsa_circ_IPCEF1/miR-3619-5p axis need to be further investigated in PTC cells to conclusively determine which target genes interact which cellular functions are affected during tumorigenesis. Taken together, the present study identified distinctive circRNA expression patterns in PTC. hsa_circ_IPCEF1 may bind to hsa-miR-3619-5p through sponge action, thus regulating the target genes of cell proliferation, migration and invasion and triggering the initiation and progression of PTC. These results provide insight into the pathogenesis of PTC and suggest the potential roles of hsa_circ_IPCEF1 in PTC diagnosis and treatment.

There are still some potential limitations and shortcomings in the present study. First, in PTC diagnosis, diagnostic performance was primarily based on neck ultrasonic diagnosis rather than expensive fine-needle aspiration cytology. PTC tumors were obtained from aggressive surgical approaches on patients with high-risk PTC, who were diagnosed and evaluated by neck ultrasonic imaging and immunohistochemistry. Limited by the small size of PTC tumors, the present study did not perform diagnostic assays on these samples to clarify the PTC subtype and genotype. However, the PTC tumors obtained for microarray were conventional PTC, exhibiting fibrovascular cores and enlarged nuclei with nuclear clearing, which were sufficiently accurate for the present study. Second, the correlation between the expression level of hsa_circ_IPCEF1 and the clinical parameters in patients with PTC was analyzed only in 57 pairs of subjects. A new set of patients with PTC should be recruited to test the predictive value and the utility of hsa_circ_IPCEF1 as a prognostic biomarker in the near future. Nevertheless, the present study provided new insight into the pathogenesis of PTC at the molecular level and proposes one potential candidate PTC biomarker. Third, the analysis of the present study was limited to the tissue and whole blood of patients with PTC and the biological function of the maladjusted circRNAs was predicted only through bioinformatics. The effects of the dysregulated circRNAs on the proliferation and migration of PTC tumor cells still need to be studied in vivo and in vitro. Also, the construction of the ceRNA network is based on the sponge function of circRNAs; other regulatory functions still need to be further studied. The mechanisms and functions of the hsa_circ_IPCEF1/hsa-miR-3619-5p axis in PTC will be studied in the near future.

In conclusion, the present study reported seven significantly differentiated circRNAs in patients with PTC compared to healthy subjects, suggesting a key role of dysregulated circRNAs in the pathogenesis of PTC. As a downregulated circRNA, hsa_circ_IPCEF1 may show promise as a potential biomarker for PTC. Moreover, the bioinformatics analysis suggested that the hsa_circ_IPCEF1/hsa-miR-3619-5p axis may be involved in the pathogenesis of PTC, providing a new theoretical basis for the future diagnosis and treatment of PTC (Fig. 9).

Figure 9. Schematic diagram of the proposed role of hsa_circ_IPCEF1 in PTC. PTC, papillary thyroid cancer; circRNA, circular RNA.

Supplementary Material

Supporting Data

Supporting Data

Acknowledgements

Not applicable.

Funding

This study was supported by the National Natural Science Foundation of China (Contract grant no. 81770054); National Natural Science Foundation of China (Contract grant no. 8197052197); Chinese Natural science foundation of Heilongjiang province (grant no. LH2019H017); Tianjin Municipal Science and Technology Project (grant no. 18JCZDJC44900); The Fundamental Research Funds for the Central Universities-Nankai University (grant no. 63191150); Spark Research Fund from The Fourth Affiliated Hospital of Harbin Medical University (grant no. HYDSYXH201907).

Availability of data and materials

The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.

Authors' contributions

SSL and XJ designed the study and interpreted the data. MG analyzed and interpreted the patient data regarding papillary thyroid carcinoma. SSL and YS performed the histological examination of PTC, molecular studies and were major contributors in writing the manuscript. JD, YL and SY performed RT-qPCR experiments. YZ and XJ performed bioinformatics analysis. SSL and XJ confirm the authenticity of all the raw data. All authors read and approved the final manuscript.

Ethics approval and consent to participate

The study procedure was approved by Medical Ethics Committee of Harbin Medical University (approval no. IRB3011619).

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

References