Clinical Significance and Potential Regulatory Mechanisms of Serum Response Factor in 1118 Cases of Thyroid Cancer Based on Gene Chip and RNA-Sequencing Data.

BACKGROUND Thyroid cancer (TC) is one of the most prevalent endocrine malignancies and there may be many unclarified molecular events and gene types involved in TC. The objective of this study was to assess the clinical implications and potential mechanisms of serum response factor (SRF) in TC. MATERIAL AND METHODS RNA-sequencing and gene chip data with TC expression were collected from The Cancer Genome Atlas/Genotype-Tissue Expression, Gene Expression Omnibus, ArrayExpress, Sequence Read Archive, and Oncomine. SRF expression of all TC and adjacent non-cancerous tissue were calculated using the t test, STATA, and Meta-DiSc. The related pathways of the potential SRF target genes and target miRNAs were explored. Dual-luciferase reporter assay was performed to validate the association between SRF and its putative miRNA. RESULTS One RNA-sequencing and 15 gene chips were collected, and the pooled standardized mean difference of SRF was -1.00. Furthermore, the area under the curve of sROC of SRF in TC was 0.8251, indicating a dramatic decreased expression of SRF in TC tissues based on 1118 cases. The intersection of differentially expressed genes in TC, SRF co-expressed genes, and SRF potential target genes achieved from Cistrome Cancer led to 169 overlapped genes. miR-330-5p was predicted to target SRF, which was further confirmed by dual-luciferase reporter assay. CONCLUSIONS The reduction of SRF appears to play a crucial role in the origin of TC. These properties are accomplished by the target genes of SRF, as a transcription factor, or by the axes with the associated miRNAs.

Background

Thyroid cancer (TC) is one of the most prevalent endocrine malignancies, and was responsible for 567 000 cases and 41 000 deaths worldwide in 2018. The global incidence of TC for both sexes was 3.1%, ranking ninth among all malignant tumors. The global incidence rate in females was 10.2 per 100 000, which is 3 times higher than that in males, and the disease represents 5.1% of the total cancer burden in females and ranking the fifth in the female incidence of cancer. The incidence of TC in recent decades has increased rapidly, especially in China, South Korea, Japan, and other countries in Asia [1-3], which has aroused wide public concern.

Papillary thyroid cancer (PTC) and follicular thyroid cancer (FTC), categorized as differentiated thyroid cancer (DTC), are the major histopathological types of TC [4-7]. There are no evident clinical manifestations in the early stage of TC. Even though most of the cases can be confirmed by ultrasonography, thyroid function test, and fine-needle aspiration cytology via pathology diagnosis, some patients still do not receive an accurate and timely diagnosis [8-11]. Elucidation of the molecular mechanism underlying the occurrence and development of TC will improve understanding of the biological behaviors of TC in clinical practice, which would improve the clinical diagnosis and treatment of TC. However, there are many unclarified molecular events and gene types involved in TC [12-15] and further research is urgently needed.

Among all published potential dysregulated genes involved in TC, there are contradictory reports concerning the serum response factor (SRF), which is a member of the highly conserved MADS (MCM1, Agamous, Deficiens, SRF) box family of transcription factors encoding an immanent nuclear protein. There are only 2 publications concerning the clinical role of SRF in TC. Kim et al. [16] in 2009 was the first to report that the SRF protein expression level showed an obvious increase in 63 cases of PTCs compared to 30 cases of nodular hyperplasia (33%) as detected by immunohistochemistry and Western blot. However, the sample size was small, the study was conducted at a single institution, and the results have not been verified. At the same time, it is aimed at the protein expression level of SRF, the mRNA expression level, and the amplification and mutation status of SRF were not been considered. Thus, many questions remain regarding the clinical role and mechanism of SRF in TC. In contrast, in 2018 Wang et al. [17] reported downregulation of SRF mRNA in 78 cases of TC tissues, as compared to 4 cases of normal thyroid epithelial cells, based on a gene chip dataset (GSE27155). Similarly, the sample size of the study was small and did not exceed 80. The results of the study were based on gene chip data, which are very unlikely to have bias, because there are a large number of high-quality RNA-sequencing and gene chip data available in public databases for wide-ranging analysis. Therefore, the currently available results on SRF in TC are contradictory. The clinical expression of SRF and its specific regulatory mechanisms in TC have not been comprehensively investigated. Therefore, it is crucial to study the expression of SRF and its essential molecular mechanism in TC, to make advances in clinical practice for TC.

In the present study, to comprehensively and objectively evaluate the clinical significance of SRF in TC, we collected all available gene chips and RNA-sequencing datasets, and extracted and integrated SRF expression data. For the first time, 16 RNA-sequencing and gene chip data were summarized. We found that the mRNA expression of SRF showed a tendency to decrease in TC tissues. We then attempted to elucidate the possible molecular mechanisms of this phenomenon. Because SRF is a transcription factor, its potential target genes are key to understanding its biological function. We integrated the putative targets of SRF, the differentially expressed genes in TC, and the co-expressed genes with SRF, and deciphered the hub genes and related pathways of SRF in TC. Finally, we also predicted and preliminarily verified the upstream microRNA of SRF. The above findings will help refine the parameters of further research on the molecular characterization of SRF in TC.

Material and Methods

Collection of RNA-sequencing and gene chip data related to TC expression profile

The Cancer Genome Atlas (TCGA) provides huge comprehensive roadmaps of the essential genomic variations in 33 classes of malignancies, including TC. The UCSC Xena Browser supports functional genomics data from TCGA, including gene-, transcript-, exon-, miRNA-, lncRNA-, and protein-expressions. RNA-sequencing data of patients in TCGA-THCA (Thyroid Cancer) were downloaded through the UCSC Xena Browser. These data included the gene expression data from RNA-sequencing, and various clinical features of 513 TC tissues and 59 TC-adjacent controls, which contains expression value data of SRF. The Genotype-Tissue Expression (GTEx) project has established an integrated public resource with tissue-specific gene expression in non-cancerous tissues from autopsies. Thus, we obtained the expresstion data from normal thyroid tissues from GTEx. All extracted data from TCGA and GTEx were caculated into TPM and presented in log2 (TPM+0.001) using RSEM and Kallisto methods.

