Integrative Analysis of CircRNA Expression and Associated ceRNA Network in Graves’ Disease

Objective Graves’ disease (GD) is the most common subtype of autoimmune thyroid disease which have the involvement of circular RNAs (circRNAs) in the pathogenesis. However, it is largely unknown about the role of circRNAs in GD. In present study, we performed a comprehensive study of the circRNAs in GD by bioinformatics analyses. Methods CircRNAs were detected in plasma of 3 newly diagnosed GD patients and 3 healthy controls. We selected 2 up-regulated and 1 down-regulated circRNAs with different expression in GD for validation by transcriptase-quantitative polymerase chain reaction in both GD and healthy control subjects. The GTRD base was used for predicting the transcript factors of the different expression of circRNAs. Then the competing endogenous RNAs(ceRNAs) network was assembled and the analysis of the Gene Ontology and Kyoto Encyclopedia of Genes and Genomes were performed to predict the functions of different expression of circRNAs.


Introduction
Graves' disease (GD) is the most common subtype of autoimmune thyroid disease, and manifests as thyrotoxicosis, diffuse goiter and extra-thyroidal physical signs such as ophthalmopathy and localized dermopathy. The annual incidence of this disease has been reported to range from 20 to 50 cases per 100,000 persons [1]. To date, research into the treatment of GD lacks a signi cant breakthrough.
Therefore, it is important to explore the molecular pathogenesis in order to determine optimal biomarkers and molecular targets of GD, which might open new therapeutic avenues for patients with this disease.
Due to the common effects of genetics and environment, organismal immune balance in GD is disrupted, causing the breakdown of immune tolerance and production of thyrotropin receptor antibody (TRAb). It is considered that these might be the core parts of GD's pathogenesis. The overproduction of thyroid is induced by the stimulation of TRAb to the thyroxine receptor (TSHR) located in the thyroid gland, which is the primary cause of the clinical symptoms of GD. However, little is known about the exact operation of this mechanism. Recent study indicates that epigenetic regulation that includes circular RNA (circRNA) plays an important role in the mechanism of GD [2].
Differing from the traditional linear RNA comprised of 5' and 3' ends, circRNAs have characteristic structures formed by the covalent combination of the downstream 3' end and the upstream 5' end 3' [3].
Computer predictions and experiments have shown that under the action of regulatory factors such as transcription factors the parental gene produced circRNA after variable splicing. According to the parental genes, they are divided into four categories: exonic circRNAs (ecirc RNA), intronic circRNAs (ciRNA), exonintron circRNAs (EIci RNA) and intergenic circRNAs or fusion circRNAs (f-circRNAs) [4]. In the eukaryotic transcriptome, circRNAs are characterized by abundant, stable, conserved and tissue-speci c expression [5]. Moreover, compared to circRNAs from traditional tissue, plasma circRNAs as precise body uid specimens have the advantages of being non-invasive and the ability to be applied in real-time [6]. An increasing number of studies have observed multiple functions of circRNAs, but the most noteworthy is the role of microRNA sponges: the competitive endogenous RNA (ceRNA) hypothesis proposes that circRNAs competitively bind to microRNA through microRNA response elements (MERs) to regulate downstream gene expression. It has been reported that circRNAs are involved in regulating the pathogenesis of various autoimmune diseases, such as rheumatoid arthritis (RA) and systemic lupus erythematosus (SLE) [7,8]. Nevertheless, the exploration of circRNAs in the pathogenesis of GD is still in its infancy.
In the present study, circRNA microarray was used to investigate the different expression of plasma circRNAs between GD and healthy control subjects. Through the reverse transcriptase-quantitative polymerase chain reaction (RT-qPCR), the circRNAs with statistically differential expression were screened out and veri ed. By predicting the transcription factors of different expression of circRNA-parental genes, the transcription factors regulating the production of circRNAs could be obtained indirectly. Finally, we used the circRNA-miRNA-mRNA co-expression network to explore the potential role of circRNAs in the pathogenesis of GD.

