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

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%.¹ 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.² 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.³ However, the tumorigenesis of PTC has not yet been fully elucidated, alongside other microcarcinomas which also lead to severe symptoms.⁴ 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.⁷ 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.⁸ The special circular structure of circRNAs, closed circular sequences, is more stable than linear RNAs, which have 5ʹ and 3ʹ ends.⁹ Moreover, longer circRNAs may also locally appear as double stranded RNAs in order to increase their stability.¹⁰ 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.¹³ 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.¹⁴ Accumulating research evidence has opined that circRNAs are usually aberrantly expressed in human cancers, and participate in oncogenesis via multiple approaches.¹⁵ 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

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

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.¹⁶ 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

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.2¹⁷ 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.¹⁸ We then collected the unmapped reads for circRNA identification using Find_circ v1,¹⁹ followed by blasting against circBase²⁰ 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.0²¹ 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.²²

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

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. Abbreviations: NC, normal cancer tissues; T, tumor 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).

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

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.

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.

circRNA-miRNA-mRNA Network Construction

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,²³ 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.²⁴ 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.²⁵ 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.³⁰ 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.³¹ One study reported that circRNAs are involved in the pathogenesis of PTC through modulating miRNAs and their target mRNA.³² 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,³³ prostate cancer,³⁴ hepatocellular carcinoma,³⁵ and cholangiocarcinoma.³⁶ In addition, Hsa-mir-22-3p regulates lung adenocarcinoma progression,³⁷ while Wnt9A activation is related to repression of human colorectal cancer cell proliferation.³⁸ Furthermore, FOXA2 suppresses the metastasis of hepatocellular carcinoma,³⁹ and ZEB1 activation induces epithelial-mesenchymal transition.⁴⁰ 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.