ERAP2 as a potential biomarker for predicting gemcitabine response in patients with pancreatic cancer

Drug sensitivity and transcripts in pancreatic cell lines

Firstly, we downloaded the sensitivity data of pancreatic cancer cell lines treated with gemcitabine from the public GDSC data, which included two data sets: GDSC1 and GDSC2. In this part, GDSC1 data set was adopted to screen genes related to gemcitabine sensitivity, which included IC50 values of 30 pancreatic cancer cells treated with gemcitabine (Table 1). Then, we conducted correlation analysis between gene expression value and IC50 value. We found that the IC50 of gemcitabine was positively relevant to 359 genes and negatively correlated with 456 genes (p < 0.05).

Univariate Cox survival analysis of gemcitabine sensitivity-related genes

We have identified a total of 815 gemcitabine sensitivity-related genes and annotated them with Gencod.v22 gene annotation file, but only 807 were successfully annotated. Then the TCGA-PAAD data set of pancreatic cancer patients was downloaded from the UCSC XENA database, including gene expression sequencing data (FPKM format), survival time, survival status, and gemcitabine medication information. 68 cases of patients with pancreatic cancer treated with gemcitabine were identified through data screening, and their corresponding gene expression sequencing data, survival time and survival status were obtained. After data cleaning and collating, 815 gemcitabine sensitivity relevant genes were selected for univariate Cox survival analysis using the Survival package of R software to determine the influence of gene expression on the patients. In order to establish the multi-factor Cox regression model, we selected the most significant 20 genes for variable screening (Table 2).

A multiple factors Cox proportional hazards regression model was established based on Lasso regression analysis

Next, we excluded OR2B2 for its barely expression, and finally 19 genes entered the variable selection of the model. Additionally, we utilized Lasso algorithm to screen out the independent variables. 1000 models were calculated by setting the parameter nlambda = 1000, and the coefficients of independent variables were obtained (Figure 1A). Besides, we filtered the lambda value (λ value) through 10-fold cross-validation, and obtained two models. It should be mentioned that one model was based on lambda.min (the error mean is the minimum corresponding lambda, dotted line on the left); the other was based on lambda.1se (the maximum lambda of the error mean within 1 standard deviation of the minimum value, dotted line on the right). As shown in Figure 1B, we adopted 9 genes as model variables for model construction. After determining the model variables, we employed the R package survival to construct the multi-factor Cox risk ratio model for these 9 genes, and identified the final model, which included the coefficients, risk values, and P values of each variable. Meanwhile, the coefficients of the best model of Lasso were also shown (Table 3). We adopted R package to plot the Forest plot of this model, and found that the P values of ERAP2 and DOCK11 were both less than 0.05, suggesting that they could be independent prognostic factors for pancreatic cancer patients treated with gemcitabine, respectively (Figure 1C). Moreover, we proved the Concordance Index of this model reaches 0.79, indicating that the model is accurate and reliable.

Figure 1. Lasso regression screening survival model variables. (A) Variable coefficient corresponding to the model with different number of variables. (B) 1000 models obtained through ten-fold cross validation. (C) Risk coefficient and P value of corresponding genes in multivariate Cox proportional risk regression model.

Risk factor association analysis

To further illustrate the validity of the Lasso-Cox regression model, we described the survival time and survival status of patients by a risk factor correlation diagram as their risk values changed. First of all, patients were assigned to high-risk group and low-risk group with cut off value equal to -0.3 (dotted line in the middle), which was calculated by Xtile software. We also mapped the distribution of survival time and status in pancreatic cancer patients based on the distribution of risk values. Patients in the high-risk group of the model group and the validation group showed shorter survival time and lower survival rates compared with the lower risk group (Figure 2A, 2B). We found ERAP2 and DOCK11 were important independent prognostic factors, and the relationship between their expression levels and risk values was also analyzed. As shown in the heat map, ERAP2 expression as a risk factor was higher in the high-risk group, while DOCK11 expression as a protective factor was lower in the low-risk group. The trend of model group and validation group is consistent, indicating a strong universality of the model (Figure 2C, 2D).

