2.1 Data acquisition and pre-processing

The mRNA expression data and corresponding phenotypic data of 233 primary TNBC specimens were extracted from TCGA database and pre-processed for subsequent analysis. All the GSE48390, GSE58812, GSE20685, GSE162228, GSE42568, GSE20711, GSE88770, GSE61304, GSE21653, GSE97342 and GSE9893 datasets were downloaded from GEO database for analysis. The prediction of immunotherapy responsiveness and the validation of the multi-centre cohort in this study were mainly obtained from a local and web-based software tool called BEST, developed by a team from Zhengzhou University in China and some members of our team. The URL is as follows: https://rookieutopia.com/app_direct/BEST/. The gene expression data were subjected to log2 transformation (normalised RSEM count of 1), and genes with low or no expression were removed (genes with an average count of > 1 and expressed in > 75% of patients). For quality control, relative log expression (RLE) and standardised scale-free standard error (NUSE) were computed using the affyPLM package. The original gene expression data were background-corrected using the robust multi-array average (RMA) algorithm, standardised using a quantitative method and summarised using the topic polishing method. MSigDB was used to obtain genes related to cell biology. All patients with primary TNBC in TCGA and GEO datasets were selected to screen characteristic genes and establish a TNM staging prognostic/risk model. Bayes test was performed (FDR < 0.05) to determine differentially expressed genes (DEGs) among multiple samples in the GEO dataset.

2.2 Screening of differentially expressed genes

The Limma package was used to identify DEGs in TCGA and GEO datasets, and the results of mRNA sequencing of DEGs were visualised on a volcano plot. A Venn diagram and volcano plot were constructed using ggplot2. Adjusted (Adj.) p-values of < 0.05 were considered statistically significant.

2.3 Gene ontology and kyoto encyclopedia of genes and genomes pathway enrichment analyses

The cluster Profiler package was used to perform Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses. Adj. p-values of < 0.05 were considered statistically significant. DEGs were subjected to KEGG/GO enrichment analyses, and the results were visualised on a histogram and bubble plot. The gene set enrichment analysis (GSEA) method was used to identify pathways associated with the expression of DEGs, and the groups were arranged 10,000 times in each analysis. Adj. p-values of < 0.05 and FDR values of < 0.25 indicated significant differences, and ‘c2.cp.kegg.v7.0.symbols.gmt’ was used as the reference gene set. The results of enrichment analysis were characterised based on adj. p-values and NESs.

2.4 Prognostic markers screened via Lasso–Cox regression

The glmnet package was used to implement Lasso–Cox regression analysis to determine genes associated with the prognosis of TNBC. Genes identified as covariates in univariate analysis were included in multivariate Cox regression analysis to determine their impact on overall survival (OS). A prognostic model was established by integrating the expression levels of the selected genes and the multivariate Cox regression coefficients. The model was used to stratify patients with TNBC using the median risk score as the threshold. Patients with a risk score greater than the median risk score were assigned to the high-risk group, whereas the remaining patients were assigned to the low-risk group. The log-rank test was performed to evaluate differences in the survival rate between the two groups, and a calibration curve was plotted using the fitplot package.

2.5 Estimation of the risk model

A time-dependent receiver operating characteristic (ROC) curve was plotted to evaluate the predictive accuracy of the risk model. The sensitivity and specificity of predicting the survival rate of patients with TNBC were calculated at 1-, 3- and 5-year follow-up intervals. Decision curve analysis (DCA) was used to assess the predictive accuracy of the risk model, and the results were corroborated using a calibration curve. The risk score and other clinical indicators, such as age, sex, TNM stage and histological grade, were also observed. The distribution of death events was demonstrated on a point map based on the increase in the risk score. A heatmap was constructed to demonstrate the distribution of the expression of each characteristic gene in the high- and low-risk groups. To assess the independent predictive performance and reliability of the prognostic model, common indicators such as age, TNM stage and histological grade were used to predict the prognosis of TNBC. Univariate and multivariate Cox regression analyses were performed based on the clinical characteristics and risk scores of patients with TNBC in TCGA cohort to determine clinical factors associated with survival. The log-rank test was performed to verify whether the risk score was related to other survival-related clinical characteristics. Subsequently, variables that were identified as independent prognostic factors were used to establish a nomogram. By intersecting the above prognostic genes in triple-negative breast cancer, we analyzed OS, RFS, PFS, DSS from GSE58812, GSE88770, GSE20685, GSE42568, GSE97342, GSE162228, GSE48390, GSE58812, GSE88770, GSE20711, GSE42568, GSE162228, GSE9893 datasets in GEO database.

