Dihydroartemisinin suppresses pancreatic cancer cells via a microRNA-mRNA regulatory network

Despite improvements in surgical procedures and chemotherapy, pancreatic cancer remains one of the most aggressive and fatal human malignancies, with a low 5-year survival rate of only 8%. Therefore, novel strategies for prevention and treatment are urgently needed. Here, we investigated the mechanisms underlying the anti-pancreatic cancer effects dihydroartemisinin (DHA). Microarray and systematic analysis showed that DHA suppressed proliferation, inhibited angiogenesis and promoted apoptosis in two different human pancreatic cancer cell lines, and that 5 DHA-regulated microRNAs and 11 of their target mRNAs were involved in these effects via 19 microRNA-mRNA interactions. Four of these microRNAs, 9 of the mRNAs and 17 of the interactions were experimentally verified. Furthermore, we found that the anti-pancreatic caner effects of DHA in vivo involved 4 microRNAs, 9 mRNAs and 17 microRNA-mRNA interactions. These results improve the understanding of the mechanisms by which DHA suppresses proliferation and angiogenesis and promotes apoptosis in pancreatic cancer cells and indicate that DHA, an effective antimalarial drug, might improve pancreatic cancer treatments.


INTRODUCTION
Despite decades of effort to improve treatments, pancreatic cancer remains one of the most aggressive and deadly human tumors, with a five-year survival rate of only 8% [1]. Pancreatic cancer is currently the eighth and ninth leading cause of worldwide cancer-related death in men and women, respectively [2]. However, patients who are treated with multimodal therapy, including surgical resection, have a 5-year survival rate of more than 20% [3]. This unfavorable prognosis for pancreatic cancer is primarily due to the poor therapeutic effects of chemotherapy and radiotherapy, the aggressive biological behavior of pancreatic tumors, and late diagnoses. Alternative therapeutic strategies are urgently needed to more effectively treat this deadly disease.
MicroRNAs, a family of non-coding RNAs approximately 21-23 nucleotides in length, regulate gene expression post-transcriptionally by targeting the 3'-untranslated regions (3'-UTR) of target genes through complementary base pairing, and play significant roles in a variety of biological processes, particularly in cancer cells [4,5]. MicroRNAs negatively regulate their targets by degrading mRNA or inhibiting its translation. microRNAs may target more than 30% of human genes [6] and regulate survival, apoptosis, proliferation, angiogenesis, invasion, metastasis, and chemotherapy resistance during the initiation and progression of human cancers [7][8][9]. Therefore, strategies that modulate microRNAs to regulate tumor amplification, apoptosis, metastasis, and invasion may provide promising novel treatments for pancreatic cancer.

Research Paper
G0/G1 cell cycle arrest [10] and apoptosis [17], inhibits growth [17] and angiogenesis [15], and increases sensitivity to chemotherapy [18,19] in pancreatic cancer. However, the details of the mechanisms underlying DHA's anti-pancreatic cancer effects have not been fully described; here, we further investigated these mechanisms.
Many experiments have confirmed that microRNAs play important roles as either promotors or inhibitors of cancer [20][21][22]. However, few experiments have demonstrated that drugs directly influence the ability of microRNAs to modulate their target genes and thus affect survival, apoptosis, proliferation, angiogenesis, invasion, metastasis, and chemotherapy resistance in cancer. In this study, we examined differentially expressed microRNAs from microarray profiles, confirmed microRNA-mRNA interactions using experimentally validated databases, and conducted a functional enrichment analysis of pathways from the Kyoto Encyclopedia of Genes and Genomes (KEGG) using quantitative real-time polymerase chain reactions and western blots, to investigate the mechanisms underlying DHA's anti-pancreatic cancer effects.

Aberrant microRNA expression profiles in PANC-1 and BxPC-3 pancreatic cancer cells
MicroRNA expression profiles in DHA-treated and control PANC-1 and BxPC-3 pancreatic cancer cells were examined using a microRNA microarray. microRNA expression results were normalized using median normalization and used to screen for differentially expressed microRNAs. Differentially expressed microRNAs were identified via fold-change filtering between the two groups in each cell line; a fold-change ≥ 1.5 identified microRNAs as up-or down-regulated. The identified microRNAs were then filtered between the two cell lines such that microRNAs for which expression was similarly up-or down-regulated in both PANC-1 and BxPC-3 cells were considered differentially expressed microRNAs. Seventy microRNAs that were up-regulated and 109 that were down-regulated (fold change ≥ 1.5) in both PANC-1 and BxPC-3 pancreatic cancer cells were identified using these filters ( Figure 1).

