Identification of Altered Circular RNA Expression in Serum Exosomes from Patients with Papillary Thyroid Carcinoma by High-Throughput Sequencing.

BACKGROUND This study aimed to identify altered exosome circular RNA (circRNA) in the serum of patients with papillary thyroid carcinoma using high-throughput sequencing. MATERIAL AND METHODS Serum was collected from three patients with papillary thyroid carcinoma and three patients with a benign thyroid goiter. Exosomes were isolated using an exosome isolation kit and confirmed by transmission electron microscopy. Exosome circRNAs were analyzed by high-throughput sequencing using the HiSeq 4000 sequencer. The differentially expressed circRNAs were confirmed by fluorescence quantitative real-time polymerase chain reaction (qRT-PCR). RESULTS Twenty-two differentially expressed circRNAs were screened, which included three that were upregulated and 19 that were down-regulated in serum from patients with papillary thyroid carcinoma compared with controls. Gene Ontology (GO) enrichment analysis showed that these differentially expressed circRNAs were associated with 16 signaling pathways, including the thyroid hormone signaling pathway, the PI3K-Akt signaling pathway, and the AMPK signaling pathway. Three differentially regulated circRNAs included hsacirc_007293, hsacirc_031752, and hsacirc_020135 were confirmed by qRT- PCR. The expression trends were consistent between the high-throughput sequencing technique and qRT-PCR. CONCLUSIONS The findings from this study have shown that gene regulation can be studied from exosomes obtained from serum samples in patients with papillary thyroid carcinoma, and supports the need for further studies on the role of exosome circRNAs in thyroid cancer.

Background

Thyroid cancer is the most common type of malignant endocrine tumor of the head and neck, with an incidence of 1–2% of all types of cancer [1,2]. Clinically, thyroid cancer is divided into differentiated thyroid carcinoma, medullary carcinoma, and undifferentiated carcinoma [3,4]. The most common thyroid cancer, differentiated thyroid carcinoma, includes papillary thyroid carcinoma and follicular thyroid carcinoma. Papillary thyroid carcinoma accounts for 75% of all types of thyroid cancer [5]. With the improvement in living standards and the increase in population, the global incidence of papillary thyroid carcinoma is increasing annually, and the number of patients who present with lymph node metastasis is also increasing in frequency [6,7]. Also, the incidence of papillary thyroid carcinoma is increasing, particularly in women, which highlights the importance of studies on the molecular pathogenesis of this malignancy and the identification of diagnostic biomarkers that may improve early diagnosis and treatment [2,4,8].

Circular RNA (circRNA) is a non-coding RNA without a 5′-cap and 3′-polyadenylate tail. In 1976, the first circRNA was identified in a virus [9]. With the development of RNA sequencing and single-cell RNA sequencing technology, increasing numbers of circRNAs have been screened and identified in different species, including humans [10-15]. According to their origin, circRNAs can be divided into three categories, exon sequence, intron sequence, and combined exon and intron sequence [16,17]. Also, circRNA has been reported to regulate the expression of microRNA (miRNA) and to function as microRNA sponges [18,19]. Further research has shown that circRNAs have important roles in the pathogenesis of human diseases, including cancer [20,21].

Exosomes are small endocrine vesicles that exist in most cells and contain transporters, mRNAs, and miRNAs [22], and in exosomes, these molecules have a potential role as diagnostic biomarkers of human disease [23]. Recent studies have also shown that a large number of circRNAs exist in exosomes [24,25]. RNA sequencing data showed that the concentration of circRNA in exosomes was greater than in normal cells [24,25], and a further study has indicated that exosome circRNAs could reflect circRNA levels in cells and tissues [26]. Therefore, exosome circRNAs might be potential biomarkers for thyroid carcinoma.

This study aimed to identify altered exosome circRNA in the serum of patients with papillary thyroid carcinoma using high-throughput sequencing and to compare the findings with control patients diagnosed with benign nodular goiter. To avoid the false-positive results, fluorescence quantitative real-time polymerase chain reaction (qRT-PCR) was used to confirm the data obtained by high-throughput sequencing technique.

Material and Methods

Patients

Three patients with papillary thyroid carcinoma were enrolled in the study, who were diagnosed and treated at the Second Affiliated Hospital of Nanchang University and three patients with benign thyroid goiter were included as the study controls. Written informed consent was obtained from each patient. This study was approved by Ethics Committee of the Second Affiliated Hospital of Nanchang University. The clinical and demographic characteristics of the study participant are shown in Table 1.

Table 1

The clinical information of the patients included in the study.

