Solute carrier family 31 member 1 (SLC31A1) mRNA expression data and clinicopathological characteristics in breast cancer were downloaded from the Cancer Genome Atlas (TCGA). The median expression of SLC31A1 was used to separate high/low expression groups. Human breast cancer tissues (n =  27) and normal breast tissues adjacent to the tumor (n =  27) were obtained from our breast center, Fujian Medical University Union Hospital. We perform prognostic analysis through the Kaplan–Meier Plotter website¹.

STRINGS² is an online website, which was used for protein–protein interaction (PPI) network analysis on SLC31A1. The result revealed that there were 10 proteins binding with SLC31A1, including ATP7B, MTF1, SLC22A2, ZBED3, CCS, CP, COX17, ATP7A, SLC11A2, and ATOX1. Furthermore, we, respectively, obtained the top 50 significant genes that are positively correlated and negatively correlated with SLC31A1 through the GEPIA2 website³ and the UALCAN website⁴. We took an intersection of 100 related genes between these two data sets and obtained 65 common genes (Figure 1A). Combined with 10 binding proteins and 65 related genes of SLC31A1, the functional enrichment analysis was performed by the Metascape website⁵.

Our study attempted to explore the relationship between the SLC31A1 expression and the expression of m6A-related genes in BRCA, and a correlation heat map was used to demonstrate them. All m6A-related genes analyzed are shown in Figure 2.

The starBase website⁶ is an open-source platform, which was used for decoding miRNA–ceRNA and miRNA–ncRNA interaction networks in our study (10, 11). The website integrates seven databases, containing PITA, RNA22, miRmap, microT, miRanda, PicTar, and TargetScan, and was employed to forecast upstream miRNAs potentially binding to SLC31A1. The CLIP data of the website for this study was set as strict stringency (≥5), and the other basic settings as default. In this study, at least three database-predicted miRNAs were included in the analysis. Then, we established the miRNA-SLC31A1 regulatory network with Cytoscape software. In addition, we used this website for the correlation analysis between target lncRNA and miRNA.

The website miRNet⁷ is an integrated open-source platform linking miRNAs from several miRNA-linked databases (TarBase, miRTarBase, miRecords, and miRanda). In this study, we predicted potential target lncRNAs of miRNA via this website. Subsequently, lncRNAs predicted by miRNet and starBase database were intersected to obtain the most potential regulatory lncRNAs. A comprehensive analysis of integrative miRNA–mRNA and miRNA–lncRNA was conducted with a negative correlation between expression levels to establish a key lncRNA–miRNA–mRNA (SLC31A1) ceRNA network for BRCA.

The Tumor Immune Estimation Resource (TIMER)⁸ was used for a comprehensive analysis of tumor-infiltrating immune cells in our study. Several researchers have suggested that tumor-infiltrating lymphocytes(TILs) were evaluated as a prognostic feature in multiple molecular subtypes of BRCA (12–14). Our study attempted to analyze the associations and prognosis between immune infiltrate cells (CD8+ T cells, CD4+ T cells, B cells, dendritic cells, macrophages, and neutrophils) and SLC31A1 expression by the TIMER website. In addition, we analyzed the correlation between the expression level of SLC31A1 in BRCA and the immune marker. A recent study has shown that cancer-associated fibroblast(CAF)plays a vital role in cancer progression and immunity (15). CAF is a stromal component that constitutes the tumor microenvironment (TME). To evaluate the role of SLC31A1 in the tumor microenvironment, we analyzed the correlation between the SLC31A1 expression and CAF in pan-cancer by the TIMER 2.0 website⁹.

In this study, breast cancer cell lines (MCF-7 and BT-549) and the non-tumorigenic epithelial breast cell line (MCF-10A) were purchased from the Chinese Academy of Sciences (Shanghai, China) and cultured at 37°C with 5% CO₂ in a humidified incubator. The DMEM medium (Gibco) with 10% fetal bovine serum (FBS; Gibco) was used to culture MCF-7 and BT-549 cell lines. DMEM/Ham’s F-12 (1,1; Gibco) with 5% FBS (Gibco) was used to culture MCF-10A.