Figure 1:
The multi-step approach adopted to identify microRNAs differentially expressed after DHA treatment and their target mRNAs in two pancreatic cancer cell lines. www.impactjournals.com/oncotarget Identification of DHA-associated microRNA-mRNA interaction networks Four high-confidence microRNA-mRNA interaction databases (miR2Disease, miRecords, TarBase and miRTarBase) were used together to screen the experimentally identified differentially expressed microRNAs and their targets in DHA-treated cells. As shown in Figure 2A, 5,092 microRNA-mRNA interactions were identified between 62 of the 179 differentially expressed microRNAs (117 had no significant interactions) and 4,362 mRNAs. The node degree distribution of this network revealed a scale-free structure, with R 2 =0.86. These results indicate that the DHA-associated microRNA-mRNA interaction network is well-characterized by a core set of post-transcriptional regulatory relationships that are distinct from those of randomly generated networks (25800746).

Identification of DHA-associated microRNApathway regulation networks
To determine which molecular pathways are regulated by these differentially expressed microRNAs in response to DHA treatment, a functional enrichment analysis was performed based on the microRNA-mRNA interaction network. Neighbor mRNA nodes in the network for each microRNA were used as input gene sets to examine relationships with hundreds of KEGG pathways. The hypergeometric test was used to identify significantly enriched pathways (p<0.05). As shown in Figure 2C, 204 microRNA-pathway interactions were identified, involving 33 microRNAs and 59 pathways. In the microRNA-pathway interaction network, several important cancer-related pathway categories, such as "pathways in cancer," "cell cycle," "P53 signaling," and "pancreatic cancer pathways," were connected to more To identify microRNAs and their targets (mRNAs) that were regulated by DHA, four experimentally verified microRNA-mRNA interaction databases were examined. 5,092 microRNA-mRNA pairs were identified that involved 62 of the differentially expressed microRNAs. B. The selected microRNAs and their targets (mRNAs) were used to construct the DHA-regulated microRNA-mRNA regulatory network. C. A global view of the DHA-associated microRNA-pathway regulatory network. Functional enrichment analysis was performed based on the microRNA-mRNA interaction network to identify biological functions of microRNAs that were regulated by DHA. Ultimately, 204 microRNA-pathway interactions were constructed between 33 microRNAs and 59 pathways. D. The microRNAs were significantly associated with the KEGG pathway database pancreatic cancer pathway (p<0.05). microRNA nodes than others. These observations indicate that cancer risk-related pathways tend to be hub nodes in the microRNA-pathway interaction network. Cytoscape software (version 3.2.1) (14597658) was used to construct and illustrate the network.

Confirmation of relevant mRNA protein levels with western blots
To determine protein levels of the mRNAs modulated by DHA via microRNAs, we next examined Cdk4, Cdk6, VEGF, IKKα, MEK1, E2F3, Rac1, E2F1, and CDC42 protein levels, which affect growth, inhibit angiogenesis, and promote apoptosis in DHA-treated and vehicle-treated pancreatic cancer cells.

microRNAs identified by microarrays and systematic analysis contribute to the antipancreatic cancer effects of DHA
To confirm that the microRNAs identified by microarray and systematic analysis contribute to the antipancreatic cancer effects of DHA, western blots were used to measure target mRNA protein levels in the presence of microRNA inhibitors. As shown in Supplementary Figure S4, DHA reduced Cdk4, Cdk6, VEGF, IKKα, and Rac1 levels compared to the control group. However, transfection of AMO-195-5p, AMO-34a-5p, or AMO-30c-5p reversed this DHA-induced inhibition compared to the group treated with DHA alone. These results indicate that the anti-pancreatic cancer effects of DHA were mediated at least partly by miR-195-5p, mir-34a-5p, and mir-30c-5p.