2.6 Assessment of immune infiltration

The CIBERSORT algorithm (URL: https://cibersort.stanford.edu/) was used to evaluate the infiltration of immune cells based on the RNA-Seq data of TNBC tissues. CIBERSORT provides gene expression signatures of 24 immune cell types, such as B cells, CD4+T cells, CD8+T cells, neutrophils, macrophages and dendritic cells.

2.7 Analysis of tumor burden in triple-negative breast cancer

The cBioPortal database is a molecular dataset that helps to understand genetics, epigenetics, gene expression and proteomics for oncology and cytology research [17]. It provides an interactive interface for custom data that allows researchers to explore the relationship between genetic changes and clinical research. The mutation frequency and percentage of key prognostic genes in triple-negative breast cancer were analysed, and clinical characteristics, such as clinical stage, pathological grade and expression of tumor-specific markers, were shown through the public database (http://www.cbioportal.org/).

2.8 Histochemical analysis of key prognostic genes

The Human Protein Atlas (HPA) database is based on proteomics, transcriptomics and systems biology data. It can be used to map tissues, cells and organs and contains the protein expression data of not only tumor tissues but also normal tissues. In addition, it provides survival curves for patients with tumors. In this study, the HPA database was used to evaluate the expression of the identified prognostic genes in TNBC based on histological staining.

2.9 WGCNA

The gene expression data of patients with TNBC were extracted from TCGA. The median absolute deviation (MAD) of each gene was calculated, and the top 50% of genes with the smallest MAD were removed. Subsequently, the good Samples Genes tool in the R package WGCNA was used to eliminate outliers and construct a scale-free co-expression network. To construct the network, all pairwise genes were analysed via Pearson correlation analysis and the average linkage method. Subsequently, a weighted adjacency matrix was constructed using a power function: A_mn =|C_mn|^β (C_mn = Pearson correlation coefficient between Gene_m and Gene_n; A_mn = adjacency between Gene m and Gene n). In this equation, β represents a soft-thresholding parameter that emphasises strong correlations between genes and penalises weak correlations. After selecting the power of 7, the adjacency matrix was transformed into a topological overlap matrix (TOM) that measured the network connectivity of a gene defined as the sum of its adjacency with all other genes for generating a network, and the corresponding dissimilarity (1-TOM) was calculated. To classify genes with similar expression profiles into gene modules, average linkage hierarchical clustering was performed according to TOM-based dissimilarity, with a minimum number (gene group) set to 30 genes for constructing a dendrogram. For further analysis of these modules, the dissimilarity of module eigengenes was calculated, a cut-off line was selected for the module dendrogram and some modules were merged. Additionally, modules with distances less than 0.25 were merged, and 9 co-expression modules were eventually obtained. Notably, a gene set that could not be assigned to any module was included in the grey module. To screen for hub genes, the correlation between the module feature vector and gene expression was evaluated to obtain MM. Highly correlated genes within a clinically significant module were identified as hub genes based on the cut-off criteria (|MM|> 0.8).

2.10 Cell culture and transfection

The murine TNBC cell line 4T1 and the human TNBC cell line MDA-MB-231 were obtained from the Type Culture Collection Committee of the Chinese Academy of Sciences (Shanghai, China). The cell lines were authenticated through short tandem repeat analysis. Mycoplasma contamination was not observed in the cell lines. 4T1 and MDA-MB-231 cells were cultured in RPMI-1640 medium and Dulbecco’s modified Eagle’s medium (DMEM) supplemented with 10% foetal bovine serum (FBS), respectively. The cultures were maintained in a humidified atmosphere of 95% air and 5% CO2 at 37 °C.

The shRNA sequences of the ZDHHC9 gene were synthesised and cloned into the pLV-H1-EF1α-puro vector. The lentivirus was generated in 293 T cells and transfected into 4T1 cells. Of the two stable cell lines, the most efficient cell line was used for subsequent experiments. shRNA sequences are listed in Supplementary Table 2.

2.11 Patients and tissue samples

A total of 114 TNBC tissues were retrospectively collected from patients with a histological diagnosis of TNBC who underwent modified radical mastectomy between January 2010 and January 2012 at the Tianjin Medical University Cancer Institute and Hospital, China. Data regarding the clinicopathological features of patients were collected retrospectively. In addition, eight pairs of fresh TNBC and para-cancerous tissues were collected during surgery between June and December 2022. Total protein was extracted from the tissues for detecting ZDHHC9 expression via western blotting. The use of these samples and patient information was approved by the Tianjin Medical University Cancer Institute and the ethics committee of the hospital in accordance with the ethical standards as laid down in the 1964 Declaration of Helsinki and its later amendments.

2.12 Immunohistochemical analysis

After de-paraffinisation and hydration, TNBC tissue sections were subjected to antigen retrieval under alkaline conditions (EDTA, pH 9.0) and incubated with a primary antibody against ZDHHC9 (Proteintech, 24046-1-AP, 1:200) overnight at 4 ℃. The following day, the sections were treated with HRP-conjugated secondary antibodies at 37 ℃ for 1 h, and the bound secondary antibodies were stained with DAB (ZSGB-BIO, ZLI-9018). The final staining scores were calculated by multiplying the score for staining intensity by the score for staining range. Staining intensity was scored as follows: 0, negative; 1, low; 2, medium; 3, high. Staining range was scored as follows: 0, 0% staining; 1, 1–25% staining; 2, 26–50% staining; 3, 51–100% staining. The final score ranged from 0 to 9: < 2, negative (−); 2–3, low staining ( +); 4–6, moderate staining (+ +); > 6, high staining (+ +  + +). Five areas were randomly examined using a light microscope (100 × magnification). The samples were independently analysed by two pathologists who were informed of the clinical features and examination findings of the patients.

2.13 Cellular immunofluorescence analysis

Cells (MDA-MB-231 and 4T1) were fixed with 4% fresh paraformaldehyde, permeabilised with 0.2% Triton X-100 and incubated with 5% goat serum for 1 h to block non-specific interactions. Subsequently, the cells were incubated with primary antibodies against ZDHHC9 (Proteintech, 24046-1-AP, 1:100) and F-actin (Abcam, ab130935, 1:200) overnight at 4 ℃. The following day, the cells were incubated with fluorescently labelled secondary antibodies (1:200) (Thermo Fisher Scientific, USA) for 1 h, mounted with a DAPI-containing antifade mounting medium (Vector Laboratories) and sealed with clear nail polish. The cells were extensively washed with PBS after each step. Finally, images were captured using a confocal fluorescence microscope (ZEISS, Germany).

2.14 Tissue immunofluorescence analysis

Two sets of 114 PDAC tissues were used for immunological assessment of ZDHHC9 and CD8/NCR1 expression. The tissues were stained using the Opal 7-Color Manual IHC Kit according to the manufacturer’s instructions. Briefly, after de-paraffinisation and hydration, TNBC tissue sections were fixed via microwave irradiation for 20 min under appropriate repair conditions and incubated with primary antibodies diluted in a blocking solution overnight at 4 °C (Supplementary Table 3). The following day, the tissues were incubated with HRP-conjugated secondary antibodies for 10 min at room temperature and stained with a fluorescent dye for 10 min at room temperature. The steps of microwave-assisted fixation, primary and secondary antibody incubation and cross-linking of fluorescent signals were repeated until the sections were stained for all indicators. The nuclei were stained with DAPI, and the tissue slices were sealed. The samples were washed with TBST after each step. Finally, images were captured using a fluorescent microscope (KEYENEC, Japanese).

Tissue immunofluorescence (IF) values were independently determined by two pathologists who were under a duty of confidentiality regarding the clinical features, examination findings and IHC analysis results of patients. Five randomly selected areas were evaluated microscopically (200 × magnification). CD8 + and NCR1 + cells in the tumor stroma were used to visualise CD8 + T cells, and NK cells, respectively. Based on the results of IHC analysis, patients were divided into the high- and low-ZDHHC9-expression groups.

2.15 Western blotting

Whole-cell extracts were prepared by lysing cells in RIPA buffer supplemented with a protease inhibitor cocktail (Sigma-Aldrich, USA). The extracted protein was quantified using a BCA protein assay kit (Thermo Fisher Scientific, USA). A total of 20 μg of extracted protein was separated via SDS-PAGE, transferred onto polyvinylidene fluoride (PVDF) membranes and incubated with 5% skim milk in TBST for 1 h at room temperature to block non-specific interactions. Subsequently, the membranes were incubated with primary antibodies against ZDHHC9 (1:1000; Proteintech, 24046-1-AP), GAPDH (1:5000; Proteintech, 60004-1-Ig) and β-actin (1:5000; Santa Cruz Biotechnology, SC-47778) overnight at 4 ℃. The following day, the membranes were incubated with secondary antibodies (goat anti-rabbit or mouse antibodies, 1:5000; Abmart), and protein bands were visualised using an ECL system (Millipore, USA).

2.16 RT-PCR

Total mRNA was extracted using the TRIzol reagent. The quality and quantity of the extracted RNA were measured on an ND-1000 spectrophotometer (Nanodrop Technologies). cDNA was synthesised using the qScript cDNA synthesis kit (Quanta Biosciences) according to the manufacturer’s instructions. Quantitative PCR was performed using the SYBR Green Master Mix, and the amplified products were detected via agarose gel electrophoresis. β-actin was used as a loading control. The primers used for PCR are listed in Supplementary Table 2.

2.17 Cell viability assay

Cell viability was assessed using the 3-(4,5-Dimethylthiazol-2-yl)-2,5-Diphenyltetrazolium Bromide (MTT) assay, following the described methods [18].

2.18 Wound-healing and transwell assays

Wound healing and migration assays were conducted using the specified 4T1 cells, following an established protocol as previously described [19]. For the Transwell assay, a sixfold diluted Matrigel was added to the upper chamber of the Boyden cell to evaluate invasion. The cells that migrated to the bottom of the filter were stained using a three-step staining kit from Thermo Fisher Scientific. Subsequently, a microscope was utilized to count the migrated cells, and statistical analysis was performed.

2.19 Animals and tumor model

Immunologically sound 4–6-week old female C57BL/6 mice were used for in vivo analysis. All mice were housed under specific pathogen-free conditions. All experimental procedures involving animals were approved by the Ethics Committee of the Tianjin Medical University and Hospital Cancer Institute and were performed in accordance with the policies and procedures recommended by the NIH Guide for the Care and Use of Laboratory Animals.

A total of 1 × 106 tumor cells (4T1-scramble and 4T1-shZDHHC9) in phosphate-buffered saline (PBS) were subcutaneously injected into the lateral abdomen of C57BL/6 mice. The width and length of tumors were measured every 4 days using an electronic calliper, and tumor volume was calculated using the following formula: volume = length × width2/2. The tumor volumes of the mice in this experiment did not exceed the maximum size or burden required by the ethics committee. The maximum diameter of tumours in mice was set at 15 mm, and the maximum volume was limited to 1500 mm3.

2.20 Flow cytometry

Tumors were treated with 1 mg/mL collagenase, 2.5 U/mL hyaluronidase and 0.1 mg/mL DNase to produce single-cell suspensions. The abundance of tumor-infiltrating CD8 + T cells (CD3 + and CD8 +), CD4 + T cells (CD3 + and CD4 +), Treg cells (CD4 + , CD25 + and FOXP3 +), NK cells (NK1.1 +) and exhausted CD8 + T cells (PD-1 + and Tim3 +) was evaluated. Isotypic controls were used as negative controls. Data were analysed using the FCS Express 7 software. Relevant antibodies are listed in Supplementary Table 3.

2.21 Statistical analysis

The R (version 3.4.0.3) and SPSS Statistics (version 22) software were used for statistical analysis. The Coxph function was used to implement univariate and multivariate Cox regression analyses, and the glmnet package in R was used to perform Lasso–Cox regression analysis. Additionally, the log-rank test was performed using the survdiff function in the survival package, and time-dependent ROC analysis was performed using the timeROC package. The characteristic gene expression profiles were demonstrated on heatmaps created using the ggplot2 package. The rms package in R was used to establish and use a nomogram. Pathway enrichment analysis was performed using the clusterProfiler package in R. Each experiment was performed independently at least three times. Data were expressed as the mean ± squared difference (SD). Student’s t-test was used to compare the means of two groups. The median survival time was estimated using the Kaplan–Meier method, and differences in survival between groups were examined using the log-rank test. Spearman rank correlation analysis was used to examine the correlation between parameters. Two-way analysis of variance and post hoc analysis with repeated measures (time × tumor volume) were performed to test for differences in tumor growth in mice among different groups. Statistical tests were performed using the GraphPad Prism 8 software. All statistical tests were two-tailed, and p-values of < 0.05 were considered statistically significant. Significant values were marked with an asterisk and indicated as follows: ns, p ≥ 0.05; *p < 0.05; **p < 0.01; ***p < 0.001; ****p < 0.0001.