Participants
The present study consisted of 3 newly diagnosed GD and 3 healthy control subjects whose age and gender matched them. All of them visited the endocrinology of our hospital from March to August 2020.
The 3 newly diagnosed GD were diagnosed using the criteria recommended by American Thyroid Association (ATA) [9]: (i) a palpable goiter and syndrome comprising heat intolerance, increased sweating, accelerated heart, (ii)elevated levels of serum thyrotrophin receptor antibody (TRAb) and free thyroxine estimates (FT4), decreased levels of serum thyrotropin (TSH). The normal values of TRAb, FT4 and TSH were 0-1.75 IU/L, 12-22 pmol/L and 0.27-4.2 mIU/L, respectively. These participants had never been treated for hyperthyroidism or ophthalmopathy.
Exclusion criteria were as follows: hepatic disfunction, renal dysfunction, infection, malignant disease, immune disease, a history of using of immunosuppressive drugs within one year, such as lentinan, mannotide, glucocorticoids, interferon, and cyclophosphamide, and pregnant or lactating women.
People eligible for the study were informed about the research, its purposes and requested to provide written consent of participation. The research protocol was approved by the Local Bioethics Committee.
RNA extraction and quality control A 10ml peripheral venous blood sample was obtained from each participant. Withdrawn samples were centrifuged at 4°C for 10 min at 3000 rpm, which were stored at -80℃until analysis. Total RNAs were extracted from the plasma of all participants by Trizol reagent (Invitrogen life technologies). According to the manufacturer's manuals, the quality of the RNA samples was assayed by a Nanodrop ND-1000 equipment (Thermo Fisher Science, Inc.).

Microarray assay
We used the microarray assay to detect the difference between GD and healthy control subjects. Samples labelling and microarray hybridization were processed according to the standard procedures of Arraystar (Arraystar Inc.). The speci c process was as follows: Firstly, RNAs were puri cated, ampli cated, and transcripted into uorescent cRNA. Secondly, the labeled cRNAs were hybridized onto the assembled expression microarray slide (Arraystar). After cleaning, Agilent G2505C scanner was used to scan the array. Finally, the analysis of images from acquired array was performed by Agilent Feature Extraction software (version 11.0.1.1).
We used the R software LIMMA package (version 3.22.7) (25) and Gene Spring GX v12. 1 (Agilent Technologies) to process the quantile normalization and subsequent data. The different expression of circRNAs between the GD and healthy control subjects were identi ed by fold change ltering. The signi cantly different expression of circRNAs were de ned as those met the standards of fold change > 1.3 and P value < 0.05.

Real-time qPCR validation
The RNAs samples were reverse transcribed into cDNAs for real-time quantitative PCR analysis. We placed the annealing mixture which was consisted of 3.0 µgRNA, 1µl Random N9(0.5ug/ul), 1.6µl dNTPs Mix(2.5mM) and H2O without RNAase in a water bath with a temperature at 65°C for 5 minutes and then put them on ice for 2 minutes.
The reverse transcription reaction system was prepared sequentially, composing 4µl 5X First-Strand Bufferr, 1µl DTT (0.1 M), 0.3µl RNase inhibitor, and 0.2 µl Superscript III RT, and the above mixture. This system underwent incubation in water successively with a temperature of 37˚C for 1 minutes, 50˚C for 60 minutes, and 70˚C for 15 minutes up until the transcription was completed. The PCR reaction mixture was performed according to the following condition: the rst denaturation at 95°C for 10 minutes, followed by 40*PCR cycles (at 95°C for 10 s and 60°C for 60 s), and at the end of the ampli cation reaction, the mixture was processed at 95°C for 10 s ,60°C for 60 s, nally 95°C for 15 s.
Each sample was assayed in triplicate. The Primer sequences were shown in Table 1. We determined the relative expression level of circRNAs using 2-ΔΔCq method. An unpaired t-test was applied to compare the expression of circRNAs between GD patients and healthy control subjects. The statistically signi cance is de ned as a P-value < 0.05.