All available gene chip data related to SRF in TC tissues and adjacent tissues or normal thyroid tissues were obtained from Gene Expression Omnibus (GEO), ArrayExpress, Sequence Read Archive (SRA), and Oncomine. The following key words were used for data searching: (malignan* OR neoplas* OR cancer OR tumor OR tumor OR carcinoma OR adenoma OR adenocarcinoma) AND (Thyroid OR Thyroidea). The gene chips were screened and selected according to the following criteria: (1) gene chip datasets with information and expression data of SRF; (2) gene chips simultaneously containing TC and adjacent tissues or healthy thyroid tissues; (3) no fewer that 3 samples in each experimental group and control group; (4) gene chips data based on human TC tissues or cancer cell lines rather than other animals; (5) the preoperative patients did not receive radiotherapy and chemotherapy or any other adjuvant treatment.

We concurrently searched the relevant literature in PubMed, EBSCO, Cochrane Central Register of Controlled Trials, Wiley Online Library, Web of Science, Google Scholar, Ovid, EMBASE, and LILACS to August 2019. The key words were: (malignan* OR neoplas* OR cancer OR tumor OR tumor OR carcinoma OR adenoma OR adenocarcinoma) AND (Thyroid OR Thyroidea) AND (SRF OR “serum response factor” OR MCM1 OR “c-fos serum response factor” OR “cFos Serum Response Factor” OR “p67 Serum Response Factor” OR p67-srf). Included studies had to meet the following criteria: SRF expression data was based on human TC tissues, and SRF expression data could be extracted.

Statistical analysis for the SRF level in TC tissues

SRF expression data of all TC tissue and adjacent non-cancerous tissue collected from TCGA, GTEx, and the gene chip database were calculated using the t test using IBM SPSS Statistics 19.0, and are presented as a violin plot using GraphPad Prism v8.0 software. A P value <0.05 was set to be statistically significant. A receiver-operating characteristic (ROC) curve was generated using GraphPad Prism to distinguish TC from non-cancerous thyroid tissues. To enhance the reliability of the data, TC samples whose purity was equal to or greater than 75% were obtained for further verification of differences in expression between TC and non-cancerous tissue. For the clinicopathological role of SRF in TCGA, IBM SPSS 19.0 was used to analyze differences in expression of SRF related to each relevant clinical parameter.

To ensure the reliability of the findings, all the datasets obtained were combined in a comprehensive meta-analyses. STATA 12.0 (StataCorp LP, College Station, TX, USA) was employed for pooling the standardized mean difference (SMD) with a 95% confidential interval (CI). Heterogeneity was assessed by I2 test. When I2 <50%, a fixed-effects model was used, and a random-effects model was used when I2 ≥50%. Begg’s and Egger’s tests were carried out to evaluate whether publication bias existed in the studies. Sensitivity analysis was used to assess the influence of individual chip data on the overall pooling results via omission of each study one at a time. Moreover, we used Meta-DiSc v.1.4 to assess sensitivity, specificity, positive and negative likelihood ratios (LR), diagnostic odds ratio (OR), and summary ROC (sROC) based on the TP, FP, FN, and TN value of all studies.

Caculation of differentially expressed genes (DEGs) in TC using limma, edgR, and Robust rank aggregation (RRA) methods

The DEGs between TC and non-cancer tissues were analyzed by R-based limma method for all gene chip data. To select DEGs from TCGA and GTEx, edgR was used. The cut-off value was set with |log2FC|>1, FDR<0.05. Since there are multiple datasets involved in selecting DEGs in TC, we gathered those upregulated or downregulated genes appearing in at least 2 independent datasets as the first part of DEGs. RRA method was used for the collection for the second part of DEGs, which uses a probabilistic model for aggregation that is robust to noise and calculates the significance probabilities for all the elements in the final ranking [18-23]. The P value cut-off was set as 0.05 for RRA. The second part of DEGs from all gene chip datasets and TCGA/GTEx were obtained by RRA method. Finally, we merged the first and second part for the final candidates for DEGs in TC, which was used in the next step of calculation for the potential target genes of SRF in TC.

Analysis of SRF co-expressed genes in TC using Pearson correlation coefficient calculation

The correlation coefficient indexes between SRF and all other genes in each dataset were calculated by using the correlation coefficient calculation method based on R. The SRF co-expressed genes were selected if |Pearson’s r|≥0.5 and P<0.05. Since this study contained multiple datasets, we selected those co-expressed genes that appeared more than 2 times as candidates for further study.

Collection of SRF target genes in TC using the Cistrome Cancer database

Given that SRF functions as a transcription factor, the target genes of SRF are the direct regulators of the function of SRF. The Cistrome Cancer database compiled by Liu et al. [24] provided predicted target genes of TFs in each cancer type by comprehensively analyzing TCGA molecular expression profile data in cancer and public TFs ChIP-seq data. The predicted target genes of TFs not only showed close correlations with TFs, but were also supported by TF binding information derived from ChIP-seq data. The candidate target genes with both high expression correlation and high regulate potential score in TFs ChIP-seq data were utilized, and the ChIP-seq dataset that met the requirements was selected by using the generalized multiple regression model. After selecting up to 10 relevant ChIP-seq datasets in the model, linear combination of regression coefficients was used to generate the adjusted regulated potential score. The final target genes were those that exceeded the adjusted regulated potential score threshold. We downloaded the target gene list of TC for the follow-up research. Then, the final selected target genes of SRF were those overlapped genes among DEGs, SRF co-expressed genes, and predicted SRF target genes. Next, these overlapped target genes were entered into KEGG pathway and GO enrichment analysis for discovering the underlying molecular mechanisms of SRF in TC. ClueGo in Cytoscape was applied for visualzing the network of GO enrichment results and KEGG pathways. Protein–protein interaction (PPI) enrichment analysis of the SRF target genes was carried out, and a PPI network was constructed using Metascape.

