We used two independent data sources in this study. The breast cancer patients from the TCGA-BRCA database were considered as a constructing database. At the same time, the breast cancer patients from the GSE162228 chip were treated as a testing database.

For the constructing database, we downloaded the gene information profile and corresponding clinical information of breast cancer patients from The Cancer Genome Atlas (TCGA, https://tcga-data.nci.nih.gov/tcga/) database, which contains 1,082 breast cancer patients with complete survival information. In addition, for the testing database, we downloaded a dataset from the GEO database (GEO, https://www.ncbi.nlm.nih.gov/geoprofiles/?term=GSE162228), which contains 109 breast cancer patients with complete survival information. The genetic information of the breast cancer patients in the GEO dataset was analysed using the Affymetrix Human Genome U133A Array platform.

In this study, the identification and selection of TME-relevant prognostic genes followed the previously established protocol (14). Essentially, the complete list of genes was obtained from 10 published studies providing transcriptomic signatures for multiple immune and stromal cell populations (15–25).

The final differentially expressed TME-relevant prognostic genes for breast cancer patients were selected and verified in two steps: first, a single-factor Cox regression analysis for breast cancer prognosis was developed to analyse the expression values of each TME-relevant prognostic gene in breast cancer samples, with a threshold of p< 0.05; second, bootstrapping was performed to test the genes that passed the initial filtering for robustness as follows: 70% of patients randomly selected from the cohort were tested for the survival impact of their genes.

To avoid over-fitting the model, redundant prognosis-related molecules were removed from the dimension using LASSO Cox regression analysis. The model’s penalty parameter (lambda value) was also determined using 10-fold cross-validation, with the smallest lambda value selected to remove redundant prognosis-related molecules.

For this purpose, the glmnet function package of the R language was performed for LASSO Cox regression analysis to further infer the TME prognostic genes related to breast cancer prognosis. The following formula was established to calculate the risk score for each individual breast cancer sample:

In this formula, Coefi represents the risk coefficient of each factor estimated by the LASSO-Cox model. X_i indicates the expression activity of each TME prognostic gene. The breast cancer patients of the constructing database and the testing database could be subdivided into a high-risk group and a low-risk group based on the median of the corresponding calculated risk score.

The survival and survminer packages in the R language were used to analyse the survival of each patient. The package was based on the Kaplan–Meier method (26).

The targets of potential miRNAs were achieved via miRBase (www.http://www.mirbase.org), TargetScan (http://www.targetscan.org), miRanda (http://www.microrna.org/microrna/home.do), miRWalk (http://www.ma.uni-heidelberg.de/apps/zmf/mirwalk), TarBase (http://diana.cslab.ece.ntua.gr/tarbase), and miRecords (http://mirecords.biolead.org).

The mda-MB-453 cell line, purchased from ATCC Co., was selected as a breast cancer cell line in this study and cultured in L15 culturing medium with 10% fetal bovine serum plus 1% P/S.

The 3’ non-coding region of each TME gene was synthesised and inserted into the XhoI and NotI sites of the pCheck2 reporter luciferase vector downstream of the luciferase gene. The wild type (WT) or mutant plasmid and negative control or miRNA were transfected together into mda-MB-453 cells. Final luciferase analysis was performed using a dual luciferase reporter analysis system (Promega, USA).

The corresponding TME prognostic gene siRNAs (50 nmol/L) (OriGene) and the corresponding negative control transfection were obtained by Lipofectamine-2000 (Sigma-Aldrich, USA).

The siRNA-transfected cells were harvested and the proliferation ability was measured using the CCK-8 method (Fisher, China). At the same time, the apoptosis ability of the cells was examined by the method of flow cytometry after Annexin V FITC/PI double staining. For transwell migration assays, 2.5×10⁴ cells were seeded in the appropriate serum-free medium into the pre-coated upper chamber. The 500-μl complete medium was used as chemoattractant in the lower chambers. The incubation time was set at 48 h, after which cells without migration or invasion were removed. The final number of migrated or invaded cells was examined using Image-Pro Plus version 6.0 software. All experiments were repeated three times independently.

The multi-factor Cox regression model was built to analyze whether risk score could predict the survival of patients with breast cancer independently of all other factors. The statistical analysis was established by R software, with version number v4.2.2.

First, the breast cancer patients from the TCGA-BRCA database were used as a constructing database. We focused on TME genetic profiling in breast cancer patients. The selection scope was based on a complete list of 760 TME genes, which are listed in Supplementary Table S1 . Differential analysis was performed using single-factor Cox regression analysis. The 760 TME-relevant prognostic gene expression values were treated as continuous variables in the regression analysis. The 55 differentially expressed TME prognostic genes were finally screened out as p-value< 0.05 ( Table 1 ).

In the list, the genes with HR value less than 1 were favourable for breast cancer prognosis (protective TME prognostic genes). On the other hand, the genes with HR value greater than 1 were unfavourable for breast cancer prognosis (risk TME prognostic genes). For this purpose, 31 of the 55 genes were considered as protective genes, while the remaining 24 genes were considered as risk genes.

To remove redundant prognosis-related molecules, the 55 differentially expressed TME prognostic genes were then plotted in a LASSO Cox regression analysis. As shown in Figure 1A , the optimal number of differentially expressed TME prognostic genes was determined to be 15 ( Figure 1A , with the lowest lambda value). Therefore, the forest plot of the top 15 genes with the smallest p-value among the 55 genes is shown in Figure 1B . The 15 TME prognostic genes were ENPEP, CCDC102B, FEZ1, NOS2, SCG2, RPLP2, RELB, RGS3, EMP1, PDLIM4, EPHA3, PCDH9, VIM, GFI1, and IRF1.

The expression values of 15 TME prognostic genes were then weighted with the regression coefficients of the LASSO Cox regression model to generate a risk score for predicting survival in breast cancer patients. Each patient’s risk score was calibrated individually. Using the risk score as a further selection criterion, the 1,082 patients were divided into a high-risk group and a low-risk group based on the median of the risk score. According to the survival analysis, the patients in the high-risk sample showed an overall poor survival rate compared to the low-risk sample ( Figure 1C ).

In addition, it could be shown that the AUC of the 1-year, 3-year, and 5-year survival period of the breast cancer patients were 0.666, 0.721, and 0.681, respectively, obtained from the time-dependent ROC ( Figure 1D ). These results indicated that the risk model could accurately predict the prognosis of breast cancer patients. Meanwhile, the 15 TME prognostic genes were remarkably different when comparing the high-risk and low-risk groups ( Figure 1E ), which further confirmed the specificity of the 15 TME prognostic genes as well as the efficiency of the risk score constructed by them.

Previously, we selected primary TME prognostic genes from the constructing database and obtained the corresponding risk factor using the TME prognostic genes. We then sought to confirm the results using an independent testing database (GSE database). Out of 109 breast cancer patients (stage 0 patients were excluded) with complete clinical information, 102 were processed for the study. The age, TNM stage, and risk score of each individual breast cancer patient were all subjected to multivariate Cox regression analysis to decide whether the risk score was an independent prognostic indicator for breast cancer patients in the test database. As shown in Figure 2A , it could be demonstrated that the risk score was dramatically associated with the overall survival of the testing database, and the samples with a high risk score had a higher risk of death and were unfavourable for prognosis (HR = 3.8, 95% CI: 2.34–6, p< 0.001).

To investigate the prognostic value of the risk score established by 15 potential TME-relevant prognostic genes, we further regrouped patients from the testing database and performed Kaplan–Meier survival analysis. Patients were divided into group A ( Figure 2B , <60 years old) and group B ( Figure 2C , ≥60 years old). Patients were defined as high risk or low risk based on the risk score determined by 15 potential TME-relevant prognostic genes. Regardless of age, patients in the high-risk group had a significantly lower overall survival rate than those in the low-risk group ( Figures 2B, C ). These results concluded that the risk score constructed by primary TME prognostic genes was an independently accurate indicator for predicting the prognosis of breast cancer patients. Furthermore, the risk factor was shown to be independent of TNM stage ( Supplementary Figure 1 ).

To investigate the relationships between 15 TME prognostic genes and the prognosis of breast cancer patients in depth, the expressions of these TME prognostic genes were re-entered into the regression model with the survival probability of each individual patient from the construction database. As shown in Figure 3 , the risk TME prognostic genes were negatively associated with prognosis ( Figures 3A–E , ENPEP, CCDC102B, FEZ1, NOS2, and SCG2), while the protective TME prognostic genes were positively associated with prognosis ( Figure 3F , RPLP2). The remaining TME prognostic genes are shown in Supplementary Figure 2 . Based on the results of forest plot ( Figure 1B ) and regression model analysis ( Figure 3 ), ENPEP, CCDC102B, FEZ1, NOS2, SCG2, RGS3, EMP1, EPHA3, and PCDH9 were suggested as risk TME prognostic genes. In contrast, RPLP2, RELB, PDLIM4, VIM, GFI1, and IRF1 were suggested as protective TME prognostic genes.

Among them, there was an apparent linear internal relationship for CCDC102B ( Figure 3B ), but non-linear relationships for others ( Figures 3A, C–F ).

The 15 TME prognostic genes and associated potential target miRNAs were searched using six well-established search tools, namely, miRBase, TargetScan, miRanda, miRWalk, TarBase, and miRecords. The search results were sorted according to TarGetScore (the complete list is shown in Table 2 ).

In the list, the higher the score, the higher the confidence. The primary miRNAs with prediction scores above 80 were considered relatively reliable, while those with prediction scores below 60 were considered less reliable. The final network of 15 TME prognostic genes and targeting miRNAs was constructed using all primary miRNAs with prediction scores greater than 80 ( Figure 4A ). Notably, all miRNAs with high scores interacted with EPHA3, suggesting that the miRNA-manipulated signalling cassette may play a pivotal role in the functions of this TME prognostic gene.

To further validate our predictive results, the direct interactions between potential miRNAs and targeted TME prognostic genes were supported by the luciferase reporter assay. As shown in Figure 4B , the WT-specific fragments of the 3’ UTR of EPHA3 and mutant (MUT) were cloned into the luciferase reporter. Subsequently, hsa-miR-559 (with the highest predictive score for EPHA3) and non-targeting miRNA control were transfected into mda-MB-453 cells. According to the luciferase reporter assay results, hsa-miR-559 could directly bind to the 3’UTR of EPHA3 and inhibit its function (WT; Figure 4C ). On the other hand, the inhibition was significantly reduced when the sequence of the identified interaction site was mutated (mutant; Figure 4C ).

To date, it has been suggested that promoting cell proliferation, destroying cell apoptosis, promoting cell invasion, and promoting cell migration are the primary causes of cancerogenesis. Based on these hypotheses, we sought to investigate the functions of 15 TME prognostic genes. We used siRNA technology to knock down the expression of each TME prognostic gene. Knockdown of risk TME prognostic genes (ENPEP, CCDC102B, FEZ1, NOS2, SCG2, RGS3, EMP1, EPHA3, and PCDH9) could induce decreased cell proliferation, increased cell apoptosis, and impaired cell invasion or migration. In contrast, deletion of protective TME prognostic genes (RPLP2, RELB, PDLIM4, VIM, GFI1, and IRF1, shown in dark colour in each graph) showed the opposite pattern.

Among the 15 TME prognostic genes, most of them (NOS2, SCG2, RGS3, EMP1, PDLIM4, PCDH9, and GFI1) were involved in all three entries, while 6 of them (CCDC102B, RPLP2, RELB, EPHA3, VIM, and IRF1) were favourable for maintaining cell invasion or migration. Furthermore, only two of them (ENPEP and FEZ1) were favourable for the processes of cell proliferation and apoptosis ( Figures 5A–C ).