ATRX is a predictive marker for endocrinotherapy and chemotherapy resistance in HER2-/HR+ breast cancer through the regulation of the AR, GLI3 and GATA2 transcriptional network

Function and pathway enrichment analysis of DEGs for drug resistance in BC

The workflow of this study is shown in Figure 1A. To investigate the differentially expressed genes (DEGs) involved in the drug resistance of BC, we downloaded and merged three relevant GEO datasets to obtain a GSE1 dataset. Moreover, on eliminating the batch effect between different sequencing platforms, a total of 1,461 DEGs were screened, wherein 818 DEGs were downregulated and 650 DEGs were upregulated (Figure 1B, 1C). Gene Ontology (GO) and Kyoto Encyclopaedia of Genes and Genomes (KEGG) analysis on DEGs (Figure 1D, 1E) revealed that the TOP10 biological process (BP) of GO was focused on transport, regulation of response to stimulus and cellular component assembly. KEGG revealed the enrichment of immune-related signalling pathways like PD-L1 expression and PD-1 checkpoint pathway in cancer, Toll-like receptor signalling pathway and T cell receptor signalling pathway. Moreover, the PI3K-Akt signalling pathway, PPAR signalling pathway and ABC transporters, which are known as drug resistance signalling pathways, were also enriched, thereby further validating our results.

Screening of hub genes using WGCNA and PPI

A cluster dendrogram of 1,461 DEGs was constructed using WGCNA based on the criteria of soft-threshold β = 9 and scale-free R² = 0.9 (Figure 2A, 2B). Subsequently, 14 gene modules were identified in the hierarchical clustering, based on a merge cut height of 0.25 and a minimum module size of 30 (Figure 2C). Among them, the brown module (correlation coefficient = -0.60, P = 7.88×10^-5) exhibited a significantly strong negative correlation with drug resistance. Therefore, the brown module was selected for further analysis as it contained 286 DEGs, which was the highest number of genes in a module, despite not having the strongest correlation with drug resistance (Figure 2D and Supplementary Table 2). Furthermore, 286 DEGs were selected for breast mammary tissue protein-protein interaction networks (PPI) analysis on a Network Analyst Database and 45 DEGs were identified using a Degree filter ≥ 15 (Figure 3A). Based on the intersection of DEGs from the three GEO datasets, we further identified 17 common DEGs that were downregulated and three common DEGs that were upregulated (Supplementary Figure 1A, 1B). Furthermore, ATRX and CDC5L were identified as hub genes using a Venn diagram of 45 DEGs from the brown module and 20 common DEGs (Figure 3B).

Exploring the role of ATRX in the drug resistance of BC

ATRX, a chromatin remodelling element, plays an important role in transcriptional regulation. Hence, an intersection analysis of ATRX-related genes in BC from TCGA was performed using the LinkedOmics database, Human TFs and DEGs in GSE1, consequently identifying 36 common DEGs (Figure 3C, 3D). PPI analysis highlighted that ATRX combined with AR, GLI3, EMX2, CUX1, FOXI1 and GATA2 formed a regulation circle (Figure 3E). Cox univariate analysis also revealed that, apart from EMX2 and CUX1, ATRX (Hazard Ratio (HR) = 0.063, 95% confidence interval (CI): 0.008-0.478) and FOXI1 (HR = 0.523, 95% CI: 0.285-0.963) showed a negative correlation with drug resistance, while AR (HR = 2.044, 95% CI:1.06-3.943), GLI3 (HR = 4.673, 95% CI:1.319-16.558) and GATA2 (HR = 2.229, 95% CI:1.032-4.812) showed a positive correlation with drug resistance (Figure 3F).

