NcRNAs-mediated P2RX1 expression correlates with clinical outcomes and immune infiltration in patients with breast invasive carcinoma

P2RX1 is down-regulated in a variety of tumors

The mRNA expression of P2RX1 was analyzed to identify the differential expression in normal and tumor tissues. We first used the Cancer Genome Atlas (TCGA) database and found the mRNA levels of P2RX1 are comparatively lower in tumor tissues than paired adjacent normal tissues in bladder urothelial carcinoma (BLCA), BC, cervical squamous cell carcinoma and endocervical adenocarcinoma (CESC), colon adenocarcinoma (COAD), esophageal carcinoma (ESCA), kidney chromophobe (KICH), lung adenocarcinoma (LUAD), lung squamous cell carcinoma (LUSC), prostate adenocarcinoma (PRAD), rectum adenocarcinoma (READ), stomach adenocarcinoma (STAD), and uterine corpus endometrial carcinoma (UCEC), while higher in head and neck squamous cell carcinoma (HNSC) and kidney renal clear cell carcinoma (KIRC) (Figure 1A). Then we analyzed the P2RX1 expression between normal and tumor tissues using the TCGA database and the Genotype-Tissue Expression (GTEx) database, and similar results showed that P2RX1 is lowly expressed in a variety of tumors (Figure 1B). The Gene Expression Profiling Interactive Analysis (GEPIA) website further confirmed the differential expression (Figure 1C). This indicates that P2RX1 may be associated with an underlying carcinoma.

Figure 1. Expression analysis of P2RX1 in pan-cancer. (A) P2RX1 expression in 33 types of tumor tissues and paired adjacent normal tissues. (B) Differential expression analysis of P2RX1 based on TCGA and GTEx database. (C) P2RX1 expression in GEPIA database.

The para-cancerous tissues are complex tissues composed of multiple cell types, including normal mammary epithelial cells, immune cells, and fibroblasts. To evaluate whether P2RX1 is highly expressed in mammary epithelial cells, we examined the matched breast tumor/normal pairs by RT-qPCR and Western blot. RT-qPCR demonstrated that the expression of P2RX1 in normal human mammary epithelial cells (MCF 10A) was significantly downregulated, compared with that in human breast cancer cells (MDA-MB-231, MCF-7) (Figure 2A). The Western blot result showed the same trend (Figure 2B).

Figure 2. P2RX1 expression in normal human mammary epithelial cell line and two human breast cancer cell lines. (A) Analysis of P2RX1 expression by RT-qPCR. (B) Analysis of P2RX1 expression by Western blot.

P2RX1 is significantly associated with prognosis in BC

We evaluated the impact of P2RX1 on cancer prognosis among these 12 human cancers with differential expression of P2RX1 mRNA. Analysis of the GEPIA website manifested that high P2RX1 expression was correlated with better overall survival (OS) in BC, CESC and LUAD patients (p = 0.026, 0.035, and 0.0037, respectively) (Figure 3A). Besides, P2RX1 expression was also strongly associated with disease-free survival (DFS) in BC patients (p = 0.032) (Figure 3B). According to the survival analysis results, BC patients with low P2RX1 expression had 1.45 times the risk of OS and 1.52 times of risk of DFS than those with high P2RX1 expression. Further setting the cut-off to the 25th and 75th percentile expression levels of P2RX1 still showed that the high expression of P2RX1 was associated with better prognoses in BC patients (p = 0.025) (Figure 3C). Therefore, a low expression level of P2RX1 may be an unfavorable prognostic biomarker for BC patients.

Figure 3. Survival analysis. (A) Overall survival analysis of P2RX1 expression in BC, CESC, and LUAD using the GIPIA platform. (B) Disease-free survival analysis of P2RX1 expression in BC using the GIPIA platform. (C) Overall survival analysis of P2RX1 expression in BC with the cut-off of the 25th and 75th percentile expression levels.

Prediction of upstream miRNA of P2RX1

