Nine Genes Mediate the Therapeutic Effects of Iodine-131 Radiotherapy in Thyroid Carcinoma Patients

Background Thyroid carcinoma (THCA) is one of the most common malignancies of the endocrine system, which is usually treated by surgery combined with iodine-131 (I131) radiotherapy. Aims This study is aimed at exploring the potential targets of I131 radiotherapy in THCA. Methods The RNA-sequencing data of THCA in The Cancer Genome Atlas database (including 568 THCA samples) was downloaded. The differentially expressed genes (DEGs) between the tumour samples whether or not subjected to I131 radiotherapy were identified using edgeR package. Using the WGCNA package, the module that was relevant with I131 radiotherapy was selected. The intersection genes of the hub module nodes and the DEGs were obtained as hub genes, followed by the function and pathway enrichment analyses using the clusterProfiler package. Moreover, the protein-protein interaction (PPI) network for the hub genes was constructed using Cytoscape software. In addition, more important hub genes were analysed with function mining using the GenCLiP2 online tool. The qPCR analysis was used to verify the mRNA expression of more important hub genes in THCA tissues. Results There were 500 DEGs (167 upregulated and 333 downregulated) between the two groups. WGCNA analysis showed that the green module (428 nodes) exhibited the most significant correlation with I131 radiotherapy. A PPI network was built after the identification of 53 hub genes. In the PPI network, CDH5, KDR, CD34, FLT4, EMCN, FLT1, ROBO4, PTPRB, and CD93 exhibited higher degrees, which were mainly implicated in the vascular function. The relative expression of nine mRNAs in the THCA tissues treated with I131 was lower. Conclusion I131 radiotherapy might exert therapeutic effects by targeting CDH5, KDR, CD34, FLT4, EMCN, FLT1, ROBO4, PTPRB, and CD93 in THCA patients.


Introduction
As one of the most common malignancies of the endocrine system, thyroid carcinoma (THCA) is classified into papillary (PTC), follicular (FTC), medullary (MTC), and undifferentiated types [1]. Although THCA only accounts for 2% of all malignant tumours, it has attracted attention due to a younger target population and higher incidence [2]. In particular, anaplastic THCA usually has a poor outcome and is responsible for 14-39% of THCA-related deaths [3,4].
There have been several guidelines for the diagnosis and initial disease management for the treatment of thyroid car-cinoma [5]. The most recent American Thyroid Association guidelines (2009) recommend total thyroidectomy for tumours > 1 cm and possible lobectomy for tumours ≤ 1 cm [6]. Nevertheless, when THCA invades a relatively wider area, the tumours sometimes cannot be completely removed surgically [7]. Therefore, some patients need to be treated with iodine-131 (I 131 ) radiotherapy postoperatively to destroy any remaining neoplastic cells that were not removed during surgery [8][9][10]. I 131 is a radioactive isotope of iodine, which can kill human tissue cells [11]. To avoid environmental and human contamination, radioiodine treatment usually uses an I 131 dose of 100-200 mCi with the patient in total isolation for 48 h. This therapy has shown positive results in terms of reducing recurrence rates in patients with welldifferentiated THCA [12]. However, the molecular mechanisms underlying the effects of I 131 -related radiotherapy on THCA remain poorly understood.
In this study, we intended to explore the key genes that respond to I 131 radiotherapy. RNA-sequencing data of THCA patients whether or not subjected to I 131 radiotherapy were downloaded from a public database, followed by a series of bioinformatics analyses. The qPCR analysis was used to verify the mRNA expression of hub genes in THCA tissues.