To further explore the role of these genes, we first detected their expression. The resistant samples revealed a reduction in ATRX and FOXI1 expression and an increase in AR, GLI3 and GATA2 expression (Figure 4A), which was consistent with the data of BC from TCGA. Notably, most of them had significant differences in both DFS and OS analysis (Supplementary Figure 2). Following this, Network Analyst Database analysed these genes in breast mammary tissue PPI, revealing a regulation network composed of ATRX, AR, GLI3 and GATA2 (Figure 4B). Furthermore, GO analysis for the network showed that the AR signalling pathway, negative regulation of transcription by RNA polymerase II and chromosome organization were enriched. KEGG analysis also revealed that the PI3K-Akt signalling pathway, Basal transcription factors, Endocrine resistance and Platinum drug resistance were enriched. Moreover, immunity-regulation signalling pathways were particularly enriched: such as the B cell receptor signalling pathway, Natural killer cell-mediated cytotoxicity and T cell receptor signalling pathway (Figure 4C).

As the above findings suggested the importance of immune regulation in drug resistance, the immunity analysis for GSE22513, which contained tissue-related data and contained more information about immunity than the one from the cell level, was considered for analysis. CIBERSORT analysis revealed that the proportion of 22 immune cells in each case, wherein naïve B cells, naïve CD8 T cells, CD4 T cells and activated Dendritic cells were lost or reduced while resting CD4 memory T cells and resting Dendritic cells were increased in the resistance group (Figure 5A). Furthermore, resting CD4 memory T cells and resting Mast cells were strongly upregulated while activated Dendritic cells were significantly suppressed in resistance samples, which validated that immunity-regulation was inhibited in the resistance group (Supplementary Figure 1C). Moreover, ssGSEA revealed a significantly lower immune status in the ATRX-Low group than in the ATRX-High group (Figure 5B). Furthermore, ATRX showed a negative correlation with resting mast cells and M0 Macrophages, indicating that elevated ATRX levels could promote immunity. While AR, GLI3 and GATA2 were negatively correlated with follicular helper T cells, activated Mast cells, CD8 and CD4 memory resting T cells, suggesting that elevated TFs may inhibit immunity (Figure 5C). Therefore, these findings indicate that the role of ATRX was opposite to AR, GLI3 and GATA2 in immunity.

Elucidating the potential mechanism of ATRX-mediated drug resistance in BC

To further explore the potential mechanism of drug resistance in BC, we explored the relationship between ATRX and the three TFs in the network. We predicted and combined the target genes from JASPAR, ENCODE, CHEA and the MotifMap database. Notably, ATRX was significantly negatively correlated with certain DEGs, which were the three TFs’ target genes (BPGM, DZIP1 and CORIN as AR’s target genes; CLGN, DNALI1 and LRRC6 as GATA2’s target genes; PDE4A, KIAA1462 and ADAMTS6 as GLI3’s target genes). However, ATRX revealed no relationship with the three TFs at the RNA level (Figure 6A–6C). Moreover, GSEA illustrated that the low expression of the target genes promoted immunity regulation (Figure 6D).

ATRX can bind to H3K9me3 to facilitate heterochromatin formation, which inhibits its transcription. Therefore, we compared H3K9me3 ChIP-seq data of normal breast epithelium tissues and MCF-7 from ENCODE Database, wherein both the signal intensities and the number of signal peaks on ATRX, AR, GLI3 and GATA2 were more abundant in normal tissues than that in MCF-7 (Figure 7). Furthermore, we determined the simple nucleotide variation data of BC samples from TCGA. Overall, missense mutations were the most common type. Moreover, a horizontal histogram revealed the frequent mutant genes in BC, namely TP53 (89.7%), ATRX (7.2%), BRCA2 (6.4%), BRCA1 (6.1%) and EGFR (3.4%) (Supplementary Figure 3A), indicating the downregulation of ATRX in BC drug resistance could be attributed to mutation. As ATRX can bind with AR directly, we predicted their potential interaction patterns, such as acetylation, ubiquitination, SUMOylation and so on (Supplementary Figure 3B). Moreover, decreased ATRX reduced the recruitment of H3K9me3 to AR, GATA2 and GLI3 to reduce the inhibition of transcription, resulting in the inhibition of immunity, promotion of the AR and PI3K-Akt signalling pathways, endocrine resistance and platinum drug resistance.