Potential target miRNAs and validation with dual-luciferase reporter assay

miRNAs can perform complex functions and regulatory by forming a feedback or feed-forward loop network, in which the TFs and miRNA regulate the expression of each other, or together upregulate or downregulate a list of target genes directly or indirectly. As a TF, SRF may be regulated or regulate miRNAs to affect the progression of TC. miRWalk, starBase, and TargetScan v.7.2 were used to predict the potential miRNA interacting with the 3′untraslated regions (UTR) of SRF.

We then selected 1 interesting miRNA as an example and performed validating experiment with dual-luciferase reporter assay. Dual-luciferase reporter assay was performed in 96-well plates to validate the miRNA-gene interaction. HEK293T cells were transfected with SRF 3′UTR luciferase reporter plasmids containing wild-type (WT) or mutant (MUT) putative binding sites of the predicted miRNA, together with miRNA mimics or miRNA mimics negative control (GeneChem, Shanghai, China). After 48-h transfection, firefly and renilla luciferase activities were detected with the dual-luciferase reporter assay system (Promega, WI, USA). The relative luciferase activity (Renilla luciferase activity/firefly luciferase activity) was calculated as mean±SD using IBM SPSS Statistics 19.0. P value <0.05 was regarded statistically significant.

Results

Expression and clinical significance of SRF based on RNA-sequencing data

A total of 513 TC cases and 59 non-TC thyroid tissues were obtained from the TCGA database. Since the number of non-cancer controls was much smaller than that of the TC patients, we added 257 samples with thyroid tissues obtained from GTEx to increase the reliability of RNA-sequencing results, and the non-cancerous thyroid controls reached a comparable number of 316 cases. The relative expression value of SRF in the TC group was 3.9320±0.65829, which was clearly lower than that in the non-cancer group (4.6614±0.69757, P<0.0001, Figure 1A). The ROC curve showed that SRF could be used as an indicator to distinguish cancer tissues from non-cancer thyroid tissues (AUC=0.8257, P<0.0001, Figure 1B). Because our findings were inconsistent with a previous report [16], in order to exclude the effect of non-tumor mesenchymal cells on sequencing results, we extracted high-purity samples for in-depth analysis to determine the accurate expression of SRF in TC cells. A total of 196 TC samples with purity ≥75% were obtained, and the expression value of SRF was 3.9263±0.72374, which was also signifantly lower compared to the non-cancer group (Figure 1C) and demostrated a significant downregulation of SRF in TC cells, but not in mesenchymal cells. As is shown in Table 1, the SRF expression level was significantly different in the particular groups according to age, sex, and T stage (Figure 1D–1F).

Figure 1

Downregulation of SRF mRNA level based on RNA-sequencing data. (A) Violin plot for RNA-sequencing data from TCGA and GTEx database; (B) The ROC curve of SRF; (C) Scatter plot for RNA-sequencing data from TCGA database shows a significant down-expression in 196 TC tissues (purity ≥75%) compared with non-cancerous tissues; (D) Age; (E) Gender; (F) T stages.

Table 1

Expression of SRF mRNA in TC based on RNA-sequencing data.

Characteristic		n	Relevant expression of SRF (log2 (TPM+0.001))
			Mean±SD	t/F value	p value
Tissue	TC	512	3.9320±0.65829	−15.137	<0.0001****
	Normal	316	4.6614±0.69757
Sex	Male	136	3.8573±0.60684	−1.637	0.102
	Female	368	3.9656±0.67685
Age	≤45	238	4.0255±0.67526	2.888	0.004**
	>45	266	3.8566±0.63651
Race	Asia	51	3.9990±0.49463	0.473	0.623
	Black/African American	27	4.0761±0.67443
	White	333	3.9557±0.68235
Pathologic T	T1	142	3.9894±0.69378	4.360	0.0048**
	T2	166	4.0292±0.66935
	T3	171	3.8404±0.62304
	T4	23	3.6303±0.42139
Pathologic T	T1–T2	308	4.0109±0.67990	3.268	0.0012**
	T3–T4	194	3.8155±0.60564
Pathologic N	N0	229	3.9463±0.65191	0.744	0.457
	N1	225	3.9016±0.62822
Pathologic M	M0	281	3.9451±0.61552	0.159	0.874
	M1	9	3.9123±0.31240
Pathologic TNM stages	I	283	3.9902±0.69224	1.837	0.139
	II	52	3.8698±0.54982
	III	112	3.8848±0.68050
	IV	55	3.7993±0.48891
Pathologic TNM stages	I–II	335	3.9715±0.67280	1.847	0.065
	III–IV	167	3.8566±0.62373

Altogether, 507 cases in TCGA contained complete pathological clinical information. However, 3 cases (TCGA-ET-A2N1, TCGA-DJ-A13W, TCGA-DJ-A2Q8) did not have SRF expression data. Hence, 504 cases were included in the analysis concerning the clinical role of SRF in TC. Race: Not Evaluated 1 case. Unknown 28 cases. Not Available 64 cases. Pathologic T: TX 2 cases. Pathologic M: MX 213 cases. TCGA-FY-A2QD: no information provided. Pathologic N: NX 50 cases. Pathologic stages: TCGA-FY-A2QD and TCGA-EL-A3CP: no information provided.

Consolidation of SRF expression in data from gene chips

