CKI (Batch No: 20181034, total alkaloid concentration of 25 mg/mL) was provided by Zhendong Pharmaceutical Co., Ltd (China). The CKI was diluted ten-fold in ultrapure water, and 2 µl of the solution was used for further analysis. Five control compounds were obtained, including sophocarpine (Batch No: 20052711, purity ≥ 99.84%), oxymatrine (Batch No: 20041315, purity ≥ 98.76%), matrine (Batch No: 20200820, purity ≥ 98%), hesperidin (Batch No: 200621, purity ≥ 98%) and sophoridine (Batch No: 19113001, purity ≥ 96.97%). Matrine was purchased from Beijing North Weiye Metrology Institute Co., Ltd (Beijing, China), and the other four compounds were purchased from Beina Biotechnology Institute Co., Ltd (Beijing, China).

The CKI was separated by applying a Waters ACQUITY UPLC BEH C18 column (2.1 mm × 100 mm, 1.7 µm, Waters Corporation, Milford, MA, United States) at 35°C. Mobile phases comprised 0.1% aqueous formic acid and acetonitrile. A gradient elution with the flow rate of 0.2 mL/min was executed as follows: 6% acetonitrile at 0-2 min, 6-15% acetonitrile from 2 to 4 min, 15-25% acetonitrile from 4 to 8 min, 25-45% acetonitrile from 8 to 14 min, 45-60% acetonitrile from 14 to 16 min, 60% acetonitrile from 16 to 18 min, 60-6% acetonitrile from 18 to 19 min, and 6% acetonitrile from 19 to 22 min. The MS analysis was performed using the electrospray ionization (ESI) source in both positive and negative ion modes, and leucine enkephalin was utilized for mass accuracy correction. The temperatures of ion source and desolvation gas were set at 120°C and 350°C, respectively. The flow rates of the cone and the desolvation gas were 50 L/h and 600 L/h, respectively. The capillary voltages were set to 3.0 kV and 2.5 kV in positive and negative ion modes, respectively. The cone and extraction cone voltages were set to 40kV. MS/MS analysis was performed with a low collision energy of 4 eV and a high collision energy of 20-35 eV. The scan area was set at m/z 50-1200. Data acquisition and analysis were conducted with MassLynx™ v4.1 (Waters Co., Ltd) and UNIFI R Scientific Information System v1.7 (Waters Co., Ltd).

The RNA-seq dataset of CKI on breast cancer MDA-MB-231 cells was download from European Nucleotide Archive (35) with the accession number PRJNA517432 (36, 37), in which 12 samples at 48-hour (three untreated in batch 1, three untreated in batch 2, three CKI-treated in batch 1 and three CKI-treated in batch 2) were included in our study ( Supplementary Table 1 ). FastQC (version 0.11.9, Babraham Bioinformatics) was used to check the quality of raw reads before proceeding with downstream analysis. Trim_galore (version 0.6.6, Babraham Bioinformatics) was used to trim adaptors and low-quality sequences. STAR (version 2.7.7a) (38) was then applied to construct a reference genome index based on the DNA sequence and gene transfer format (GTF) files of the reference genome [GRCh38, Ensembl Release 103 (39)], and thereby the trimmed reads were further aligned to the above index by the STAR software and Binary Alignment/Map (BAM) format files were sorted by samtools (version 1.10) (40). Then, read counts data was generated after preparing reference sequences (reference genome: GRCh38, Ensembl Release 103) and calculating expression values by RSEM (v1.3.3) (41).

Differential expression analysis was performed by the DESeq2 (42), edgeR (43) and limma (44, 45) packages, respectively. The official pipelines of the three R packages that consider batch effects in RNA sequencing data in Bioconductor (46) were referenced in our study. The “filterByExpr” function in edgeR was used to screen genes with sufficiently large counts for a statistical analysis and scaling factors were calculated with the trimmed mean of M values (TMM) method (47) to convert raw library sizes into effective library sizes. The apeglm method (48) in the “lfcShrink” function was used to shrink log2 fold changes (FCs) when applying the DESeq2 method to perform differential expression analysis. The significance threshold for differential gene expression screening was set as adjusted P < 0.05 and |log₂FC| > log₂1.5, and the overlapping genes from the three methods were considered as DEGs. Normalized expression values of RNA-seq data were obtained with the voom algorithm (49), and batch effects were removed by the sva package (50) for further analysis.