Validating the expression of ATRX using the TMA of BC

IHC staining validated ATRX expression on the TMA. A representative TMA stained for ATRX is shown in Figure 8A, 8B. ATRX expression was localised in the nucleus of tumour cells but was abundant in the 136 tissues of the TMA, which were defined as ATRX-High group (Figure 8A and Table 1). However, ATRX staining was low in the other 108 tissues, which were defined as ATRX-Low group (Figure 8B and Table 1). Interestingly, certain tissues appeared the same tendency of ATRX and TILs (Figure 8C–8E). Next, we analysed the relationship between ATRX expression and clinicopathological factors, revealing a significant correlation between ATRX expression and the ER status (P = 0.008), PR status (P = 0.024), histologic grade (P = 0.035), TILs (P = 0.045), chemotherapy (P = 0.041), OS (P = 0.013) and DFS (P = 0.033) (Table 1).

The mean follow-up time for OS and DFS were 66.93 months (range, 1–103 months) and 64.07 months (range, 1-103 months), respectively. During the follow-up period, 16.8% (41 of 244) of the patients had recurrence and/or metastasis whereas 13.5% (33 of 244) died. Univariate and Multivariate analyses of ATRX associated with OS in patients revealed that ER, T, N and combination ATRX with TILs significantly correlated with OS (Table 2). Similarly, PR, T and combination ATRX with TILs correlated significantly with DFS (Table 3). Survival analyses were performed using the Kaplan-Meier method with a log-rank test. The results showed that ATRX expression was positive with OS (P = 0.013) and DFS (P = 0.037) in patients with BC while the combination of ATRX with TILs showed a better correlation with both OS (P = 0.002) and DFS (P = 0.014). Furthermore, ATRX was significantly correlated with OS (P = 0.003) and DFS (P = 0.032) of patients with HER2-/HR+ BC. Additionally, the combination of ATRX with TILs showed a correlation with OS (P = 0.025) in patients with HER2-/HR+ BC (Figure 9).

In particular, the low expression of ATRX predicted poor DFS of the HER2-/HR+ BC subgroup who underwent endocrine (P = 0.034) or chemotherapy (P = 0.031) treatment. Similarly, suppressed ATRX revealed a bad OS for the same subgroup who received endocrine (P = 0.006) or chemotherapy (P = 0.001) treatment (Figure 10). Univariate analysis was performed for screening of risk factors, wherein combination ATRX with TILs (HR = 0.718, P = 0.016), PR (HR = 0.384, P = 0.007), ER (HR = 0.396, P = 0.006) were negatively correlated with OS while N (HR = 1.456, P = 0.004) was positively correlated with OS (Figure 8F).

Constructing a nomogram to predict survival in patients with BC

Nomogram survival prediction plots are commonly used to predict patient survival, with scores reflecting the values of several prognostic variables. Thus, we constructed a nomogram to estimate the probability of survival at 3, 5 and 7 years. The C-index value of the nomogram was 0.79. The calibration curves depicting the actual and nomogram-predicted survival at 3, 5 and 7 years were within the reference limits, reminding the nomogram based on our prognostic signature is precise and reliable (Figure 11). Moreover, receiver operating characteristics (ROC) analysis evaluated the predictive ability of ARTX in the TMA data. The AUC of ATRX for predicting the 3-year OS (AUC =0.842, P = 0.042) and the 5-year OS (AUC =0.848, P = 0.021) of the HER2-/HR+ BC patients who underwent endocrinotherapy were ideal (Supplementary Figure 4A). The similar good performance of ATRX was also found for the HER2-/HR+ BC patients who underwent chemotherapy (Supplementary Figure 4B).

Drug-gene interactions

