A total of 50 patients who underwent surgical treatment with a final diagnosis of PTC at Peking Union Medical College Hospital (Beijing, China) between January 2017 and December 2018 were enrolled in this study, and the clinicopathological characteristics of the patients are listed in Tables I and III. The recruited patients did not receive any adjunctive treatment prior to surgery, such as radiotherapy, chemotherapy or targeted therapy. Clinical information, including patient age, sex, multifocality, subtype (classical/follicular variant PTC), Ki67, vascular endothelial growth factor (VEGF), pT, pN, metastasis, BRAF mutation, capsular invasion and thyroiditis, was collected from the clinical records and stratified as characteristics according to the 8th edition of the American Joint Committee on Cancer (AJCC) Cancer Staging Manual (27).

A total of 50 paired specimens (tumour and normal tissues) were obtained from the same patients during surgery. Following gross examination, the tissues were resected and washed with PBS (pH 7.2). One part of each specimen was immediately snap-frozen using liquid nitrogen within 30 min and cryopreserved in liquid nitrogen (−196°C) until RNA sequencing and gene expression analysis. The remaining part of the specimen was fixed in formalin at room temperature for 48 h and embedded in paraffin, and then 4-µm sections were cut from each block and stained with haematoxylin and eosin (H&E). The morphology of H&E slides were reviewed by 2 experienced pathologists to confirm the final diagnosis (28).

The present study was approved by the Institutional Review Board of Peking Union Medical College Hospital, and written informed consent was obtained from all patients regarding the use of their tissues for research purposes.

A total amount of 5 µg RNA per sample was used as input material for RNA sample preparation. First, ribosomal RNAs were depleted using a Ribo-Zero™ rRNA Removal Kit (EPIcentre) to obtain rRNA-depleted RNAs. The rRNA-depleted RNAs were further treated with RNase R (EPIcentre) and were then subjected to TRIzol extraction (Thermo Fisher Scientific, Inc.). Subsequently, sequencing libraries were generated using the rRNA-depleted and RNaseR-digested RNAs by NEBNext Ultra™ Directional RNA Library Prep Kit for Illumina^® (New England Biolabs) following the manufacturer's recommendations. Then, clustering of the index-coded samples was performed on a cBot Cluster Generation System using HiSeq PE Cluster Kit v4 cBot (Illumina, Inc.) according to the manufacturer's instructions. After cluster generation, the library preparations were sequenced on an Illumina Hiseq 2500 platform and 125-bp paired-end reads were generated.

Subsequently, the circRNAs were identified according to the manual (10). Differential expression in normal tissues and PTC was evaluated by using DESeq2 from Dr. Michael Love (version 1.6.3) (29). P-values were adjusted as q (Padj) values by the Benjamini & Hochberg method (30); both q (Padj) <0.05 and fold-change in expression of >2 was set as the threshold for differential expression by default.

Thyroid tissues were processed using the IKA ULTRA-TURRAX T10 basic Disperser/Homogenizer (IKA Works, Inc.), and then total RNA was extracted using TRIzol reagent (Thermo Fisher Scientific, Inc.) according to the manufacturer's instructions. The total RNA concentrations and quality were determined using Nano Drop One (Thermo Fisher Scientific, Inc.). The RNA integrity was assessed through denaturing agarose gel electrophoresis.

