- Research
- Open access
- Published:
Identification of SUMOylation-related biomarkers in papillary thyroid carcinoma
Cancer Cell International volume 24, Article number: 149 (2024)
Abstract
Background
Small ubiquitin-like modifier (SUMO) modification is increasingly recognized as critical in tumorigenesis and progression. This study identifies biomarkers linked to SUMOylation in papillary thyroid carcinoma (PTC), aiming to advance therapeutic and prognostic strategies.
Methods
Employing PTC datasets and SUMO related genes (SRGs), we utilized univariate Cox regression for prognosis-related SRGs, conducted differential expression analyses, and integrated findings to pinpoint candidate genes. These genes underwent further validation through survival, gene set enrichment, immune infiltration, and drug sensitivity analyses, including external validation via quantitative RT-qPCR. In our final step, we conducted immunohistochemical staining on tumor samples from PTC patients at our center and integrated this with their clinical data to validate BMP8A’s effectiveness in predicting recurrence in PTC.
Results
Three biomarkers—BMP8A, RGS8, and SERPIND1—emerged as significant. Gene Set Enrichment Analysis (GSEA) showed their involvement in immune-related pathways, with differential immune infiltration patterns and drug response correlations observed, underscoring their potential for targeted therapy. Lastly, we validated the efficacy of BMP8A in predicting the recurrence of PTC in patients using clinical and pathological data from our center.
Conclusion
The study identifies BMP8A, RGS8, and SERPIND1 as key biomarkers associated with SUMOylation in PTC. Their linkage to immune response and drug sensitivity highlights their importance as targets for therapeutic intervention and prognosis in PTC research.
Background
Within the global context, the incidence of thyroid cancer has been incrementally rising, establishing it as the most prevalent malignancy in the head and neck region [1, 2]. Papillary thyroid carcinoma (PTC), constituting 80% of all thyroid cancer cases [3], is noted for its lower malignancy and generally favorable prognosis. However, PTC frequently undergoes lymph node metastasis in its early stages, posing significant treatment challenges [4, 5]. The primary treatment regimen includes surgery, supplemented by postoperative thyroid-stimulating hormone (TSH) suppression therapy and radioactive iodine treatment, which, while effective, does not prevent recurrence or mortality in all cases [6,7,8]. Consequently, identifying new biomarkers and understanding their mechanisms in PTC are crucial for improving early diagnosis and treatment strategies.
Small ubiquitin-like modifier modification (SUMOylation) is a pivotal post-translational modification that involves the conjugation of substrate proteins to small ubiquitin-like modifier (SUMO) proteins via a cascade of enzymatic reactions. This modification significantly influences the activity and function of substrate proteins, impacting various cellular functions [9]. SUMOylation is instrumental in regulating a myriad of biological processes such as cell proliferation, apoptosis, carcinogenesis, and DNA damage repair [10, 11]. Notably, it is intricately linked with the oncogenesis and progression of diverse cancers [12]. Studies have revealed that disruptions in SUMOylation lead to pronounced mitotic anomalies in mice [13] and increased genomic instability due to perturbations in the SUMO signaling pathway [14]. These aspects are critical in tumorigenesis. Moreover, SUMOylation is implicated in the metabolic reprogramming and growth of tumor cells; deficiencies in this pathway are associated with diminished mitochondrial oxidative phosphorylation, enhanced glycolysis, and increased lactate production, all of which contribute to tumor proliferation [15]. Recent investigations have also illuminated a significant association between SUMOylation and the anti-tumor immune response, particularly involving macrophages [16, 17]. Thus, SUMOylation emerges as a potential therapeutic target, offering new perspectives for cancer diagnosis, treatment, and prognosis. Despite its recognized importance, the specific role of SUMOylation in papillary thyroid carcinoma (PTC) remains underexplored and warrants further investigation.
Consequently, the principal objective of our current research is to identify and characterize potential SUMOylation-related biomarkers in papillary thyroid carcinoma (PTC). This involves a comprehensive analysis of the biomarkers’ functions, their implications on the immune microenvironment, as well as an examination of their mutations and drug sensitivity. By delving into these aspects, our study aims to contribute novel theoretical insights into the diagnosis and therapeutic strategies for PTC, potentially leading to more effective and personalized treatment approaches.
Materials and methods
Data source
The gene expression matrix and survival data of the The Cancer Genome Atlas (TCGA)-thyroid cancer (THCA) dataset were retrieved from the UCSC Xena platform (https://xenabrowser.net/datapages/), and then the tumor data containing only PTC was obtained by filtering the tumor data from the THCA data. A total of 560 samples were obtained from the TCGA-PTC gene expression matrix, including 502 PTC tumor samples and 58 control samples. The GSE60542 dataset was sourced from Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/), including 33 PTC tumor tissue samples and 30 normal thyroid tissue samples. In addition, totally 209 SUMOylation-related genes (SRGs) were obtained by searching “SUMOylation” in the Molecular Signatures Database (MsigDB, http://www.broadinstitute.org/gsea/msigdb/index.jsp).
Differential expression analysis
In TCGA-PTC, differential expression analysis was completed to screen differentially expressed genes (DEGs) between the PTC tumor tissue and control tissue samples by the ‘DESeq2’ R package (version 1.38.0) [18], with the threshold values of |log2FC| > 0.5 and P.adj < 0.05. The volcano plot and heatmap of DEGs were drawn using the ‘ggplot2’ R package (version 3.4.1) [19].
Identification of SUMOylation related-subtypes
To obtain prognosis-related SRGs, SRGs expression was extracted in TCGA-PTC, and then prognosis-related SRGs were screened by univariate Cox regression analysis (P < 0.05). Consensus clustering was an extensively utilized unsupervised method for classifying subtype. In this study, PTC tumor samples in TCGA-PTC were clustered by consensus clustering analysis based on prognosis-related SRGs using the ‘ConsensusClusterPlus’ R package (Version1.62.0) [20] to identify SUMOylation-related subtypes.
The dimensionality reduction analysis was finished using t-distributed stochastic neighbor embedding (t-SNE) to plot the point map of the different subtypes. Subsequently, to explore survival differences between the different subtypes, survival analysis of different subtypes was completed in TCGA-PTC. Expression heatmap of prognosis-related SRGs was drawn for different subtypes based on different clinical characteristics (age, gender, race and survival status) to demonstrate different clinical information and gene expression. Differential expression analysis was performed between the two subtypes with the most significant survival differences using the ‘DESeq2’ R package, and DEGs were screened by setting |log2FC| > 1.0 and P.adj < 0.05. To better visualize the differences in gene expression between the two subtypes, volcano plot was drawn using the ‘ggplot2’ R package.
Identification of differential genes associated with SRGs
The overlap analysis of DEGs between the PTC tumor tissue and control tissue samples and DEGs between two subtypes was performed, and the intersected genes were regarded as differential genes associated with SRGs. Subsequently, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were conducted on differential genes associated with SRGs by using ‘clusterProfiler’ R package (Version4.7.1.3) [21].
Identification of biomarkers
Differential genes associated with SRGs were screened by utilizing two machine learning algorithms. The feature genes were screened by incorporating the differential genes associated with SRGs into machine learning algorithms, including least absolute shrinkage and selection operator (LASSO) and support vector machine recursive feature elimination (SVM-RFE). The LASSO and SVM-RFE methods were performed using the ‘glmnet’ (Version4.1.4) [22] and ‘caret’ (Version6.0.93) R packages, respectively. Then, the feature genes identified from the machine learning algorithms were overlapped using ‘ggvenn’ R package (Version0.1.9) to obtain candidate genes.
Subsequently, the receiver operating characteristic (ROC) curves of the candidate genes were plotted in the TCGA-PTC and GSE60542 datasets using the ‘pROC’ R package (version1.18.0) [23], respectively. Candidate genes with area under the curve (AUC) of ROC greater than 0.8 were utilized as candidate biomarkers. Next, the different biomarkers were grouped according to their optimal cutoff values for expression in TCGA-PTC, and then the Kaplan-Meier (K-M) survival curves of different groups of PTC patients were plotted. Candidate biomarkers with P < 0.05 were screened as the biomarkers for subsequent analysis.
Gene set enrichment analysis (GSEA) of biomarkers
To further explore the potential biological functions of the biomarkers, GSEA of the biomarkers was performed. First, the samples in TCGA-PTC were categorized into high and low expression groups based on the biomarker expression. Next, differential expression analysis between high and low expression groups was accomplished using ‘DESeq2’ R package, and log2FC values was calculated. Subsequently, log2FC were sorted from large to small. Then GSEA was conducted with ‘clusterProfiler’ R package, and the reference gene set was the HALLMARK gene set in the MSigDB database. Thresholds were set at |NES| > 1 and P < 0.05.
Immune infiltration analysis of biomarkers
To further explore the relationship between biomarkers and the immune microenvironment, the relative abundance of immune-infiltrating cells was assessed for all samples in TCGA-PTC using CIBERSORT algorithm [24]. Then, Wilcoxon-test was utilized to evaluate the difference in immune cell infiltration between different samples (P < 0.05), and the immune scores of immune cells in different samples were shown by drawing box plots. Meanwhile, the correlation analysis among immune infiltrating cells was completed using the ‘corrplot’ R package. Finally, the Spearman correlation analysis was finished between immune cells and biomarkers.
Drug sensitivity analysis
Calculate the 50% inhibitory concentration (IC50) values of 138 commonly used chemotherapy and targeted therapy drugs from the Genomics of Drug Sensitivity in Cancer (GDSC) database for all TCGA-PTC patients using the “prophytic” R package (version 0.5) [25]. The Spearman correlation analysis was performed to screen drugs associated with biomarkers (|cor| > 0.5). Meanwhile, the samples in TCGA-PTC were categorized into high and low groups according to the expression of different biomarkers, and the IC50 differences of the drugs between different biomarker high and low expression groups were compared by Wilcoxon-test.
Validation of biomarker expression
To assess the robustness of biomarker expression across diverse datasets, we initially compared the expression levels of biomarkers between PTC and control samples in both the TCGA-PTC and GSE60542 datasets.
Subsequently, validation was performed using clinical specimens obtained from 20 patients who underwent surgical procedures at our institution. This included both PTC tissue and adjacent normal tissue, which were promptly stored at -80 °C.
Extract total RNA from the collected tissue using TRIzol reagent, followed by reverse transcription using PrimeScript-RT-Master Mix. Analyze the expression of target genes using the SYBR Premium Ex Taq II PCR detection system, with GAPDH serving as an internal control reference. Calculate the normalized expression level using the 2^-ΔΔCT method to determine the relative fold change. Refer to Table 1 for the primer sequences.
Validate the relationship between BMP8A expression and recurrence risk stratification based on patient clinical data
A retrospective analysis was conducted on 120 patients diagnosed with PTC at the Affiliated Hospital of Jiujiang University between September 2022 and August 2023. The inclusion criteria for the study were as follows: patients who had undergone surgery at the hospital, confirmed PTC through postoperative pathology, and had complete clinical data. Patients with other malignant tumors, postoperative pathology confirming different types of thyroid cancer, missed diagnosis, or incomplete clinical information were excluded. Subsequently, clinical data were collected, and immunohistochemical staining of tumor tissues was performed on these 120 patients.
In the immunohistochemical process, all tissue specimens were sectioned into continuous slices of 4 μm. The sections underwent a series of procedures including dewaxing, rehydration, antigen retrieval, and blocking with goat serum. Subsequently, they were incubated overnight at 4 °C with a mouse anti-human BMP8A polyclonal antibody. The sections were then washed with sterile Phosphate Buffered Saline (PBS) and incubated with the secondary antibody at room temperature for one hour. Following the addition of the substrate and counterstaining with hematoxylin, the sections were observed under a microscope. The immunohistochemical results were independently recorded and reviewed by two experienced pathologists. The presence of brown or brownish-yellow granules in the cytoplasm or on the cell membrane was marked as positive. Cells without staining were marked as negative.
The pathological staging and recurrence risk stratification of PTC patients were based on the latest AJCC guidelines for TNM staging and the ATA guidelines for recurrence risk stratification. Accordingly, we categorized the patients into low-risk and intermediate-high-risk groups.
Statistical analysis
The analysis was conducted using R software (version 3.6.1; available at https://www.r-project.org/). The relationship between categorical variables was verified using the Chi-square test. Variables deemed potentially significant in univariate analysis were included in multivariate analysis, and their significance was validated using logistic regression analysis. A p-value of less than 0.05 was considered statistically significant.
Results
DEGs1 and prognosis-related SRGs were obtained
A total of 2,984 DEGs1 (PTC vs. control) were identified in TCGA-PTC, including 1,783 up-regulated genes and 1,201 down-regulated genes, and the top 10 up-regulated and down-regulated genes were shown by volcano maps and heat maps (Fig. 1A-B). By univariate Cox regression analysis, 22 prognosis-related SRGs were screened with the threshold at P < 0.05 (Supplementary Table S1), and the genes with the top 15 h values were displayed (Fig. 1C).
SUMOylation related-subtypes were obtained through consensus clustering analysis
The results of consistent clustering analysis suggested that the PTC samples in TCGA-PTC were clustered into three SUMOylation-related subtypes based on 22 prognosis-related SRGs (cluster 1, cluster 2 and cluster 3; Fig. 2A). The point map of the different subtypes showed that the 3 different subtypes could be better distinguished from each other (Fig. 2B). The results of the survival analysis suggested that cluster 1 and cluster 2 showed more significant survival differences among the three subtypes (Fig. 2C). The heatmap of the 22 prognosis-related SRGs expression in different subtypes of clinical features was shown in Fig. 2D.
By analyzing the differentially expressed genes of two subtypes, Cluster 1 and Cluster 2, we obtained DEGs2, which consists of a total of 24 genes with differential expression. Among these, 5 genes were found to be upregulated, while 19 genes were downregulated (Fig. 2E).
Differential genes associated with SRGs were obtained
Totally 14 differential genes associated with SRGs (SFTPA2, HOXD1, TNNI3, LCN10, BMP8A, RGS8, IL1RAPL1, PENK, KIF19, GRIA3, SERPIND1, CADM2, NPY, and MROH2A) were identified by taking intersection of 2,984 DEGs1 (PTC vs. control) and 24 DEGs2 (cluster 1 vs. cluster 2), and the Venn diagram was shown in Fig. 3A.
Moreover, GO and KEGG pathway analyses of the 14 differential genes associated with SRGs were achieved to discover the potential biological functions. The GO results indicated that these genes were mainly involved in “neuronal cell body”, “positive regulation of behavior”, “receptor ligand activity” and etc. biological processes (Fig. 3B). Meanwhile, the KEGG results showed that the most significantly enriched pathways contained “cAMP signaling pathway”, “neuroactive ligand-receptor interaction”, “motor proteins”. (Figs. 3C).
BMP8A, RGS8, and SERPIND1 were biomarkers
A total of 14 differential genes associated with SRGs were submitted into LASSO, and the results showed that 12 feature genes (BMP8A, RGS8, HOXD1, TNNI3, LCN10, IL1RAPL1, PENK, GRIA3, SERPIND1, CADM2, MROH2A and KIF19) were obtained, as shown in Fig. 4A. The SVM-RFE method identified the top 13 genes with accuracy as feature genes (BMP8A, SERPIND1, RGS8, NPY, IL1RAPL1, CADM2, LCN10, HOXD1, PENK, KIF19, SFTPA2, MROH2A and TNNI3), as shown in Fig. 4B. Totally 11 candidate genes, BMP8A, RGS8, HOXD1, TNNI3, LCN10, IL1RAPL1, PENK, SERPIND1, CADM2, MROH2A and KIF19, were secured by intersecting the feature genes screened through machine learning algorithms (Fig. 4C).
Subsequently, the ROC curves of the candidate genes in the TCGA-PTC and GSE60542 datasets showed that the AUC values of BMP8A, RGS8, HOXD1, and SERPIND1 were all greater than 0.8 (Fig. 4D). Next, BMP8A, RGS8, and SERPIND1 were identified as biomarkers through the K-M survival curves, with the threshold of P < 0.05 (Fig. 4E).
Pathways of biomarkers are acquired by GSEA
The results of GSEA suggested that both BMP8A and RGS8 were enriched in “allograft rejection”, “interferon alpha response” and “interferon gamma response” (Fig. 5A and B). SERPIND1 was enriched in “kras signaling dn” (Fig. 5C) (Supplementary Table S2).
Biomarkers were associated with immune infiltrating cells
Assessment of the relative abundance of immune cells showed that abundances of M0 macrophage, resting memory CD4 + T cell and M2 macrophage were higher, abundances of neutrophil, resting mast cell and eosinophil were lower in PTC samples (Fig. 6A); while abundances of resting memory CD4 + T cell, CD8 + T cell and M0 macrophage were higher, and abundances of neutrophil, activated myeloid dendritic cell and activated memory CD4 + T cell were lower in control samples (Fig. 6B). The results of the immune score showed significant differences in 15 types of immune cells in PTC and control samples (Fig. 6C), including M1 macrophage, M2 macrophage, activated mast cell, etc. Correlation analysis among immune cells suggested that the highest positive correlation between M0 macrophage and monocyte (cor=-0.42), and the highest negative correlation between M1 macrophage and plasma B cell (cor = 0.416) were found (Fig. 6D). Meanwhile, correlation analysis between biomarkers and immune cells showed that BMP8A and RGS8 were positively correlated with eosinophil and CD8 + T cell, and negatively correlated with resting myeloid dendritic cell and regulatory T cell (Tregs); and SERPIND1 was positively correlated with M1 macrophage and resting mast cell, and negatively correlated with M0 macrophage and M2 macrophage (Fig. 6E).
The IC50 of the drugs showed differences
The correlation analysis between PTC biomarkers in the TCGA-PTC dataset and 138 drugs from the GDSC database revealed that GDC0941 exhibited the strongest correlation between BMP8A and RGS8 (Fig. 7A). Furthermore, drug sensitivity analysis results suggested that IC50s of the drugs in the high and low expression groups of each biomarker showed differences (Fig. 7B-D).
BMP8A, RGS8, and SERPIND1 were down-regulated in PTC samples
The expression validation results of the biomarkers suggested that BMP8A, RGS8, and SERPIND1 were significantly down-regulated in PTC samples and showed consistent trends in TCGA-PTC and GSE60542 datasets (Fig. 8A-B).
Subsequently, we validated the expression of three prognostic genes in 20 PTC tissues and their adjacent normal tissues through Quantitative Real-time PCR (qRT-PCR) analysis. The results also indicate that the expression of these three prognostic genes is significantly downregulated in PTC tissue (Fig. 8C-E).
Clinical data validation of the relationship between BMP8A expression and recurrence risk stratification
We validated the relationship between BMP8A expression and recurrence risk stratification in patients with papillary thyroid carcinoma (PTC) by integrating other clinicopathological indicators. The expression of BMP8A in PTC samples was confirmed using immunohistochemistry (Fig. 9). Out of 120 PTC patient samples, 97 showed negative expression of BMP8A, while 21 were positive. Subsequently, we analyzed the relationship between BMP8A expression and patient recurrence risk stratification using corresponding clinicopathological data. The results indicated a strong negative correlation between BMP8A expression and recurrence risk stratification (p = 0.005) (Table 2). This finding is consistent with our analysis results from the TCGA database, confirming the role of BMP8A expression in predicting patient disease recurrence.
Discussion
Papillary thyroid carcinoma (PTC) is generally characterized as a relatively indolent tumor, boasting long-term survival rates exceeding 95%. Nonetheless, certain aggressive forms of PTC lead to diminished disease-free and overall survival rates [26]. Notably, prevalent treatment approaches like total thyroidectomy and radioactive iodine therapy might result in overtreatment for many patients with a low-risk profile [27]. Consequently, there is an increasing emphasis on identifying novel biomarkers for PTC. A pertinent study highlighted that mRNA levels of SUMOylation-related proteins are significantly reduced in most PTC tissues compared to normal thyroid tissues [28]. Despite these initial findings, research specifically focusing on the role of SUMOylation-related biomarkers in PTC remains scant, underscoring the need for more targeted investigation in this area.
In our study, we identified three genes associated with SUMOylation modification—BMP8A, RGS8, and SERPIND1—as significant to the prognosis of papillary thyroid carcinoma (PTC) patients. Notably, patients with high expression of these genes exhibited a significantly reduced prognosis compared to those in the low-expression group. Further analysis confirmed that BMP8A, RGS8, and SERPIND1 are intricately involved in the initiation and progression of tumors, indicating their potential as prognostic markers in PTC [29,30,31].
BMP8A is part of the bone morphogenetic protein (BMP) gene family, which is known to encode secreted ligands of the transforming growth factor-β (TGF-β) protein superfamily. These ligands interact with various TGF-β receptors and activate SMAD family transcription factors, thereby influencing gene expression critical for a range of cellular functions from embryonic development and postnatal growth to tissue repair, regeneration, and organismal homeostasis [32]. The BMP signaling pathway is also implicated in inducing epithelial-mesenchymal transition, a key process in cancer progression [33]. Studies have established a strong correlation between BMP8A expression and prognosis in renal cell carcinoma and breast cancer, where high BMP8A expression typically signals a poorer prognosis [34, 35]. This effect may be due to BMP8A’s role in stimulating Nrf2 phosphorylation, thus affecting ROS balance and promoting tumor cell proliferation. In thyroid cancer research, Liu et al. [36] reported abnormal overexpression of BMP8A in thyroid cancer tissues. While our study corroborates the link between high BMP8A expression and poor prognosis in PTC patients, it interestingly notes a decreased expression of BMP8A in PTC tissues compared to normal thyroid tissues. Therefore, the precise relationship between BMP8A expression and PTC prognosis necessitates further investigation.
Genes encoding regulators of G-protein signaling (RGS) produce RGS proteins that are pivotal in modulating G-protein coupled receptors (GPCRs). The interaction between GPCRs and RGS proteins is well-established in the context of cancer onset and progression [37]. While various RGS family members have been implicated in malignant tumor development [38], information on RGS8 remains scarce. Current research has only identified its abnormal expression in brain and breast cancers [30], with the detailed mechanisms of its involvement in tumorigenesis still unclear. In the realm of thyroid cancer, Bai et al. [39] observed that RGS8 is typically downregulated in thyroid cancer tissues, a finding that aligns with the results of our study. Similar to BMP8A, there are no comprehensive studies yet that explore the relationship between RGS8 expression and the prognosis of papillary thyroid carcinoma, indicating an area ripe for further research.
Recent years have seen increasing research into the relationship between SERPIND1 and various malignant tumors. Specifically, SERPIND1 has been implicated in promoting the development of ovarian cancer by augmenting the phosphorylation of the PI3K/AKT pathway. This enhancement leads to increased proliferation, migration, invasion, and cell cycle regulation of ovarian cancer cells while concurrently inhibiting apoptosis [31]. Despite these insights, the association between SERPIND1 and thyroid cancer, particularly papillary thyroid carcinoma (PTC), remains unexplored, marking a significant gap in current thyroid cancer research.
To investigate the related pathways and biological functions of BMP8A, RGS8, and SERPIND1, we performed a comprehensive biological functional correlation analysis. The results from Gene Set Enrichment Analysis (GSEA) indicate that BMP8A and RGS8 share enrichment in several key pathways, including allograft rejection, interferon alpha response, and interferon gamma response. Conversely, SERPIND1 is predominantly associated with the kras signaling down pathway. Based on these findings, we postulate that these three genes may be intricately linked to the immune mechanisms contributing to tumor development and progression.
Immune cell infiltration is a pivotal component of the tumor microenvironment and significantly influences the initiation and progression of tumors [40]. Given the critical role of immune factors in the progression of papillary thyroid carcinoma (PTC) [41, 42], our study aimed to elucidate the correlations between BMP8A, RGS8, and SERPIND1 and immune cell infiltration. We discovered positive correlations of these genes with immune cells known for their anti-tumor properties, such as eosinophils, CD8 + T cells, M1 macrophages, and resting mast cells. Conversely, negative correlations were observed with immune cells typically associated with tumor promotion, including resting myeloid dendritic cells, regulatory T cells (Tregs), M0 macrophages, and M2 macrophages. Drawing on existing knowledge about immune cells’ roles in PTC, it is evident that eosinophils, CD8 + T cells, and M1 macrophages generally exhibit anti-tumor effects, whereas Tregs, M0 macrophages, and M2 macrophages contribute to tumor development [43]. Therefore, it can be inferred that downregulation of BMP8A, RGS8, and SERPIND1 might lead to a shift in the immune microenvironment towards a pro-tumor state, thus facilitating the occurrence and advancement of PTC.
To comprehensively assess the relationship between identified biomarkers and responses to chemotherapy and molecular targeted therapy, we determined the half-maximal inhibitory concentration (IC50) values of widely used drugs in relation to BMP8A, RGS8, and SERPIND1. Our comparative analysis identified a pronounced correlation between GDC0941, also known as Pictilisib, and the expression levels of BMP8A and RGS8. Notably, higher expression of these genes correlated with significantly improved efficacy of GDC0941. Pictilisib, a PI3K inhibitor, is instrumental in regulating key cellular functions including proliferation, differentiation, and apoptosis by suppressing PI3K activation and thus, inhibiting tumor cell proliferation and inducing apoptosis [44,45,46]. While it is predominantly used in treating breast cancer [47, 48], emerging evidence suggests its broader application across various malignancies, including thyroid cancer [49,50,51]. Notably, studies by Burrows et al. [52, 53] and Kandil et al. [54] have demonstrated Pictilisib’s efficacy in inhibiting metastasis and reducing the growth of thyroid cancer cells, respectively. These insights underscore the potential of targeting the PI3K/AKT signaling pathway in thyroid cancer treatment, offering promising avenues for future therapeutic strategies.
Subsequently, we selected BMP8A, the gene with the highest expression among the three we identified, to validate its role in predicting the prognosis of PTC patients. Since the mortality rate of PTC patients is quite low, our primary evaluation focused on the relationship between BMP8A expression and tumor recurrence in PTC patients. Currently, the risk of recurrence in PTC patients is clinically assessed based on the ATA guidelines, which categorize patients into low, intermediate, and high-risk groups. Due to the increased public awareness of health and advancements in medical imaging technology, most patients are diagnosed and treated for PTC at an early stage. Consequently, there are few high-risk patients in our center; therefore, we divided the patients into low-risk and intermediate-high-risk groups [55, 56]. We discovered that the expression of BMP8A in the tumor tissues of patients in the intermediate-high-risk group was significantly higher than that in the low-risk group. Hence, BMP8A may serve as a valuable supplement in predicting the recurrence risk of PTC patients.
This study successfully identified SUMOylation-related biomarkers in papillary thyroid carcinoma (PTC) using a series of bioinformatics approaches. It delved into their relationships with tumor subtypes, clinical features, and the immune microenvironment, laying the groundwork for new theoretical approaches in PTC research and treatment.
In subsequent clinical practice, we can initially conduct BMP8A testing on PTC patients and tailor postoperative management accordingly based on the test results. For instance, for BMP8A-positive patients, we can consider shortening their follow-up intervals, implementing more effective TSH inhibition treatment, and assessing outcomes through long-term follow-up.
However, this study is not without limitations. While the focus was specifically on PTC, other types of thyroid cancer warrant similar in-depth exploration to understand the broader implications of SUMOylation in thyroid malignancies. In this study, we only validated the association between BMP8A expression and prognosis in clinical patients, however, further validation is required for the other two genes.
Conclusions
Conclusively, SUMOylation-related biomarkers, particularly BMP8A, RGS8, and SERPIND1, demonstrate promising prognostic value, offering insights into the immune microenvironment of PTC and potential therapeutic pathways. These markers show great potential in diagnosing and predicting the prognosis of PTC patients.
Data availability
No datasets were generated or analysed during the current study.
Abbreviations
- BMP:
-
Bone morphogenetic protein
- DEGs:
-
Differentially expressed genes
- GEO:
-
Gene Expression Omnibus
- GO:
-
Gene Ontology
- GPCRs:
-
G-protein coupled receptors
- GSEA:
-
Gene set enrichment analysis
- IC50:
-
50% inhibitory concentration
- K-M:
-
Kaplan-Meier
- KEGG:
-
Kyoto Encyclopedia of Genes and Genomes
- LASSO:
-
Least absolute shrinkage and selection operator
- MSigDB:
-
Molecular Signatures Database
- PTC:
-
Papillary thyroid carcinoma
- RGS:
-
G-protein signaling
- ROC:
-
Receiver operating characteristic
- SRGs:
-
Small ubiquitin-like modifier modification -related genes
- SUMOylation:
-
Small ubiquitin-like modifier modification
- t-SNE:
-
T-distributed stochastic neighbor embedding
- TCGA:
-
The Cancer Genome Atlas
- THCA:
-
Thyroid cancer
- Tregs:
-
Regulatory T cell
- TSH:
-
Thyroid stimulating hormone
References
Siegel R, Miller K, Jemal A. Cancer statistics, 2020. CA: Cancer J Clin. 2020;70:7–30.
Chen D, Lang B, McLeod D, Newbold K. M. Haymart. Thyroid cancer. Lancet (London England). 2023;401:1531–44.
Derwich A, Sykutera M, Bromińska B, Andrusiewicz M, Ruchała M. Sawicka-Gutaj. Clinical implications of mTOR expression in papillary thyroid Cancer-A systematic review. Cancers. 2023;15:1665.
Bogdanova T, Saenko V, Brenner A, Zurnadzhy L, Rogounovitch T, Likhtarov I, Masiuk S, Kovgan L, Shpak V, Thomas G, Chanock S, Mabuchi K, Tronko M, Yamashita S. Comparative histopathologic analysis of radiogenic and sporadic papillary thyroid carcinoma: patients born before and after the chernobyl accident. Thyroid: Official J Am Thyroid Association. 2018;28:880–90.
Hou D, Xu H, Yuan B, Liu J, Lu Y, Liu M, Qian Z. Effects of active localization and vascular preservation of inferior parathyroid glands in central neck dissection for papillary thyroid carcinoma. World J Surg Oncol. 2020;18:95.
Zhang K, Li C, Liu J, Tang X, Li Z. DNA methylation alterations as therapeutic prospects in thyroid cancer. J Endocrinol Investig. 2019;42:363–70.
Saravana-Bawan B, Bajwa A, Paterson J, McMullen T. Active surveillance of low-risk papillary thyroid cancer: a meta-analysis. Surgery. 2020;167:46–55.
Hartl D, Hadoux J, Guerlain J, Breuskin I, Haroun F, Bidault S, Leboulleux S, Lamartina L. Risk-oriented concept of treatment for intrathyroid papillary thyroid cancer. Best Pract Res Clin Endocrinol Metab. 2019;33:101281.
Hendriks I, Vertegaal A. A comprehensive compilation of SUMO proteomics. Nat Rev Mol Cell Biol. 2016;17:581–95.
Eifler K, Vertegaal A. SUMOylation-mediated regulation of cell cycle progression and cancer. Trends Biochem Sci. 2015;40:779–93.
Zeng M, Liu W, Hu Y, Fu N. Sumoylation in liver disease. Clin Chim Acta. 2020;510:347–53.
Han Z, Feng Y, Gu B, Li Y, Chen H. The post-translational modification, SUMOylation, and cancer (review). Int J Oncol. 2018;52:1081–94.
Nacerddine K, Lehembre F, Bhaumik M, Artus J, Cohen-Tannoudji M, Babinet C, Pandolfi P, Dejean A. The SUMO pathway is essential for nuclear integrity and chromosome segregation in mice. Dev Cell. 2005;9:769–79.
Bergink S, Jentsch S. Principles of ubiquitin and SUMO modifications in DNA repair. Nature. 2009;458:461–7.
Shangguan X, He J, Ma Z, Zhang W, Ji Y, Shen K, Yue Z, Li W, Xin Z, Zheng Q, Cao Y, Pan J, Dong B, Cheng J, Wang Q. Xue W. SUMOylation controls the binding of hexokinase 2 to mitochondria and protects against prostate cancer tumorigenesis. Nat Commun. 2021;12:1812.
Zhang L, Xie F, Zhang J, Dijke P, Zhou F. SUMO-triggered ubiquitination of NR4A1 controls macrophage cell death. Cell Death Differ. 2017;24:1530–9.
Wang K, Zhou W, Cai Q, Cheng J, Cai R, Xing R. SUMOylation of KLF4 promotes IL-4 induced macrophage M2 polarization. Cell Cycle (Georgetown Tex). 2017;16:374–81.
Love M, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:550.
Coleman A, Bose A, Mitra S. Metagenomics data visualization using R. Methods in molecular biology. (Clifton N J). 2023;2649:359–92.
Wilkerson M, Hayes D. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinf (Oxford England). 2010;26:1572–3.
Yu G, Wang L, Han Y, He Q. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012;16:284–7.
Friedman J, Hastie T, Tibshirani R. Regularization paths for generalized Linear models via coordinate descent. J Stat Softw. 2010;33:1–22.
Robin X, Turck N, Hainard A, Tiberti N, Lisacek F, Sanchez J, Müller M. pROC: an open-source package for R and S + to analyze and compare ROC curves. BMC Bioinformatics. 2011;12:77.
Sturm G, Finotello F, List M. Immunedeconv: an R Package for unified access to computational methods for estimating immune cell fractions from bulk RNA-sequencing data. Methods in molecular biology. (Clifton N J). 2020;2120:223–32.
Geeleher P, Cox N, Huang R. pRRophetic: an R package for prediction of clinical chemotherapeutic response from tumor gene expression levels. PLoS ONE. 2014;9:e107468.
Donaldson L, Yan F, Morgan P, Kaczmar J, Fernandes J, Nguyen S, Jester R, Day T. Hobnail variant of papillary thyroid carcinoma: a systematic review and meta-analysis. Endocrine. 2021;72:27–39.
Harries V, Wang L, McGill M, Xu B, Tuttle R, Wong R, Shaha A, Shah J, Ghossein R, Patel S. Ganly I. Should multifocality be an indication for completion thyroidectomy in papillary thyroid carcinoma? Surgery. 2020;167:10–7.
Tuccilli C, Baldini E, Sorrenti S, Di Gioia C, Bosco D, Ascoli V, Mian C, Barollo S, Rendina R, Coccaro C, Pepe M, Catania A, Bononi M, Tartaglia F, De Antoni E, D’Armiento M. Ulisse S. Papillary thyroid cancer is characterized by altered expression of genes involved in the sumoylation process. J Biol Regul Homeost Agents. 2015;29:655–62.
Tang J, Tan M, Liao S, Pang M, Li J. Recent progress in the biology and physiology of BMP-8a. Connect Tissue Res. 2023;64:219–28.
Sethakorn N, Dulin N. RGS expression in cancer: oncomining the cancer microarray data. J Recept Signal Transduct Res. 2013;33:166–71.
Guo Q, Zhu L, Wang C, Wang S, Nie X, Liu J, Liu Q, Hao Y, Li X, Lin B. SERPIND1 affects the malignant biological behavior of epithelial ovarian cancer via the PI3K/AKT pathway: a mechanistic study. Front Oncol. 2019;9:954.
Sagorny K, Chapellier M, Laperrousaz B. Maguer-Satta. [BMP and cancer: the Yin and Yang of stem cells]. Med Sci (Paris). 2012;28:416–22.
Gonzalez D, Medici D. Signaling mechanisms of the epithelial-mesenchymal transition. Sci Signal. 2014;7:re8.
Yu Y, Cai L, Wang X, Cheng S, Zhang D, Jian W, Wang T, Yang J, Yang K, Zhang C. BMP8A promotes survival and drug resistance via Nrf2/TRIM24 signaling pathway in clear cell renal cell carcinoma. Cancer Sci. 2020;111:1555–66.
Katsuta E, Maawy A, Yan L, Takabe K. High expression of bone morphogenetic protein (BMP) 6 and BMP7 are associated with higher immune cell infiltration and better survival in estrogen receptor–positive breast cancer. Oncol Rep. 2019;42:1413–21.
Liu K, Gao M, Qin D, Wang H, Lu Q. Serous BMP8A has clinical significance in the ultrasonic diagnosis of thyroid cancer and promotes thyroid cancer cell progression. Endocr Metab Immune Disord Drug Targets. 2020;20:591–8.
Hurst J, Hooks S. Regulator of G-protein signaling (RGS) proteins in cancer biology. Biochem Pharmacol. 2009;78:1289–97.
Li L, Xu Q, Tang C. RGS proteins and their roles in cancer: friend or foe? Cancer Cell Int. 2023;23:81.
Bai M, Ke S, Yu H, Xu Y, Yu Y, Lu S, Wang C, Huang J, Ma Y, Dai W, Wu Y. Key molecules associated with thyroid carcinoma prognosis: a study based on transcriptome sequencing and GEO datasets. Front Immunol. 2022;13:964891.
Gajewski T, Schreiber H, Fu Y. Innate and adaptive immune cells in the tumor microenvironment. Nat Immunol. 2013;14:1014–22.
Bergdorf K, Ferguson D, Mehrad M, Ely K, Stricker T. Weiss V. Papillary thyroid carcinoma behavior: clues in the tumor microenvironment. endocrine-related cancer. 2019;26:601–14.
Yang Z, Wei X, Pan Y, Xu J, Si Y, Min Z. Yu B. A new risk factor indicator for papillary thyroid cancer based on immune infiltration. Cell Death Dis. 2021;12:51.
Xie Z, Li X, He Y, Wu S, Wang S, Sun J, He Y, Lun Y. Zhang J. Immune cell confrontation papillary thyroid carcinoma microenvironment. Front Endocrinol. 2020;11:570604.
Vanhaesebroeck B, Guillermet-Guibert J, Graupera M, Bilanges B. The emerging mechanisms of isoform-specific PI3K signalling. Nature reviews. Mol cell Biology. 2010;11:329–41.
Yue Q, Khojasteh S, Cho S, Ma S, Mulder T, Chen J, Pang J, Ding X, Deese A, Pellet J, Siebers N, Joas L, Salphati L, Ware J. Absorption, metabolism and excretion of pictilisib, a potent pan-class I phosphatidylinositol-3-Kinase (PI3K) inhibitor, in rats, dogs, and humans. Xenobiotica. 2021;51:796–810.
Sarker D, Ang J, Baird R, Kristeleit R, Shah K, Moreno V, Clarke P, Raynaud F, Levy G, Ware J, Mazina K, Lin R, Wu J, Fredrickson J, Spoerke J, Lackner M, Yan Y, Friedman L, Kaye S, Derynck M, Workman P, de Bono J. First-in-human phase I study of pictilisib (GDC-0941), a potent pan-class I phosphatidylinositol-3-kinase (PI3K) inhibitor, in patients with advanced solid tumors. Clin cancer Research: Official J Am Association Cancer Res. 2015;21:77–86.
Pictilisib stalls. Advanced ER+/PR + breast cancer. Cancer Discov. 2015;5:OF5.
Krop I, Mayer I, Ganju V, Dickler M, Johnston S, Morales S, Yardley D, Melichar B, Forero-Torres A, Lee S, de Boer R, Petrakova K, Vallentin S, Perez E, Piccart M, Ellis M, Winer E, Gendreau S, Derynck M, Lackner M, Levy G, Qiu J, He J, Schmid P. Pictilisib for oestrogen receptor-positive, aromatase inhibitor-resistant, advanced or metastatic breast cancer (FERGI): a randomised, double-blind, placebo-controlled, phase 2 trial. Lancet Oncol. 2016;17:811–21.
Khalafi S, Zhu S, Khurana R, Lohse I, Giordano S, Corso S, Al-Ali H, Brothers S, Wahlestedt C, Schürer S, El-Rifai W. A novel strategy for combination of clofarabine and pictilisib is synergistic in gastric cancer. Translational Oncol. 2022;15:101260.
De Wolf E, De Wolf C, Richardson A. ABT-737 and pictilisib synergistically enhance pitavastatin-induced apoptosis in ovarian cancer cells. Oncol Lett. 2018;15:1979–84.
Shapiro G, LoRusso P, Kwak E, Pandya S, Rudin C, Kurkjian C, Cleary J, Pilat M, Jones S, de Crespigny A, Fredrickson J, Musib L, Yan Y, Wongchenko M, Hsieh H, Gates M, Chan I, Bendell J. Phase ib study of the MEK inhibitor cobimetinib (GDC-0973) in combination with the PI3K inhibitor pictilisib (GDC-0941) in patients with advanced solid tumors. Investig New Drugs. 2020;38:419–32.
Burrows N, Babur M, Resch J, Ridsdale S, Mejin M, Rowling E, Brabant G, Williams K. GDC-0941 inhibits metastatic characteristics of thyroid carcinomas by targeting both the phosphoinositide-3 kinase (PI3K) and hypoxia-inducible factor-1α (HIF-1α) pathways. J Clin Endocrinol Metab. 2011;96:E1934–1943.
Burrows N, Telfer B, Brabant G, Williams K. Inhibiting the phosphatidylinositide 3-kinase pathway blocks radiation-induced metastasis associated with Rho-GTPase and hypoxia-inducible factor-1 activity. Radiotherapy Oncology: J Eur Soc Therapeutic Radiol Oncol. 2013;108:548–53.
Kandil E, Tsumagari K, Ma J, Abd Elmageed Z, Li X, Slakey D, Mondal D, Abdel-Mageed A. Synergistic inhibition of thyroid cancer by suppressing MAPK/PI3K/AKT pathways. J Surg Res. 2013;184:898–906.
Ywata de Carvalho A, Kohler HF, Gomes CC, Vartanian JG, Kowalski LP. Predictive factors for recurrence of papillary thyroid carcinoma: analysis of 4,085 patients. Acta Otorhinolaryngol Ital. 2021;41:236–42.
Hassan A, Razi M, Riaz S, Khalid M, Nawaz MK, Syed AA, Bashir H. Survival analysis of papillary thyroid carcinoma in relation to stage and recurrence risk: a 20-year experience in Pakistan. Clin Nucl Med. 2016;41:606–13.
Acknowledgements
Not applicable.
Funding
The authors declare that no funds, grants, or other support were received during the preparation of this manuscript.
Author information
Authors and Affiliations
Contributions
All authors contributed to the study conception and design. XL designed this work. XL and ZD integrated and analyzed the data. XL wrote this manuscript. YT provided guidance for data analysis and revised the manuscript. All authors read and approved the final manuscript.
Corresponding author
Ethics declarations
Ethics approval and consent to participate
The Affiliated Hospital of Jiujiang University’s Ethics Committee approved this study in accordance with the Declaration of Helsinki.
Consent for publication
All the authors provided consent for publication.
Competing interests
The authors declare no competing interests.
Additional information
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Electronic supplementary material
Below is the link to the electronic supplementary material.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.
About this article
Cite this article
Li, X., Ding, Z. & Tong, Y. Identification of SUMOylation-related biomarkers in papillary thyroid carcinoma. Cancer Cell Int 24, 149 (2024). https://doi.org/10.1186/s12935-024-03323-3
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s12935-024-03323-3