Datasets in the GEO meeting the following criteria were included in our study and formed the discover and validation datasets. First, gene expression data of patients was acquired using the Affymetrix HG-U133 A (GPL96) or HG-U133 Plus 2.0 (GPL570) microarray platform (Affymetrix, Santa Clara, USA). Next, breast cancer patients were with a full record of their ER status and treated with tamoxifen adjuvant treatment. Finally, distant metastasis-free survival (DMFS) records were available for each patient in these datasets. Since all the patients' data were obtained from public available GEO, ethic approvals of study and informed consent were already handled when they were submitted to GEO.

The lncRNA expression data from gene expression chips was obtained and analyzed as previously reported [111213]. Briefly, the Robust Multichip Average package and Combat function in the Surrogate Variable Analysis package were utilized to normalize raw data among different datasets [1415]. Moreover, the Guided Principal Components Analysis package was utilized to assess the batch effect before and after normalization [16]. The statistical variance and the first two principal components from each batch were compared before and after normalization for the discover and validation datasets, respectively. The p-value of batch variance was analyzed accordingly.

To determine lncRNA expression, probes from the arrays were aligned to the human genome (GRCh38/hg38) using SeqMap [17] such that probes matching the lncRNA chromosomal positions from GENCODE (http://www.gencodegenes.org; GRCh38, release 25) were identified [18].

To identify lncRNAs correlated with DMFS, univariate Cox proportional hazards regression analysis was firstly performed to evaluate the relationship between each lncRNA and DMFS in the discover cohort. Only lncRNAs related to DMFS with p<0.002 were considered statistically significant. Each lncRNA was evaluated by stratifying patients in the discover cohort with Cutoff Finder [19]. Multivariate Cox proportional hazards regression analysis was conducted by inputting the significant lncRNAs. By integrating these prognosis-related lncRNAs, a Cox proportional hazards regression prediction model was established. The predictive risk score of each patient was then calculated from the linear combination of lncRNA expression with its regression coefficients generated in multivariate Cox proportional hazards regression analysis. Separated by the optimized cutoff score value, patients were classified into a low-risk group with relatively good survival and a high-risk group with poor survival. Survival curves were derived using the Kaplan-Meier method with log-rank tests to evaluate differences in DMFS between the low- and high-risk groups with R package “survminer”. Receiver operating characteristic (ROC) analysis was performed to assess the predictive capability of the model with the R package “survivalROC” [20].

To evaluate the biological processes driving tamoxifen resistance in the tamoxifen efficacy-related lncRNA signature (TLS) high-risk patients, Gene Set Enrichment Analysis (GSEA; http://www.broadinstitute.org/gsea) was performed to detect the biological differences between TLS high-risk and low-risk patients using MSigDB c2: curated gene sets: all canonical pathways [2122]. Gene sets associated with TLS high-risk patients and identified with a false discovery rate (FDR)<0.01 and p<0.005 were considered as statistically significant. Furthermore, we attempted to annotate the potential function of each lncRNA in the TLS. Because GSEA failed to show significant results for individual lncRNAs in our analysis, a previously described method was adopted [111323]. Briefly, mRNAs highly correlated with each lncRNA were identified in the discover dataset by Pearson correlation analysis (top 1.0%). Next, the positively or negatively correlated mRNAs were input into the Database for Annotation, Visualization, and Integrated Discovery (DAVID; version 6.8; https://david.ncifcrf.gov/) [2425]. Finally, DAVID functional annotations with FDR<0.01 and p<0.005 were visualized with the Enrichment Map plugin in Cytoscape (version 3.2.0; http://www.cytoscape.org) for each lncRNA [26].

All data in this study was analyzed with R software (version 3.3.1; http://www.r-project.org/) and Bioconductor (http://www.bioconductor.org/). PAM50 intrinsic subtypes, 21-gene recurrence score (Oncotype DX®; Genomic Health, Redwood City, USA) and 70-gene prognosis signature (MammaPrint®; Agendia, Amsterdam, the Netherlands) were obtained by using the “genefu” package in R [27]. A p-value less than 0.05 was considered significant.

A total of 1,056 tamoxifen-treated patients with ER-positive breast cancer from GSE6532, GSE9195, GSE17705, GSE19615, GSE26971, and GSE45255 datasets were enrolled in this study. In detail, 197 ER-positive patients from GSE6532 (GPL96 platform part) and 298 ER-positive patients from GSE17705 (GPL96) treated only with 5 years of tamoxifen were combined as the discovery dataset for lncRNA-based model construction. Additionally, 88 patients from GSE6532 (GPL570 platform part), 77 patients from GSE9195 (GPL570), 62 patients from GSE19615 (GPL96), 258 patients from GSE26971 (GPL96), and 74 patients from GSE45255 (GPL96) who were also treated with 5 years of tamoxifen were combined as the validation dataset. Principal component analysis revealed no significant variance among batches for both the discover and validation datasets after normalization (Supplementary Figures 1 and 2, available online).

In the discover cohort, there were 315 and 164 patients with negative and positive lymph nodes. A total of 147 patients were classified as luminal A subtype, while 348 patients had other intrinsic subtypes of tumors. For the 21-gene recurrence score, 253 tumors were classified as recurrence score low and medium, while 242 patients had high recurrence score tumors. In terms of the 70-gene prognosis signature, there were 48 and 447 patients with low-risk and high-risk signatures (Table 1). The medium follow-up period was 7.4 (0.0–16.3) years and 126 patients (25.5%) experienced distant metastasis in the discover cohort.

A total of 558 patients were included in the validation cohort. Two hundred thirty-nine patients had tumors no larger than 2.0 cm, while 308 patients had tumors larger than 2.0 cm. A total of 285 patients were lymph node-negative and 231 patients were lymph node-positive. Additionally, 180 patients were classified as luminal A subtype and 378 as other subtypes. In terms of 21-gene recurrence classification, there were 244 and 314 patients classified as having low to medium and high recurrence scores, respectively. Meanwhile, 94 and 464 patients had low and high 70-gene prognosis scores, respectively (Table 2). After a medium follow-up of 6.3 (0.0–17.6) years, 117 patients (21.0%) had distant metastasis events.

After establishing the discover and validation datasets, we re-annotated the probes corresponding to lncRNAs for the HG-U133A and HG-U133 Plus 2.0 Affymetrix platform. A total of 598 probes corresponding to 452 lncRNAs were obtained for the HG-U133A microarray, while 5,654 probes were matching with 3,793 lncRNAs in the HG-U133 Plus 2.0 microarray. There were 365 lncRNAs overlapping between the two platforms. Figure 1 shows a diagram of the data analysis and model construction.

In the discover group of 495 ER-positive patients treated with tamoxifen, six of these 365 overlapping lncRNAs were significantly associated with DMFS in univariate Cox proportional hazard regression analysis (Table 3). Each of these six lncRNAs was capable of classifying patients into high- and low-risk groups, which could predict disease outcomes in this cohort of patients (Figure 2).

These six lncRNAs were then integrated in a multivariate Cox proportional hazard regression model to construct TLS, which assessed each tamoxifen-treated patient with an individual risk score. TLS scores were calculated as follows: TLS score=0.3753×expression value of RP11-189B4.7+0.1399× expression value of RP11-59H7.3+0.3604×expression value of CTD-2090|13.1+0.1562×expression value of LINC01399+ 0.4007×expression value of RP11-119F7.5+0.0561×expression value of RP11-193F5.1. The ROC curve of TLS was subsequently obtained with an area under the ROC curve (AUC) of 0.764 for 5 years of DMFS in the discover cohort (Figure 3A). In the ROC curve, by selecting the point nearest to the perfect prediction point, the cutoff value was set to 0.362 for TLS. Patients in the discover dataset were then separated into 133 high-risk cases and 362 low-risk cases (Table 1). The heat map revealed distinct lncRNA expression and DMFS in the TLS high- and low-risk groups (Figure 3B). Moreover, patients in the high-risk group had worse DMFS than those in the low-risk group by Kaplan-Meier analysis (Figure 3C). Univariate Cox proportional hazard regression demonstrated that lymph node status (hazard ratio [HR], 1.99; 95% confidence interval [CI], 0.62–6.37; p=0.001), TLS (HR, 4.04; 95% CI, 2.83–5.77; p<0.001), intrinsic subtype (HR, 1.68; 95% CI, 1.10–2.58; p=0.016), 21-gene recurrence score (HR, 2.27; 95% CI, 1.56–3.29; p<0.001) and 70-gene prognosis signature (HR, 2.63; 95% CI, 1.08–6.44; p=0.028) were related to DMFS in the discover cohort. Multivariate Cox proportional hazard regression analysis demonstrated that only TLS (HR, 3.40; 95% CI, 2.34–4.96; p<0.001) and 21-gene recurrence score (HR, 1.72; 95% CI, 1.15–2.58; p=0.008) were independent factors associated with DMFS in the discover cohort.

In the validation cohort, TLS separated patients into 178 high-risk cases and 380 low-risk cases (Table 2). This six lncRNA-based signature was also found to be associated with DMFS (Figure 4). In univariate Cox proportional hazard regression analysis, we found that the TLS (HR, 1.66; 95% CI, 1.15–2.40; p=0.006), 21-gene recurrence score (HR, 1.49; 95% CI, 1.02–2.17; p=0.038), and 70-gene prognosis signature (HR, 1.90; 95% CI, 1.05–3.45; p=0.032) were related to DMFS. Multivariate analysis showed that only TLS (HR, 1.54; 95% CI, 1.06–2.24; p=0.023) was independently correlated with DMFS.

Overall, a total of 1,053 patients were in the discover and validation cohorts. There were 600, 395, and 58 patients with lymph node-negative, -positive, and unknown disease. Patients with known lymph node status were included in the following analysis. According to the PAM50 classification criteria, 327 cases were identified as luminal A subtype, 496 cases as luminal B subtype, 61 cases as HER2-enriched subtype, 44 cases as basal-like subtype and 125 cases as normal-like subtype. Patients with luminal A and B disease were included in further subgroup analysis. In terms of the 21-gene recurrence score, there were 232, 265, and 556 patients classified as having low-, medium-, and high-risk recurrence scores, respectively. For the 70-gene prognosis signature, 142 and 911 patients were classified as having low-risk and high-risk signatures. Subgroups analysis showed that TLS performed better in the lymph node-negative subgroup (HR, 3.28; 95% CI, 2.26–4.75; p<0.001) than in the positive subgroup (HR, 2.02; 95% CI, 1.39–2.93; p<0.001), better in the luminal B subgroup (HR, 2.41; 95% CI, 1.73–3.37; p<0.001) than in the luminal A subgroup (HR, 1.96; 95% CI, 1.07–3.58; p=0.026), better in the 21-gene recurrence score high-risk subgroup (HR, 2.60; 95% CI, 1.89–3.56; p<0.001) than in the low-risk subgroup (HR, 1.91; 95% CI, 0.96–3.83; p=0.062) and medium-risk subgroup (HR, 1.85; 95% CI, 1.00–3.44; p=0.047) and better in the 70-gene prognosis signature high-risk subgroup (HR, 2.35; 95% CI, 1.81–3.06; p<0.001) than in the low-risk subgroup (HR, 2.55; 95% CI, 0.73–8.92; p=0.129) (Figure 5).

The predictive accuracy of TLS was then compared with the 21-gene recurrence score and 70-gene prognosis signature in all patients with AUC values of 0.656, 0.635, and 0.631, respectively (Figure 6). In all lymph node-negative patients, the AUC values were 0.693, 0.676, and 0.623 for the TLS, 21-gene recurrence score and 70-gene prognosis signature in terms of DMFS outcome prediction. Similar results were observed in the lymph node-positive subgroup (AUC values of 0.627, 0.573, and 0.611, respectively). For different intrinsic subtypes of all patients, TLS was not superior to the 21-gene recurrence score and 70-gene prognosis signature in the luminal A subtype (AUC values of 0.548, 0.688, and 0.534, respectively). However, TLS performed much better than the 21-gene recurrence score and 70-gene prognosis signature in the luminal B subtype with AUC values of 0.659, 0.580, and 0.580, respectively. Furthermore, integration of TLS into the current 21-gene recurrence score and 70-gene prognosis signature improved the disease outcome prediction power (Figure 6).

GSEA was performed to determine the biological difference between the TLS high-risk and low-risk group patients and biological processes driving tamoxifen resistance in TLS high-risk patients. In GSEA, gene sets significantly associated with TLS high-risk patients were identified in the discover cohort (FDR<0.01, p<0.005) (Supplementary Table 1, available online). Our result indicated that most gene sets associated with TLS high-risk patients were mainly correlated with cell cycle (most gene sets were correlated with G0/G1 phase, some were related to S phase and G2/M phase) and nucleotide metabolism (such as transcription, DNA and RNA synthesis). We also attempted to annotate the function of each lncRNA in the TLS by DAVID annotation analysis whose results are shown in Supplementary Figure 3 (available online).