The RT-qPCR experiments were conducted as previously described (31,32). Total RNA was digested with RNaseR at 37°C for 30 min and reverse-transcribed to synthesise cDNA using M-MLV reverse transcriptase (Promega Corporation; cat. no. M1701) and random primers at 70°C for 10 min followed by a 60-min incubation at 37°C. Then the cDNA was analysed by qPCR with Fast SYBR^® Green Master Mix (Applied Biosystems; Thermo Fisher Scientific, Inc.; cat. no. 4385612) and a 7500 Fast Real-Time PCR System (Applied Biosystems; Thermo Fisher Scientific, Inc.). Amplification conditions were 95°C (20 sec), 40 cycles of 95°C (3 sec) and 60°C (30 sec) according to the manufacturer's instructions. GAPDH was used as the internal control. The divergent primers used for circRNA validation are listed in Table SI, the specificity of primers was assessed in NCBI BLAST (https://blast.ncbi.nlm.nih.gov/Blast.cgi), and the backsplice junction of circRNA in the PCR product was verified with Sanger sequencing. The relative expression level of each circRNA with GAPDH was calculated using the 2^−ΔΔCq method (24,33).

IHC staining was performed on 4-µm sections using the DAKO Autostainer Link 48 (Dako; Agilent Technologies, Inc.). The tissue epitopes were repaired using the automated water bath heating process in Dako PT Link (Dako; Agilent Technologies, Inc.); the sections were incubated in Tris-EDTA retrieval solution (10 mM Tris, 1 mM EDTA; pH 9.0) at 98°C for 20 min. The sections were subsequently incubated for 50 min at room temperature with the primary anti-Ki67 (cat. no. MIB1; Dako; Agilent Technologies, Inc.) and anti-VEGF (cat. no. VG1; Dako; Agilent Technologies, Inc.) antibodies, diluted at 1:200 and 1:100, respectively, in Dako Envision™ Flex antibody diluents, followed by anti-rabbit immuno-peroxidase polymer (Envision FLEX/HRP) for 20 min at room temperature according to the manufacturer's instructions, and developed with freshly prepared 0.05% 3,3′-diaminobenzidine tetrahydrochloride. Finally, the slides were counterstained with hematoxylin, dehydrated, and mounted.

Genomic DNA was extracted from fresh-frozen tissues using the TIAamp Genomic DNA Kit (Tiangen Biotech Co., Ltd.; cat. no. DP304). Subsequently, the genomic DNA was used to detect the V600E mutation in the BRAF oncogene through the AmoyDx BRAF^V600E Mutations Detection Kit (AmoyDx; cat. no. ADx-BR04).

Statistical analyses were performed using SPSS 19.0 (IBM Corp.). The significance of the RT-qPCR validation between the PTC and the normal tissue groups was assessed by the paired Student's t-test (18,34,35). The Mann-Whitney U test was used to analyse the associations between the mean ± standard deviation of circRNA expression and the clinicopathological characteristics of PTC, including age, sex, multifocality, Ki67, pN stage, BRAF mutation status, capsular invasion and thyroiditis, and the Kruskal-Wallis test followed with post-hoc Dunn's multiple comparison test was used to analyse the association between circRNA expression and pT-stage (35,36,37). P<0.05 was considered to indicate statistically significant differences (38). Receiver operating characteristics (ROC) curves were established to evaluate the diagnostic value of circRNAs for PTC. An area under the ROC curve (AUC) of 0.5 was considered to represent a test with no discriminating ability (i.e., no better than chance), while an AUC of 1.0 was considered to represent a test with perfect discrimination (39).

The clinicopathological characteristics of the 5 patients are listed in Table I. The majority were male (4 cases, 75%), 3 cases were classical PTCs, 1 case was follicular PTC and 1 case was mixed classical and follicular PTC, with different tumour-node-metastasis (TNM) stages. Two cases were multifocal and all displayed capsular invasion. The histomorphology of tissues are shown in Fig. 1. Among them, Fig. 1A and C shows the histomorphology of normal tissues and with lymphocyte infiltration. Fig. 1B and represent the classical PTC and follicular variant of PTC, respectively.

Subsequently, RNA-seq demonstrated that the circRNA expression profile of the 5 paired normal and PTC tissue samples exhibited significant differences. The hierarchical clustering (Fig. 2A) and volcano plots (Fig. 2B) revealed that there were 53 circRNAs differentially expressed in carcinoma tissues, with both q (Padj) <0.05 and fold-change in expression of >2 set as the threshold. Among those, 45 circRNAs were upregulated and 8 were downregulated in tumour tissues.