DHA suppresses growth, inhibits angiogenesis, and promotes apoptosis in vivo
To examine the anti-tumor effects of DHA in vivo, we established a xenograft pancreatic cancer model in nude mice that did or did not receive 18 days of DHA treatment ( Figure 7A). Tumor sizes were measured every three days, and volumes were estimated according to the formula V=(π/6)×(larger diameter)×(smaller diameter) 2 ( Figure 7B). DHA treatment decreased tumor volumes over time compared to the control group. Mice were sacrificed and tumors excised after treatment ( Figure 7C).
To evaluate the anti-proliferation, anti-angiogenesis, and apoptotic effects of DHA in tumor tissues, Ki-67 (a cell proliferation marker) and CD31 (a microvessel density marker) levels and TUNEL staining were examined using immunohistochemistry. As shown in Figure 7D, tumors from DHA-treated mice had fewer Ki-67-positive cells than those from the vehicle-treated group. Similarly, DHA treatment reduced microvessel density compared to the vehicle-treated group ( Figure 7D). DHA treatment also increased the percentage of TUNEL positive cells compared to the vehicle-treated group ( Figure 7D).

DISCUSSION
The American Cancer Society estimates that approximately 53,070 Americans will be diagnosed with pancreatic cancer in 2016. Despite significant efforts to improve treatments, 41,780 Americans will die from pancreatic cancer, which has a five-year survival rate of only 8% [1]. Gemcitabine, a broad spectrum drug used to treat solid tumors, is well-tolerated in pancreatic cancer patients. However, patients who show initial sensitivity to gemcitabine chemotherapy rapidly acquire resistance, and the efficacy of gemcitabine in treating pancreatic cancer remains at only 20-30% [9,23]. Therefore, novel and effective chemotherapeutic agents for the treatment of pancreatic cancer are urgently needed.
DHA is an effective anti-malarial drug, and many studies have revealed that DHA also has anticancer effects in a variety of cancers. In a previous study, we found that DHA suppressed pancreatic cancer cell growth both in vitro and in vivo [17]. We also found that DHA downregulated cdks and cyclins, such as Cdk4, Cdk6, and cyclin E, which play essential roles in the regulation of cell cycle progression; DHA thus increased G0/G1 cell cycle arrest in pancreatic cancer cells. More importantly, DHA inhibited the DNA-binding activity of NF-κB in pancreatic cancer cells [10]. Furthermore, DHA increased the antipancreatic cancer effects of gemcitabine by inactivating NF-κB both in vitro and in vivo [18] and suppressed angiogenesis by regulating the NF-κB pathway [15]. Finally, we found that DHA enhanced Apo2L/TRAILmediated apoptosis via ROS-mediated up-regulation of death receptor 5 [19]. In summary, we have demonstrated that DHA, in addition to curing malaria, is a promising chemotherapeutic agent for treating pancreatic cancer. However, the mechanisms underlying the anti-pancreatic

