Sage Journals: Discover world-class research

Abstract

BACKGROUND:

The incidence of thyroid cancer has risen rapidly over the last decades. Although mortality rates are relatively low compared to other cancers, the rate of new cases started to increase in the early 2000s. While tumor suppressors and oncogenes were recently identified in thyroid cancer, the potential roles of these genes in thyroid cancer remain unclear.

OBJECTIVE:

Analyze the roles and functions of tumor suppressors and oncogenes in thyroid cancer.

METHODS:

Thyroid cancer data were collected from public databases, such as the UCSC Xena database of TCGA thyroid cancer, TISIDB, and UALCAN. The genes frequently associated with unfavorable thyroid cancer were examined and validated. The association of these target genes with thyroid tumorigenesis, stages, subtypes, and survival rates were analyzed. Additionally, the genes aberrantly expressed in thyroid cancer and significantly involved in thyroid tumorigenesis, stages, subtypes, and survival rates were identified.

RESULTS:

Female sex was identified as a risk factor for thyroid cancer. The expression of PAPSS2, PDLIM3, COPZ2, ALDH1B1, ANTXR1, GUF1, and SENP6 negatively correlated with thyroid cancer prognosis.

CONCLUSION:

Female sex was a risk factor for thyroid cancer. In addition, our analysis suggested that PAPSS2, PDLIM3, COPZ2, ALDH1B1, ANTXR1, GUF1, and SENP6 are negatively correlated with the prognosis of thyroid cancer. The expression of ANTXR1, GUF1, and PDLIM3 was weakly associated with thyroid cancer’s immune and molecular subtypes.

Keywords

Database thyroid cancer sex differences prognosis gene anticancer drugs

1. Introduction

The incidence of thyroid cancer is rapidly increasing in the United States. According to the National Cancer Institute (NCI) report, there were around 915,664 individuals with thyroid cancer in the United States in 2019, and an

Table 1
The 5- and 10-year trends of thyroid cancer during the last decade

Sex	Average annual percent change (AAPC) estimates
	Year range	AAPC (%)	Lower 95% CI.	Upper 95% CI.	$P$ -Value	Direction
Female	2010–2019	$-$ 0.4	$-$ 1.3	0.4	0.31	Not Significant
	2015–2019	$-$ 2.4	$-$ 3.5	$-$ 1.3	$<$ 0.01	Falling
Male	2010–2019	0.5	$-$ 0.5	1.5	0.30	Not Significant
	2015–2019	$-$ 1.0	$-$ 2.3	0.2	0.10	Not Significant

Data source: SEER 22 2015–2019. By Sex, Delay-adjusted SEER Incidence Rate, All Races, All Ages, All Stages.

estimated 43,000 new cases will occur in 2022. Thyroid cancer is an abnormal growth of thyroid cells, forming a lump within the thyroid. Most of the thyroid nodules are non-cancerous (benign), and only 5% are cancerous. Based on the histological examinations, four types of thyroid cancer were identified: papillary, follicular, medullary, and anaplastic. As the most common type of thyroid cancer, papillary thyroid cancer (PTC) accounts for 90% of thyroid cancers and has a better prognosis than the other types [1]. Genes associated with thyroid cancer have been identified, including BRAF, RET/PTC, and RAS [2, 3, 4]. However, no significant association between these genes and PTC subtypes and sex was observed [5].

Thyroid cancer is one of the few cancers that are more common in females, with a 2.9-fold higher risk of developing thyroid follicular cancer in women than in men (female: male ratio 16.3:5.7) [6]. Based on SEER’s 2019 cancer prevalence proportions in 12 Areas and the 2019 US population estimates from the US Bureau of the Census, 85714 females (60%) and 31,673 males (30%) were living with cancer. The estimates included white individuals, including Hispanics, 65–74 years old, and less than 27 years since the initial diagnosis. Gender was associated with thyroid cancer incidence (Table 1), as well as aggressiveness and prognosis [7, 8]. Recent studies have reported that both chromosomal and gonadal factors contribute to sex bias in disease traits such as obesity, nonalcoholic fatty liver disease, cardiovascular diseases, diabetes, and COVID-19 infection [9, 10, 11, 12, 13]. It seems that sex chromosomes play a relatively more influential role in these age-dependent diseases because older individuals undergo menopause and reduced hormone production [14, 15]. Most sex chromosome effects result from X-chromosome inactivation (XCI) escape, whereby some genes escape from XCI during aging. XCI is the inactivation of one of the X-chromosomes in females to balance X dosage in males [18]. Some genes on the X chromosome can escape XCI under specific conditions, resulting in imbalanced gene expression between males and females [16, 17]. It has been proposed that the sex hormones and the X chromosome could be involved in the sex-associated thyroid cancer risk [19, 20]. However, the mechanisms underlying the male preponderance are not entirely known yet.

Recently, bioinformatics analysis using pan-cancer next-generation sequencing (NGS) database has been used for cancer research, which could provide more information about molecular features than conventional histological classification. To better understand the gene functions within the thyroid tumor and improve survival, we screened and explored the function of tumor suppressors and oncogenes, then used these genes to identify potential thyroid cancer predictor genes. We aimed to identify gene expression aberrations and interpret their significance on the growth and development of the thyroid tumor.

In this study, we performed a comprehensive bioinformatics analysis of thyroid cancer-associated genes to identify potential tumor suppressors and oncogenes associated with serous thyroid cancer. In addition, the expression data analysis provided specific information on the molecular basis and subtypes of thyroid cancer and elucidated the genes driving thyroid tumors.

2. Materials and methods

2.1 Ethical statement

The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki (Seventh revision (2013)). The study was approved by the Ethics Committee of the Children’s Hospital of Fudan University, and informed consent was obtained from all patients.

2.2 Data source and pretreatment

The cancer stat facts of the National Institute of Health (National Cancer Institute) were used to obtain the thyroid cancer case and death numbers in the US, including the rates of new cases and deaths, 5-year relative survival, new cases by age, race/ethnicity, and sex.