In the dysregulated circRNAs (Table SII), the Log₂ (FC) of abnormally highly expressed circRNAs ranged from 1.48 to 7.34, and the Log₂ (FC) of abnormally lowly expressed circRNAs ranged from-1.31 to −5.593. In addition, the magnitude of fold change was highest for chr5:161330882-161336769- in the upregulated circRNAs (fold change, 162; P<0.01), and chr20:20425608-20472956- for downregulated circRNAs (fold changes, 48; P=0.01). As presented in Fig. 2C, the transcription of the differentially regulated circRNAs was found to be widely distributed among all chromosomes, except chr10, chr13, chr14, chr17, chr18, chr19, chr22 and chrM. Moreover, 30 were identified as new circRNAs that have not been annotated in the circBase or circ2 Traits database. Furthermore, multiple circRNAs can be produced from a single gene and alternative splicing events, so circRNAs can arise from exons, introns, intergenic regions, UTRs, ncRNA loci and from locations antisense (40,41). In our research, of the 53 dysregulated circRNAs (Table SII), 29 (54.7%) were from the antisense locus, 10 (18.9%) were formed from exon-intron junction, 8 (15.1%) were from intron, 4 (7.5%) were from exon, and the other 2 were from intergenic regions, respectively.

A total of 8 differentially expressed circRNAs were selected according to transcripts per million reads (TPM), fold change, q (Padj) value, junction read, features (arise from exonic/intronic), distribution and length, and then were verified using 45 pairs of tissues by RT-qPCR. The circRNA relative expression level of GAPDH was calculated using the 2^−ΔΔCq method. The results confirmed that the upregulation of chr5: 161330882-161336769- (P<0.001), chr12:129699809-129700698- (P<0.05), chr9:22046750-22097364+ (P<0.001), chr20:17456347-17465553+ (P<0.001), chr7:116699070-116700284+ (P<0.01), chr8:18765448-18804898- (P<0.001) and chr7:22308338-22318037- (P<0.001) and the downregulation of chr5:38523418-38530666- (P<0.001) in carcinoma (Fig. 3) were in accordance with the RNA-seq data (Fig. 2D).

Furthermore, the mRNA expression of these 8 circRNAs was accessed with ROC curve analysis (Fig. 4A), which revealed that all the AUCs of circRNAs were >0.7 (P<0.001), and chr5:161330882-161336769- with the highest AUC (0.878) could be considered as a powerful marker for clinical PTC diagnosis (Table II).

The clinical characteristics of 45 PTC patients are summarized in Table III. The median age of the patients at diagnosis was 45 years (range, 17–64 years), while the majority of the patients were female (n=33, 73.3%); 32 cases (71.1%) were classical PTC, 12 cases (26.7%) were the follicular variant of PTC and only 1 case (2%) was a mixture of the classical and follicular variant. Moreover, 31 cases (68.9%) were unifocal, and 14 cases (31.1%) were multifocal. In regards to pT stage, 29 cases (64.4%) had a diameter <2 cm or localized, 4 (8.9%) ranged between 2 and 4 cm, and 11 (24.4%) were >4 cm or invading the soft tissues beyond the thyroid. Only one-quarter of the patients had regional lymph node metastasis. None of the patients had metastasis. The Ki67 index of 33.3% patients was very high (≥3), and VEGF was positive for 91.1% patients. BRAF^V600E mutation, capsular invasion and thyroiditis with diffuse lymphocytic infiltration were observed in 80, 75.6 and 20% of the cases, respectively.