Fifteen GSE datasets were retrieved from GEO, SRA, ArrayExpress and Oncomine databases: GSE3467, GSE3678, GSE6004, GSE6339, GSE9115, GSE27155, GSE29265, GSE29315, GSE33630, GSE53157, GSE50901, GSE35570, GSE58545, GSE60542, and GSE65144 (Table 2). The expression levels of SRF in TC and non-cancer tissues in each of the included datasets are presented in Figure 2. Among these datasets, GSE3678, GSE27155, GSE29265, GSE33630, GSE50901, GSE35570, GSE58545, GSE60542, and GSE65144 showed a significant downregulated expression pattern of SRF in TC compared to that in non-cancerous tissues (P<0.05). Figure 3 shows the ROC curves of SRF in each GSE dataset. Among them, GSE3678, GSE27155, GSE29265, GSE29315, GSE33630, GSE35570, GSE58545, GSE60542, and GSE65144 illustrated moderate-to-high distinguishing capacity of SRF level to separate TC from non-cancerous tissues (P<0.05). The forest plot in Figure 4A shows the SRF expression data from the 15 gene chips with the pooled SMD of SRF being −1.00 (95% CI: −1.30 to −0.71) by the random-effects model. The I-squared value was 71.0% and the p value was <0.001. The SRF expression level was remarkably decreased in 606 cases of TC tissues. We also calculated the publication bias using a funnel plot (Figure 4B). The P value from Begg’s test was 1.000 (>0.05) and that from Egger’s test was 0.679 (>0.05), which showed no significance of publication bias. Sensitivity analysis (Figure 4C) showed no significant difference among all the included datasets. In the combined analysis in the random-effects model, the sensitivity was 0.79 (95% CI=0.76–0.83; I2=80.7%) (Figure 5A), specificity 0.65 (95% CI=0.61–0.70; I2=80.1%) (Figure 5B), the positive likelihood ratio was 2.59 (95% CI=1.81–3.70; I2=73.0%) (Figure 6A), the negative likelihood ratio 0.33 (95% CI=0.23–0.48; I2=76.7%) (Figure 6B), and the diagnostic OR was 9.37 (95% CI=4.30–18.27; I2=69.7%) (Figure 7A). Most importantly, the area under the sROC curve of SRF in TC was 0.8231 (Figure 7B). These results indicate that SRF has a potential role in distinguishing TC from non-cancerous tissues based on a TC population of 606 cases.

Table 2

The basic information of included RNA-sequencing and gene chip data.

ID	Authors	Year	Country	Citation	N1 (Cancer group)	M1	SD1	N2 (Normal control)	M2	SD2
GSE3467	Sandya L et al.	2005	USA	PMID: 16365291 [45]	9	7.2327	0.60427	9	7.6055	0.50707
GSE3678	Ismael R et al.	2006	USA	/	7	7.201	0.22267	7	8.355	0.46224
GSE6004	Sandya L et al.	2006	USA	PMID: 17296934 [46]	14	6.083	0.30475	4	6.49	0.66481
GSE6339	Fontaine J et al.	2007	France	PMID: 17968324 [47]	48	1.023	0.53493	134	1.2965	1.12104
GSE9115	Salvatore G et al.	2007	USA	PMID: 17981789 [48]	15	0.513	0.53191	4	0.4713	0.38562
GSE27155	Rork K et al.	2011	USA	PMID: 16609007 [49]	156	2.595	0.20521	42	2.7301	0.3127
GSE29265	Gil T et al.	2012	Belgium	/	29	5.088	0.34102	20	5.6747	0.45155
GSE29315	Gil T et al.	2012	Belgium	/	31	7.8073	0.48499	40	8.0665	0.58535
GSE33630	Gil T et al.	2012	Belgium	PMID: 22266856 [50]PMID: 22828612 [51]	60	6.6565	0.21228	45	7.1179	0.49957
GSE53157	Branca M et al.	2013	Portugal	PMID: 19809427 [52]	24	7.5897	0.22492	3	7.7751	0.02563
GSE50901	Barros-Filho MC et al.	2014	Brazil	PMID: 25867809 [53]	51	−1.4673	0.52346	4	−0.0917	1.47755
GSE35570	Swierniak M et al.	2015	Poland	PMID: 26810418 [54]	65	4.4749	0.23213	51	5.106	0.72787
GSE58545	Swierniak M et al.	2015	Poland	PMID: 26625260 [55]	27	4.1394	0.34757	18	4.6134	0.48307
GSE60542	Maxime T et al.	2015	Belgium	PMID: 25965298 [56]	58	5.7927	0.14496	34	6.2147	0.51585
GSE65144	John A et al.	2015	USA	PMID: 25675381 [57]	12	7.8081	0.67202	13	8.5873	0.84682
TCGA and GTEX	/	2019	/	/	512	3.932	0.65829	316	4.6614	0.69757

Figure 2

Scatter plots of SRF expression in the included gene chip data. (A) Scatter plots of SRF expression in GSE3467; (B) GSE3678; (C) GSE6004; (D) GSE6339; (E) GSE9115; (F) GSE27155; (G) GSE29265; (H) GSE29315; (I) GSE33630; (J) GSE53157; (K) GSE50901; (L) GSE35570; (M) GSE58545; (N) GSE60542; (O) GSE65144.

Figure 3

ROC curves of SRF in TC tissues from all the included gene chips. (A) ROC curves of SRF in GSE3467; (B) GSE3678; (C) GSE6004; (D) GSE6339; (E) GSE9115; (F) GSE27155; (G) GSE29265; (H) GSE29315; (I) GSE33630; (J) GSE53157; (K) GSE50901; (L) GSE35570; (M) GSE58545; (N) GSE60542; (O) GSE65144

Figure 4

The comprehensive SRF expression level in TC generated by 15 gene chips. (A) Forest plot of SRF expression in TC based on included gene chips. TC vs. normal, fixed-effects model. (B) Funnel plot of different datasets related to SRF. (C) Sensitivity analysis of included gene chips.

Figure 5

Sensitivity and Specificity of SRF in TC tissues in relevant gene chips. (A) Sensitivity; (B) Specificity

Figure 6

Positive LR and Negative LR of SRF in TC tissues in relevant gene chips. (A) Positive LR; (B) Negative LR

Figure 7

Diagnostic OR and SROC Curve of SRF in TC tissues in relevant gene chips. (A) Diagnostic OR; (B) SROC Curve

Overall SRF expression level in all gene chip and RNA-sequencing data

