Abstract

   Breast cancer is a leading malignancy in women, with mortality
   disparities between developed and underdeveloped regions. Accumulating
   evidence suggests that the competitive endogenous RNA (ceRNA)
   regulatory networks play paramount roles in various human cancers.
   However, the complexity and behavior characteristics of the ceRNA
   network in breast cancer progression have not been fully elucidated.
   The expression profiles of three RNAs (long non-coding RNAs [lncRNAs],
   microRNAs [miRNAs], and mRNAs) were extracted from breast cancer and
   adjacent samples were sourced from the TCGA database. The SLC26A4-AS1-
   hsa-miR-19a-3p-NTRK2 ceRNA network related to the prognosis of breast
   cancer was obtained by performing bioinformatics analysis. Importantly,
   we identified the SLC26A4-AS1/NTRK2 axis within the ceRNA network
   through correlation analysis and found it to be a potential prognostic
   model in clinical outcomes based on Cox regression analysis. Moreover,
   methylation analysis suggests that the aberrant downregulation of the
   SLC26A4-AS1/NTRK2 axis might be attributed to hypermethylation at
   specific sites. Immune infiltration analysis indicates that the
   SLC26A4-AS1/NTRK2 axis may have implications for the alteration of the
   tumor immune microenvironment and the emergence and progression of
   immune evasion in breast cancer. Finally, we validated the expression
   of SLC26A4-AS1-hsa-miR-19a-3p-NTRK2 in breast cancer cell lines. In
   summary, the present study posits that the SLC26A4-AS1/NTRK2 axis,
   based on the ceRNA network, could be a novel and significant prognostic
   factor associated with breast cancer diagnosis and outcomes.

Supplementary Information

   The online version contains supplementary material available at
   10.1007/s12672-025-02080-7.

   Keywords: Breast cancer; LncRNA; MiRNA; MRNA; CeRNA network;
   Prognosis,; Methylation; Immune infiltration

Introduction

   Breast cancer (BC) stands as the predominant malignant tumor in women,
   securing a spot within the top three most diagnosed cancers globally.
   According to GLOBOCAN 2022, BC remains the most frequently diagnosed
   cancer in women, with an estimated 2.3 million new cases and 685,000
   deaths globally in 2022, accounting for 11.6% of all new cancer cases
   and 6.9% of all cancer-related deaths. While developed regions have
   observed a notable decline in breast cancer-related deaths, primarily
   due to early detection and effective systemic treatments, breast cancer
   remains a leading cause of cancer-related deaths in underdeveloped
   regions, where late-stage diagnoses and limited access to treatment are
   significant challenges [[34]1]. Therefore, the persistent advancement
   in breast cancer diagnosis and therapy becomes paramount to curbing its
   mortality. Thus, continuous updates in breast cancer diagnosis and
   treatment are essential for reducing mortality rates in breast cancer.

   The public availability of comprehensive cancer genomic datasets,
   notably The Cancer Genome Atlas (TCGA), offers an abundance of clinical
   and molecular insights for diverse cancer types, serving as a pivotal
   resource in pinpointing potential tumor biomarkers. The burgeoning use
   of high-throughput RNA sequencing (RNA-seq) has catalyzed the creation
   of enhanced tools, recognized for their non-invasiveness,
   reproducibility, sensitivity, and precision, benefiting early
   screenings, diagnoses, patient evaluations, and continual monitoring
   [[35]2].

   Interestingly, while roughly 75% of the human genome can generate
   functional transcripts, a mere 2% encode proteins, directing the
   spotlight to non-coding RNAs (ncRNAs) which encompass lncRNAs, miRNAs,
   circular RNAs, and pseudogenes [[36]3]. Despite mRNA research being
   relatively mature, investigations on lncRNAs, miRNAs, and
   pseudotranscripts remain nascent, with vast uncharted territories
   awaiting exploration [[37]4]. Defined as non-coding RNA variants
   spanning over 200 nucleotides with little to no protein-coding
   capabilities, lncRNAs engage in modulating numerous cellular activities
   through a myriad of channels, encompassing transcriptional governance
   via interactions with DNA [[38]5], RNA [[39]6], and proteins [[40]7],
   and more. Their influence spans epigenetic modifications, protein/RNA
   stability, translation, post-translational modifications, and even
   direct interactions with signaling receptors [[41]8]. The pivotal roles
   of lncRNAs in cancer biology, especially their frequent anomalies in
   cancers [[42]9], including BC, underline their clinical significance.
   Some lncRNAs have been associated with cancer relapses, metastasis
   [[43]10], and post-radiotherapy complications [[44]11]. Therefore,
   investigating the role of lncRNAs in the prognosis and diagnosis of BC
   holds significant clinical importance.

   On another note, microRNAs (miRNAs) are diminutive, intrinsic, and
   highly conserved ncRNAs, orchestrating nearly 30% of gene expressions
   post-transcriptionally. Typically, they hinder mRNA translation and
   stability, thereby governing genes pivotal in cellular functions like
   inflammation, cell cycle, stress reactions, differentiation, apoptosis,
   and migration [[45]12]. Moreover, the intricate play of miRNA-mediated
   regulation leans heavily on an assortment of RNA-binding proteins,
   granting miRNAs the capacity to influence tumor onset and development.

   Within the cytoplasmic realm, Salmena et al. [[46]13] proposed the
   intriguing competing endogenous RNA (ceRNA) hypothesis, which portrays
   a dynamic RNA transcript communication network. By competing for shared
   miRNA sequences, these transcripts adjust one another's expressions.
   This theory underscores pivotal discoveries concerning the interplay
   between coding and non-coding RNAs mediated by miRNA. Earlier research
   illustrated that lincRNA-RoR (long intergenic noncoding RNA regulator
   of reprogramming) could bolster the invasive prowess of triple-negative
   breast cancer via the miR-145/ARF6 route by destabilizing E-cadherin
   positioning [[47]14], a key epithelial adhesion molecule that maintains
   intercellular junctions and tissue integrity. The loss of E-cadherin
   destabilizes cell–cell adhesion, promoting cytoskeletal rearrangement
   and enhancing cellular motility, processes central to
   epithelial-mesenchymal transition (EMT). EMT not only drives breast
   tumor invasion and metastasis but also contributes to chemoresistance
   and radio-resistance by enabling tumor cells to acquire
   mesenchymal-like traits and evade therapeutic interventions,
   underscoring the therapeutic potential of targeting EMT-regulating
   factors like lncRNAs and miRNAs [[48]15]. Furthermore, the PI3K-Akt
   signaling pathway, a key regulator of proliferation, survival, and
   migration, is frequently dysregulated in breast cancer and has been
   associated with resistance to various therapies, making it a crucial
   target for potential interventions [[49]16]. In addition, Chou et al.
   [[50]17] highlighted that the lncRNA MALAT1 escalates breast cancer
   cell migration and intrusion by competing for miR-1, catalyzing the
   surge of cell division control protein 42 (CDC42). Contrarily, Kim et
   al. [[51]18] discerned an inverse relation between MALAT1
   concentrations and breast cancer advancement, with MALAT1 acting more
   as a metastasis suppressor. It acts as a suppressor of tumor metastasis
   rather than a promoter. Therefore, unraveling the complex interactions
   of specific lncRNA-miRNA-mRNA axes through high level bioinformatics
   and network analysis is crucial for deciphering these pathways and
   identifying specific targets for breast cancer diagnosis and treatment.

   In this study, we hypothesize that specific ceRNA axes, such as
   SLC26A4-AS1-hsa-miR-19a-3p-NTRK2, play critical roles in the
   progression of breast cancer by influencing gene expression through
   epigenetic modifications and modulating the tumor immune
   microenvironment. To test this hypothesis, we constructed a
   tri-regulatory lncRNA-miRNA-mRNA ceRNA network specific to BC and
   identified SLC26A4-AS1-hsa-miR-19a-3p-NTRK2 as a key axis associated
   with breast cancer prognosis. We further explored the functional
   mechanisms of this axis through methylation assessments, immune
   infiltration analysis, and pathway enrichment studies, integrating
   bioinformatics approaches and experimental validation using the GEO
   database and breast cancer cell lines. Our findings reveal that the
   ceRNA-based SLC26A4-AS1/NTRK2 axis may serve as a novel and significant
   prognostic factor in breast cancer, providing potential therapeutic
   targets for future research.