Total RNA was isolated from cells using Trizol (Takara, Dalian, China) and reverse transcribed by using the PrimeScriptTM RT reagent kit (TaKaRa). Real-time PCR(RT-qPCR)was performed using the SYBR Premix Ex Taq II kit (Takara). GAPDH was used as an internal reference, and the relative expression of SLC31A1 mRNAs was quantified using the 2 − ΔΔCt method. Primer sequences were:

In this study, the Xiantao Academic website¹⁰ was used to perform most of the statistical analysis. Furthermore, we also conducted mapping and statistical analysis through GraphPad Prism software (version 7.0). Spearman’s rank correlation coefficient was applied for correlation analysis. A value of p of <0.05 was considered statistically significant. The t-test was used to compute statistical significance and annotated by the number of stars (*, value of p <0.05; **, value of p <0.01; and ***, value of p <0.001).

To determine whether SLC31A1 plays a role in human cancers, we first examined the expression level of SLC31A1 in the TCGA database. Our result suggested that the expression of SLC31A1 mRNA was upregulated in 22 tumor types while downregulated in four tumor types compared with the normal tissue (Figure 3A). In a further pairwise comparison, SLC31A1 mRNA expression was also upregulated in BRCA (Figure 3B). Furthermore, we examined the expression of SLC31A1 using real-time PCR in human breast samples, breast tumor cell lines (MCF-7 and BT-549), and human mammary epithelial cells (MCF10A). The result revealed that the mRNA expression level of SLC31A1 was higher in breast cancer tissues than in normal breast tissues (Figure 3C). Similarly, the mRNA expression level of SLC31A1 was also higher in the luminal (MCF-7) cell line and triple-negative breast tumor cell line (BT-549) than in normal breast cell lines (MCF-10A) (Figure 3D). Furthermore, we found that the expression level of SLC31A1 in the tumor tissue of all breast cancer subtypes was higher than that in the normal tissue through the UALCAN website (Figure 3E).

In our study, the correlation between the SLC31A1 expression and clinicopathological characteristics in breast cancer patients based on the TCGA database was analyzed. As shown in Table 1, the differential expression of SLC31A1 was significantly related to race (p = 0.004), ER status (p < 0.001), HER2 status (p = 0.026), and PAM50 (p < 0.001). Furthermore, the Kaplan–Meier survival analysis showed that the high expression of SLC31A1 in breast cancer patients was related to the worse relapse-free survival (RFS) (HR = 1.48, p = 2.7e-14) (Figure 4A), distant metastasis-free survival (DMFS) (HR = 1.24, p = 0.0061) (Figure 4B), and worse overall survival (OS) trend (HR = 1.19, p = 0.075) (Figure 4C).

As shown in Figure 1B, we constructed a co-expression network with 10 interacting proteins and 65 common co-expressed genes. Functional enrichment analysis indicated that these genes were significantly related to copper ion transport, body fluid secretion, positive regulation of lamellipodium assembly, and so on (Figure 1C). These results suggested that SLC31A1 interaction with these proteins may affect the transfer of copper ions and induce cuproptosis.

More recently, m6A modifications are associated with breast cancer progression (16, 17). To further assess whether SLC31A1 is related to m6A modification, we assessed the association between the expression of SLC31A1 and 20 m6A-related genes in breast cancer. Results revealed that the SLC31A1 expression was significantly positively associated with all m6A-related genes in breast cancer (Figure 2A). Importantly, the SLC31A1 expression was strongly positively associated with YTHDF3 (r = 0.479, p < 0.001), YTHDF2 (r = 0.375, p < 0.001), RBM15 (r = 0.373, p < 0.001), and YTHDF1 (r = 0.333, p < 0.001) (Figures 2B–E).

Many studies have proven that miRNA plays its biological function by participating in the regulation of the translation process of its downstream genes. A recent study showed that an ATPase copper transporter promotes breast cancer cells’ cisplatin resistance by being negatively regulated by miR-148a-3p (18). Therefore, our study attempted to predict the upstream regulatory miRNA genes of SLC31A-the copper death-related gene. As shown in Figure 5A, we finally screened out 15 predicted miRNAs that could bind to SLC31A1 through the starBase website. Furthermore, we performed a prognostic analysis of predicted miRNAs and an expression correlation analysis between miRNA and SLC31A1. Hsa-miR-28-5p, hsa-miR-29a-3p, hsa-miR-31-5p, hsa-miR-98-5p, hsa-miR-29b-3p, hsa-miR-124-3p, hsa-miR-29c-3p (Figure 5D), hsa-miR-196b-5p, and hsa-miR-543 were favorable prognostic biomarkers for OS of patients with breast cancer, while hsa-miR-105-5p, hsa-miR-219a-5p, and hsa-miR-193b-5p were unfavorable prognostic biomarkers (Figure 5B). The expression correlation between predicted miRNAs and SLC31A1 in breast cancer using the starBase database is shown in Figure 5C and Supplementary Table S1. As shown in Figure 5E, hsa-miR-29c-3p had the strongest negative correlation with SLC31A1 expression. Combined with these results, hsa-miR-29c-3p was identified to be the most potential miRNA regulating SLC31A1 in breast cancer.

