Discovery of plasma messenger RNA as novel biomarker for gastric cancer identified through bioinformatics analysis and clinical validation

Genome-wide expression analysis by Oncomine

Expression profiling, including 304 GC cancer samples and 174 normal controls, was obtained from Oncomine microarray database (http://www.oncomine.org) (Rhodes et al., 2004). In order to analyze the expression pattern of cancer vs. normal tissue mRNA, we focused on primary tumors and the following cut-offs were employed p-value ≤ 10⁻⁴, fold change ≥ 2 and gene rank ≤ 10%. Heat maps of overexpressed and underexpressed mRNAs in GC were available for each study.

Clinical specimens

Peripheral blood samples from 56 patients with gastric adenocarcinoma were addressed before therapeutic intervention by venipuncture and processed within 2 h at the First Affiliated Hospital of Xiamen University. We also collected blood samples from 14 healthy volunteers. The healthy volunteers presented neither a history of cancer nor other diseases. All patients were pathologically diagnosed as having gastric cancer using surgical specimens and biopsies. Plasma was isolated from 4 ml blood specimens after centrifugation at 1,600× g for 10 min at 4 °C and 10,000× g for 10 min, and then stored at −80 °C until the next step. Demographic, clinical and histopathological parameters of all these cases were summarized in Table 1. All experimental protocols were approved by the Clinical Research Ethics Committee of the First Affiliated Hospital of Xiamen University. All methods were performed in accordance with the Declaration of Helsinki. Written informed consent was obtained from all human participants after complete description of the study (Fig. S1 and Table S2).

RNA extraction

RNA was extracted from 500 µl plasma using TRIzol LS reagent (cat#10296018, Thermo Fisher Scientific Inc.) following manufacturer’s instructions as previously described (Pucciarelli et al., 2012; Zhang et al., 2012). In brief, 500 µl plasma was mixed with 500 µl TRIzol reagent. After 5 min incubation at 4 °C, 500 µl chloroform was added to the mixture, and violently shaken for 30 s. The mixture was immediately centrifuged at 12,000× g for 5 min at 4 °C. The above aqueous layer was transferred into a fresh tube containing 800 µl isopropyl alcohol. Next, the mixture was centrifuged at 12,000× g for 5 min at 4 °C and washed with 1 ml 70% ethanol for twice. Lastly, the RNA pellet was dissolved in 15 µl RNase-free water. Qubit RNA HS Assay Kit and Qubit 3.0 fluorometer (ThermoFisher Scientific Inc. cat#Q32851) were employed to quantify the concentration of RNA solution through spectrophotometry. The concentration of RNA isolated from plasma ranged was 35.24–168.18 ng/µl.

Quantitative PCR (qPCR)

The extracted total RNA was reverse-transcription into cDNA using PrimeScript RT reagent Kit (TAKARA cat#RR047A) according to the manufacturer’s protocol in triplicate. The resulting cDNA was stored at −80 °C for next PCR amplification. Primer sequences were designed through web-based version 4.1.0 of Primer 3 and were shown in Table S1. For clinical validation of the bioinformatics analysis results, qPCR was conducted by ABI ViiA 7 Real-Time PCR System (Applied Biosystems) with melting curve analysis. QPCR was carried out in triplicate at 50 °C for 2 min, denaturing at 95 °C for 5 min, followed by 40 cycles at 95 °C for 30 s, 58 °C for 1 min. Glyceraldehyde-3-phosphate dehydrogenase (GAPDH) was selected as an internal control. Water negative controls contained all components for the qPCR reaction without target RNA. Positive controls of RNA were extracted from SGC7901 cells obtained from the Cancer Center of Xiamen University (Xiamen, China). The 2^−ΔCT algorithm (ΔCT = Ct. target − Ct. reference) was employed for data analysis (Maru et al., 2009).

The Kaplan–Meier plotter analysis

The prognostic value of candidate circulating mRNAs was analyzed using the Kaplan–Meier plotter database, an online database containing 54,675 genes on survival based on 1065GC samples with a mean follow-up of 33 months (Szasz et al., 2016). Overall survival of patients with high and low expression levels of target genes was displayed using Kaplan–Meier survival curves. Hazard ratios (HR) with 95% confidence intervals (CI), and log-rank p-values were also calculated and summarized.

Functional enrichment network

Gene functional network was performed using gene ontology enrichment (GO) analysis. Enrichment map were created using the Cytoscape (v3.6.0). FDR < 0.05 was considered to be significant.

Statistical analysis

The Mann–Whitney U test was used to compare the expression status of circulating mRNAs in normal and GC groups and calculate the relationship between clinicopathologic characteristics and expression levels of relevant mRNA. Data was shown as median and range. Receiver operating characteristic (ROC) curves and the area under the curve (AUC) was used to identify the diagnosis value of selected mRNA (Brumback, Pepe & Alonzo, 2006). Cut-off values were assessed at different sensitivities and specificities and at the maximum Youden’s index = (sensitivity + specificity − 1) (Youden, 1950). Then, the logistic regression model was performed to obtain a combined ROC curve. GraphPad Prism 6.0 (GraphPad Soft-ware Inc., La Jolla, CA) and SPSS (version 22.0, IBM SPSS co., USA) was used for these statistical analyses. A two-sided p value of less than 0.05 was defined as statistically significant.

Identification of candidate mRNAs from Oncomine database

We compared the mRNA expression level in gastric cancer vs. normal samples obtained from Oncomine database. In total, 304 GC samples including 50 diffuse gastric adenocarcinoma, 21 gastric adenocarcinoma, 113 gastric intestinal type adenocarcinoma, 22 gastric mixed adenocarcinoma, 6 gastrointestinal stromal tumor and 92 other GC and 174 normal controls from 5 selected microarray datasets were analyzed (Chen et al., 2003; Cho et al., 2011; Cui et al., 2011; D’Errico et al., 2009; Wang et al., 2012). As shown in Fig. 1, expression profiling analysis generated 40 altered mRNA expressions in GC. The top 10 overexpressed and underexpressed genes embodied important genes involved in carcinogenesis and progression as well as several uncharacterized candidates (collagen type VI alpha 3 chain (COL6A3), serpin family H member 1 (SERPINH1), pleckstrin homology and RhoGEF domain containing G1 (PLEKHG1), collagen type I alpha 2 chain (COL1A2), claudin 1 (CLDN1), metallopeptidase inhibitor 1 (TIMP1), NOP56 ribonucleoprotein (NOP56), collagen type IV alpha 1 chain (COL4A1), immunoglobulin like domain containing receptor 1 (ILDR1); cadherin 11 (CDH11), pepsinogen 4, group I pepsinogen A (PGA4), potassium voltage-gated channel subfamily E regulatory subunit 2 (KCNE2), gastric intrinsic factor (GIF), solute carrier family 9 member A4 (SLC9A4), estrogen related receptor gamma (ESRRG), ATPase H+/K+ transporting subunit alpha (ATP4A), ATPase H+/K+ transporting subunit beta (ATP4B), PRDM16 divergent transcript (FLJ42875), MAS related GPR family member D (MRGPRD) and tripartite motif containing 50 (TRIM50), respectively).

Experimental validation of the potential noninvasive biomarkers

To validate our expression profiling analyses, six selected genes including top three overexpressed (COL6A3, SERPINH1 and PLEKHG1) and underexpressed (PGA4, KCNE2 and GIF) mRNAs were experimental verified using circulating mRNA extracted from 56 GC patients and 14 healthy subjects.

The expression of COL6A3, PLEKHG1, PGA4 and KCNE2 was detected in 55/56, 41/56, 32/56 and 55/56 GC plasma samples respectively, whereas the expression of SERPINH1 and GIF was detected in all GC samples. As shown in Fig. 2, the expression level of COL6A3 (P = 0.0128), SERPINH1 (P = 0.0217) and PLEKHG1 (P = 0.0064) was significantly higher in GC plasmas than those in normal controls. However, the expression levels of PGA4, KCNE2 and GIF had no significant change.

Diagnostic value of plasma COL6A3, SERPINH1 and PLEKHG1 for GC

We performed ROC curve analysis to assess the diagnostic value of these three circulating mRNAs (Fig. 3). The plasma level of COL6A3 had a sensitivity of 100% and a specificity of 46.2% to differentiate the GC patients from the healthy subjects, with an AUC of 0.720 at the optimal cut-off point (95% CI [0.522–0.919]). SERPINH1 had a sensitivity of 58.9% and a specificity of 78.6%, with an AUC of 0.698 (95% CI [0.543–0.852]). PLEKHG1 had a sensitivity of 68.3% and a specificity of 100%, with an AUC of 0.833 (95% CI [0.699–0.968]). Combining any two of the biomarkers had values for ROC area ranging from 0.676 to 0.870, sensitivities from 60% to 82.9%, specificities from 76.9% to 100% (Fig. 4). Furthermore, we used logistic regression analysis to combine these three circulating mRNAs and obtained an increased AUC value of 0.907 (95% CI [0.820–0.993]), with a sensitivity of 82.9% and a specificity of 100%.

Associations between clinicopathological characteristics and plasma COL6A3, SERPINH1 and PLEKHG1

The associations between clinicopathological characteristics and these three circulating mRNAs were investigated. As shown in Table 2, the higher COL6A3 level was significantly associated with increased lymph node metastasis (P = 0.0233), whereas the elevated expression of SERPINH1 was associated with advanced age (P = 0.0034) and poor differentiation (P = 0.0231) in GC patients. No significant association between PLEKHG1 and any clinicopathological characteristics including age, sex, stage, lymph node metastasis or differentiation was found.

Increased COL6A3, SERPINH1 and PLEKHG1 is associated with poor prognosis

To further evaluate whether the expression levels of COL6A3, SERPINH1 and PLEKHG1 can predict prognosis, we performed a survival analysis based on publicly gene expression datasets from the Kaplan–Meier Plotter resource. As shown in Fig. 5, the higher expression of COL6A3 (HR = 1.32, 95% CI [1.11–1.58], p = 0.0018), SERPINH1 (HR = 1.97, 95% CI [1.61–2.41], p = 3.1e⁻¹¹) and PLEKHG1 (HR = 1.34, 95% CI [1.07–1.69], p = 0.012) were all significantly correlated with poor OS in GC. These results indicated that GC patients with high COL6A3, SERPINHI or PLEKHG1 tend to have unfavorable outcome.

GO functional enrichment analysis

Functional enrichment network of COL6A3, SERPINH1 and PLEKHG1 was constructed. As shown in Fig. 6A, COL6A3 was predicted to have the main functions: extracellular matrix organization, cell adhesion, multicellular organism development, animal organ development and system development. SERPINH1 played roles in extracellular matrix organization, collagen fibril organization, skeletal system development, collagen metabolic process, and animal organ morphogenesis (Fig. 6B). PLEKHG1 was found associated with Rho guanyl-nucleotide exchange factor activity (Fig. 6C).