Figure 2. Correlation analysis of risk factors between model group and validation group. (A) The survival time and survival status of pancreatic cancer patients in the model group as a function of their risk values. (B) The survival time and survival status of pancreatic cancer patients in the validation group as a function of their risk values. (C) Heat map of gene expression as a function of risk values in the model group. (D) Heat map of gene expression as a function of risk values in the validation group.

Kaplan-Meier survival analysis based on risk score of pancreatic cancer patients

According to the variable coefficients of the optimal model of Lasso regression mentioned above (Table 3), the Lasso regression coefficients corresponding to 9 genes were obtained. Then we obtained the Risk Score of each pancreatic cancer patient according to the expression levels of 9 genes, and obtained the value of risk for 68 pancreatic cancer patients with gemcitabine treatment. After calculating the cut off value by Xtile software, patients with a risk value above -0.3 were considered to be in the high-risk group, while the rest were considered to be in the low-risk group. The corresponding Kaplan-Meier survival curve was shown in Figure 3A. At the same time, 113 pancreatic cancer patients without gemcitabine medication information were applied to validate the model. Similarly, the risk value of each patient was calculated by the expression levels of 9 genes, and the corresponding Kaplan-Meier survival curve was plotted (Figure 3B). It was obvious that the prognosis of low-risk patients in the model group (n = 68) and validation group (n = 113) was better than that of high-risk patients. We further draw the ROC curves of the prediction results of the model group and the validation group, and the 1-, 3- and 5-year ROC curves of pancreatic cancer patients, respectively. In the model group, the 1-, 3- and 5-year AUC area was 0.83, 0.83, and 0.76, respectively, indicating the good predictive ability on the survival of pancreatic cancer patients treated with gemcitabine of the model, and we observed that the Lasso-Cox regression model composed of these 9 genes had a very high accuracy (Figure 3C). In the validation group, the 1-, 3- and 5-year AUC areas were 0.6, 0.7 and 0.68, respectively. In the case of using external data, the model still maintained a high accuracy, indicating that the Lasso-Cox regression model constituted by these 9 genes had good universality (Figure 3D).

Figure 3. Evaluation and validation of Lasso-Cox regression model for pancreatic cancer. (A) Risk-based Kaplan-Meier survival curve for 68 pancreatic cancer patients treated with gemcitabine. (B) Risk-value Kaplan-Meier survival curve for 113 pancreatic cancer patients with no medication information. (C) AUC area of pancreatic cancer patients in the model group at 1, 3 and 5 years. (D) AUC area of pancreatic cancer patients in the validation group at 1, 3 and 5 years.

The expression of ERAP2 was associated with gemcitabine sensitivity and response

From the above Cox proportional risk regression model, we identified two genes, ERAP2 and DOCK11, associated with gemcitabine sensitivity. As a risk factor, ERAP2 expression is positively relevant to the risk value of pancreatic cancer patients, and is more suitable as a target for further research. First, we analyzed the relationship between ERAP2 expression and gemcitabine IC50 in pancreatic cancer cells in the GDSC1 and GDSC2 datasets. In GDSC1, the correlation between ERAP2 and gemcitabine log2 (IC50) value reached 0.5. In GDSC2, the correlation between ERAP2 and gemcitabine log2 (IC50) value also reached 0.47 (Figure 4A). In addition, we downloaded the expression chip GSE78982 of pancreatic cancer fibroblast from the GEO database, and analyzed the data by the GEO2R online tool. ERAP2 expression was significantly up-regulated in the gemcitabine resistant group, with differential multiple LogFc value of 1.19 and p-value of 0.0075 (Figure 4B). More importantly, we screened gemcitabine treated patients in the TCA-PAAD pancreatic cancer data set and measured their treatment efficacy, and found significantly lower ERAP2 expression values in the effective group (PR/CR) than in the ineffective group (PD) (Figure 4C). These results suggested that ERAP2 expression was significantly related to gemcitabine sensitivity in pancreatic cancer.