Pathway analyses were performed on the 50 hallmark gene sets described in the Molecular Signatures Database (MSigDB, version 7.4) (51), and were then accomplished on biological process signatures in the Gene Ontology (GO) (52) and six kind of pathways comprising metabolism, genetic information processing, environmental information processing, cellular processes, organismal systems and human diseases deposited in the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (Last updated: May 1, 2021) (53). The gene sets with at least two genes found in the RNA-seq data were retained, and the GSVA algorithm was called from within the GSVA package (54) to calculate the enrichment score of each biological pathway in each sample with transcriptomic data. Subsequently, the empirical Bayesian approach within the limma package was applied to determine significantly changed pathways, with adjusted P < 0.05 as a significant cutoff criterion.

The METABRIC (55) cohort (a total of 1,980 patients, 100% female; including 320 patients with ER-, PR- and HER2- status) of breast cancer patients was included as the training set in our study. Non-primary breast cancer cases or patients without expression profiles were excluded, so 298 primary TNBC patients (100% female, average age = 56 ± 14 years) were reserved for further analysis ( Supplementary Table 2 ). The event was defined as died of disease and the patients who were living or died of other causes were censored. Furthermore, another cohort with 107 TNBC patients (56) was included as the test set to externally validate the performance of our survival model (GSE58812; 100% female, average age = 57 ± 13 years), and the metastasis-free survival (MFS) time of these cases was extracted ( Supplementary Table 3 ). The CEL format files of this microarray data were downloaded and were normalized with the Robust Multi-array Average (RMA) method in the affy package (57).

Feature gene sets of 28 subpopulations of tumor-infiltrating immune cells (58) were referenced, which include 28 immune cell types (13 innate and 15 adaptive immune cells), namely, activated dendritic cells (DCs), CD56bright natural killer (NK) cells, CD56dim natural killer (NK) cells, eosinophils, immature dendritic cells (DCs), macrophages, mast cells, myeloid-derived suppressor cells (MDSCs), monocytes, natural killer (NK) cells, Natural killer T (NKT) cells, neutrophils, plasmacytoid dendritic cells (pDCs), activated B cells, activated CD4⁺ T cells, activated CD8⁺ T cells, CD4⁺ central memory T (Tcm) cells, CD8⁺ central memory T (Tcm) cells, CD4⁺ effector memory T (Tem) cells, CD8⁺ effector memory T (Tem) cells, gamma delta T (γδ T) cells, immature B cells, memory B cells, regulatory T (Treg) cells, T follicular helper (Tfh) cells, type 1 T helper (Th1) cells, type 17 T helper (Th17) cells and type 2 T helper (Th2) cells. In the METABRIC cohort, the ssGSEA algorithm was called from within the GSVA package to quantify the relative infiltration of the 28 immune cell types in the tumor microenvironment (TME). The normalized enrichment score (NES) of these immune cell signatures calculated by ssGSEA was utilized to indicate the relative abundance of each immune cell in each TNBC sample. The prognostic value of each cell subset was evaluated with the univariate Cox proportional hazards model. ESTIMATE (59) was applied to evaluate the infiltration level of immune cells (immune score), the level of stromal cells (stromal score) and tumor purity for each patient.

The unsupervised clustering Pam method based on Euclidean and Ward’s linkage was called from within the ConsensuClusterPlus package (60) to determine the optimal number of stable immune-based TNBC clusters, and was repeated 100,000 times to ensure the stability of classification. The hierarchical clustering and k-means clustering were performed to confirm the robustness of the clustering. Principle component analysis (PCA) and t-Distributed Stochastic Neighbor Embedding (t-SNE) were further applied for dimensional reduction of TNBC patients based on tumor-infiltrating immune cells. Kaplan-Meier survival analysis with the log-rank test was implemented to compare the differences in prognosis among distinct immune subgroups. To identify TME-relevant genes, the patients in the METABRIC cohort were grouped into three different clusters based on immune-cell infiltration. The empirical Bayesian approach within the limma package was applied to determine DEGs among distinct immune phenotypes, and DEGs with adjusted P < 0.05 and |log₂FC| > log₂1.5 were considered as TME-related genes. The correlation between the relative expression value (Z-Score transformed) of each gene and DSS was evaluated by the univariate Cox regression analysis, and genes with P < 0.05 were considered as prognosis-related genes. The genes associated with both TME and prognosis were finally reserved for further analysis.

The CKI-perturbed genes were intersected with the TME- and prognosis-relevant genes, and subsequently the down-regulated genes with HR > 1 and the up-regulated genes with HR < 1 were reserved. Functional annotation for immune-related genes regulated by CKI was performed by the hypergeometric test. The formula is as follows:

where N is the number of all genes annotated by GO; n is the number of target genes in N; M is the number of all genes annotated in a specific GO term; m is the number of target genes that can be annotated in a specific GO term. After the P value was adjusted by the Benjamini-Hochberg correction, GO terms that met adjusted P < 0.05 were defined as significantly enriched GO terms. To obtain the most useful immune-related prognostic markers in the training set, the least absolute shrinkage and selection operator (LASSO) penalty within the glmnet package (61) was implemented to reduce dimensionality. Genes represented by an optimal value of the penalty parameter λ (the λ with minimum mean cross-validated error) that was determined by five-fold cross-validation constituted the immune-related gene signature in this study. Finally, these genes were fitted into a multivariate Cox regression model, and the linear combination of the expression value of the gene multiplied by its regression coefficients derived from the multivariate Cox regression model generated a prognostic risk score with the genes. The formula is as follows:

where t is the survival time, h₀ (t) is the baseline hazard, n is the number of genes, x_i is the expression value of the ith gene and β_i is the regression coefficient of the ith gene.

The patients were divided into either low- or high-risk group according to the optimal cut-off point determined via the maximally selected rank statistics method in the survminer package, and meanwhile the minimal proportion of observations per group was set to 30% to prevent the problem of too few patients in a certain group. The survival curves were generated by the Kaplan-Meier method, and the log-rank test was used to assess the differences between the low- and high-risk groups. Survival predictive accuracy of the prognostic model was estimated through receiver operating characteristic curve (ROC) and Harrell’s concordance index (C-index) analyses. The predictive efficacy of the gene signature for patients’ prognosis was further validated in the test set mentioned above. The univariate and multivariate Cox proportional hazards models were utilized to analyze whether the risk scores had a prognostic significance. Age, tumor size, the number of positive lymph nodes and PAM50 subtypes were firstly evaluated in the univariate Cox proportional hazards model, and all statistically significant variables were then added as covariates in the multivariate Cox proportional hazards model. Further validation on the prognostic value of risk scores was performed by comparing the predictive ability between two models: one comprising tumor size and the number of positive lymph nodes as covariates; and the other including risk scores as the third covariate. The time-dependent C-index and the area under curve (AUC) of time-dependent ROC curve were set as the indicators of prognostic efficacy.