A growing body of evidence has revealed that lncRNA promotes breast cancer progression and functions as ceRNA to target genes by affecting mi-RNA (6, 19, 20). Therefore, we predicted potential lncRNAs that might bind to hsa-miR-29c-3p by starBase and miRNet websites and then took the intersection (Figure 6A). We finally found 52 common predicted lncRNAs, and then the lncRNA-hsa-miR-29c-3p regulatory network was constructed using Cytoscape software (Figure 6B). The TCGA database was used to determine whether these lncRNAs were differentially expressed in breast cancer tissues compared to normal tissues. As shown in Figure 6C, a total of eight lncRNAs were significantly downregulated or upregulated in breast cancer compared to the normal tissues. The prognostic value of these eight lncRNAs in breast cancer was then evaluated. However, all eight lncRNAs show no statistically significant difference in survival(data not shown). Subsequently, the starBase database was used to evaluate the expression correlation between these eight lncRNAs and hsa-miR-29c-3p. The strongest negative expression correlation between LINC00511 and hsa-miR-29c-3p (R = −0.368, p = 4.62e-36) is shown in Figure 6D and Supplementary Table S2. Combined with the expression analysis, survival analysis, and correlation analysis, LINC00511 was considered the most significant lncRNA upstream of the hsa-miR-29c-3p/SLC31A1 axis in breast cancer.

A recent study showed that intratumoral and stromal CD4+ T cell density is an independent predictor of pCR in triple-negative patients (21). Our study attempted to investigate whether the expression of SLC31A1 is related to the level of immune infiltration in breast cancer through the TIMER website. The SLC31A1 expression showed a positive correlation with the levels of B cells (cor = 0.186, p = 4.44e-6), CD8+ T cells (cor = 0.341, p = 5.58e-28), CD4 + T cells (cor = 0.103, p = 1.32e-3), macrophages (cor = 0.24, p = 2.12e-14), neutrophils (cor = 0.291, p = 5.07e-20), and dendritic cells (DC, cor = 0.231, p = 5.84e-13) (Figure 7A). B cells, macrophages, and dendritic cells have been reported to be associated with tumor immune escape. Therefore, SLC31A1 may positively regulate B cells, macrophages, and dendritic cells to promote tumor immune escape. Furthermore, we analyzed the differential expression of 24 immune cells between the high SLC31A1 expression group and the low SLC31A1 expression group in breast cancer. Compared with the low expression group, aDC, DC, iDC, macrophage, neutrophil, T helper cells, Tcm, Tgd, Th1, and Th2 increased in the high expression group of SLC31A1, while CD8 + T cell, NK CD56bright cells, and pDC decreased (Figure 7B). In addition, we explored the correlation between SLC31A1 and immune markers of various tumor immune infiltration cells in breast cancer. The results showed that the expression of SLC31A1 was significantly positively correlated with immune markers in breast cancer, especially for STAT3 of Th17 (r = 0.382, p < 0.001), STAT1 of Th1 (r = 0.358, p < 0.001), and CCR8 of Treg (r = 0.304, p < 0.001) (Figure 7C).

Several studies have shown that CAF plays an essential role in TME and correlates with resistance and tumor progression (22–24). Therefore, we analyzed the correlation between the SLC31A1 expression and CAF by the TIMER 2.0 website. We found that the SLC31A1 expression was positively related to CAF in most tumors (Figure 8A). For BRCA-LumA disease, the SLC31A1 expression was positively associated with CAF by TIDE (Figure 8B), MCPCOUNTER (Figure 8C), and EPIC (Figure 8D) algorithms on the TIMER 2.0 website.