Sample	Gender	Age (years)	Clinical diagnosis
Control 1	Female	55	Nodular goiter
Control 2	Female	62	Nodular goiter
Control 3	Female	46	Nodular goiter
Cancer 1	Female	44	Papillary thyroid carcinoma
Cancer 2	Female	48	Papillary thyroid carcinoma with lymphocytic thyroiditis
Cancer 3	Female	30	Papillary thyroid carcinoma

Isolation of serum exosomes

Serum samples were centrifuged at room temperature for 10 min (2000×g) to remove residual cells and debris, and the supernatant was collected from each and transferred to a new centrifuge tube, and one-third of the volume consisted of GSTM Exosome Isolation Reagent A (Thermofisher Scientific, Waltham, MA, USA). The samples were completely mixed and centrifuged at room temperature for 10 min (13000×g). The supernatant was discarded, and GSTM Exosome Isolation Reagent B with the same volume as Reagent A was added, mixed, and centrifuged for 10 min at 13 000×g, and the supernatant was discarded. The precipitates that included the isolated exosomes were stored at −80°C. The exosomes were identified by transmission electron microscopy. Briefly, the precipitates were fixed in 2.5% glutaraldehyde, dehydrated, embedded, and stained with 3% uranium acetate and lead citrate, and observed and photographed by transmission electron microscopy using a JEM-1230 electron microscope (JEOL Ltd., Akishima, Tokyo, Japan).

RNA extraction

Exosome RNA was extracted using TRIzol® reagent (Invitrogen, Carlsbad, CA, USA), as previously described [26]. The concentration and purity of the RNA were detected using the Nanodrop 2000 spectrophotometer (Thermofisher Scientific, Waltham, MA, USA). The integrity of the RNA was confirmed by agarose gel electrophoresis, and the RNA integrity number was measured using an automated Agilent 2100 Bioanalyzer system (Agilent, Santa Clara, CA, USA). The total RNA level was >5 μg, the concentration was ≥200 ng/μL, and the optical density (OD) at 260/280 nm was between 1.8 and 2.0. The ribosomal RNA (rRNA) and the linear RNA were removed using the Ribo-Zero Magnetic Kit and RNase R (Epicentre Technologies, Madison, WI, USA), respectively. The paired-end sequencing library was constructed according to TruSeqTM Stranded Total RNA Library Prep Kit assay (Illumina, San Diego, CA, USA). RNA sequences were identified using the HiSeq 4000 sequencer platform (Illumina, San Diego, CA, USA).

Data analysis

SepPrep and Sickkle software were used to measure the quality of the data, and the data were compared to reference genome data by Bowtie. KNIEF and CIRCE Explorer 2 were used to predict circular RNAs (circRNAs), and the expression of circRNAs in the samples was calculated and classified.

Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis

Differentially expressed gene enrichment from GO and the KEGG were analyzed by GOATOOLS and included molecular function, cell components, and biological processes. The KOBAS web server was used for gene and protein functional annotation and to enrich the KEGG pathway analysis. The target genes of circRNAs were predicted using the DESeq package.

Fluorescence quantitative real-time polymerase chain reaction (qRT-PCR)

Fluorescence quantitative real-time polymerase chain reaction (qRT-PCR) was used to verify the differential expression of circRNAs. According to the results of high-throughput sequencing, three differentially expressed circRNAs, hsacirc_020135, hsacirc_007293, and hsacirc_031752, were detected by qRT-PCR. U6 was used as an internal reference, as previously described [27]. The sequences of the primers used are listed in Table 2.

Table 2

Primer sequences of the circular RNAs (circRNAs).

	Sequence (5′-3′)	Primer length (bp)	Product length (bp)	Annealing temperature (°C)
Hsa-circ-020135 F	ATCCAGATAATGTATGGCTGCG	22	19215	58.6
Hsa-circ-020135 R	GGAGGAGGGCAAGAGCTGGA	20
Hsa-circ-007293 F	TGCCGGAAGAGGGGTTCCATG	21	299	60.5
Hsa-circ-007293 R	GCGAAGATTCTGCCCATCATGT	22
Hsa-circ-031752 F	CGTTTCTGCATTATTGCTCCC	21	233	56.2
Hsa-circ-031752 R	GGGACAGAAAGGATCCTCGG	20
U6 F	CTCGCTTCGGCAGCACA	17	94	60.0
U6 R	AACGCTTCACGAATTTGCGT	24

Statistical analysis

Data were analyzed using SPSS version 19.0. The differences between groups were compared using a t-test. A P-value <0.05 was considered to be statistically significant.

Results

Exosome confirmation and ultrastructure

Identification of the exosomes is shown in Figure 1. Transmission electron microscopy (TEM) showed that the exosomes were cystic with double membranes, and they were about 100 nm in diameter.

Figure 1

Structure of the exosomes shown by transmission electron microscopy (TEM). The exosome is cystic with double membranes, with a diameter of about 100 nm. Based on the findings from this study, circular RNAs (circRNAs) appear to be concentrated between the exosome membranes.