Figure 7: DHA decreases proliferation and angiogenesis and increases apoptosis in pancreatic cancer cells in vivo.
A. Schematic representation of the experimental protocol described in the Methods. B. Tumors were measured using Vernier calipers and volumes were calculated using the formula V=(π/6)×(larger diameter)×(smaller diameter) 2 . C. Mice were sacrificed and tumors were excised at the end of treatment. D. Cell proliferation was measured by immunohistochemical analysis of Ki-67, apoptotic cells were measured by TUNEL staining, and microvessel density was measured by CD31 immunohistochemistry in tumor tissues. E. Ki-67-positive cells were counted to estimate the proliferation index. F. TUNEL-positive cells were counted to estimate the apoptosis index. G. CD31stained microvessels were counted to estimate microvessel density. *p<0.05, **p<0.01 compared to the control. cancer effects of DHA are not fully understood. Therefore, we investigated the mechanisms by which DHA exerts anti-pancreatic cancer effects using microarray profiles and systematic analyses.
To our knowledge, this is the first study to examine differences in microRNA expression after DHA treatment using microarrays and systematically analyze DHAassociated microRNA-mRNA interaction networks to identify the mechanisms by which DHA exerts its anticancer effects. Surprisingly, we found that four crucial microRNAs (miR-34a-5p, miR-195-5p, miR-30c-5p, and miR-130b-3p) regulated the expression of many mRNAs (Cdk4, Cdk6, VEGF, IKKα, MEK1, E2F3, Rac1, E2F1, ERK1, and CDC42) and their proteins, and thus were crucial to the anti-pancreatic cancer effects of DHA.
Cdk4, Cdk6, E2F3, and E2F1 play key roles in the regulation of cell cycle progression in pancreatic cancer. Down-regulation of Cdk4, Cdk6, E2F3, and E2F1 expression increases G0/G1 cell cycle arrest in pancreatic cancer cells [24,25]. Here, we found that DHA treatment up-regulated miR-34a-5p, miR-195-5p, miR-130b-3p, and miR-30c-5p expression and down-regulated the expression of the target mRNAs Cdk4, Cdk6, E2F3, and E2F1, respectively; DHA treatment also decreased protein levels translated from these mRNAs. VEGF plays a key role in angiogenesis [26], and down-regulation of VEGF expression suppresses angiogenesis in pancreatic cancer. Here, we found that DHA treatment down-regulated VEGF mRNA expression and protein levels by upregulating the expression of miR-34a-5p and miR-195-5p.
The Ras-Raf-MEK-ERK signalling pathway, which is one of best-characterized kinase cascades in cancer cell biology, influences various processes in tumors, including cancer cell survival, proliferation, migration, and differentiation [27]. MEK1 and ERK1 play key roles in the Ras-Raf-MEK-ERK signalling pathway [28]. In this study, DHA treatment up-regulated miR-34a-5p and miR-30c-5p and down-regulated MEK1 mRNA expression and protein levels, which both microRNAs target. These results indicate that the Ras-Raf-MEK-ERK signalling pathway is involved in the anti-pancreatic cancer effects of DHA. CDC42, a Rho family GTPase, also plays an B. Differential expression of the 9 mRNA targets in pancreatic tumor tissues. C. Cdk4, Cdk6, VEGF, IKKα, MEK1, E2F3, Rac1, E2F1, and CDC42 protein levels in pancreatic cancer tissues were detected by western blot. β-actin was used as a protein loading control D. The density of each band was measured and compared to β-actin. *p<0.05, **p<0.01 compared to the control. important role in various cancers [29,30] and promotes proliferation and metastasis [31,32]. Here, the DHAinduced increase in miR-195-5p down-regulated CDC42 expression, likely contributing to the ability of DHA to inhibit proliferation and suppress metastasis. Rac1, another Rho family GTPase, is also involved in various cancers, and down-regulation of Rac1 expression inhibits proliferation, viability, and migration in pancreatic cancer cells [33]. Here, DHA-induced, miR-30c-5p-mediated down-regulation of Rac1 expression might be another novel mechanism by which DHA inhibits cancer cell proliferation, viability, and migration.
In a previous study, we found that, although total NF-κB levels did not change, DHA treatment downregulated nuclear NF-κB levels in a dose-dependent manner [10]. Therefore, DHA prevented the transfer of NF-κB into the nucleus, and down-regulated gene products downstream of NF-κB, thus promoting apoptosis, suppressing growth, inhibiting angiogenesis, and increasing sensitivity to chemotherapy in pancreatic cancer cells. However, the mechanisms by which DHA prevents transport of NF-κB into the nucleus are unknown. Surprisingly, in the present study, DHA treatment upregulated miR-195-5p, which down-regulated IKKα mRNA expression. IKK phosphorylates and degrades IκB in a ubiquitination-dependent manner; this prevents the translocation of the RelA(P65) and NF-κB1(P50) dimer, which is sequestered in the cytoplasm in its inactive state via interaction with the endogenous inhibitor IκB [34,35]. This mechanism may explain the ability of DHA to inhibit the DNA-binding activity of NF-κB in pancreatic cancer cells.
In summary, we found that the DHA-induced microRNA-mRNA regulatory network suppressed growth, inhibited angiogenesis, and promoted apoptosis in pancreatic cancer. Our results provide mechanistic evidence that the anti-proliferative, pro-apoptotic, and anti-angiogenesis effects of DHA were associated with the up-regulation of miR-34a-5p, miR-195-5p, miR-30c-5, and miR-130b-3p. While some of these results confirmed our previous findings, most of them identified novel mechanisms underlying the anti-pancreatic cancer effects of DHA, which may be effective not only in treating malaria but also as a chemotherapeutic drug for treating pancreatic cancer.