The expression data of SRF was unable to be extrated from the literature, so only aquired gene chip and RNA-sequencing data was entered into analysis. The forest plot (Figure 8A) included the SRF expression data from TCGA database and gene chips. The random-effects model was used in the analysis due to the I-squared being 71.2%. The I-squared value may be produced by patient variation, samples processing methods, and the statistical model used. The pooled SMD of SRF was −1.00 (95% CI −1.24 to −0.75) using the random-effects model, indicating dramaticly decreased expression of SRF in TC tissues based on 1118 cases. The funnel plot in Figure 8B reveals an absence of publication bias among these studies, since the P value gained from Begg’s test was 0.685 (>0.05) and that from Egger’s test was 0.763 (>0.05). The sensitivity analysis (Figure 8C) showed no significant difference among studies. The pooled sensitivity and specificity were 0.77 (95% CI: 0.74–0.79) and 0.68 (95% CI: 0.63–0.72), respectively (Figure 9A, 9B), positive likelihood ratio was 2.77 (95% CI=1.92–4.00; I2=77.2%) (Figure 10A), the negative likelihood ratio was 0.34 (95% CI=0.25–0.46; I2=76.3%) (Figure 10B), the diagnostic OR was 9.58 (95% CI=5.19–17.66; I2= 70.0%) (Figure 11A), and the area under the curve of sROC was 0.8251 (Figure 11B). These results further confirm that SRF can differentiate TC from non-cancerous tissues, as evidenced by the large sample number of 1118 cases. We also attempted to analyze the prognostic value of SRF in TC. Unfortunately, there was insufficient data a systematic analysis to calculate the pooled hazard ratio. Based on the RNA-sequencing data, there was no evidence supporting that SRF has a predictive value for patient survival status of TC (data not shown). Larger cohorts are required to investigate the prognostic effect of SRF in TC.

Figure 8

Overall expression level of SRF in TC with all data from gene chips and RNA-sequencing. (A) Forest plot of SRF expression in TC based on included gene chips and TCGA/GTEx RNA-sequencing data. TC vs. non-cancerous, random-effects model. (B) The funnel plot showing the publication bias of gene chips and TCGA/GTEx data. (C) Sensitivity analysis of gene chips and TCGA/GTEx data.

Figure 9

Sensitivity and Specificity of SRF in TC tissues by data from relevant gene chips and TCGA/GTEx RNA-sequencing. (A) Sensitivity; (B) Specificity

Figure 10

Positive LR and Negative LR of SRF in TC tissues by data from relevant gene chips and TCGA/GTEx RNA-sequencing. (A) Positive LR; (B) Negative LR

Figure 11

Diagnostic OR and SROC Curve of SRF in TC tissues by data from relevant gene chips and TCGA/GTEx RNA-sequencing. (A) Diagnostic OR; (B) SROC Curve

GO and KEGG pathway analysis for selected target genes from DEGs, SRF co-expressed genes and predicted SRF target genes

After analyzing the DEGs in 15 gene chips and 1 RNA-sequencing data, a total of 9243 DEGs were obtained, and the volcano plots for each study are presented in Figure 12. RRA method was used to merge all the DEGs from each study, after which 107 genes were significantly upregulated and 105 genes were significantly downregulated. The top 20 significantly upregulated and downregulated genes are represented in Figure 13. Finally, we selected the union set of 2248 DEGs that appeared more than 2 times in 16 datasets and the 212 DEGs from RRA as final DEGs. Interestingly, RRA did not increase the gene number; in fact, all 2248 DEGs were determined for the following step. In terms of SRF co-expressed genes, 2637 genes appeared in more than 2 studies among the 16 datasets. As for predicted SRF target genes, 3766 genes were obtained from Cistrome Cancer. The intersection of above DEGs, SRF co-expressed genes, and SRF potential target genes led to 169 overlapped genes (Figure 14A), which were entered into GO enrichment and KEGG pathway analysis. The GO enrichment included 3 parts: biological processes (BP), cellular component (CC), and molecular function (MF). For BP, genes were mainly enriched in negative regulation of ossification, actin filament bundle assembly, response to corticosteroid, skeletal muscle organ development, metanephros development, positive regulation of epithelial cell proliferation, and respiratory system development (Figure 14B). For MF, the target genes were involved in promoter-specific chromatin binding, activating transcription factor binding, heparin binding, peptidase activator activity, insulin-like growth factor binding, fibronectin binding, protein kinase C binding, and phosphatidylserine binding (Figure 14C). In CC, genes were highly enriched in platelet dense granule, A band, desmosome, and intercalated disc (Figure 14D). In KEGG pathway analysis, AGE-RAGE signaling pathway in diabetic complications, breast cancer, heparin binding, thyroid hormone synthesis, thyroid hormone signaling pathway, protein kinase C binding, TNF signaling pathway, Apelin signaling pathway, fibronectin binding, and activating transcription factor binding were significant. Thyroid cancer is one of the pathways in which SRF target genes are mainly enriched, which confirms the crucial role of SRF in TC (Figure 14E, Table 3).

Figure 12

(A–P) Volcano plots of each dataset in the current study. Red dots represent upregulated expression genes, and green dots represent downregulated ones. Grey dots represent stable ones.

Figure 13

Dysregulated genes in TC tissue by RRA method. Top 20 genes with red backgrounds were significantly upregulated, while the other top 20 genes with green were significantly downregulated expressed in TC based on 16 included datasets.

Figure 14

GO and KEGG analyses of 169 overlapped target genes using ClueGo in Cytoscape. (A) Venn diagram of the overlapped target genes of SRF. (B) Biological process (BP); (C) Molecular function (MF); (D) Cellular component (CC); (E) KEGG pathway annotations.

Table 3

The GO annotation and KEGG pathway analysis of the potential targets of SRF in TC.