Quality analysis

The quality of the data was initially analyzed. As shown in Table 3, the sequence numbers in the control group were 89274794, 92765116, and 93254666, and the sequence numbers in the papillary thyroid cancer group were 868446436, 85036516 and 84788014. The bases with a Phred quality score >20 accounted for about 97% of the total bases, indicating that the original sequencing data was of good quality and could be used for subsequent analysis.

Table 3

Quality control of the original data.

Sample	No. of reads	Q30 (bp)	N (%)	Q20 (%)	Q30 (%)
Control 1	89274794	12735895669	0.006518	97.8	94.47
Control 2	92765116	13171999297	0.009318	97.59	94.03
Control 3	93254666	13307202465	0.009464	97.88	94.5
Cancer 1	86846436	12307924955	0.006775	97.51	93.85
Cancer 2	85036516	12132210290	0.006969	97.84	94.48
Cancer 3	84788014	11950214828	0.00877	97.02	93.33

As shown in Table 4, the total number of reference genome sequences was >90%, indicating that the quality control was good. Principal component analysis (PCA) showed that there was a significant difference between the sample groups, and no differences between sample repeats (Figure 2). These data were determined to be appropriate for use in subsequent analysis.

Table 4

Map of the RNA sequences identified by high-throughput sequencing.

Sample	Clean reads	Total mapped	Multiple mapped	Uniquely mapped
Control 1	88787808	81755230 (92.08%)	19625139 (24.00%)	62130091 (76.00%)
Control 2	92182006	83727339 (90.83%)	17432824 (20.82%)	66294515 (79.18%)
Control 3	93002378	85246689 (91.66%)	24293743 (28.50%)	60952946 (71.50%)
Cancer 1	85950894	79377677 (92.35%)	19018369 (23.96%)	60359308 (76.04%)
Cancer 2	84481058	79139236 (93.68%)	54246616 (68.55%)	24892620 (31.45%)
Cancer 3	83232062	75664229 (90.91%)	17226243 (22.77%)	58437986 (77.23%)

Figure 2

Results of principal component analysis (PCA). Data shows that there was a significant difference between the sample groups, and between the sample repeats.

Differentially expressed gene (DEG) analysis

Using DEGseq software, 22 differentially expressed circular RNAs (circRNA) were screened from the two study groups using the screening criteria of multiple expression of differentially expressed genes (DEGs) (change fold >1, and P<0.05). Three circRNAs were upregulated and 19 circRNAs were down-regulated in the serum from patients with papillary thyroid carcinoma compared with the control group (Figure 3).

Figure 3

Volcano map of the differentially expressed circular RNAs (circRNAs). The left indicates the downregulated circular RNAs (circRNAs) in the papillary thyroid carcinoma group, and the right indicates the upregulated circRNAs in the papillary thyroid carcinoma group.

Gene Ontology (GO) functional enrichment analysis

First, all genes were mapped to terms in the Gene Ontology (GO) database and the number of differentially expressed circRNAs were calculated. GO analysis of biological processes showed that differentially expressed genes were enriched in GO: 0000435 positive regulation of transcription from RNA polymerase II promoter by galactose, GO: 0036494 positive regulation of translation initiation in response to endoplasmic reticulum, and GO: 0043666 regulation of phosphoprotein phosphatase activity. GO analysis of cell types showed that differentially expressed genes were enriched in GO: 0030891 VCB complex, GO: 0042582 azurophil granule, GO: 0044666 MLL3/4 complex, and GO: 0032937 SREBP-SCAP-Insig complex. GO analysis of molecular function showed that differentially expressed genes were enriched in GO: 0008440 inositol-1,4,5-trisphosphate 3-kinase activity, GO: 0004416 hydroxyacyl glutathione hydrolase activity, GO: 0043813 phosphatidylinositol-3,5-phosphate 5-phosphatase activity, GO: 0043812 phosphatidylinositol-4-phosphate phosphatase activity, GO: 0047453 ATP-dependent NAD(P)H-hydrate dehydratase activity, GO: 0052855 ADP-dependent NAD(P) H-hydrate dehydratase activity and GO: 0072542 protein phosphatase activator activity.

Enrichment analysis

The Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis showed that the differentially expressed circRNAs were mainly concentrated in 16 signaling pathways that included the thyroid hormone signaling pathway, the phosphatidylinositol 3′-kinase (PI3K)-Akt signaling pathway, and the adenosine 5′-monophosphate (AMP)-activated protein kinase signaling pathway (Figure 4).

Figure 4