In order to verify the role of ATRX in the first-line drugs sensitivity to HER2-/HR+ BC, the half maximal inhibitory concentration (IC₅₀) of PTX, DOX and TMX on MCF-7 cell and MCF-7 cell with lentivirus-mediated ATRX overexpression were detected. Westernblot showed when the multiplicity of infection (MOI) of lentivirus was 50 in MCF-7 cell, the expression of ATRX was increased significantly (Figure 12A). Therefore, we chose this MOI value for the following experiment. Interestingly, we found overexpression of ATRX significantly inhibited the IC₅₀ of the three first-line drugs on MCF-7 cell, especially the IC₅₀ of PTX and TMX dropped sharply below 50% (Figure 12B).

As ATRX, AR, GATA2 and GLI3 are vital to gene mediated drug resistance, the CTDbase database was used to obtain potential drugs. A total of four drugs were identified to promote ATRX and inhibit AR, GLI3 and GATA2. Among the identified drugs, bisphenol and tretinoin were considered ideal drugs (Figure 12C).

GEO data

As doxorubicin, paclitaxel and antiestrogen therapy are first-line anti-tumour drugs in BC, three microarray expression profiles on drug resistance in BC (GSE24460, GSE22513 and GSE8562) were downloaded from the Gene Expression Omnibus (GEO) database. The details are presented in Supplementary Table 1. The GSE24460 dataset (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE24460) consists of two MCF-7 doxorubicin Resistance cells vs two MCF-7 cells (Sensitive). In GSE22513 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE22513), 20 non-pCR breast biopsies (paclitaxel) vs eight pCR breast biopsies (paclitaxel) samples were present, and we defined non-pCR breast biopsy as the resistant group and pCR breast biopsy as the sensitive group. The GSE8562 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE8562) dataset consisted of three MCF-7/X-box binding protein 1 (XBP1) (Estrogen Resistance) cells vs three MCF-7 cells (Sensitive), where XBP1 is a well-known endocrine resistant gene in BC [41, 42]. Due to the limited sample size, we merged the three datasets into a singular large dataset (named GSE1) and corrected the batch effect using the SangerBox 3.0 platform.

Bioinformatics analysis

