Identification of circRNA/lncRNA-miRNA-Transcription Factor-mRNA Networks to Identify Potential Biomarkers in Response to Neoadjuvant Therapy for Breast Cancer

Introduction

Breast cancer (BC) is the most frequently diagnosed malignancy in among women globally.^1,2 Recent studies indicate a significant increase in the incidence of BC.³ BC can be classified into several subtypes depending on the expression levels of estrogen receptor (ER), progesterone receptor (PR), and human epidermal growth factor receptor 2 (HER2) in tumor cells: luminal A (ER+, PR+, HER2-), luminal B (ER+, PR±, HER2 +), HER2+ (ER-, PR-, HER2+), and triple-negative breast cancer (TNBC) (ER-, PR-, HER2-).^4,5

Traditionally, therapeutic options for treating BC include surgery, chemotherapy, radiation therapy, targeted therapy, and endocrine therapy.³ In recent years, neoadjuvant therapy (NAT), followed by definitive surgery, has become the standard option for early-stage or locally advanced BC, as it enhances the likelihood of breast-conserving surgery.^6,7 NAT, also known as preoperative therapy, offers several potential benefits for BC patients: first, it can facilitate tumor downstaging, enabling previously inoperable patients to meet surgical criteria; second, for patients who initially require total mastectomy, NAT enables the possibility of breast-sparing surgery, thereby enhancing surgical outcomes and improving the patient′s quality of life.⁸ However, it is important to note that not all patients may benefit from NAT. Consequently, identifying biomarkers that may predict responsiveness to NAT is crucial for guiding treatment decisions.

Non-coding RNAs (ncRNAs) comprise a category of RNA molecules that typically do not encode proteins but play crucial roles in regulating cellular signal transduction and biological processes.⁹ NcRNAs, including long noncoding RNA (lncRNAs), circular RNA (circRNAs), and microRNAs (miRNA), are pivotal in cancer progression.^10–12 Yuan et al indicated that BC patients with low circRNA CDR1-AS expression exhibited increased sensitivity to paclitaxel-cisplatin-based neoadjuvant chemotherapy.¹³ Meanwhile, Sha et al reported a significant association between the expression level of lncRNA LOC100505851 and a higher pathological complete response rate after NAT as well as improved relapse-free survival and overall survival in BC patients.¹⁴ These findings suggest a potential relationship between circRNA/lncRNAs and the response to NAT treatment in BC patients.

Mechanistically, both lncRNAs and circRNAs can function as competitive endogenous RNAs (ceRNAs) by binding to miRNAs, thus influencing the effects of miRNAs on target genes.¹⁵ Nevertheless, the involvement of circRNA/lncRNA-miRNA axes in the antitumor efficacy of NAT for BC is still unknown. Therefore, this study aims to identify and validate circRNA/lncRNA–miRNA–mRNA ceRNA networks associated with the response to NAT in BC, thereby providing a novel perspective for BC research.

Materials and Methods

Clinical Samples

Eighteen paired tumor tissues (pre-therapy tumor tissues) and adjacent normal tissues were obtained from eighteen patients with BC who underwent surgery at the Affiliated Hospital of Inner Mongolia Medical University. The inclusion criteria were as follows: (1) all participants met the diagnostic criteria for BC, which was confirmed through pathological examination; (2) complete patient clinicopathological information was provided; (3) complete clinical and imaging data were available; and (4) no patient had a history of other malignant tumors. The exclusion criteria were as follows: (1) patients who were clinically diagnosed with BC without biopsy for pathological diagnosis, and (2) patients who had undergone treatment for cancer.

Furthermore, eighteen tumor samples (post-NAT tumor tissues) were collected from three BC patients who received neoadjuvant therapy before tissue resection. The inclusion criteria were as follows: (1) patients diagnosed clinically and pathologically as locally advanced (American Joint Committee on Cancer [AJCC] Stage III) (except T3N1M0); (2) patients who could undergo surgery but did not meet the requirements for preservation of the breast and armpit; (3) patients with clinically positive lymph nodes and tumor size > 2 cm; and (4) no family history of malignant tumors and no organ failure. The exclusion criteria were as follows: (1) inability to tolerate chemotherapy, (2) patients who had uncontrollable and serious infections, and (3) patients who had serious underlying disease. This study was approved by the Ethics Committee of the Affiliated Hospital of Inner Mongolia Medical University (approval number: YJS2024055), in accordance with the Declaration of Helsinki. Written informed consent was obtained from all patients.

NAT Plan

HER2 negative/ER positive BC patient received 8 cycles (21 days/cycle) of chemotherapy consisting of cyclophosphamide, epirubicin, paclitaxel, and docetaxel. During the first four cycles of chemotherapy, patients received chemotherapy with epirubicin (100 mg/m²) and cyclophosphamide (600 mg/m²). During the final four chemotherapy treatment cycles, patients received docetaxel (100 mg/m²) or paclitaxel (80 mg/m²). After 8 cycles of neoadjuvant chemotherapy, the patients underwent surgery.

