The Involvement of the hsa_circ_0088494-miR-876-3p-CTNNB1/CCND1 Axis in Carcinogenesis and Progression of Papillary Thyroid Carcinoma

Recently, growing studies have demonstrated that circular RNAs (circRNAs) function as critical players in multiple human tumors, including papillary thyroid carcinoma (PTC). However, the expression and underlying potential mechanism of circRNAs in PTC are still not fully elucidated. In this study, 14 candidate differentially expressed circRNAs (DECs) between normal thyroid tissues and benign thyroid tissues or PTC were first screened using the GSE93522 dataset by the GEO2R online tool. Then, the structural loop graphs of these 14 circRNAs were obtained through the CSCD database. After performing miRNA co-prediction by combination of CSCD and CRI databases, a potential circRNA-miRNA sub-network, consisting of 9 circRNAs and 21 miRNAs, was successfully constructed. Subsequently, the expression and prognostic values of these miRNAs were further determined by starBase, and two miRNAs, namely, miR-605-5p and miR-876-3p, were identified as key miRNAs in PTC. Then, their downstream target genes were predicted by the miRNet database. CTNNB1 and CCND1 were found to be two most potential targets of miR-876-3p by combination of multiple in silico analyses, including protein–protein interaction (PPI), hub gene screening, correlation analysis, and expression analysis. Conclusively, we established a key hsa_circ_0088494-miR-876-3p-CTNNB1/CCND1 axis linked to carcinogenesis and progression of PTC, which may provide promising therapeutic targets in treating PTC in the future.