Figure 4. Relationship between expression of ERAP2 and sensitivity to gemcitabine in pancreatic cancer. (A) Correlation between the expression value of ERAP2 in GDSC1 and GDSC2 pancreatic cancer cell lines and the IC50 value of gemcitabine. (B) Expression of ERAP2 in two types of pancreatic cancer fibroblasts. (C) Relationship between the effect of gemcitabine treatment and ERAP2 expression in TCGA pancreatic cancer patients.

ERAP2 promoted pancreatic cancer progression

The above findings show that ERAP2 is a key gene affecting gemcitabine resistance in patients with pancreatic cancer, which enable us to further study whether there is a relationship between ERAP2 and the diagnosis and prognosis of pancreatic cancer. We downloaded the expression microarray GSE62452 from GEO database of pancreatic tumors and paracancerous tumors, and found that ERAP2 was highly expressed in pancreatic tumors (Figure 5A). In addition, ERAP2 was also significantly overexpressed in blood extracellular vesicle samples of pancreatic cancer patients from the BBCancer database (Figure 5B). Next, we analyzed the ERAP2 expression in pancreatic cancer by the GEPIA database and found that ERAP2 was significantly overexpressed in pancreatic cancer samples, which further demonstrated the cancer-promoting effect of ERAP2 (Figure 5C). In addition, we performed a Kaplan-Meier survival analysis in the TCGA pancreatic cancer dataset, using the median expression of ERAP2 as the threshold. The results are consistent with previous analysis, indicating that high ERAP2 expression is a risk factor for poor prognosis (Figure 5D), while DOCK11 failed to differentiate patient prognosis (Supplementary Figure 1). Finally, we adopted the immunohistochemical data of pancreatic tissues and pancreatic tumors in ProteinAtlas database to verify the expression of the protein level of ERAP2, and further confirmed the high expression of ERAP2 in pancreatic tumor (Figure 5E). The above results suggest that ERAP2 is overexpressed in pancreatic cancer, resulting in a poor prognosis.

Figure 5. Expression and prognosis of ERAP2 in pancreatic cancer. (A) Expression of ERAP2 in GSE45452 chip pancreatic tumor and paracancerous tumor. (B) The expression of ERAP2 in extracellular vesicles of blood origin. (C) Expression of ERAP2 in the pancreatic cancer dataset from GEPIA database. (D) Kaplan-Meier survival curve of ERAP2 in pancreatic cancer data set from GEPIA database. (E) Immunohistochemical images of ERAP2 in pancreatic tissues and pancreatic tumors from ProteinAtlas database.

Expression of ERAP2 was positively correlated with gemcitabine resistance

Next, we examined the viability of BxPC-3 and PANC-1 cells exposed to gemcitabine for 0, 24, 48 and 72 h. We discovered that gemcitabine reduced the proliferation of BxPC-3 and PANC-1 cells with IC50 of 10.49 μM and 17.7 μM (72 h) (Figure 6A). We then examined the ERAP2 levels in pancreatic cancer cell lines with gemcitabine treatment. As shown in Figure 6B, 6C, ERAP2 expression and transcription level were significantly increased after gemcitabine treatment. We conclude that ERAP2 expression may play an important role in the sensitivity of pancreatic cancer cells to gemcitabine.

Figure 6. The relationship between ERAP2 expression and gemcitabine resistance. (A) Cell viability assay in BxPC-3 and PANC-1cell lines treated with gemcitabine. (B) The Protein levels of ERAP2 in BxPC-3 and PANC-1 cell lines treated with gemcitabine were assessed by western blot. (C) The mRNA levels of ERAP2 in BxPC-3 and PANC-1cell lines treated with gemcitabine were quantified by qRT-PCR. Data were shown as the mean ± SD of three independent experiments (*p < 0.05, **p < 0.01 and ***p < 0.001).

Inhibition of ERAP2 attenuated the tumorigenesis and increased sensitivity to gemcitabine of pancreatic cancer cells