HER2 positive BC patient received 6 cycles of chemotherapy consisting of paclitaxel (80 mg/m²) or docetaxel (100 mg/m²), carboplatin (carboplatin dose was based on the area under the plasma concentration–time curve of 6), trastuzumab (8 mg/kg for the first time and 6 mg/kg for the subsequent time points), and pertuzumab (840 mg for the first time and 420 mg for the subsequent time points). After receiving 6 cycles of neoadjuvant chemotherapy, the patients underwent surgery.

Patient with TNBC received eight cycles (21 days/cycle) of chemotherapy consisting of docetaxel, paclitaxel, epirubicin, and cyclophosphamide. During the first four cycles of chemotherapy, patients received chemotherapy with epirubicin (100 mg/m²) and cyclophosphamide (600 mg/m²). During the final four chemotherapy treatment cycles, patients received docetaxel (100 mg/m²) or paclitaxel (80 mg/m²). After 8 cycles of neoadjuvant chemotherapy, the patients underwent surgery.

Ribosomal RNA (rRNA)-Depleted RNA Sequencing (RNA-Seq)

Total RNA was extracted from the tumor tissues using the mirVana™ miRNA ISOlation Kit (Ambion-1561, ThermoFisher Scientific). RNA integrity was determined using an Agilent 2100 Bioanalyzer. Next, TruSeq Stranded Total RNA with Ribo-Zero Gold (RS-122-2301, Illumina) was used to construct libraries according to the manufacturer′s instructions. Subsequently, libraries were sequenced on the Illumina sequencing platform.

Trimmomatic software was initially employed to remove adapters, low-quality bases, N-bases, and low-quality reads, resulting in high-quality clean reads.¹⁶ The clean reads were then aligned to the reference genome of the experimental species using the hisat2.¹⁷

For linear RNAs, including lncRNAs and mRNAs, sequencing reads from each sample were aligned with the mRNA transcript sequences, known lncRNA sequences, and predicted lncRNA sequences using bowtie2. Subsequently, eXpress was used to perform quantitative gene analysis, resulting in the acquisition of FPKM values and counts (the number of reads for each gene in each sample). For circRNA prediction, we utilized BWA software (version 0.7.5a) to align the sequencing reads of each sample with the reference genome, which facilitated the generation of SAM files. The default parameters of CIRI software (version v2.0.3) was employed to detect paired chiastic clipping signals, allowing for the prediction of circRNA sequences based on junction reads and GT-AG cleavage signals.¹⁸

To identify differentially expressed mRNAs (DEmRNAs), differentially expressed lncRNAs (DElncRNAs) and differentially expressed circRNAs (DEcircRNAs), the estimateSizeFactors function in the DESeq (version 1.42.1) R package were employed to normalize the counts.¹⁹ Subsequently, p values and fold-change values were calculated using the nbinomTest function. The False Discovery Rate (FDR) method was used to correct for multiple hypothesis testing. Genes with |log2FoldChange|>1 and FDR<0.05, were considered differentially expressed.

The Cancer Genome Atlas (TCGA)-BC Dataset

