This article is mentioned in:

Introduction

Breast cancer (BC) is one of the most common malignancies worldwide, with 2.26 million new cases in 2022, representing 23.8% of all cases of cancer in women (1). BC is the leading cause of cancer-associated mortality in women, having caused 666,103 mortalities in 2022 worldwide (1). In Mexico, 31,043 new cases and 8,195 mortalities from BC were reported in 2022 (1), making it one of the most prevalent and fatal types of cancer in Mexico. The estimated national incidence of BC age-standardized rate (ASR) is 40 cases per 100,000 person-years (1). Likewise, the first population-based cancer registry in Southeastern Mexico revealed an ASR incidence of 49.3 cases per 100,000 person-years (2). Despite the high incidence of BC, emerging therapeutic targets, diagnostic and prognostic biomarkers have managed a 5-year overall survival rate of 90% in the USA (3,4). However, long-term chemotherapy and targeted therapy can induce resistance and cancer progression (5). Nevertheless, for individuals without social security in emerging countries, such as Mexico, 55% BC cases are typically diagnosed already at locally advanced stages, with a 5-year survival of 69.9% (6). By contrast, patients with social security, who often have greater access to comprehensive medical care and timely treatments, have a 5-year overall survival rate of 90.4% (7); however, upon analyzing data according to stages of the disease, a progressive decrease in survival can be observed. Stage I has a survival rate of 98.8%, stage II has a survival rate of 97.5%, stage III has a survival rate of 89.4% and stage IV has a survival rate of 22.7% (7).

BC is highly heterogeneous and is classified based on expression of the estrogen receptor (ER), progesterone receptor (PR) and human epidermal growth factor receptor 2 (HER2) (8). Molecular classification of BC based on global expression analysis has contributed to determining prognosis and appropriate treatment options (9–13). BC tumors can be divided into luminal A, luminal B, HER2, basal, claudin-low and normal-like (14). Basal and claudin-low subtypes are clinically known as triple-negative BC (TNBC), an aggressive subtype of BC characterized by the lack of ER, PR and HER2 expression (15) Clinically, TNBC presents additional challenges compared with other subtypes of BC, including high aggressiveness and a high recurrence rate due to the lack of specific therapeutic options and molecular targets for targeted therapy, exerting a considerable negative impact on the survival of patients with TNBC (14).

Consequently, efforts have been focused on characterizing TNBC to identify predictive biomarkers and therapeutic targets. Sub-classification of TNBC has been proposed, namely basal-like 1, basal-like 2, mesenchymal, mesenchymal stem-like, immunomodulatory and luminal androgen receptor (16). Molecular classification, biomarkers and therapeutic approaches have been designed based on the profile of protein-coding RNAs. However, these transcripts only represent 2% of the total transcriptome in a human cell (17). The remaining 98% of transcripts consist of non-coding RNAs (ncRNAs) that may exert relevant clinical and biological functions. A diverse range of non-coding RNAs have been identified in tumor cells, among which long non-coding RNAs (lncRNAs) represent a tissue-specific group of transcripts exhibiting diverse functions in tumor progression and metastasis (18). LncRNAs are transcripts consisting of >200 nucleotides that lack protein coding potential (19). These transcripts can interact with DNA, mRNA, microRNAs (miRNA) and protein molecules to control gene and protein expression on various levels, such as epigenetic, transcriptional, post-transcriptional, translational and post-translational levels (20). LncRNAs can be located in the nucleus or cytoplasm of a cell (21). In the cytoplasm, lncRNAs can modulate mRNA stability, regulate translation, mediate protein modifications, act as competing endogenous RNAs, function as a sponge of miRNA by binding miRNAs and sequestering them and can also be precursors of miRNAs (22,23). By contrast, nuclear lncRNAs can control chromatin states, transcriptional and posttranscriptional gene expression, RNA stability and RNA processing and splicing (24). LncRNAs exert a plethora of functions in a cell, implying roles in a variety of normal and pathological cellular processes, including tumorigenesis, invasion and metastasis (18). Expression of lncRNAs is commonly dysregulated during tumorigenesis, which suggests oncogenic or tumor-suppressive functions and chemo-resistance (25,26). Several studies have previously focused on describing the role of lncRNAs in the development and progression of cancer, such as breast, colon, melanoma, pancreas and liver cancers (27), which have been attributed to oncogenic or tumor-suppressor functions (28–30), such as metastasis-associated lung adenocarcinoma transcript 1, homeobox transcript antisense RNA, urothelial cancer associated 1, growth arrest-specific transcript 5, nuclear paraspeckle assembly transcript and lincRNAp21.

The present study aimed to identify the molecular profile of lncRNAs in a series of BC tumors from Mexican patients, according to the expression of routine biomarkers ER, PR, HER2. TNBC (ER-, PR- and HER2-) was then compared with luminal B HER2+ (ER+, PR+ and HER2+) BC. Previous studies have reported an increased prevalence of TNBC tumors in Mexico compared with USA populations (31,32). In these studies, it was reported that the prevalence of TNBC in Mexico was 23.1%, compared with the 10–13% in USA, making it essential to find suitable prognostic biomarkers and therapeutic options. The aim of the present study was to find a lncRNA signature associated with TNBC prognosis, since TNBC is known to be the most aggressive form of BC (15). Differentially expressed lncRNAs were validated using a cohort from The Cancer Genome Atlas (TCGA), through a cohort obtained from the Gene Expression Omnibus database (accession no. GSE134359), and an independent cohort of Mexican patients with BC. These samples were different from those used in the microarray analyses of the present study; although they were selected based on the same characteristics, they were independent samples collected from UMAE Oncology Hospital. It is hoped that this approach can allow for the identification of lncRNA profiles with a potential impact on the overall survival of patients with TNBC.

Materials and methods

Sample collection

Approval for the present study was obtained from the Research and Ethics committee of the Mexican Institute of Social Security under approval nos. R-2013-785-045 and R-2020-785-154. A total of 192 BC samples with histopathological reports were collected from pathology records at the Unidad Médica de Alta Especialidad Oncology Hospital, XXI National Medical Center of the Mexican Institute of Social Security (Mexico City, Mexico). These samples were collected retrospectively from cases registered between January 2009 and December 2013, following project authorization in July 2013. Clinical data, including age, sex, grade and stage, were compiled from the medical records at the UMAE Oncology Hospital (Table SI).

Histological samples processing

BC tissues were processed at the Pathology department, UMAE Oncology Hospital (Mexico City, Mexico) as follows: Tumor tissue biopsies were previously fixed in 10% neutral formaldehyde buffer for 24 h at room temperature. Subsequently, the tissue was processed using a tissue processor. The tissues were dehydrated for 1 h at room temperature in each of the following solvents: Ethanol 70%, two changes of 96% ethanol, two changes of absolute ethanol, one change of a 50% mixture of absolute ethanol and absolute xylene, two changes of xylene for clearing and three changes of paraffin at 58°C, concluding with the embedding process. Afterwards, tissue sections were cut to 5-m thickness using a microtome and deparaffinized at 58°C for 20 min and immersed in xylene. The tissues were then rehydrated by immersion in alcohol baths of 100, 96 and 70% alcohol, followed by water.

H&E staining

Following the deparaffinization and rehydration of the tissue sections, they were stained with hematoxylin at room temperature for 6 min, rinsed with water and differentiated with acid alcohol. Afterward, they were rinsed in tap water, blued in mildly alkaline water and rinsed with tap water. Subsequently, eosin staining was performed at room temperature for 1 min. The slides were dehydrated in baths of 70 and 96% alcohol, absolute alcohol and xylene. The slides were then mounted using synthetic resin. Light microscopic evaluation of the samples was performed by a pathologist (author AM) for diagnosis.

Tumor microarray construction