Enrichment analysis of differentially expressed genes of the Kyoto Encyclopedia of Genes and Genomes (KEGG) signal pathways. The ordinate is the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway entry, and the abscissa is the ratio of the number of differentially enriched genes to the number of genes annotated on the pathway. The size of the dot in the figure represents the number of differentially enriched genes and the number of all genes annotated. The intensity of the color represents the prominence of the pathway.

Quantitative real-time polymerase chain reaction (qRT-PCR)

The results from qRT-PCR showed that the expression of hsacirc_007293 and hsacirc_031752 genes were upregulated in the patients with papillary thyroid carcinoma, and the expression of the hsacirc_020135 gene was upregulated (Figure 5). These data were consistent with the results from high-throughput sequencing.

Figure 5

Quantitative real-time polymerase chain reaction (qRT-PCR) verified the differential circular RNAs (circRNAs). Compared with the control group, * P<0.05.

Discussion

In this study, alterations of exosome circular RNAs (circRNAs) were identified from the serum of patients with papillary carcinoma of the thyroid using high-throughput sequencing. Compared with the control group, who were patients diagnosed with benign goiter, the serum-derived exosomes contained three upregulated circRNAs and 19 down-regulated circRNAs in the patients with papillary carcinoma of the thyroid. The differentially expressed circRNAs were found to be related to 16 signaling pathways.

The circRNAs have a role as miRNA sponges to regulate the expression levels of other related RNAs, by binding to miRNAs to block their inhibitory actions on the expression of their target genes [28]. Alterations in circRNA expression have also been identified during the onset of renal injury [29], and in gastrointestinal diseases [30]. Reduced levels of circRNAs have also previously been reported in tumor tissues from patients with colorectal cancer when compared with adjacent normal colorectal tissue [31]. Also, overexpression of circTCF25 (hsa_circ_0041103) has been shown to down-regulate miR-103a-3p and miR-107, increase the expression of cyclin-dependent kinase 6 and to promote cell proliferation and migration of bladder cancer cells [32]. Also, the findings from a recent study have shown that circRNAs are enriched and stable in exosomes and can be monitored to diagnose traumatic brain injury [33]. In the present study, three exosome circRNAs were upregulated and 19 were down-regulated in the serum of patients with papillary carcinoma of the thyroid.

The formation of malignant tumors includes the activation of proto-oncogenes and/or the inactivation of tumor suppressor genes, together with changes in apoptosis-regulating genes and DNA repair genes, resulting in cell transformation that leads to the development of malignant tumors [34]. The development of malignant tumors usually involves multiple signaling pathways. Enrichment analysis of the Kyoto Encyclopedia of Genes and Genomes (KEGG) signaling pathway is based on the whole genome, and by counting the number of differentially expressed genes (DEGs) at different levels, the metabolic pathways and signaling pathways that involve DEGs can be determined. In the present study, KEGG analysis showed that the DEGs were mainly concentrated in 16 signaling pathways, including the thyroid hormone signaling pathway, the PI3K-Akt signaling pathway, the AMPK signaling pathway, and the phosphatidylinositol signaling pathway. Both the PI3K-Akt signaling pathway and the AMPK signaling pathway have been shown to be involved in the development of thyroid carcinoma [35-37]. The data from the present study showed that serum exosome circRNAs might be potential diagnostic molecular biomarkers for papillary thyroid carcinoma.

In this study, interaction network analysis between circRNAs and miRNA was performed. To determine the reliability of the high-throughput sequencing data, the expression of hsacirc_007293, hsacirc_031752, and hsacirc_020135 were found to be consistent with the high-throughput sequencing results by fluorescence quantitative real-time polymerase chain reaction (qRT-PCR). These results suggest that the expression of circRNA in the exosomes of patients with simple goiter differed from that of patients with papillary thyroid cancer, indicating that circRNAs may play a role in the development of papillary thyroid carcinoma.

This study had several limitations in this study. Although circRNAs were shown to be abnormally regulated in serum exosomes in patients with papillary thyroid carcinoma, these preliminary findings require support from further in vitro and in vivo studies. Controlled clinical studies should be designed to determine whether circRNAs might have a role as diagnostic, prognostic, or therapeutic biomarkers in patients with papillary thyroid carcinoma. Also, although the findings of the present study are supported by the findings from studies that have shown that circRNAs have previously been reported to promote the progression of papillary thyroid carcinoma [38,39], the signaling pathways regulated by these circRNAs remain to be confirmed.

Conclusions

The findings from this study have shown that gene regulation can be studied from exosomes obtained from serum samples in patients with papillary thyroid carcinoma, and supports the need for further studies on the role of exosome circular RNAs (circRNAs) in thyroid cancer. This study screened out some differentially expressed circRNAs using high-throughput sequencing technology and enriched them into tumor-related signaling pathways. The potential functions of these differentially regulated circRNAs as diagnostic, prognostic, or therapeutic biomarkers should be further investigated.