Subsequently, we analysed the association between the mean ± standard deviation of circRNA expression in the tumour and the clinical characteristics in 45 PTC patients, including age, sex, multifocality, Ki67, capsular invasion, VEGF, BRAF^V600E mutation, TNM stage and thyroiditis (Table III). Among the 8 dysregulated circRNAs, chr5:161330882-161336769- (P=0.015), chr9:22046750-22097364+ (hsa_circ_0008796, P=0.041) and chr8:18765448-18804898- (hsa_circ_0002111, P=0.036) were associated with the BRAF^V600E mutation, and chr12:129699809-129700698- was associated with capsular invasion (P=0.025). Additionally, chr5:38523418-38530666- (hsa_circ_0072309) was associated with pT (P=0.037) and pN (P=0.002) stage. This association between circRNA expression and clinicopathological characteristics indicated that those circRNAs may participate in tumourigenesis and progression of PTC. The other three dysregulated circRNAs exhibited no statistically significant association with clinicopathological characteristics.

Furthermore, the association between the mRNA expression of these 5 circRNAs and the clinical characteristics was assessed with ROC curve analysis (Fig. 4B), revealing that the AUCs of all circRNAs were >0.7 (P<0.05), with chr5:38523418-38530666-(associated with pN stage) having the highest AUC (0.8209, P<0.01; Table IV). This result was consistent with the data mentioned above, and indicated that those 5 circRNAs associated with BRAF^V600E mutation, capsular invasion, advanced pT stage and lymph node metastasis may be putative biomarkers for the diagnosis and evaluation of the progression of PTC.

GO enrichment analysis for the host genes of the identified differentially expressed circRNAs was performed, and the significantly enriched GO terms in the biological process, cellular component and molecular function categories, e.g., carbohydrate derivative catabolic process (GO: 1901136, P=0.00021557), membrane (GO: 0016020, P=0.0012195) and GTPase regulator activity (GO: 0030695, P=0.00021632), are presented in Fig. 5A. A total of 28 host genes of the differentially expressed circRNAs were associated with the membrane GO term. The results of the Kyoto Encyclopaedia of Genes and Genomes pathway analysis of the 20 highly related circRNA genes are presented in Fig. 5B. Those circRNA genes were mainly associated with important pathways (Rap1, VEGF and Ras signalling pathways), tumourigenesis (transcriptional misregulation and proteoglycans in cancer), cell signal transduction (axon guidance, retrograde endocannabinoid signalling, cytokine-cytokine receptor interaction), infection and immune defence (bacterial invasion of epithelial cells, chemokine signalling pathway). Furthermore, the Rap1 and Ras signalling pathway contained more related circRNA genes, which suggested that circRNAs may play a key role in PTC tumourigenesis and progression.

circRNAs have been shown to contain miRNA-binding sites and negatively regulate the inhibitory effects of miRNAs on their target mRNAs, thereby essentially acting as miRNA sponges (14). To evaluate the potential functions of the identified differentially expressed circRNAs, miRNAs that were potentially associated with those circRNAs were investigated using the RegRNA2 platform and Circular RNA Interactome (42,43). The first 5 predicted miRNAs for the dysregulated circRNAs are listed in Table V. For example, circRNA chr5:161330882-161336769- was predicted to harbour hsa-miR-1273g-3p, hsa-miR-498, hsa-miR-1268a, has-miR-363-5p and hsa-miR-566.

To further explore which miRNAs may be cancer-related, the starBase v2.0 (44) and miRPathDB v1.1 (45) platforms were used, and a network map of circRNA-miRNA interactions was constructed to identify the dysregulated circRNAs in PTC that potentially interacted with cancer-related miRNAs. From the initial list of 40 potential miRNAs identified by the RegRNA2 and circRNA Interactome analysis, several cancer-related miRNAs were identified following a literature research, including hsa-miR-498 (46), hsa-miR-766-3p (47,48) and hsa-miR-661 (49,50), which are associated with cell proliferation in several cancers. The results indicated that those dysregulated circRNAs may act as miRNA sponges, leading to the progression of PTC, but the specific underlying mechanism requires further investigation.