For subsequent microarray (n=30) and reverse transcription-quantitative PCR (RT-qPCR) (n=21) analysis, only samples that met the following criteria were used: Sufficient tumor tissue, samples belonging to the TNBC or luminal B HER2 subtypes and samples with a high RNA quality. A high RNA quality was classed as an absorbance value at 260/280 of 1.9–2.2, as well as an RNA Integrity Number (RIN) >6. Histological evaluation was performed by a pathologist (author AM) who confirmed diagnosis and selected representative fragments for the construction of tissue microarrays (TMAs). TMAs were built using Chemicon ATA-100 equipment (Leica Biosystems) and representative cylinders of 1 mm in diameter. Histological sections that were 4-mm thick were cut and mounted on electrocharged slides (cat. no. AMS90-Color; Cancer Diagnostics, Inc.) for H&E staining and marker detection by immunohistochemistry.

Immunohistochemistry

Immunohistochemistry was performed on Ventana BenchMark XT using the UltraView Universal DAB Kit (cat. no. 760-500; Roche Tissue Diagnostics) following the provider's protocol. Key steps included a deparaffinization process using EZprep at 76°C in 4-min intervals over a total of 16 min, followed by antigen exposure in cell conditioning solution 1 (cat. no. 950-124) at 95°C for 30 min. Endogenous peroxidase activity was subsequently blocked with UV INHIBITOR for 4 min at 37°C. The antibodies used were intended for in vitro diagnostic use, were ready to use and it was not necessary to perform any dilution. The samples were then incubated with Ventana primary antibodies specific to the target antigens: CONFIRM ER (SP1) (cat. no. 790-4324), CONFIRM PR (1E2) (cat. no. 790-2223) and PATHWAY anti-HER2 (4B5) (cat. no 790-4493), each incubated at 37°C for 20 min. Next, the samples were incubated with Amplification Kit (cat. no. 760-080; Roche Diagnostics) for 20 min at room temperature and treated with UV HRP UNIV MULT for 8 min at room temperature, followed by detection using UV DAB and UV DAB H2O2 for another 8 min. Additionally, they were incubated with UV COPPER for 4 min. Finally, the samples were counterstained with hematoxylin incubated at room temperature for 6 min, dehydrated and mounted in synthetic resin.

Cases were classified based on the expression of ER, PR and HER2 into the following subtypes: Luminal A (ER+, PR+ and HER2-); luminal B HER2+ (ER+, PR+ and HER2+); non-luminal HER2+ (Er-, PR- and HER2+); and triple-negative (ER-, PR- and HER2-; Fig. S1).

Image analysis

Immunohistochemistry images were scanned using Aperio equipment (Leica Microsystems, Inc.) and quantified using the ImageJ software (version 1.47; National Institutes of Health) following parameters described by Tuominen et al (33). The ImageJ plugin ‘ImmunoRatio’ (version 1.0c; 14.2.2011) and ‘ImmunoMembrane’ (version 1.0i; 8.7.2011) were used to analyze nuclear and membrane location, respectively, of the analyzed molecules. Receptor status was considered positive when ≥1% cells expressed the corresponding proteins, according to the American Society of Clinical Oncology/College of American Pathologists guideline recommendations (34,35). Chromogenic in situ hybridization (Ventana Her2 dual ISH DNA; cat. no. 760-6072) was performed on Her2+ samples with an equivocal result of >2. This suggested that the HER2 expression results obtained through immunohistochemistry were inconsistent as the presence of complete or incomplete membrane staining with weak to moderate intensity or complete and strong membrane staining was observed in <10% of the cells. Therefore, the evaluation of HER2 amplification using ISH evaluation was required. The HER2 detection process involves a cocktail of two DNA probes (cat. no. 800-4422; Ventana), with one targeting the HER2 gene, labeled with dinitrophenol, and another targeting the centromere of chromosome 17, labeled with digoxigenin (DIG). Probe hybridization were visualized using the Ventana UltraView SISH DNSP kit (cat. no. 800-098) and the Ventana UltraView Red ISH DIG detection kit (cat. no. 800-505).

Expression analysis

A transcriptome analysis of both luminal B HER2+ (ER+, PR+ and HER2+; n=15) and triple-negative (n=15) BC samples was performed to analyze critical molecular differences between aggressive tumors and improve understanding of the tumor biology of TNBC. RNA was extracted using TRIzol® Reagent (cat. no. 15596026; Invitrogen; Thermo Fisher Scientific, Inc.). The RT kit used was the GeneChip® Human Transcriptome Array 2.0 (cat. no. 902310; Applied Biosystems; Thermo Fisher Scientific, Inc.). A total of 30/197 samples were chosen for subsequent expression analysis as they contained sufficient tumor tissue, the appropriate molecular subtype [luminal B HER2+ (ER+, PR+ and HER2+) and TNBC (ER-, PR- and HER2-)] and high RNA quality (absorbance value at 260/280 of 1.9–2.2 and RIN >6). This analysis was conducted using the Human Transcriptome Array 2.0 platform (cat. no. 902162) developed by Affymetrix (Thermo Fisher Scientific, Inc.), with the capacity to identify a total of 44,699 transcripts associated with coding genes, as well as 22,829 transcripts corresponding to non-coding genes.

