The mRNA expression of SAE1 and clinical data for TNBC patients were extracted from the cancer genome atlas (TCGA)¹ and gene expression omnibus (GEO) database² including GSE21653 (Sabatier et al., 2011b), GSE31448 (Sabatier et al., 2011a), GSE45827 (Gruosso et al., 2016), GSE53752, GSE65216 (Maire et al., 2013), GSE36295 (Merdad et al., 2014), GSE37751 (Terunuma et al., 2014), GSE38959 (Komatsu et al., 2013), and GSE61724 (Mathe et al., 2015). The method for extracting microarray gene expression values is based on our previous research (Mathe et al., 2015; Wang Q. S. et al., 2019). First, the probe ID was converted into a gene symbol. When a gene was mapped to different probes, the genic expression value was calculated by the average expression value. Next, the median normalization was performed using the robust multichip averaging method. Kaplan-Meier survival curve was conducted to estimate the overall survival (OS) and disease-free survival (DFS) by using TCGA and GSE21653, respectively.

Two hundred seven patients from Fujian Medical University Union Hospital between August 2013 and October 2017 were retrospectively enrolled in this study. All 207 patients were histologically confirmed TNBC, with a median age of 55 years (range, 27–72 years). Relevant clinicopathological information, including tumor size, nodal status, tumor grade, lymphovascular invasion (LVI), type of surgery, chemotherapy, and radiotherapy status, were obtained from patients’ medical records. OS was defined as the time from the date of diagnosis until death from any cause. DFS was defined as the time from the date of diagnosis to the date of clinical relapse (with histopathology confirmation or radiological evidence of tumor recurrence). The last follow-up was December 31, 2020. All patients met the following inclusion criteria: (1) no history of other malignant tumors or bilateral breast cancer; (2) received post-operative adjuvant chemotherapy with at least three cycles of anthracycline-based and three cycles of taxane-based regimens; (3) primary tumor size was pT1c-pT2 (1 cm < T ≤ 5 cm); and (4) demographic, clinicopathological, and follow-up information were complete. This study was approved by the Research Ethics Committee of Fujian Medical University Union Hospital and informed consent was obtained from each participant.

We conducted immunohistochemistry (IHC) staining analysis to measure the protein expression of SAE1 in TNBC tissues and adjacent normal breast tissues according to the standard immunoperoxidase staining procedure. Slides were incubated with anti-SAE1 (ab97523, Abcam, Cambridge, United Kingdom, diluted 1:300) according to the manufacturer’s instructions. A negative control was also prepared via substitution of a primary antibody with phosphate-buffered saline (PBS). The IHC staining scores of SAE1 were evaluated by two independent pathologists. The percentage of stained positive cells was calculated from 1 to 4: 1, 0–25%; 2, 26–50%; 3, 51–75%; and 4, 75–100%. The staining intensity score was ranged from 0 to 3: 0, no staining; 1, weak staining; 2, moderate staining; and 3, strong staining. The percentage of positive tumor cells and the staining intensity were multiplied to produce a weighted score for each case. A score of 8–12 was defined as high expression level and a score of 0–7 was defined as low expression.

The methylation site and beta values of SAE1 promoter in TNBC was obtained from UCSC Xena Browser³. The beta values of CpG islands had been acquired using the Illumina HumanMethylation450K platform and had been preprocessed by TCGA by using standard protocols. One of the most widely used techniques to measure DNA methylation is the Illumina Infinium HumanMethylation450 BeadChip array, which covers approximately 450,000 CpG sites. At each CpG site, methylation is quantified by the beta value beta.

Weighted gene coexpression network (WGCNA) is a common algorithm for building gene coexpression networks and performed through WGCNA R package (Langfelder and Horvath, 2008). The WGCNA hypotheses that the coexpression gene network follows the scale less distribution firstly, defines the adjoining function of the gene coexpression correlation matrix and gene network formation, calculates the difference coefficient of different nodes, and constructs a hierarchical cluster tree accordingly finally.