For two-group comparisons, the Shapiro-Wilk test was applied to assess the assumption of normal distribution, and statistical significance of differences between non-normally distributed variables was estimated with the Wilcoxon test. For comparisons among more than two groups, the Kruskal-Wallis test was utilized as a non-parametric method. Correlation coefficients (ρ) were generated by Spearman’s and distance correlation analyses. Survival analysis was performed via the Kaplan-Meier method, and differences in survival distributions were evaluated using the log-rank test. The “surv-cutpoint” function in the survminer package, which repeatedly tests all potential cut points to find the maximum rank statistics, was applied to determine the best cut-off level for each prognostic marker. Univariate and multivariate analyses were conducted employing Cox proportional hazard models. All the tests were two sided, and statistical significance was set at P < 0.05, unless otherwise stated. The Benjamini-Hochberg correction was applied in multiple tests to reduce false positive rates. All statistical analyses in this study were performed with the R software (version 4.0.3, http://www.R-project.org).

In this study, we identified a total of 23 chemical components from CKI by using UPLC-Q-TOF-MS analysis ( Figure 1 and Table 1 ). Consistent with previous findings (37, 62), we also identified the six alkaloids that have been considered as the major active ingredients of CKI, namely, matrine, oxymatrine, N-methylcytisine, sophoridine, sophocarpine and oxysophocarpine.

The results of RNA-seq data analysis for the CKI-perturbed cell line samples were shown in Figure 2 and Supplementary Tables 4–6 . The PCA showed the difference in the sample distribution before or after removing batch effects ( Figures 2A, B ). The hierarchical clustering presented a strong intra-group correlation and a relatively low inter-group correlation ( Figure 2C ). We found 1925 consistent up-regulated genes and 1767 consistent down-regulated genes between CKI-treated and control groups ( Figure 2D ), and Pearson’s correlation analysis exhibited high correlation of log₂FC calculated by three methods. The heatmap and volcano plots for the DEGs were displayed in Figures 2F, G .

Analysis of hallmark pathway gene signatures highlighted that cell cycle-related pathways like MYC targets, G2M checkpoint and E2F targets were significantly down-regulated in the CKI treatment group compared with the control group ( Figures 3A, B ). According to the GSVA results of KEGG pathways ( Figure 3C ), 113 biological pathways (including 53 activated and 60 suppressed) were significantly changed between CKI-treated and control groups, which exhibited a comprehensive influence of CKI on metabolism, cellular processes, genetic and environmental information processing, and organismal systems. Metabolic pathway analysis showed CKI decreased pathways involved in nucleotide metabolism (purine and pyrimidine metabolism) and lipid metabolism (fatty acid biosynthesis) ( Figure 3D ). Especially, the KEGG pathways associated with immune system such as natural killer cell mediated cytotoxicity, T cell receptor signaling pathway, Th1 and Th2 cell differentiation, Th17 cell differentiation and B cell receptor signaling pathway were activated in the CKI-treated group versus the control group, which suggested its potential effects on immunoregulation ( Figures 4A, B ). Analysis of GO gene signatures highlighted that most of biological processes that promote immune response were activated after CKI treatment ( Figures 4C, D ). Compared with the samples in the control group, most biological processes that improve T cell proliferation, differentiation and T cell receptor signaling pathway or induced cell death of T cells were activated, whereas those that inhibit T cell activation, proliferation, differentiation and T cell receptor signaling pathway were inactivated in the CKI-perturbed samples ( Figure 4C ). Meanwhile, the activation and proliferation of NK, NKT and B cells were up-regulated after CKI treatment. Notably, for the Th1 cells, which execute anti-tumor immunity functions, multiple processes that promote their activities such as positive regulation of T-helper 1 cell cytokine production and positive regulation of T-helper 1 type immune response were positively regulated in the CKI-treated group. However, for the Th2 cells, which executing pro-tumor, immune suppressive functions, the cell cytokine production of Th2 cells was down-regulated in the CKI-treated group ( Figure 4D ).

Spearman’s correlation analysis showed that the abundance of cells executing anti-tumor reactivity was positively associated with the abundance of cells delivering pro-tumor suppression ( Figures 5A , 8A ). The universal landscape of immune cell interaction in TME was visualized in Figure 5B , which demonstrated that the relative abundance of most infiltrating immune cell populations was positively correlated with each other, stromal scores and immune scores and negatively correlated with tumor purity. Prognostic significance of each microenvironment cell in TNBC was shown in Figures 5C, D . Based on the enrichment scores of 28 gene signatures for the 298 TNBC samples, we performed unsupervised clustering to classify these patients into three independent subtypes (we named them as immunity low, immunity medium and immunity high) ( Figures 6A–F , 8B ). The three clusters witnessed a remarkable difference in the relative infiltration of immune cell populations ( Figures 7 , 6G , 8C, D ), stromal scores, immune scores and tumor purity ( Figures 6H , 8E ). Considering the vital role of TME in prognosis, we investigated the clinical relevance of the immune clusters. Kaplan-Meier analysis exhibited significant associations between immune clusters and DSS (log-rank test, P = 0.012), and the immunity high or medium cluster had a longer survival time than the immunity low cluster ( Figure 6I ). The multivariate Cox proportional hazards model also disclosed that the immune clusters independently predicted a better DSS in TNBC (Immunity medium: HR, 0.55; 95% confidence interval, 0.36-0.82; P = 0.004; Immunity high: HR, 0.47; 95% confidence interval, 0.26-0.84; P = 0.011; Table 2 ). We also analyzed the three crucial immune checkpoints PD-L1, PD1 and CTLA4 in each immunity subtype. The immunity high cluster was featured by a significantly higher PD-L1/PD1/CTLA4 expression levels while the immunity low cluster with a lower PD-L1/PD1/CTLA4 expression levels ( Figure 9 ). After conducting differential expression analysis, we found 1591 DEGs between the immunity high and immunity low groups, and 313 between the immunity medium and immunity low groups ( Supplementary Table 7 ). We aggregated these genes and finally identified 1593 unique genes as immune-related genes in TNBC.

For the 1593 immune-related genes in TNBC, univariate Cox analysis showed 304 of them were correlated with DSS time (P < 0.05), in which two genes (DYNC2I2 and MARVELD2) with HR > 1 can be down-regulated by CKI and 26 genes (PRKCH, TNFRSF1B, PAG1, LAT2, ARHGAP4, SEL1L3, RASSF2, LPXN, IL23A, ALDH2, ST8SIA4, HSD11B1, ARHGAP9, STX11, SLCO2B1, STAT4, FERMT3, GBP2, CTSW, CD7, SLCO3A1, SEMA4D, SERPINA1, MICA, SIPA1 and RASSF5) with HR < 1 can be up-regulated by CKI ( Figure 10A ). In the present study, these 28 genes were considered as key TME- and prognosis-related genes regulated by CKI. Clearly, the 26 genes were positively associated with the relative abundance of most infiltrating immune cell populations, while the other two genes showed negative associations ( Figure 10B ). GO functional annotation presented that the 28 CKI-perturbed immune-related genes were significantly enriched in immune system process, immune response, regulation of immune response, adaptive immune response, T cell activation and T cell mediated immunity ( Figure 10C ). The expression of the 26 genes that up-regulated by CKI were positively correlated with each other, PD-L1, PD-1 and CTLA4, while the expression of MARVELD2 and WDR34 that down-regulated by CKI were negatively correlated with the other 26 genes ( Figure 10D ). These genes were enrolled to construct a prognostic signature.

After performing LASSO and multivariate Cox regression analysis, a five-gene signature correlated with prognosis was constructed, in which two genes (DYNC2I2 and MARVELD2) with HR > 1 were considered as risky genes and three genes (RASSF2, FERMT3 and RASSF5) with HR < 1 were considered as protective genes ( Figures 11 – 13 ). The regression coefficient for each gene was also generated, and the survival risk score was computed as follows: risk score = (0.1202 × expression level of DYNC2I2) + (0.2372 × expression level of MARVELD2) + (-0.3560 × expression level of RASSF2) + (0.4986 × expression level of FERMT3) + (-0.5813 × expression level of RASSF5). For the training set, the 206 patients with risk scores higher than the cutoff value (1.464) were classified into the high-risk group, while the rest 92 patients were classified into the low-risk group ( Figure 11A ). The Kaplan-Meier survival analysis exhibited that patients in the high-risk group had shorter survival time in comparison with patients in the low-risk group (Log-rank test P < 0.001), suggesting the expression levels of the five genes could effectively discriminate the survival risks of these TNBC patients ( Figure 11C ). The C-index was 0.646 in the training set, and the AUC of the ROC curve was 0.70, 0.68, 0.68, 0.68, 0.68 and 0.66 for 1-year, 2-year, 3-year, 5-year, 7-year and 10-year DSS, respectively, confirming a good predictive efficacy of the prognostic gene signature ( Figure 11E ). For the test set, the 83 patients with risk scores higher than the cutoff value (1.464) were classified into the high-risk group, while the rest 24 patients were included into the low-risk group ( Figure 11B ). The Kaplan-Meier survival analysis showed that patients in the high-risk group had shorter survival time in comparison with patients in the low-risk group (Log-rank test P = 0.001), suggesting the expression levels of the five genes could effectively discriminate the survival risks of these TNBC patients ( Figure 11D ). The C-index was 0.696 in the test set, and the AUC of the ROC curve was 0.72, 0.73, 0.72, 0.71 and 0.72 for 2-year, 3-year, 5-year, 7-year and 10-year MFS, respectively, confirming a good predictive efficacy of the prognostic gene signature ( Figure 11F ). The multivariate Cox proportional hazards model also revealed that the risk score independently predicted a worse DSS in TNBC (HR, 1.51; 95% confidence interval, 1.27-1.81; P < 0.001; Table 3 ). The time-dependent C-index and the time-dependent ROC curve demonstrated that the addition of risk scores into the Cox proportional hazards model enhanced the prognostic efficacy ( Figures 11G, H ). Higher expression of DYNC2I2 and MARVELD2 was correlated with worse prognosis, whereas higer expression of RASSF2, FERMT3 and RASSF5 was correlated with better prognosis ( Figures 12A, B ). As shown in Figures 12C, D , DYNC2I2 and MARVELD2 were up-regulated in the high-risk group versus the low-risk group, and RASSF2, FERMT3 and RASSF5 were down-regulated in the high-risk group versus the low-risk group. The expression of RASSF2, FERMT3 and RASSF5 that up-regulated by CKI were positively correlated with each other, PD-L1, PD-1 and CTLA4, whereas the expression of MARVELD2 and WDR34 that down-regulated by CKI were negatively correlated with RASSF2, FERMT3 and RASSF5 ( Figures 12E, F ). We also analyzed the expression of PD-L1, PD-1 and CTLA4 in the high- and low-risk groups. The low-risk group was characterized by a significantly higher PD-L1/PD1/CTLA4 expression levels, while the high-risk group with a lower PD-L1/PD1/CTLA4 expression levels ( Figures 12G, H ).