Ontology	GO term	P value corrected with Bonferroni step down	Count	Involved genes
BP (P<0.01)	Metanephros development	0.01	7	BASP1, EGR1, FAT4, FOXJ1, ID3, LIF, OSR1
	Respiratory system development	0.00	10	ALDH1A3, BASP1, CRISPLD2, CTGF, FGF7, FOXJ1, ID1, LIF, SRF, WNT11
	Positive regulation of epithelial cell proliferation	0.01	10	CCND1, CDH3, FGF7, ID1, IQGAP3, JUN, NR4A1, NR4A3, OSR1, PTK2B
	Actin filament bundle assembly	0.01	9	CENPJ, CTGF, FAM107A, FHOD1, ID1, LPAR1, PTK2B, SRF, WNT11
	Actin filament bundle organization	0.01	9	CENPJ, CTGF, FAM107A, FHOD1, ID1, LPAR1, PTK2B, SRF, WNT11
	Response to corticosteroid	0.00	11	CCND1, CTGF, DUSP1, FAM107A, FOS, FOSB, FOXO1, GHR, PPARGC1A, S100B, SERPINF1
	Response to glucocorticoid	0.00	10	CCND1, DUSP1, FAM107A, FOS, FOSB, FOXO1, GHR, PPARGC1A, S100B, SERPINF1
	Skeletal muscle organ development	0.00	10	BASP1, EGR1, EGR2, EPHB1, FOS, GPX1, HDAC4, RBM24, RCAN1, S100B
	Skeletal muscle tissue development	0.01	9	EGR1, EGR2, EPHB1, FOS, GPX1, HDAC4, RBM24, RCAN1, S100B
	Osteoblast differentiation	0.00	11	CYR61, FHL2, GDF10, HDAC4, ID1, ID3, PTK2B, SNAI1, TOB1, TP53INP2, WNT11
	Regulation of ossification	0.00	10	CYR61, EGR2, GDF10, HDAC4, ID1, ID3, IER3, OSR1, PTK2B, TOB1
	Negative regulation of ossification	0.00	7	GDF10, HDAC4, ID1, ID3, IER3, PTK2B, TOB1
CC (P<0.05)	Intercalated disc	0.02	3	FHOD1, FXYD1, SCN4B
	Desmosome	0.01	3	EVPL, RCAN1, SRPX
	A band	0.02	3	CRYAB, FHL2, HDAC4
	Platelet dense granule	0.00	3	CLEC3B, ITPR1, TIMP3
MF (P<0.05)	Phosphatidylserine binding	0.04	3	CAVIN2, ICAM5, THBS1
	Fibronectin binding	0.02	3	CTGF, FBLN1, THBS1
	Protein kinase C binding	0.03	3	CAVIN2, DACT1, PTK2B
	Insulin-like growth factor binding	0.03	3	CTGF, CYR61, IGFBPL1
	Heparin binding	0.00	9	CLEC3B, COL13A1, CRISPLD2, CTGF, CYR61, FGF7, KLK10, PCOLCE2, THBS1
	Peptidase activator activity	0.04	3	FBLN1, PCOLCE2, RPS27L
	Promoter-specific chromatin binding	0.04	3	EGR1, HDAC4, PPARGC1A
	Activating transcription factor binding	0.01	5	EGR2, FOS, HDAC4, JUN, WFS1
	RNA polymerase II activating transcription factor binding	0.04	3	EGR2, FOS, JUN
KEGG (P<0.05)	Longevity regulating pathway	0.14	3	CRYAB, FOXA2, FOXO1
	Apelin signaling pathway	0.04	6	CCND1, CTGF, EGR1, HDAC4, ITPR1, PPARGC1A
	GnRH signaling pathway	0.12	4	EGR1, ITPR1, JUN, PTK2B
	Thyroid hormone synthesis	0.04	3	GPX1, GPX3, ITPR1
	Thyroid hormone signaling pathway	0.06	5	CCND1, FOXO1, PLN, RCAN1, RCAN2
	AGE-RAGE signaling pathway in diabetic complications	0.01	6	CCND1, EGR1, FOXO1, JUN, SELE, THBD
	TNF signaling pathway	0.07	5	FOS, JUN, LIF, MAP3K8, SELE
	Amphetamine addiction	0.11	3	FOS, FOSB, JUN
	Colorectal cancer	0.03	5	CCND1, FOS, GADD45B, JUN, TCF7L1
	p53 signaling pathway	0.06	4	BID, CCND1, GADD45B, THBS1
	Endometrial cancer	0.14	3	CCND1, GADD45B, TCF7L1
	Thyroid cancer	0.07	3	CCND1, GADD45B, TCF7L1
	Basal cell carcinoma	0.12	3	GADD45B, TCF7L1, WNT11
	Melanoma	0.08	3	CCND1, FGF7, GADD45B
	Breast cancer	0.01	7	CCND1, FGF7, FOS, GADD45B, JUN, TCF7L1, WNT11

GO – gene ontology; KEGG – Kyoto Encyclopedia of Genes and Genomes; BP – biological process, CC – cellular component; MF – molecular function.

Construction of the PPI network

A PPI network (Figure 15A) was constructed and the MCODE algorithm was then leveraged to the network to determine neighborhoods where proteins are densely connected (Figure 15B). GO enrichment analysis was used in the MCODE network to assign “meanings” to the network components. MCODE components include proteins SRF, FHL1, and FHL2. Interestingly, both FHL1 and FHL2 own binding sites with SRF, as evidenced by public ChIP-seq sequencing data (Figure 15C, 15D). GO data showed that these proteins are associated with positive regulation of cell death, muscle structure development, and positive regulation of apoptotic process.

Figure 15

Protein–protein interaction network of 169 overlapped target genes of SRF in TC. (A) Center genes from the protein–protein interaction network. Nodes represent gene-encoded proteins. Connections between nodes show the regulatory association between proteins. (B) PPI MCODE components, including proteins SRF, FHL1 and FHL2. (C) The binding site of FHL1 with SRF: GEO or ENCODE: GSM1505773, Cell Line: HUES64, binding score: 0.602, Coordinate: chrX: 136147399–136211359, visualized by UCSC. (D) The binding sites of FHL2 with SRF: GEO or ENCODE: GSM803425, Cell Line: H1, binding score: 0.657, Coordinate: chrX: 105360825–105399118, visualized by UCSC.