Materials and Methods
2.1. Data Source. The University of California Santa Cruz's (UCSC's) Xena platform (http://xena.ucsc.edu/) contains the public data from databases, such as The Cancer Genome Atlas (TCGA, https://cancergenome.nih.gov/), International Cancer Genome Consortium (http://www.icgc.org), and Therapeutically Applicable Research to Generate Effective Treatments (http://target.cancer.gov). The data on this platform was standardized for subsequent analysis [13]. The RNA-sequencing data and relevant clinical data of THCA in the TCGA database (including 568 THCA samples and 615 phenotypes) were downloaded using the UCSC Xena platform [13].

Data Preprocessing.
Using the R package biomaRt (version 2.38.0, http://bioconductor.org/packages/biomaRt/) [14], the serial numbers of genes were converted into gene symbols. The serial numbers that could not be matched to gene symbols were deleted. When multiple serial numbers corresponded to the same gene symbol, the average value of the serial numbers was taken as the final expression value of this gene symbol and was used for subsequent analysis.
The clinical expression matrices that matched to the gene symbols were extracted according to gene types, and only mRNAs were extracted for subsequent analysis. The mRNAs with low expression were difficult to be verified and could interfere with the results of differential expression analysis, and thus, these mRNAs (expressing in at least 25% samples) were removed. The remaining genes were included in the subsequent analysis.
Normal samples were removed from the clinical expression profile of THCA. From the remaining tumour samples, the samples with or without clear clinical information of receiving I 131 radiotherapy were screened and the samples without such records were deleted. The eligible tumour samples were used for the subsequent analysis.

Differential Expression Analysis.
To analyse the gene expression differences between the tumour samples whether or not subjected to I 131 radiotherapy, differential expression analysis was conducted using the R package, edgeR (version 3.24.3, http://bioconductor.org/packages/edgeR/) [15]. The threshold for screening differentially expressed genes (DEGs) was set at p value ≤ 0.01.

Weighted Gene Coexpression Network Analysis
(WGCNA). WGCNA is an algorithm that can be used to cluster genes with similar expression patterns and construct a gene coexpression network [16]. First, the WGCNA algorithm assumed that the gene network obeys a scale-free distribution and defined the correlation matrix of the gene coexpression and the adjacency function formed by the gene network. Then, the dissimilarity coefficient of different nodes was calculated, and the hierarchical clustering tree was constructed accordingly. Finally, the associations between modules and specific phenotypes or diseases were explored, and the target genes and gene networks for disease treatment were ultimately identified. Using the R package, WGCNA (version 1.66, https://cran.r-project.org/web/ packages/WGCNA/) [16], WGCNA analysis was carried out. In detail, the expression matrix obtained through data preprocessing was subjected to log 2 ðx + 1Þ standardization, and the genes with median absolute deviation in the top 75% and larger than 0.01 were screened. After screening and removing the outlier, the one-step method was used to build the WGCNA network, and each module was visualized in the hierarchical clustering tree. Finally, receiving I 131 treatment or not was used as the external biological parameter for performing correlation analysis of the modules, and the module that was the most highly correlated with the external biological parameter was selected for further study.   Figure 3: The hierarchical clustering tree of WGCNA. The module that was relevant with I 131 radiotherapy in THCA was analysed by using the WGCNA. The expression matrix obtained through data preprocessing was subjected to log 2 ðx + 1Þ standardization, and the genes with median absolute deviation in the top 75% and larger than 0.01 were screened. After screening and removing the outlier, the one-step method was used to build the WGCNA network, and each module was visualized in the hierarchical clustering tree. The table showed the number of clustered genes in the functional modules. The figure showed the hierarchical clustering tree of the weighted gene coexpression network. Grey indicates a cluster of genes that are not included in any of the functional modules.
3 Disease Markers 2.5. Functional and Pathway Enrichment Analyses. After the module that was most highly correlated to I 131 treatment was identified, the module nodes were subjected to Gene Ontology (GO) [17] function and Kyoto Encyclopedia of Genes and Genomes (KEGG) [18] pathway enrichment analyses using the clusterProfiler package (version 3.10.1, http:// bioconductor.org/packages/clusterProfiler/) [19] in R. The thresholds for screening the enrichment results were defined at adjusted p value < 0.01 and q value < 0.05.
Using Cytoscape software (http://www.cytoscape.org/) [20], the key module was analysed and the degree of each module node was calculated. The module nodes with degree > 200 were selected as hub nodes. By comparing the hub nodes with the DEGs, the intersection genes that were statistically significant and played key roles in the module were identified as hub genes. Furthermore, enrichment analysis for the hub genes was conducted using the clusterProfiler package [19]. The screening threshold was set at adjusted p value < 0.05.

Protein-Protein Interaction (PPI) Network Analysis.
Using the STRING database (http://string-db.org/) [21], PPI analysis for the hub genes was conducted. The threshold for PPI analysis used the default values in the database. The obtained PPI pairs were visualized in the PPI network using Cytoscape software [20], and the network nodes with degree > 5 were selected as more important hub genes.

Function
Mining of More Important Hub Genes. Using the GenCLiP2 online tool (http://ci.smu.edu.cn/GenCLiP2/ analysis.php) [22], the functions of the more important hub genes were analysed through literature mining. The enrichment results with p value < 1e-4 and involving at least four genes were selected.

Collection of Tissue Samples.
The study protocol was approved by the Medical Ethics Committee of Huzhou Central Hospital (China). Informed consent was obtained from patients or guardians. In total, 47 THCA patients subjected to I 131 radiotherapy and 50 THCA patients without preoperative therapy at the Huzhou Central Hospital from March 2018 to March 2019 were included in the present study. The pathological diagnosis of THCAs was confirmed by two senior pathologists. The collected tissue samples were snap-frozen and stored at -80°C.
2.9. qPCR Analysis. A TRIzol reagent (Invitrogen, Carlsbad, USA) was used to isolate the total RNA. The First-Strand cDNA Synthesis kit (K1622; Thermo Fisher, USA) was used to synthesize the first-strand cDNA. Quantitative real-time

Data Preprocessing and Differential Expression Analysis.
After data preprocessing, a total of 17003 genes and 298 tumour samples were selected. Therefore, a 17003 × 298 expression matrix was obtained for the subsequent analysis.
There were 500 DEGs between the tumour samples whether or not subjected to I 131 radiotherapy, including 167 upregulated genes and 333 downregulated genes. The volcano plot of the DEGs is presented in Figure 1.

WGCNA Analysis.
After preprocessing, a 12751 × 298 matrix was obtained for WGCNA analysis. In the process of outlier screening, one outlier was found in all samples. After removing the outlier, the 12751 × 297 matrix was used for constructing the WGCNA network. In WGCNA analysis, the soft threshold (power) was selected as 12 ( Figure 2(a)). The mean connectivity tended to zero as the soft threshold got larger, indicating that the calculation of the soft threshold was reliable (Figure 2(b)). The weighted gene coexpression network was built and visualized in the hierarchical clustering tree. There were a total of 10 functional modules (black module, 133 nodes; blue module, 906 nodes; brown module, 833 nodes; green module, 428 nodes; magenta module, 54 nodes; pink module, 92 nodes; purple module, 34 nodes; red module, 231 nodes; turquoise module, 5487 nodes; and yellow module, 692 nodes) (Figure 3).

Disease Markers
Both correlation analysis among the functional modules (Figure 4(a)) and correlation analysis between each functional module and phenotype (whether or not subjected to I 131 radiotherapy) (Figure 4(b)) were conducted. Our results showed that the green module exhibited the most significant correlation with I 131 radiotherapy. The genes in the green module were listed in the supplementary material (Table S2). Therefore, the green module was selected for the subsequent analysis. In addition, the correlation between each functional module and gene ( Figure 5(a)) and the correlation between I 131 radiotherapy and gene ( Figure 5(b)) were separately analysed. The genes involved in the green module had the strongest correlations with all genes. These results suggested that the green module was the most relevant and significant functional module for THCA patients subjected to I 131 radiotherapy. The results of the correlation analysis indicated that the genes in the green module were significantly correlated not only with the corresponding module but also with I 131 radiotherapy ( Figure 5(c)). Thus, the genes in the green module deserved further exploration.

Functional and Pathway Enrichment Analyses.
To explore the biological functions and pathways involving the components of the green module, the 428 genes in the green module were subjected to enrichment analysis. The enriched functions mainly included circulatory system processes (GO biological progress (BP), adjusted p value = 1.72E-19, and q value = 6.23E-16), extracellular matrix (GO cellular component (CC), adjusted p value = 1.20E-09, and q value = 4.34E-07), and channel activity (GO molecular function (MF), adjusted p value = 4.03E-07, and q value = 1.24E -04). In addition, the pathway enrichment analysis showed that the 428 genes were mainly involved in focal adhesion (adjusted p value = 7.10E-08 and q value = 1.14E-05), cGMP-PKG signaling pathway (adjusted p value = 1.04E -07 and q value = 1.14E-05), and Rap1 signaling pathway (adjusted p value = 5.92E-07 and q value = 3.33E-05) ( Table 1).
A total of 70 hub nodes with a degree > 200 were identified from the green module and were then compared with the 500 DEGs. Finally, 53 intersection genes were obtained as hub genes. These hub genes were implicated in 15 GO biological progress (BP) functional terms (such as endothelium development, adjusted p value = 1.69E-08), two GO cellular component (CC) functional terms (such as cell-cell junction, adjusted p value = 6.60E-03), four GO molecular function   Table 2).

Function
Mining of More Important Hub Genes. The functions of the nine more important hub genes were mined from the PubMed database. The enrichment analyses showed that these more important hub genes were mainly related to vascular function (such as angiogenesis, cell adhesion, cell activation, and necrosis) (Figure 7).
3.6. qPCR Analysis of Hub Genes in Tissues. The nine more important hub genes were verified in the clinical tissue samples. The relative mRNA expression of nine hub genes in THCA tissues was shown in Figure 8. The experimental group and control group were the 47 THCA samples treated with I 131 and 50 THCA samples not treated with I 131 , respec-tively. The mRNAs of the nine genes including CDH5, KDR, CD34, FLT4, EMCN, FLT1, ROBO4, PTPRB, and CD93 showed a lower expression compared to the control group (p < 0:05, independent sample t-test).

Discussion
Our study, for the first time, investigated the key genes that respond to I 131 radiotherapy in THCA by analysis of the RNA-sequencing data. A total of 500 DEGs (167 upregulated and 333 downregulated genes) between the tumour samples whether or not subjected to I 131 radiotherapy were obtained. Among the 10 functional modules identified by WGCNA, the green module was the most relevant and significant functional module for THCA patients receiving I 131 radiotherapy. After comparing the 70 hub nodes in the green module with the 500 DEGs, 53 intersection genes were selected as hub genes. For the 53 hub genes, 15 GO BP functional terms, two GO CC functional terms, four GO MF functional terms, and three KEGG pathways were enriched. From the PPI network, CDH5, KDR, CD34, FLT4, EMCN, FLT1, ROBO4, PTPRB, and CD93 were selected as more important hub genes, which were mainly implicated in the vascular function. Vascular functions are the basis of normal metabolism of tumour cells [23]. Thus, I 131 radiotherapy might affect the normal metabolism of tumour cells and kill tumour cells by influencing their vascular function.
Snail and N-cadherin exhibit inducible and constitutive expression in PTC and may affect the development and progression of PTC. Therefore, Snail and N-cadherin may serve   The color scale at the top right represents the relative expression of mRNA. The histogram (b) showed that the relative expression of nine mRNAs in the THCA treated with I 131 group was lower than that in the THCA treated without I 131 group, and the differences were statistically significant ( * represents p < 0:05, independent sample t-test).

10
Disease Markers as effective markers for the tumour [24]. N-cadherin contributes to thyroid tumourigenesis by regulating epithelialmesenchymal transition and certain signaling pathways and may be considered a promising therapeutic target for THCA [25]. The plasma membrane protein, ROBO4, serves as a therapeutic target in tumour endothelial cells; ROBO4 conjugates suppress tumour angiogenesis and growth [26,27]. Protein tyrosine phosphatase receptor type J (PTPRJ) has been reported to be a tumour susceptibility gene and acts as a tumour suppressor gene during the carcinogenesis of THCA [28]. These findings suggested that CDH5, ROBO4, and PTPRB might be associated with the mechanisms of THCA. Thus, they may serve as target genes of I 131 radiotherapy. The enrichment analysis for the 428 genes in the green module shows that the hub genes are associated with the vasculature of tumour tissue but not the tumour cells themselves (Table 1). Tumour microenvironment may be involved in the regulation of the iodine-131 radiotherapy in THCH. For example, the HIF1a or VHL-HIF1a axis is targeted by iodine therapy of THCH. The iodine target in tumour cells is HIF1a which could lead to microvessel decrease. Therefore, further studies to evaluate the relationship between the hub genes and HIF1a or VHL-HIF1a axis may provide new research directions for iodine-131 radiotherapy of THCH.
Vascular endothelial growth factor (VEGF), vascular endothelial growth factor receptor 1 (VEGFR1, also called FLT1), and VEGFR2 (also known as KDR) are upregulated in PTC, and inhibiting VEGFR offers a novel approach for managing the development of THCA [29]. Both VEGFR2 and RET protooncogenes in MTC are targeted by the superior dual inhibitor, AGN-PC-0CUK9P, and thus, AGN-PC-0CUK9P may be used for improving the therapy of MTC [30]. The activities of genes in tumour vessels (such as VEGFR2 and VEGFR3 (also called FLT4)) and tumour cells can be suppressed by the multikinase inhibitor, sorafenib, and sorafenib exerts antiangiogenic and antiproliferative effects in advanced THCA [31]. Through a cell-mediated approach, soluble FLT1 has been found to be able to repress tumour angiogenesis and growth in patients with FTC [32]. Therefore, I 131 might also function in THCA treatment by targeting KDR, FLT4, and FLT1.
CD34 is implicated in cell adhesion and signal transduction and is an indicator of neoplastic behaviour, and CD34 + stromal cells are related to the carcinogenesis of PTC [33]. Precursor CD34 + stromal cells or monocytes can develop into dendritic cells (DCs), and a higher number of DCs may contribute to the favourable prognosis of PTC [34]. The transmembrane receptor, CD93, is highly expressed in tumour vessels of multiple cancers, indicating that CD93 functions by organizing the extracellular matrix and mediating vascular maturation [35,36]. The expression of five genes (EMCN, engulfment and cell motility 1, solute carrier organic anion transporter family member 2A1, potassium voltagegated channel subfamily A member regulatory beta subunit 1, and inter-alpha-trypsin inhibitor heavy chain 5) decreased in FTC compared with benign follicular thyroid adenoma (FTA), and the 5-gene classifier has high specificity and sensitivity and can be used to distinguish FTC from FTA [37].
Taken together, the differential expressions of CD34, CD93, and EMCN after I 131 radiotherapy indicated that the three genes may be targets of I 131 radiotherapy.
Although nine key genes were identified to be associated with I 131 radiotherapy in THCA, these genes were not validated in multicenter and large samples, which was a limitation of this study. Whereas I 131 is radioactive and cell experiment is difficult to be carried out, the cell experiments could not be performed. In the follow-up clinical treatment, the weight of genes would be comprehensively considered and a gene screening model based on the nine genes may be constructed to provide guidance for the selection of therapeutic options.

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.

Consent
Written informed consent was obtained from the patients for the publication of this paper.