The proportion of peripheral DNT cells is increased in BC patients

The preoperative percentage of DNT cells in peripheral blood samples from BC patients and healthy controls was analyzed by flow cytometry. As shown in Fig. 1, the percentage of DNT cells was significantly higher in the BC patients compared to the controls (P < 0.05). There were no significant differences between the two groups in terms of demographic parameters like age and gender.

Fig. 1

Gating strategies and expression profiles of double-negative T cells (DNT cells). DNT cells express CD3 but lack CD4 and CD8. A FSC and SSC were used for the lymphocyte gating. B FSC-A and FSC-H were used to distinguish cell aggregates from single cells. C Non-viable cells were excluded through 7-AAD staining. D The B cells were excluded on the basis of CD19 and CD3 expression. E Representative flow cytometry dot plots of peripheral blood mononuclear cells (PBMCs) isolated from a patient with breast cancer (BC). F Representative dot plots of PBMCs isolated from a healthy control (HC). G Percentage of DNT cells in BC patients and HCs. P values were calculated by two-tailed unpaired Student’s t-test (**, P < 0.01)

Identification and functional characterization of DEGs

The DEGs in the DNT cells isolated from BC patients were screened using Smart-seq2. Using |log2FC|> 2 and P < 0.05 as the criteria, we obtained 289 DEGs, of which 137 were upregulated and 152 were downregulated in the BC group. The volcano map of the DEGs is shown in Fig. 2.

Fig. 2

Volcano plot showing the differentially expressed genes. The criteria for the DEGs were |log2FC|> 2 and P < 0.05. Red dots represent upregulated genes, turquoise dots represent downregulated genes, and grey dots indicate genes without significantly altered expression

To further explore the potential biological functions of these DEGs, we performed GO and KEGG enrichment analyses using the clusterProfiler package of R software. The categories of biological processes (BP), molecular functions (MF), and cellular components (CC) were included in the GO analysis. The significantly enriched BP terms for the DEGs included immunoglobulin mediated immune response, complement activation, classical pathway, humoral immune response mediated by circulating immunoglobulin, complement activation, and B cell receptor signaling pathway. In addition, the DEGs were significantly associated with CC terms including external side of plasma membrane, immunoglobulin complex, and postsynaptic membrane. Regarding MF terms, DEGs were enriched in antigen binding and immunoglobulin receptor binding. Finally, KEGG analysis revealed that the DEGs were mainly enriched in pathways related to protein digestion and absorption, hematopoietic cell lineage, B cell receptor signaling, ATP-binding cassette transporters, and complement and coagulation cascades (Fig. 3).

Fig. 3

Gene Ontology and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses of the differentially expressed genes (DEGs). A The top 10 enriched biological process (BP), cellular component (CC), and molecular function (MF) terms. B Significantly enriched KEGG pathways. The sizes of the dots represent the number of genes in each term

Identification of core genes related to the DNT cells

To investigate possible associations between the DEGs, we constructed a PPI network using the STRING website. As shown in Fig. 4A, the PPI network consisted of 183 nodes with 121 edges. Using the MCODE plug-in of Cytoscape, we filtered five modules and obtained two pivotal modules. One module included IFIT1, IFIT3, IFI27, IFI44L, RSAD2, and EGR1 (Fig. 4B), and the second module included C1QB, TMEM176A, C1QC, and TMEM176B (Fig. 4C). The top 10 hub genes were identified by CytoHubba (ranked in MCC), and included nine upregulated genes (IFIT1, RSAD2, IFI27, IFIT3, EGR1, IFI44L, C1QB, C1QC, NGFR) and one downregulated gene (VCAM1, Fig. 4D).

Fig. 4

Protein–protein interaction (PPI) network and modular analysis of differentially expressed genes (DEGs). A PPI network. B–C The top two modules of the DEGs according to MCODE. D Selection of hub genes using the CytoHubba plug-in of Cytoscape

Survival analysis of the identified core genes

We identified the potential core genes as those screened by both MCODE and CytoHubba. The prognostic values of these genes were estimated through survival analysis using the Kaplan–Meier method. As shown in Fig. 5, IFIT3, RSAD2, TMEM176B, C1QB and C1QC were significantly associated with survival rates, indicating their potential as prognostic indicators based on expression level.