To explore the potential mechanism of the P2RX1 regulation, we firstly predicted possible upstream miRNAs and 11 upstream miRNAs that might be involved in the regulation of P2RX1 expression were discovered. Then we performed correlation analysis to determine the correlation between P2RX1 and these 11 miRNAs expression levels, and the results demonstrated that the expression of hsa-let-7c-5p (r = 0.223, p = 1.08E-13), hsa-let-7f-5p (r = 0.199, p = 3.88E-11), and hsa-let-7i-5p (r = 0.197, p = 5.96E-11) were positively correlated with the expression of P2RX1, while the expression of hsa-let-7b-5p (r = -0.169, p = 2.13E-08), hsa-let-7d-5p (r = -0.15, p = 7.07E-07), and hsa-let-7e-5p (r = -0.133, p = 1.03E-05) were inversely correlated with the expression of P2RX1 (Table 1). Previous studies have proved that miRNA usually negatively regulate the expression of targeted genes expression [16]. Therefore, the expression of hsa-let-7b-5p, hsa-let-7d-5p, and hsa-let-7e-5p in BC was assessed continuously. The results exposed that hsa-let-7b-5p was significantly downregulated in BC tissues, while hsa-let-7d-5p and hsa-let-7e-5p were significantly up-regulated which were positively correlated with the BC patients prognosis (p < 0.001) (Figure 4). Hsa-let-7d-5p and hsa-let-7e-5p were the most likely candidate miRNAs that might regulate P2RX1 expression.

Figure 4. Hsa-let-7d-5p and hsa-let-7e-5p are potential upstream binding miRNAs of P2RX1 in BC. (A) Differential expression analysis of hsa-let-7b-5p, hsa-let-7d-5p and hsa-let-7e-5p in BC based on TCGA and GTEx database. (B) The prognostic value of hsa-let-7d-5p and hsa-let-7e-5p in BC evaluated by Kaplan-Meier plotter.

Prediction of upstream lncRNA of hsa-let-7d-5p and hsa-let-7e-5p

The Starbase database was used to predict the upstream lncRNA of hsa-let-7d-5p and hsa-let-7e-5p. Hsa-let-7d-5p and hsa-let-7e-5p belong to the same miRNA family with similar upstream regulatory lncRNAs. Among the 40 predicted lncRNAs (Figure 5), 19 lncRNAs were lowly expressed in BC tissues (Supplementary Figure 1). We further analyzed the prognostic values of these 19 lncRNAs in BC and revealed that low STAG3L5P-PVRIG2P-PILRB expression was associated with a significantly worse OS in BC patients (Figure 6). However, there was no statistical difference in DFS. STAG3L5P-PVRIG2P-PILRB might be the upstream lncRNA of hsa-let-7d-5p and hsa-let-7e-5p in BC.

Figure 5. The lncRNA-miRNA regulatory network. The 40 lncRNAs predicted by the Starbase database were used to construct the lncRNA-miRNA regulatory network.

Figure 6. Overall survival and disease-free survival analysis of STAG3L5P-PVRIG2P-PILRB expression in BC using the GIPIA platform. High STAG3L5P-PVRIG2P-PILRB expression predicted better OS in BC patients (p = 0.014). But the two groups showed no statistical difference in DFS.

The construction and analysis of protein-protein interactions (PPI) network

The PPI network complex constructed by the STRING database consisted of 11 nodes and 44 edges (Figure 7). Important nodes included P2RX6, P2RX7, and P2RX5. The average local clustering coefficient was 0.893. P2RX1 was involved in calcium-mediated signaling using extracellular calcium sources and G protein-coupled ADP and ATP receptor activity. KEGG analysis indicated that P2RX1 mainly contributed to neuroactive ligand-receptor interaction, Platelet activation, and Calcium signaling pathway.

Figure 7. The protein-protein interaction (PPI) network of P2RX1 based on the STRING database. PPI network analysis showed the interaction of hub genes, including P2RX6, P2RX7, and P2RX5.

The correlation between P2RX1 and 6 types of immunocytes infiltration in BC

