The GEO database was exhaustively searched for breast-cancer-related microarray datasets. Four GEO datasets, with accession numbers GSE182471, GSE101410, GSE111504, and GSE159980, were selected for the screen of circRNAs dysregulated in tamoxifen-resistant and metastatic luminal breast cancer. Dataset GSE182471 summarizes circRNA expression profiles from five pairs of breast cancer tissue samples and adjacent non-tumor tissue samples by microarray and RT-qPCR. GSE101410 comprises RNA-seq Luminal subtype-specific circRNAs in breast cancer cells. GSE111504 identified the circRNA expression signatures for breast cancer metastasis. GSE159980 is composed of microarray analysis of tamoxifen-sensitive and -resistant cells. Raw expression data were scored and extracted/normalized using the GEO-online tool GEO2R after applying log2 transformation. The overlapped circRNAs among these four datasets were available by VENNY diagram.

GSE62022, GSE110204, and GSE121172 were analyzed using GEO2R for the luminal-specific miRNAs. The dysregulated miRNAs were subjected to a volcano-diagram presentation. The potential miRNAs for the obtained circRNAs were collected using Circular RNA Interactome (https://circinteractome.nia.nih.gov/index.html). The dysregulated miRNAs from GSE datasets were integrated with the circRNA-targeting miRNAs. GSE48390 was normalized and log2 transformed using GEO2R for differentially expressed genes between the high- and low-risk groups. Genes with absolute logFC values more than 2 and adjusted p-value less than 0.05 were screened as significantly differentially expressed.

Significantly differentially expressed genes were subjected to STRING for the prediction of protein interactions. Networks were constructed using Cytoscape software and the hub genes in the network were identified using cyto-Hubba, a plug-in Cytoscape.

Gene symbols were transformed for the analysis of gene-centric functional and canonical pathway enrichment analysis using OmicShare Tools.

Combining RNA interactions and expression levels, the circRNA-miRNA-mRNA networks were established as two distinct sets using Cytoscape software. The circRNA-miRNA-mRNA axes containing hub genes were used for the visualization of final biological networks.

To be eligible for this study, women had to be at their first diagnosis of luminal-subtype breast cancer (Engstrøm et al., 2013). The candidates received no cancer-related therapy before surgery and immediately underwent breast surgery, followed by tamoxifen treatments. The tumor core biopsies, along with the adjacent normal tissues, were obtained from all patients and preserved in RNAlater for circRNA, miRNA, and mRNA expression evaluation. All women included in this study provided written informed consent. The study was approved by the Research Ethics Committee at Shanghai Tenth People’s Hospital, Tongji University School of Medicine. Information on survival time and distant metastasis was obtained from electronic or paper medical records in our institution.

The cell lines used in the analysis include MCF7, ZR-75-1, MCF7 Tam1, and ZR-75-1 Tam1, purchased from American Type Culture Collection (Rockville, MD, United States). Whereas ZR-75-1 and MCF7 were respectively cultured in 10% FBS-supplemented RPMI-1640 and EMEM, culture media for MCF7 was additionally contained and 10 μg/mL human insulin (Sigma. St. Louis, MO, United States). The tamoxifen-resistant cell lines, MCF-7 Tam1 and ZR-75-1 Tam1, used the same culture media as the parental cell lines, additionally containing 1 µM 4-hydroxytamoxifen (Sigma. St. Louis, MO, United States).

Hsa_circ_0086735 siRNA (si-circ), STAT1 siRNA (si-STAT1), and non-silencing control plasmids (si-NC and ctr-STAT1) were purchased from Guangzhou RiboBio. Mi-1296-5p inhibitor (in-miR) and non-silencing control (in-NC) were purchased from Guangzhou Genecopoeia. The above vectors were transfected into cells by transient transfection using Lipofectamine 3,000 transfection reagents from Invitrogen.

First, total RNA was extracted from tissue homogenate and cells using TRIzol (Life Technologies, Carlsbad, CA, United States). Then, total RNA was purified using the miRNeasy Mini kit (QIAGEN, Hilden, Germany) and assessed for purity using a NanoVue spectrophotometer (Ge BioSciences, United States). Total RNA, combined with anchored oligo (dT)18 primer and reverse transcriptase from the Transcriptor High FidelitycDNA Synthesis Kit (Roche, Germany), was subjected to the first strand synthesis for cDNA. The cDNA samples were examined by qRT-PCR using FastStart™ universal SYBR^® Green (Roche) for determination of miR-1296-5p and STAT1 mRNA, using U6 and beta-actin as the endogenous control. For hsa_circ_0086735, total RNA was prepared to the linear RNA using RNase R digestion (Epicentre Technologies, Madison, WI, United States). RNase R-treated RNA was reversely transcribed to cDNA using SuperScript III First Strand Synthesis System (Invitrogen). qRT-PCR was performed using Power Up Sybr Green Master Mix (Thermo Fisher Scientific). Each sample was assayed in triplicate determinations at least.

Luciferase reporter plasmids, including WT-circ, MUT-circ, WT-STAT1, and MUT-STAT1, were constructed by LandM Biotechnology (Guangzhou, China). These reporter plasmids were singly transfected into MCF7 and ZR-75-1 cells with miR-1296-5p mimic or inhibitor. After incubation, each cells lysate was prepared using two kinds of buffers, Firefly luciferase substrate and Firefly luciferase inhibitor/Renilla luciferase substrate, corresponding to the first measure of the Firefly luciferase activity and second measure of the Renilla luciferase activity. Finally, the normalized luciferase activity was calculated by Firefly luciferase activity/Renilla luciferase activity. Data were collected from at least three independent experiments for each subgroup.

The combination of miR-1296-5p and hsa_circ_0086735, as well as miR-1296-5p and STAT1, was determined by RNA pull-down assay. The miR-1296-5p probe (biotin-miR) and its negative control (biotin-NC) were constructed by Ribobio. Harvested MCF7 and ZR-75-1 were lysed using lysis buffer and then incubated with a biotin-labeled miT-1296-5p probe and streptavidin agarose magnetic beads. Hsa_circ_0086735 and STAT1 mRNA in miR-1296-5p probe pull-down complex were detected by qRT-PCR. Data were collected from at least three independent experiments for each subgroup.

An appropriate number of cells (1 × 10⁵ of MCF7 and 8 × 10⁴ of ZR-75-1 or ZR-75-1 Tam1) were seeded into a 96-well plate. Different periods after plating, 0, 12, 24, 48, and 72 h for MCF7 and ZR-75-1 and 0, 2, 4, 6, 9 days for ZR-75-1 Tam1, 10 μL of CCK-8 solution (Millipore-Sigma) was added to each well, followed by a 2-h incubation. The optical density of each well was measured at 450 nm using a BioTek spectrophotometer.

For tamoxifen-sensitive studies with the ZR-75-1 tamoxifen-resistant cell line, cells were seeded with increasing tamoxifen concentrations (0–120 μg/mL). After 6 days, cell viability was evaluated by CCK-8 solution (Millipore-Sigma). Data were collected in biological duplicates run in triplicate for each time point.

To quantify apoptosis, transfected cells cultured for 48 h were harvested, digested, and pipetted to achieve a uniform single-cell suspension. Then cells were doubly stained with the Annexin V-FITC Apoptosis Detection Kit (Sigma-Aldrich). To score Annexin V and/or PI-positive cells, cells were subjected to analysis on a BD Accuri C6 Plus Flow Cytometer (Franklin Lakes, NJ, United States). All data were obtained from at least three independent culture experiments.

Hsa_circ_0086735, miR-1296-5p, and STAT1 mRNA expression levels were calculated from qRT-PCR CT values using 2^−ΔΔCT. The paired t-test was conducted for comparison of the expression differences between the tumor and normal tissues. Based on the mean value of hsa_circ_0086735 in the enrolled patients’ tissues, patients were classified into a low-expression group and a high-expression group. Chi-square test or Spearman coefficient of correlation was used to compare tumor characteristics across categories of hsa_circ_0086735 expression levels. The prognostic significance of relevant variables was assessed by Cox’s proportional-hazards regression from multivariate analysis. Overall and distant metastasis-free survival was assessed by Kaplan-Meier analysis with the log-rank test. Whereas normally distributed variables were compared using two-sided t-tests, non-normally distributed variables were compared using Mann-Whitney U tests. A p-value less than 0.05 was the cutpoint to be considered significant.

The flowchart of construction for the circRNA-miRNA-mRNA network was shown in Figure 1. To filter differentially expressed circRNAs, GSE18247 was analyzed by GEO2R (Figure 2A). To determine differentially expressed circRNAs in tamoxifen-resistant luminal breast cancer, the circRNAs obtained from GSE182471 were crossed and overlapped with those from GSE101410, GSE111504, and GSE159980. There were eight circRNAs significantly dysregulated in tamoxifen-resistant luminal breast cancer (Figure 2B). The differential expression of these eight circRNAs, hsa_circ_0086735, hsa_circ_0008267, hsa_circ_0001921, hsa_circ_0000458, hsa_circ_0008817, hsa_circ_0022601, hsa_circ_0004036, and hsa_circ_0006314, was presented as a heat map in Figure 2C.

The luminal-related miRNAs were obtained from GSE32022, GSE110204, and GSE121172, and presented as a Volcano diagram (Figure 3A). The downstream miRNAs of the eight dysregulated circRNAs were collected and overlapped with the luminal-subtype-related differentially expressed miRNAs (DE-miRNAs) (Figures 3B,C). In total, 26 miRNAs were classified as present across the panel of eight dysregulated circRNAs. These target miRNAs were used for further differential analysis.

The prognostic gene signature from GSE48390 for Han Chinese breast cancers was analyzed by GEO2R. A total of 508 differentially expressed and prognosis-related genes were identified (Figure 4A). Subsequently, 70 genes were screened to cover the downstream genes of the above 26 miRNAs, and their PPI was shown in Figure 4B with orphaned nodes deleted. The top 10 hub nodes in the PPI network were identified as PTPRC, IL7R, STAT1, LCK, SLAMF1, CCR5, ITK, CD38, CD3D, and CD24. We further conducted the GO and KEGG pathway enrichment analysis of 70 genes. The GO-enriched results showed that the genes mainly enriched in the main immune system process, developmental process, response to stimulus, and multicellular organismal process (Figure 4C). KEGG enrichment analysis showed that they were enriched for 124 pathways (Figure 4D). In KEGG, the top 20 enriched pathways included PD-L1 expression and PD-1 checkpoint pathway in cancer, cAMP signaling pathway, and JAK-STAT signaling pathway (Figure 4E).

A total of 96 circRNA-miRNA-mRNA axis, including 47 correctly sequenced ones (Figure 5A) and 49 invalid-ordered ones (Figure 5B), were obtained with the above-obtained RNAs. As for the hub genes (Supplementary Table S1), the circRNA-miRNA-mRNA-pathway network was mapped (Figure 5C). Given STAT1 was one of the hub genes that ranked second (Score = 10) in the interaction network by Degree method (Supplementary Table S1), we chose the hsa_circ_0086735-miR-1296-5p-STAT1 axis for further validation.

qRT-PCR was used to validate the expression levels of hsa_circ_0086735, miR-1296-5p, and STAT1 in luminal-subtype breast cancer. Hsa_circ_0086735 was significantly upregulated in luminal breast tumor tissues (p < 0.001, Figure 6A), which corresponded to the array data from GSE182471. In contrast, miR-1296-5p was significantly downregulated in the tumor tissues (p < 0.001, Figure 6B). STAT1 mRNA was an upregulated mRNA in luminal cancerous tissues (p < 0.001, Figure 6C). Luminal-subtype cell lines showed the same altered tends correspondingly (p < 0. 01, Figures 6D–F).

To confirm the interaction between hsa_circ_0086735 and miR-1296-5p, as well as miR-1296-5p and STAT1, Luciferase reporter assay and RNA pull-down assay, were introduced. The binding sites between hsa_circ_0086735 and miR-1296-5p were shown in Figure 7A, to gain the specific sequence of the combination of the two and thus to construct Luciferase report plasmids and sequence used in RNA pull-down assay. The expression of hsa_circ_0086735 was inversely correlated to that of miR-1296-5p (Figure 7B). Then the plasmids containing mutant- or wild-type hsa_circ_0086735 were used to study the influence of miR-1296-5p dysregulation on their luciferase activity. As expected, the reporter that carried binding sites was significantly influenced by miR-1296-5p mimics or inhibitors in MCF7 cells, whereas a luciferase reporter without binding was negligible (p < 0.001, Figure 7C). Then, biotin-labeled RNA pull-down assays demonstrated that biotin-labeled specific miR-1296-5p probes could pull down abundant hsa_circ_0086735 in ZR-75-1 cells (p < 0.001, Figure 7D). We simultaneously tested the binding of miR-1296-5p to STAT1 (Figure 7E). MiR-1296-5p was inversely expressed with STAT1 mRNA (Figure 7F). With Luciferase reporter assay, we assessed the binding of the wild-STAT1-containing reporter to miR-1296-5p and found that wild-type STAT1 resulted in repression of fluorescence intensity by miR-1296-5p mimics or induction by a miR-1296-5p inhibitor (p < 0.001, Figure 7G). RNA pull-down assay confirmed their binding (p < 0.001, Figure 7H).

To determine whether hsa_circ_0086735-miR-1296-5p-STAT1 axis was associated with luminal breast cancer cell function or patients’ survival, we assessed cell proliferation, cell apoptosis, and patients’ overall survival. Cell transfection was verified by the expression of STAT1 mRNA (Figures 8A,B). We observed a robust reduction in the proliferation of luminal breast cancer cells if hsa_circ_0086735 was individually inhibited, but increased again if both hsa_circ_0086735 and miR-1296-5p were inhibited until there was another reversal if hsa_circ_0086735, miR-1296-5p, and STAT1 were all inhibited (Figures 8C,D). The opposite trend occurs in cell apoptosis suggesting that hsa_circ_0086735 and STAT1 deficiency induced the apoptotic capacity, while miR-1296-5p inhibition reduced the apoptotic capacity (Figures 8E,F). The mean level of hsa_circ_0086735 (1.66) was defined as the cutoff for grouping of high expression and low expression. Next, high hsa_circ_0086735 was significantly correlated with patients’ histological type, tumor grade, molecular phenotype, and baseline Ki67 (Table 1). Additionally, high hsa_circ_0086735 was associated with poor overall survival in luminal breast cancer patients by Kaplan-Meier analysis (Figure 8G) and Multivariate Cox proportional hazards analysis (Table 2; HR 11.889, p = 0.007).

To determine whether hsa_circ_0086735-miR-1296-5p-STAT1 axis contributes to tamoxifen resistance, the area under the ROC curve (AUC) for STAT1, along with the other three genes in circRNA-miRNA-mRNA-pathway network, discriminating the non-response and response groups, was obtained from CTR-DB datasets (Supplementary Figure S1). Then, the expression of hsa_circ_0086735 was determined as an upregulated status in tamoxifen-resistant cell lines compared to that in the parental cell lines (Figure 9A). ZR-75-1 Tam1 cell line with a higher hsa_circ_0086735 expression was transfected for the subsequent experiments (Figure 9B). ZR-75-1 Tam1 developed a vulnerability to tamoxifen after hsa_circ_0086735 knockdown or hsa_circ_0086735/miR-1296-5p/STAT1 inhibition but become resistant to tamoxifen when hsa_circ_0086735 and miR-1296-5p were simultaneously inhibited, presented as the change in IC50 values (Figure 9C). Subsequent functional assays suggested that inhibited hsa_circ_0086735 reduced the cell proliferation potential, while miR-1296-5p can weaken this effect but STAT1 can consolidate it (Figure 9D). The induction of hsa_circ_0086735 on cell apoptosis rate can be abolished by miR-1296-5p and reproduced with the inhibition of STAT1 (Figure 9E). In patients treated with tamoxifen, high hsa_circ_0086735 was associated with poor distant metastasis-free survival (DMFS) (Figure 9F; Table 2, HR 6.945, p = 0.004).