The human breast carcinoma cell lines MCF-7 (expressing ERα and wild-type p53) and MDA-MB-231 [ERα-negative, bearing mutated p53 (R280K)] were obtained from the European Collection of Cell Cultures (ECACC) and maintained in Dulbecco's modified Eagle's medium (DMEM) (Sigma-Aldrich; Merck KGaA) supplemented with 10% fetal bovine serum (FBS; Gibco; Thermo Fisher Scientific, Inc.) and 1% penicillin/streptomycin (Lonza Group, Ltd.) at 37°C in a humidified atmosphere containing 5% CO₂ until they reached 70% confluency. Where indicated, cells were treated with 10 ng/ml IFN-γ (Sigma-Aldrich; Merck KGaA) for 24 h or 10 µM etoposide (ETO) for 24 h (Sigma-Aldrich; Merck KGaA).

Concentration of 5 µM of the siGENOME PDIA1 siRNA and 5 µM of the siGENOME non-targeting siRNA pool was added to each well containing 2×10⁵ cells in DMEM and incubated for 72 h according to the suppliers' instructions (Dharmacon, UK) as described previously (4). The sequences of the siRNAs used are indicated as follows: siGE-NOME PDIA1-targeting siRNA pool: ACA GGA CGG UCA UUG AUU A, GGA CGG UCA UUG AUU ACA A, CCA AGAG UGU GUC UGA CUA, and CAGAGAGGAUCACAGAGUU; siGENOME non-targeting (scramble) siRNA pool: UAG CGA CUA AAC ACA UCA A, UAA GGC UAU GAA GAG AUA C, AUG UAU UGG CCU GUA UUA G, and AUG AAC GUG AAU UGC UCA A.

Cellular extracts from MCF-7 and MDA-MB-231 and 30 µg of total protein per sample were loaded on a 20% precast polyacrylamide gel and transferred to a PVDF membrane. The membranes were then incubated in 5% milk in PBS-0.1% Tween-20 (v/v) with anti-P4HB anti-body (Santa Cruz Biotechnology, sc-74551) (dilution 1:500) or β-actin (Sigma-Aldrich; Merck KGaA; A1978) (dilution 1:10,000) overnight at 4°C. After incubation with secondary anti-mouse immunoglobulin G conjugated to horseradish peroxidase (GE Healthcare) (dilution 1:1,000) in 2.5% milk in PBS-0.1% Tween-20 (v/v) for 1 h at 25°C, the protein bands were visualized using the ChemiDoc MP imaging system (Bio-Rad Laboratories) The blots were quantified using ImageJ version 1.51 (National Institutes of Health).

MCF-7 and MDA-MB-231 breast cancer cell lines were seeded in 6-well plates at a concentration of 3×10⁵ cells per well for 24 h after treatment. RNA was extracted according to the RNeasy Mini Kit supplier instructions (Qiagen, UK). After collecting RNA, the samples were maintained at −80°C before sequencing carried out by the Next Generation Sequencing Facility (University of Leeds, UK).