Gene Ontology (GO) analysis and kyoto encyclopedia of genes and genomes (KEGG) pathway enrichment were calculated by functional enrichment tool DAVID⁴ (Huang da et al., 2009). DAVID bioinformatics resources provide an integrated biological database and a repository of analytic tools for systematic exploration of biological meaning of gene set DAVID. The default parameters in the tool were used, and enriched pathways were ranked according to their enrichment scores. A P-value of <0.05 was identified as enriched functions.

CancerSEA⁵ (Yuan et al., 2019) depicts single-cell functional status maps that contain 14 functional states (including angiogenesis, apoptosis, invasion, EMT, differentiation, proliferation, DNA damage, metastasis, hypoxia, inflammation, cell cycle, DNA repair, stemness, and quiescence) from 25 types of cancer including 41,900 individual cells. CancerSEA was used to explore the potential roles of SAE1 and its coexpression genes in breast cancer.

Cistrome DB Toolkit database⁶ (Zheng et al., 2019) is a website tool that allows users to query transcription factors (TFs) that might regulate genes. We used the Cistrome DB Toolkit to predict which TFs are most likely to increase SAE1 expression in breast cancer. The correlation coefficient R > 0.3 is considered relevant.

Gene set variation analysis (GSVA) provides increased power to detect subtle pathway activity changes over a sample population in comparison with corresponding methods. In the study, the pathway activity of protein SUMOylation in TNBC were analyzed using R package “GSVA.”

In the study, t-test was applied to calculate the statistical correlation. Correlations between SAE1 expression and clinicopathological characteristics were performed by the chi-squared test. DFS and OS were calculated by the Kaplan-Meier method, and survival curves were compared via log-rank test. We conducted Cox regression analysis for the univariate and multivariate survival analyses. All P-values < 0.05 were considered statistically significant.

To examine the mRNA expression pattern of SAE1 in TNBC, publicly available datasets from TCGA and nine related GEO databases, including GSE21653, GSE31448, GSE45827, GSE53752, GSE65216, GSE36295, GSE37751, GSE38959, and GSE61724, were applied (Figures 1A–J). The mRNA level of SAE1 was found to be significantly increased in TNBC tissues (TCGA, p < 0.001; GSE21653, p < 0.001; GSE31448, p < 0.001; GSE45827, p < 0.001; GSE53752, p < 0.01; GSE65216, p < 0.001; GSE36295, p < 0.01; GSE37751, p < 0.001; GSE38959, p < 0.001; and GSE61724, p < 0.001). Moreover, meta-analysis containing 847 tissues from 10 TNBC databases mentioned above further demonstrated that the level of SAE1 mRNA was increased in TNBC (p < 0.001; Figure 1K). Next, we found that the average mRNA expression of SAE1 was higher in TNBC patients compared with luminal A or luminal B breast cancer subtypes based on GEO database GSE31448, GSE45827, and GSE65216 (Supplementary Figure 1). However, no significant differences were detected between TNBC- and HER2-enriched subtypes. These results could indicate that the increased expression of SAE1 might be more associated with ER-negative and PR-negative status.

In order to explore whether the mRNA expression of SAE1 could predict poor survival for TNBC patients, Kaplan-Meier survival analysis was performed for the OS and DFS. TNBC patients with high SAE1 expression had a significantly detrimental prognosis compared with those with low SAE1 expression (OS, p > 0.05; DFS, p < 0.001) (Figures 2A,B). Meanwhile, we performed GSVA to conduct GO analysis and assign protein SUMOlyation pathway activity estimates to individual samples and found that the activity of SUMOylation pathway was higher in TNBC patients compared with all other breast cancer subtypes (Supplementary Figure 2). The result of survival analysis indicated that high activity of SUMOylation pathway had a significantly detrimental DFS prognosis compared with those with low activity of SUMOylation pathway (Supplementary Figure 3).