P2RX1 is highly expressed on the surface of immune cells, which may play an important role in immune cell infiltration. The significant differences in immunocytes infiltration levels under various copy numbers of P2RX1 were observed in Figure 8, which were mainly correlated with arm-level gain. Correlation analysis results indicated that the expression of P2RX1 positively correlates with 6 immune cells, including B cell, CD8+ T cell, CD4+ T cell, macrophage, neutrophil, and dendritic cell in BC. Further investigation of immune biomarkers shown in Supplementary Figure 2 validated that P2RX1 expression was significantly positively correlated with B cell’ biomarkers (CD27, CD19), CD8+ T cell’ biomarkers (CD8A and CD8B), CD4+ T cell’ biomarker (CD4), M1 macrophage’ biomarkers (CD68, CXCL10), M2 macrophage’ biomarkers (CCL18, CD163), neutrophil’ biomarkers (ITGAM), and dendritic cell’ biomarkers (CD83, CD1c). These findings support the positive correlation between P2RX1 and immunocytes infiltration in BC.

Figure 8. The relationship between P2RX1 and immune cell infiltration. (A) Infiltration levels of 6 immune cells with different copy number variations of P2RX1 in BC. (B) P2RX1 expression is positively correlated with 6 immune cells types evaluated by the TIMER database.

The correlation between P2RX1 and immune checkpoints in BC

At present, the main immune checkpoints are PD-1 (PDCD1), PD-L1 (CD274), and CTLA4 which are important biomarkers for predicting the efficacy of immunotherapy. As shown in Figure 9, P2RX1 expression was significantly positively correlated with PDCD1 (r = 0.507, p = 3.84e-66), CTLA4 (r = 0.455, p = 6.07e-52) and CD274 (r = 0.351, p = 3.26e-30) in BC after adjusted by tumor purity. Data from GIPIA showed results similar to those from the TIMER database which suggested P2RX1 may participate in the process of tumor immunomodulation.

Figure 9. The expression of P2RX1 is significantly positively correlated with PDCD1, CTLA4, and CD274 in BC evaluated by TIMER and GIPIA database.

Data acquisition and pan-cancer analysis