Initial alignment of RNA transcript counts, sequence short reads, and normalization were carried out by the Next Generation Sequencing Facility (University of Leeds, UK). Sequence data in Fastq format were quality-checked using FastQC software (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Cutadapt version 1.16 software (https://cutadapt.readthedocs.org/en/stable/) was used to trim poor quality bases (Phred quality score <20) and contaminating adapter sequences from raw reads. Reads trimmed to fewer than 30 were discarded. Reads were aligned to a human (hg38) genome reference sequences, obtained from the UCSC database (41) using the splicing-aware STAR aligner (42). STAR aligner was run with known splice junctions supplied in GTF file format, obtained from the UCSC database using Table Browser tool (43,44). The resulting alignments in BAM file format were checked for quality using QualiMap software (44) and Picard tools version 1.90 (http://picard.sourceforge.net). Picard tools also used to mark PCR/Optical duplicate alignments. BAM files were indexed using Samtools software (45) and visualised using IGV browser (46) to check for genomic DNA contamination and the presence of PCR duplicates. Bioconductor R package RSubread (47) was used to extract raw sequenced fragment counts per transcript using the RefSeq hg38 annotation dataset used by STAR during alignment. Multi-mapping read pairs were counted as a fraction of all equivalent alignments. Read count data were generated with the inclusion of reads marked as PCR/optical duplicates.

Significantly downregulated or upregulated genes (Benjamini-Hochberg <0.01) in the PDIA1-silenced MCF-7 and MDA-MB-231 treated with IFN-γ or ETO compared to MCF-7 and MDA-MB-231 cells expressing PDIA1 treated with IFN-γ (39) or ETO (40) were selected. To visualize the overall changes in gene expression between PDIA1-silenced and scramble-transfected MCF-7 and MDA-MB-231 cells, the log value of the change in gene expression was graphed against the normalized read count for each transcript using the plotMA function of the DeSeq2 R package (47). The read count data required for each pairwise analysis was imported into the R package DeSeq2 and the effect size of each sample's library was defined and used to normalize the read count data.

To understand the clinical importance of the genes identified and filtered from the RNA-seq data, the Kaplan-Meier Plotter website (http://kmplot.com/analysis/) was utilized. For this purpose, Kaplan-Meier (KM) plots for ERα-positive and ERα-negative patients were generated for the genes identified to be upregulated or downregulated in the RNA-seq data (48). The following parameters were selected: Survival: OS, Split patients by: upper quartile, Follow-up threshold: 240 months, Probe set option: only JetSet best probe sets (49), ER status-IHC: ER-positive (for the genes upregulated or down regulated in MCF-7 cells) or ER negative (for the genes upregulated or downregulated in MDA-MB-231 cells) ER status-array: ER positive (for the genes upregulated or downregulated in MCF-7 cells) and ER negative (for the genes upregulated or down regulated in MDA-MB-231 cells). Log rank P-values <0.05 for the KM plots of all genes was considered statistically significant.

The Pearson correlation coefficient for the comparison of PDIA1 mRNA levels with those of the genes upregulated or downregulated in the PDIA1-silenced MCF-7 cells treated with etoposide (ETO) or interferon γ (IFN-γ) in breast cancer patients was explored using the cBioPortal platform (http://cbioportal.org) in combination with data provided in the Breast Cancer (SMC 2018) database using the mRNA expression z-scores relative to all sample (RNA seq TPM-168 samples) options.

In order to understand the biological importance of the gene lists generated from the previous section, the GeneCards tool (50) was used to identify linked genes to the genes listed in Tables I and II. After this, Metascape tool (51) was employed. Gene lists were entered in the metascape tool using the Official Gene Symbol as the identifier and Homo sapiens as the species.

Following the identification of genes important for breast cancer patient survival, the potential clinical targeting of each gene was explored via DRUGSURV (52). Genes exerting oncogenic function that were identified as statistically significant for patient OS via KM Plotter were queried in turn, and the approved drugs that targeted the genes directly or indirectly were recorded.

Total mRNA obtained from MCF-7 (Fig. 1A and B) or MDA-MB-231 (Fig. 1C and D) cells transfected with either scramble or specific siRNA to silence PDIA1 gene expression and treated with either ETO (Fig. 1A and C) or IFN-γ (Fig. 1B and D) was submitted to RNA-seq analysis. The transfection efficiency for the PDIA1-siRNA transfected MCF-7 cells compared to the transfection efficiency for the scramble siRNA was 67.13% (Fig. 1E). The transfection efficiency for the PDIA1-siRNA transfected MDA-MB-231 cells compared to the transfection efficiency for the scramble siRNA was 96.67% (Fig. 1F).

The total number of genes exhibiting statistically significant downregulation in ETO-treated MCF-7 cells in which PDIA1 gene expression had been silenced compared with ETO-treated MCF-7 cells transfected with scramble siRNA was 841, and the total number of genes upregulated under these conditions was 863 (Fig. 2A). The total number of genes exhibiting statistically significant (P-adj<0.01) downregulation in IFN-γ-treated MCF-7 cells in which PDIA1 gene expression had been silenced compared with IFN-γ-treated MCF-7 cells transfected with scramble siRNA was 773, and the total number of genes upregulated under these conditions was 919 (Fig. 2B). The total number of genes exhibiting statistically significant downregulation in ETO-treated MDA-MB-231 cells in which PDIA1 gene expression had been silenced compared with ETO-treated MDA-MB-231 cells transfected with scramble siRNA was 55, whereas the total number of genes upregulated in these cells under the same conditions was 27 (Fig. 2C). In the IFN-γ-treated MDA-MB-231 cells in which PDIA1 gene expression had been silenced compared with IFN-γ-treated MDA-MB-231 cells transfected with scramble siRNA the total number of the downregulated genes was 11, whereas the total number of genes upregulated in these cells under the same conditions was 27 (Fig. 2D). The numbers of upregulated and downregulated genes are presented in Fig. 2E and F, respectively.

To shed light on the clinical significance of the aforementioned findings, KM curves for the upregulated or downregulated genes in the PDIA1-silenced compared with the PDIA1-expressing and either IFN-γ- or ETO-treated MCF-7 and MDA-MB-231 cells were plotted in the patients with ERα-positive and ERα-negative breast cancer, respectively, as described in the summary of the experimental design (Fig. 3). The upregulated or downregulated genes in MCF-7 cells under the conditions described above that exhibited statistically significant (P<0.05) KM results are shown in Table I. The KM curves plotted for these genes for ERα-positive patients are shown in Fig. 4A-E. High expression of aurora kinase-A (AURKA), mitotic checkpoint serine/threonine-protein kinase BUB1 (BUB1), cyclin-B2 (CCNB2), cell division cycle 25C (CDC25A) (Fig. 4A), transcription factor E2F8 (E2F8), junD proto-oncogene (JUND), kinesin family member 11 (KIF11), L1 cell adhesion molecule (L1CAM) (Fig. 4B), melanoma antigen gene D1 (MAGED1), midkine (MDK), pleckstrin homology-like domain family A member 2 (PHLDA2), membrane-associated tyrosine- and threonine kinase 1 (PKMYT1) (Fig. 4C), peroxiredoxin-2 (PRDX2), protein tyrosine phosphatase receptor type S (PTPRS), pyrroline-5-carboxylate reductase 1 (PYCR1) (Fig. 4D), uncoupling protein 2 (UCP2), SMAD family member 6 (SMAD6), transketolase (TKT) and tubulin β 2A class IIa (TUBB2A) (Fig. 4E) genes was associated with statistically significant (P<0.05) shorter OS probability in the patients with ERα-positive breast cancer. On the other hand, high expression of the S100P-binding protein (S100PBP) gene was associated with statistically significant (P<0.05) longer OS probability when comparing the 75th quartile with the 25th quartile in KM plotter in patients with ERα-positive breast cancer (Fig. 4D).

To explore whether PDIA1 gene expression in patients with ERα-positive breast cancer was correlated with the upregulated or downregulated genes in either IFN-γ- or ETO-treated PDIA1-silenced MCF-7 cells, the cBioPortal bioinformatics tool was used to investigate the Pearson correlation coefficient between PDIA1 mRNA levels and that of the genes shown in Table I. This investigation was carried out by analyzing the data provided in the Breast Cancer (SMC 2018) database using the mRNA expression z-scores relative to all samples (RNA seq TPM, 168 samples) option. Positive correlations were observed between PDIA1 gene expression and that of AURKA, CDC25A, centromere protein F (CENPF), cytoskeleton-associated protein 5 (CKAP5), E2F8, extra spindle pole bodies like 1, separase (ESPL1), flap structure-specific endonuclease 1 (FEN1), G2 and S phase-expressed protein 1 (GTSE1), JUND, KIF11, kinesin-like protein KIF20A (KIF20A), MAGED1, MDK, PKMYT1, PRDX2, PTPRS, PYCR1, RAN binding protein 1 (RANBP1), S100PBP, SMAD6, transforming acidic coiled-coil-containing protein 3 (TACC3), TKT, thyroid hormone receptor interactor 13 (TRIP13), tubulin b class 1 (TUBB), ubiquitin-conjugating enzyme E2 C (UBE2C) and UCP2 genes. In other words, in the presence of PDIA1, the expression of these genes would be upregulated, whereas in the absence of PDIA1, the expression of these genes would be downregulated in patients with ERα-positive breast cancer. Among these genes, the correlation between PDIA1 (shown as P4HB) gene expression and that of the CKAP5, FEN1, PYCR1, TACC3, TKT and TUBB genes exhibited higher Pearson correlation coefficients (Fig. 5).

Positive correlations between PDIA1 gene expression and that of MAGED1, MDK, PRDX2, PTPRS, PYCR1, SMAD6, TKT and UCP2 genes were found, which were demonstrated to be upregulated in the PDIA1-silenced MCF-7 cells (Table I), thus indicating that in the presence of PDIA1 in patients with ERα-positive breast cancer these genes would be downregulated. According to the KM plots, the downregulation of MAGED1, MDK, PRDX2, PTPRS, PYCR1, SMAD6, TKT and UCP2 mRNA levels was associated with higher OS probability in the patients with ERα-positive breast cancer (Fig. 4). This suggested that PDIA1 can associate with MAGED1, MDK, PRDX2, PTPRS, PYCR1, SMAD6, TKT and UCP2 to induce anti-oncogenic effects in patients with ERα-positive breast cancer. To explore the pathways involved in the PDIA1-mediated anti-oncogenic effects in patients with ERα-positive breast cancer, the Metascape bioinformatics tool was employed (51). The results of this analysis showed that the pathways by which PDIA1 exerts tumor-suppressive effects include cell-substrate adhesion, regulation of growth and response to ROS (Table III).

Direct correlations between PDIA1 gene expression and that of the AURKA, CDC25A, CENPF, CKAP5, E2F8, ESPL1, FEN1, GTSE1, JUND, KIF11, KIF20A, PKMYT1, RANBP1, S100PBP, TACC3, TRIP13, TUBB and UBE2C genes were also observed, which were found to be downregulated in the PDIA1-silenced MCF-7 cells (Table I). This suggested that in the presence of PDIA1 in patients with ERα-positive breast cancer these genes would be upregulated. The KM plots indicated that high levels of these genes were associated with shorter OS in patients with ERα-breast cancer (Fig. 4A-C and E), apart from S100PBP (Fig. 4D). Therefore, the association of PDIA1 with these genes is a potential route by which PDIA1 exerts oncogenic effects in patients with ERα-positive breast cancer. The pathways affected by the networks formed by these genes were investigated using the Metascape bioinformatics tool and they included control of cell cycle transition to the mitotic phase as well as regulation of the cell cycle checkpoints (Table IV).

Similar analysis to that described for the MCF-7 cells was carried out for the upregulated or downregulated genes in the IFN-γ- or ETO-treated MDA-MB-231 cells with PDIA1 silenced compared with IFN-γ- or ETO-treated cells transfected with scramble siRNA that exhibited statistically significant (P<0.05) KM results are shown in Table II, and the KM curves plotted for these genes in ERα-negative patients are presented in Fig. 6. KM survival curves were plotted for the upregulated or downregulated genes in MDA-MB-231 cells under various oxidative stress conditions in patients with ERα-negative breast cancer (Fig. 6A and B). High expression of the coiled-coil-helix -coiled-coil-helix domain containing 3 (CHCHD3), heat shock protein 90 β family member 1 (HSP90B1), lysosomal-associated transmembrane protein 5 (LAPTM5) (Fig. 6A), protein tyrosine phosphatase receptor type J (PTPRJ), ring finger protein 213 (RNF213), spectrin repeat containing nuclear envelope protein 1 (SYNE1), vacuolar ATPase assembly factor VMA21 (VMA21) and X-linked inhibitor of apoptosis (XIAP) (Fig. 6B) genes was associated with statistically significant (P<0.05) longer OS in the patients with ERα-negative breast cancer. On the other hand, high gene expression of insulin-like growth factor-binding protein 3 (IGFBP3) was associated with lower OS probability in these patients (Fig. 6A).

The Pearson correlation coefficients between PDIA1 gene expression and the upregulate or downregulated genes in the PDIA1-silenced MDA-MB-231 cells under oxidative stress conditions (Table II) were investigated using the cBioPortal platform. To carry out this analysis, only samples originating from patients with ERα-negative breast cancer were included, but no statistical significance or data were available to manually calculate the correlation between PDIA1 gene expression and these genes in the patients with ERα-negative breast cancer. Thus, it was hypothesized that the Pearson correlation coefficients between PDIA1 gene expression with that of the genes shown in Table II were positive, meaning that in the presence of PDIA1, HSP90B1, IGFBP3, LAPTM5, RNF213 and SYNE1 genes would be downregulated, while CHCHD3, PTPRJ, VMA21 and XIAP would be upregulated, in patients with ERα-negative breast cancer.

Direct associations between PDIA1 and CHCHD3, PTPRJ, VMA21 and XIAP indicated that these genes can induce tumor-suppressive effects in patients with ERα-negative breast cancer since the KM plots in these patients showed that the high expression of these genes was associated with longer OS (Fig. 6A and B). PDIA1 would mediate tumor-suppressive effects in patients with ERα-negative breast cancer patients via IGFBP3 as well since this gene in the presence of PDIA1 is downregulated (Table II) and the KM plot in these patients showed that low expression of IGFBP3 was associated with longer OS (Fig. 6A). In addition, direct association between PDIA1 and HSP90B1, LAPTM5, RNF213 and SYNE1 gene expression suggested that these genes confer oncogenic effects to patients with ERα-negative breast cancer since the KM plots for these patients indicated that low expression was associated with lower OS probability (Fig. 6A and B).

The Metascape bioinformatics tool was used to identify the pathways affected by the genes modulated by PDIA1 to confer anti- or pro-carcinogenic effects in patients with ERα-negative breast cancer. This analysis revealed that PDIA1 may exert tumor-suppressive effects in patients with ERα-negative breast cancer by interfering with the NF-kB signaling pathway (XIAP), glycolysis (IGFBP3), mitochondrial biogenesis (CHCHD3) and negative regulation of cell migration (IGFBP3, P4HB, PTPRJ) (Table V). Metascape analysis also showed that PDIA1 may trigger pro-oncogenic effects by forming signaling networks with calmodulin-3 (CALM3), HSP90B1, LAPTM5 and RNF213 genes, which could affect pathways such as the adaptive immune system and cytokine signalling in the immune system, as well as activating signal transduction in the immune response (Table VI).

To investigate whether approved or experimental drugs targeting the genes involved in the pathways conferring oncogenic activity in the patients with ERα-positive (Table VII) or ERα-negative (Table VIII) breast cancer, the genes associated with PDIA1 and its oncogenic effects were investigated in each subgroup of patients with breast cancer using the DRUGSURV (52) bioinformatics tool. This investigation revealed that in the patients with ERα-positive breast cancer, drugs targeting the function of the AURKA, FEN1, KIF11, PKMYT1 and TUBB genes, which are involved in the regulation of cell cycle progression, would be selectively effective for this group of patients. On the other hand, drugs targeting the function of the HSP90B1 gene, which is involved in the regulation of toll-like receptor translocation through the endoplasmic reticulum to the plasma membrane (53), would be preferential for the patients with ERα-negative breast cancer.