Cell viability assay

GSP-treated PANC-1 cells were plated in 96-well cell culture plates at 5×10³ cells/well and incubated for 24 h at 37°C. Subsequently, 50 µl of MTT solution (5 mg/ml; Sigma-Aldrich; Merck KGaA) was added to each well and the cells were incubated for a further 4 h at 37°C. Following 3 min centrifugation (5500 × g) at 4°C, the supernatant was removed from each well. The formazan crystals produced from MTT in each well were dissolved in 150 µl DMSO and the optical density values were measured at 490 nm.

Flow cytometry

GSP-induced apoptosis in PC cells was determined by flow cytometry using the Annexin V-fluorescein isothiocyanate (FITC) Apoptosis Detection kit (BD Biosciences, Franklin Lakes, NJ, USA). Following treatment with GSPs for 48 h, 2×10⁵ cells were harvested, washed twice with PBS and incubated with Annexin V-FITC and propidium iodide for 10 min in the dark at room temperature. The stained cells were then detected and analyzed by the MoFlo XDP flow cytometer (Beckman Coulter, Inc., Brea, CA, USA) and Cell Quest 3.3 software (BD Biosciences, Franklin Lakes, NJ, USA).

RNA extraction and small RNA sequencing

Total RNA was extracted from PC cells using TRIzol reagent (Invitrogen; Thermo Fisher Scientific, Inc.), according to the manufacturer’s protocol. Subsequently, RQ1 DNase (Promega Corporation, Madison, WI, USA) was used to remove DNA. The quality and quantity of the purified RNA was monitored by absorbance at 260 and 280 nm, and the A260: A280 ratio, using a SmartSpec Plus Spectrophotometer (Bio-Rad Laboratories, Inc., Hercules, CA, USA). RNA integrity was further verified by electrophoresis of a 1.5% agarose gel.

Total RNA (3 µg) from each sample was used for small RNA cDNA library preparation with a Balancer NGS Library Preparation kit (Gnomegen, San Diego, CA, USA), according to the manufacturer’s protocol. The whole library was subjected to 10% native polyacrylamide gel electrophoresis and the bands corresponding to miRNA insertion were cut and eluted. Following ethanol precipitation and washing, the purified libraries were quantified using the QubitFluorometer (Invitrogen; Thermo Fisher Scientific, Inc.). These small RNA libraries were applied to NextSeq 500 (Illumina, Inc., San Diego, CA, USA) 76 nt pair-end sequencing.

Identification of conserved and novel miRNA