We next identified the tumorigenesis and sensitivity of pancreatic cancer cells to gemcitabine after ERAP2 knockdown. It was indicated that ERAP2 knockdown significantly restrained the growth of pancreatic cancer cells (Figure 7A, 7B). Interestingly, ERAP2 knockdown also significantly increased sensitivity to gemcitabine in these two pancreatic cancer cell lines (Figure 7A, 7B). We also tested the effects of ERAP2 on cell migration and invasion by wound healing assay and Transwell invasion assay. As suggested in Figure 7C, 7D, knockdown of ERAP2 remarkably attenuated the migration ability of pancreatic cancer cell lines, and the invasion ability of pancreatic cancer cell lines was also inhibited. Compared with gemcitabine alone, knockdown ERAP2 combined with gemcitabine further reduced the migration ability of pancreatic cancer cells, and the invasion ability was more significantly inhibited (Figure 7C, 7D). Overall, our results suggest that down-regulation of ERAP2 significantly blocks the tumorigenic ability of pancreatic cancer cells and significantly enhances the anti-pancreatic cancer activity of gemcitabine.

Figure 7. Effect of ERAP2 knockdown on tumorigenesis and gemcitabine sensitivity in pancreatic cancer cell lines. (A) Cell viability assay in two pancreatic cancer cell lines treated with ERAP2 knockdown, gemcitabine alone or gemcitabine with ERAP2 knockdown. (B) Cell proliferation assay in two pancreatic cancer cell lines treated with ERAP2 knockdown, gemcitabine alone or gemcitabine with ERAP2 knockdown. (C) Wound healing assay in two pancreatic cancer cell lines treated with ERAP2 knockdown, gemcitabine alone or gemcitabine with ERAP2 knockdown. (D) Transwell invasion assay in two pancreatic cancer cell lines treated with ERAP2 knockdown, gemcitabine alone or gemcitabine with ERAP2 knockdown. Data were shown as the mean ± SEM of three independent experiments (*p < 0.05, **p < 0.01 and ***p < 0.001).

Identification of genes related to gemcitabine sensitivity

Drug sensitivity data and transcript data were obtained from the GDSC database. The drug response data includes two datasets, GDSC1 and GDSC2. GSDC1 was a discovery set to obtain genes related to gemcitabine sensitivity, and the final results were validated in GDSC2. The IC50 values of pancreatic cancer cells and their corresponding transcripts were analyzed, and finally 359 genes positively as well as 456 genes negatively correlated with IC50 were obtained.

Construction of drug-sensitive prognosis model