Transcription factor prediction
The dates of human transcription factor (TF) binding sites were downloaded from GTRD base (http://gtrd.biouml.org/) to predict transcription factors of parental genes and used them as transcription factors for circRNAs. Subsequently, the TFs of circRNAs were further screened according to the following standards: (i) The TFs located in the promoter region of a parent gene (ii) Those veri ed in more experimental data sets were preferred. Cytoscape software (v3.8.1) was used to establish the TF-parental genes-circRNAs network based on the predicted TFs and their downstream targets.
the circRNA-miRNA-mRNA co-expression network construction The circRNA-miRNA-mRNA co-expression network was applied to analyze the different expression of circRNAs which had been veri ed by qPCR. Studies have identi ed that MREs was the foundation of this network [10]. The potential MREs and their downstream mRNAs were predicted by Miranda (http://www.microrna.org/) and TargetScan (http://www.targetscan.org/). The potential microRNAs and gene targets were screened according to the following standards: (i) the veri ed microRNA by literature and MREs 1. (ii) downstream mRNAs related to the pathogenesis of GD. Cytoscape software (v3.8.1) was used to establish the ceRNA network based on data of circRNAs and their predicted downstream targets.

Bioinformatics analysis
Gene ontology (GO) analysis (http://www.geneontology.org) was performed to explore the functions of differentially expressed circRNAs in terms of biological processes (BP), cellular components (CC) and molecular functions (MF). Kyoto Encyclopedia of Genes and Genomes (KEGG) (http://www.genome.jp/keg/) was used to explore those pathways related to circRNAs.

Result
Microarray analysis showed that there were signi cant differences in expression of circRNAs between GD and healthy control subjects.
There were signi cant differences in the expression pro le of circRNAs between GD and healthy control subjects (Fig. 1A). In addition, the volcano plot showed differences in the expression of circRNAs between the two groups ( Compared with the control group, hsa_circ_0090364 was the different expression of circRNA in GD. To further prove the result of microarray analysis, we selected 3 circRNAs including 2 up-regulated circRNAs (hsa_circ_0001228, hsa_circ_0090364) and 2 down-regulated circRNAs (hsa_circ_0008339) randomly from 366 different expression of circRNAs by RT-qPCR ( Table 1). As shown in Fig. 2, only hsa_circ_0090364 was consistent in microarray analysis and RT-qPCR results.

Transcription factor prediction
The circRNA was derived from parental genes, and its transcription factors could be obtained indirectly by predicting of transcription factors of parental genes. JADE3 was the parental gene of hsa_circ_0090364 (Table 1). In this study, we found a total of 3608 transcription factors of JADE3 which included 82 transcription factors in its core promoter region. As shown in Fig. 3, transcription factors were in the promoter region and had been veri ed in more than 3 experimental data sets were acquired. Together with JADE3 and hsa_circ_0090364, a total of 24 of the TFs met criterias that were melted into the network. The transcription factors that we focus on had the function of regulating cell proliferation, differentiation and cell metabolism.
Prediction of the circRNA-miRNA-mRNA co-expression network To explore the regulatory effect of circRNAs in the pathogenesis of GD, we integrated hsa_circ_0090364 into a circRNA-miRNA-mRNA network (Fig. 4A). The network was consisted of 1 circRNA nodes, 5 miRNA nodes and 188 mRNA nodes. As shown in Fig. 4B, the 20 most signi cant hub genes were screened out according to the criterion of node degrees > 5 by cytoHubba and integrated them into a circRNA-miRNA-mRNA network.
The enriched functional analysis in the circRNA-microRNA-mRNA network Based on the circRNA-miRNA-mRNA co-expression network described above, GO functional analysis was performed to predict the potential functions of the effectively different expression of circRNA. Generally, 52 BP, 7 CC, and 12 MF GO terms were predicted to be signi cantly enriched (P < 0.05), and then the top-10 enriched GO terms of BP and MF,7 CC were listed (Fig. 5A). BP showed that hsa_circ_0090364 was involved in MAPK/ERK signaling pathway, epidermal growth factor receptor signaling pathway and phosphatidylinositol-mediated signaling and it has the function of regulating cell proliferation, differentiation and cell metabolism. On the other hand, the GO analysis of CC was indicated that hsa_circ_0090364 was signi cantly linked to the synthesis of IGF-1-IGF1R complex. For molecular function, hsa_circ_0090364 was associated with protein kinase binding, phosphotyrosine binding, transforming growth factor, beta receptor binding, calmodulin binding, calmodulin-dependent protein kinase activity and insulin receptor binding (Fig. 5A). The top 20 KEGG pathways suggested that hsa_circ_0090364 involved in TGF-beta signaling pathway, MAPK signaling pathway, T cell receptor signaling pathway and Toll-like receptor signaling pathway. Hsa_circ_0090364 might be associated with Prostate cancer, Glioma, Small cell lung cancer, Non-small cell lung cancer, Endometrial cancer, Hepatitis B and Viral carcinogenesis. (Fig. 5B) hsa_circ_0090364 was found to interact with miR-378f, miR-378c and miR-378a-3p (Fig. 6). The complementarity between the 1 candidate circRNA and the 3 microRNAs was perfect according to 8mer matching types.

Discussion
Our study indicated that a total of 366 circRNAs were differentially expressed in plasma of GD compared with that of healthy control subjects. Among them, there were 195 circRNAs up-regulated and 171 downregulated, suggesting that circRNAs might play a crucial regulatory role in pathogenesis of GD.
In the microarray analysis of all different expression of circRNAs, we found that exonic circRNAs accounted for the largest proportion (84.3%) which might be related to their function. The exonic circRNA is a loop containing a 3'-5' phosphodiester bond which was described as a covalent link between the 3' end of an exon and the 5' end of either the same exon or the upstream exon. The speci c circle structure protects circRNAs from 'exon skipping' and leads them to become more inclined to regulate the linear coding RNAs from which they are derived [3].
The chromosomal distribution analysis of circRNAs in this study showed that the different expression of circRNAs were mainly localized on chr17, chr2 and chr3. Some studies have indicated that X chromosome disorder was associated with female predominance in adult autoimmune thyroid disease (ATD) [11,12], however a similar result was not re ected in this chromosomal distribution analysis. This difference may be due to the microarray analysis of small samples to lead to the poor reproducibility of the experimental results. It is therefore necessary to perform a high-quality and large-scale study to further explore this point. Ecsedi et al. showed that chr17 imbalance may lead to cell proliferation activation, supporting our ndings [13]. We screened three signi cantly different expressions of circRNAs that included 2 up-regulated circRNAs (hsa_circ_0090364, hsa_circ_0001228) and 1 down-regulated circRNA (hsa_circ_0008339) for RT-qPCR validation. The results showed that hsa_circ_0090364 was consistent in microarray analysis and RT-qPCR, which indicated that the microarray analysis was highly reliable. In summary, circRNAs may have regulatory functions in the epigenetic regulation mechanism of GD, which deserves further attention.
Hsa_circ_0090364 is a high-expression circular exonic RNA in GD, which is mainly generated by the selective splicing of its parental gene Jade3 exon. Therefore, the higher levels of JADE3 in patients implies the potentially elevated expression of hsa_circ_0090364. According to Circbase's annotation, JADE3 was involved in inducing histone acetylation during transcription. Although the relationship between histone acetylation and GD has not been discussed, a study has revealed that the high level of histone H3 acetylation (H3ac) increased the level of IL-6 in Rheumatoid Arthritis Synovial Fibroblast (RASFs) to promote rheumatoid arthritis (RA) [14]. In addition, some studies have identi ed that the high levels of IL-6 were correlated with increased risk for GD [15,16]. The ndings above suggest that JADE3 plays an important role in the pathogenesis of AITD. Until now, the nature of the speci c mechanisms between JADE3 and hsa_circ_0090364 remain unclear. To further investigate the upstream-regulatory mechanism of circRNAs, this study indirectly obtained the transcription factors of hsa_circ_0090364 by predicting the transcription factors of JADE3. As the predicted results show, a total of 3608 transcript factors were predicted. We focused on the 24 transcript factors which have been con rmed by three or more experiments. Functional annotations of these transcript factors suggested that they could promote cell proliferation, differentiation and immune in ammation. Furthermore, the pathogenesis of GD was also related to these biological processes [2]. Given the above studies, we speculated that there was a regulatory relationship among transcript factors, parental gene and circRNAs in the development of GD, with further studies required to illustrate this in detail.
In the ceRNA network of hsa_circ_0090364, we found that hsa_circ_0090364 interacted with multiple miRNAs in the miR-378 family (miR-378h, miR-378f, miR-378b, miR-378c, miR-378a-3p). The 2D structure of microRNA and circRNA (Fig. 6) showed that the binding site of miR-378a-3p, miR-378c and miR-378f on hsa_circ_0090364 had the same microRNA seed sequence CUGGAC, which was able to match perfectly with hsa_circ_0090364. The microRNA seed sequence is de ned as nucleotides 2-7 of the 5' region on the microRNA, and is believed to be the most conserved section of the microRNAs, which is of great importance for miRNA recognition [17]. The same microRNA seed sequence further con rmed that hsa_circ_0090364 and each member of miR-378 family had the same binding ability. Therefore, we presume that hsa_circ_0090364 acts as an 'endogenous sponge' for many microRNAs from the miR-378 family by using the MRE as 'RNA language'.
The miR-378 family has not been reported in studies of GD patients, but was reported by Shi Haizhong et al [18]. It is suggested that the low expression of the miR-378 family enhanced the levels of its downstream molecule IL-33 and played a role in the proin ammatory response of the autoimmune disease ulcerative colitis (UC). At the same time, the investigators screened for miR-378a-3p to verify it. In the present study, hsa_circ_0090364 was detected to be signi cantly up-regulated in GD, and it was established that pathogenesis of GD was the process involved in the autoimmune in ammatory response leading to thyroid cell proliferation resulting in thyroid hormone secretion. Therefore, we hypothesized that hsa_circ_0090364 may act as an endogenous sponge of miR-378h, miR-378f, miR-378b, miR-378c and miR-378a-3p. When hsa_circ_0090364 was up-regulated in GD, miR-378h, miR-378f, miR-378b, miR-378c and miR-378a-3p were repressed at the transcriptional level, thereby regulating the target genes and prompting the occurrence and development of GD. However, this hypothesis needs to be veri ed by future study.
The GO functional analysis of the ceRNA network showed that multiple biological processes of GO enrichment in the network were related to multiple metabolic processes, such as activation of the MAPK pathway, phosphatidylinositol-mediated signal transduction and IGF-1-IGF1R complex synthesis. The MAPK pathway plays a key role in inducing protein synthesis and cell differentiation. Thus, stimulating antibodies in uences thyroid cell activation and proliferation by modulating MAPK signalling cascades in GD. PI3Ks are downstream signal molecules of phosphatidylinositol. The PI3K pathway has the function of regulation of thyroid cell proliferation and hormone synthesis [19]. The IGF-1-IGF1R complex is the core component of the development of thyroid-associated ophthalmopathy. Thyroid stimulating immunoglobulins activate the IGF-1-IGF1R complex through TSHR-IGF1R-cross-talk and increase the level of glycosaminoglycan synthesis and in ammatory molecules. Moreover, immuno-globulins against the IGF1R can activate signals in orbital broblasts, leading to the production of hyaluronic acid and cytokines [20]. In the KEGG pathways analysis of circRNA-miRNA-mRNA network, activating the T cell receptor pathway was one of multiple enriched pathways. In GD, excessive thyroid hormones are secreted by thyroid cells due to being stimulated by thyroid stimulating hormone receptor (TSHR) autoantibodies, which are produced by plasmacyte and local B cells controlled by activated T cells. In this process, activating the T cell receptor pathway is crucial [21,22].
It must be considered that limitations still exist. Firstly, with a greater number of study subjects patient samples used for detection may be increased, and the range of differential expression of circRNAs in GD patients may be narrowed. Secondly, this study con ned its analysis to the different mapping of circRNAs and explored their potency to act as internal competitive endogenous RNAs (ceRNAs), while the functional analysis of circRNAs was not comprehensive. In future studies, more functional validation shall be required.

Conclusions
In this study, the expression pro le of circRNAs in GD was elucidated by microarray analysis. Different expression of circRNAs was found in plasma of GD compared with normal controls. A transcription factor-parental gene circRNA network was structured to explore the upstream regulatory mech-anisms of circRNAs. The ceRNA network was constructed to explore downstream microRNAs and target genes. GO annotation and KEGG pathway analysis were used to clarify the potential of circRNAs to play a regulatory role in the pathogenesis of GD.

Declarations
Funding sources  The Transcription factors-JADE3-hsa_circ_0090364 co-expression network.