To obtain reliable clean reads, raw reads were processed by the FASTX-Toolkit (version 0.0.13)(http://hannonlab.cshl.edu/fastx_toolkit/). During this procedure, adaptor sequences and low quality tags were discarded. Based on the length of the mature miRNA and adapter lengths, RNAs <16 or >30 nts in length were also excluded from further analysis. The high-quality clean reads were subsequently searched against the Rfam database (version 12.0) (http://rfam.xfam.org/) using Bowtie (19). Matches to ribosomal RNAs and transfer RNAs were excluded. Subsequently, the remaining unique sequences were aligned against the miRBase database (20) using Bowtie, with one mismatch allowed. The matched small RNA sequences were considered to be conserved miRNAs and the unaligned sequences were potential candidates for novel miRNAs. To identify novel miRNAs, the unique sequences were aligned to the human reference genome sequence (GRCH38) using the miRDeep algorithm (21).

Bioinformatics analysis

To investigate the expression profiles of identified miRNAs, the frequency of miRNA counts were normalized to transcripts per million (TPM) using the following formula: Normalized expression=actual read count/total read count ×10⁶. DE miRNAs were analyzed using the edgeR (v3.22) package of Bioconductor software (22). A fold change (FC) ≤2 or <0.5 and P≤0.01 indicated a statistically significant DE miRNA. The miRanda algorithm (23) was used on human miRNA and transcript sequences of miRBase and TargetScan (7) for the prediction of putative miRNA targets.

To predict the gene function and calculate the frequency distribution of functional categories, GO and KEGG analysis were employed using the DAVID bioinformatics database (24). Networks were constructed by calculating the Pearson’s correlation coefficient (PCC) for the expression levels of DE miRNAs and target genes. Cytoscape (version 3.0.2) was used to display the co-expression network (25).

Validation of the target genes by reverse transcription-quantitative polymerase chain reaction (RT-qPCR)

To validate the RNA-seq data, RT-qPCR was performed for selected target genes and normalization was achieved with the human reference gene GAPDH. The primers used are presented in Table I. The same RNA samples for RNA-seq were used for RT-qPCR, and RNA extraction was performed with the same protocol and materials. In each pooled sample, l µg of total RNA was reverse transcribed using the PrimeScript™ RT Reagent kit (Takara Bio, Inc., Otsu, Japan), according to the manufacturer’s protocol. qPCR was performed with the S1000 Thermal Cycler (Bio-Rad Laboratories, Inc.) and Bestar SYBR Green RT-PCR Master mix (DBI Bioscience, Shanghai, China). The following thermocycling conditions were used: 95°C for 10 min, 38 cycles of 95°C for 15 sec and 60°C for 1 min. PCR amplifications were performed in triplicate for each sample, and the results were quantified by 2^−∆∆Cq methods (26).

Statistical analysis

All data are presented as the mean ± standard deviation. For comparisons between two groups, statistically significant differences between means were identified by Student’s t-test. For multiple comparisons, the significance was determined by one-way analysis of variance followed by a Tukey’s post hoc test. P<0.05 was considered to indicate a statistically significant difference.

miRNA expression profiles in PC cells treated with GSPs

To avoid high levels of necrotic cells and investigate how GSPs modulate the expression of miRNAs in PC cells, PANC-1 cells were treated with 20 µg/ml GSP for 3, 12 and 24 h, and the miRNA samples termed SS3, SS12 and SS24 were prepared, respectively. Cells treated with DMSO for 3, 12 and 24 h served as controls, and the control miRNA samples termed SC3, SC12 and SC24 were prepared, respectively. Two biological replicates were generated at each time point, therefore, a total of 12 small RNA libraries (SS3-A, SS3-B, SC3-A, SC3-B; SS12-A, SS12-B, SC12-A, SC12-B; SS24-A, SS24-B, SC24-A and SC24-B) were constructed for miRNA-seq.

Using Illumina NextSeq 500, >89.38 million reads were generated, corresponding to ~7.44 million sequence reads per sample and ~70% of all reads were successfully mapped against the current human reference genome (Table II). In addition, the proportion of clean reads mapped to the Rfam database for each sample is presented in Table III. In total, >79.6% of reads were matched to Rfam and multiple mapped reads were >84.2%.

The clean reads were mapped to miRBase and it was identified that ~39.3% reads were matched to mature miRNAs (Table IV). The read counts were normalized to TPM and a total of 2,578 miRNAs were revealed, consisting of 502 novel miRNAs (data not shown).

Determination of DE miRNAs

DE miRNAs were determined using edgeR (22) and the following criteria: P≤0.01 and FC ≥2 or ≤0.5. A total of 24, 83 and 83 DE miRNAs were identified in SS3 vs. SC3, SS12 vs. SC12 and SS24 vs. SC24, respectively (Fig. 1D and E). The highest and lowest FC values were 2×10^8.59and 2×10^−8.12, respectively (data not shown). These findings indicated that treatment with GSPs is associated with the expression levels of numerous miRNAs in PANC-1 cells. In addition, only three common DE miRNAs were identified between all three comparisons, suggesting that treatment with GSPs can regulate the expression levels of different miRNAs depending on the duration of treatment.

Prediction of DE miRNA target genes

miRNAs commonly exert their functions by binding to complementary target sites in the mRNAs of their target genes. Using the miRanda algorithm (23) on human miRNA and transcript sequences of miRBase and TargetScan (7), 74, 598 and 1,204 target genes were identified for the DE miRNAs in SS3 vs. SC3, SS12 vs. SC12 and SS24 vs. SC24, respectively (Fig. 2). This suggests that the DE miRNAs have the capacity to regulate the expression levels of multiple genes.

In our previous study, PANC-1 cells were also treated with 20 µg/ml GSPs for 3, 12 and 24 h, and 12 cDNA libraries (S3-A, S3-B, C3-A, C3-B; S12-A, S12-B, C12-A, C12-B; S24-A, S24-B, C24-A and C24-B) were constructed for RNA-seq. The current study analyzed the expression levels of the predicted target genes using the transcriptome database generated in our previous study (GSE85610). This revealed that all predicted target genes of the DE miRNAs were downregulated or upregulated following treatment with GSPs, suggesting that GSPs may modulate gene expression by the regulation of miRNAs. The expression levels determined by RT-qPCR of certain target genes, including MutS homolog 6, epidermal growth factor receptor, cyclin dependent kinase 6 and DNA methyltransferase 1, were in accordance with the RNA-seq results (Fig. 3).

Functional analysis of target genes

To identify the pathways associated with the target genes of the DE miRNAs, GO enrichment analysis was first conducted (Fig. 4). It was demonstrated that 7, 77 and 124 significant functional terms were identified for the target genes of the DE miRNAs in SS3 vs. SC3, SS12 vs. SC12 and SS24 vs. SC24, respectively (P<0.05; Fig. 4A, C and E). Multiple target genes were enriched in ‘positive regulation of cell proliferation’ (GO: 0008284) and ‘negative regulation of apoptotic process’ (GO: 0043066), which indicates that the DE miRNAs induced by treatment with GSPs may serve a regulatory role in the proliferation and apoptosis of PANC-1 cells. Targeting factors associated with cell cycle is a potential approach for cancer therapy; therefore, GSPs may exhibit inhibitory effects on the proliferation of cancer cells.

In addition, KEGG enrichment analysis was performed. It was revealed that 15, 42 and 47 significant functional terms were identified for the target genes of the DE miRNAs in SS3 vs. SC3, SS12 vs. SC12 and SS24 vs. SC24, respectively (P<0.05). Although ‘microRNAs in cancer’ (ID: hsa05206) and ‘pancreatic cancer’ (ID: hsa05212) were not presented in Fig 4 (only top10 terms were presented), numerous target genes were enriched in them and ‘pathways in cancer’ (ID: hsa05200). This suggests that treatment with GSPs may exhibit anticancer effects by disrupting the functional pathways of PC cells.