Serum microRNAs as potential new biomarkers for cisplatin resistance in gastric cancer patients

Cell culture

The MGC803 cell line was obtained from the Chinese Academy of Medical Sciences (Beijing, China). DDP-resistant GC cells were established by culturing MGC803 cells in a continuous stepwise fashion with gradually increasing concentrations of DDP from 0.05 to 5 µg/ml over a period of 15 months. This resulted in a shift towards cell proliferation and an apoptosis phenotype (Hong et al., 1988; Yu, Ma & Chang, 2000). The cell line was continuously cultured with 2 µg/ml DDP to maintain its resistance. These two cell lines were cultured following standard culturing procedures and identified by short tandem repeat (STR) profile data compared with the ATCC, DSMZ, or JCRB databases (Cobioer Biosciences Co., Ltd, Los Altos, CA, USA).

Patient enrollment and ethics statement

Serum specimens were collected from 74 GC patients who met the following criteria: (1) histologically confirmed gastric adenocarcinoma; (2) DDP-based chemotherapy used as first-line treatment; (3) received at least 2 cycles of neoadjuvant chemotherapy or palliative treatment; (4) chemotherapy efficacy could be evaluated by unenhanced and enhanced computed tomography (CT) after 2 or 3 cycles; (5) serum samples were collected before the first chemotherapy. All patients provided informed consent prior to the collection of a blood sample. This study was approved by the Ethics Committee (Approval Number: 2018-P2-045-01) of Beijing Friendship Hospital, Capital Medical University and met the ethical requirements of the Declaration of Helsinki. Prechemotherapy serum samples were collected between January 2015 and January 2019 from patients who received neoadjuvant chemotherapy or palliative treatment. The serum samples were collected and processed following the standard operating procedure of the Early Detection Research Network. The chemotherapy principle was executed according to the guidelines of the American Society of Clinical Oncology (ASCO). The chemotherapy response effect was evaluated according to CHOI’s principle (Choi et al., 2004). Treatment resistance was evaluated by the existing progressive disease (PD) in <3 treatment cycles, and the treatment response was evaluated by the existing complete response (CR), stable disease (SD), or partial response (PR) lasting for 2 or 3 cycles. The chemotherapy response-sensitive group included CR and PR patients, whereas the chemotherapy response-resistant group included SD and PD patients.

RNA extraction

Total RNA was isolated from the GC cells using Trizol (Invitrogen, USA). The same amount of Caenorhabditis elegans cel-39-3p miRNA was spiked into each serum sample as an external calibration for RNA extraction, reverse transcription, and miRNA amplification. Total RNA was extracted and purified from serum using the miRNeasy Serum/Plasma kit (Qiagen, cat. 217184). The quantity and integrity of the RNA yield was assessed using the Qubit2.0 fluorometer (Life Technologies, USA) and Agilent 2200 TapeStation (Agilent Technologies, USA), separately.

sRNA-seq and data analysis