SRF may be a direct target of miR-330-5p

Seven (miRWalk, Microt4, miRMap, PITA, RNA22, RNAhybrid, and TargetScan) of the 12 target prediction algorithms of miRWalk indicated that SRF is one of the target genes of miR-330-5p. Predicted consequential pairing of the target region on SRF and miR-330-5p was identified in TargetScan (Figure 16A). In starBase v.3.0, we next discovered a relationship between miR-330-5p and SRF expression, but it was not significant (Figure 16B). Dual-luciferase reporter assay was also leveraged to confirm whether SRF was a direct biological target of miR-330-5p. In the transfection of SRF-3′UTR-WT luciferase reporter plasmids, the luciferase activity of the miR-330-5p group was lower than that of the miR-330-5p NC group (P<0.0001), indicating a significant difference. The luciferase activity of SRF-3′UTR-MUT showed no significant alteration between the 2 groups (P>0.05), which confirmed that miR-330-5p particularly binds to the 3′UTR of SRF (Figure 16C).

Figure 16

Potential target miRNA of SRF in TC. (A) Predicted consequential pairing of SRF 3′UTR and miR-330-5p; (B) Correlation analysis of hsa-miR-330-5p and SRF; (C) Correlation between miR-330-5p with SRF with a dual-luciferase reporter assay. Relative luciferase activity in cells co-transfected with SRF-3′UTR-WT or SRF-3′UTR-MUT and miRNA negative control or miR-330-5p mimic (**** P<0.0001).

Discussion

In this research, we integrated available public RNA-sequencing and gene chip data of TC tissues to comprehensively elucidate the clinical implications of SRF. We found clear downregulation of SRF mRNA levels in 16 independent datasets with 1118 TC patients. As a transcription factor, SRF fulfills its biological function via targeting a certain group of genes. Combining the DEGs of TC tissues, SRF co-expressed genes, and SRF predicted target genes; 169 genes had a high probability of being the real target genes of SRF in TC. Moreover, SRF also acts as a target of miR-330-5p, as verified by dual-luciferase reporter assay. These findings help to demonstrate the critical role of SRF in the regulatory mechanism in TC.

Overexpression of SRF has been noted in several cancer types, including gastric cancer [25], cervical cancer [26], prostate cancer [27], and hepatocellular carcinoma [28]. However, there are only 2 publications concerning the clinical role of SRF expression level in TC [16,17], which have reported an opposite expression pattern of SRF in TC. The first study on SRF expression in TC was published in 2009 by Kim et al. [16] and examined the expression of SRF protein in 63 cases of PTCs, 9 cases of anaplastic TCs, 30 cases of follicular adenoma, and 30 cases of adenomatous hyperplasia. The results showed that 50 of 63 PTCs were SRF-positive, accounting for 79% of the total PTCs; while 6 of 9 anaplastic TCs were SRF-positive, accounting for 67%. For the non-cancerous controls, 18 of 30 follicular adenomas were SRF-positive (60%), and 10 of 30 nodular hyperplasia were SRF-positive (33%). Overall, the positive expression of SRF protein in TCs tended to increase, but the degree of increase was not obvious (79% for PTCs vs. 60% for adenomas). However, the sample size of this study was small, which could lead to data deviation. This is the only study so far assessing the clinical significance of SRF in TC by immunohistochemistry and Western blot, and the results have never been verified by other research groups.

Interestingly, Wang et al. [17] reported the opposite trend of SRF expression level in TC in 2018. They analyzed data from a gene chip, including 17 cases of thyroid adenoma, 78 cases of TCs, and 4 cases of normal thyroid epithelia. They found that SRF mRNA levels were clearly downregulated in TCs, but the sample size of this study was also small, and the results were probably biased. To obtain comprehensive SRF expression data, we carried out several levels of research in the present study. Firstly, we performed a recalculation of the SRF mRNA expression from the RNA-sequencing data. To balance the numbers of the experimental and control groups, and to reduce the potential statistical deviation, we added RNA-seq data from the GTEx autopsy database of non-cancer thyroid tissues. The results showed that SRF mRNA expression in 513 TCs was significantly lower than that in non-cancer thyroid tissue. However, this result is only from the data of a single cohort. To determine SRF gene expression in a larger sample size, we then collected various data from multiple databases. From GEO, ArrayExpress, SRA, Oncomine, and various literature databases, we obtained data from another 15 gene chips with SRF gene expression. We conducted a statistical analysis for each dataset and found that in most of the gene chip data, SRF mRNA was expressed at low levels, consistent with the RNA-sequencing data, but some were not consistent. To gain a relatively comprehensive representation of SRF mRNA expression in these gene chips, we performed SMD and sROC calculations. The summarized results based on 606 cases showed that the SMD was −1.00 (95% CI: −1.30, −0.71), and the AUC was 0.8231. Finally, we also merged all of the available data and found that the overall SMD was −1.00 (95% CI: −1.24, −0.75) and the AUC was 0.8251, which was identical to that found by gene chips alone. The above results suggest that SRF expression is decreased or absent in TC. Because the above gene chips and RNA-sequencing experimental materials are from human tissues that contain a variety of cell types, in order to exclude the influence of non-tumor cells on SRF expression values, we separately extracted the cases with tumor purity information for analysis, and found that the expression of SRF in TCs was indeed lower, and there was no contrary situation of higher SRF expression. In summary, we can conclude that SRF is obviously downregulated in TCs in a total of 1118 cases, which unsurprisingly coincides with the previous study [17] included in our current integration.

However, another study on SRF in TC did not find lower SRF expression [16]. The reasons for the different results of this study [16] compared with Wang et al. [17] and our current work may be due to the following reasons. Firstly, the sample size in the Kim et al. study was relatively small, having no more than 100 cases, and the number of non-cancerous controls was even smaller, which is insufficient to represent the exact overall situation of TCs. In contrast, our study is based on 1118 cases and the sample size is over 16 times more that of Kim et al., which lends more credibility to our results. Secondly, Kim et al. measured the protein levels, while Wang et al. [17] and the present study examined the mRNA levels of SRF.