Limma analysis screened DEGs of the GEO data. The absolute value was identified |FC| > 1.2 and Benjamini–Hochberg-adjusted P < 0.05. GO and KEGG analyses were performed using the cluster profiler package in R in the SangerBox 3.0 platform. P < 0.05 was set as the cut-off criterion. Weighted Correlation Network Analysis (WGCNA) was performed in the SangerBox 3.0 platform. Additionally, interactive relationships and PPI networks of the DEGs were evaluated using the Network Analyst Database and STRING database. Human transcription factors (TFs) were downloaded from the human transcription factors database (http://humantfs.ccbr.utoronto.ca/allTFs.php). CIBERSORT algorithm analysed the proportion of 22 immune cells in resistant and sensitive samples, and the difference was compared using a t-test in the SangerBox 3.0 platform. Pearson’s correlation analysis calculated the correlation coefficient between immune cell infiltration and TFs. ssGSEA detected the different distribution of the 22 immunity cells between the ATRX-Low and ATRX-High groups using the SangerBox 3.0 platform. Furthermore, GSEA was used to associate genes with possible pathways in the SangerBox 3.0 platform. A false discovery rate < 0.5 and P < 0.05 were used as the criteria for statistical significance. Additionally, the trimethylation of histone H3 at lysine 9 (H3K9me3) ChIP-seq of normal breast epithelium tissue (female adult breast epithelium tissue) (https://www.encodeproject.org/experiments/ENCSR936LAH/) and MCF-7 (https://www.encodeproject.org/experiments/ENCSR999WHE/) in ENCODE were analysed, wherein ENCODE3 GRCh38 processed data were downloaded and normalised to compare the signal intensities and visualised using IGV.2.14.1 software. The mutation profiling in BC samples obtained from TCGA was presented in the SangerBox 3.0 platform. The potential interaction mechanism of AR and ATRX was predicted using the cBioPortal database. Drug-gene interactions were performed on the Network Analyst Database.

Patient information and tissue microarray

Human experiments were performed according to the ethical standards of the Helsinki Declaration and the China Ministry of Health’s ‘Ethical Review of Human Biomedical Research (Tentative, 2007)’. Written informed consent was obtained from all participants. Primary BC samples were obtained from the Hospital Affiliated with Nantong University between 2012 and 2017. We collected formalin-fixed paraffin-embedded surgery tissues from 244 patients with BC along with their complete clinicopathological data for preparation of tissue microarray (TMA). The diagnosis had given by two pathologists who were blinded to the clinicopathological data. Samples were classified into the following BC subtypes: HER2+/HR+, HER2+/HR-, HER2-/HR+ and TNBC.

Immunohistochemical (IHC) analysis

TMA was incubated with a primary antibody against a 1:200 dilution of ATRX (Rabbit-anti-human, Abcam) overnight at 4° C after the dewaxing and blocking of endogenous peroxides of the tissues. Visualisation of the antibody complex was achieved through a diaminobenzidine reaction, resulting in the brown staining of the cell nucleus. TMA was counterstained by Meyer’s haematoxylin. IHC staining was scored by pathologists based on the intensity of staining in the tumour cell nucleus from three hot spots of each tissue. The intensity of staining was scored as 0, negative; 1, weak; and 3, strong. We defined IHC staining score 0 and 1 as low ATRX expression and score 3 as high ATRX expression.

Assessment of tumour-infiltrating lymphocytes (TILs)

Hematoxylin eosin staining was used to stain the TMA. At x200 magnification, the percentage of the area containing TILs in the tumour nest and stroma to the total tumour area was calculated, and its mean value was defined as the TILs ratio. The proportions of TILs ≤ 5%, between 6% and 49% and ≥50% were considered as negative, weak and high TILs proportion groups, respectively. Moreover, we defined the negative and weak TILs groups as low TILs while the high TILs group stayed the same.

Western blot

The radio immunoprecipitation assay (Beyotime Biotechnology, Shanghai, China) lysis buffer supplemented with PMSF and protease inhibitor was used to extract total cellular protein. Protein was separated by 7.5% SDS-PAGE (Beyotime Biotechnology, China), transferred to polyvinylidene fluoride membrane, and sealed with 5% skim milk at room temperature. The membrane was sealed with diluted ATRX specific antibody (ab97508, Abcam, UK) at 4° C overnight, washed and incubated with secondary antibody at room temperature. After cleaning the membrane, enhanced chemiluminescence (Tanon, Shanghai, China) was used to detect protein bands.

Cell viability assay

To determine the IC₅₀, cell viability was assessed by CCK-8 (Vazyme, Nanjing, China) analysis. Cell suspensions of MCF-7 and MCF-7 with lentivirus-mediated ATRX overexpression, at a concentration of 5×10⁴ cell/100 μL, were seeded in a 96-well plate. The next day, cells were treated with different concentrations of DOX, PTX or TMX for 48 hours. Absorbances were assessed by detecting OD450 with an automatic spectrophotometer.

Statistical analysis

SPSS V17 software was used to analyse GEO data or experimental data. IC₅₀ was analysed by Graphpad Prism. Correlations between ATRX expression and clinicopathological characteristics were analysed using the χ² test. We defined OS as the time from the date of the primary surgery until the date of death or the last follow-up. DFS was defined as the time from the date of the primary surgery to the date of recurrence, which indicated locoregional recurrence, distant metastasis or death from any cause. Kaplan-Meier detected survival analysis, and the log-rank test estimated associations between variables and survival. Univariable and Multivariable Cox regression models identified significant prognostic factors. Validation of the nomogram was performed based on the primary group. The concordance index (C-index) assessed the discrimination of the nomogram. A calibration curve detected the calibration. Notably, the curve was corroborated with 1000 resamples conducted for validation. P < 0.05 was considered significant.