To identify members of a putative molecular network that might affect SAE1 expression, we examined TFs that might regulate SAE1 gene transcription using the Cistrome DB Toolkit. We checked the TNBC cell line MDA-MB-231 and found that jun proto-oncogene (JUN), tumor protein p53 (TP53), methyl-CpG-binding domain protein 3 (MBD3), E2F transcription factor 1 (E2F1), v-myc avian myelocytomatosis viral oncogene homolog (MYC), and forkhead box M1 (FOXM1) possessed the regulatory potential in MDA-MB-231 cell (Figure 3A). Subsequent analyses revealed that SAE1 mRNA and FOXM1 mRNA expression were positively correlated in TNBC tissues using TCGA database (Figure 3B, r = 0.34, p < 0.001). However, the expression of SAE1 mRNA was not significantly correlated with JUN, TP53, MBD3, E2F1, and MYC (Figures 3C–G). The transcription factor FOXM1 coordinates the expression of cell cycle-related genes and plays a pivotal role in tumorigenesis and cancer progression. To understand the role of FOXM1 in breast cancer, gene expression data (GSE55204) from breast cancer cells before and after FOXM1 knockdown to highlight specific FOXM1 transcriptional networks were performed by Bergamaschi et al. (2014). Based on the GSE55204 database, we found that when FOXM1 was knocked down (Supplementary Figure 4A), SAE1 expression decreased (Supplementary Figure 4B). Meanwhile, the expression of FOXM1 and SAE1 was positively correlated (Supplementary Figure 4C). These experimental results further support our conclusion.

The methylation of promoter is one of the critical mechanisms regulating gene expression which has been shown to control transcription in mammals. Previous studies have revealed that increased expressions of many genes were associated with hypomethylation of promoter. The association between the SAE1 expression and methylation of SAE1 for TNBC patient was detected via UCSC Xena. The associations between methylation of 11 SAE1 CpG sites and SAE1 expression are displayed in Figure 3H. The Pearson’s correlation coefficient between methylation site (cg14042711) and SAE1 is −0.32 (p < 0.05). However, we did not detect any significant correlation between 10 other SAE1 CpG sites and SAE1 expression. These results indicate that TFs and DNA methylation modifications might play a considerable role in TNBC processes by regulating SAE1 expression.

It is well-known that tumor-infiltrating immune cells could play an important role during cancer development and are of prognostic value in TNBC patients. Therefore, we performed a correlation analysis of SAE1 expression with the quantity of 10 tumor-infiltrating immune cells (T cells, CD8 T cells, cytotoxic lymphocytes, B lineage, NK cells, monocytic lineage, myeloid dendritic cells, neutrophils, endothelial cells, and fibroblasts). However, it is indicated that the expression of SAE1 was not correlated with different immune cell infiltrates (Figures 4A–J).

To better understand the relevance and underlying mechanisms of SAE1 expression in breast cancer, we investigated the functional state of SAE1 in single-cell level via the CancerSEA database. Figure 5A displays the correlation between SAE1 and 14 functional states (including angiogenesis, apoptosis, invasion, EMT, differentiation, proliferation, DNA damage, metastasis, hypoxia, inflammation, cell cycle, DNA repair, stemness, and quiescence). Results indicated that SAE1 was positively correlated with cell cycle, DNA repair, DNA damage, and proliferation (Figures 5B–F).

To reveal the highly correlated genes and coexpression networks of SAE1 gene in TNBC patients, WGCNA analysis was conducted. Fifty-one TNBC samples obtained from GSE53752 database was used to build the WGCNA network (Figures 6A–E). We calculated the network topology for soft-thresholding powers from 1 to 30 to choose the best threshold. One of the most critical parameters in WGCNA network construction is the power value, which affects the average connectivity and independence of each coexpression module. Power value 8 was the lowest power for the scale-free topology. The coexpression similarity matrix was transformed into adjacency matrix, and the topological overlap matrix (TOM) was calculated. The dynamic tree cut analysis gave rise to 46 modules with different colors, among which the genes coexpressed with SAE1 belong to the green module containing 1,176 genes.