Inconsistent expression abundance of mRNA and protein levels of a gene is also possible. The following are possible reasons for this. First, if the protein level of a gene is elevated in tissues, it likely reduces its transcription level through a mechanism of negative feedback regulation; on the contrary, if the protein level is reduced, the cells may promote its transcriptional level to keep the balance. Secondly, the time and space at which transcription and translation of eukaryotic gene expression occur have spatial and temporal intervals. After transcription, there are several stages of post-transcriptional processing, degradation of transcription products, translation, post-translational processing, and modification. Therefore, it is understandable that the level of transcription and the level of translation are not completely consistent. Due to the different time points of detection, the mRNA may have degraded when the protein reaches its peak, or it may occur when the amount of protein is still increasing, but the mRNA reaches its peak. Thirdly, it is of course possible that the process of SRF translation may also be regulated by other factors, such as its target genes, which can also modulate the expression of SRF, or post-translational modifications may also result in different expression of mRNA and protein levels. In conclusion, the protein level and clinical significance of SRF in TCs, as well as their specific mechanisms, need further exploration.

In the 2 previous studies of SRF in TCs [16,17], the molecular mechanism of SRF action has not been investigated. Therefore, we also explored the reason and prospective mechanism for the decrease of SRF mRNA level in TCs. Since SRF is a transcription factor, its target genes play a crucial role in the biological function of SRF. To gather the potential target genes of SRF in TC, we collected 3 parts of the gene to narrow the range. Firstly, the differentially expressed gene profiles of TC were sorted out. We synthetically analyzed all gene expression data from 16 datasets and obtained 2248 differentially expressed genes, which may play a role in the incidence and progression of TC. Secondly, we summarized the co-expressed genes of SRF in TC. We calculated the correlation coefficients of SRF with other genes based on the information in the same 16 datasets mentioned above, because the genes associated with SRF are better target gene candidates. Altogether, 2637 co-expressed genes of SRF were taken into the next step. Thirdly, with the assistance of the Cistrome Cancer database [24], which provides the predicted targets of a TF, we achieved a set of potential target genes (n=3766) of SRF based on ChIP-sequencing experiments. An integration and modeling of more than 10 000 cancer molecular profiles from TCGA, as well as more than 23 000 ChIP-seq and chromatin accessibility profiles, were systematically assessed by Cistrome Cancer. The reconstruction of functional enhancer profiling, “super-enhancer” targets, and predictions of TF and their targets are available in Cistrome Cancer. Hence, the final 169 genes intersected by the above 3 parts became the potential target genes of SRF in TC. Among these potential targets, FHL1 and FHL2 were selected based on the PPI network to show the binding sites of SRF. Furthermore, the relative pathways are essential for the elucidation of the underlying mechanisms of SRF in TC. Several thyroid-related pathways were enriched by KEGG analysis, including thyroid hormone synthesis, thyroid hormone signaling pathway, and thyroid cancer. Interestingly, the pathway of activating transcription factor binding was also refined, which indicates that SRF plays a pivotal role in the transcription modulation in TC.

Transcription factors may play a role by targeting their downstream DNA targets. At the same time, they themselves may also be targets of some non-coding RNAs. In many subtypes of non-coding RNAs, miRNAs are short non-coding RNAs with a length of about 22 nucleotides, involved in post-transcriptional modulation of gene expression. SRF has been found to exert its biological functions in different diseases by acting as a target gene for certain miRNAs. For instance, in some non-cancerous diseases, an SRF/miR-1 axis has been documented in heart failure [29]. SRF acts as a target of miR-22 in human umbilical vein endothelial cells (HUVECs) [30]. miR-483-3p regulates endothelial progenitor cells dysfunction in deep vein thrombosis patients via SRF [31]. The cardiac fibrosis is also modulated by an SRF/miR-133a axis [32]; for example, the target miRNAs miR-647 [33], miR-101-3p [34], and miR-199a-5p [35] have all been confirmed to target SRF in human gastric cancer. Furthermore, the SRF-miR-29b axis obstructs the infiltration and metastasis of non-small cell lung cancer [36]. In TC, the relationship between SRF and miRNAs has not been reported. Through target gene prediction and differentially expressed gene profile, we predicted the potential upstream miRNAs of SRF in TC. An unpublished miR-330-5p was selected to carry out preliminary dual-luciferase reporter assay, which confirmed that miR-330-5p can bind the 3′UTR of SRF. The expression level and mechanism of miR-330-5p have been studied in various cancers, including melanoma [37] hepatocellular carcinoma [38,39], pancreatic cancer [40,41], cervical cancer [42], prostate cancer [43], and non-small cell lung cancer [44]. However, the expression and function of miR-330-5p in TC and its molecular mechanism have not been clarified, and these topics warrant further study.

Some of the work that was not been done in the present study can be considered in the future, and these are also the limitations of our study. Non-invasive detection may be of clinical value to discover the potential biomarkers used to diagnose the disease or to assess the disease progression and the prognosis, as well as the therapeutic effect. Hence, the expression of SRF in the body fluids of TC patients may be detected in the future. Moreover, the different functions and roles of SRF, especially the prognostic value of SRF in each subtype of TC, remain to be explored. Lastly, the cellular biological function of SRF and the potential molecular mechanism, including the downstream target genes and the upstream miRNAs, also need to be elucidated by in vitro and in vivo experiments.

Conclusions

For the first time, we have combined the RNA-sequencing and gene chip data of TCs. Based on the expression data of 1118 cases of TC patients, we confirmed that the SRF mRNA level was clearly downregulated in TC tissues compared to non-cancerous thyroid controls. The reduction or absence of SRF mRNA may play a crucial role in the origin of TC. These functions may be accomplished by the target genes of SRF, as a transcription factor, or by the axes with the associated miRNAs. Some target genes and potential miRNAs selected in this study provide preliminary suggestions for further study of the molecular mechanism of SRF in TC, but more experiments need to be carried out to clarify these hypotheses.