Identification of PIMREG as a novel prognostic signature in breast cancer via integrated bioinformatics analysis and experimental validation

Analysis of PIMREG expression level

The gene expression data of 1109 BC tissues and 113 adjacent normal tissues were obtained from the Breast Invasive Carcinoma Project (TCGA-BRCA) of the TCGA database. The clinicopathological characteristics of the patient’s data were summarized. Wilcoxon rank-sum test was used to assess the expression level of PIMREG between BC and adjacent normal tissue; visualization was performed using the ‘ggplot2’ package in R software v3.6.3 (R Core Team, 2020; Ding et al., 2022). Further analysis was conducted on information obtained from 112 pairs of breast cancer tissue samples and adjacent normal breast tissues. Moreover, subgroup analyses were conducted based on T, N, M stage, pathologic stage, PR status, HER-2 status, and ER status. Statistical significance was considered when p <0.05. Findings were verified through the GSE42568 data set, which including104 BC tissues and 17 normal tissues. Data obtained from Gene Expression Omnibus (GEO) database set aimed to analyze the expression difference of PIMREG between BC tissues and normal breast tissue (Clarke et al., 2013; Zhou et al., 2020). All figures were visualized using the R package ggplot2. To evaluate and verify the expression of PIMREG in BC, images of immunohistochemical (IHC) staining for the PIMREG protein obtained from the Human Protein Atlas (HPA) (http://www.proteinatlas.org/) were used (Pontén et al., 2011; Zhu et al., 2022).

Prognosis analysis of PIMREG expression within BC

Using RNAseq data with Log2 transformed TPM (transcripts per million reads) of TCGA-BRCA, the prognostic significance of PIMREG expression was first investigated. To compare the survival difference, as well as overall survival (OS), disease-specific survival (DSS), and progression-free interval (PFI), between high PIMREG and low PIMREG expression groups, we applied the Kaplan–Meier (KM) survival analysis with the log-rank test. Log-rank tests and univariate Cox proportional-hazards regression were performed to create the KM curves, which included p-values and hazard ratio (HR) with a 95% confidence interval (CI). The KM curve of OS based on the GEO database (GSE42568) was also plotted with the same method. R software with ‘survminer’ and ‘survival’ packages was used to perform these calculations. Furthermore, using the “plotROC” package, the diagnostic property was verified through receiver operating characteristic (ROC) analysis of BC. Meanwhile, in order to explore and further clarify the risk factors in patients with BC, the univariate analysis of the correlation between clinicopathological characteristics and OS in BC was also performed. A forest plot was developed to visualize these univariate and multivariate analyses. Based on overexpressed PIMREG and other potential predictors, a predictive model based on TCGA-BRCA data was also developed to predict the mortality risk (Balachandran et al., 2015; Iasonos et al., 2008; Liu et al., 2018). A nomogram using ‘rms’ and ‘survival’ R packages was constructed to create a graphical representation of potential predicting factors for calculating mortality risk for an individual patient and predicting the 1-, 3-, and 5-year overall survival. All figures were visualized using the R package ‘ggplot2’.

GSEA

Based on transcriptional sequences obtained from TCGA, GSEA was used to identify gene sets and pathways associated with PIMREG. Gene expressions were classified into high-expression and low-expression groups. To assess the potential functions of these groups, we conducted a comparative analysis using GSEA via the Broad Institute website (https://www.gsea-msigdb.org/gsea/index.jsp), utilizing the R package cluster profiler (Subramanian et al., 2005; Yu et al., 2012). The results of this analysis were visually presented using the R package ‘ggplot2’.

Explore other genes/proteins that interacted with PIMREG

The protein-protein interaction (PPI) networks of the PIMREG were constructed via GeneMANIA (http://www.genemania.org) (Vlasblom et al., 2015). In addition, the MaxLink web server (https://funcoup5.scilifelab.se/maxlink/) was conducted to find and rank other cancer gene candidates by their connectivity to the PIMREG and the two most commonly expressed genes for BC, i.e., BRCA1 and BRCA2 (Guala, Sjölund & Sonnhammer, 2014).

Infiltration of immune cells

We utilized TIMER2.0, a web-based tool for comprehensive molecular characterization of tumor-immune interactions, to examine the infiltration of different immune cells and their clinical effects (Li et al., 2017). Using the Gene module in TIMER2.0, we created plots to explore the correlation between the expression of PIMREG members and immune infiltration level in BC. A significant correlation is set for the cutoff value of Cor >0.2 and p<0.05 (Li et al., 2016; Liu et al., 2021; Zhong et al., 2021).

To investigate interactions and relationships between the immune system and PIMREG, we employed the TISIDB database (http://cis.hku.hk/TISIDB). TISIDB is an online platform for evaluating immune system and tumor interactions, including almost 1,000 reported immune-related anti-tumor genes, etc., and immunological data aggregated from seven public databases (Ding et al., 2021; He et al., 2021; Xu et al., 2021). This study used TISIDB to explore highly associated immunomodulators with PIMREG in BC.

Relationship between PIMREG and immune checkpoints in BC

To investigate the potential oncogenic role of PIMREG in breast cancer (BC), we assessed its associations with immune checkpoints using TIMER and GEPIA. In this study, we focused on the immune checkpoints programmed death 1 (PD1)/PD1-L1 and cytotoxic T-lymphocyte-associated protein 4 (CTLA-4) as they are well-recognized critical checkpoints in tumor immune evasion (De Velasco et al., 2017; Li et al., 2017; Lou et al., 2021; Xu et al., 2018).

Investigate potential drugs that demonstrate interaction with PIMREG in BC

Q-omics (version 1.0) was employed to investigate drugs that exhibit a potentially favorable response to PIMREG. The score, developed by Q-omics, represents the log-fold change in drug response between samples with high PIMREG expression and those with low PIMREG expression. This scoring system was utilized to identify drugs demonstrating the strongest response to high PIMREG expression (Jeong et al., 2020; Lee et al., 2021).

Human tissue samples collection

We collected five surgically removed breast cancer (BC) tissues and their matched adjacent normal tissues from Northern Jiangsu People’s Hospital, affiliated with Yangzhou University, between October 2022 and February 2023. The clinical information of the patients is provided in Table 1. All participants in this study provided written informed consent, and the study protocol was approved by the ethics committee of Northern Jiangsu People’s Hospital, affiliated with Yangzhou University (Approval ID: 2022KY316).

Cell lines

The human mammary epithelial cell line MCF-10A (#CL-0525) and the human breast cancer cell line MDA-MB-231 (#CL-0150B) were obtained from Procell (Wuhan, China). MCF-10A cells were cultured in DMEM/F12 (Gibco, CA, USA) supplemented with 5% horse serum, 0.5 μg/ml hydrocortisone, 20 ng/ml human epidermal growth factor, 10 μg/ml insulin, and 1% penicillin/streptomycin. MDA-MB-231 cells were cultured in Dulbecco’s modified Eagle’s medium (DMEM, Hyclone, UT, USA), supplemented with 10% fetal bovine serum and 1% penicillin/streptomycin. Before further experiments, these cells were incubated in a humidified atmosphere containing 5% CO₂ at 37 °C.

Quantitative real-time reverse transcription PCR (qRT-PCR)

To assess the mRNA levels in cells, quantitative real-time polymerase chain reaction (qRT-PCR) was performed following a previously published protocol (Xie et al., 2022). Total RNA was extracted from the cells using the TRizol reagent (Takara, Shiga, Japan). Subsequently, the PrimeScriptTM RT reagent kit (Takara, Shiga, Japan) was utilized to synthesize single-stranded cDNA from the RNA, following the manufacturer’s instructions. The 2 −ΔΔCt value method, which involves determining the cycle at which the fluorescence level reaches the threshold limit set by the qPCR system, was employed to calculate the relative gene expression in the target and reference samples (Rao et al., 2013). This study used it to calculate the relative mRNA expression of PIMREG or HDAC2. Glyceraldehyde 3-phosphate dehydrogenase (GAPDH) was used as an endogenous control. The sequences of RT-PCR primers were used as follows: PIMREG: forward 5′-CCTGGAAACGCCTGGAAAC-3′and reverse 5′-CAAAGCACTCTTAGCTGAGCG-3′; HDAC2: forward 5’- ATGGCGTACAGTCAAGGAGG-3’ and reverse 5’-TGCGGATTCTATGAGGCTTCA-3’; GAPDH: forward 5′-ATTCCACCCATGGCAAATTC-3′and reverse 5′-TGGGATTTCCATTGATGACAAG.

Western blot assay

To examine the protein expression of PIMREG in cells and BC tissues, a western blot assay was performed based on a previously published protocol (Xie et al., 2022). Briefly, the cells or BC tissues were lysed using RIPA reagent (Beyotime, Shanghai, China) and boiled for 10 min. Subsequently, the total proteins were separated by 10% sodium dodecyl sulfate-polyacrylamide gel electrophoresis (SDS-PAGE) and transferred onto PVDF membranes (Millipore, Billerica, MA, USA). The membranes were then blocked with 5% nonfat milk diluted with Tris-buffered saline-Tween (TBST) at room temperature for 1 h and incubated overnight at 4 °C with the primary antibody: anti-PIMREG (1:500; abs152589; Absin Bioscience, Inc., Shanghai, China). The next day, the membranes were incubated with the HRP-conjugated secondary antibodies at room temperature for 1 h. The Western blots were visualized using an ECL detection kit (Proteintech Group, Inc., Rosemont, IL, USA) for chemiluminescence. The density of the protein bands was quantified using Image J software (Pierce, Rockford, IL, USA), and the values was normalized to β-Actin.

Immunohistochemistry assay

The immunohistochemistry (IHC) assay followed a previously published protocol (Jin et al., 2018). Briefly, paraffin-embedded breast tissue slides with a thickness of 5 µm were incubated overnight at 4 °C with anti-PIMREG (1:50; abs152589; absin, Shanghai, China). Subsequently, the slides were incubated with the secondary antibody at room temperature for 1 h. Then, the slides were incubated in diaminobenzidine (DAB) and counterstained with hematoxylin. Images of the slides were captured under a microscope.

Transfection assay

To knockdown the expression of the PIMREG, the siRNA transfection assay was conducted based on the previously published paper (Yao et al., 2019). Briefly, cells were seeded into a 6-well plate and transfected with Lipo8000™ Transfection Reagent (#C0533, Beyotime, Shanghai, China) following the manufacturer’s instructions. After 48 h of siRNA transfection, the knockdown efficiency was confirmed by conducting a western blot assay or qRT-PCR. The siRNA targeting human PIMREG or HDAC2 mRNA was designed by Genepharma Technologies Inc (Shanghai, China), and the siRNA sequences were used as follows: siPIMREG: 5 ′-UCCUGGAAACGCCUGGAAATT-3 ′, si-HDAC2: 5 ′-GGAUUACAUCAUGCU AAGA-3 ′.

Colony formation assay

To investigate the impact of HDAC2 on breast cancer, a colony formation assay was performed following a previously published protocol (Kim & Ma, 2018). Briefly, MDA-MB-231 cells were transfected with si-HDAC2 and subsequently seeded in 6-well plates at a density of 1000 cells per well. The cells were then incubated for 7 days to allow colony formation. Afterward, the colonies were fixed with 4% paraformaldehyde for 30 min and stained with 0.5% crystal violet for 30 min. The number of colonies was quantified using ImageJ software.

Cell counting kit-8 (CCK-8) assay

To assess the cell proliferation of MDA-MB-231 cells transfected with siRNA-PIMREG, a CCK-8 assay was performed following a previously published protocol (Xie et al., 2022). Briefly, after transfection with siRNA-PIMREG, MDA-MB-231 cells were seeded in 96-well plates at a volume of 100 μL/well, with 5000 cells per well. The cells were then incubated for 0 h, 24 h, 48 h, and 72 h. Subsequently, 10 μL of CCK8 solution (Beyotime, Shanghai, China) was added to each well, and the plate was further incubated for 30 min. Finally, the cell proliferation was measured by recording the absorbance at 450 nm using a microplate reader.

Transwell migration assay

To assess the migration ability of cells transfected with siRNA-PIMREG, a transwell migration assay was conducted following a previously published protocol (Yao et al., 2019). Briefly, after transfection with siRNA-PIMREG, 5 × 10⁴ cells were seeded in the upper chamber with an 8 µm pore size (Corning Costar). The lower chamber was filled with DMEM containing 10% FBS to stimulate cell migration. After incubation for 12 h, the migrated cells were fixed with 4% paraformaldehyde for 30 min and stained with 0.1% crystal violet for 30 min. Images of the migrated cells were captured under a microscope.

Statistical analysis

The bioinformatics analysis in this study was conducted using the online database mentioned above or the R software. The experimental results are presented as mean ± standard deviation (SD), and the differences between the two groups were compared using the Student’s t-test. Each experiment was repeated three times. Statistical significance was defined as a P value <0.05 or a log-rank P value <0.05.

Correlation analysis of PIMREG expressions in BC

Among the data downloaded from TCGA, a total of 1109 samples were available. However, the demographic and clinical data were only accessible for 1083 patients with breast cancer (Table 2). The analysis revealed significant associations between the expression of PIMREG and various clinicopathological parameters. These parameters include race (P < 0.001), T stage (P < 0.001), pathologic stage (P = 0.002), histological type (P < 0.001), PR status (P < 0.001), ER status (P < 0.001), and PAM50 status (P < 0.001) in breast cancer patients. These findings suggest that PIMREG may serve as a promising novel biomarker in the context of breast cancer.

Upregulation of PIMREG in BC

The analysis of mRNA expression data from the TCGA database revealed a significantly higher average expression level of PIMREG in BC tissues compared to normal tissues (P < 0.001) (Fig. 1A). This finding was further confirmed by analyzing matched tumors and surrounding samples, which demonstrated higher expression of PIMREG in BC tumor tissues (Fig. 1B). Consistent results were obtained from the expression difference analysis based on the GEO database (P < 0.001) (Fig. 1C). Subsequent analysis showed a significant difference in PIMREG expression between pathologic stage I and stage II, with higher expression observed in stage II (P < 0.001) (Fig. 2A). Moreover, more advanced T stage was associated with higher expression of PIMREG (P < 0.001) (Fig. 2B). However, there was no association between the M and N stages with PIMREG expression (Figs. 2C–2D). Regarding hormone receptor status, PIMREG expression was found to be elevated in estrogen receptor (ER) negative BC tissues compared to ER-positive tissues, and it was even more increased in progestin receptor (PR) negative BC tissues compared to PR-positive tissues (P < 0.001) (Figs. 2E–2F). However, no statistically significant difference was observed between human epidermal growth factor receptor 2 (HER2) positive and HER2-negative BC tissues (Fig. 2G). Additionally, PIMREG expression was significantly higher in infiltrating ductal carcinoma compared to infiltrating lobular carcinoma (P < 0.001) (Fig. 2H). To further support these findings, immunohistochemistry (IHC) images obtained from datasets of the Human Protein Atlas (HPA) were used to demonstrate PIMREG protein expression in BC and breast tissues. The intensity of PIMREG staining was weak in normal breast tissue but moderate in BC tissues, confirming the elevated expression of PIMREG in BC compared to normal breast tissue (Fig. 3). Overall, these results indicate that PIMREG is likely overexpressed in breast cancer.

Prognostic property of PIMREG in BC

The TCGA database was used to assess the impact of the upregulation of PIMREG on BC prognosis. Its upregulation was revealed to be associated with worse OS (HR: 1.48; 95%CI [1.17–2.03]; P = 0.018) and DSS (HR: 1.78; 95%CI [1.15–2.76]; P = 0.009), as well as shorter PFI (HR: 1.64; 95%CI [1.18–2.28]; P = 0.003) in patients with BC (Figs. 4A–4C). In addition, receiver operating characteristic (ROC) analysis was also performed with an area under the curve (AUC) value of 0.940 (95%CI [0.919–0.961]), indicating a high prognostic value of PIMREG (Fig. 4D). Further analysis of subgroups found that the upregulation of PIMREG was associated with worse OS in infiltrating lobular carcinoma with HR of 2.79 (95%CI [1.13–6.75]; P = 0.026), but not in infiltrating ductal carcinoma. (Figs. 5A–5B) Based on the GEO database, the higher expression of PIMREG was also found to be associated with worse OS in BC (HR: 2.28; 95%CI [1.13–4.59]; P = 0.021). (Fig. 5C) As PAM50 classifies BC into five molecular intrinsic subtypes that differ in their biological properties and prognoses, a subgroup analysis based on PAM50 subtypes was also performed, but did not indicate statistical significance (Fig. 6).

To further explore and identify the risk factors in patients with BC, a univariate Cox proportional-hazards regression analysis was conducted. The results showed that T stage, M stage, N stage, pathologic stage, age, and menopause status were associated with worse OS (P = 0.046, P < 0.001, P < 0.001, P = 0.003, P < 0.001, and P = 0.003, respectively) (Fig. 7). In addition, the higher expression of PIMREG was detected (HR = 1.475, 95%CI [1.070–2.034], P = 0.018) and was determined to be a significant predictor of worse OS.

The above results suggested that the T stage, N stage, M stage, age, pathologic stage, and menopause status might be linked to prognosis. As the pathologic stage of breast cancer is based on the TNM stage system, thus only the TNM stage was applied here for nomogram development. Therefore, a prognostic nomogram based on TNM stage, age, menopause, and PIMREG expression level was developed to predict individual survival probability (Fig. 8A). The C index of the nomogram was 0.741(0.713−0.769), indicating its potential for the prediction of overall survival. The calibration curve of this prediction model showed that the established lines of 1-, 3- and 5-year survival matched the ideal line at a relatively high degree (Fig. 8B). Collectively, these results indicated that PIMREG may be a prognostic signature for BC.

Identities in PIMREG-related signaling pathway identified by GSEA

GSEA was conducted to identify the differences in signaling pathways between low- and high-PIMREG expression datasets (Fig. 9). The result revealed that PIMREG was related to the neuroactive ligand–receptor interaction, NABA-secreted factors, cell cycle checkpoint, Class A1 Rohdopsin-like receptor, DNA repair, mRNA splicing, and translation pathways.

Immune cell infiltration of PIMREG in BC

The relationship between PIMREG and immune cell infiltration was studied using the TIMER database, as immune cell levels often correspond to cancer cell proliferation and progression (Fig. 10) (Grenier, Yeung & Khanna, 2018). A statistically significant correlation was found between PIMREG expression and the infiltration of macrophages, neutrophils, and dendritic cells. However, the correlation coefficients did not quite reach conventional levels of a positive Spearman’s rho value (<0.3).

To further explore the association between PIMREG and immunomodulators, three immunoinhibitors that were most correlated with PIMREG were Lymphocyte Activating 3 (LAG3), Indoleamine 2,3-dioxygenase 1 (IDO1), and CTLA-4 (Figs. 11A–11C). On the other hand, the three immunostimulators that showed the strongest correlation with PIMREG were poliovirus receptor (PVR), human UL16-binding protein 1 (ULBP1), and tumor necrosis factor (TNF)-receptor superfamily 13C (TNFRSF13C) (Figs. 11D–11F).

Regarding major histocompatibility complex (MHC) molecules, transporter associated with antigen processing 2 (TAP2), transporter associated with antigen processing 1 (TAP1), and human leukocyte antigen-A (HLA-A) showed significant associations with PIMREG (Figs. 11G–11I). These findings suggest that PIMREG may be involved in immune cell infiltration and immune modulation in breast cancer.

Relationship between PIMREG and immune checkpoints in BC

The relationship between PIMREG and immune checkpoints, including PD1 (PDCD1), PD1-L1 (CD274), and CTLA-4, which play a crucial role in tumor immune evasion, was assessed (Zhang et al., 2021). According to the analysis of TIMER data, PIMREG expression showed a relatively strong positive correlation with CTLA-4 only in BRCA (invasive breast carcinoma). This correlation was further adjusted by tumor purity using partial Spearman’s correlation (Fig. 12A). However, the analysis using GEPIA data did not find a significant positive correlation between PIMREG and PD1 (PDCD1), PD1-L1 (CD274), and CTLA-4 (Figs. 12B–12D).These findings suggest that the relationship between PIMREG and immune checkpoints may vary depending on the specific analysis method and dataset used. Further studies are needed to explore the potential interactions between PIMREG and immune checkpoints in breast cancer.

Other genes/proteins interacted with PIMREG in patients with BC

GeneMANIA analysis revealed that PIMREG is implicated in various functions, including protein deacetylation, CHD-type complex, transcription repressor complex, histone deacetylase complex, ATPase complex, and SWI/SNF superfamily-type complex. The top 20 interactors associated with PIMREG included phosphatidylinositol binding clathrin assembly protein (PICALM), Tudor domain-containing 7 (TDRD7), RB binding protein 7 (RBBP7), RB binding protein 4 (RBBP4), metastasis-associated 1 family member 2 (MTA2), histone deacetylase 2 (HDAC2), suppressor APC domain-containing 2 (SAPCD2), kinesin family member 15 (KIF15), chromatin licensing and DNA replication factor 1 (CDT1), histone deacetylase 1 (HDAC1), prion protein (PRNP), PHD finger protein 19 (PHF19), thymidylate synthetase (TYMS), zyg-11 family member B (ZYG11B), cyclin E1 (CCNE1), BUB1 mitotic checkpoint serine/threonine kinase (BUB1), ubiquitin-associated and SH3 domain-containing B (UBASH3B), SPC24 component of NDC80 kinetochore complex (SPC24), BLM RecQ like helicase (BLM), and transforming acidic coiled-coil-containing protein 3 (TACC3) (Fig. 13A).

Network analysis using MaxLink of BRCA1, BRCA2, and PIMREG revealed that MTA2, Cell Division Cycle 20 (CDC20), HDAC2, and Cyclin-Dependent Kinase 1 (CDK1) were important links connecting PIMREG with BRCA1/2 (Fig. 13B). Among these links, HDAC2 was the only direct connection between PIMREG and both BRCA1 and BRCA2, and MTA2 and HDAC2 were identified as interactors via GeneMANIA.

Investigate potential drugs that demonstrate interaction with PIMREG in BC

Using Q-omics (version 1.0), we analyzed 493 potential drugs. Among them, we identified 46 drugs that positively responded to increased expression of PIMREG (Fig. 14). The selection of these drugs was based on the Y-score, which represents the log (fold-change) of the drug response between samples with high PIMREG expression and samples with low PIMREG expression. The top three responsive drugs, based on the Y-score, were Mitoxantrone, Bleomycin, and AZD7762. For more detailed information about these potential drugs, please refer to Supplementary File 1.

Knockdown of PIMREG decreased the cell proliferation and migration in MDA-MB-231 cells

As anticipated, both qRT-PCR and western blot assays revealed the overexpression of PIMREG in BC cells when compared to normal human mammary epithelial cells (Figs. 15A–15B, P < 0.01). The western blot assay conducted on tissue samples (Fig. 15C, P < 0.01) and the experimental IHC assay (Fig. 15D) further confirmed the overexpression of PIMREG in BC tissue samples. To assess the impact of PIMREG on cell proliferation and migration, siRNA was employed to knockdown PIMREG expression (Figs. 15E–15F, P < 0.01). Following transfection with PIMREG-siRNA, both the CCK8 assay and transwell migration assay demonstrated a significant reduction in cell proliferation and migration in MDA-MB-231 cells when compared to the untreated group (Figs. 15G–15H, P < 0.01). These findings suggest that PIMREG may play a crucial role in BC cell proliferation and migration. Additionally, the colony formation assay was conducted to verify the effects of HDAC2 (Fig. 15I). The results demonstrated that the colony formation ability of MDA-MB-231 cells was reduced after the knockdown of HDAC2, indicating that HDAC2 may serve as a potential therapeutic target for breast cancer.