In order to gain further insights into the function of 1,176 genes in green module, we used DAVID database to analyze GO and KEGG pathway (Figures 6F,G). The terms of biological processes (BP) were cell-cell adhesion, intracellular protein transport, covalent chromatin modification, vesicle-mediated transport, cell division, protein autophosphorylation, protein transport, chromatin remodeling, DNA repair, mRNA processing, and protein phosphorylation. The terms of KEGG pathway terms were acute myeloid leukemia, protein processing in endoplasmic reticulum, insulin signaling pathway, endocytosis, non-small cell lung cancer, progesterone-mediated oocyte maturation, regulation of autophagy, chronic myeloid leukemia, thyroid hormone signaling pathway, and cell cycle.

Meanwhile, we analyzed the differentially expressed genes (DEGs) between TNBC tissues and adjacent normal tissues by using GSE45827, GSE65216, GSE31448, and GSE21653 databases. Figure 7A shows that the 1,176 genes in green module and DEGs in GSE45827, GSE65216, GSE31448, and GSE21653 have 10 overlap genes including minichromosome maintenance complex component 5 (MCM5), polo-like kinase 1 (PLK1), PHD finger protein 19 (PHF19), ubiquitin-conjugating enzyme E2S (UBE2S), secretion associated, Ras related GTPase 1A (SAR1A), extra spindle pole bodies homolog 1 (ESPL1), coactosin-like F-actin binding protein 1 (COTL1), antisilencing function 1B histone chaperone (ASF1B), apolipoprotein O (APOO), and transcription factor 19 (TCF19). Next, single-cell analysis was performed to investigate the functions of the 11 genes including the 10 overlap genes and SAE1 gene in breast cancer via CancerSEA website. Results indicated that these 11 genes were positively correlated with cell cycle, DNA repair, DNA damage, and proliferation (Figures 7B–F).

The prognostic analysis of the 10 overlap genes for TNBC patients were further evaluated (Figure 8). TNBC patients with low APOO expression exhibited significantly decreased DFS. TNBC patients with low MCM5, PHF19, PLK1, and UBE2S expression exhibited significantly increased DFS. However, there was not any significant relationship between the expression of ASF1B, COTL1, ESPL1, SAR1A, and TCF19 expression and the DFS of TNBC patients.

In order to construct a risk score model for predicting DFS of TNBC, LASSO Cox regression model was used to build a prognostic classifier, which included SAE1, APOO, MCM5, PHF19, PLK1, and UBE2S (Figures 9A,B). Using the LASSO Cox regression models, we calculated a risk score for each patient based on those six gene expression: risk score = (PLK1 × 0.1122) + (SAE1 × 0.3116). Survival analysis revealed that the survival time of TNBC patients in the high-risk group was significantly shorter than those patients in the low-risk group (Figure 9C). Then, the model reliability was verified through the ROC curves, and the results exhibited that the AUC value was 0.76 (Figure 9D). The above results indicated that the SAE1 and PLK1 risk assessment model had predictive value for the prognosis of TNBC patients.

Finally, the protein expression of SAE1 in TNBC was detected via IHC staining analysis. SAE1 protein expression was demonstrated significantly increased in 207 TNBC tissues compared with the paired adjacent normal breast tissues (Figures 10A–C). In our study, we found that SAE1 mainly localized to the cytoplasm but not the nucleus in TNBC patients (Supplementary Figure 5). The clinicopathological characteristics of 207 TNBC patients are presented in Table 1, and the correlations between relevant clinicopathological factors and SAE1 protein expression were analyzed. However, no associations between clinicopathological features and SAE1 expression were observed.

Survival analysis was also applied to analyze the relationship between SAE1 protein expression and the prognosis of these 207 TNBC patients. The results of Kaplan-Meier curve implied that TNBC patients with lower protein expression of SAE1 have better prognosis in both OS and DFS than those with higher expression (Figures 10D,E). Meanwhile, univariate and multivariate Cox analyses were performed to identify the independent prognostic factor for poor survival in TNBC patients (Tables 2, 3). Lymph node metastasis, tumor grade, LVI, and SAE1 expression were determined to be independent prognostic factors for OS of TNBC patients by multivariate Cox analysis. For DFS, tumor size, lymph node metastasis, tumor grade, LVI, and SAE1 expression were identified to be independent prognostic factors in TNBC patients.