Strategy
As shown in Figure 1, a multi-step approach was used to identify differentially expressed microRNAs that were up-regulated by DHA and, in turn, the mRNAs regulated by these miRNAs, in two pancreatic cancer cell lines. First, PANC-1 and BxPC-3 pancreatic cancer cells were treated with DHA or vehicle, and the microRNAs that were differentially expressed between the treatment groups were identified using microarrays. 179 differentially expressed microRNAs were identified using two standards: first, the microRNAs were up-or downregulated with a fold-change in expression ≥ 1.5; second, the microRNAs were differentially expressed similarly in both PANC-1 and BxPC-3 cell lines. Next, microRNA-mRNA interactions were identified using the miR2Disease, miRecords, TarBase, and miRTarBase databases; only microRNAs with confirmed target relationships in these databases were examined further. Based on the microRNA-mRNA interaction data, a DHA-associated microRNA-pathway regulatory network was identified, and the "pancreatic cancer pathway" specifically was selected for further investigation. Five microRNAs, 11 mRNAs, and 19 microRNA-mRNA interactions identified in the pancreatic cancer pathway were investigated to determine their roles in the mechanisms underlying the anti-pancreatic cancer effects of DHA. The differentially expressed microRNAs and mRNAs identified in this way were measured in vitro and in vivo using qRT-PCR, and levels of the protein products of these mRNAs were measured in vitro and in vivo using western blots. ThemicroRNA-target mRNA regulatory relationships selected by microarray and systematic analysis were confirmed with western blots in vitro. Finally, the antitumor effects of DHA were examined in vivo by measuring tumor sizes and via immunohistochemistry.

Cell culture
The human pancreatic cancer cell lines PANC-1, BxPC-3, SW1990, and CFPAC-1 were obtained from the American Type Culture Collection (Rockville, USA) and were cultured in DMEM or RPMI 1640 medium supplemented with fetal bovine serum (10%), penicillin (100 U/mL), and streptomycin (100 mg/mL) (Irvine Scientific, Irvine, CA). All cells were maintained at 37°C in humidified air with 5% CO 2 . All reagents were from HyClone China Ltd. (China). Mycoplasma contamination was tested using the Mycoplasma Stain Assay Kit (Beyotime Institute of Biotechnology, Beijing, China); none of the cell cultures were contaminated with mycoplasma. www.impactjournals.com/oncotarget

Treatment of cells
Subconfluent cells (60-70%) were treated with 50 μM DHA in DMSO in complete cell culture medium; control group cells were treated with 0.1% DMSO. Subsequent experiments were repeated three times.

RNA extraction
Total RNA was extracted and isolated using TRIzol reagent (Invitrogen Life Technologies) according to the manufacturer's instructions. RNA quantity and quality were measured using a NanoDrop ND-2000 spectrophotometer. Acceptable OD A260/A280 ratios were between 1.8 and 2.1 (the target is approximately 2.0 for pure RNA), and OD A260/A230 ratios were greater than 1.8. RNA integrity was assessed by standard denaturing agarose gel electrophoresis.

MicroRNA microarray procedure
After RNA was isolated from the samples, the miRCURY™ Hy3™/Hy5™ Power labeling kit (Exiqon, Vedbaek, Denmark) was used to label microRNAs according to the manufacturer's guidelines. One microgram of each sample was 3'-end-labeled with the Hy3 TM fluorescent label using T4 RNA ligase and the following procedure. RNA was suspended in 2.0 μL of water and combined with 1.0 μL of CIP buffer and CIP (Exiqon). The mixture was incubated for 30 min at 37°C, and the reaction was terminated by incubation for 5 min at 95°C. Then, 3.0 μL of labeling buffer, 1.5 μL of the fluorescent label (Hy3 TM ), 2.0 μL of DMSO, and 2.0 μL of labeling enzyme were added to the mixture. The labeling reaction was incubated for 1 h at 16°C and terminated by incubation for 15 min at 65°C. The Hy3 TM -labeled samples were then hybridized on the miRCURY TM LNA Array (v.18.0) (Exiqon) according to the array manual. The entire samples (Hy3 TM -labeled samples in 25 μL of hybridization buffer) were denatured for 2 min at 95°C, incubated on ice for 2 min, and then hybridized to the microarray for 16-20 h at 56°C in a 12-Bay Hybridization System (Nimblegen Systems, Inc., Madison, WI, USA), which provides active mixing action at a constant incubation temperature to improve hybridization uniformity and enhance signals. Following hybridization, the slides were washed several times using the Wash buffer kit (Exiqon) and dried by centrifugation for 5 min at 400 rpm. The slides were then scanned using the Axon GenePix 4000B microarray scanner (Axon Instruments, Foster City, CA). The scanned images were imported into GenePix Pro 6.0 software (Axon) for grid alignment and data extraction. Replicate microRNAs were averaged, and microRNAs with intensities ≥30 in all samples were used to calculate the normalization factor. The expression data were normalized using median normalization. After normalization, differentially expressed microRNAs were identified through fold-change filtering. Hierarchical clustering was performed using MEV software (v4.6, TIGR).