We downloaded transcriptome and clinical data of pan-cancer from the TCGA database (https://tcga-data.nci.nih.gov/tcga/) and GTEx database (http://www.gtexportal.org/). After sample normalization by Toil [30], differential expression analysis was performed on pan-cancers. Patients were classified into P2RX1 high and low expression groups based on the 50th percentile.

Differential expression and prognosis analysis

GEPIA is a website based on TCGA and GTEx data, using the output of a standard processing pipeline for high-throughput RNA sequencing data (http://gepia.cancer-pku.cn/) [31]. We used GEPIA to evaluate the expression difference of P2RX1 between tumor (TCGA samples) and normal (GTEx samples) tissues in 33 tumor tissues, as well as prognostic profiles including OS and DFS. GEPIA was also used to analyze the prognostic value of the candidate lncRNAs in BC patients. In addition, the correlation between P2RX1 and immune checkpoint inhibitors (PDCD1, CD274, and CTLA4) was determined using GEPIA.

Real-time quantitative reverse transcription PCR

Normal human mammary epithelial cells (MCF 10A) and human breast cancer cells (MDA-MB-231 and MCF-7) were bought from Procell Life Science and Technology Co., Ltd (Wuhan, China) and were authenticated with STR assays and tested negative for mycoplasma infection. Total RNA was extracted from MCF 10A, MDA-MB-231, and MCF-7 cells with an RNA extraction kit (Baidai, China), and RNA concentration was determined with Nanodrop (Thermo Fisher, USA). Reverse transcription was performed using ReverTra Ace qPCR RT Master Mix (TOYOBO, Japan). RT-qPCR was performed on a QuantStudio 6 Pro Real-Time PCR Systems (Thermo Fisher, USA) using the SYBR Green PCR kit (TOYOBO, Japan). The normalized internal reference was GAPDH mRNA. The relative expression of P2RX1 was calculated using the 2-ΔΔCT method. The primers sequence was as follows: forward, ATGGTGCTGGTGCGTAATAAG, and reverse, GGAAGACGTAGTCAGCCACA.

Western blot analysis

Normal human mammary epithelial cell line MCF 10A was cultured in DMEM containing 5% horse serum, 20ng/ml Epidermal Growth Factor, 0.5μg/ml Hydrocortisone, 10μg/ml Insulin, and 1% Non-Essential Amino Acids (CM-0525, Procell, China). The human breast cancer cell line MDA-MB-231 was cultured in Leibovitz's L-15 medium (PM151010, Procell, China) containing 10% fetal bovine serum (Gibco, USA). The human breast cancer cell line MCF-7 was cultured in DMEM medium (Gibco, USA) containing 10% fetal bovine serum (Gibco, USA). MCF 10A and MCF-7 were maintained in a humidified incubator at 37° C and 5% CO2, and MDA-MB-231 was maintained in a humidified incubator at 37° C.

Proteins from normal human mammary epithelial cells and human breast cancer cells were collected using RIPA buffer (Beyotime, China) and quantified using the BCA protein assay kit (Beyotime, China). Proteins were separated using 10% sodium dodecyl sulfate-polyacrylamide gel electrophoresis and transferred to polyvinylidene fluoride membranes. Membranes were blocked with 5% nonfat dry milk and washed with Tris-buffered saline containing 0.1% Tween 20.

The following primary antibody was used for overnight incubation at 4° C: rabbit anti-P2RX1 (1:2000, DF9728, Affinity, China) and mouse anti-GAPDH (1:4000, T0004, Affinity, China). Secondary antibodies were goat anti-rabbit IgG (H+L) HRP conjugate and goat anti-mouse IgG (H+L) HRP (horseradish peroxidase) conjugate (anti-rabbit, 1:4000, Affinity, China, and anti-mouse, 1:20000, Beijing Zhongshan). ECL hypersensitive luminescent liquid (Millipore, USA) and Multiplex Fluorescence Western blot imaging system (ProteinSimple, USA) were used to detect protein strips.

Associated miRNA prediction

In this study, we used 5 target gene prediction programs to predict potential upstream binding miRNAs of P2RX1, namely miRWalk, TargetScan, miRDB, RNA22, and miRmap. The miRNAs predicted by these five programs at the same time were considered candidate miRNAs and were used for further analysis.

Starbase database analysis

Starbase (http://starbase.sysu.edu.cn/), an online database for ncRNA functions and their coordinated regulatory networks, was used to predict candidate lncRNA, and analyze the correlation between miRNA and P2RX1, as well as candidate lncRNA and P2RX1 [32].

Kaplan-Meier plotter database analysis

Kaplan-Meier plotter database (http://kmplot.com/analysis/) was used to analyze the effect of potential miRNAs on the survival of BC patients who were divided into two groups based on the median expression of miRNAs. Kaplan-Meier survival plots were constructed and the log-rank p value was calculated.

STRING database analysis

The STRING tool (https://string-db.org) is a database of known and predicted PPI. STRING database was used to investigate the PPI network of P2RX1 and analyze the functional interactions. Nodes and edges of different colors indicate the relationship between proteins and proteins.

TIMER database analysis

TIMER database (https://cistrome.shinyapps.io/timer/) is a tool for the analysis of immune infiltrates. TIMER database was used to analyze the relationship between P2RX1 expression and tumor immune infiltrating cells including B cells, CD4+ T cells, CD8+ T cells, neutrophils, macrophages, and dendritic cells. Then the immune checkpoints were also evaluated to explore their correlation with P2RX1 expression.

Statistical analysis

The expression difference of P2RX1 in pan-cancer was compared by the Wilcoxon rank sum test. The survival differences were analyzed by log-rank test and Kaplan-Meier method. The significant difference was determined at p < 0.05. These analyses were conducted in R Version 3.6.1 and Stata/MP Version 13.1.