We collected the mRNA expression profiles and clinical data of patients with breast cancer from the Cancer Genome Atlas (TCGA) database ¹ and the Genotype Tissue Expression Project (GTEx) ² database (n = 179). The level 3 HTSeq-FPKM format data were normalized as transcripts per million reads (TPM). We also obtained the RNA-sequencing data in TPM format from the UCSC Xena ³ database and the GTEx database for pan-cancer analysis.

A total of 31 samples of paraffin-embedded breast cancer tissues and their matched paracancerous tissues were collected between April 2016 and April 2018 at the Pathology Department of Affiliated Hospital of Guilin Medical University. This study was approved by the medical ethics committees of Affiliated Hospital of Guilin Medical University (Approval Number: QTLL202136) and was conducted in line with the Declaration of Helsinki.

Immunohistochemistry (IHC) staining was carried out as previously described (DeRycke et al., 2009). Briefly, tumor tissues and paracancerous tissues were fixed in 10% formalin, paraffin-embedded, sliced into 4∼6 μm sections, and placed onto slides. After deparaffinization, rehydration and microwave antigen retrieval, the slides were incubated with MCTS1 (Abcam, Cat #ab238825) antibody at 1:800 dilution at 4°C overnight. Afterwards, the slides were incubated with secondary antibody at room temperature for 30 min and stained with DAB substrate, followed by haematoxylin counterstaining.

According to the median score of MCTS1 expression, patients with breast cancer in TCGA were divided into high and low MCTS1 expression groups. The R package DESeq2 was used to perform the differentially expressed gene (DEG) analysis between these two groups (Love et al., 2014), and adjusted p value <0.05, and |log₂-fold-change (FC)|>1 were set as the thresholds of DEGs. The correlation between the expression of the top 10 DEGs and MCTS1 was evaluated using Spearman’s correlation analysis.

Functional enrichment analyses, including Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis, were implemented for the DEGs using the R package GOplot (version 1.0.2) (Walter et al., 2015). Gene set enrichment analysis (GSEA) was carried out using the R package clusterProfiler (Subramanian et al., 2005; Yu et al., 2012), and an adjusted p value <0.05 and false discovery rate (FDR) < 0.25 were regarded as statistically significantly enriched function or pathway terms.

Based on the DEGs, a protein-protein interaction (PPI) network was established using the online STRING database ⁴ with a confidence score >0.7 and other parameters left as default, and the PPI network was visualized using the Cytoscape software (version 3.5.1) ⁵ (Szklarczyk et al., 2019). Subsequently, CytoHubba, a plugin in the Cytoscape software, was utilized to identify the top 10 hub genes of these DEGs (Chin et al., 2014).

A total of 24 immune cells were used to calculate the level of immune infiltration, and the relative enrichment score of these immune cells in breast cancer was assessed by single-sample GSEA, which was accomplished using the R package GSVA (Bindea et al., 2013). The correlation between the expression of MCTS1 and these immune cells was investigated using the Spearman’s correlation analysis, and the differences in the level of immune infiltration between the high and low MCTS1 expression groups were evaluated using the Wilcoxon rank-sum test.

To explore the underlying mechanism of MCTS1 on breast cancer, the UALCAN database ⁶ was used to investigate the status of MCTS1 promoter methylation (Chandrashekar et al., 2017). Additionally, the prognostic value of the MCTS1 methylation level was assessed using the MethSurv database ⁷ , which is an online tool for multivariable survival analysis based on DNA methylation data (Modhukur et al., 2018).

The Kaplan-Meier method with the log-rank test was used for survival analysis, and the cut-off value was set at the median expression level of MCTS1. Univariate and multivariate Cox regression analyses were used to assess the effect of clinical variables on patient outcomes. The prognostic variables p < 0.1 in the univariate Cox regression analysis were entered into multivariate Cox regression analysis. The R package ggplot2 was used to visualize the forest map.

To predict the overall survival probability, a nomogram was established based on independent prognostic factors in multivariate Cox analysis. Calibration plots were then used to assess the performance of the nomogram, and the concordance index (C-index) was used to quantify the discrimination of the nomogram. The nomogram and calibration plots were created using the R package RMS (version 5.1–4) ⁸ . The time-dependent receiver-operating characteristic (ROC) curve was performed to evaluate the predictive accuracy using the timeROC package.

All statistical analyses were conducted using R (version 3.6.3) ⁹ . The Wilcoxon rank-sum test and paired sample t-test were used to assess statistical significance for the expression of MCTS1 in the non-paired and paired tissues, respectively. The Wilcoxon rank-sum test and logistic regression were used to assess the correlations between clinical features and MCTS1 expression. All of the tests were two-sided, and p values <0.05, were regarded as statistically significant.

Our cohort included 1,065 breast cancer patients with clinical information and RNA-sequencing data, of which 110 patients had matched adjacent normal tissue samples, which were retrieved from TCGA. In addition, to increase the sample size of normal breast tissues, we obtained the gene expression data of normal breast tissues (n = 179) from the GTEx database. The clinicopathological characteristics of patients with breast cancer are shown in Supplementary Table S1.

The pan-cancer analysis showed that the expression of MCTS1 was highly expressed in most types of cancers, such as adrenocortical carcinoma, bladder urothelial carcinoma, cervical squamous cell carcinoma, adenocarcinoma, and cholangiocarcinoma (Figure 1A). The expression of MCTS1 was significantly higher in breast cancer samples than in normal breast tissues (p < 0.001) (Figure 1B). In addition, MCTS1 was highly expressed in 110 paired breast cancer tissues (p < 0.001) (Figure 1C). Furthermore, the ROC curve indicated that MCTS1 expression had good predictive power with an area under the curve (AUC) of 0.894 (95% confidence interval [CI] = 0.877–0.911) to discriminate breast cancer tissues from normal tissues (Figure 1D).

To further identify the expression of MCTS1, IHC staining was carried out on a cohort comprising 31 cases of primary breast cancer tissues paired with noncancerous tissues. There were 31 females with a mean age of 50.4 years (range 26–78 years) involved in this cohort. The expression of MCTS1 in 80.6% (25/31) of breast cancer tissues was upregulated via IHC. Representative images are presented in Figure 2.

As shown in Table 1 and Figure 3, high expression of MCTS1 was significantly associated with pathologic stage (stage IV vs. stage I, p = 0.013), M stage (p = 0.007), histological type (p < 0.001), PAM50 (Luminal B vs. Luminal A, p = 0.005; human epidermal growth factor receptor 2 [HER2] vs. Luminal A, p = 0.028), age (p = 0.002), overall survival (OS) (p < 0.001), and disease-specific survival (DSS) events (p = 0.003). Meanwhile, the results of the univariate logistic regression analyses showed that there were certain clinicopathological differences between the groups with high and low expression of MCTS1, including M stage (odds ratio [OR] = 3.048, 95% CI = 1.170–9.435, p = 0.032), age (OR = 1.281, 95% CI = 1.006–1.632, p = 0.045), race (OR = 0.667, 95% CI = 0.496–0.894, p = 0.007), PAM50 (OR = 1.463, 95% CI = 1.143–1.874, p = 0.003), and histological type (OR = 0.439, 95% CI = 0.317–0.605, p < 0.001) (Table 2).

A total of 226 genes were differentially expressed between the groups with high and low expression levels of MCTS1, including 89 upregulated DEGs (39.4%) and 137 downregulated DEGs (60.6%) (adjusted p value <0.05, |Log₂-FC| > 1) (Figure 4A and Supplementary Table S2). Next, the relationship between the top 10 DEGs (including CSN2, LALBA, SMR3B, SMR3A, FGF4, KLHL1, IAPP, CHGA, MAGEA12, and MAGEB16) and MCTS1 are presented in Figure 4B. To explore the potential interactions among all identified DEGs, we constructed a PPI network using the online STRING tool, and then identified the hub genes. As shown in Supplementary Figure S1, the network of the DEGs was complex, and the top 10 hub genes were CYP2A6, CYP2A13, ADH1B, CD19, PLIN1, LCE3A, LCE3D, LCE1B, CASP14, and LEP.

GO enrichment analysis, including biological processes, cellular compositions, and molecular functions revealed that DEGs were enriched in different GO terms such as glucuronate metabolic process, glucuronate metabolic process, cornified envelope, glucuronosyltransferase activity, and hormone activity (Figure 4C and Supplementary Table S3). Additionally, KEGG pathway analysis showed that significantly DEGs-enriched pathways included pentose and glucuronate interconversions, chemical carcinogenesis, metabolism of xenobiotics by cytochrome P450, and drug metabolism (Figure 4D and Supplementary Table S4). Subsequently, GSEA was applied between the high- and low-MCTS1 expression groups, and more immune-related biological processes were found to be significantly enriched in the low MCTS1 expression group, suggesting that the high expression of MCTS1 conferred a decreased immune phenotype in breast cancer (Figures 5A–D).

To clarify the underlying mechanisms of MCTS1 overexpression in breast cancer tissues, we also investigated the correlation between MCTS1 expression levels and methylation status using online tools. First, we observed that the breast cancer tumor tissues exhibited a significantly lower level of DNA methylation at the promoter than that in the normal breast tissues using the UALCAN database (p < 0.001) (Figure 6A). We found that the majority of methylation sites in the DNA sequences of MCTS1 were hypomethylated in breast cancer, and the degree of methylation was correlated with patient outcomes (i.e., patients with low MCTS1 methylation had poorer overall survival than patients with high MCTS1 methylation) (Figure 6B). Finally, several methylation sites were indicative of poor prognosis, including cg03900860, cg21978299, cg24931094, cg25622910, cg17246352, and cg19911179 (Figures 6C–K).

The expression of MCTS1 was significantly negatively correlated with the levels of immune cell infiltration of natural killer (NK) cells (r = –0.240, p < 0.001), CD8⁺ T cells (r = –0.220, p < 0.001), effector memory T (TEM) cells (r = –0.210, p < 0.001), and plasmacytoid dendritic cells (pDCs) (r = –0.210, p < 0.001) (Figure 7A). In addition, the enrichment scores of NK cells, CD8⁺ T cells, TEM cells, and pDCs in the MCTS1 high expression group were markedly lower than those in the MCTS1 low expression group (all p < 0.001) (Figures 7B–I).

The correlation between MCTS1 expression and the prognosis of patients with breast cancer was calculated using the Kaplan-Meier method. The median value of MCTS1 expression was used as a cut-off score, and the patients were divided into high and low MCTS1 expression groups. Compared with the low MCTS1 expression group, both the OS and DSS of the high MCTS1 expression group exhibited a significantly worse prognosis (OS: hazard ratio [HR] = 2.32, 95% CI = 1.65–3.25, p < 0.001; DSS: HR = 2.56, 95% CI = 1.61–4.06, p < 0.001) (Figures 8A,B). Additionally, in our hospital cohort, breast cancer patients with high MCTS1 expression exhibited a low disease-free survival (p = 0.023) (Supplementary Figure S2). Next, the relationships between the expression of MCTS1 and prognosis in different subgroups were evaluated.

Regardless of OS or DSS, the prognosis of patients with high MCTS1 expression was notably more unfavorable in several subgroups, including T1 and T2, N2 and N3, M0, stage II and III, age >60 years, estrogen receptor (ER)-positive, progesterone receptor (PR)-positive, HER2-positive, Luminal B, and infiltrating ductal carcinoma (IDC) subgroups (all p < 0.05) (Figure 9).

Univariate and multivariate Cox regression analyses were conducted to identify prognostic indicators. The results of the multivariate analysis demonstrated that MCTS1 expression (adjusted HR = 2.564, 95% CI = 1.346–4.883, p = 0.004), M stage (adjusted HR = 13.072, 95% CI = 1.016–168.232, p = 0.049), age (adjusted HR = 2.716, 95% CI = 1.289–5.722, p = 0.009), and menopausal status (adjusted HR = 3.345, 95% CI = 1.137–9.845, p = 0.028) were independent factors of OS in patients with breast cancer (Figure 8C and Supplementary Table S5). Similarly, for DSS, MCTS1 expression (adjusted HR = 3.267, 95% CI = 1.723–6.195, p < 0.001), M stage (adjusted HR = 10.452, 95% CI = 1.230–88.798, p = 0.032), and T3 stage (adjusted HR = 0.219, 95% CI = 0.049–0.986, p = 0.048) proved to be prognostic indicators (Supplementary Table S6).

To predict the prognosis of patients with breast cancer, a nomogram based on the independent factors of OS was generated. On the nomogram, a higher total number of points was associated with a worse prognosis (Figure 10A). Additionally, calibration curves were used to assess the prediction efficacy of the nomogram (Figures 10B–D). The bootstrap corrected C-index of the nomogram was 0.715 (95% CI = 0.687–0.743), indicating that the model had a moderate predictive accuracy for OS of patients with breast cancer. Furthermore, the time-dependent ROC curve was applied to evaluate the discriminative ability of MCTS1 expression and the constructed nomogram model, respectively (Supplementary Figures 3A–B). These results indicated that the nomogram was appropriate.