Materials and methods

Data preparation and preprocessing

   We obtained mRNA sequencing (mRNA-seq) data for 1113 cases of BC and
   113 adjacent non-cancerous samples, as well as miRNA sequencing
   (miRNA-seq) data for 1103 cases of BC and 104 adjacent non-cancerous
   samples, from the TCGA database ([52]https://portal.gdc.cancer.gov/).
   Survival information for the TCGA-BC dataset was sourced from UCSC Xena
   ([53]http://xena.ucsc.edu/). We divided the publicly available data
   into cancerous and adjacent non-cancerous groups. We conducted
   differential analysis on the raw count matrix of the selected public
   data using the DESeq2 package, following the standard procedure.
   Furthermore, we applied the Variance Stabilizing Transformations (VST)
   method provided by the DESeq2 package to normalize the raw count
   matrix. For protein expression analysis, we utilized the UALCAN portal
   ([54]http://ualcan.path.uab.edu/analysis-prot.html). In addition, we
   obtained gene expression profiles from the Gene Expression Omnibus
   (GEO, [55]https://www.ncbi.nlm.nih.gov/GEO/) database, specifically
   [56]GSE45827 (Platform: [57]GPL570) and [58]GSE12790 (Platform:
   [59]GPL570), to further validate our findings.

Screening of differentially expressed genes (DEGs)

   Differential expression lncRNAs (DElncRNAs) and differential expression
   mRNAs (DEmRNAs) were identified with thresholds of |logFC|> 1 and
   p.adj < 0.05. The differential expression miRNAs (DEmiRNAs was
   identified with thresholds of |logFC|> 0.5 and p.adj < 0.05. Volcano
   plots for visualizing DEGs were generated using the R package
   ‘‘EnhancedVolcano’’. Heatmap clustering of DEGs was constructed using
   the R package ‘‘pheatmap’’.

Establishment of the ceRNA network in BC

   In accordance with the hypothesis that lncRNAs act as natural sponges
   in the cytoplasm, competitively regulating mRNA expression indirectly
   through miRNAs, the construction steps for the ceRNA network are as
   follows: (1) Utilizing miRNet ([60]https://www.mirnet.ca) to predict
   potential miRNAs targeted by DElncRNAs and lncRNA-miRNA interactions;
   (2) Using TargetScan ([61]http://www.targetscan.org) and miRDB
   ([62]http://www.mirdb.org) to predict target genes by overlapping the
   predicted miRNAs with DEmiRNAs; (3) Employing the R package
   ‘‘VennDiagram’’ to facilitate the intersection between target genes and
   DEmRNAs. selecting the resulting mRNA intersection for subsequent
   analysis; (4) Integrating the lncRNA-miRNA pairs with the miRNA-mRNA
   pairs to construct the lncRNA-miRNA-mRNA triple-regulation network. The
   generated network was visualized using Cytoscape software (version
   3.9.1).

   LncRNA sequences were obtained from LNCipedia
   ([63]https://lncipedia.org), and lncRNA cellular localization was
   determined using the LncLocator database
   ([64]http://www.csbio.sjtu.edu.cn/bioinf/lncLocator/).

Functional enrichment analysis

   We conducted functional enrichment analysis on DEmRNAs in the
   lncRNA-miRNA-mRNA triple-regulation network using the R software
   package ‘‘Goplot’’ and ‘‘clusterProfiler’’. Significance thresholds for
   GO terms and KEGG pathways were set at p.adj < 0.05. Additionally, we
   obtained a set of 20 NTRK2-binding proteins based on the STRING
   database ([65]https://cn.string-db.org) and supported by experimental
   evidence. We analyzed the correlation among these genes pairwise in BC,
   then visualized the correlation results using the ‘‘circlize’’ package
   and performed GO enrichment analysis (Molecular fuction; MF, Cellular
   component; CC, Biological process; BP) and KEGG pathway analysis.
   Visualization of the results was carried out using the R package
   ‘‘ggplot2’’.

Survival analysis and construction of survival-related ceRNA in BC

   Survival status and time in BC patients were sourced from TCGA clinical
   dataset, filtered for female gender and no history of neoadjuvant
   treatment. Overall survival (OS) was defined as the time from diagnosis
   to death from any cause. Disease-specific survival (DSS) was defined as
   the time from diagnosis to death specifically caused by breast cancer.
   Progression-free interval (PFI) was defined as the time from diagnosis
   to disease progression or recurrence. Patients without these events at
   the end of follow-up were censored. We conducted batch survival
   regression analysis on the ceRNA network using the Cox proportional
   hazards model in the “survival” package in R to candidate lncRNAs,
   miRNAs, and mRNAs associated with OS, DSS, and PFI. Kaplan–Meier
   survival curves were generated to visualize survival trends, and hazard
   ratios (HRs) with 95% confidence intervals (CIs) were calculated to
   assess the impact of gene expression on survival outcomes.

   Based on the expression levels of candidate genes, we performed
   time-dependent receiver operating characteristic curve (ROC) analysis
   using timeROC package to assess the predictive capability of
   biomarkers, calculating the area under the curve (AUC). The results
   were visualized using ‘‘ggplot2’’. When the value of a molecule (risk
   factor) trends towards promoting event occurrence, the AUC for that
   molecule is > 0.5. Larger areas (AUC value closer to 1) indicate better
   prognostic efficacy. Conversely, when a molecule (protective factor)
   trends in the opposite direction of event occurrence, the AUC for that
   molecule is < 0.5. Smaller areas (AUC values closer to 0) indicate
   better prognostic efficacy.

Methylation and expression analysis of NTRK2

   DNA methyltransferases (DNMT1, DNMT3A and DNMT3B) are critical players
   in epigenetic regulation, with roles ranging from the maintenance of
   existing DNA methylation patterns to the establishment of new ones.
   Their impact on gene expression and epigenetic memory makes them
   important subjects of study in the field of molecular biology and
   genetics. First, we investigated the expression levels of three DNA
   methyltransferases in NTRK2-high and NTRK2-low groups using the TCGA
   database. Subsequently, we assessed the methylation levels of NTRK2
   between BC and adjacent normal tissues using the DiseaseMeth 3.0
   version ([66]http://bio-bigdata.hrbmu.edu.cn/diseasemeth/) and UALCAN.
   Furthermore, we analyzed the methylation levels of NTRK2 CpG sites in
   BC using MethSurv ([67]https://biit.cs.ut.ee/methsurv/) and examined
   the correlation between these CpG sites and the survival of breast
   cancer patients. Finally, we employed SMART
   ([68]http://www.bioinfo-zs.com/smartapp/) to investigate the
   methylation values of different CpG sites of NTRK2 in BC and normal
   tissues and analyzed the relationship between these CpG sites and NTRK2
   expression.

Immune infiltration levels and expression analysis of NTRK2

   We employed the ssGSEA algorithm provided in R package GSVA to
   calculate the immune infiltration status in BC using markers for 24
   immune cell types as provided by Bindea et al. [[69]19]. Subsequently,
   we analyzed the correlation between NTRK2 and these immune cells and
   visualized using ‘‘ggplot2’’. Additionally, based on the median
   expression level of NTRK2 (4.435, derived from the log2-transformed
   expression values) in BC samples, we stratified the samples into two
   groups: high-NTRK2 expression group (≥ 4.435) and low-NTRK2 expression
   group (< 4.435). This stratification allowed us to analyze the
   differences in immune infiltration between these two groups.
   Furthermore, we used the TIDE ([70]http://tide.dfci.harvard.edu/query/)
   to evaluate the relationship between NTRK2 expression and CTL and
   assess the relationship between NTRK2 and T-cell exclusion score,
   T-cell dysfuction score, response to immune checkpoint blockade (ICB)
   therapy, and gene knockout phenotype in CRISPR screens.

Cell culture

   Human MCF 10A (American Type Culture Collection, USA) and MDA-MB-468
   cells (Cell Bank of the Chinese Academy of Sciences, China) were
   propagated in Dulbecco’s modified Eagle’s medium (DMEM; Corning)
   supplemented with 10% fetal bovine serum. Cells were cultured in
   standard T-25 or T-75 culture flasks (Corning) and maintained in an
   incubator at 37 °C with 5% CO₂ and saturated humidity. The medium was
   replaced every 2–3 days, and cells were subcultured at 70–80%
   confluence. Cell morphology and viability were monitored under an
   inverted phase-contrast microscope before each experiment, and cells
   were regularly tested for Mycoplasma contamination using a Mycoplasma
   Detection Kit (Biyuntian, C0303S).

Quantitative reverse transcription-polymerase chain reaction (qRT-PCR)

   Cells were washed twice with PBS, collected by centrifugation, and
   lysed with Trizol Reagent (Shanghai Pufei Biotech Co., Ltd., Shanghai,
   China) to extract total RNA, including the miRNA fraction, following
   the manufacturer's instructions. RNA purity and concentration were
   determined using a NanoDrop 2000 spectrophotometer (Thermo Fisher
   Scientific), and integrity was verified by agarose gel electrophoresis.

   For reverse transcription, 2 μg of total RNA was used with M-MLV
   reverse transcriptase (Promega) at 43.5 °C for 60 min, followed by
   enzyme inactivation at 73.5 °C for 3 min. qPCR reactions were carried
   out using LightCycler 480 SYBR Green (Roche) on a LightCycler 480 II
   (Roche). The reaction mix (12 μL) included 6.0 μL SYBR Green premix,
   0.5 μL each of forward and reverse primers (5 μM), 1.0 μL cDNA, and
   4.0 μL RNase-free water.

   qPCR cycling conditions were as follows: initial denaturation at 95 °C
   for 5 min, followed by 40 cycles of denaturation at 95 °C for 15 s,
   annealing at 60 °C for 20 s, and extension at 72 °C for 20 s. Melt
   curve analysis confirmed amplification specificity.

   Primers for qPCR were purchased from Ribobio and Shanghai Sangon
   Biotech. The specific primers for SLC26A4-AS1 were
   5′-TTTCTGCTCCTTGTTGCTT-3′ (forward) and 5′-TGCCTCCTGTTGTTACTCA-3′
   (reverse); for hsa-miR-19a-3p, they were ssD809230917 (forward;
   Ribobio), SSD089261711 (reverse; Ribobio), and ssD809230225 (RT;
   Ribobio); for NTRK2, they were 5′-TCACCTTGACTTGTCTGAACTG-3′ (forward)
   and 5′-CTGGACTGGATTTAGCCTCTT-3' (reverse). The primers for ACTB were
   5′-GCGTGACATTAAGGAGAAGC-3′ (forward) and 5′-CCACGTCACACTTCATGATGG-3′
   (reverse); for U6, they were ssD0904071008 (forward; Ribobio),
   ssD0904071006 (reverse; Ribobio), and ssD0904071007 (RT; Ribobio). The
   ΔΔCt method was used for data analysis, with ACTB and U6 serving as
   internal controls.

Western blotting

   The cells were lysed in RIPA buffer (Biyuntian, P0013), and total
   proteins were quantified using the BCA Protein Assay Kit (Biyuntian,
   P0010S). Equal amounts of protein (20–30 µg) were separated on a 10%
   denaturing polyacrylamide gel and transferred onto a PVDF membrane
   (Millipore) at 4 °C, 200 mA for 120 min. Membranes were blocked with 5%
   non-fat milk in TBST for 1 h at room temperature and incubated
   overnight at 4 °C with primary antibodies: anti-NTRK2 (Proteintech,
   13,129-1-AP, rabbit, 1:3000) and anti-β-actin (Santa Cruz, sc-69879,
   mouse, 1:13,000). After washing, membranes were incubated with
   HRP-conjugated secondary antibodies (CST, #7074 for anti-rabbit IgG;
   #7076 for anti-mouse IgG, 1:8000) for 1.5 h at room temperature.
   Protein bands were visualized using the LumiGLO® Reagent and Peroxide
   Kit (CST, #7003), and images were captured on X-ray films. NTRK2
   (90–95 kDa) and β-actin (43 kDa) served as the target and loading
   control, respectively.

Statistical analysis

   In this study, all data analysis and visualization were performed using
   GraphPad Prism software (version 9.0.0) and R software (version 4.2.1;
   [71]https://www.r-project.org/) and the corresponding packages. A
   p-value of less than 0.05 was considered statistically significant.

Results

   Identification of differentially expressed genes (DEGs) between BC and
   adjacent normal samples.

   A total of 1113 BC samples and 113 adjacent normal samples for
   identifying DElncRNAs and DEmRNAs, and 1103 BC samples and 104 adjacent
   normal samples for identifying DEmiRNAs. In BC samples, we identified
   4177 DElncRNAs, of which 64.3% (2686) were upregulated and 35.7% (1491)
   were downregulated. Similarly, 5065 DEmRNAs were identified, with 62.1%
   (3146) upregulated and 37.9% (1919) downregulated. For DEmiRNAs, we
   identified 118 differentially expressed miRNAs, including 62.7% (74)
   upregulated and 37.3% (44) downregulated, using stringent criteria
   (|log2 FC|> 0.5 and p.adj < 0.01). Volcano plots and heatmaps provided
   a visual representation of the distribution of DEGs between BC and
   adjacent normal samples, as well as the expression of 20 significant
   variable genes (Fig. [72]1).

Fig. 1.

   [73]Fig. 1
   [74]Open in a new tab

   Visualization of Differentially Expressed RNAs and Corresponding
   Heatmaps. A The volcano plot highlights the DElncRNAs. Accompanying the
   plot, a heatmap is presented, showcasing the expression patterns of the
   top 20 most variable DElncRNAs. B Similarly, a volcano plot for the
   DEmiRNAs is displayed. Adjacent to it, a heatmap elucidates the
   expression profiles of the top 20 most variable DEmiRNAs. C A volcano
   plot captures the DEmRNAs, with an adjoining heatmap revealing the
   expression trends of the top 20 most variable DEmRNAs

Construction of the lncRNA-miRNA-mRNA triple regulatory network

   DelncRNAs were subjected to miRNet database, which predicted 631
   potential miRNAs targeting them. Among these, 14.1% (89) overlapped
   with the identified DEmiRNAs. These 89 miRNAs were further used to
   predict mRNA targets through TargetScan and miRDB databases, yielding
   2309 mRNAs by intersecting with the 5065 DEmRNAs, accounting for 45.6%
   of the total DEmRNAs. (Fig. [75]2A). GO and KEGG enrichment analyses
   were conducted on these mRNAs, with a significant threshold of
   corrected p < 0.05. GO enrichment analysis revealed enrichment in
   development and structural organization GO terms in Biological
   Processes (BP), including urogenital system development (GO: 0001655),
   external encapsulating structure organization (GO: 0045229), nuclear
   division (GO: 0000280), renal system development (GO: 0072001), and
   extracellular matrix organization (GO: 0030198). Cellular Component
   (CC) included providing structural organization and cellular
   compartmentalization enriched GO terms, mainly associated with
   collagen-containing extracellular matrix (GO: 0062023), synaptic
   membrane (GO: 0097060), basement membrane (GO: 0005604), condensed
   chromosome (GO: 0000793), and chromosome, centromeric region (GO:
   0000775). Molecular Functions (MF) were enriched in fundamental
   biological processes and molecular interactions GO terms, with notable
   terms such as extracellular matrix structural constituent (GO:
   0005201), transmembrane receptor protein kinase activity (GO: 0019199),
   transmembrane receptor protein tyrosine kinase activity (GO: 0004714),
   glycosaminoglycan binding (GO: 0005539), and DNA-binding transcription
   activator activity (GO: 0001216). Additionally, fundamental cellular
   processes and signaling KEGG pathways were identified, including Cell
   cycle including Cell cycle (hsa04110), Focal adhesion (hsa04510),
   PI3K-Akt signaling pathway (hsa04151), ECM-receptor interaction
   (hsa04512), and Axon guidance (hsa04360) (Fig. [76]2B). Finally, a
   triple regulatory network comprising 192 lncRNAs, 89 miRNAs, and 2309
   mRNAs (Supplementary Fig. 1) was constructed using Cytoscape. Moreover,
   we conducted batch survival regression analysis on the network, and a
   triple regulatory network associated with BC prognosis-related genes
   was generated (Fig. [77]2C).

Fig. 2.

   [78]Fig. 2
   [79]Open in a new tab

   Construction of the lncRNA-miRNA-mRNA Network and Its Functional
   Enrichment Analysis. A Venn diagrams display the overlap between
   DEmiRNAs and the miRNAs targeted by DElncRNAs, as well as the
   intersection between DEmRNAs and the mRNAs targeted by shared miRNAs
   from the previous overlap. B Visualization of the Gene Ontology (GO)
   and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment
   analyses for the intersecting DEmRNAs and their predicted mRNA
   counterparts. C A highlighted subset of genes within the ceRNA
   co-regulatory network that have associations with breast cancer
   prognosis

Identification of a prognostic ceRNA network in BC

   To identified key ceRNA with significant prognosis value in BC, we
   analyzed the status and prognostic value of hub triple regulatory
   network RNAs associated with prognosis in BC. Our results showed 6
   upregulated and 3 downregulated lncRNAs, 6 upregulated and 4
   downregulated miRNAs, 33 upregulated mRNAs, and 18 downregulated mRNAs
   in BC compared to adjacent normal breast tissues (Fig. [80]3). Survival
   analysis revealed the high expression of 1 lncRNA (LINC01588), 2 miRNAs
   (hsa-miR-9-5p and hsa-miR-1307-3p), and 7 mRNAs (FAM83D, SDC1, ANLN,
   CCNE1, CENPA, LRP8, and CENPN) was associated with the worst OS, DSS,
   and PFI prognosis, while high expression of 5 lncRNAs (MAPT-IT1,
   LINC01224, SLC26A4-AS1, LINC02568, and DIO3OS) and 14 mRNAs (CCL19,
   ABCB1, NEFH, IL2RB, KRT17, FCRLA, GREB1L, NTRK2, NR3C2, ETNPPL, ABAT,
   ZNF703, DMD, and CCL5) was associated with best prognosis (Fig. [81]4).

Fig. 3.

   [82]Fig. 3
   [83]Open in a new tab

   Distribution of Expression Patterns for 73 Prognostically-Associated
   Hub-RNAs Derived from the Triple Regulatory Network in the TCGA BRCA
   Dataset. A Expression patterns of 9 pivotal long non-coding RNAs
   (hub-lncRNAs) in BRCA samples. B Expression distribution of 10 key
   microRNAs (hub-miRNAs) within BRCA samples. C Expression trends of 54
   central messenger RNAs (hub-mRNAs) observed in BRCA samples

Fig. 4.

   [84]Fig. 4
   [85]Open in a new tab

   Prognostic Efficacy of the 73 Hub-RNAs Derived from the Triple
   Regulatory Networks. A Forest plots displaying the overall survival,
   disease-specific survival, and progression-free survival of lncRNAs. B
   Forest plots presenting the survival outcomes for miRNAs, encompassing
   overall survival, disease-specific survival, and progression-free
   survival. C Forest plots illustrating the survival profiles of mRNAs
   with regards to overall survival, disease-specific survival, and
   progression-free survival

   Considering the upregulated LINC001588 interacted with the upregulated
   hsa-miR-9-5p, and SLC26A4-AS1 was mainly localized in the cytoplasm,
   whereas LINC001588 was mainly distributed in the cytosol (Fig. [86]5A).
   Thus, SLC26A4-AS1 may act as a ceRNA in BC, sequestering hsa-miR-19a-3p
   to regulated mRNAs expression. Furthermore, expression correlation
   analysis showed that ADRB1, ABCB1, MDFIC, NTRK2, NR3C2, GABRA2, and
   SOCS3 were positively correlated with SLC26A4-AS1 and negatively
   correlated with hsa-miR-19a-3p expression (Fig. [87]5B), among which,
   we observed a significant decrease in NTRK2 total protein levels in BC,
   consistent with mRNA expression levels using data from CPTAC
   (Fig. [88]5C). Therefore, we constructed the
   SLC26A4-AS1-hsa-miR-19a-3p-NTRK2 ceRNA network. Target prediction of
   SLC26A4-AS1 and NTRK2 3' UTRs paired with miR-19a-3p was conducted
   using the starBase version 2.0 database
   ([89]http://starbase.sysu.edu.cn) (Fig. [90]5D). We also observed a
   significant downregulation of NTRK2 in BC patients and breast cancer
   cell lines in the GEO databases [91]GSE45827 and [92]GSE12790, further
   validating NTRK2 expression. Thus, we selected the SLC26A4-AS1/NTRK2
   axis in the ceRNA network as a potential prognostic model for further
   analysis (Fig. [93]5E).

Fig. 5.

   [94]Fig. 5
   [95]Open in a new tab

   Development and Comprehensive Analysis of the ceRNA Network. A
   Inference of the cellular localization for SLC26A4-AS1 and LINC01588. B
   Correlational analysis involving SLC26A4-AS1, hsa-miR-19a-3p, and
   associated mRNAs in BRCA. C Quantitative representation of the protein
   expression of genes interacting with SLC26A4-AS1 in normal and primary
   tumor tissues of BRCA; other mRNAs were undetected in the Breast
   dataset. D Conceptual model delineating the ceRNA network, accompanied
   by predicted binding site indications. Overexpressed genes are marked
   in red, and underexpressed genes in blue. E Empirical validation of
   NTRK2 expression within breast tissues and breast cancer cell lines,
   utilizing the GEO database.onstruction and integrated analysis of the
   ceRNA network. Notations: *p < 0.05, **p < 0.01, ***p < 0.001

   Clinical relevance of the SLC26A4-AS1- hsa-miR-19a-3p-NTRK2 ceRNA
   network in BC patients.

   SLC26A4-AS1, hsa-miR-19a-3p, and NTRK2 are candidate genes for
   constructing BC specific prognostic biomarkers. Grouping gene
   expression by median, K-M survival analysis showed that patients with
   high expression of SLC26A4-AS1 (Fig. [96]6A) and NTRK2 (Fig. [97]6B)
   had higher OS than those with low expression, while patients with high
   expression of hsa-miR-19a-3p had lower OS than those with low
   expression (Fig. [98]6C). It is noteworthy that SLC26A4-AS1, NTRK2, and
   hsa-miR-19a-3p do not significantly predict prognosis in the early
   stages (I-II) of BC but have a strong predictive capability for the
   prognosis in late-stage (III-IV; Supplementary Fig. 2). Time-dependent
   ROC curve analysis confirmed that SLC26A4-AS1, hsa-miR-19a-3p, and
   NTRK2 have significant prognostic value (Fig. [99]6). We also analyzed
   correlations between SLC26A4-AS1- hsa-miR-19a-3p-NTRK2 ceRNA network
   and the pathological stage of BC patients, and we found that low NTRK2
   expression was significantly correlated with earlier tumor stage
   (I-II), while expression levels of SLC26A4-AS1 and hsa-miR-19a-3p were
   not (Supplementary Fig. 2A–C).

Fig. 6.

   [100]Fig. 6
   [101]Open in a new tab

   Prognostic Assessment of Predictive RNAs in BRCA. A Comprehensive
   analysis encompassing expression profiling, Kaplan–Meier survival
   estimation, and time-dependent receiver operating characteristic (ROC)
   curve evaluation for SLC26A4-AS1. B Parallel evaluations for NTRK2. C
   Analogous analyses for hsa-miR-19a-3p

   Furthermore, we explored the correlation between the expression levels
   of SLC26A4-AS1, hsa-miR-19a-3p, and NTRK2 with clinical factors. Our
   study showed that the expression levels of SLC26A4-AS1, hsa-miR-19a-3p,
   and NTRK2 were significantly associated with PR status and ER status
   (p < 0.001), as shown in Supplementary Table 1. However, no significant
   correlation was observed between the expression levels of these genes
   with other clinical factors, including age, menopause status,
   pathologic T stage, pathologic N stage, pathologic M stage, pathologic
   stage, HER2 status, and anatomic neoplasm subdivisions (p > 0.05). In
   addition, in the univariate cox regression analysis models of
   SLC26A4-AS1, hsa-miR-19a-3p, and NTRK2, some prognostic factors (Age,
   Pathologic T stage, Pathologic N stage, Pathologic M stage, Pathologic
   stage, Menopause status) were closely associated with OS in BC patients
   in the TCGA cohort (p < 0.05). Importantly, in the multivariate Cox
   regression analysis, aberrant expression of SLC26A4-AS1 remained
   associated with OS in breast cancer patients, while aberrant expression
   of hsa-miR-19a-3p and NTRK2 was not related to poor prognosis
   (Supplementary Fig. 3A-B). Therefore, SLC26A4-AS1 may serve as an
   independent prognostic factor in BC patients.

The methylation status of NTRK2 in BC

   To further elucidate the mechanism underlying the downregulation of
   NTRK2 in BC tissues, we explored the correlation between NTRK2
   expression levels and its methylation status through various methods.
   Firstly, differential expression analysis of three DNA
   methyltransferases (DNMT1, DNMT3A, and DNMT3B) between breast cancer
   NTRK2-high and NTRk2-low groups showed a significant increase in DNMT3B
   expression in the NTRK2-low group compared to the NTRK2-high group,
   although DNMT1 and DNMT3A exhibited no significant differences
   (Fig. [102]7A). Secondly, UALCAN analysis revealed that the methylation
   levels of NTRK2 were higher in breast cancer tissues than in normal
   breast tissues (Fig. [103]7B). Similarly, analysis using DiseaseMeth
   2.0 showed that the methylation levels of NTRK2 were significantly
   higher in BC tissues compared to normal tissues (Fig. [104]7C).

Fig. 7.

   [105]Fig. 7
   [106]Open in a new tab

   Methylation analysis of NTRK2. A Comparative expression of three DNA
   methyltransferases (DNMT1, DNMT3A, and DNMT3B) between high-NTRK2 and
   low-NTRK2 expressing BRCA samples. B Analysis of NTRK2 methylation
   patterns in BRCA via the UALCAN platform. C Examination of NTRK2
   methylation status in BRCA utilizing the DiseaseMeth version 2.0 tool.
   D Discrepancies in methylation intensities of CpGs within NTRK2 between
   BRCA and normal tissue samples. Correlational investigation between
   methylation values across multiple CpG sites and NTRK2 expression
   magnitude in BRCA. E Heatmap representation highlighting the clustering
   of CpG methylation degrees within NTRK2 in BRCA, processed through the
   Methsurv database. Methylation levels are denoted continuously, with a
   spectrum ranging from blue (fully unmethylated, 0) to red (fully
   methylated, 1). Each row signifies a specific CpG, whereas columns
   represent individual patients. F Survival implications depicted via a
   forest plot, focused on diverse CpG sites within NTRK2 in BRCA.
   Notations: ns indicates no significance; *p < 0.05, **p < 0.01,
   ***p < 0.001

   Interestingly, through analysis using the SMART App, we found that
   several CpGs within the NTRK2 gene exhibited higher methylation levels
   in BC samples compared to normal sample, including cg01009697,
   cg03628748, and cg08470639. Furthermore, the methylation (beta-value)
   of these CpGs was negatively correlated with NTRK2 expression. In
   addition, when compared to normal tissues, cg13620631, cg14273545, and
   cg13654445 exhibited lower methylation levels in BC, but except for
   cg13620631, cg14273545, and cg13654445 methylation was positively
   correlated with NTRK2 expression (Fig. [107]7D). Subsequently, we
   utilized MethSurv to present the methylation status of multiple CpGs
   within the NTRK2 gene in BC in the form of a heatmap (Fig. [108]7E). It
   is noteworthy that high methylation of CpGs cg13620631, cg13698224,
   cg14273545, and cg01009697 was associated with higher overall survival
   rates in TCGA breast cancer patients, while high methylation of
   cg01292475 and cg13965062 associated with lower overall survival rates
   in breast cancer patients (Fig. [109]7F).

Correlation between NTRK2 expression and immune infiltration in BC

   Tumor-infiltration cells are an important component of the tumor
   microenvironment and are closely associated with the occurrence,
   development, or metastasis of cancer. We assessed the potential
   correlation between NTRK2 expression and immune infiltration levels in
   breast cancer. The results showed that NTRK2 was negatively correlated
   with the infiltration levels of Th2 cells, Treg, aDC, and NK CD56dim
   cells, but positively correlated with CD8 T cells, DC, Eosinophils,
   iDC, Mast cells, Neutrophils, NK CD56bright cells, NK cells, pDC, T
   cells, T helper cells, Tcm, Tem, Th17 cells, and TReg (Fig. [110]8A).
   Additionally, breast cancer tissues with low NTRK2 expression showed
   lower infiltration levels of B cells, CD8 T cells, DC, Eosinophils,
   iDC, Mast cells, Neutrophils, NK CD56bright cells, NK cells, pDC, T
   cells, T helper cells, Tcm, Tem, and TFH, while aDC, TReg, and Th2
   cells showed higher infiltration levels (Fig. [111]8B).

Fig. 8.

   [112]Fig. 8
   [113]Open in a new tab

   Immune Landscape Pertaining to NTRK2 in BRCA. A Association of NTRK2
   expression with the extent of immune cell infiltration in BRCA tissues.
   B Discrepancies in immune cell infiltration between NTRK2
   high-expressing and low-expressing BRCA cohorts. C Elucidation of the
   interrelation between NTRK2 expression and cytotoxic T-lymphocyte (CTL)
   infiltration. D Exploration of the interplay between NTRK2 expression
   and various parameters: T-cell exclusion metric, T-cell dysfunction
   index, responsiveness to immune checkpoint blockade (ICB)
   interventions, and the phenotypic manifestation in CRISPR-based gene
   ablation studies. Notations: ns indicates no significance; *p < 0.05,
   **p < 0.01, ***p < 0.001

   There are two main mechanisms of immune evasion in cancer: (1) inducing
   T-cell dysfunction in tumors with high cytotoxic T lymphocyte (CTL)
   infiltration, and (2) preventing T-cell infiltration in tumors with low
   CTL levels. To predict whether NTRK2 was involved in BC immune evasion,
   we explored the correlation between NTRK2 and CTL using the TIDE
   database. The results showed that NTRK2 was negatively correlated with
   CTL infiltration in BC (Fig. [114]8C). Figure [115]8D showed that NTRK2
   was associated with a positive T-cell dysfunction score in the Metabric
   dataset. Additionally, NTRK2 expression levels were higher in
   tumor-associated M2 macrophages (TAM M2) and lower in CAF FAP, cell
   types known to promote T-cell exclusion. These findings suggest that
   NTRK2 may contribute to immune dysregulation and facilitate immune
   evasion in BC patients through mechanisms involving T-cell exclusion.

Exploring NTRK2-associated proteins and pathways in BC

   To further screen NTRK2-targeted binding proteins and investigate the
   molecular mechanism of NTRK2 in BC, we retrieved 40 NTRK2 binding
   proteins (medium confidence) from the STRING database. Figure [116]9A
   displays the protein–protein interaction network. KEGG pathway
   enrichment analysis revealed that neurotrophin signaling pathway, Ras
   signaling pathway, and PI3K-Akt signaling pathway are the most
   significantly involved in NTRK2's impact on BC (Fig. [117]9B).
   Furthermore, the expression level of NTRK2 remained unaffected in
   breast cancer cell lines treated with a PI3K-Akt signaling pathway
   inhibitor (Fig. [118]9C), suggesting that NTRK2 may be upstream of the
   PI3K-Akt signaling pathway. Additionally, GO enrichment analysis
   related to NTRK2, including biological processes (Fig. [119]9D),
   cellular components (Fig. [120]9E), and molecular functions
   (Fig. [121]9F), mainly enriched in signaling, synaptic function, and
   kinase activity regulation. Finally, the results of SLC26A4-AS1,
   hsa-miR-19a-3p, and NTRK2 expression were validated in breast cancer
   cell MDA-MB-468 and breast epithelial cell line MCF-10A (Fig. [122]9G).
   However, the observed changes in NTRK2 protein in breast cancer cell
   were not parallel changes in NTRK2 RNA levels (Fig. [123]9H).

Fig. 9.

   [124]Fig. 9
   [125]Open in a new tab

   Exploration of NTRK2-Associated Gene Enrichment and Expression
   Validation. A Correlational scrutiny of 20 experimentally ascertained
   NTRK2-binding proteins in BRCA samples. B KEGG pathway exploration
   anchored on NTRK2-associated genes. C Evaluative analysis of the
   influence of PI3K-Akt pathway inhibitors on NTRK2 expression in JIMT1
   and HCC1954 breast cancer cell lines. D–F GO enrichment analysis for
   NTRK2-associated genes distributed across Biological Processes D,
   Cellular Component E, and Molecular Function (F). G Quantitative
   assessment of mRNA expression levels for SLC26A-AS1, hsa-miR-19a-3p,
   and NTRK2 employing RT-qPCR. H The protein expression magnitude of
   NTRK2 ascertained through Western blotting analysis. Notations: ns
   indicates no significance; *p < 0.05, **p < 0.01, ***p < 0.001

Discussion

   Breast cancer is one of the three most common cancers worldwide. While
   early-stage breast cancer is considered curable, advancements in
   treatment have shifted the focus towards reducing both overtreatment
   and undertreatment. Identifying reliable biomarkers is crucial for
   developing new therapeutic targets and improving patient prognosis.
   Increasing evidence has revealed that non-coding RNAs (ncRNAs) play
   pivotal roles in breast cancer progression. However, the intricate
   regulatory interplay between different types of ncRNAs remains poorly
   understood. In this study, we established a lncRNA-miRNA-mRNA
   tri-regulatory network in BC, identifying 192 lncRNAs, 89 miRNAs, and
   2,309 mRNAs. Enrichment analysis highlighted the network’s focus on key
   processes such as "extracellular matrix organization" and
   "extracellular matrix receptor interaction." This suggests that these
   ncRNA interactions may influence tumor progression by modulating the
   extracellular matrix, contributing to either tumorigenic or
   tumor-suppressive effects. Subsequent prognostic analysis identified a
   core ceRNA network of 9 lncRNAs, 10 miRNAs, and 54 mRNAs, among which
   the SLC26A4-AS1-hsa-miR-19a-3p-NTRK2 axis emerged as significantly
   associated with prognosis in BC. Expression profiling and survival
   analyses further confirmed that low expression of this axis correlates
   with poor clinical outcomes. The roles of SLC26A4-AS1 and NTRK2 in
   cancer have been reported previously. Yuan et al. [[126]20] found that
   lncRNA SLC26A4-AS1 was significantly downregulated in thyroid cancer
   and lymph node metastasis specimens, and low expression was linked to
   poor prognosis. They also demonstrated that SLC26A4-AS1 inhibits the
   MRE11/RAS50/NBS1 complex-mediated DNA repair signal, thereby
   suppressing thyroid cancer metastasis by disrupting DDX5. In addition,
   SLC26A4-AS1 has been shown to promote autophagy and mitigate thyroid
   papillary carcinoma progression by recruiting ETS1 to enhance ITPR1
   expression [[127]21]. Similarly, Li et al. [[128]22] identified
   SLC26A4-AS1 downregulation in gastric cancer, where it was also
   significantly associated with poor prognosis. Consistent with these
   findings in other cancers [[129]23, [130]24], our study demonstrates
   that low expression of SLC26A4-AS1 in BC is linked to poor prognosis,
   and for the first time, we validated this in the breast cancer cell
   line MDA-MB-468.

   Currently, research on hsa-miR-19a-3p mainly focuses on non-cancer
   domains. For instance, reduced expression of hsa-miR-19a-3p was found
   in premenopausal women with genital arousal disorder [[131]25] and in
   patients with spinal cord injury-induced neuropathic pain [[132]26],
   while its expression was aberrantly increased in patients with vitiligo
   [[133]27]. We are the first to identify and validate the abnormally
   high expression of hsa-miR-19a-3p in BC, which is associated with poor
   prognosis in BC patients. This suggests that the aberrant expression of
   hsa-miR-19a-3p is related to the onset and development of human
   cancers.

   The Tropomyosin Receptor Kinases (TRKs) were first identified in 1986
   when Neurotrophic Receptor Tyrosine Kinase 1 (NTRK1) was discovered to
   be a part of an oncogenic fusion gene in colorectal cancer. Shortly
   thereafter, NTRK2 and NTRK3 were also identified. NTRK2, also known as
   Trk-β, encodes the TRKB protein, which binds to neurotrophins and
   regulates neuronal development and function [[134]28]. It has been
   reported that NTRK fusions are common in some rare types of tumors,
   such as secretory carcinoma of adult breasts and salivary glands, as
   well as infantile fibrosarcoma (IFS) in pediatrics, with reported
   fusion positivity rates often exceeding 75% [[135]29, [136]30].
   Furthermore, drugs targeting NTRK fusions are currently in clinical
   trials, with some, like TRK inhibitors larotrectinib and entrectinib,
   being approved as tumor-agnostic treatments for solid tumors harboring
   NTRK fusions [[137]31]. In addition, Gomez et al. [[138]32] found that
   overexpression of NTRK2 in non-small cell lung cancer is associated
   with in vitro proliferation and invasion. Hu et al. [[139]33] concluded
   that NTRK2, as a target oncogene of miR-22, is involved in the in vitro
   development of gastric cancer. However, current research on NTRK2 in
   breast cancer is limited, and its expression in breast cancer differs
   from its expression in other tumors. A recent study showed that NTRK2
   is abnormally downregulated in breast cancer cell lines MCF-7 and
   MDA-MB-231 [[140]34]. Similarly, based on the TCGA database, our
   results show that the expression of SLC26A4-AS1 and NTRK2 in BC is
   significantly lower than in normal tissue. Survival analysis revealed
   that low expression of SLC26A4-AS1 and NTRK2 is associated with poor
   prognosis. In contrast, the expression of hsa-miR-19a-3p is
   upregulated, and its high expression level indicates a poorer OS
   prognosis for BC. Additionally, we validated the mRNA levels of NTRK2
   in the [141]GSE12790 and [142]GSE45827 datasets.

   DNA methylation is the most extensively studied epigenetic
   modification, playing a pivotal role in promoting embryonic
   development, genomic imprinting, and X-chromosome inactivation, among
   other vital biological processes. Aberrations in DNA methylation can
   lead to alterations in the cellular microenvironment, impacting gene
   expression patterns, ultimately giving rise to various pathological
   conditions, including carcinogenesis [[143]35]. In recent years,
   emerging non-invasive DNA methylation biomarkers have shown significant
   clinical value in predicting cancer prognosis and drug responses
   [[144]36]. Therefore, we consulted various online sources to explore
   DNA methylation patterns that could elucidate the aberrant expression
   of NTRK2 in BC. We discovered that, compared to adjacent normal
   samples, the methylation level of NTRK2 in BC samples is elevated. This
   is consistent with the downregulation of NTRK2 in BC observed when
   combined with DNA methyltransferase analysis (DNMT1, DNMT3A, and
   DMNT3B; Supplementary Fig. 4). We also noticed that the low expression
   group of NTRK2 significantly co-occurred with high expression of
   DNMT3B, providing insights into the downregulation of NTRK2 in BC.
   Additionally, compared to normal tissue, we identified specific
   methylation sites of NTRK2 in BC samples that displayed higher
   methylation levels and showed a negative correlation with NTRK2
   expression. In contrast, some sites with low methylation in BC
   positively correlated with NTRK2 expression and negatively correlated
   with the prognosis of BC patients. Lastly, we observed that high
   methylation in the TSS1500 region (g01009697), body region (cg13698224
   and cg14273545), and the 3'UTR region (cg13620631) of NTRK2 resulted in
   a higher overall survival rate for TCGA breast cancer patients.
   Conversely, elevated methylation in the TSS200 region (cg01292475) and
   5'UTR (cg13965062) resulted in a lower overall survival rate. These
   findings suggest that aberrant methylation may be associated with BC
   prognosis, and more evidence is needed to validate the potential role
   of NTRK2 DNA methylation in the onset of BC tumors.

   The tumor microenvironment (TME) is a principal regulator of tumor
   initiation, progression, and therapeutic responses. The TME comprises
   stromal cells (including myoepithelial cells, fibroblasts, and vascular
   endothelial cells), immune cells (such as T cells, B cells,
   macrophages, and natural killer cells), and the extracellular matrix.
   Most cancer pathobiology research has predominantly focused on changes
   and drivers within cancer cells, with relatively less attention paid to
   the surrounding milieu. Compared to other cancers like melanoma or lung
   cancer, the breast cancer microenvironment exhibits significant
   cellular and spatial heterogeneity [[145]37]. Tumor-infiltrating
   lymphocytes (TILs) are immune cells that migrate to malignant lesions
   as part of the immune response promoting tumor elimination. They
   encompass various immune cells, primarily T cells (CD8 + cytotoxic T
   cells and CD4 + helper or regulatory T cells (Tregs)), as well as B
   cells and NK cells. A systematic review encompassing 15 studies
   involving nearly 14,000 breast cancer patients revealed the highest
   level of TIL infiltration in triple-negative breast cancer (TNBC),
   followed by HER2-positive breast cancer, with the luminal subgroups
   exhibiting the lowest TIL infiltration [[146]38]. Moreover, in a
   retrospective analysis of over 5,000 female breast cancer samples, low
   Treg abundance correlated with pCR post-neoadjuvant chemotherapy in
   TNBC women but was unrelated to ER-positive/HER2-negative disease
   [[147]39]. These findings indicate the presence of immune infiltrating
   cells varies considerably within breast cancer subtypes and between
   different subtypes. Dendritic cells (DCs) play a pivotal role in
   antigen presentation during immune surveillance and orchestrating
   anti-tumor immune responses. DCs isolated from breast cancer tumors
   [[148]40] and peripheral blood [[149]41] of breast cancer patients have
   been shown to exhibit impaired maturation and reduced T cell
   stimulatory activity. Tregs in the TME exhibit immunosuppressive
   actions, promoting the growth and spread of malignant tumors. Through
   secreting inhibitory cytokines, Tregs can foster intra-tumoral T cell
   exhaustion by modulating the expression of inhibitory receptors in
   CD8 + TILs and exhaustion-associated transcriptomic features [[150]42].
   In our study, we found low NTRK2 expression correlates with reduced TIL
   and DC infiltration, but a higher degree of TReg infiltration. Tumors
   with high cytotoxic T lymphocyte (CTL) infiltration can induce T cell
   dysfunction, leading to tumor immune evasion. Our results also show
   NTRK2 negatively correlates with CTL infiltration in BC, with elevated
   expression levels in tumor-associated M2 macrophages (TAM M2) and
   reduced expression in CAF FAP. These findings suggest aberrant NTRK2
   expression relates to immune dysregulation in BC patients and may
   promote immune evasion via T cell exclusion. Consequently, the role and
   functionality of NTRK2 in immunotherapy for BC warrant further
   investigation.

   To further explore the biological functions of NTRK2, we conducted GO
   and KEGG enrichment analyses for NTRK2 and its associated proteins. The
   enriched pathway related to NTRK2 is the PI3K-Akt signaling pathway,
   and analysis of the [151]GSE69189 dataset suggests that NTRK2 might be
   positioned upstream in the PI3K-Akt signaling pathway, consistent with
   earlier research results [[152]43]. GO enrichment analyses associated
   with NTRK2 (including BP, CC, and MF) were primarily enriched in
   signaling, synaptic function, and regulation of kinase activity.
   Lastly, we confirmed the mRNA expression of SLC26A4-AS1,
   hsa-miR-19a-3p, and NTRK2 in breast cancer cell line MDA-MB-468 and
   breast epithelial cell line MCF-10A. However, we did not observe any
   differential expression of NTRK2 protein between breast cancer cell
   lines and normal breast cells. The inconsistency between NTRK2 mRNA and
   protein levels suggests that reduced NTRK2 RNA levels might be common,
   but they may not correlate with the actual total protein expression or
   the expression of NTRK2 in BC might involve complex transcriptional and
   post-transcriptional regulatory mechanisms involving non-coding RNA, an
   issue that deserves further investigation.

   Although we have constructed the SLC26A4-AS1/NTRK2 axis based on the
   ceRNA mechanism, which appears to be a potential prognostic biomarker
   for clinical application, certain limitations must be acknowledged.
   First, the binding affinities of lncRNAs, miRNAs, and mRNAs retrieved
   from databases require further experimental validation. Secondly, the
   function and mechanism of the SLC26A4-AS1/NTRK2 axis in the initiation
   and progression of BC need to be further explored through in vitro and
   in vivo studies.

   In conclusion, we have established a low-expression ceRNA network
   related to BC prognosis (SLC26A4-AS1-hsa-miR-19a-3p-NTRK2), which
   offers deeper insights into the correlations among lncRNA-miRNA-mRNA.
   Additionally, we found that the SLC26A4-AS1/NTRK2 axis based on the
   ceRNA mechanism might be a novel and critical prognostic factor
   involved in BC. This prognostic model contributes to our understanding
   of the pathogenic mechanism of BC.

Supplementary Information

   [153]Supplementaty material 1.^ (21.8MB, docx)
   [154]Supplementaty material 2.^ (22.3KB, docx)

Acknowledgements