The mRNA expression profiles, which included 1077 breast cancer (BC) tissues and 106 normal tissues, along with the clinical information of BC patients, were obtained from the TCGA-BC dataset. The Wilcoxon rank-sum test was used to examine the differences in gene expression between the normal and BC groups. R package “survival” (https://CRAN.R-project.org/package=survival) and R package “survminer” (https://CRAN.R-project.org/package=survminer) were used to assess the overall survival in different groups. Differences in survival rates were evaluated using the Log rank test.

Construction of the Co-Expression Network for lncRNA/circRNA/mRNA

Using the assembled sequences of circRNAs and lncRNAs, along with miRNA sequences obtained from the miRBase database, miRanda software was used to predict the degree of complementary matching (total score) and free energy (total energy) of the composite structure formed by their interactions.²⁰

Functional Analysis

The “clusterProfiler” (version 4.8.3) package in R language was employed for Gene ontology [GO, including biological process (BP), molecular function (MF), cellular component (CC)] and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses.²¹ The threshold for significance was set at p-adjust < 0.05.

Real-Time Quantitative Reverse-Transcription PCR (RT-qPCR)

Total RNA was extracted from tissue samples using the Redzol reagent (FTR-50, SBS Genetech Co., Ltd). Complementary DNA (cDNA) was synthesized from the extracted RNA using the Reverse Transcription kit (QP056, GeneCopeia). The cDNA was subsequently amplified by quantitative PCR using the 2× SYBR Green qPCR Master Mix (None ROX) (MPC2203026, Servicebio) on an iQ5Real-TimePCR system. The relative gene expression was calculated using the 2^–ΔΔCt method by normalization to U6 or GAPDH. The primers used in this study were listed in Table 1.

Table 1 Primers Used in This Study

Statistical Analysis

The Pearson correlation method was used for correlation analysis. Statistical significance was set at p < 0.05. All statistical analyses were performed using the R software version 4.3.3.

For the analysis of the experimental data, the Student′s t-test was employed to assess differences between the groups. Data were presented as mean ± SD, with p values < 0.05 considered statistically significant.

Results

Identification of DElncRNAs Between Groups

R language was used to identify DElncRNAs between pre-therapy tumor and adjacent normal groups, and between post-NAT tumor and pre-therapy tumor groups. Compared to the adjacent normal group, 134 lncRNAs were downregulated and 128 lncRNAs were upregulated in the pre-therapy tumor group (Figure 1A, B and Table S1). Compared with the pre-therapy tumor group, 82 lncRNAs were downregulated and 105 lncRNAs were upregulated in the post-NAT tumor group (Figure 1C, D and Table S1). Next, a Venn diagram was used to distinguish the common lncRNAs between the two comparison groups (pre-therapy tumor vs adjacent normal and post-NAT tumor vs pre-therapy tumor). A total of 25 DElncRNAs were found to be dysregulated across these two comparison groups (Figure 1E and F). As shown in Figure 1E, compared to the adjacent normal group, the levels of 15 lncRNAs were remarkably elevated in the pre-therapy tumor group; however, these levels were notably reduced in the post-NAT tumor group compared to the pre-therapy tumor group. As shown in Figure 1F, the levels of 10 lncRNAs in the pre-therapy tumor group were considerably lower than those in the adjacent normal group, whereas these levels were significantly higher in the post-NAT tumor group than in the pre-therapy tumor group.

Figure 1 Identification of DElncRNAs between groups. (A and B) Volcano and heatmap of DElncRNAs between pre-therapy tumor (n=3) and tumor adjacent normal (n=3) samples. Blue color, downregulated lncRNAs; red color, upregulated lncRNAs. (C and D) Volcano and heatmap of DElncRNAs between post-NAT (n=3) and pre-therapy tumor (n=3) samples. Blue color, downregulated lncRNAs; red color, upregulated lncRNAs. (E) Venn diagram of 15 common DElncRNAs between the upregulated genes in the pre-therapy tumor vs adjacent normal comparison group and the downregulated genes in the post-NAT vs pre-therapy group. (F) Venn diagram of 10 common DElncRNAs between the downregulated genes in the pre-therapy tumor vs adjacent normal comparison group and the upregulated genes in the post-NAT vs pre-therapy group.

Of these 25 DElncRNAs, 20 were detected across all samples in the three groups (Table S1). Next, the miRanda database was used to screen for miRNAs targeted by these 20 DElncRNAs. As shown in Figure S1, six lncRNAs (ADAMTS9-AS2, MIAT, LINC0067, LINC02884, LOC105372310, and LOC10099662) and 149 potential miRNAs were identified using the miRanda database.

Identification of DEcircRNAs Between Groups

Next, we screened DEcircRNAs between adjacent normal and pre-therapy tumor groups, and between pre-therapy tumor and post-NAT tumor groups. Compared to the adjacent normal group, 1384 circRNAs were notably increased and 3782 circRNAs were notably decreased in the pre-therapy tumor group (Figure 2A, B and Table S2). Compared to the pre-therapy tumor group, 1582 circRNAs were significantly reduced, and 3732 circRNAs were remarkably elevated in the post-NAT tumor group (Figure 2C, D and Table S2). Next, a Venn diagram was used to distinguish the common circRNAs between the two comparison groups (pre-therapy tumor vs adjacent normal and post-NAT tumor vs pre-therapy tumor). A total of 2693 DEcircRNAs were dysregulated across these two comparison groups (Figure 2E and F). As shown in Figure 2E, compared to the adjacent normal group, the levels of 1047 circRNAs increased in the pre-therapy tumor group; however, compared to the pre-therapy tumor group, these levels decreased in the post-NAT tumor group. As shown in Figure 2F, the levels of 1646 circRNAs in the pre-therapy tumor group were notably lower than those in the adjacent normal group, whereas these levels were significantly higher in the post-NAT tumor group than in the pre-therapy tumor group.

Figure 2 Identification of DEcircRNAs between groups. (A and B) Volcano and heatmap of DEcircRNAs between pre-therapy tumor (n=3) and tumor adjacent normal (n=3) samples. Blue color, downregulated circRNAs; red color, upregulated circRNAs. (C and D) Volcano and heatmap of DEcircRNAs between post-NAT (n=3) and pre-therapy tumor (n=3) samples. Blue color, downregulated circRNAs; red color, upregulated circRNAs. (E) Venn diagram of 1047 common DEcircRNAs between the upregulated genes in the pre-therapy vs adjacent normal comparison group and the downregulated genes in the post-NAT and pre-therapy group. (F) Venn diagram of 1646 common DEcircRNAs between the downregulated genes in the pre-therapy vs adjacent normal comparison group and the upregulated genes in the post-NAT and pre-therapy group.

In these 2693 DEcircRNAs, a total of 104 circRNAs were detected across all samples in the three groups (Table S2), and the target miRNAs of these DEcircRNAs were predicted using the miRanda database. Based on sequence complementarity (total score>200) and free energy of formation (total energy<80), 77 miRNAs that putatively targeted by eight DEcircRNAs (circRNA_31003, circRNA_34064, circRNA_42934, circRNA_40273, circRNA_15026, circRNA_42276, circRNA_50033, and circRNA_46221) were identified using the miRanda database, and the miRNA-circRNA pairs were constructed (Figure S2 and Table S3).

Identification of Common miRNAs Targeted by 6 DElncRNA and 8 circRNAs

Next, the common miRNAs targeted by the six DElncRNAs and eight circRNAs were screened using Venn diagrams. As shown in Figure 3A, 39 miRNAs that potentially interact with both circRNAs and lncRNAs were identified. Subsequently, a circRNA/lncRNA-miRNA network was constructed, containing eight circRNAs, six lncRNAs, and 39 miRNAs (Figure 3B).

Figure 3 Identification of common miRNAs targeted by 6 DElncRNA and 8 circRNAs. (A) Venn diagram of 39 common miRNAs between DElncRNA-target miRNAs and DEcircRNA-target miRNAs. (B) The circRNA/lncRNA-miRNA network. Blue color represent lncRNAs/circRNAs that were downregulated in the pre-therapy tumor vs adjacent normal comparison but upregulated in the post-NAT vs pre-therapy tumor comparison; red color represents lncRNAs/circRNAs that were upregulated in the pre-therapy tumor vs adjacent normal comparison but downregulated in the post-NAT vs pre-therapy tumor comparison; green color, miRNAs.

Identification of DEmRNAs Between Groups

DEmRNAs were then identified between the adjacent normal and pre-therapy tumor groups, as well as between the pre-therapy tumor and post-NAT tumor groups. In comparison to the adjacent normal group, 1373 genes exhibited a significant increase, while 1014 genes demonstrated a notable decrease in expression within the pre-therapy tumor group (Figure 4A, B and Table S4). When comparing the post-NAT tumor group to the pre-therapy tumor group, 794 genes were significantly reduced, and 454 genes were remarkably elevated in the post-NAT tumor group (Figure 4C, D and Table S4). Furthermore, a Venn diagram was used to distinguish the common DEmRNAs between the two comparison groups (pre-therapy tumor vs adjacent normal and post-NAT tumor vs pre-therapy tumor). A total of 936 DEmRNAs were found to be dysregulated across these two comparison groups (Figure 4E and F).

Figure 4 Identification of DEmRNAs between groups. (A and B) Volcano and heatmap of DEmRNAs between pre-therapy tumor (n=3) and adjacent normal (n=3) samples. Blue color, downregulated mRNAs; red color, upregulated mRNAs. (C and D) Volcano and heatmap of DEmRNAs between post-NAT (n=3) and pre-therapy tumor (n=3) samples. Blue color, downregulated mRNAs; red color, upregulated mRNAs. (E) Venn diagram of 619 common DEmRNAs between the upregulated genes in the pre-therapy vs adjacent normal comparison group and the downregulated genes in the post-NAT vs pre-therapy comparison group. (F) Venn diagram of 317 common DEmRNAs between the downregulated genes in the pre-therapy vs adjacent normal comparison group and the upregulated genes in the post-NAT vs pre-therapy comparison group. (G and H) GO and KEGG analyses were performed on 619 common DEmRNAs screened in (E). (I and J) GO and KEGG analyses were performed on 317 common DEmRNAs screened in (F).

As illustrated in Figure 4E, compared to the adjacent normal group, the levels of 619 genes were found to be elevated in the pre-therapy tumor group; however, when compared to the pre-therapy tumor group, these levels decreased in the post-NAT tumor group. The results of GO and KEGG analyses showed that these 619 DEmRNAs were significantly enriched in 754 GO-BP terms (eg “Chromosome segregation”), 69 GO-CC terms (eg “External side of plasma membrane”) and 49 GO-MF terms (eg “Extracellular matrix structural constituent”), and 52 KEGG pathways (eg “Cell cycle”) (Figure 4G, H and Table S5).

As shown in Figure 4F, the levels of 317 genes in the pre-therapy tumor group were notably lower than those in the adjacent normal group. In contrast, these levels were significantly elevated in the post-NAT tumor group when compared to the pre-therapy tumor group (Figure 4F). The results of GO and KEGG analyses indicated that these 317 DEmRNAs were significantly enriched in 20 GO-BP terms (eg “Regulation of membrane potential”), 35 GO-CC terms (eg “Basolateral plasma membrane”) and 24 GO-MF terms (eg “Glycosaminoglycan binding”), and 5 KEGG pathways (eg “PPAR signaling pathway”, “AMPK signaling pathway”) (Figure 4I, J and Table S6).

Screening Hub Differential Expressed-Transcription Factors (DE-TFs)

To investigate the potential TFs involved in the regulation of these DEmRNAs, Venn diagrams were employed to analyze the overlap between 1562 TFs and 936 DEmRNAs. This analysis led to the identification of 58 DE-TFs (Figure 5A and Table S7). To further refine the selection of candidate DE-TFs, the top five DE-TFs (HOXC11, HOXD11, PITX1, NKX2-2, and PRAME), which exhibited the most significant log2 fold change (log2 FC) in comparisons between pre-therapy tumor and adjacent normal tissue samples, were selected for subsequent analysis (Table S7).

Figure 5 Screening hub differential expressed-transcription factors (DE-TFs). (A) Venn diagram of 58 common DE-TFs between the 1562 TFs and 936 DEmRNAs. (B) Heatmap of expression levels of HOXC11 and its five target DEmRNAs among adjacent normal, pre-therapy and post-NAT groups. The color gradient from blue to red represents that the gene expression level ranges from low to high. (C) Heatmap of expression levels of NKX2-2 and its five target DEmRNAs among adjacent normal, pre-therapy and post-NAT groups. The color gradient from blue to red represents that the gene expression level ranges from low to high. (D) Heatmap of expression levels of PRAME and its five target DEmRNAs among adjacent normal, pre-therapy and post-NAT groups. The color gradient from blue to red represents that the gene expression level ranges from low to high.

Next, the ChIp-seq data for these five TFs were downloaded from the Cistrome Data Browser (http://cistrome.org/db/) to predict their target genes. Of the 936 DEmRNAs analyzed, 58 represented DE-TFs and were excluded from downstream targeting analysis. The remaining 878 DEmRNAs were thus considered potential targets. Utilizing ChIP-seq data from the Cistrome database, we identified three DE-TFs (HOXC11, NKX2-2, and PRAME) with potential downstream targets (Table S7).

To predict the target DEmRNAs of HOXC11, we interrogated its ChIP-seq data (GEO: GSM1306052) with a screening threshold greater than 0.1.²² As presented in Table 2, the binding sites and binding strengths between five target DEmRNAs (CCL5, HBA2, HBA1, HOXC12, and ELAVL4) and HOXC11 were predicted. Specifically, HOXC11 was found to have one binding region with HBA2, HBA1, CCL5, and HOXC12, whereas HOXC11 and ELAVL4 had three binding regions (Table 2). The expression levels of HOXC11 and its five target genes across the three groups were illustrated in a heatmap (Figure 5B).

Table 2 The Binding Regions and Binding Score Between HOXC11 and Its Target Genes

The ChIp-seq data for NKX2-2 (GEO: GSM1208780) were analyzed to identify its target DEmRNAs using a screening threshold set to greater than 1.²³ As presented in Table 3, the binding sites and binding strengths between five target DEmRNAs (NEK2, KIF20A, SPC24, LRRC49, and CDC45) and NKX2-2 were predicted. Specifically, NKX2-2 interacts with KIF20A and CDC45 through one binding region each, and with LRRC49 through three binding regions. For NEK2, NKX2-2 binds to each canonical transcript through one binding region and also interacts with one protein-coding sequence. Additionally, NKX2-2 binds to six canonical transcripts of SPC24. The expression levels of NKX2-2 and its five target genes across the three groups were displayed in a heatmap (Figure 5C).

Table 3 The Binding Regions and Binding Score Between NKX2-2 and Its Target Genes

Chip-seq data for PRAME (GEO: GSM648586) were used to predict the targeted DEmRNAs of PRAME, with a screening threshold set to greater than 1.²⁴ As shown in Table S7, 281 target DEmRNAs of PRAME were predicted, encompassing 1373 predicted binding regions. Thus, to further refine hub target genes of PRAME, these 281 candidates were subjected to survival and differential expression analyses based on data from the TCGA-BC cohort. Among them, the levels of RAD54L, PGM2L1, and RPGRIP1L exhibited significantly higher expression in tumor tissues compared to controls, and their high expression was associated with poor prognosis in BC patients (Figure S3A–S3F). Conversely, FBXO2 and FAM181B showed significantly lower expression in tumor tissues, and their reduced expression was also associated with poor prognosis in patients with BC (Figure S3G–S3J). The predicted binding sites and strengths between these five genes and PRAME are detailed in Table 4. PRAME interacts with FBXO2 and FAM181B through one binding region each. For RAD54L, PRAME binds to each of its three canonical transcripts via one binding region and also interacts with two protein-coding sequences. In the case of PGM2L1, PRAME binds to one canonical transcript and one protein-coding sequence. For RPGRIP1L, PRAME interacts with six canonical transcripts and nine protein-coding sequences. The expression levels of PRAME and its five target genes across the three groups were illustrated in a heatmap (Figure 5D).

Table 4 The Binding Regions and Binding Score Between PRAME and Its Target Genes

Construction of circRNA/lncRNA–miRNA–TF-mRNA Network

Furthermore, the miRWalk database was conducted to predict the interactions between 39 miRNAs and 3 hub TFs (HOXC11, NKX2-2, and PRAME), and miRNAs with a binding prediction score exceeding 0.9 were designated as hub miRNAs (Table S8).

Finally, based on the identified circRNA/lncRNA-miRNA pairs (Figure 3B) and miRNA-TF-mRNA pairs (Tables 2–4, Table S8), circRNA/lncRNA-miRNA-TF-mRNA ceRNA networks were constructed (Figure S4A–S4C).

According to the ceRNA mechanism, circRNAs and lncRNAs are positively correlated with target gene expression.^25,26 Given that the expression levels of three TFs (HOXC11, NKX2-2, and PRAME) were potentially elevated in the pre-therapy tumor group compared to the adjacent normal group, we mainly focused on the upregulated circRNAs/lncRNAs in pre-therapy tumor samples (Figure 6A–C).

Figure 6 Construction of upregulated circRNA/upregulated lncRNA–miRNA–TF-mRNA network. (A) The upregulated circRNA/upregulated lncRNA–miRNA–HOXC11-mRNA network. (B) The upregulated circRNA/upregulated lncRNA–miRNA–NKX2-2-mRNA network. (C) The upregulated circRNA/upregulated lncRNA–miRNA–PRAME-mRNA network. Blue color represent mRNAs that were downregulated in the pre-therapy tumor vs adjacent normal comparison but upregulated in the post-NAT vs pre-therapy tumor comparison; pink color represents lncRNAs/circRNAs/TFs/mRNAs that were upregulated in the pre-therapy tumor vs adjacent normal comparison but downregulated in the post-NAT vs pre-therapy tumor comparison; green color, miRNAs.

As shown in Figure 6A–C, circRNA_31003 and MIAT may indirectly regulate all three TFs, suggesting that circRNA_31003 or MIAT may be the hub circRNA or lncRNA in BC. Additionally, circRNA_40273, circRNA_42276, and LOC105372310 may indirectly regulate two TFs (HOXC11 and NKX2-2), suggesting that these genes may play a significant role in BC.

Validation of Expression Levels of Key circRNAs/lncRNAs, miRNAs, TFs and mRNAs Among Adjacent Normal, Pre-Therapy Tumor and Post-NAT Tumor Groups

RT-qPCR was performed to validate the expression levels of key circRNAs, lncRNAs, miRNAs, TFs, and mRNAs across the adjacent normal, pre-therapy tumor, and post-NAT tumor groups in the circRNA/lncRNA-miRNA-TF-mRNA ceRNA networks (Figure 7A-L). The results indicated that the expression levels of circRNA_31003, circRNA_42276, MIAT, HOXC11, NKX2-2, PRAME, CCL5, NEK2, and RAD54L were significantly upregulated in pre-therapy tumor tissues compared with adjacent normal tissues; however, these levels were decreased in the post-NAT tumor group compared to pre-therapy tumor group (Figure 7A–C, G–L). In contrast, the expression levels of miR-1225-3p, miR-661, and miR-143-5p were markedly downregulated in pre-therapy tumors relative to adjacent normal tissues; however, these levels were significantly elevated in the post-NAT tumor group compared to pre-therapy tumor group (Figure 7D–F).

Figure 7 Validation of expression levels of key circRNAs/lncRNAs, miRNAs, TFs and mRNAs among adjacent normal, pre-therapy tumor and post-NAT tumor groups. (A–L) RT-qPCR was employed to assess the expression levels of key molecules, including two circRNAs (circRNA_31003, circRNA_42276), one lncRNA (MIAT), three miRNAs (miR-1225-3p, miR-661, miR-143-5p), three TFs (HOXC11, NKX2-2, PRAME), and three mRNAs (CCL5, NEK2, RAD54L) across adjacent normal, pre-treatment tumor, and post-NAT tumor tissues. ***P<0.001, n=15.

Discussion

In this study, we conducted rRNA-depleted RNA-seq to identify DElncRNAs, DEcircRNAs, and DEmRNAs across adjacent normal tissues, pre-therapy tumor tissues and post-NAT tumor tissues with the aim of constructing lncRNA-miRNA-mRNA networks that respond to NAT in BC. A total of six DElncRNAs (ADAMTS9-AS2, MIAT, LINC0067, LINC02884, LOC105372310, and LOC10099662) and eight DEcircRNAs (circRNA_31003, circRNA_34064, circRNA_42934, circRNA_40273, circRNA_15026, circRNA_42276, circRNA_50033, and circRNA_46221) were identified between the adjacent normal and pre-therapy tumor groups and between the pre-therapy and post-NAT tumor groups. We also explored their downstream miRNA targets to construct the DE circRNA/DE lncRNA-miRNA axis. These lncRNAs and circRNAs may serve as candidate biomarkers for the response of BC patients to NAT treatment. Previous studies have highlighted the roles of ADAMTS9-AS2 and MIAT in BC.^27,28 For instance, Ni et al found that DAMTS9-AS2 levels were decreased in TNBC tissues, and that DAMTS9-AS2 overexpression suppressed TNBC progression.²⁷ Conversely, MIAT is overexpressed in BC and its overexpression has been shown to promote BC progression.^28,29 These findings suggest that ADAMTS9-AS2 may function as a tumor suppressor and MIAT may act as an oncogene in BC. Our results showed that ADAMTS9-AS2 levels were reduced, whereas MIAT levels were elevated in pre-therapy tumor tissues compared with those in adjacent normal tissues, which is consistent with previous research. Conversely, in post-NAT tumor tissues, ADAMTS9-AS2 levels were found to be increased, whereas MIAT levels were decreased when compared to the pre-therapy tumor group. These results implied that ADAMTS9-AS2 and MIAT may serve as potential biomarkers for identifying BC patients who respond to NAT treatment. However, the roles of other 4 lncRNAs and eight circRNAs in BC progression have not yet been investigated, warranting further studies.

TFs regulate the expression of downstream genes by binding to specific DNA sequences, thereby promoting or inhibiting transcription.³⁰ They play dual roles in tumorigenesis, functioning as both tumor suppressor genes and oncogenes, which underscores their significant importance in cancer development.^31,32 In this study, we identified 58 common TFs between the adjacent normal and pre-therapy tumor groups, as well as between the pre-therapy tumor and post-NAT tumor groups. To focus on these genes, we selected the top five based on their expression ratios, identifying three DE-TFs (HOXC11, NKX2-2, and PRAME) that have potential downstream targets. These three DE-TFs (HOXC11, NKX2-2, and PRAME) were all found to be elevated in pre-therapy tumor tissues compared with adjacent normal tissues. In contrast, the levels of HOXC11, NKX2-2, and PRAME were reduced in post-NAT tumor tissues when compared to the pre-therapy tumor group. Consequently, we further examined their upstream miRNA targets and downstream DEmRNA targets to construct the DEcircRNA/DElncRNA-miRNA-DE-TF-DEmRNA axis.

HOXC11 has been identified as an oncogene in various cancers, such as colon adenocarcinoma and pancreatic adenocarcinoma, with elevated levels of HOXC11 correlating with poor overall survival in cancer patients.³³ Furthermore, nuclear HOXC11 has been associated with poor disease-free survival in patients with BC,³⁴ suggesting its relevance to BC prognosis. Furthermore, evidence has shown that HOXC11 may serve as a transcriptional activator and repressor,^35,36 acting as a dual transcriptional regulator in tumorigenesis. Using the Cistrome Data Browser database, we predicted that CCL5, HOXC12, ELAVL4, HBA1, and HBA2 are potential downstream targets of HOXC11. Meanwhile, the levels of CCL5, HOXC12, and ELAVL4 were elevated, and HBA1 and HBA2 levels were reduced in pre-therapy tumor tissues compared with those in adjacent normal tissues. Evidence has shown that CCL5 and ELAVL4 may play oncogenic roles, whereas HBA1 may serve as a potential tumor suppressor gene in different cancers.^37,38 Additionally, BC patients exhibiting high CCL5 expression had poorer disease-free survival.³⁹ Thus, we speculated that HOXC11 might promote BC progression by potentially promoting the transcription of oncogenes (CCL5 and ELAVL4) and suppressing the transcription of the tumor-suppressor gene (HBA1). However, this assumption requires further investigation in future studies.

NKX2.2 plays a dual role in cancer progression.^40,41 Several studies have indicated that NKX2.2 is reduced in several malignancies, such as colorectal cancer and osteosarcoma, where it functions as a tumor suppressor.^40,42 Conversely, a study by Lawson et al revealed that overexpression of NKX2.2 can accelerate cell survival in small cell lung cancer, implying a potential oncogenic role for NKX2.2.⁴¹ Using the Cistrome Data Browser database, we predicted that SPC24, KIF20A, NEK2, and CDC45 are potential downstream targets of NKX2.2. Evidence has shown that these four target genes are upregulated in various tumor tissues, and their elevated levels correlate with poor prognosis.^43–46 Notably, SPC24, KIF20A, and NEK2 were upregulated in BC tissues when compared to normal controls,^47–49 aligning with our findings. NKX2.2 possesses both transcriptional activation and transcriptional repression capabilities.⁵⁰ Thus, we speculate that NKX2.2 may promote BC progression by activating the expression of its downstream genes; nevertheless, these hypotheses warrant further investigation in future studies.

Al-Khadairi et al discovered that PRAME can facilitate TNBC progression by enhancing epithelial-to-mesenchymal transition.⁵¹ Additionally, high PRAME expression has been correlated with unfavorable outcomes in patients with BC.^52,53 A previous study demonstrated that PRAME suppresses the expression of its downstream targets.⁵⁴ In contrast, PRAME activates meiotic genes in uveal melanoma.⁵⁵ These findings indicate that PRAME may play a dual role in regulating the expression of downstream genes. Our results showed that RPGRIP1L, RAD54L, and PGM2L1 are potential targets of PRAME. These three genes have been identified as oncogenes in cancers.^56–58 Consistent with prior studies, our results showed that PRAME, RPGRIP1L, RAD54L, and PGM2L1 levels were elevated in pre-therapy tumor tissues compared to those in adjacent normal tissues. Thus, we speculate that PRAME may promote BC progression through upregulation of these three genes; however, these assumptions warrant further investigation in future studies.

LncRNAs or circRNAs can function as ceRNAs by competitively occupying the shared binding sequences of miRNAs. This interaction diminishes the capacity of miRNAs to interact with target genes, thereby altering their expression.^59–61 Luan et al demonstrated that lncRNA MIAT acts as a ceRNA by sponging miR-155-5p, thereby enhancing the expression of DUSP7 in breast cancer.²⁹ Zhang et al reported that MIAT competitively binds to miR-411-5p, leading to the restoration of STAT3 (a TF) expression, which subsequently upregulates the levels of its downstream target, PD-L1, in hepatocellular carcinoma cells.⁶² Additionally, Cao et al indicated that circRNA circRNF20 plays a role in BC progression via the miR-487a/HIF-1α/HK2 axis.⁶⁰ These findings suggest that lncRNAs and circRNAs may influence tumor progression by sponging miRNAs, thereby modulating the TF-mRNA axis. In this study, our results indicated that circRNA_31003 and MIAT may indirectly regulate the expression of these three TFs (HOXC11, NKX2-2, and PRAME) via sponging miRNAs, and circRNA_40273, circRNA_42276, and LOC105372310 may indirectly regulate two TFs (HOXC11 and NKX2-2) via similar mechanisms, suggesting that these three circRNAs and two lncRNAs may play significant roles in BC. Furthermore, all these three TFs, three circRNAs and two lncRNAs were all found to be elevated in pre-therapy tumor tissues compared to adjacent normal tissues. In contrast, these levels were reduced in the post-NAT tumor group compared to the pre-therapy tumor group, suggesting that these genes may be potential factors in the response to NAT in BC.

In general, as illustrated in Figure 6A–C, circRNA_31003, circRNA_40273, circRNA_42276, LOC105372310, and MIAT may affect HOXC11 expression by sponging miRNAs (such as miR-1285-5p, miR-1225-3p, or miR-4739), which in turn affects the levels of CCL5, ELAVL4, and HBA1. These circRNAs and lncRNAs may also modulate HKX2-2 expression by sponging miRNAs (such as miR-939-5p, miR-6797-3p, or miR-661), thereby affecting the levels of SPC24, KIF20A, NEK2, and CDC45. Additionally, circRNA_31003 and MIAT may regulate PRAME expression by sponging miR-143-5p, which subsequently influences the RPGRIP1L, RAD54L, and PGM2L1 levels. These circRNA/lncRNA-miRNA-TF mRNA regulatory networks may play a potential role in the response to NAT in BC.

In this study, we employed RT-qPCR to assess the expression levels of key molecules, including two circRNAs (circRNA_31003, circRNA_42276), one lncRNA (MIAT), three miRNAs (miR-1225-3p, miR-661, miR-143-5p), three TFs (HOXC11, NKX2-2, PRAME), and three mRNAs (CCL5, NEK2, RAD54L), within proposed circRNA/lncRNA-miRNA-TF-mRNA regulatory networks across adjacent normal, pre-therapy tumor, and post-NAT tumor tissues. The RT-qPCR results revealed that, compared with adjacent normal tissues, the expression levels of circRNA_31003, circRNA_42276, MIAT, HOXC11, NKX2-2, PRAME, CCL5, NEK2, and RAD54L were significantly upregulated in pre-therapy tumor samples. However, these levels were downregulated in the post-NAT tumor group compared to the pre-therapy tumor group. Conversely, the levels of miR-1225-3p, miR-661, and miR-143-5p were markedly downregulated in pre-therapy tumors relative to normal tissues, but showed a significant increase after NAT. Based on the data presented in Figure 6, we constructed three regulatory axes: circRNA_42276/MIAT-miR-1225-3p-HOXC11-CCL5 axis, circRNA_31003-miR-661-NKX2-2-NEK2 axis, circRNA_31003/MIAT-miR-143-5p-PRAME-RAD54L axis. These axes may play potential roles in BC pathogenesis. Previous studies have identified CCL5, NEK2, and RAD54L as oncogenic factors in various cancers.^39,57,63,64 Inhibition of the CCL5-mediated Akt/NF-κB pathway has been shown to induce apoptosis in ovarian cancer cells.⁶⁵ Similarly, downregulation of NEK2 has been found to inhibit BC cell proliferation and induce apoptosis through the modulation of ERK/MAPK signaling.⁶⁴ Additionally, Wang et al reported that Rad54L facilitates the progression of bladder cancer by influencing cell cycle and senescence.⁵⁷ Our findings suggest that circRNA_31003, circRNA_42276, and MIAT may contribute to BC development via their roles within a ceRNA network. By acting to modulate the expression of CCL5, NEK2, and RAD54L, these circRNAs/lncRNAs may subsequently influence key cellular processes including apoptosis, cell cycle progression, and senescence. Thus, it is plausible that these three axes may contribute to tumorigenesis in BC. Their disruption following NAT may consequently impede key molecular mechanisms essential for tumor development. Targeting these pathways may represent a promising therapeutic strategy to enhance the efficacy of NAT in BC.

A limitation of this study is that the ceRNA network construction was based on bioinformatic analyses utilizing only three samples per group, which may limit statistical power and introduce bias. While the expression of potential key biomarkers has been validated by RT-qPCR in an independent cohort of 15 samples per group, it remains necessary to validate these preliminary observations in future studies involving larger and more diverse cohorts. In addition, the predicted interactions within the ceRNA networks, such as those between circRNAs/lncRNAs and miRNAs, as well as miRNAs and TFs, require experimental confirmation in future studies. Furthermore, although our current findings have identified potential biomarkers that may respond to NAT treatment in BC, it remains to be validated in future large-scale clinical studies with long-term follow-up whether these molecules can serve as independent predictors of patient outcomes. Additionally, functional investigations using cellular and animal models are imperative to definitively confirm the mechanistic roles of these molecules in tumorigenesis and treatment responsiveness in the future.

Conclusion

Our study identified several potential DElncRNAs, DEcircRNAs, DE-TFs, and DEmRNAs that may respond to NAT. Additionally, we constructed associated regulatory networks to explore potential biomarkers that may indicate the response of BC to NAT treatment. These results may have important clinical implications for guiding personalized treatment strategies in the neoadjuvant setting.