To further screen out key genes, PAAD transcript and clinical data were collected from the TCGA database for subsequent analysis through the UCSC Xena browser (http://xena.ucsc.edu), and we totally obtained 68 pancreatic cancer patients who received gemcitabine. Finally, 815 candidate genes were analyzed by univariate Cox analysis (359 + 456), and the top 20 genes were taken into the Lasso algorithm to create the final multivariate Cox model. Additionally, we included other 113 pancreatic cancer patients as validation set. The R package glmnet (version 2.0.18) and ggrisk (version 1.2) were used for visualization of these results. The cut off value of risk score was obtained by Xtile software (version 3.6.1).

Functional annotation and pathway analysis

To investigate a comprehensive set of functionally annotated hub genes, David online analysis tool was utilized to carry out GO function annotation and KEGG pathway analysis on 359 positive correlation genes and 456 negative correlation genes, respectively (https://david.ncifcrf.gov). The regulatory network of ERAP2 is constructed using IPA software.

The analysis of expression and prognosis

The GEPIA2 online tool (http://gepia2.cancer-pku.cn) was used to analyze the expression and prognosis of the target gene, and the protein level of the gene was further verified in The Human protein atlas database (https://www.proteinatlas.org).

Cell culture

Pancreatic cancer cell lines PANC-1 and BxPC-3 were grown in Dulbecco's modified Eagle (DMEM, Gibco) with the addition of 10% fetal bovine serum (FBS; ExCell Bio).

Reagents

Gemcitabine was purchased from Selleck. Anti-GAPDH (1:5000) and ERAP2 (1:1000) antibodies were purchased from Proteintech.

shRNA transfection

6-well plates were inoculated with PANC-1 and BxPC-3 cells. Add 1 mL of lentiviral supernatant containing sh-ERAP2 (Genechem) lentiviral construct or control vector and 1 mL DMEM media with10% FBS to cells when the cell concentration reaches 50%-60%. Two days after infection, 2 μg/mL puromycin was utilized to screen cells that had been successfully infected. Two days later, the inhibition of ERAP2 was investigated by western blot.

Cell viability assay

The cells were evenly plated in a 96-well culture plate with 5000 cells per well. The cells were treated with 0, 0.1, 1 and 10 μM gemcitabine on the second day. The cytotoxicity of gemcitabine was determined at 0, 24, 48 and 72 h, respectively. 10 μL CCK-8 reagent (Biotool) was supplemented to each well, and the absorbance was determined at 450 nm after incubation at 37° C for 2 h.

Quantitative real-time PCR

Total RNA was extracted with MgZol reagent (Magen) according to the instructions. Quantitative real-time PCR was conduct to detect the mRNA level of ERAP2. Primers were designed and synthesized by Sangon Biotech. The primer sequences:

ERAP2:

Forward 5’- GAGGCGGAGTCTTGCTCTGTTG-3’

Reverse 5’- GAGGCAGGAGAATGGCGTGAAC-3’

Western blot

Cells were washed twice with PBS and collected. The mixture of RIPA buffer (Beyotime) and protease inhibitor (Bimake) was added. The cells were lysed for 30 min and the supernatant was collected by centrifugation. Protein concentration was determined by BCA protein assay kit (Beyotime). Protein was isolated with 10% SDS-PAGE and transferred to PVDF membrane after sufficient protein separation. The PVDF membrane was sealed with 5% skim milk at room temperature for 1 h, cleaned with PBS for 15min, and incubated overnight with primary antibody at 4° C. Subsequently, the PVDF membrane was incubated with secondary antibody at room temperature for 1 h. After cleaning with PBS for 30 min, ECL reagent was used for chemiluminescence detection.

Colony formation assay

Cells were seeded in 6-well culture plates with 1,000 cells per well. On the second day, cells were treated with 0, 0.1, 1, and 10 μM gemcitabine, respectively. After 24 h, the fresh medium was replaced. The fresh medium was replaced every 3 days thereafter. The cell status was observed after 10-14 days, and the medium was discarded after the cells formed suitable clones. After washed with PBS 2 times, the cells were fixed by 500 μL paraformaldehyde. After 15 min, paraformaldehyde was discarded. After washed with PBS 2 times, the cells were stained by 500 μL crystal violet. After 10 min, crystal violet was discarded. After washed with PBS 2 times, colonies were counted after the six-well plate becomes dry.

Wound healing assay

Cells were plated in 6-well plates until the density reached 80%-100%. The cells were scratched using a sterile 200μL spear tip. Wash the cells gently with PBS three times without breaking the drawn lines. Then 10 μM gemcitabine was added and cell microscope images were taken at 0, 24 and 48 h. Three different areas were selected from each well and the ability of cells to migrate was assessed according to the degree of healing.

Cell invasion assay

Cells invasion ability were evaluated using 24-well chemotaxis chambers (Costar, #3422). After pretreatment at 37° C for 12 h with a mixture of 60 ul serum-free medium (7): Matrigel (1) in the upper chamber, cells (50000 cells/well) suspended on 100 μL serum-free medium were added, followed by 500 μL culture base containing 30% FBS in the lower chamber. 24 h later, after washed with PBS 2 times, the cells were fixed by 500 μL paraformaldehyde. After 15 min, paraformaldehyde was discarded. After washed with PBS 2 times, the cells were stained by 500 μL crystal violet. After 10 min, crystal violet was discarded. After washed with PBS 2 times, photographed and counted.

Data analysis

Data were statistically analyzed by the unpaired two-tailed Student’s t-test of Graphpad Prism software (version 8.0), and P < 0.05 was considered as statistically significant difference between groups.