INTRODUCTION
Papillary thyroid carcinoma (PTC) is one of the most prevalent malignant tumors all over the world, which accounts for more than 80% of thyroid cancers (Lim et al., 2017). Among all malignant cancer types, PTC's prognosis is generally good and the 10-year survival rate of PTC has reached about 90% due to its low rate of tumor growth and high degree of tumor differentiation FIGURE 1 | Identification of candidate circRNAs associated with progression of papillary thyroid carcinoma (PTC). (A) The volcano plot of differentially expressed circRNAs (DECs) between 6 benign thyroid lesions and 6 normal thyroid tissues. (B) The volcano plot of differentially expressed circRNAs (DECs) between 6 PTC and 6 normal thyroid tissues. The red dots and green dots respectively represent upregulated DECs and downregulated DECs with significance (P < 0.05 and |log 2 FC| > 1). The black dots are those DECs without significance. (C) The interaction analysis for upregulated DECs in both compared groups. (D) The interaction analysis for downregulated DECs in both compared groups. (Gui et al., 2020). However, the total survival rate of PTC patients with some aggressive features, such as lymph node and distant metastatic potentials, is still unsatisfactory (Xue et al., 2020). Thus, it is of great significance to uncover the underlying molecular mechanism of PTC carcinogenesis and progression, which may provide novel therapeutic targets for PTC treatment.
Circular RNAs (circRNAs) are a class of endogenous noncoding RNAs (ncRNAs) with covalently closed loops (Qu et al., 2015). Increasing evidences have showed that deregulation of circRNAs widely appears in a series of human diseases and dysregulated circRNAs play key roles in pathogenesis of these diseases, especially cancers (Ding et al., 2020;Guo et al., 2020;. In PTC, several circRNAs have been validated to be involved in development of PTC. For example, Wang et al. (2018) suggested that circ-ITCH inhibited progression of PTC through the miR-22-3p/CBL/beta-catenin pathway; Bi et al. (2018) showed that circRNA_102171 promoted the progression of PTC by regulating CTNNB1P1-dependent activation of the beta-catenin pathway; Liu et al. (2018) also indicated that hsa_circ_0060060 enhanced cisplatin resistance of human thyroid carcinoma cells by autophagy regulation. However, the current knowledge regarding circRNAs in PTC is still not enough, and the expression profile and detailed mechanism in PTC need to be further explored.
In this study, we first identified differentially expressed circRNAs (DECs) between normal thyroid tissues and benign thyroid tissues (normal vs. benign comparison) or PTC (normal vs. PTC comparison). Those identified DECs that were only significantly upregulated or downregulated in normal vs. PTC comparison were screened out and were regarded as candidate circRNAs in PTC. Next, a potential hsa_circRNA/miRNA/mRNA regulatory axis linked to PTC carcinogenesis and progression was established through a series of in silico analyses, including miRNA co-prediction,

Selection of Potential circRNAs in PTC
To find some potential circRNAs associated with carcinogenesis and progression of PTC, the GSE93522 dataset from the GEO database was finally used by a series of criteria selection, after performing differential expression analysis between normal thyroid tissues and benign thyroid tissues or PTC using the GEO2R tool. As shown in Figure 1A and Table 1, 29 and 7 circRNAs were significantly upregulated and downregulated in benign thyroid tissues when compared with normal thyroid tissues, respectively. Moreover, a total of 19 significant DECs, consisting of 18 upregulated and 1 downregulated DECs, were discovered in normal vs. PTC comparison ( Figure 1B and Table 2). In order to improve the analytic accuracy, integration analysis between normal vs. PTC and normal vs. benign comparisons was carried out. We screened out those circRNAs that were only significantly differentially expressed in normal vs. PTC comparison (not in normal vs. benign comparison). Consequently, as presented in Figures 1C,D, we found that 13 upregulated (red box) and 1 downregulated (green box) circRNAs were identified. These 14 candidate circRNAs and their detailed information, including circBase ID, parental gene, and genome location, are listed in Table 3. To further understand the 14 circRNAs, the structural loop graphs were described by the CSCD database. Finally, 11 of 14 candidate circRNAs were accessed (Figure 2), including hsa_circ_0003892 ( Figure 2K). All of them were upregulated circRNAs in PTC. The other 3 circRNAs were not available in the CSCD database. Therefore, the 11 circRNAs were regarded as the final candidate circRNAs in PTC.

Prediction and Analysis of Binding miRNAs of circRNAs in PTC
Our team and other groups have previously documented that circRNAs can act as miRNA sponges. Thus, we predicted the binding miRNAs of 11 candidate circRNAs. Two online prediction databases, containing CSCD and CRI, were employed to predict binding miRNAs of circRNAs (data were not shown). Only these miRNAs that were commonly appeared in both CSCD miRNA set and CRI miRNA set were selected for the following investigation, and these miRNAs were considered as candidate miRNAs. Finally, 9 of 11 circRNAs had candidate miRNAs. To better visualization, a potential circRNA-miRNA network, consisting of 9 circRNAs and 21 miRNAs, was established as presented in Figure 3. Based on the ceRNA hypothesis, the binding miRNAs of upregulated circRNAs should play opposite effects of miRNAs in PTC. Subsequently, the expression of 21 miRNAs in thyroid carcinoma was analyzed using starBase (Figure 4). The result suggested that 8 miRNAs (miR-1179, miR-335-5p, miR-605-5p, miR-145-5p, miR-432-5p, miR-876-3p, miR-545-3p, and miR-532-3p) were significantly downregulated in thyroid carcinoma compared with normal thyroid tissues. Moreover, the prognostic values of these 21 miRNAs in thyroid carcinoma were also assessed through the Kaplan-Meier plotter (Figure 5). Among the 21 miRNAs, only thyroid carcinoma patients with high expression of 8 miRNAs (miR-637, miR-605-5p, miR-661, miR-1272, miR-604, miR-876-3p, miR-1184, and miR-646) had favorable prognosis. As shown in Figure 6A, among these miRNAs, only two miRNAs (miR-605-5p and miR-876-3p) were significantly downregulated in thyroid carcinoma and were positively correlated with survival time of patients with thyroid carcinoma. Collectively, miR-605-5p and miR-876-3p might be two most potential miRNAs of candidate circRNAs in PTC.

Prediction and Analysis for Target Genes of miRNAs in PTC
Next, the miRNet database was used to predict downstream target genes of miR-605-5p and miR-876-3p, and 687 targets were finally found as listed in Supplementary Table S1. To better understand the interaction relationship among all these targets, PPI network analysis was performed with the help of the STRING database. As presented in Supplementary Table S2, 3459 PPI pairs were obtained. According to node degree, the top 20 hub genes (TP53, GAPDH, CTNNB1, ACTB, EEF2, CCND1, VEGFA, JUN, HNRNPK, XPO1, MDM2, HNRNPA2B1, NCBP2, DHX9, DICER1, HNRNPD, NCL, DDX5, YBX1, and PCBP2) were screened, with TP53 being the highest node degree, which was 144 ( Figure 6B). Then, GO functional annotation for the target genes of miR-605-5p and miR-876-3p was also conducted. As shown in Figures 6C-E, the enriched GO items for these target genes included regulation of the metabolic process, regulation of the primary metabolic process, and regulation of the macromolecule metabolic process in the biological process (BP) category; nucleus, nuclear part, and nucleoplasm in the cellular component (CC) category; and heterocyclic compound binding, binding, and organic cyclic compound in the molecular function (MF) category. Furthermore, KEGG pathway enrichment analysis was also conducted. The top 5 enriched KEGG pathways were microRNAs in cancer, proteoglycans in cancer, spliceosome, viral carcinogenesis, and focal adhesion ( Figure 6F). Construction of a Potential hsa_circ_0088494-miR-876-3p-CTNNB1/CCND1 Axis in PTC Subsequently, a miRNA-hub gene network, consisting of miR-605-5p, miR-876-3p, and the top 20 hub genes, was established using Cytoscape ( Figure 7A). The expression correlation between miRNAs and their respective target genes in thyroid carcinoma were evaluated by the starBase database ( Table 4). As presented in Figures 7B-D, only three miRNA-gene pairs (miR-876-3p-CTNNB1, miR-876-3p-ACTB, and miR-876-3p-CCND1) revealed a statistically negative expression relationship in thyroid carcinoma. Next, the expression levels of CTNNB1, ACTB, and FIGURE 3 | The potential binding miRNAs of 9 circRNAs predicted by CSCD and circular RNA interactome (CRI) databases.
CCND1 in thyroid carcinoma were also determined (Figures 7E-G). Among the three genes, only CTNNB1 and CCND1 expression was remarkably increased in thyroid carcinoma tissues when compared with normal thyroid tissues. These findings suggested that CTNNB1 and CCND1 might be two downstream key targets of miR-876-3p in PTC. Taken together, dysregulation of a potential pathway, namely, hsa_circ_0088494-miR-876-3p-CTNNB1/CCND1, might play important roles in the process of PTC carcinogenesis and development ( Figure 8A). Finally, 25 pairs of clinical PTC tissues and adjacent normal tissues were employed to preliminarily validate RNA expression and RNA-RNA correlation of the established potential pathway in PTC. As shown in Figures 8B-E, expression levels of hsa_circ_0088494, CTNNB1, and CCND1 were significantly increased but miR-876-3p was markedly downregulated in PTC tissues when compared with normal tissues. These findings were identical with our previous analytic results. Furthermore, despite that no statistical expression correlation was observed between the miR-876-3p/CTNNB1 RNA pair, miR-876-3p was negatively correlated with hsa_circ_0088494, CTNNB1, or CCND1 and hsa_circ_0088494 was positively linked to CTNNB1 or CCND1 in PTC (Figures 8F-J), which was in accordance with competing endogenous RNA (ceRNA) hypothesis (Salmena et al., 2011).

DISCUSSION
Papillary thyroid carcinoma is one of the most common human cancer types worldwide. Despite that the overall outcome of patients with PTC is good, the detailed molecular mechanism of PTC pathogenesis is not fully elucidated and needs to be investigated to further improve prognosis of patients with PTC.
To uncover the expression and mechanism of circRNA in PTC, GEO2R was used to perform differential expression analysis for the GSE93522 dataset from the GEO database. A total of 14 candidate circRNAs were identified, which were only significantly dysregulated in PTC. Some of these "p < 0.05" represents significant difference, and "ns" represents no significance.   circRNAs have been documented to be closely associated with initiation and progression of human cancers. For example, hsa_circ_0001955 was linked to development of colorectal cancer and breast cancer (Afzali and Salimi, 2019;Ding et al., 2020). Salmena et al. (2011) proposed that ncRNAs, including circRNAs, could "talk" to messenger RNAs (mRNAs) using microRNA response elements (MREs). Recent studies have also experimentally validated that circRNAs could modulate gene expression through sponging shared miRNAs, thus influencing PTC progression. For example, Wei et al. (2018) indicated that circZFR facilitated cell proliferation and invasion of PTC by upregulating c8orf4 through sponging miR-1261; Li et al. (2018) suggested that circNUP214 sponged miR-145 to enhance ZEB2 expression in PTC cells.
Thus, the binding miRNAs of these candidate circRNAs were predicted by combination of two databases, namely, CSCD and CRI. Finally, 21 miRNAs targeting 9 circRNAs were found. The 9 circRNAs were all significantly upregulated in PTC when compared with normal thyroid tissues. In accordance with ceRNA hypothesis, there is a negative relationship between circRNA and its binding miRNA, indicating that miRNAs should be tumor-suppressive modulators in PTC. Therefore, expression analysis and survival analysis were performed to improve the analytic accuracy. Consequently, two miRNAs, miR-605-5p and miR-876-3p, were identified as the most potential binding miRNAs of candidate circRNAs. Despite that the roles of the two miRNAs in PTC have not been studied, miR-605-5p and miR-876-3p have been found to act as tumor-suppressor genes in human cancers. For example, miR-605-5p was markedly downregulated in melanoma tissues and cells and inhibited melanoma progression and glutamine catabolism through targeting GLS ; miR-876-3p functioned as a tumor suppressor and correlated with cell metastasis in pancreatic adenocarcinoma via targeting JAG2 .
To understand the downstream action mechanism of miR-605-5p and miR-876-3p, their target genes were predicted using a comprehensive target gene prediction database, namely, miRNet. It has been widely acknowledged that genes usually exert their functions by interacting with each other. Thus, these predicted target genes were entered into the STRING database to conduct a PPI network analysis, and the top 20 hub genes in this PPI network were screened, after which a miRNA-hub gene network was established. Based on the action mechanism of miRNA, there should be a negative expression relationship between miRNAs and target genes, and target genes of tumor suppressive miRNAs should function as oncogenes. At the end, CTNNB1 and CCND1 were screened out, which were both the potential targets of miR-876-3p. CTNNB1 is part of a complex of proteins that constitute adherens junctions, which is correlated with cancer pathogenesis (Ilina et al., 2020;Ye et al., 2020;Zhang et al., 2020;Zhitnyak et al., 2020). CCND1 belongs to the highly conserved cyclin family, whose members are featured by a dramatic periodicity in protein abundance throughout the cell cycle. Recent studies have also suggested that CCND1 is involved in the process of cancer initiation and progression in multiple cancers, including breast cancer (Montaudon et al., 2020), lung cancer (Yang Y. et al., 2020), and bladder cancer (Chen et al., 2019;Ying et al., 2020). Besides, CCND1 and CTNNB1 were also reported to function as oncogenes in thyroid cancer (Jung et al., 2009;Guo et al., 2019;Wong et al., 2019;Li et al., 2020). The previous analytic results together with these reports indicated that hsa_circ_0088494-miR-876-3p might influence CTNNB1 and CCND1 expression and function, thereby exerting their roles in mediating cancer carcinogenesis and progression.
In conclusion, our study elucidated a potential pathway, namely, hsa_circ_0088494-miR-876-3p-CTNNB1/CCND1 axis, in PTC by combination of a series of bioinformatic analyses, which may provide key clues for developing effective therapeutic targets in treating patients with PTC. However, the current findings need to be further validated by basic experiments and clinical trials in the future.

Inclusion of Dataset
In this study, we aimed to identify some potential functional circRNAs in initiation and progression of PTC. The NCBI GEO database 1 was used to screen possible datasets. Selection criteria are as follows: (1) The selected datasets should contain normal thyroid tissues, benign thyroid lesions, and PTC tissues; (2) the sample count of selected datasets should be more than 10; and (3) datasets about human PTC cells or animal PTC tissues should be excluded. Finally, we found that only GSE93522 met all these criteria. GSE93522 contained a total of 18 thyroid samples, consisting of six normal thyroid tissues, six benign thyroid lesions, and six PTC tumors. GSE93522, based on the platform of GPL19978 Agilent-069978 Arrarystar Human CircRNA microarray V1, studied the circRNA expression profile in different thyroid tissues.

Differential Analysis
To obtain the differentially expressed circRNAs (DECs) associated with carcinogenesis and progression of PTC, differential analysis was performed using GEO2R 2 which is an online tool provided by the NCBI GEO database as previously described (Ding et al., 2020). | log 2 FC| > 1 and P < 0.05 were set as the thresholds for identifying DECs. circBase Analysis circBase 3 is a database for investigating circRNAassociated studies, which provides scripts to identify was introduced to acquire location and parental genes of candidate circRNAs.

Cancer-Specific circRNA Database (CSCD) Analysis
Cancer-specific circRNA database 4 is a database for cancerspecific circular RNAs, contributing to the research for the function and regulation of cancer-associated circRNAs . The CSCD database was employed to obtain structural loop graphs of candidate circRNAs and predict the potential microRNA response element (MRE) sites for each candidate circRNAs.

Circular RNA Interactome (CRI) Analysis
Circular RNA Interactome 5 , a web tool for exploring circular RNAs and their interacting proteins and microRNAs, was also introduced to predict binding miRNAs of candidate circRNAs (Dudekula et al., 2016). The predicted miRNAs that potentially bind to candidate circRNAs were directly exported from the CRI database.
starBase Analysis starBase database 6 is an open-source platform for decoding miRNA-ceRNA, miRNA-ncRNA, and protein-RNA interaction networks from 108 CLIP-seq data sets generated from 37 independent studies (Yang et al., 2011;Li et al., 2014). In this study, starBase was utilized to analyze the expression levels of miRNAs and target genes in thyroid carcinoma. Moreover, this database was used to assess the expression correlation of miRNAs with their corresponding target genes in thyroid carcinoma. R < -0.1 and P < 0.05 were considered as statistically significant.

Kaplan-Meier Plotter Analysis
As we previously described (Lou et al., 2019a,b), the prognostic values of miRNAs in thyroid carcinoma were evaluated by the Kaplan-Meier plotter database 7 , which is an online tool capable of accessing the effect of 54,000 genes on survival in 21 cancer types, and its miRNA subsystems include 11,000 samples from 20 distinct cancer types, containing thyroid carcinoma.