Transcriptome Analysis Console (version 4.0.1; Applied Biosystems; Thermo Fisher Scientific, Inc.) was subsequently used to identify differentially expressed (DE) genes (DEGs) between TNBC and luminal samples. Quality control and normalization were conducted using Robust Multichip Average (36). Transcripts were considered as DEGs if their fold change (FC) was found to be >2 or <-2 with P-values and a false discovery rate (FDR) <0.05. Transcript annotation was verified with the ‘BioMart-Ensembl’ tool (version GRCh38p.14; http://www.ensembl.org/biomart/martview/4cb2962c2a96efbcf3aa813eb2716e42). Raw data was deposited into gene expression omnibus (GEO) database with submission number GSE270721.

To visualize expression levels of DEGs in Mexican patients with BC, clustering analysis using ‘clustVis’ (https://biit.cs.ut.ee/clustvis/) (37), a variant of the R heatmap package (version 0.7.7). Heatmaps were generated for both coding and non-coding transcripts, comparing samples between luminal B, HER2+ and TNBC groups. Sample grouping was conducted using an unsupervised approach based on Euclidean correlation distances (38).

In silico validation and survival analysis using public databases

To validate DEGs found between TNBC and luminal B HER2+ BC from Mexican patients, data from a BC cohort from TCGA was used. This was carried out through the Genomic Data Commons Data Portal, using the term ‘breast cancer’ for the data search. Expression levels of mRNA transcripts in both TNBC and luminal B HER2+ BC were analyzed using the UCSC Xena browser platform, (version 2.0; http://xena.ucsc.edu/). Additionally, expression of lncRNA from the same tumor samples was evaluated using ‘The Atlas of Noncoding RNAs in Cancer’ (version: 2.2.1; http://www.tanric.org/), an open platform for exploring lncRNAs associated with TCGA data. LncRNAs not validated in TCGA cohort were considered exclusive to the Mexican patients of the present study. In addition, the 14 lncRNAs exclusive to Mexican patients were validated in a second cohort of Mexican patients from the GEO database (accession no. GSE134359) (39). To compare the present data with the GSE134359 study, transcripts were considered DE if they showed a FC >2 or FC <-2, with P-values and FDR values <0.05.

Overall survival of patients with BC

The impact of lncRNA expression levels (low or high) on the overall survival of breast cancer patients was evaluated. Furthermore, the impact of lncRNA expression on overall survival of patients was categorized according to molecular subtype (luminal A. luminal B and basal subtypes), through the use of tanric (40). The overall survival analysis was conducted by grouping patients according to lncRNA expression into low or high categories. For this purpose, the samples were divided using the median as a cutoff point. Subsequently, cox regression and log-rank P-values were calculated to assess the differences between the groups.

Key pathway analysis

DEGs were used to infer biological processes enriched in TNBC tumors. Ingenuity pathway analysis (IPA) software (version 122103623; Qiagen GmbH) was used to identify enriched pathways and to infer upstream regulators. Ensemble gene IDs, FC and FDR values were used as input parameters for IPA analysis. Results with z-scores ≥2 and a P<0.05 were considered statistically significant.

In silico exploration of the long intergenic non-coding RNA (LINC)01087 interactome

To explore the interaction of LINC01087 with other molecules, the online RNA interactome database from RNA inter (version 3.0; http://www.rnainter.org/) was used, which integrates RNA interaction data, including RNA-RNA, RNA-protein, RNA-DNA and RNA-histone modification. This method relies on confidence scores, which are derived from three different levels of supporting evidence: Interactions backed by strong experimental evidence, those supported by weak experimental evidence and predicted evidence. Each interacting molecule is assigned a confidence score, meaning that molecules with higher confidence scores are more likely to interact with LINC01087. RNA inter integrates information from literature and databases 10.1093/nar/gkab997, such as LncTarD (https://lnctard.bio-database.com/), NPInter (http://bigdata.ibp.ac.cn/npinter5), NoncoRNA (http://www.ncdtcdb.cn:8080/NoncoRNA/), miRDB (https://mirdb.org/), oRNAment (https://rnabiology.ircm.qc.ca/oRNAment) and tRFtarget (http://trftarget.net/). An overrepresentation analysis of the top molecules interacting with LINC01087 was also conducted with g:profiler software (version e112_eg59_p19_25aa4782; ELIXIR infrastructure). GProfiler performs gene set enrichment analysis interrogating DEGs using Gene Ontology, KEGG, Reactome and WikiPathways pathways; miRNA targets from miRTarBase and regulatory motif matches from TRANSFAC.

Cell culture

To validate the clinical findings using an experimental model, the present study examined whether some lncRNAs, previously identified in tumor biopsies, were differentially expressed and similarly regulated in a controlled environment, such as cell lines. Luminal (MCF7, ZR75 and T47D) and basal (MDA-MB-231, MDA-MB-468 and BT20) cell lines were acquired from American Type Culture Collection (ATCC). MCF7, T47D, ZR75 and BT20 cell lines were cultured in RPM1 medium (cat. no. 10-040-CV; Corning Inc) and MDA-MB-231 and MDA-MB-468 cell lines were cultured in DMEM medium (cat. no. 10-013-CV; Corning Inc). Cell cultures were supplemented with 5% FBS from ATCC (cat. no. 30-2020) and maintained at 37°C with 5% CO2.

RNA extraction

Total RNA was obtained and purified from tissue sections embedded in paraffin and fixed with formalin for 24 h at room temperature with the RNeasy FFPE Kit (cat. no. 73504; Qiagen GmbH), following the manufacturer's instructions. Cellularity of tumors was analyzed unblinded by a pathologist and only samples with ≥70% tumor cells were selected for subsequent experiments. RNA was quantified using spectrophotometry on the EPOCH system (Norgen BioTek Corp) and PicoGreen (cat. no. DNAQF; Sigma-Aldrich; Merck KGaA). Samples with high RNA quantity and purity were processed using the Sensation Plus kit (Affymetrix; Thermo Fisher Scientific, Inc.). After RNA extraction, samples underwent DNase digestion using the RQ1 RNase-Free DNase kit (cat. no. M6101; Promega Corporation). A single unit of DNase was used to degrade 1 µg of DNA in 10 min at 37°C.

RT-qPCR

RNA was reverse transcribed into cDNA using the High-Capacity cDNA Reverse Transcription Kit (Applied Biosystems; Thermo Fisher Scientific, Inc.), following the manufacturer's instructions. The following thermocycling conditions were used: 25°C for 10 min, 50°C for 120 min and 85°C for 5 min. Finally, qPCR was performed using the SYBR-select Master Mix kit (cat. no. 4472908; Applied Biosystems; Thermo Fisher Scientific, Inc.) and the LightCycler 480 system [Roche Diagnostics (Shanghai) Co., Ltd.]. The expression of several lncRNAs that showed differential expression were subsequently analysed. The primers used for RT-qPCR were as follows: LINC01087 forward, 5′-CGGTCTTTGTCATCGAGGCA-3′ and reverse, 5′-GGCAAAAGGATGGCTTGGAC-3′; Lnc-peroxidasin (Lnc-PXDN)-3:1 forward, 5′-CAGCGTCTCGTTCACCTCTT-3′ and reverse, 5′-GCTCCCTGAAGCACTGACAT-3′; SOX9-AS1 forward, 5′-AGCTGTCCCATGAGTGAAGC-3′ and reverse, 5′-CACTGGATGTCAAGGCTGGT-3′; Lnc-SYDE forward, 5′-AGAGAGGCTAGGTCGTCAGA-3′ and reverse, 5′-ACATGAGGCTGCTTAGTGACA-3′ and GAPDH forward, 5′-AGCCACATCGCTCAGACAC-3′ and reverse, 5′-GCCCAATACGACCAAATCC-3′.

The PCR amplification conditions for lncRNAs were: 50°C for 2 min, 95°C for 2 min, followed by 40 cycles of 95°C for 15 sec, 60°C for 15 sec and 60°C for 1 min. For GAPDH, the same initial conditions were used, followed by 40 cycles of 95°C for 15 sec, 58°C for 15 sec and 60°C for 1 min.

The primers designed aligned 100% with the corresponding target sequences. RT-PCR was performed and a single product was verified using melting curves. Product size was verified by gel electrophoresis. The expression level of each gene was normalized to GAPDH and calculated using the 2−ΔΔCq method (41).

Statistical analysis

Statistical analyses were performed in GraphPad Prism version 10.3.1 (Dotmatics). Before conducting the ANOVA, the normality of the data within each group was assessed using the Shapiro-Wilk test. To calculate statistical significance between two groups, an unpaired two-tailed student's t-test with Welch correction was used in case both groups had unequal variances. For analysis involving two groups or more, a one-way ANOVA parametric or non-parametric test was performed. and the mean rank of each column was compared with the mean rank of every other column. The normality of the data within each group was assessed using the Shapiro-Wilk test. If data were normally distributed, a parametric one-way ANOVA was used and the Brown-Forsythe test applied. If the normality assumption was violated, a non-parametric alternative ANOVA was used, followed by the Kruskal-Wallis test. Data were presented as mean ± SEM.

Results

Clinical and pathological characteristics

In the present study, of the initial 192 samples selected, 98.4% (n=189) of all analyzed cases were found in women. Age range varied from 25 to 88 years, with a mean of 57±12.6 years. Notably, 34.4% (n=66) patients studied were aged <50 years. In total, 97.9% (n=188) of the cases were pathologically diagnosed as invasive ductal carcinoma, whereas 73% (n=140) were classified as high-grade histological cases, based on standardized histological criteria (Scarff-Bloom Richardson scheme) (42,43) that evaluates cellular characteristics such as the degree of cellular differentiation, mitotic index and irregularities in the size and shape of the nuclei of tumor cells. Analysis of expression biomarkers revealed that 58% (n=111) exhibited luminal A characteristics, 21% (n=40) were identified as luminal B, 7% (n=13) were non-luminal HER2+ and 14% (n=26) were classified as TNBC.

TNBC exhibits different expression profiles to luminal B HER2+ tumors

To improve understanding of the biological characteristics of TNBC, 15 samples were selected and their transcriptomic profiles (mRNA and lncRNA) were compared with those in luminal B HER2+ tumors (n=15), which have one of the most favorable prognosis among aggressive breast tumors (44,45). Most of the selected samples corresponded to stages ranging from IIB to IIIC, with an average age of 58.4±12.6 years for the luminal B HER2+ type and 52.8±15.4 years for TNBC samples (Table SI).

The transcriptional portrait of TNBCs compared with luminal B HER2+ tissues revealed a total of 746 DEGs. DE RNAs were defined as transcripts with a fold change >2 or <-2 harboring P-values and an FDR <0.05. Among these DE transcripts, 555 were protein-coding (mRNAs), whilst 191 were non-coding RNAs (ncRNAs). The DE ncRNAs were then classified into 102 long non-coding RNAs (lncRNAs), 16 miRNAs, 36 small nucleolar RNAs (snoRNAs), 29 pseudogenes and 8 rRNAs (Fig. 1A). Amongst the DE mRNAs, 139 show increased expression, whereas 416 presented decreased expression in the TNBC samples, compared with those in luminal B HER2+ samples. Analysis of the DE lncRNAs revealed 27 upregulated DEGs and 75 downregulated DEGs (Fig. 1B).

Figure 1. Transcriptional changes in coding genes of patients with TNBC vs. luminal B HER2+. (A) Classification of DE transcripts in tissue samples from TNBC vs. luminal B HER2+ samples. (B) Bar graph representing the number of DE protein coding transcripts in TNBC vs. luminal B HER2+ samples. (C) Heatmap displaying the gene expression of coding transcripts between TNBC and luminal B HER2+ tissues. Columns represent expression of each sample. Range of gene expression spans from red (increased) to blue (decreased). (D) Volcano plot shows the classification of transcripts based on their expression and significance from patients with TNBC vs. luminal B HER2+ (left) and top DE mRNAs (right). Sankey diagram with a moon plot based on interpretative phenomenological analysis, showing significantly enriched (E) activated and (F) inhibited signaling pathways in DE coding and non-coding transcripts from samples of Mexican patients with TNBC. Moon plot includes z-score (left), P-value (right) and gene ratio. DE, differentially expressed; miRNA, microRNA; lncRNA, long non-coding RNA; TNBC, triple-negative breast cancer; HER2, human epidermal growth factor receptor 2.

Unsupervised hierarchical clustering analysis revealed that the samples segregated into two distinct groups, luminal B HER2+ and TNBC. Heatmap illustrating expression patterns of protein coding DE transcripts are depicted in Fig. 1C. The expression profile of lncRNAs provided a clearer distinction in transcriptional profiles between the two groups compared with the coding RNA profile.

The mRNAs with the highest change were GABRP (FC=28.19), ELF5 (FC=9.1), MUC5B (FC=14.34), CALML5 (FC=8.93), PROM1 (FC=5.52), SCGB2A2 (FC=−345.93), ERBB2 (FC=−81.24), PIP (FC=−78.39), GATA3 (FC=71.2), ANKRD30A (FC=−39.66), CDK12 (FC=−28.95) and TBC1D9 (FC=−26.19; Fig. 1D). A total of 57 DE transcripts are involved in processes that contribute to invasion and metastasis of BC, by promoting cell migration and remodeling of the extracellular matrix (data not shown).

A pathway enrichment analysis was also conducted on the DE transcripts to compare TNBC with luminal B HER2+ tumors using the IPA software. Analysis revealed activation of several signaling pathways, such as ‘Wnt/β-catenin’, ‘Rho GDP Dissociation Inhibitor (RHOGDI)’, ‘MAPK’ and ‘NF-κB’ (Fig. 1E). These pathways have been previously described to be the signaling pathways that can promote proliferation, migration and invasion of triple negative breast cancer cell lines, which are characteristics of TNBC tumors (46,47). Additionally, pathways and processes were identified that were considerably inhibited in TNBC, such as ‘extracellular matrix organization’, ‘integrin cell surface interactions’, ‘estrogen receptor signaling’, ‘Interleukin-4 (IL-4) signaling’, ‘mesenchymal-epithelial transition’ and ‘homeobox antisense intergenic RNA regulatory pathways’. These processes and pathways revealed the nature of TNBC tumors. Unlike luminal tumors, which exhibit strong adhesion properties that influence their movement, TNBC tumors possess weaker adhesion traits. This difference affects how triple-negative cancer cells migrate, interact with tumor microenvironment and metastasize (Fig. 1F).

To confirm the gene expression results retrieved from the microarrays, expression changes of DE mRNAs were next analyzed on the UCSC Xena platform. This software integrates accurate and consistent data sources and clinical information, including TCGA, information that has been validated by various independent research teams (48). A cohort of 758 TCGA BC samples was used, where only samples with a basal or luminal molecular profile were considered. Expression of DEGs was evaluated in basal (n=142), luminal B (n=194) and luminal A (n=422) tumors. The data of the present study were successfully validated in luminal B samples and luminal A samples. Figs. 2 and 3 display the top 25 transcripts with the highest rate of change, regulated both negatively (Fig. 2A-C) and positively (Fig. 3A-C), in patients with basal BC compared with those with the luminal subtype. As expected, reduced expression of ESR1, PGR and ERBB2 was identified in basal tumors. The genes showing the highest fold change between basal and luminal samples were anterior gradient protein 3, ankyrin repeat domain-containing protein 30a (ANKRD30A) and arylamine n-acetyltransferase 1 (downregulated in basal tumors), whereas phosphoserine aminotransferase 1, prominin-1 (PROM1) and forkhead box c1 were upregulated in basal tumors.

Figure 2. Validation of inhibited mRNA expression profiles in basal vs. luminal TCGA samples. (A) Heatmap displaying genetic expression across different tissue subtypes of patients with breast cancer. Columns represent expression of each validated transcript in TCGA samples, whilst rows indicate each evaluated sample. Genetic expression is visualized on a gradient ranging from red (increased) to blue (decreased). (B) Expression levels of differentially expressed mRNAs among subtypes of breast cancer. (C) Expression levels of differentially expressed mRNAs among subtypes of breast cancer. Statistical analysis was conducted using a one-way ANOVA followed by the Brown-Forsythe test for non-parametric data or Kruskal Wallis test for parametric data (****P<0.0001, ***P<0.001 and *P<0.05). TCGA, The Cancer Genome Atlas.

Figure 3. Validation of highly expressed mRNA expression profiles in basal vs. luminal TCGA samples. (A) Heatmap displaying genetic expression across different tissue subtypes of patients with breast cancer. Columns represent the expression of each validated transcript in TCGA samples, whilst rows indicate each evaluated sample. Genetic expression is visualized on a gradient ranging from red (increased) to blue (decreased). (B) Expression levels of differentially expressed mRNAs among subtypes of breast cancer. (C) Expression levels of differentially expressed mRNAs among subtypes of breast cancer. Statistical analysis was conducted using a one-way ANOVA followed by the Brown-Forsythe test for non-parametric data or Kruskal Wallis test for parametric data (****P<0.0001 and ***P<0.001). TCGA, The Cancer Genome Atlas.

Expression profiles of lncRNAs in TNBC and luminal subtypes of BC

Heatmap illustrating expression patterns of non-coding DE transcripts are depicted in Fig. 4A. LncRNAs exhibiting the most significant changes between TNBC and luminal B HER2+ were as follows: AL157387.1 (FC=−39.39), AC093001.1 (FC=−28.74), LINC01087 (FC=−13.22), AC044784.1 (FC=−12.83), SOX9-AS1 (FC=12.9), mucin 5B-AS1 (FC=8.45), AC092168.2 (FC=5.46) and small nucleolar RNA host gene 16 (FC=4.97). Although some of these molecules have been described in certain types of tumors, they have not been fully characterized and their functions therefore remain undefined (Fig. 4B). For example, it has been reported that AC093001.1 may serve an important role in clear cell renal cell carcinoma (49) and the lncRNA Sox9-AS1 may be important in hepatocellular carcinoma and breast cancer (50). LINC01087 has been reported to increase the aggressiveness of breast cancer and to predict patient outcome, however the mechanism of action of this LINC01087 is poorly understood (51), Mucin 5B-AS1 has been involved in promoting cell migration and invasion processes in lung cancer cells through interaction with MUC5B (52).

Figure 4. Analysis and validation of lncRNAs in patients with triple-negative breast cancer versus luminal B (HER2+) using Mexican and TCGA datasets. (A) Heatmap showing DE lncRNAs in TNBC tissue samples vs. luminal B HER2+ samples. Range of gene expression spans from red (increased) to blue (decreased). (B) Volcano plots and top DE lncRNAs. (C) Overexpressed LncRNAs in basal tumors vs. luminal B (HER2+). Each point on the graph represents a patient, either with Basal (orange points) or luminal B (HER2+) tumors (green points). (D) Inhibited LncRNAs in basal tumors compared with luminal B (HER2+). (E) LncRNA signature found only in the cohort of Mexican patients. Changes in expression levels in TNBC samples are shown in blue (decreased) and red (increased), whereas green rectangles highlight lncRNAs validated in a second cohort of Mexican patients (accession no. GSE134359). Additionally, the font size represents the magnitude of the change in expression. Larger font sizes indicate more pronounced changes, while smaller font sizes correspond to more subtle changes. (F) Overall survival of LINC00174 in TNBC tumors using TCGA data. ****P<0.0001, ***P<0.001, **P<0.01 and *P<0.05. LINC, long intergenic non-coding RNA; LncRNA, long non-coding RNA; TNBC, with triple-negative breast cancer; HR, hazard ratio; DE, differentially expressed; FDR, false discovery rate; HER2, human epidermal growth factor receptor 2.

Next, the expression of DEGs in patients with luminal BC and TNBC tumors from TCGA were assessed. Amongst the 102 lncRNAs, only 48 well-annotated lncRNAs had information available. Of these lncRNAs, 40/48 (83.3%) were successfully validated with statistical significance. Representative validated lncRNAs in Fig. 4C and D were those exhibiting the most pronounced change. The integration of in silico and experimental analysis enabled the identification of potential lncRNAs associated with basal BC.

These DE lncRNAs were validated using TCGA samples originating from non-Mexican patients. The majority of these differences were made to emphasize the discrepancies found in the analysis of the present study of Mexican patients compared with TCGA data. Additionally, a set of 14 lncRNAs were identified which exhibited DE patterns in the cohort of Mexican patients that were not identified in the European TCGA cohorts (Fig. 4E). Additionally, a second cohort of Mexican patients (accession no. GSE134359) was analyzed in the same manner as the DEGs in basal and luminal B (HER2+) samples (cut-off criteria: FC=>2 or <-2/P-values and FDR <0.05). Fig. 4E revealed that AC011676.1, Lnc-lactalbumin α-1, Lnc-zinc finger protein 200-1 and LINC00174 were successfully validated, meaning the same expression pattern (increased or decreased according to the established criteria) was observed in this second cohort. Although the remaining 10 lncRNAs showed a similar trend, they failed to meet established cut-off criteria.

To assess the influence of these lncRNAs in the prognosis of patients with BC, TCGA data was used to determine whether any of the lncRNAs validated in the second cohort were associated with overall survival. The findings revealed that the increased expression of LINC00174 was associated with worse overall survival in patients harboring basal breast tumors. However, this association was not statistically significant (Fig. 4F). High and low LINC00174 expression was determined by dividing the expression data into percentiles.

Additionally, four lncRNAs were selected and their expression was analyzed in luminal and TNBC cell lines. To validate the clinical findings using an experimental model, the expression of the lncRNAs was examined to determine whether they are expressed and upregulated or downregulated in a controlled environment, such as in cell lines, in a similar manner to how they behave in previously identified tumor biopsies from Mexican patients in the present study. Fig. 5A, D, G and J revealed the expression levels of DE lncRNAs between luminal B HER2+ and TNBC subtypes derived from Mexican patients with BC. These data originated from the expression analysis using Human Transcriptome Array 2.0 microarrays. Subsequently, Fig. 5B, E, H and K illustrate the expression levels of these lncRNAs in a panel of BC cell lines belonging to luminal (MCF7, ZR75 and T47D) and triple-negative (MDA-MB-231 and MDA-MB-468, BT20) subtypes. Specifically, SOXA9_AS1 exhibited increased expression levels in TNBC cell lines compared with those in the luminal cell lines, whilst Lnc-PXDN-3:1, Lnc-SYDE and LINC01087 exhibited decreased expression levels in TNBC cell lines. These findings are consistent with those found by the microarray analysis. Expression of downregulated lncRNA SYDE was found to be reduced in TNBC cell lines, but the results were not statistically significant (Fig. 5H).

Figure 5. Analysis of lncRNA expression in cell lines across several breast cancer subtypes. Relative expression of the (A) lncRNA SOX9_AS, based on microarray data obtained from Mexican patients with BC, (B) between luminal and TNBC cell lines evaluated by RT-qPCR and (C) in luminal vs. basal, Claudin-low and HER2 breast cancer cell lines evaluated by RT-qPCR. (D) Relative expression of the lncRNA PXDN-3:1, based on microarray data obtained from Mexican patients with BC, (E) between luminal and TNBC cell lines evaluated by RT-qPCR and (F) in luminal vs. basal, Claudin-low and HER2 breast cancer cell lines evaluated by RT-qPCR. (G) Relative expression of the lncRNA SYDE, based on microarray data obtained from Mexican patients with BC, (H) between luminal and TNBC cell lines evaluated by RT-qPCR and (I) in luminal vs. basal, Claudin-low and HER2 breast cancer cell lines evaluated by RT-qPCR. (J) Relative expression of the LINC01087, based on microarray data obtained from Mexican patients with BC, (K) between luminal and TNBC cell lines evaluated by RT-qPCR and (L) in luminal vs. basal, Claudin-low and HER2 breast cancer cell lines evaluated by RT-qPCR.. Each data point on the graph represents a biological replicate, which is composed of the mean of two technical replicates. Statistical analysis was conducted using two groups with an unpaired student's t-test followed by the Welch correction in case of unequal variances or a one-way ANOVA followed by the Brown-Forsythe test for non-parametric data. ****P<0.0001, ***P<0.001, **P<0.01 and *P<0.05. lncRNA or lnc, long non-coding RNA; LINC, long intergenic non-coding RNA; PXDN, peroxidasin; TNBC, triple-negative breast cancer; RT-qPCR, reverse transcription-quantitative PCR.

Since TNBC can harbor both basal and claudin-low molecular subtypes (53), cell lines were next grouped according to their molecular subtypes, namely Basal (MDA-MB-468 and BT20), claudin-low (MDA-MB-231), luminal (MCF7, T47D and ZR75) and HER2 (SKBR3 and MDA-MB-453) (Fig. 5C, F, I and L). Expression of lncRNA SOX9-AS was revealed to be increased in basal and claudin-low cell lines compared with that in luminal and HER2 cell lines (Fig. 5C). LINC01087 and SYDE exhibited different expression patterns in basal and claudin-low cell lines, with elevated levels observed particularly in basal cell lines. In the case of lncPXDN, significant differences were identified when comparing luminal cell lines with basal and claudin-low cell lines. Furthermore, HER2 cell lines also showed significant differences compared with claudin-low and basal cell lines (Fig. 5F). Specifically, LINC01087 exhibited decreased expression in basal, but it was absent in claudin-low cells (Fig. 5L), suggesting that it may only be expressed in the basal molecular subtype of TNBC. Therefore, expression of LINC01087 could be used to differentiate between basal and claudin-low subtypes. Furthermore, SYDE was expressed at increased levels in luminal cell lines compared with the other cells. The expression of SYDE was significantly decreased in claudin-low and HER2 cells compared with luminal cels. Basal cell lines had a slightly decreased expression of SYDE compared with luminal cells, although this was not statistically significant (Fig. 5I).

DE lncRNAs were also evaluated in an independent validation cohort of 21 samples obtained from formalin-fixed paraffin-embedded tumors from Mexican patients. Samples were divided into TNBC (n=8) and luminal B HER2+ tumors (n=13). Results revealed a significant decrease in the expression of lncRNAs PXDN-3:1 and LINC01087 in TNBC samples compared with those in luminal B HER2+ tumors. Furthermore, a decrease in the expression of lncRNA SYDE was observed in TNBC tumors, whilst that of lncRNA-SOX9-AS1 was also increased, though no statistical significance was reached (Fig. 6).

Figure 6. Analysis of lncRNA expression in an independent cohort of triple-negative breast cancer vs. luminal breast cancer tissues. The lncRNAs were evaluated in an independent validation cohort of 21 samples obtained from FFPE tumors from Mexican patients. The graphs show expression of lncRNAs using reverse transcription-quantitative PCR. (A) lncRNA Sox9-AS, (B) lnc-PXDN-3:1, (C) lnc-SYDE and (D) LINC01087. Each point on the graph represents a sample from an independent patient. Statistical analysis was conducted using two groups with an unpaired student's t-test followed by the Welch correction. **P<0.01, ***P<0.001. LncRNA, long non-coding RNA; LINC, long intergenic non-coding RNA.long intergenic non-coding RNA; PXDN, peroxidasin; SYDE, synapse defective rho gtpase homolog.

In summary, the validation graphs depicting the expression of lncRNAs across the different cell lines are consistent with findings from the Mexican patient samples and TCGA data. These findings suggest that the differential expression patterns of identified lncRNAs can serve as a key biomarker for the molecular classification of BC, while possibly serving a fundamental functional role in TNBC progression.

LncRNAs with clinical significance in patients with BC

Clinical implication of lncRNAs in the survival of patients with BC holds importance (54). A previous study has demonstrated that certain lncRNAs may be associated with diagnosis, prediction or treatment response in patients with BC (54). Therefore, overall survival of DE lncRNAs was next assessed using TCGA BC data.

The relationship between 102 DE lncRNAs and overall survival were evaluated using TCGA BC data. Validation of lncRNAs in both cell lines and a second set of tissue samples different to the samples used for microarray, revealed that the decreased expression of genes, such as LINC01087, LINC02568, ACO22196.1 and eosinophil granule ontogeny transcript (EGOT), are associated with poor overall survival in samples of BC, regardless of molecular subtype (Fig. 7). Furthermore, overall survival analysis revealed a significant impact of LINC01087 on the survival of patients with the luminal B subtype of BC (Fig. 8C), where the decreased expression of this lncRNA was found to be associated with worse prognosis. No significant associations were found between overall survival of patients with BC and the expression of lncRNAs SOX9-AS1, lncRNA-PXDN-3:1 and SYDE (data not shown). When the samples were classified by molecular subtype, lncRNA-PXDN-3:1 revealed a significant association in ER+ samples, where worse survival was associated with the high expression of lncRNA-PXDN-3:1 (Fig. 8A).

Figure 7. LncRNAs with clinical significance in patients with breast cancer. Kaplan-Meier plots from TCGA breast cancer data demonstrate that lncRNAs (A) LINC01087, (B) LINC02568, (C) ACO22196.1 and (D) EGOT may be associated with overall survival in patients with breast cancer, regardless of molecular classification. The green line on Kaplan-Meier plots indicates high expression of lncRNA, whilst the brown line indicates low expression. LncRNA, long non-coding RNA; LINC, long intergenic non-coding RNA; EGOT, eosinophil granule ontogeny transcript.

Figure 8. LncRNAs involved in the survival of patients with breast cancer with specific molecular subtypes. Kaplan-Meier plots displaying lncRNAs involved in the overall survival of patients with (A) ER+, specifically the PXDN-3.1 (B) Luminal A, specifically the Lnc-GATA-3-1 (C) Luminal B, the LINC01087 (D) Luminal B, the Lnc- BX679671.1 (E) Luminal B, the Lnc-FSIP1.6:3 (F) Luminal B, the LINC00182 (G) basal, the LINC00504 and (H) basal, Lnc-ARHGEF38-IT1. Blue lines indicate low expression of lncRNA, whilst red lines indicate high expression. Blue lines indicate low expression of lncRNA, whilst red lines indicate high expression. LA, luminal A, LB, luminal B; LncRNA, long non-coding RNA; LINC, long intergenic non-coding RNA; PXDN, peroxidasin; GATA-3-1, GATA binding protein 3-1; FSIP1, fibrous sheath interacting protein 1; ARHGEF38-IT1, rho guanine nucleotide exchange factor 38 intronic transcript 1.

In addition, the overall survival of patients with distinct molecular subtypes of BC according to the low or high expression of other lncRNAs was assessed. In luminal A samples, low expression of lncRNA-GATA-3-1 was found to associate with worse overall survival (Fig. 8B). In luminal B samples, low expression of lncRNA BX679671.1 was associated with worse overall survival (Fig. 8D), whilst high expression of lncRNAs fibrous sheath interacting protein 1-6:3 and LINC00182 was associated with reduced survival (Fig. 8E and F). In basal samples, low expression of lncRNAs LINC00504 and ARHGEF38-IT1 were associated with worse survival. Notably, these lncRNAs exhibited decreased expression in TNBC samples compared with luminal B HER2 samples, consistent with the aggressive characteristics of TNBC (Fig. 8G and H).

LINC01087 interacts with histones and transcription factors to possibly regulate ESR-mediated signaling

In silico analysis was next conducted to characterize SOX9-AS1, PXDN-3.1 and LINC01087 using the RNA interactome repository (55), to investigate potential interactions of this lncRNA with other molecules (Fig. 9A-C). Analysis revealed that the majority of interactions with LINC01087 strongly associated with chromatin marks, such as H3K27me3, H3K4me3, H327ac and H3K4me1/2 (Fig. 9C). Additionally, interactions of LINC01087 with several transcription factors were identified, including CCCTC-binding factor, estrogen receptor 1, Forkhead box (FOX) A1, SOX2, STAT3 and Krüppel-like factor. Additionally, analysis of LINC01087 revealed interactions with other lncRNAs, such as non-coding RNA activated by DNA damage and AP000526.1 (Fig. 9C). Analysis of co-regulated gene patterns to explore any possible functional relationships is shown according to gene color and size (Fig. 9D).

Figure 9. Prediction of SOX-AS1, PXDN-3.1 and LINC01087 interactions with associated biological processes and signaling pathways. Circus plot indicating interactions of (A) SOX-AS1, (B) PXDN-3.1 and (C) LINC01087 with molecules, such as histone marks (light green dots), transcription factors (yellow dots), lncRNAs (dark blue dots), G protein-coupled proteins (red dots), chromatin remodeling factors (light blue dots), RNA-binding proteins (dark green dots). (D) Gene co-regulation analysis, with the font size reflecting the degree of correlation between molecules and the color reflects the degree of expression, with red as the highest degree of expression. (E) Processes and signaling pathways enriched by the presence of LINC1087. LINC, long intergenic non-coding RNA.

Since LINC01087 has been associated with aggressiveness in several types of cancers, including breast cancer (56–59), the possible function of this lncRNA was explored by investigating the involvement of biological processes and signaling pathways associated with LINC01087. Functional enrichment analysis successfully identified relevant signaling pathways where LINC01087 may participate. Results revealed engagement of this transcript in various biological functions, such as ‘signaling by nuclear receptors’, ‘ESR-mediated signal’, ‘transcription factor binding’, ‘estrogen-dependent gene expression’, ‘chromatin binding’ and ‘protein dimerization activity’, suggesting that LINC01087 is involved in nuclear receptor signaling. These analyses suggest that molecules interacting with LINC01087 participate in ER-mediated signaling, which suggests that this lncRNA may regulate the expression of ER-related genes and have a role in luminal tumors (Fig. 9E).

Discussion

BC is a notably heterogeneous disease (60) and the TNBC subtype is particularly distinguished for its aggressive nature (44). Since the different molecular subtypes of breast cancer have distinct clinical behaviors and specific therapies, research has been focused on the identification of biomarkers to stratify patients both therapeutically and prognostically (15). Results of the data generated in the present study are comparable with other studies conducted worldwide, due to the expression profiles of these tumors and the signaling pathways and processes dysregulated in these tumors (61,62). However, the observed frequency of TNBC was lower when compared with that reported in previous studies in the Mexican population (31,32). This discrepancy may be attributed to the cut-off points in the present study used to positively diagnose TNBC. The official cutoff value was used to assess the expression of estrogen and progesterone receptors in BC using IHC, following American Society of Clinical Oncology/College of American Pathologists guidelines (35). A result was considered positive if ≥1% of tumor cell nuclei showed staining. The transcriptomic analysis, using cutoff values based on fold change and statistical analysis, confirmed the expression findings. No hormonal receptors were observed in TNBC, whereas they were expressed in luminal tumors.

In TNBC, lncRNAs have emerged as promising markers for diagnosis, prognosis, treatment resistance and monitoring the spread of cancer. These lncRNAs interact with a variety of molecules, such as DNA, RNA and proteins, to activate the metastatic transcriptional network (18). Although luminal subtypes are generally characterized by a more favorable prognosis in comparison with TNBC, it is important to highlight that the luminal B subgroup is also associated with a poor prognosis (45,54,63). Given that transcriptomic profiles and molecular classification panels for BC rely heavily on mRNA, the present study focused on the investigation of lncRNAs.

The present study aimed to elucidate differences in the gene expression levels of lncRNAs and mRNAs among tissue samples between patients with TNBC and luminal B HER2+ BC. To ensure robustness, reliability and accurate interpretation of the microarray data, gene expression changes were validated in patient samples using TCGA database. Results confirmed the consistency of both coding and non-coding RNA data. In Fig. 4C and D, the representative lncRNAs, that is, those with the most pronounced and significant changes between luminal B and TNBC samples, were shown. Thus, the graphs only show 19 lncRNAs, chosen for their higher statistical robustness and clearer representation of the observed trends.

Regarding AC011676, no significant difference was found between Luminal B and TNBC Samples. Since only 26 lncRNAs were upregulated, information was available for only 8/26 lncRNAs in the TCGA cohort. As a result, AC011676 was one of the 8 lncRNAS that failed validation.

Expression of lncRNAs in TNBC has been previously described in other studies, which have identified specific lncRNAs differentially expressed across the different molecular subtypes (54,63). However, some transcripts are difficult to identify due to their low expression levels and tissue-specific patterns. In particular, some transcripts may not be abundant enough in the sample to be reliably detected, especially in heterogeneous tissues such as breast tumors. Additionally, variability in expression patterns across different tissue types makes it difficult to distinguish signals from non-cancerous tissues within the sample. Therefore, the present study provided specific information from Mexican patients and contributed support for future clinical considerations. A total of 14 lncRNAs were identified that were not DE in TCGA cohorts, which mainly included data from European patients. Among these, 4 lncRNAs were validated in a different Mexican cohort, including TNBC and luminal B HER2+ tumors. However, since several lncRNAs remain unidentified and their functions remain poorly defined, the presented study integrated coding and non-coding data to investigate potential signaling pathways associated with the TNBC tumor phenotype.

Several signaling pathways were found to be activated in TNBC that are known for regulating various processes during tumor progression and metastasis (18)Specifically, the present study identified that the NF-κB, Wnt/β-catenin, MAPK and RHOGDI signaling pathways are activated in TNBC. NF-κB and Wnt/β-catenin are established signaling axes that serve key roles in cell proliferation, maintenance of stemness and development of drug resistance in TNBC (64). Their activation considerably affects the survival of patients with TNBC. MAPK pathways are involved in signal transduction cascades that promote tumor progression, invasion and metastasis of TNBC (65). Constitutive activation of MAPKs has been observed in aggressive TNBC, which is associated with worse overall survival (66). By contrast, RHOGDIs regulates the Rho family of GTPases, which control several signal transduction pathways, such as the actin cytoskeleton remodeling pathway and the cell cycle regulatory pathway. RHOGDI can both inhibit or activate Rho GTPases, potentially serving a dual role in cancer (67). Notably, RhoGDI-2 is overexpressed in TNBC and its inhibition can induce apoptosis whilst sensitizing BC cells to cisplatin treatment (46).

The present study revealed several processes that are inhibited, such as extracellular matrix organization, integrins and syndecan interactions, ESR and IL-4 signaling. The altered pathways reflected changes in adhesion properties of TNBCs. TNBC cells have fewer adhesive properties, which facilitates cell migration and invasion (18). In this regard, several clinical trials are currently underway to evaluate the efficacy of pharmacological agents that can regulate components of the extracellular matrix in TNBC (68), For example, Lucitanib, which inhibits angiogenesis and reduces matrix metalloproteases and collagen, or Reparixin, an interleukin-8 receptor C-X-C chemokine receptor type 1 and 2 inhibitor (69,70). LncRNAs have emerged as key components in BC biology, serving a key role in the regulation of gene expression and disruption of key biological pathways in the development and progression of BC (71,72). For example, the AKT/PI3K/mTOR pathway is mediated by lncRNAs such as HOX transcript antisense intergenic RNA (HOTAIR), metastasis associated lung adenocarcinoma transcript 1 and lnc-urothelial carcinoembryonic antigen 1. Additionally, HOTAIR and lncRNA maternally expressed gene 3 (lncRNA-MEG3) regulate NF-κB signaling (72) and lncRNA growth arrest-specific transcript 5 influences apoptosis through the miR-378a-5p/suppressor of fused homolog pathway in TNBC (73). In addition, the diversity and specificity of lncRNAs in different molecular subtypes of BC suggest potential applications in patient stratification and personalized therapies (25,74)

It is noteworthy that the transcription profile of lncRNAs is more consistent compared with that of mRNAs. This is reflected by the reduced variability of data among molecular subgroups, in addition to a marked difference in the diversity of molecular subtype expression. Congruence in validation results using RT-qPCR and expression microarray data serves to underscore the steadfastness and robustness of the findings. Reduction in the expression of PXDN-3:1 and LINC01087 in TNBC cell lines supports their potential role in characterizing this subtype. LncRNA PXDN-3:1 was also highly expressed in the luminal and HER2 molecular subtypes compared with that in the basal and claudin-low cell lines, where its expression was scarce. Similarly, increased expression of lncRNA-SOX9-AS1 in this subtype suggests its potential contribution to the molecular characteristics of TNBC.

Although SYDE results were not initially significant in TNBC cell lines, subgroup analysis highlighted its potential role in specific subtypes. This detailed approach could provide an improved reflection of biological variability in TNBC cell lines. Certain lncRNAs may exhibit such sensitivity in their expression levels to enable classification of subtypes, particularly within the TNBC subtype. This sensitivity can lead to more precise stratification in terms of classification and survival outcomes. Claudin-low samples have been observed to display a more aggressive behavior compared with basal samples, where this variability has been reported to be associated with different mutations that can impact drug response (75). In the present study, LncRNA SYDE exhibited reduced expression in claudin-low and HER2 cell lines compared with that in basal cell lines. This suggests the possibility that SYDE may be involved in proliferation processes, given that basal lines tend to display a higher proliferation rate compared with claudin-low lines (76). However, the functional role of lncRNA SYDE remains unknown. Therefore, additional studies are required to deepen the understanding of this lncRNA. In addition, a significant decrease in the expression of PXDN-3:1 and LINC1087 transcripts was observed in an independent set of FFPE TNBC tissues compared with luminal tissues. These results support the validity of these lncRNAs as potential biomarkers for characterizing tumor subtypes.

The decreasing trend in SYDE expression observed in TNBC tumors underscores the importance of using a larger sample size to detect more pronounced differences. Similarly, the increasing trend in lncRNA-SOX9-AS1 expression in TNBC tumors, whilst not reaching statistical significance in the present study, merits further exploration. Previous studies have also revealed that SOX9-AS1 is highly expressed in TNBC tumors compared with other BRCA subtypes (77,78). Among the various TNBC molecular subtypes, the overexpression of SOX9-AS1 is associated with a favorable prognosis (77). Furthermore, other studies suggest that SOX9-AS1 may be involved in migration, invasion and lipid metabolic reprogramming in TNBC cells (77,79). Previously, SOX9-AS1 was identified to regulate cellular senescence resistance in TNBC cells through the Wnt pathway (78). SOX9-AS1 knockdown promotes senescence and immune cell infiltration, suggesting its involvement in modulating the immune microenvironment (78). However, databases report different isoforms of SOX9-AS1 that are susceptible to more robust biological analysis to determine if they have any clinical implications.

Several studies have revealed the potential of the association of LINC01087 with carcinogenesis and suggest that it may possess significant diagnostic value (56–58). This lncRNA has been previously investigated in esophageal, gastric, ovarian and breast cancers (59,80). It has been reported that LINC01087, located in the cytoplasm, functions as an endogenous competitor (ceRNA) by binding to specific miRNAs to prevent interaction with their mRNA targets. In the case of gastric cancer (GC), it has been observed that LINC01087 modulates tumorigenesis through the miR-135a-5p/caspase activity and apoptosis inhibitor 1 (CAAP1) signaling pathway (56). LINC01087 is highly expressed in GC and acts as a molecular sponge for miR-135a-5p, allowing CAAP1 expression to increase the migration and invasion capacity of tumor cells (56). The analysis of LINC01087 was deepened owing to the growing evidence of its relevance in the context of BC (51,58,59,80). In TNBC, in silico experiments previously revealed that LINC01087 can form a ceRNA network with LINC01087/miR-135b-5p by co-expressing with sulfite oxidase, E6-AP carboxyl terminus domain E3 ubiquitin protein ligase 2 and solute carrier family 39 member 6 (79).

However, interaction analysis conducted in the present study suggested that besides interacting with miRNAs, LINC01087 may interact with nuclear molecules, such as transcription factors ER2, PR, Runt-related transcription factor 1 and GATA. This finding may provide novel insights into the function of this lncRNA in the nucleus and its role in BC. The data from our study support findings, indicating a decrease in LINC01087 in TNBC, in addition to upregulation in luminal subtypes (51). Therefore, LINC01087 may have a considerable role in the classification of luminal and TNBC samples. Furthermore, data of the present study indicated that increased LINC01087 expression is associated with a favorable prognosis in BC. Consequently, low expression levels of LINC01087 were associated with reduced survival, as seen in TNBC cases (51). No information was found regarding the lncRNA lnc-Syde, which could present a novel research opportunity.

In the present study, low expression of LINC01087, ACO22196.1 and EGOT is associated with a reduction in 5-year survival. Furthermore, it was demonstrated that the decreased expression of FOXP1-IT1, LINC02568 and Z9330.2 predicts worse 20-year survival, highlighting their impact on long-term prognosis. These data suggest that these lncRNAs may have a role in controlling tumor progression. Survival analysis specific to molecular subtypes provides additional information on the relationship between lncRNA expression and the nature of each molecular subtype, reinforcing the potential clinical relevance of these findings.

Association between low expression levels of several lncRNAs and worse overall survival, along with specific associations in molecular subtypes, emphasizes the significance of considering these lncRNAs as potential biomarkers for patient stratification and development of more personalized therapeutic approaches in BC. In this sense, it is relevant to classify and identify biological subtypes of clinical importance using different approaches, such as mutational profiles of mitochondrial DNA, in addition to both radiological and histopathological image analysis (81,82). Nevertheless, further studies are warranted. Despite the use of additional tissue samples to validate results, the sample size remains limited, calling for the expansion of analysis to larger sample sizes to confirm these findings and achieve more robust statistical significance, thereby substantiating and further exploring these prognostic associations.

The present study identified a set of lncRNAs that are DE in Mexican patients with TNBC and luminal BC. However, some limitations hindered the ability to obtain a unique profile of lncRNAs specifically associated with TNBC. Analysis focused on only two of the five molecular subtypes, meaning that confirmation on whether the DE lncRNAs reported in the present study are predominantly associated with TNBC remains unconfirmed. Another limitation is the small sample size that was analyzed in the microarray, which can be less representative and may not allow generalization to a broader population. A strength of the present study is the representative selection of tumor tissues in the presence or absence of immunohistochemical biomarkers for transcriptomic analysis. Another strength is that the transcriptomic analysis was performed in a group of tumors from Mexican patients, providing more specific information about transcriptomic signatures of clinical relevance in Mexico. Additionally, the present analysis identified transcriptomic data that, regardless of the population, were consistent with those of European tumors according to TCGA data. Likewise, the lncRNA signature enriched in the Mexican population was identified, which was validated in another independent Mexican cohort (GSE134359). Additionally, previous studies have compared transcriptomic profiles of TNBC and HER2-low or lncRNA expression profiles in luminal B associated with neoadjuvant chemoresistance (83,84). However, to the best of our knowledge, studies that analyzed the DEG profile between TNBC and luminal B HER2+ could not be identified. Therefore, the present study provided information to identify the signature of lncRNAs associated with prognosis.

Further genomic and epigenomic studies must be conducted on Latin American populations to achieve a more representative landscape, since exposure to environmental factors may vary between populations, in addition to genetic differences.

The present study identified transcripts associated with tumor aggressiveness, including lncRNAs linked to poorer prognosis in patients with TNBC. Furthermore, it identified molecules involved in signaling pathways that could drive the more aggressive phenotype of TNBC compared with the luminal B HER2+ subtype, contributing to the discovery of potentially useful biomarkers for patient prognosis. Since TNBC is the most aggressive subtype with limited therapeutic options, generating insights into its biology is key for improving clinical management (85).

There are studies comparing TNBC and luminal B HER2+ or triple-positive (TPBC) subtypes, analyzing clinicopathological differences, such as tumor stage, metastatic sites and lymphovascular invasion (LVI) (61,62). In a study by Mandor et al (61), these comparisons aimed to evaluate possible associations between these clinical characteristics and patient recurrence or survival. The results showed no significant differences between the two groups in terms of disease progression. However, significant changes in progression, such as LVI, were observed within each group separately (61).

A study previously conducted in Korea highlighted differences in the clinical courses of TNBC and TPBC. The authors found significant variations in histologic grade, nuclear grade, lymphatic invasion and long-term survival after treatment in patients with these BC subtypes (62). However, to the best of our knowledge, studies specifically focusing on transcriptomic evaluations that compared only these two molecular subtypes could not be identified, underscoring the relevance of the work in the present study.

To conclude, BC subtypes are highly heterogeneous, such that comparing the transcriptomes of each subtype may facilitate understanding into how molecular differences can influence their clinical behavior. The present study focused on comparing TNBC and luminal B HER2+ subtypes because whilst luminal B HER2+ is an aggressive subtype, it shows improved outcomes compared with TNBC (44,86) due to the availability of targeted therapies. By contrast, the TNBC subtype has a worse prognosis because of its aggressiveness and the lack of specific treatments.

Collectively, these discoveries underscore the intricate involvement of mRNAs and lncRNAs in the landscape of BC biology as a whole and in the context of each one of the molecular subtypes. The findings reveal that specific mRNAs may have a key influence on BC progression, impacting key signaling pathways. Furthermore, emphasis is placed on the significance of select lncRNAs as biomarkers for patient stratification and their potential use as prognostic indicators for individuals with BC. Indeed, quantitative assessment of lncRNAs is considered to hold clinical relevance for the diagnosis and prognosis of patients with TNBC. Whilst promising associations have been uncovered, a deeper exploration into these prognostic relationships is necessary.

Supplementary Material

Supporting Data

Supporting Data

Acknowledgements

Not applicable.

Funding

The present study was supported by Instituto Mexicano del Seguro Social through the Health Research Coordination (grant nos. FIS/IMSS/PROT/PRIO/13/027 ‘Temas prioritarios IMSS 2013’) and funded by ‘Protocolos de investigación multidisciplinarios de Cohorte de Largo Aliento sobre Temas Prioritarios en el IMSS 2023’ (grant no. R-2020-785-154).

Availability of data and materials

The data generated in the present study may be found in the GEO under accession number GSE270721. (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE270721). The data generated in the present study may be requested from the corresponding author.

Authors' contributions

RCO conducted RT-qPCR experiments, bioinformatic analysis, analyzed results and wrote the manuscript. CVV obtained RNA from samples and conducted IHC experiments. KVS conducted RT-qPCR experiments, bioinformatics analysis and wrote the manuscript. AMM confirmed the histopathological diagnosis and selected representative areas of tumor for molecular analysis. MERT, JT, NRS and HM analyzed the results. PPS conceived the present study, conducted experiments, analyzed and revised results and wrote the manuscript. RCO and PPS confirm the authenticity of all the raw data. All authors have approved the final version of the manuscript.

Ethics approval and consent to participate

The present study was approved by the Research Ethics Committee and by the National Committee for Scientific Research from Mexican Institute of Social Security (approval nos. R-2013-785-045 and R-2020-785-154).

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

References