The RNA sequencing (RNA-seq) and survival data were obtained from the UCSC Xena database of TCGA Thyroid Cancer (THYROID CANCER) (https://xenabrowser.net/) [8]. The Transcripts per million (TPM) was calculated, followed by log2 transformation. The downloaded expression data were analyzed and processed using R (Version 4.2.1 2022-06-23 ucrt).

TISIDB (http://cis.hku.hk/TISIDB/index.php) is an online website focusing on the interaction between cancer and the immune system and containing multiple heterogeneous data types [9]. In this study, TISIDB was employed to construct the expression heatmap of prognostic genes in different immune subtypes and identify correlations with cytokines.

The UALCAN (http://ualcan.path.uab.edu/index.html) database is an interactive web portal for in-depth analyses of TCGA gene expression data [10]. The expression of prognostic genes in thyroid cancer stages and types was obtained from this database, and the survival data of prognostic genes were also validated.

2.3 Distribution analysis of thyroid cancer

We searched for thyroid cancer in the NCI’s Annual Report to the Nation on the Status of Cancer and obtained the Cancer Stat Facts: Thyroid Cancer (https://seer.cancer.gov/statfacts/html/thyro.html). The clinical characteristics (e.g., sex, age, tumor stage, and survival time) of thyroid cancer were downloaded accordingly, then loaded to GraphPad Prism (Version 9.4.1.681).

2.4 Association between oncogenes expression and thyroid cancer

The relationship of candidate oncogenes’ expression with thyroid cancer was analyzed using the database of human protein atlas, UALCAN, and the UCSC Xena web tool. The heat map, hierarchical clustering, and violin plots of oncogenes were created using R language by implementing the R packages “ggplot2”, “ggpubr”, and “survival” (with a cut-off $p$ -value $<$ 0.05).

2.5 Statistical analysis

The difference in clinicopathological features among the three subtypes was analyzed using the Chi-square test. The differences in the expression levels among the three subtypes were determined using ANOVA. Two-group comparisons were performed using the student T-test. The Pearson correlation coefficient was used for correlation analysis. R was used for statistical analyses. Statistical significance was defined as ${}^{*}p<$ 0.05, ${}^{**}p<$ 0.01, and ${}^{***}p<$ 0.001. ${}^{\#}p<$ 0.05, ${}^{\#\#}p<$ 0.01, and ${}^{\#\#\#}p<$ 0.001 ${}^{\dagger}p<$ 0.05, ${}^{{\dagger}{\dagger}}p<$ 0.01, and ${}^{{\dagger}{\dagger}{\dagger}}p<$ 0.001. ${}^{\ddagger}P<$ 0.05 .

3. Results

Figure 1.

Female sex is one of the risk factors for thyroid cancer. A. Estimated number of deaths and new thyroid and other cancer types in the US in 2019. B. The distribution of thyroid cancer cases over the clinical stages. C. 5-year relative survival rate among thyroid cancer stages (SEER 17 2012–2018, all Races, both sexes by SEER combined summary of stage). D. New cases in SEER 12. Death rates from the US Mortality sheet (All races, both sexes, age-adjusted), 1975–2019. E. The rates of thyroid cancer cases and death per 100,000 persons by age. F. Female-to-male ratios of new thyroid cancer cases based on ethnicity. G. Female-to-male ratios of the death rate among thyroid cancer cases based on ethnicity. H. Male and female thyroid cancer incidence and mortality from SEER 12. Death rates from the US Mortality sheet (All races, Age-Adjusted), 1956–2019. (Data source https://seer.cancer.gov/statfacts/html/thyro.html).

3.1 Female sex is a risk factor for thyroid cancer

Table 2
Patient samples information from the SEER database

Age group	Percent of new cases (%)	Percent of deaths (%)
$<$ 20	2.0	0
20–34	16.2	0.7
35–44	18.3	1.7
45–54	21.1	6.2
55–64	21.0	17.0
65–74	14.7	27.0
75–84	5.5	28.4
$>$ 84	1.3	19.0

Data source: SEER 22 2015–2019, Age-Adjusted.

We examined the distribution of the most common cancers in the United States using the National Cancer Institute database. Thyroid cancer ranked 13 ${}^{\text{th}}$ , with estimated 43,800 new cases and 2,230 deaths in 2022 [21]. The data in the Annual Report to the Nation on the Status of Cancer (ARN) reported that the death rate of thyroid cancer was the least common in terms of estimated deaths and new cases compared to other cancers (Fig. 1A). Considering the distribution of stages, there were 65.66% localized, 29.29% regional, 3.03% distant, and 2.02% unknown cases at diagnosis (Fig. 1B). The 5-year relative survival rates were relatively high in all the stages (99.9% localized, 98.3% regional, and 90.5% unknown) except for distant tumors (53.3%) (Fig. 1C). Figure 1D illustrates long-term trends (1975–2020) in thyroid cancer incidence rate, reflecting death rate, rates of new cases, and 5-year relative survival per 100,000 persons in the US. Globally, the pattern in the incidence rate of new cases started to rise during the early 2000s, reflecting a surge in the detection of asymptomatic thyroid cancer. The overall death rate of thyroid cancer generally decreased in the early 1990s and remained stable until 2015, then decreased again during the 2020s. In contrast, the 5-year relative survival rate increased during the last 40 years, dropping slightly in 2000–2006 (Fig. 1D).

We also analyzed new and death cases as a function of age using 10-year bins. The ratio of new cases started to increase around the age of 20, reaching a peak of 21.1% between 45–65 years, then decreasing along with aging (green circle) (Fig. 1E and Table 2). In contrast, the percentage of death cases showed a distinct trend, keeping a base level during the early years and rapidly increasing between 45 and 54 years, with the highest level of 28.4% between 75–84 years (red square). These trends could be related to the menopausal transition, which most often begins between 45 and 55 years, indicating that menopause may be a detrimental factor for thyroid cancer (Fig. 1E). A sex imbalance in the overall survival (OS) and mortality from thyroid cancer patients has been reported worldwide [8, 22]. Interestingly, when analyzing the sex imbalance of thyroid cancer based on ethnicity, the female-to-male ratio of new cases per 100,000 persons was always larger than one in all ethnicities (Fig. 1F). We further analyzed the sex bias of death rate based on ethnicity, and the overall mortality rate was higher in females compared to male patients in almost all ethnicities except for Non-Hispanic Whites (Fig. 1G). Based on the 2017–2019 WHO data, the long-term global trend of thyroid cancer since 1955 was approximately 2 folds higher in women during their lifetime compared to men. Similarly, the female population had a higher mortality rate than males until 1996–2006, when the sex gap slowly narrowed and reversed. Indeed, the female mortality rate ratio declined from 1.1 in 1955 to 0.44 in 2017, and the male mortality rate increased from 0.36 in 1990 to 0.49 in 2017 (unknown reasons) (Fig. 1H).

Figure 2.

Unfavorable prognostic genes in thyroid cancer. A. Lollipop plot displaying the unfavorable prognostic gene levels (vertical axis) and their position on the chromosomes. B–U. mRNA expression levels of candidate prognostic genes in normal and thyroid tumor samples. Normal, $n=$ 59; thyroid primary tumor, $n=$ 505.

Table 3

Selected thyroid cancer genes captured by the HPA database

Gene	mRNA	$p$ -value	Chr. location	Position	Description	Predicted location
	(cancer)
ROR2	1.9	4.28E-08	9q22.31	91722601..91950228	receptor tyrosine kinase like orphan receptor 2	Membrane
PRKAB2	2.3	6.78E-07	1q21.1	147155106..147172470	protein kinase AMP-activated non-catalytic subunit beta 2	Intracellular
SNAI1	4.1	8.82E-07	20q13.13	49982980..49988886	snail family transcriptional repressor 1	Intracellular
FIBIN	2.8	1.85E-06	11p14.2	110191270..110193338	fin bud initiation factor homolog	Intracellular
PDLIM3	1.6	2.13E-06	4q35.1	185500660..185535507	PDZ and LIM domain 3	Intracellular
COPZ2	8	2.89E-06	17q21.32	48026167..48048053	COPI coat complex subunit zeta 2	Intracellular
FAM171A1	4.3	4.90E-06	10p13	15211643..15374554	family with sequence similarity 171 member A1	Intracellular, Membrane
SYNDIG1	1.6	6.16E-06	20p11.21	149672703..149846312	synapse differentiation inducing 1	Membrane
MAP1A	1.4	1.09E-05	15q15.3	43510954..43531611	microtubule associated protein 1A	Intracellular
GAS1	1.6	1.13E-05	9q21.33	887003..888682	growth arrest specific 1	Intracellular
ALDH1B1	7.9	1.17E-05	9p13.1	38392702..38398661	aldehyde dehydrogenase 1 family member B1	Intracellular
ANTXR1	11.5	1.22E-05	2p13.3	69013144..69249327	ANTXR cell adhesion molecule 1	Intracellular,Membrane
WDR37	4.5	1.34E-05	10p15.3	1056385..1132372	WD repeat domain 37	Intracellular
GUF1	4.9	1.65E-05	4p12	44678420..44700928	GUF1 homolog, GTPase	Intracellular
TSHZ1	3.6	1.71E-05	18q22.3	84029752..84105833	teashirt zinc finger homeobox 1	Intracellular
CILP	1.5	2.24E-05	15q22.31	65194760..65211473	cartilage intermediate layer protein	Secreted
RINT1	3.6	2.39E-05	7q22.3	105532201..105567677	RAD50 interactor 1	Intracellular
SENP6	8.2	2.83E-05	6q14.1	75601880..75718281	SUMO specific peptidase 6	Intracellular
PAPSS2	2.9	3.21E-05	10q23.2-q23.31	87659878..87747705	3’-phosphoadenosine 5’-phosphosulfate synthase 2	Intracellular
BCKDHA	2.3	3.27E-05	19q13.2	41397818..41425002	branched chain keto acid dehydrogenase E1 subunit alpha	Intracellular

Data source: Human Protein Atlas proteinatlas.org.

3.2 Unfavorable prognostic genes in thyroid cancer

Previous studies suggested that the severity of thyroid cancer is associated with an increased level of prognostic genes [23, 24, 25]. Here, we referred to 185 genes associated with unfavorable prognosis in thyroid cancer in the human protein atlas database (https://www.proteinatlas.org/). We examined the top 20 most significant genes related to unfavorable prognosis (Table 3). We examined the chromosome location of these genes and their fold change ( $p$ -value) at the human genome level (Fig. 2A). A 4-fold change in the expression of 8 genes was observed in thyroid cancer compared to normal samples, of which 7 genes showed an extremely low $p$ -value ( $<$ 0.00002). Among these genes, the most significant $p$ -value between normal and tumors was observed for PAPSS2 (3’-phosphoadenosine 5’-phosphosulfate synthase 2) and PDLIM3 (PDZ and LIM domain 3). Thyroid tumors showed a significant increase in the expression of ANTXR1 (ANTXR cell adhesion molecule 1), SENP6 (SUMO specific peptidase 6), COPZ2 (COPI coat complex subunit zeta 2), ALDH1B1 (aldehyde dehydrogenase 1 family member B1) and GUF1 (GUF1 homolog, GTPase) compared to normal tissues (Fig. 2A). To identify their association with tumor prognosis in thyroid cancer patients, we examined the expression level of these top 20 differentially expressed genes (DEGs) in thyroid cancer. As shown in Fig. 2B–U, dramatically decreased expression of the prognostic genes, such as COPZ2 (Fig. 2B) ( $p$ -value $=$ 1.624E-12), ALDH1B1 (Fig. 2E) ( $p$ -value $=$ 1.624E-12), ANTXR1 (Fig. 2F) ( $p$ -value $=$ 1.624E-12), GUF1 (Fig. 2H) ( $p$ -value $=$ 1.624E-12), PAPSS2 (Fig. 2L) ( $p$ -value $<$ 1E-12), PDLIM3 (Fig. 2R) ( $p$ -value $=$ 2.0220E-16), and SENP6 (Fig. 2T) ( $p$ -value $=$ 1.647E-12), was observed in primary thyroid tumors compared to their normal controls. These findings were further verified using the UALCAN database, showing consistent results (Sup. Fig. 1). These data demonstrated that COPZ2, ALDH1B1, ANTXR1, GUF1, PAPSS2, PDLIM3, and SENP6 are negatively correlated with the prognosis of thyroid cancer.

3.3 The levels of prognostic genes are dependent on the thyroid cancer

Figure 3.

The expression of prognostic genes is dependent on thyroid cancer. A–G. Levels of ALDH1B1 (A), COPZ2 (B), ANTXR1 (C), GUF1 (D), SENP6 (E), PAPSS2 (F), and PDLIM3 (G) in Solid Tissue Normal normal, Primary Tumor Solid, and Metastatic samples based on thyroid tumor types. Solid Tissue Normal normal ( $n=$ 8), Primary Tumor Solid ( $n=$ 505), Metastatic samples ( $n=$ 59); two-way ANOVA with Sidak’s post hoc correction. (Data source https://xenabrowser.net/).

It has been reported that thyroid cancer had striking female predominance, and females have three- to four-fold higher incidence than males [26]. Gender differences in thyroid cancer are mostly confined to the detection of small PTCs, and the natural causes behind gender differences in thyroid cancer are likely due to mixed effects, including sex hormones, genetics, and the immune system [27]. These issues are often associated with thyroid disorders in female patients; however, little is known about whether prognostic genes display a sexual pattern in thyroid cancer patients. We sought to determine the expression differences of prognostic genes between male and female patients and whether thyroid tumor expands sex differences in gene expression. The top 7 genes, based on the $p$ -values, were selected from the TCGA Thyroid Cancer database (THYROID CANCER), including ALDH1B1 ( $p$ -value $=$ 1.624E-12), COPZ2 ( $p$ -value $=$ 1.624E-12), ANTXR1 ( $p$ -value $=$ 1.624E-12), GUF1 ( $p$ -value $=$ 1.624E-12), SENP6 ( $p$ -value $=$ 1.647E-12), PAPSS2 ( $p$ -value $<$ 1E-12), and PDLIM3 ( $p$ -value $=$ 2.0220E-16), (Fig. 3). As shown in Fig. 2 and Sup. Fig. 1 the expression of all these genes declined in the primary tumor compared to normal tissues. Notably, when dividing the thyroid tumor into primary solid and metastatic tumors, a significant decrease in the expression of these genes was observed: ALDH1B1 (Fig. 3A) ( $p$ -value $=$ 0.0000014), COPZ2 (Fig. 3B) ( $p$ -value $=$ 0.0007694), ANTXR1 (Fig. 3C) ( $p$ -value $=$ 0.0000033), GUF1 (Fig. 3D) ( $p$ -value $=$ 0.0000212), SENP6 (Fig. 3E) ( $p$ -value $=$ 0.0668386), PAPSS2 (Fig. 3F) ( $p$ -value $=$ 0.0000001), and PDLIM3 (Fig. 3G) ( $p$ -value $=$ 0.0000002). Moreover, normal solid tissue vs. primary solid tumor displayed similar expression levels (Fig. 3A–G), which is inconsistent with the results shown in Fig. 2. The reason could be due to the limited case numbers in these subgroups (like normal solid tissue column of female in Fig. 3D, $n=$ 31).

Then, we examined the sex-related difference in the expression of these prognostic genes (PAPSS2, PDLIM3, COPZ2, ALDH1B1, ANTXR1, GUF1, and SENP6) in normal solid tissue, primary tumor, recurrent tumor, and metastatic tumor (Fig. 3A–G). No significant difference was detected; the reason may be complicated, as these genes are all located on non-sex chromosomes. According to our observations, only genes that can escape XCI may show a sexually dimorphic expression pattern [28, 29, 30]. Even though thyroid tumors showed a dramatic sex difference in the occurrence and death rate, the hormones and life habits may contribute to this difference.

Figure 4.

Prognostic gene expression in thyroid tumor stages. A–G. Levels of ALDH1B1 (A), COPZ2 (B), ANTXR1 (C), GUF1 (D), SENP6 (E), PAPSS2 (F), and PDLIM3 (G) in different stages of thyroid cancer. Normal ( $n=$ 59), Stage 1 ( $n=$ 284), Stage 2 ( $n=$ 52), Stage 3 ( $n=$ 112), and Stage 4 ( $n=$ 55). One-way ANOVA with Sidak’s post hoc correction. (Data source http://ualcan.path.uab.edu/cgi-bin/TCGAExResultNew2.pl?genenam=ALDH1B1&ctype=THCA). ${}^{*}P<$ 005 compared to normal, ${}^{***}P<$ 0.001 compared to normal, ${}^{\#}P<$ 0.05 compared to Stage 1, ${}^{{\dagger}}P<$ 0.05 compared to Stage 2, ${}^{{\dagger}{\dagger}}P<$ 0.01 compared to Stage 2, ${}^{{\dagger}{\dagger}{\dagger}}P<$ 0.001 compared to Stage 2, ${}^{{\ddagger}}P<$ 0.05 compared to Stage 3.

3.4 Prognostic gene expression in thyroid tumor stages

The correlations of prognostic genes with the pathological stages of thyroid cancer were explored to examine their clinical relevance. Our results showed that the expression of ALDH1B1 could be a reliable predictor of thyroid tumor stages, showing significant differences among all 4 stages except Stage1 vs. stage3 ( $p$ -value $=$ 9.246200E-01) (Fig. 4A). Additionally, significant correlations of prognostic genes with all 4 pathological stages were also observed, including COPZ2 (Fig. 4B), ANTXR1 (Fig. 4C), GUF1 (Fig. 4D), SENP6 (Fig. 4E), PAPSS2 (Fig. 4F), and PDLIM3 (Fig. 4G) ( $p$ -value $<$ 0.05 in normal vs. all stages). However, no significant differences could be detected among stages for ANTXR1 (Fig. 4C), SENP6 (Fig. 4E), and PDLIM3 (Fig. 4G) ( $p$ -value $>$ 0.05 in normal vs. all stages), indicating that these genes could be used as markers to identify thyroid tumors but not to predict thyroid tumor stages.

Figure 5.

Prognostic gene expression in thyroid tumor subtypes. A–G. Levels of ALDH1B1 (A), COPZ2 (B), ANTXR1 (C), GUF1 (D), SENP6 (E), PAPSS2 (F), and PDLIM3 (G) in different subtypes of thyroid cancer based on tumor histology. Normal (N) ( $n=$ 59), classical thyroid papillary carcinoma (CTPC) ( $n=$ 358), tall thyroid papillary carcinoma (TTPC) ( $n=$ 36), follicular thyroid papillary carcinoma (FTPC) ( $n=$ 102), and others (O) ( $n=$ 9). H–N. Levels of ALDH1B1 (H), COPZ2 (I), ANTXR1 (J), GUF1 (K), SENP6 (L), PAPSS2 (M), and PDLIM3 (N) in different subtypes of thyroid cancer based on immune and molecular subtypes. C1 (wound healing, $n=$ 2); C2 (IFN-gamma dominant, $n=$ 13); C3 (inflammatory, $n=$ 459); C4 (lymphocyte depleted, $n=$ 22); and C6 (TGF- $\beta$ dominant, $n=$ 3). One-way ANOVA with Sidak’s post hoc correction. (Data source http://ualcan.path.uab.edu/cgi-bin/TCGAExResultNew2.pl?genenam=ALDH1B1&ctype=THCA, http://cis.hku.hk/TISIDB/index.php). ${}^{*}P<$ 0.05 compared to normal. ${}^{**}P<$ 0.01 compared to normal. ${}^{***}P<$ 0.001 compared to normal, ${}^{\#}P<$ 0.05 compared to CTPC, ${}^{{\#}{\#}}P<$ 0.01 compared to CTPC, ${}^{{\#}{\#}{\#}}P<$ 0.001 compared to CTPC. ${}^{\dagger}P<$ 0.05 compared to TTPC, ${}^{{\dagger}{\dagger}}P<$ 0.01 compared to TTPC, ${}^{{\dagger}{\dagger}{\dagger}}P<$ 0.001 compared to TTPC, ${}^{{\ddagger}{\ddagger}}P<$ 0.01 compared to FTPC.

3.5 Expression of prognostic genes in different thyroid tumor subtypes

Determining the PTC subtype is critical in predicting patient prognosis and designing proper treatment plans. Using TCGA Thyroid Cancer data, we showed that, compared to the control tissues (N), patient tissues with classical thyroid papillary carcinoma (CTPC), tall thyroid papillary carcinoma (TTPC), follicular thyroid papillary carcinoma (FTPC), and other subtypes (O) exhibited significantly decreased expression of ALDH1B1 (Fig. 5A), ANTXR1 (Fig. 5C, GUF1 (Fig. 5D), SENP6 (Fig. 5E), PAPSS2 (Fig. 5F), and PDLIM3 (Fig. 5G) ( $p$ -value $<$ 0.05 in normal vs. all thyroid cancer subtypes), except COPZ2 (Fig. 5B) ( $p$ -value $>$ 0.05 in normal vs. other cancer subtypes). Notably, none of these genes can distinguish papillary thyroid cancer subtypes depending on their expression, which indicates that the expression of these genes had a considerable positive link to PTC diagnosis rather than distinguishing all the subtypes.

Additionally, the association of these prognostic genes among 6 immune and molecular subtypes was analyzed. The main clinical and molecular characteristics are depicted in Figs 5H–N as C1 (wound healing); C2 (IFN-gamma dominant); C3 (inflammatory); C4 (lymphocyte depleted); C5 (immunologically quiet); and C6 (TGF-b dominant). The results showed no significant difference within 5 types of immune and molecular subtypes ( $p$ -value $>$ 0.05 among different subtypes) for ALDH1B1 (Fig. 5H), COPZ2 (Fig. 5I), SENP6 (Fig. 5L), and PAPSS2 (Fig. 5M). However, the expression of ANTXR1 (Fig. 5J) was significantly different between C2 and C6 ( $p$ -value $=$ 0.01259), C3 and C4 ( $p$ -value $=$ 0.03259), C3 and C6 ( $p$ -value $=$ 0.01675), and C4 and C6 ( $p$ -value $=$ 0.0008354). In addition, GUF1 (Fig. 5K) showed a significant difference in C3 vs. C4 ( $p$ -value $=$ 0.009842), and PDLIM3 (Fig. 5N) in C3 vs. C6 ( $p$ -value $=$ 0.02517)), indicating that ANTXR1, GUF1, and PDLIM3 could be associated with immune and molecular subtypes in thyroid cancer.

Table 4
Drugs targeting ALDH1B1 and ANTXR1 from the DrugBank database

Gene	ID	Name	Drug type	Targets
ALDH1B1	DB00157	NADH	Small molecule	144
ANTXR1	DB05945	MDX-1303	Small molecule	2

Data source: https://go.drugbank.com/.

Figure 6.

The association of prognostic genes expression with thyroid cancer survival. A–G. Kaplan-Meier survival curves comparing high and low expression of ALDH1B1 (A), COPZ2 (B), ANTXR1 (C), GUF1 (D), SENP6 (E), PAPSS2 (F), and PDLIM3 (G) in thyroid cancer. The high expression level of these genes is associated with unfavorable overall survival of thyroid cancer. n (high) $=$ 252–255, n (low) $=$ 254–255. (Data source https://xenabrowser.net/).

3.6 The association of prognostic genes expression with thyroid cancer survival

We further assessed the association of the prognostic genes with the overall survival of thyroid cancer using the Kaplan-Meier curve (Fig. 6A–G). The analysis showed that shorter overall survival was significantly correlated with declined expression of ALDH1B1 (Fig. 6A), COPZ2 (Fig. 6B), ANTXR1 (Fig. 6C), PAPSS2 (Fig. 6F), and PDLIM3 (Fig. 6G) across the whole follow-up period. However, GUF1 (Fig. 6D) and SENP6 (Fig. 6E) did not affect the overall survival rate during the first 100 days. Decreased expression of these genes was associated with a lower survival rate. These observations suggest that the expression of these prognostic genes could be a detrimental factor affecting the survival of thyroid cancer.

3.7 The correlation between prognostic genes expression and immune microenvironment and modulators in thyroid cancer

Figure 7.

Correlation Between Prognostic Genes expression and Immune Microenvironment and Modulators in Thyroid Cancer. A. Spearman correlation of prognostic gene expression with immune infiltrating cells in the Tumor-Immune System Interaction Database (TISIDB). B. Spearman correlation of prognostic gene expression with chemokine expression in TISIDB. (Data source http://cis.hku.hk/TISIDB/index.php).

The immune system plays a critical role in thyroid cancer progression [31, 32, 33]. To further assess the relationship between the expression of prognostic genes in thyroid cancer and the immune microenvironment, we explored the correlations between these prognostic genes and tumor-infiltrating lymphocytes (TILs) using the TISIDB database. As shown in Fig. 7A (Columns 6 and 7), PAPSS2 and PDLIM3 expression is correlated with multiple types of TILs. Moreover, a negative correlation existed between diverse TILs and GUF1 and SENP6 expression (Fig. 7A Columns 4 and 5). These findings suggest that PAPSS2 and PDLIM3 might be associated with immune stimulation and inhibition. To further evaluate the relationship between chemokines and these prognostic genes, we evaluated their co-expression with immunomodulators (Fig. 7B). The expression of ANTXR1, PAPSS2, and PDLIM3 significantly increased compared to the rest of the candidate genes (Fig. 7B Column 6 and 5). However, the chemokines showed lower correlations with ALDH1B1, COPZ2, GUF1, and SENP6 expression.

Altogether, these results suggest that the potential roles of PAPSS2 and PDLIM3 in thyroid cancer are likely in an immune system-dependent manner.

Figure 8.

Potential anticancer drugs for thyroid cancer. A–B. Prediction of drugs targeting ALDH1B1 and ANTXR1 from the DrugBank database. C. The structure of coenzyme NADH. (Data source http://cis.hku.hk/TISIDB/browse.php?gene=ALDH1B1).

3.8 Potential anticancer drugs for thyroid cancer

Kinase inhibitors, such as Lenvatinib (Lenvima), sorafenib (Nexavar), and cabozantinib (Cabometyx) [34, 35, 36, 37], are the most common drugs used to treat thyroid cancer. New drugs are less needed unless those treatments are ineffective or may cause severe side effects. To explore the potential mechanisms involved in thyroid cancer carcinogenesis, we used the TISIDB online tool to predict the targets of drugs collected from the DrugBank database. The results are shown in Fig. 8A and Table 4. DB00157 (NADH) has a direct interaction with ALDH1B1. On the other hand, DB05945 (MDX-1303) targeted ANTXR1 only. As a coenzyme composed of ribosyl nicotinamide 5’-diphosphate, NADH could couple to adenosine 5’-phosphate by pyrophosphate linkage (Fig. 8B and C). Some evidence suggests that NADH might be useful in treating Parkinson’s disease, chronic fatigue syndrome, Alzheimer’s disease, and cardiovascular disease. At the same time, MDX-1303 (Protein Chemical Formula unavailable) is a fully human antibody against inhalation anthrax, the most lethal form of illness in humans caused by the Bacillus anthracis bacterium. It targets the anthrax protective antigen, a protein component of these lethal toxins.

These results would accelerate the understanding and future development of drugs against thyroid cancer.

4. Discussion

Thyroid cancer represents 2.3% of all new cancer cases in the US, and its incidence has risen faster over the last decades [38, 39]. This study analyzed thyroid cancer burden using data from public databases and explored the associated prognosis genes expression with thyroid tumorigenesis, stages, subtypes, and survival rate. Consistent with previous studies, we identified the female sex as a risk factor for thyroid cancer. Women were overrepresented and had a higher incidence than males in both new cases and death rates (Fig. 1F and G). A sex difference in incidence and mortality was particularly noticeable in the last 50 years (Fig. 1H).

We screened genes with the most significant association with thyroid cancer. Out of 20 genes, 7 were identified based on their fold change and $p$ -value (Fig. 2). We explored these 7 genes in association with healthy and thyroid cancer samples (sex difference) (Fig. 3), thyroid cancer stages (Fig. 4), thyroid cancer subtypes (Fig. 5), overall survival (Fig. 6), and immunological environment (Fig. 7) of thyroid cancer. Our data suggest that PAPSS2, PDLIM3, COPZ2, ALDH1B1, ANTXR1, GUF1, and SENP6 negatively correlated with the prognosis of thyroid cancer. However, no association was observed between the sex difference and these genes across all thyroid tumor subtypes. Further, gene expression analysis did not show a significant correlation between the thyroid tumor stage and the expression of ANTXR, SENP6, and PDLIM3. We also observed that ANTXR1, GUF1, and PDLIM3 were weakly associated with the immune and molecular subtypes in thyroid cancer. Moreover, their high expression positively correlated with a low survival rate in these patients.

The “female favorable” phenotype of thyroid cancer could be partly related to the higher induction of inflammatory factors and earlier activation of adaptive immunity by sex hormones. Previous clinical data have shown that the fluctuation of sex hormones may cause gender disparity in papillary thyroid cancer. Indeed, the peak incidence of PTC has been observed in women aged 40–49 years, around menopause [40, 41]. Estrogen is linked to a higher rate of thyroid cancer development in women; therefore, female sex hormones may play an important role in thyroid cancer development [42, 43, 44]. Consistent with these reports, our results showed a rapid increase in death cases around the age of 40–55 (Fig. 1E), while new cases remained stable between the ages of 40–60 years, then decreasing after that (Fig. 1E). Some studies also suggested that the X chromosome inactivation (XCI) theory, the higher expression level of X chromosome genes in women with 2 X chromosomes, can contribute to sex-specific phenotypes with more robust immune activity in many disorders including thyroid diseases and cancer [45, 46, 47]. Another unneglectable factor is the difference in lifestyle and thyroid cancer risk factors between male and female populations, such as BMI, physical activity, and tobacco smoking status. However, the association of hormones (e.g., estrogen), lifestyle, and XCI with thyroid cancer is inconclusive. None of the genes in our analysis showed a sex difference pattern between healthy tissues and tumor subtypes, which could be because these genes are not located on the X chromosome (Fig. 2A). Moreover, the function of these genes is not associated with immune processes.

Recently, Ljubic, et al. reported a set of genes and comorbidities associated with thyroid cancer based on PubMed and the HCUP SID California database. The authors identified five genes from PubMed analysis, including BRAF, RET, SLC5A5, RAS, and PTEN, and five genes using DisGeNET, including BRAF, RET, KRAS, NRAS, and PRKAR1A. However, genes that involve mutations include only BRAF (V600E) and NRAS(Q61R). These genes are generally oncogenic in various cancers, including thyroid cancer [48].

Our analysis supports the fact that the dysregulation of these prognostic genes is closely correlated with thyroid cancer malignancies. However, it should be mentioned that our study had some limitations to consider when interpreting the data. While the data from public databases and websites may validate each other, there could be some duplicates, and the sample size is insufficient. Therefore, the conclusion from sex difference patients requires further validation. We also did not have a chance to get the information on circulating sex hormone levels (estrogen and testosterone) from the database. Therefore, the contribution of the gonadal hormones to the thyroid cancer sex difference remains unknown. Of note, the conclusion drawn from this study is limited to the current information and data, and the conclusion may change as more data are collected in the future.

5. Conclusion

The present study investigated the expression variation of 7 prognostic genes in relation to sex differences, stages, and subtypes of thyroid cancer. We also analyzed the correlation of these genes in tumor immunology and their effect on the overall survival rate. Finally, we proposed potential anticancer drugs for thyroid cancer from the DrugBank database. While our study provided inconclusive results regarding the reproductive and risk factors of thyroid cancer, it may help explain the genetic susceptibility to thyroid cancer in the future.

Ethics approval and consent to participate

This study was approved by the clinical trial ethical committee of the Children’s Hospital of Fudan University. Following the approval of the ethical committee, all the methods and information in this study were performed according to the relevant guidelines and regulations.

Author contributions

Conception: Yue Huang, Yaoxin Wang, Wen-Xia Chen, and Zhengmin Xu.

Interpretation or analysis of data: Yue Huang and Yaoxin Wang.

Preparation of the manuscript: Yue Huang and Yaoxin Wang.

Revision for important intellectual content: Sining Liu, Wen-Xia Chen, and Zhengmin Xu.

Supervision: Zhengmin Xu and Wen-Xia Chen.

Funding

Not applicable.

Availability of data and materials

The original contributions presented in the study are included in the article/supplementary material. Further inquiries can be directed to the corresponding author(s).

The datasets presented in this study can be found with the following database freely. UCSC Xena (https://xenabrowser.net/), TISIDB (http://cis.hku.hk/TISIDB/index.php), cBioPortal (https://www.cbioportal.org/), STRING (http://string-db.org), and HISTome 2 (http://www.actrec.gov.in/histome2/index.php).

Code availability

Code used throughout this study is available upon reasonable request from the corresponding author(s).

Consent for publication

Not applicable.

Disclosure

Not applicable.

Supplementary data

The supplementary files are available to download from http://dx.doi.org/10.3233/CBM-230029.

sj-docx-1-cbm-10.3233_CBM-230029.docx - Supplemental material

Supplemental material, sj-docx-1-cbm-10.3233_CBM-230029.docx

Footnotes

Acknowledgments

Not applicable.

Conflict of interest

The authors have declared that no conflict of interest exists.

References

Cabanillas

M.E.

McFadden

D.G.

and Durante

, Thyroid cancer, Lancet 388(10061) (2016), 2783–2795.

Zou

, et al., Concomitant RAS, RET/PTC, or BRAF mutations in advanced stage of papillary thyroid carcinoma. Thyroid 24(8) (2014), 1256–66.

Bellevicine

, et al., Evaluation of BRAF, RAS, RET/PTC, and PAX8/PPARg alterations in different Bethesda diagnostic categories: A multicentric prospective study on the validity of the 7-gene panel test in 1172 thyroid FNAs deriving from different hospitals in South Italy. 128(2) (2020), 107–118.

Durand

, et al., Molecular characteristics of papillary thyroid carcinomas without BRAF mutation or RET/PTC rearrangement: relationship with clinico-pathological features. Endocrine-Related Cancer 16(2) (2009), 467–481.

Ming

, et al., Association between BRAF and RAS mutations, and RET rearrangements and the clinical features of papillary thyroid cancer. Int J Clin Exp Pathol 8(11) (2015), 15155–62.

Ortega

, et al., Efficacy and cost-effectiveness of the UltraCision harmonic scalpel in thyroid surgery: an analysis of 200 cases in a randomized trial. J Laparoendosc Adv Surg Tech A 14(1) (2004), 9–12.

Rahbari

Zhang

and Kebebew

, Thyroid cancer gender disparity. Future Oncol 6(11) (2010), 1771–9.

, et al., Sex disparities in thyroid cancer: a SEER population study. Gland Surg 10(12) (2021), 3200–3210.

Arnold

A.P.

, et al., Sex Hormones and Sex Chromosomes Cause Sex Differences in the Development of Cardiovascular Diseases. 37(5) (2017), 746–756.

10.

Arnold

A.P.J.H.

and behavior, The organizational – activational hypothesis as the foundation for a unified theory of sexual differentiation of all mammalian tissues. 55(5) (2009), 570–578.

11.

Schaafsma

S.M.

and Pfaff

D.W.J.F.i.n.

, Etiologies underlying sex differences in autism spectrum disorders. 35(3) (2014), 255–271.

12.

Spiering

A.E.

and de Vries

T.J.

, Why Females Do Better: The X Chromosomal TLR7 Gene-Dose Effect in COVID-19. 12 (2021).

13.

et al., Sex differences in the immune response to acute COVID-19 respiratory tract infection. Biol Sex Differ 12(1) (2021), 66.

14.

H.I.

and Freeman

E.W.

, Hormone changes associated with the menopausal transition. Minerva Ginecol 61(6) (2009), 483–9.

15.

Burger

H.G.

, et al., A review of hormonal changes during the menopausal transition: focus on findings from the Melbourne Women’s Midlife Health Project. Human Reproduction Update 13(6) (2007), 559–565.

16.

Chaligné

and Heard

, X-chromosome inactivation in development and cancer. FEBS Lett 588(15) (2014), 2514–22.

17.

Lee

J.T.

, Gracefully ageing at 50, X-chromosome inactivation becomes a paradigm for RNA and chromatin control. Nat Rev Mol Cell Biol 12(12) (2011), 815–26.

18.

Panning

, X-chromosome inactivation: the molecular basis of silencing. Journal of Biology 7(8) (2008), 30.

19.

Suteau

, et al., Sex Bias in Differentiated Thyroid Cancer. Int J Mol Sci 22(23) (2021).

20.

Sousa

, et al., Gender Differential Transcriptome in Gastric and Thyroid Cancers. 11 (2020).

21.

Siegel

R.L.

, et al., Cancer statistics, 2022. CA Cancer J Clin 72(1) (2022), 7–33.

22.

, et al., Epidemiology of Thyroid Cancer: Incidence and Mortality in China, 2015. 10 (2020).

23.

Meng

, et al., Identification of prognostic biomarkers for papillary thyroid carcinoma by a weighted gene co-expression network analysis. Cancer Med 11(9) (2022), 2006–2019.

24.

Xing

, Genetic-guided Risk Assessment and Management of Thyroid Cancer. Endocrinol Metab Clin North Am 48(1) (2019), 109–124.

25.

Ren

, et al., Identification of a Six Gene Prognosis Signature for Papillary Thyroid Cancer Using Multi-Omics Methods and Bioinformatics Analysis. 11 (2021).

26.

Sex Differences in Differentiated Thyroid Cancer. 32(3) (2022), 224–235.

27.

LeClair

, et al., Evaluation of Gender Inequity in Thyroid Cancer Diagnosis: Differences by Sex in US Thyroid Cancer Incidence Compared With a Meta-analysis of Subclinical Thyroid Cancer Rates at Autopsy. JAMA Internal Medicine 181(10) (2021), 1351–1358.

28.

, et al., X chromosome escapee genes are involved in ischemic sexual dimorphism through epigenetic modification of inflammatory signals. Journal of Neuroinflammation 18(1) (2021), 70.

29.

Yang

, et al., Global survey of escape from X inactivation by RNA-sequencing in mouse. Genome Res 20(5) (2010), 614–22.

30.

Carrel

and Willard

H.F.

, X-inactivation profile reveals extensive variability in X-linked gene expression in females. Nature 434(7031) (2005), 400–4.

31.

Galdiero

M.R.

Varricchi

and Marone

, The immune network in thyroid cancer. Oncoimmunology 5(6) (2016), e1168556.

32.

Menicali

, et al., Immune Landscape of Thyroid Cancers: New Insights. 11 (2021).

33.

Varricchi

, et al., The Immune Landscape of Thyroid Cancer in the Context of Immune Checkpoint Inhibition. 20(16) (2019), 3934.

34.

Al-Jundi

, et al., Novel Targeted Therapies for Metastatic Thyroid Cancer-A Comprehensive Review. Cancers (Basel) 12(8) (2020).

35.

Bikas

, et al., Targeted therapies in thyroid cancer: an extensive review of the literature. Expert Rev Clin Pharmacol 9(10) (2016), 1299–1313.

36.

Ferrari

S.M.

, et al., Molecular Targeted Therapies of Aggressive Thyroid Cancer. Front Endocrinol (Lausanne) 6 (2015), 176.

37.

Porter

and Wong

D.J.

, Perspectives on the Treatment of Advanced Thyroid Cancer: Approved Therapies, Resistance Mechanisms, and Future Directions. 10 (2021).

38.

Leenhardt

Grosclaude

and Chérié-Challine

, Increased incidence of thyroid carcinoma in france: a true epidemic or thyroid nodule management effects? Report from the French Thyroid Cancer Committee. Thyroid 14(12) (2004), 1056–60.

39.

Wiltshire

J.J.

, et al., Systematic Review of Trends in the Incidence Rates of Thyroid Cancer. Thyroid 26(11) (2016), 1541–1552.

40.

Kilfoy

B.A.

, et al., Gender is an age-specific effect modifier for papillary cancers of the thyroid gland. Cancer Epidemiol Biomarkers Prev 18(4) (2009), 1092–100.

41.

Negri

, et al., A pooled analysis of case-control studies of thyroid cancer. II. Menstrual and reproductive factors. Cancer Causes Control 10(2) (1999), 143–55.

42.

Moleti

, et al., Female Reproductive Factors and Differentiated Thyroid Cancer. 8 (2017).

43.

Derwahl

and Nicula

, Estrogen and its role in thyroid cancer. Endocrine-Related Cancer 21(5) (2014), T273–T283.

44.

Zane

, et al., Estrogen and thyroid cancer is a stem affair: A preliminary study. Biomed Pharmacother 85 (2017), 399–411.

45.

Santiwatana

, et al., Skewed X chromosome inactivation in girls and female adolescents with autoimmune thyroid disease. 89(6) (2018), 863–869.

46.

Ozcelik

, et al., Evidence from autoimmune thyroiditis of skewed X-chromosome inactivation in female predisposition to autoimmunity. European Journal of Human Genetics 14(6) (2006), 791–797.

47.

Chen

, et al., Investigation of the clonal origin of multifocal papillary thyroid carcinoma according to the X-chromosome inactivation pattern. Oncol Lett 17(5) (2019), 4695–4700.

48.

Ljubic

, et al., Genes and comorbidities of thyroid cancer. Informatics in Medicine Unlocked 25 (2021), 100680.

An integrative analysis of the tumor suppressors and oncogenes from sexual dimorphism and gene expression alteration features in thyroid cancer

Abstract

BACKGROUND:

OBJECTIVE:

METHODS:

RESULTS:

CONCLUSION:

Keywords

1. Introduction

Table 1 The 5- and 10-year trends of thyroid cancer during the last decade

2.1 Ethical statement

2.2 Data source and pretreatment

2.3 Distribution analysis of thyroid cancer

2.4 Association between oncogenes expression and thyroid cancer

2.5 Statistical analysis

3. Results

Table 2 Patient samples information from the SEER database

3.3 The levels of prognostic genes are dependent on the thyroid cancer

Table 4 Drugs targeting ALDH1B1 and ANTXR1 from the DrugBank database

3.7 The correlation between prognostic genes expression and immune microenvironment and modulators in thyroid cancer

4. Discussion

5. Conclusion

Ethics approval and consent to participate

Author contributions

Funding

Availability of data and materials

Code availability

Consent for publication

Disclosure

Supplementary data

sj-docx-1-cbm-10.3233_CBM-230029.docx - Supplemental material

Footnotes

Acknowledgments

Conflict of interest

References

Table 1
The 5- and 10-year trends of thyroid cancer during the last decade

Table 2
Patient samples information from the SEER database

Table 4
Drugs targeting ALDH1B1 and ANTXR1 from the DrugBank database