MicroRNA-mRNA interaction network
The experimentally verified miR2Disease (18927107), miRecords (18996891), TarBase (22135297), and miRTarBase (24304892) databases were used to identify microRNA-mRNA interaction networks; 43,260 non-redundant, validated microRNA-mRNA pairs were identified. Among these interactions, 5,092 microRNA-mRNA pairs were associated with 62 of the 179 differentially expressed microRNAs (there were no associations for the remaining 107) and incorporated into the microRNA-mRNA interaction network. Cytoscape software (version 3.2.1) (14597658) was used to construct and illustrate the network.

Functional enrichment analysis
A hypergeometric test was used to calculate the significance of the enriched genes in KEGG pathways. If the whole genome had a total of N genes, of which K were involved in the biological pathway under investigation, and a total of M microRNA target genes, of which x were involved in the same pathway, were being analyzed, then the p value for the enrichment of that pathway was calculated as follows: The pathways with p<0.05 were considered significantly enriched.

Quantitative real-time polymerase chain reaction (qRT-PCR)
Total RNA was extracted from DHA-treated and vehicle-treated PANC-1 and BxPC-3 pancreatic cancer cells using the Ultrapure RNA Kit (CWBio Co., Ltd.) according to the manufacturer's instructions. The firststrand cDNA was synthesized using the HiFi-MMLV cDNA Kit (CWBio Co., Ltd.) for mRNA and using the miRNA cDNA Kit (CWBio Co., Ltd.) for miRNA. Realtime PCR was performed with the UltraSYBR Mixture (with ROX) (CWBio Co., Ltd.) for mRNA and the miRNA Real-Time PCR Assay kit for miRNA (CWBio Co., Ltd.) to quantify relative mRNA and miRNA expression. The reactions were run on a 7500 FAST Real-Time PCR System (Applied Biosystems). The PCR primer pairs are listed in Table 1. The reverse primers for miRNA were provided in the miRNA Real-Time PCR Assay kit. The relative expression of each gene was calculated using the comparative fold-change method (2 -ΔΔct ). β-actin and U6 were used as internal controls for mRNA and microRNA expression, respectively. www.impactjournals.com/oncotarget

Western blot
Proteins were extracted measured with western blots as previously described [36][37][38]. Cells that had been treated with 50 μM DHA or vehicle were harvested, washed twice in ice-cold PBS, sonicated in RIPA buffer (Beyotime Institute of Biotechnology, Beijing, China), and homogenized. The debris was removed by centrifugation at 12,000 g for 10 min at 4°C, and protein concentrations were determined using the BCA Protein assay according to the manufacturer's instructions. Samples containing equal amounts of protein (50 μg) were separated by electrophoresis on 10% or 15% polyacrylamide SDS gels (100 V for 1 to 2 hours) and transferred to polyvinylidene difluoride (PVDF) membranes by electroblotting (100 V for 1 hour at 4°C). Running and transfer times and voltages were altered slightly as needed to optimize results. The membranes were then blocked by incubation with 5% skim milk in TBST buffer (TBS plus 0.1% Tween 20) for 2 hours and then incubated with the appropriate primary antibody overnight at 4°C with gentle agitation. The membranes were then washed several times and incubated with the appropriate horseradish peroxidase-conjugated secondary antibody (Santa Cruz Biotechnology, Carlsbad, CA, USA) for 1 hour at room temperature. Finally, membranes were washed again and protein bands were visualized with an enhanced chemiluminescence (ECL) kit, followed by exposure to X-ray film. β-actin was simultaneously measured as a loading control.

Quantification of apoptosis in tumor sections
Apoptosis was quantified as previously described [36,39]. Briefly, paraffin-embedded tissue sections (5 μm) were stained with the TUNEL agent (Roche, Shanghai, China), and numbers of TUNEL-positive cells were counted in ten randomly selected microscopic fields at ×400 magnification.

Tumor microvessel density
Tumor microvessel densities were quantified as previously described [36,39]. Briefly, tumor sections were immunostained with anti-CD31 antobody, numbers of microvessels were counted in randomly selected ten microscopic fields at ×200 magnification, and microvessel densities were calculated.

Statistical analysis
Results are expressed as means ± standard deviations (SD). The significance of differences between histopathologic scores was assessed using the Kruskal-CWallis test. Continuous data were analyzed by ANOVA and the Student-CNewman-CKeuls test. The hypergeometric test was used to identify significantly enriched pathways. Differences were considered statistically significant when p<0.05. All statistical analyses were performed using the R 3.1.3 framework.