Fig. 5

Prognostic significance of the 10 core genes in BC. A–J Kaplan–Meier curves showing the overall survival of BC patient subgroups based on the expression of IFIT1, IFII27, RSAD2, IFIT3, EGR1, IFI44L, TMEM176A, TMEM176B, C1QB, C1QC, and EGR1. The prognostic significance of each gene was determined on the basis of the hazard ratio (HR) and P value. P < 0.05 indicated statistical significance

Immune infiltration analysis

To further explore the role of DNT cells in tumor immunity and the potential mechanisms, we analyzed the correlation between the core genes and immune cells and immune checkpoints using the TIMER 2.0 and TISIDB databases. Interestingly, TMEM167B (r = 0.307, P = 3.90e-23), C1QB (r = 0.435, P = 3.20e-47), C1QC (r = 0.4, P = 2.14e-39), RSAD2 (r = 0.194, P = 7.58e-10) and IFIT3(r = 0.259, P = 9.46e-17) were positively correlated with the infiltration of CD8+ T cells. In addition, TMEM176B (r = 0.163, P = 2.36e-07), C1QB (r = 0.277, P = 6.36e-19), C1QC (r = 0.254, P = 3.80e-16), RSAD2(r = 0.09, P = 4.54e-03) and IFIT3 (r = 0.13, P = 3.94e-0.5; Fig. 6A–E) correlated significantly with the NK cells. The expression levels of three core genes were closely associated with immune checkpoints such as PDCD1, CD274, LAG3 and CTLA4 (Fig. 7A–E).

Fig. 6

Correlation analysis of gene expression and infiltration of CD8.+ T cells and NK cells in BC (TIMER database). A TMEM176B, B C1QB, C C1QC, D RSAD2, E IFIT3

Fig. 7

Correlation analysis of gene expression and immune checkpoints (PDCD1, CD274, LAG3 and CTLA4) in BC. A TMEM176B, B C1QB, C C1QC, D RSAD2, E IFIT3

DNA methylation analysis

We next examined the correlation between gene expression and promoter methylation using the UALCAN database to determine the impact of epigenetic alterations on breast tumorigenesis and development. The methylation levels of TMEM176B, RSAD2 and IFIT3 were higher in the BC samples compared to that in the normal samples (Fig. 8A, D–E; P < 0.05), whereas C1QB and C1QC showed lower promoter methylation in the BC samples (Fig. 8B–C; P < 0.05). Moreover, analysis of the MerthSuv database showed that C1QB hypomethylation in 3’UTR-open-sea-cg18763854, RSAD2 hypomethylation in TSS1500-open-sea-cg15346781, and C1QC hypermethylation in TSS1500-open-sea-cg17097874 located in the CpG islands were associated with a favorable prognosis in BC patients (Fig. 8F–J). In contrast, the methylation status of TMEM176B and IFIT3 were not significant associated with prognosis (P < 0.05).

Fig. 8

Association between promoter methylation and prognosis in BC. A–E Methylation level of the five prognostic genes in BC patients and healthy controls. F–J Kaplan–Meier curves showing the impact of promoter methylation (from the MethSurv database) on survival

Gene expression validation by RT-qPCR

To confirm the results of Smart seq-2 RNA-seq, the expression of TMEM176B, EGR1, C1QB, C1QC, RASD2 and IFIT3 were further analyzed by RT-qPCR using GAPDH as the internal control. Compared to the controls, TMEM176B was significantly downregulated, whereas EGR1, C1QB, C1QC, RASD2 and IFIT3 were upregulated in the BC samples (P < 0.05, Fig. 9).

Fig. 9

Verification of the expression levels of six genes. HC healthy control, BC breast cancer. RNA was extracted from peripheral blood DNT cells of BC (n = 5) and HC (n = 5) groups respectively. The gene expression was analyzed using the SYBR green RT-qPCR method. Data are presented as the mean ± SD. P values were calculated by two-tailed unpaired Student’s t-test (*, P < 0.05; **, P < 0.01)