Total RNA (1 µg) of GC cells were used to prepare small RNA libraries by NEBNext Multiplex Small RNA Library Prep Set for Illumina (NEB, USA) according to the manufacturer’s instructions. The libraries were sequenced by HiSeq 2500 (Illumina, USA) with single-end 50 bp reads at RiboBio Co. Ltd (RiboBio, China). Raw reads were processed by FastQC to get clean reads by filtering out those containing adapter, poly ‘N’, were of low quality, or were smaller than 17nt reads. Mapping reads were obtained by mapping clean reads to the reference genome from the BWA software. miRDeep2 was used to identify known mature miRNA based on miRBase21 (http://www.miRBase.org) and to predict novel miRNA. The databases of Rfam12.1 (http://www.rfam.xfam.org) and piRNABank (http://www.pirnabank.ibab.ac.in) were used to identify ribosomal RNA (rRNA), transfer RNA (tRNA), small nuclear RNA (snRNA), small nucleolar RNA (snoRNA), and piwi-interacting RNA (piRNA) by BLAST. miRNA expressions were calculated in RPM (reads per million) values (RPM = (number of reads mapping to miRNA/number of reads in clean data) × 10⁶). The expression levels were normalized by RPM [(number of reads mapping to miRNA/number of reads in clean data) × 10⁶].

Bioinformatics analysis

The differential expression between the two sample sets was calculated using the edgeR algorithm according to the criteria of —log2 (Fold Change) — ≥1 and P-value <0.05. TargetScan, miRDB, miRTarBase, and miRWalk were used to predict the target genes of selected miRNA. Patients with different expression levels of miRNAs were determined to have shorter or longer overall survival (OS) times for 10-year OS using the Kaplan–Meier survival curve data according to the Kaplan–Meier Plotter database (http://kmplot.com) (Szász et al., 2016; Nagy et al., 2018), the OncoLnc database (http://www.oncolnc.org/)  (Anaya, 2016), and the OncomiR database (http://www.oncomir.org/) (Wong et al., 2018).

RT-qPCR

The primers for miRNA detection were purchased from FulenGene Company (Guangzhou, China) and their sequence is listed in Table 1. Total RNA was reverse transcribed using the All-in-One™ miRNA RT-qPCR Kit (GeneCopoeia, USA) for miRNA analysis. The All-in-One™ miRNA RT-qPCR Detection Kit (GeneCopoeia, USA) was used to measure miRNAs quantitatively. All samples were normalized by the initial biofluid input volume used for RNA extraction and were calibrated by the spike-in cel-39-3p to eliminate the minute bias caused by different RNA isolation efficiencies and PCR efficiencies among samples. U6 or let-7g-5p was used as an endogenous control to normalize the relative number of miRNAs. ABI 7500 real-time fast PCR system (Applied Biosystems, USA) was used to achieve the relative quantitation of miRNA expression and the data was analyzed using the affiliated software. The Ct attenuation value of each type of miRNA in each of the serum samples was corrected by the internal reference cel-39-3p and the two housekeeping genes, U6 and let-7g-5p. The 2^−ΔΔCt method was used to calculate the relative expression of each miRNA. Each sample per reaction was performed in triplicate.

Statistical analysis

miRNA expression levels in the serum samples were divided into high and low groups taking the median as the cutoff value by RT-qPCR. Statistical differences for in vitro experiments were analyzed using student’s unpaired t tests. Associations between the clinical parameters of the patients and their miRNA expression were analyzed using the Mann–Whitney test. The AUC was used to assess the diagnostic accuracy of the predictors. Logistic regression was used to develop a panel of combined biomarkers to predict the probability of GC chemotherapy response. All data were expressed as mean ± SD. Each experiment was repeated independently at least three times. Quantitative data were analyzed and graphed using SPSS 23.0, MedCalc, and GraphPad Prism 7. Differences were considered to be significant at ****P < 0.0001, ***P < 0.001, **P < 0.01, and *P < 0.05.

Analysis of sRNA-seq data

The results of sRNA-seq in MGC803/DDP cells (Fig. 1A) and MGC803 cells (Fig. 1B) showed that the sRNA sequence included miRNA, tRNA, rRNA, snRNA, snoRNA, piRNA, and other RNA. Thirty-five miRNAs were selected and intersected, of which 11 miRNAs were up-regulated, and 24 down-regulated at higher levels in MGC803/DDP cells than the MGC803 cells (Fig. 1C). These 35 differentially expressed miRNAs in sRNA-seq data are shown in Table 2.

Bioinformatics analysis of selected miRNA

A large cohort analysis was conducted using Kaplan–Meier survival data according to the Kaplan–Meier Plotter database in order to determine the role of the 35 differentially expressed miRNAs as a potential prognostic factor. We observed that a total of seven differentially expressed miRNAs showed statistical differences (p < 0.05) for 10-year overall survival (OS) (Table 2), which indicated that these miRNAs could be used as valid biomarkers for chemotherapy response and overall survival. The differentially expressed miRNAs observed were miR-9 (including miR-9-3p and miR-9-5p, Fig. 2A), miR-146a (Fig. 2B), miR-370 (Fig. 2C), miR-433 (Fig. 2D), miR-519a (Fig. 2E), and miR-522 (Fig. 2F). TargetScan, miRDB, miRTarBase, and miRWalk were used to predict the target genes of the seven selected miRNAs, including miR-9-3p (Fig. 3A), miR-9-5p (Fig. 3B), miR-146a-5p (Fig. 3C), miR-370-3p (Fig. 3D), miR-433-3p (Fig. 3E), miR-519a-5p (Fig. 3F), and miR-522-5p (Fig. 3G). Supplementary material from four different databases (Table S1) was added to the common target genes of the 7 miRNAs (Fig. 3). These 7 miRNAs were selected as miRNA candidates.

Validation of sRNA-seq data by RT-qPCR analysis

GC cell results showed that the levels of miR-146a-5p, miR-519a-5p, and miR-522-5p were significantly downregulated in MGC803/DDP cells. The levels of miR-9-3p, miR-9-5p, miR-370-3p, and miR-433-3p were significantly upregulated in MGC803/DDP cells (Table 3, Fig. 4A). miR-519a-5p and miR-522-5p changed consistently with sRNA-seq, but the opposite results were seen from the other five miRNAs. miR-519a-5p and miR-522-5p had the same primer and exhibited the same result.

Expression levels and functions of miRNAs in human GC clinical specimens

RT-qPCR tests were run on the serum samples from 74 GC patients to identify the potential miRNAs. 2 housekeeping genes (U6 and let-7g-5p) and a spike-in cel-39-3p served as the internal references.

The clinical parameter findings of 74 patients with neoadjuvant chemotherapy or palliative treatment are shown in Table 4. Using the median ratio as the cutoff for the relative miRNA expression (7 miRNA candidates; Fold Change = 0) in serum, patients were classified into two groups: low miRNA (7 miRNA candidates) and high miRNA (7 miRNA candidates). Clinical data was collected and efficacy was evaluated by CT in 2–3 cycles following chemotherapy. Independent measurements were made by two radiologists according to CHOI’s principle (Choi et al., 2004). Unenhanced-combined-enhanced CT was applied to accurately evaluate the tumor response using tumor density measurement. RT-qPCR analysis was used to detect miRNA (7 miRNA candidates) levels in 34 chemotherapy response-sensitive (CR and PR) and 40 chemotherapy response-resistant (SD and PD) patients to validate miRNA (7 miRNA candidates) expression levels.

Our results proved that miR-9-3p expression obviously increased in the serum samples of chemotherapy response-resistant GC patients; 67.5% (27 of 40) of the high miR-9-3p samples showed DDP resistance (p < 0.0001, Fig. 4B). miR-9-5p expression increased in the serum samples of chemotherapy response-resistant GC patients; 72.5% (29 of 40) of the high miR-9-5p samples showed DDP resistance (p < 0.0001, Fig. 4C). miR-146a-5p expression was obviously induced in the serum samples of chemotherapy response-resistant GC patients; 62.5% (25 of 40) of the low miR-146a-5p samples showed DDP resistance (p < 0.0001, Fig. 4D). miR-433-3p expression increased in the serum samples of chemotherapy response-resistant GC patients; 67.5% (27 of 40) of the high miR-433-3p samples showed DDP resistance (p < 0.0001, Fig. 4E). However, there was no significant difference between the relative levels of miR-370-3p (p =0.0611, Fig. 4F), miR-519a-5p ( p = 0.4028, Fig. 4G), and miR-522-5p (p = 0.4028, Fig. 4H) in the serum of 34 chemotherapy response-sensitive gastric cancer samples and 40 chemotherapy response-resistant gastric cancer samples.

Our results indicate that miR-9-3p, miR-9-5p, miR-146a-5p, and miR-433-3p may act as potential new biomarkers for the chemotherapy response of DDP treatment. The levels of miR-9-3p, miR-9-5p, miR-146a-5p, and miR-433-3p were not related to gender and age (Table 4).

Receiver operating characteristic (ROC) analysis of miRNAs

Four miRNAs distinguished the GC chemotherapy response-resistant (SD+PD) group from the GC chemotherapy response-sensitive group (CR+PR). The results for the area under the curves (AUC), standard deviation (SD), 95% confidence intervals (CI), P-values, sensitivities (SE), and specificities (SP) of these miRNAs were as follows (Table 5): miR-9-3p (AUC = 0.824, Fig. 5A), miR-9-5p (AUC = 0.856, Fig. 5B), miR-146a-5p (AUC = 0.799, Fig. 5C) and miR-433-3p (AUC = 0.838, Fig. 5D), respectively. These four miRNAs provided promising AUC values for differentiating the GC chemotherapy response-resistant groups from the GC chemotherapy response-sensitive groups.

Two of the four candidate miRNAs were combined in a logistic model, with a significantly improved performance compared with the individual miRNA, determined by (Table 5): miR-9-3p + miR-9-5p (AUC = 0.889, Fig. 6A), and the risk score factors (RSF) were calculated as 0.731 × miR-9-3p + 0.586 × miR-9-5p + 0.192; miR-9-3p + miR-146a-5p (AUC = 0.903, Fig. 6B), and RSF = 0.800 × miR-9-3p-0.772 × miR-146a-5p + 0.392; miR-9-3p + miR-433-3p (AUC = 0.865, Fig. 6C), and RSF = 0.574 × miR-9-3p + 0.684 × miR-433-3p + 0.173; miR-9-5p + miR-146a-5p (AUC = 0.901, Fig. 6D), and RSF = 0.575 × miR-9-5p − 0.746 × miR-146a-5p + 0.345; miR-9-5p + miR-433-3p (AUC = 0.898, Fig. 6E), and RSF = 0.551 × miR-9-5p + 0.847 × miR-433-3p + 0.264; miR-146a-5p + miR-433-3p (AUC = 0.885, Fig. 6F), and RSF = −0.637 × miR-146a-5p + 0.802 × miR-433-3p + 0.327.

Three or more of the four candidate miRNAs were combined resulting in a significantly improved performance compared with the individual miRNA (Table 5: miR-9-3p + miR-9-5p + miR-146a-5p) (AUC = 0.930, Fig. 7A), and the RSF was calculated as 0.690 × miR-9-3p + 0.515 × miR-9-5p − 0.738 × miR-146a-5p + 0.386; miR-9-3p + miR-9-5p + miR-433-3p (AUC = 0.915, Fig. 7B), RSF = 0.514 × miR-9-3p + 0.526 × miR-9-5p + 0.614 × miR-433-3p + 0.240; miR-9-3p + miR-146a-5p + miR-433-3p (AUC = 0.918, Fig. 7C), RSF = 0.646 × miR-9-3p − 0.700 × miR-146a-5p + 0.504 × miR-433-3p + 0.388; miR-9-5p + miR-146a-5p + miR-433-5p (AUC = 0.926, Fig. 7D), RSF = 0.521 × miR-9-5p − 0.650 × miR-146a-5p + 0.655 × miR-433-3p + 0.376; miR-9-3p + miR-9-5p + miR-146a-5p + miR-433-3p (AUC = 0.937, Fig. 7E), RSF = 0.540 × miR-9-3p + 0.480 × miR-9-5p − 0.681 × miR-146a-5p + 0.393 × miR-433-3p + 0.398.

Pairwise comparison of ROC curves

MedCalc software was used to display the data for the different ROC curves and the results of pairwise comparison of all ROC curves (difference between the areas (DBA), standard error (SE), 95% confidence interval for the difference (CI), and P-value) (Table 6).

The ROC curve of miR-9-5p had the best AUC value of the four miRNA candidates. A comparison between the combined group (miR-9-5p + miR-9-3p + miR-146a-5p), and the miR-9-5p group showed that the two compared areas were significantly different (p = 0.026, Fig. 8A). The compared areas of the combined group (miR-9-5p + miR-9-3p + miR-433-3p) and the miR-9-5p group were significantly different (p = 0.045, Fig. 8B). The compared areas of the combined group (miR-9-5p + miR-146a-5p + miR-433-3p) and the miR-9-5p group were significantly different (p = 0.034, Fig. 8C) and the compared areas of the combined group (miR-9-5p + miR-9-3p + miR-146a-5p + miR-433-3p) and the miR-9-5p group were also significantly different (p = 0.016, Fig. 8D). However, the compared areas of the other combined group and the miR-9-5p group were not significantly different (P > 0.05).

Validate selected miRNAs in public database

To determine the role of miR-9-3p, miR-9-5p, miR-146a-5p, miR-433-3p, or their combination as a potential prognostic factor, a large cohort analysis was conducted using Kaplan–Meier survival data according to the Kaplan–Meier Plotter database, the OncoLnc database, and the OncomiR database (Table 7).

We observed that patients with higher miR-9-3p levels had shorter OS times than those with lower miR-9-3p levels in the OncoLnc database (p = 0.00836, Fig. 9A) and the OncomiR database (p = 0.06643, Fig. 9B). Patients with higher miR-9-5p levels had shorter OS times than those with lower miR-9-5p levels in the OncoLnc database (p = 0.00181, Fig. 9C) and the OncomiR database (p = 0.0044, Fig. 9D). Patients with higher miR-433-3p levels had shorter OS times than those with lower miR-433-3p levels in the OncoLnc database ( p = 0.0311, Fig. 9E) and the OncomiR database (p = 0.05651, Fig. 9F). These observations indicate that the high expression of miR-9-3p, miR-9-5p, or miR-433-3p are a valid biomarker for chemoresistance and poor survival, especially miR-9-5p.

Patients with higher levels of the combined group (miR-9-5p + miR-9-3p) had shorter OS times than those with lower combined levels in the OncomiR database (p = 0.003421, Fig. 10A). Patients with higher levels of the combined group (miR-9-5p + miR-433-3p) had shorter OS times than those with lower combined levels in the OncomiR database (p = 0.01428, Fig. 10B). Patients with higher levels of the combined group (miR-9-5p + miR-9-3p + miR-433-3p) had shorter OS times than those with lower combined levels in the OncomiR database (p = 0.008066, Fig. 10C). Patients with higher levels of the combined group (miR-9 + miR-433) had shorter OS times than those with lower combined levels in the Kaplan–Meier Plotter database (p = 0.0047, Fig. 10D). These observations indicate that higher levels of the combined groups, especially the combined group (miR-9-5p + miR-9-3p + miR-433-3p) are a valid biomarker for chemoresistance and poor survival.