Introduction

Cyclooxygenase (COX)-2 is an inducible COX whose expression is promoted by a variety of factors such as cytokines, growth factors and tumor promoters. COX-2 is widely upregulated in many human cancers, including colorectal,1 prostate,2 breast,3 gastric,4 pancreatic,5 lung,6 and head and neck cancers.7 Recent studies have highlighted the relevance of COX-2 in human carcinogenesis using carcinogen-induced animal models, indicating its role in promoting tumorigenesis.8, 9 COX-2 was found to promote cell proliferation and inhibit apoptosis by mediating the activation of downstream oncogenic pathways.10, 11

Hepatocellular carcinoma (HCC) remains the most common primary cancer of the liver. The prognosis of HCC is poor with mortality almost equaling incidence with limited effective treatment options.12 HCC development is a complex process associated with an accumulation of molecular changes that run through the steps of initiation, promotion and progression. Identification of the central factor and molecular behavior of HCC can be of high validity in the effective chemotherapy of this malignancy. High COX-2 expression has been identified in HCC tumors, whereas no or very low COX-2 expression is observed in normal liver tissues.13 COX-2 can induce proliferation and angiogenesis via p53, p27 and vascular endothelial growth factor and also inhibit apoptosis by inducing the antiapoptotic factor Bcl-2, as well as activating antiapoptotic signaling through protein kinase B pathway in HCC.14, 15, 16, 17 Thus, overexpression of COX-2 may be one of the leading factors in hepatic carcinogenesis. However, it is not known whether COX-2 overexpression in liver is sufficient to initiate and drive the formation of spontaneous HCC in vivo. We have previously reported that constitutive COX-2 expression in hepatocytes induced hepatitis in mice.18 In this study, we determined the consequence of liver-specific overexpression of COX-2 in initiating HCC using transgenic (TG) mice and further elucidated the potential mechanisms of action by which COX-2 contributes to the development of HCC.

Results

COX-2 overexpression promoted spontaneous HCC formation

COX-2 protein was confirmed to be highly expressed in the hepatocytes from the conditional TG mice but not in the wild-type (WT) mice by western blot (Figure 1a). After 24 months of normal chow feeding, the spontaneous HCC formation was only found in male TG mice (3/21, 14.3%), but not in any of the WT male mice (0/19) at harvesting. Histological examination confirmed HCC formation (well differentiated) in the liver tissues of TG mice, along with co-occurrence of inflammatory infiltration and neovessels (Figure 1b). Moreover, liver-infiltrating macrophages were identified in all TG mice (3/3) with HCC but not in WT mice (0/3) by F4/80 immunostaining (Figure 1c). To further confirm the role of COX-2 in HCC development, three different HCC cell lines (HepG2, Huh7 and SK-Hep1) were treated with the COX-2 selective inhibitor celecoxib. HCC cells treated with celecoxib significantly inhibited cell growth in dose-dependent growth (Figure 1d).

Figure 1
figure 1

Spontaneous HCC formation in liver of COX-2 TG mice. (a) Western blot of livers from WT and TG mice. The cell lysates were blotted with anti-COX-2 or anti-β-actin antibody. The arrow indicates the position of COX-2. (b) Macroscopic photographs and hematoxylin and eosin (H&E) staining (× 100 magnification) of representative images of livers from WT and COX-2 TG mice. The liver from WT mice appeared normal. Liver from TG mice showed formation of well-differentiated HCC (gray arrow or T), inflammatory infiltration (blue arrow) and neovessels (white arrow). (c) F4/80 staining of representative pictures of livers from WT and COX-2 TG mice with HCC formation (× 200 magnification). The presence of macrophages was shown as brown-staining (black arrow). Average number of macrophages were counted in different areas (at least five) in each sample (n = 3 for both genotypes). (d) Celecoxib inhibited cell growth of different HCC cell lines in a dose-dependent manner. The cell lines were exposed to different concentrations of celecoxib for 5 days. Data are represented as mean±s.d. P-values were obtained by independent Student'st-test. (*P<0.05, **P<0.01).

Expression microarray analysis identified oncogenic pathways induced in COX-2-driven HCC

To explore transcriptome-wide signatures altered in COX-2-induced HCC, we analyzed the gene expression profiles of HCC tumors from TG mice and livers from TG mice without HCC by gene expression array covering 39 430 entrez genes. We identified 1742 differentially expressed genes (DEGs), which could discriminate HCC from non-HCC tissues by unsupervised clustering (Figure 2a). Gene set enrichment analysis (GSEA) against 50 hallmark gene sets was first performed on the DEGs. We found that inflammatory response pathway and nuclear factor (NF)-kB pathway in response to tumor necrosis factor α (TNFα) were among the most enriched gene sets in COX-2-induced HCC (Supplementary Figure S1). In addition, expression of several inflammation-related genes, such as CD14, CD40, CXCL10 and interleukin-6, were significantly increased (Supplementary Figure S1).

Figure 2
figure 2

COX-2 overexpression increase the expression of target genes of AKT, STK33 and MTOR pathways in HCC of TG mice. (a) The expression of 1742 DEGs (twofold change and adjusted P-value <0.05) in three HCCs from different TG mice compared with three livers of TG mice (w/o HCC) are shown in a heat map (the green, black and red colors represent lower than average, close to average and higher than average expression value, respectively). Unsupervised hierarchical clustering analysis was applied. (b) GSEA enrichment analysis of oncogenic signature for 1742 DEGs. Normalized enrichments score (NES) and nominal P-value are shown per each gene set analyzed. Sixteen gene sets are significantly enriched at nominal P-value <0.05. GSEA enrichment plots shown in right panel demonstrates upregulation of AKT, STK33 and MTOR pathway (top gene sets) signatures in COX-2-induced HCC. (c) Western blot of three normal livers from WT mice, three livers from TG mice (w/o HCC) and three HCC from TG mice. The tissue lysates were blotted with indicate antibody. Levels of proteins were quantified by Image J (National Institutes of Health, Bethesda, MD, USA). The bar graph represents relative expression of each protein with mean values from WT group set as 1.0. Data are represented by mean±s.d. P-values were obtained by independent Student'st-test (*P<0.05, **P<0.01). (d) The graph shows each gene and the number of subsets in which it appears by leading edge analysis. (e) The mRNA expression of Hbegf, Krt23, Pak1 and TNFRSF12A in HCCs from TG mice (n=3) compared with normal livers from WT mice (n=3) or livers from TG mice without HCC (n=3). Data are represented as mean±s.d. P-values were obtained by independent Student'st-test. (*P<0.05, **P<0.01).

To specifically investigate oncogenic pathways activated in COX-2-induced HCC, we next performed GSEA against the 189 oncogenic signature gene sets. Interestingly, AKT (protein kinase B), STK33 (Serine/Threonine kinase 33) and MTOR (mechanistic target of rapamycin) pathways signature genes were found to be significantly upregulated (Figure 2b and Supplementary Figure S2), suggesting that these three oncogenic pathways could be activated in the COX-2-induced HCC. In keeping with this, increased protein expression of phosphorylated mTOR was found in COX-2-induced HCCs compared with the non-HCC controls (Figure 2c). We further performed gene screening among the 16 significantly enriched oncogenic signature gene sets in COX-2-induced HCC by leading edge analysis (Figure 2d). We found that, of the 14 core signature genes that were shared by at least three gene sets, HB-EGF, Krt23, Pak1 and TNFRSF12A were highly expressed in HCC compared with the normal liver of WT mice or liver from TG mice without HCC by quantitative reverse transcriptase–PCR (Figure 2e). Therefore, these genes may have important oncogenic role in COX-2-driven HCC.

Genome-wide DNA methylation analysis revealed epigenetic alteration in COX-2-induced HCC

To investigate the epigenetic alternations in COX-2-induced HCC, we performed reduced representation bisulfite sequencing (RRBS) and identified 3717 differentially methylated regions (DMRs) in HCC tumors. Interestingly, the majority of them (3430/3717, 92.2%) were found to be hypermethylated DMRs (Figure 3a). The DMRs were distributed mainly at the intron (24.80%), coding sequence (CDS) (21.87%) and promoter (18.70%) regions (Figure 3b). Of note, DMRs located at promoter regions showed an approximately twofold increase in methylation level in HCC as compared with non-HCC (Figure 3c). To investigate the epigenetic regulation of signaling pathways in COX-2-induced HCC, we performed Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis for 3717 DMRs and 1742 DEGs independently. Five signaling pathways (pathways in cancer, axon guidance, focal adhesion, amebiasis, gap junction and the regulation of actin cytoskeleton) were found to be commonly altered in COX-2-induced HCC (Figure 3d and Table 1).

Figure 3
figure 3

Characterization of DNA methylation in COX-2-induced HCC of TG mice. (a) The methylation levels of 3717 DMRs (1.25-fold change and false discovery rate <5%) in three HCC vs three livers (w/o HCC) and from different TG mice are shown in a heat map (the blue, green and yellow colors represent lower than average, close to average and higher than average methylation level, respectively). Unsupervised hierarchical clustering analysis was applied. (b) DNA methylation patterns in different genomic regions, including upstream, 5’-UTR, promoter, coding sequence (CDS), 3’-UTR, intron and downstream. (c) Scatter plot of 695 DMRs at promoter region in three HCC vs three livers (w/o HCC) from different TG mice. The data are represented by median with interquartile range. (d) Enrichment of Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways analysis for 3717 DMRs and 1742 DEGs. Pathways with 10 or more genes involved and a P-value <1e-06 are considered as significantly affected. Gene number involved in different pathway is indicated above bars.

Table 1 Overlap of predicted pathways by different source (DEG and DMR)

Integrative analyses of gene methylation and mRNA expression in COX-2-induced HCC

Given that promoter hypermethylation has a key role in tumorigenesis via the repression of tumor-suppressor genes, we hypothesize that promoter methylation may contribute to COX-2-driven HCC by altering their gene expression. By integrative analyses of DNA methylation genes in correlation with their mRNA expression, we identified that 68 genes (for example,ADCY5, AXIN2, LTBP1, NDK1, NTN1, PLXNC1 and PRKCZ) were hypermethylated and transcriptionally downregulated, but only 12 genes (for example,Pygb, Pitpnm1 and Ssbp4) that were demethylated and transcriptionally upregulated, in COX-2-induced HCC (Figure 4a). We further validated the significant downregulation of mRNA expression of seven randomly selected genes in COX-2-induced HCC as compared with normal liver from WT mice or liver from TG mice without HCC (Figure 4b). By bisulfite genomic sequencing, ADCY5, LTBP1, NKD1 and PLXNC1 were confirmed to be hypermethylated at the promoter region in COX-2-induced HCC as compared with WT liver in mice (Figure 4c).

Figure 4
figure 4

Identification and validation of aberrantly methylated genes involved in COX-2-induced HCC in mice. (a) Dot plots of genes of differential expression and methylation (promoter) in three HCC vs three livers (w/o HCC) from different TG mice. Several candidate genes with differential expression (delta CT value) and altered methylation level (delta M value) were listed in the right panel. (b) The mRNA expression of seven candidate genes (Adcy5, Axin2, Ltbp1, Nkd1, Ntn1, Plxnc1 and Prkcz) in three normal livers from WT mice, three livers from TG mice (w/o HCC) and three HCC from TG mice. Data are represented by mean±s.d. P-values were obtained by independent Student's t-test. (*P<0.05, **P<0.01). (c) Bisulfite genomic sequencing of Adcy5, Ltbp1, Nkd1 and Plxnc1 in HCCs from TG mice compared with normal livers from WT mice. The location of CpG island and bisulfite genomic sequencing target regions were shown for each gene, with black bars denoting the first exons. TSS, transcription start site. Average methylation levels at each site are shown. P-values were obtained by repeated-measures analysis of variance.

DNA promoter hypermethylation was associated with the loss of TET1 expression and activity in COX-2-induced HCC

To understand the causality between COX-2 overexpression and DNA hypermethylation, we screened the key regulators involved in DNA methylation and demethylation process such as DNMT1, DNMT3a, DNMT3b, tet methylcytosine dioxygenase 1 (TET1), TET2 and TET3. Among them, TET1 was found to be significantly reduced both at mRNA levels and protein levels in COX-2-induced HCC (Figures 5a and b). TET1 catalyzes the conversion of 5-methylcytosine (5mC) into 5-hydroxymethylcytosine (5hmC) and prevents DNA hypermethylation.19 Consistent with the reduced TET1 expression, loss of TET1 activity was evidenced by the decreased genomic levels of 5hmC (Figure 5c). These results suggested that COX-2 overexpression may cause reduction of TET1 expression and subsequent loss of TET1-associated 5hmc abundance, which may impair DNA demethylation and result in DNA hypermethylation in COX-2-induced HCC. To confirm this, we constructed two stable HCC cell lines (Hep3B and SK-Hep1) with COX-2 overexpression. We found that upregulation of COX-2 could decrease TET1 mRNA expression (Figure 5d), protein expression (Figure 5e), as well as 5hmC levels (Figure 5f) in these two cell lines. Conversely, HCC cells (HepG2, Huh6 and SK-Hep1) treated with the COX-2 inhibitor celecoxib showed increased TET1 expression in dose-dependent manner (Figure 5g and Supplementary Figure S3A). In our study, enhanced hepatic macrophage infiltration and NF-kB activation are demonstrated in HCCs of COX-2 TG mice. TNFα is a proinflammatory cytokine, which activates canonical NF-kB pathway. We investigate whether the inflammatory pathway could affect the TET1 expression by treating HCC cells (Hep3B, HepG2, Huh6, Huh7 and SK-Hep1) with TNFα. We found decreased TET1 expression in a dose-dependent manner in these HCC cell lines by TNFα treatment (Supplementary Figure S3B), inferring the inflammatory pathway has a role in TET1 downregulation.

Figure 5
figure 5

COX-2 overexpression reduced the expression and activity of TET1. (a) Comparison of genes that are involved in DNA methylation and demethylation process between three HCC vs three livers (w/o HCC) from different TG mice (data was acquired from gene expression array). The fold change (FC) ratio (HCC vs non-HCC) and P-value were list in the top panel. The bottom panel showed the mRNA expression of TET1 in three normal livers from WT mice, three livers from TG mice (w/o HCC) and three HCC from TG mice by using real-time PCR. Data are represented by mean±s.d. P-values were obtained by independent Student's t-test. (*P<0.05, **P<0.01). (b) Western blot of three normal livers from WT mice, three livers from TG mice (w/o HCC) and three HCC from TG mice. The tissue lysates were blotted with indicate antibody (top panel). Levels of proteins were quantified by Image J. The bar graph (bottom panel) represents relative expression of each protein with mean values from WT group set as 1.0. Data are represented by mean±s.d. P-values were obtained by independent Student's t-test. (*P<0.05, **P<0.01). (c) Comparison of global 5hmC quantification in three normal livers from WT mice, three livers from TG mice (w/o HCC) and three HCC from TG mice. P-values were obtained by independent Student's t-test. (*P<0.05, **P<0.01). (d, f) In two different stable COX-2-overexpressing cell lines (Hep3B and SK-Hep1), the mRNA expression (d) and protein expression (e) of TET1 were measured. Data are represented by mean±s.d. (n⩾3). P-values were obtained by independent Student's t-test. (**P<0.01, ***P<0.001). (f) Global 5hmC levels were assayed by dot blot analysis for two different stable COX-2-overexpressing cell lines (Hep3B and SK-Hep1). Quantification was performed by Image J. (g) Celecoxib treatment increases TET1 expression both at mRNA (top panel) and protein level (bottom panel) in dose-dependent manner. HepG2 was deprived of serum overnight and then treated with different dosage of celecoxib for 12 h. Data are represented by mean±s.d. (n⩾3). P-values were obtained by independent Student's t-test. (***P<0.001). Protein quantification was performed by Image J. (h) Schematic diagram of the molecular events for COX-2 inducing hypermethylated genes.

Validation of candidate genes (Ltbp1, ADCY5 and Prkcz) in human HCC

The mRNA expression of LTBP1, ADCY5 and PRKCZ were significantly decreased in human HCC tumors as compared with adjacent non-tumor tissues or healthy liver (Figure 6a). In keeping with this, downregulation of LTBP1, ADCY5 and PRKCZ were demonstrated in The Cancer Genome Atlas (TCGA) HCC cohort (Figure 6b). Moreover, integrative analyses of DNA methylation and mRNA expression from TCGA data set of human HCCs confirmed that LTBP1, ADCY5 and PRKCZ mRNA expression levels were inversely correlated with methylation level of CpG islands in their promoter regions (Figure 6c). Interestingly, TET1 mRNA level was found to be inversely correlated with the DNA methylation level of these three candidates (Supplementary Figure S3).

Figure 6
figure 6

COX-2-induced hypermethylated genes (ADCY5, Ltbp1 and Prkcz) are novel tumor suppressors inactivated by DNA methylation in human HCC. (a) The mRNA expression level of ADCY5, Ltbp1 and Prkcz in 16 normal livers (NL), 30 pairs of HCC (T) and adjacent non-HCC (NT) liver tissues. Data are represented by mean±s.d. P-values were obtained by Mann–Whitney U-test or Wilcoxon matched-pairs test. (*P<0.05, **P<0.01, ***P<0.001). (b) The mRNA expression level of ADCY5, Ltbp1 and Prkcz in 371 primary HCC (T) and 50 adjacent non-HCC (NT) liver tissue from TCGA LIHC (liver HCC) data set. Data are represented by mean±s.d. P-values were obtained by Mann–Whitney U-test. (*P<0.05, **P<0.01, ***P<0.001). (c) The Spearson correlation between mRNA expression level and methylation level of ADCY5, Ltbp1 and Prkcz in 371 primary HCC samples from TCGA data. The gene methylation levels were determined as the average beta-values of different probes, which are mapped to CpG island in their promoter regions on genes.

To evaluate the functional significance of COX-2-mediated silenced genes, we performed cell growth assays on LTBP1, one randomly selected candidate. We found that LTBP1 knockdown (Figure 7a) significantly promoted cell proliferation in HCC cell lines (SK-Hep1 and Hep3B) as evidenced by cell viability assay (Figure 7b). The induction of cell proliferation by LTBP1 knockdown was further confirmed by colony formation assay in these cell lines (Figure 7c), inferring the tumor-suppressive role of LTBP1 in COX-2-induced HCC. These findings collectively suggest that overexpression of COX-2 could induce promoter hypermethylation, which directly leads to the transcriptional silencing of potential tumor-suppressive genes.

Figure 7
figure 7

Knockdown of LTBP1 promoted HCC cell proliferation. (a) Hep3B and SK-Hep1 cells were transfected with negative control small interfering RNA (siRNA-NC), Ltbp1 target siRNA1 or siRNA2 for 4 days. Ltbp1 knockdown efficiency in different cell lines were validated by q-PCR (top panel) and reverse transcriptase (RT)–PCR (bottom panel). Data are represented by mean±s.d. P-values were obtained by independent Student's t-test. (**P<0.01, ***P<0.001). (b) The cell viability was measured at day 4 by MTT assay. P-values were obtained by independent Student's t-test. (*P<0.05, ***P<0.001). (c) Knockdown of Ltbp1 by Ltbp1 target siRNA2 promoted colony formation in Hep3B and SK-Hep1 cells. (***P<0.001).

Discussion

COX-2 is an inducible COX isoform that can be regulated by various growth factors and inflammatory cytokines such as interleukin-1β and TNFα.20 Therefore, COX-2 is frequently overexpressed during inflammation. In human HCC, increased COX-2 expression had been found in the early stages (usually well differentiated), indicating its causative role in hepatocarcinogenesis.13 However, there have been several studies reporting that COX-2 was upregulated in the noncancerous liver tissues as compared with HCC tissues in patients with cirrhosis or chronic hepatitis.13, 21 In contrast, COX-2 expression was not detected in normal liver without chronic inflammatory diseases.22 COX-2 has been reported to induce tumor initiation, progression and angiogenesis in different cancer types.23 To date, there has been no report on whether COX-2 overexpression can directly drive HCC. Previously, we demonstrated that liver-specific COX-2 TG mice spontaneously developed hepatitis.18 Here, we showed for the first time that COX-2 transgene is sufficient to initiate HCC in mice.

We constructed TG mice that specifically expressed human COX-2 complementary DNA in hepatocytes using the transthyretin promoter.24 After 24 months, spontaneous HCC was identified only in male TG mice (3/21, 14.3%), but not in WT mice. Histopathological examination revealed the presence of numerous foci of HCC, inflammation and neovessels in the liver of TG mice (Figure 1b). We have demonstrated that enhanced COX-2 expression in hepatocytes in TG mice could activate NF-kB pathway, stimulate the secretion of proinflammatory cytokines and enhance the recruitment of macrophage.18 Consistently, liver-infiltrating macrophage and activation of NF-kB pathway were observed in HCC of TG mice (Figure 1c and Supplementary Figure S1). Celecoxib is a COX-2 selective inhibitor, which is used for patients with pain and inflammation. Celecoxib treatment in HCC cell lines significantly decreased the COX-2 expression and inhibited cell growth (Figures 1d and 5g) further confirmed the role of COX-2 in promoting HCC development.

To investigate and elucidate the mechanisms of COX-2 on aberrant gene expression and pathways, we analyzed the transcriptome of COX-2-induced HCC as compared with non-HCC tissues. By GSEA analysis, AKT, STK33 and MTOR pathway signature genes were found to be significantly upregulated in COX-2-driven HCC (Figure 2b and Supplementary Figure S2), indicating the activation of these oncogenic cascades by COX-2 overexpression during the pathogenesis of HCC. Enhanced activation of mTOR was further confirmed in COX-2-induced HCCs compared with non-HCC controls by western blot (Figure 2c). Consistent with our observations, COX-2-derived prostanoid was reported to exert an oncogenic effect through activation of AKT or mTOR-regulated pathways, thereby promoting cell growth and survival.14, 25 TET1 activates the transcription of PTEN via directly binding to its promoter region and causing demethylation of its CpG islands in cancer.26 In keeping with this, we found that PTEN expression was silenced in the HCC compared with TG mice livers (P=0.034; Supplementary Table 3). Downregulation of TET1 therefore may activate AKT and mTOR pathways that are suppressed by PTEN. On the other hand, it was reported that inhibition of AKT pathway using phosphatidylinositol 3 kinase inhibitor had no effect on TET1 expression in vitro.27 These findings suggested that downregulation of TET1 reduces the expression of PTEN, which subsequently activate AKT and MTOR pathway in COX-2-driven HCC.

In addition, four core signature genes (HB-EGF, Krt23, Pak1 and TNFRSF12A) that shared across the enriched gene sets were subsequently confirmed to be highly elevated in COX-2-driven HCC compared with the normal livers (Figures 2d and e). HB-EGF has been reported to stimulate cell growth in either autocrine or paracrine manner, which contributes to tumorigenicity;28 KRT23 is highly expressed in colon adenocarcinomas, and knockdown of KRT23 reduced proliferation of colon cancer cells;29 Pak1 is a potent mediator for prostate cancer cell migration and tumor growth, and knockdown of Pak1 impaired prostate tumor growth;30 By interaction with TNFRSF12A protein, Tweak protein exerted mitogenic effect in promoting proliferation of oval cells.31 These results supported that COX-2 overexpression enhanced the expression of key oncogenic genes and activated the oncogenic pathways, which contributes to the initiation and progression of HCC formation in COX-2 TG mice.

Epigenetic silencing by DNA methylation is a predominant mechanism that inactivates genes involved in HCC tumorigenesis.32 A large body of evidence supports that DNA hypermethylation surrounding the proximal promoter region leads to transcriptional inactivation of putative tumor-suppressor genes.33 In this study, we demonstrated for the first time the aberrant promoter DNA hypermethylation in COX-2-induced HCC in TG mice, which was associated with downregulation of gene expression. Methylation levels of DMRs surrounding the promoter region in COX-2-induced HCC were found to be twofold higher than that in non-HCC livers (Figure 3c). Indeed, we identified 68 genes (for example,ADCY5, AXIN2, LTBP1, NDK1, NTN1, PLXNC1 and PRKCZ) that were transcriptionally downregulated by promoter hypermethylation in COX-2-induced HCC (Figure 4). It is recognized that DNA methyltransferases (DNMT1, DNMT3A and DNMT3B) and ten-eleven translocation (TET) family (TET1, TET2 and TET3) are responsible for the regulation of the abundance and distribution of DNA methylation.34 In our study, only TET1 was found to be significantly reduced at both mRNA and protein levels in COX-2-induced HCC (Figures 5a and b). Consistent result was also obtained in the paired human HCC tissues that TET1 is significantly lower in HCC compared with their adjacent non-tumor tissues (Supplementary Figure S3C). In agreement with our data, protein level of TET1 protein, but not TET2 and TET3, was reported to be downregulated in human HCC as compared with adjacent non-tumor tissues.35 Protein level of 5hmC was significantly reduced, whereas the 5mC level was increased in HCC tissues.35 Moreover, reduced Tet1 activity was associated with aberrant hypermethylated CpG islands during hepatocarcinogenesis.36 In our study, COX-2 overexpression in HCC cell lines significantly reduced both TET1 expression and 5hmC level (Figures 5d–f). Conversely, treatment of celecoxib significantly increased TET1 expression (Figure 5g and Supplementary Figure S3A). Furthermore, HCC cells treated with TNFα showed decreased TET1 expression (Supplementary Figure S3B) indicating the correlation between inflammatory pathway activation and TET1 downregulation. These results collectively suggested that promoter hypermethylation in COX-2-induced HCC was mediated, at least in part, by reducing TET1 expression and inhibiting the conversion of 5mc to 5hmc, which results in subsequent inhibition of tumor-suppressor genes by promoter hypermethylation (Figure 5h). By analyzing TCGA-LIHC data set, TET1 mRNA expression was found to be inversely correlated with the DNA methylation level of three candidate genes (LTBP1, ADCY5 and PRKCZ) (Supplementary Figure S4), which were inactivated by DNA hypermethylation in both COX-2-driven HCC of TG mice and human HCCs (Figures 4 and 6). Of these three candidates, LTBP1 was functionally found to promote cell proliferation in different liver cell lines as a potential tumor suppressor (Figure 7). LTBP1 belongs to the family of latent TGF-β binding proteins, which could regulate the secretion and activation of TGF-β.37, 38 Loss of LTBP1 was reported to reduce TGF-β activity in primary hepatic stellate cells.39 In contrast, overexpression of LTBP1 increased TGF-β activation and reduced cell proliferation in mouse embryo fibroblasts.40 Thus, LTBP1 may regulate cancer cell proliferation by modulating TGF-β signaling.41 In hepatocytes, TGF-β induced cell growth arrest mediated by the interaction of Smad proteins with Sp1 transcription factors.42 In addition, TGF-β was also reported to induce hepatocyte apoptosis.43 In this connection, reduced LTBP1 expression in hepatocytes of COX-2 TG mice could promote the HCC development via abrogating the cytostatic effects of TGF-β signaling pathway. All these results imply that COX-2 overexpression may silence tumor suppressors by DNA methylation, which in turn, promote HCC formation.

In conclusion, our study demonstrated that enhanced COX-2 expression in hepatocytes is sufficient to induce HCC. The oncogenic effect of COX-2 is mediated by inducing important oncogenic factors (HB-EGF, Krt23, Pak1 and TNFRSF12A) and signaling cascades (AKT, STK33 and MTOR pathway) and by causing transcriptional silence of tumor-suppressive genes (LTBP1, ADCY5 and PRKCZ) through promoter hypermethylation mediated by the reduction of TET1 expression. Inhibition of COX-2 may provide a potential approach for prevention and treatment of liver cancer.

Materials and methods

Generation of TG mice

COX-2 TG mice were generated as described previously.18 COX-2 TG mice, as well as WT littermates were maintained under normal condition. At 24 months of age, the animals were killed under anesthesia, and the livers were harvested, fixed and embedded for histopathological examination. A portion of the normal liver and HCCs were stored at −80 °C for further analysis. All protocols and procedures were approved by the Animal Experimentation and Ethics Committee of the Chinese University of Hong Kong.

Histopathological analysis

Five μm thick sections from paraffin-embedded mouse liver tissues were deparaffinized, hydrated, and then stained by hematoxylin and eosin. Hematoxylin and eosin staining was performed as described previously:44 Hematoxylin solution for 5 min, hydrochloric acid alcohol solution for 15 s, eosin solution for 30 s and 90% ethanol for 30 s. Neutral balsam was used for mounting and the section was observed and photographed under the microscope. HCC was confirmed histologically in a blinded manner by two pathologists.

Human HCC samples

Thirty paired HCC samples (primary tumor and adjacent non-tumor sites) and 16 normal liver biopsies were obtained at the Prince of Wales Hospital, The Chinese University of Hong Kong from 2010 to 2014. All specimens were immediately snapped frozen in liquid nitrogen and stored at −80 °C until further processing. Informed consent was provided by all participants, and this study was approved by the ethics committee of the Chinese University of Hong Kong.

Gene expression microarray

Total RNA was extracted from the mouse liver tissues and hybridized to the SurePrint G3 Mouse GE 8x60K Microarray (Agilent Technologies, Santa Clara, CA, USA). Gene expression microarray raw data were normalized using GeneSpring GX software (Agilent Technologies). Signal values for all genes were transformed to the log base 2. Quantile normalization was applied to obtain equal distributions of the probe signal intensities. The comparative analysis between the three HCC and three non-HCC livers from different TG mice was carried out usingt-test and the Benjamini–Hochberg false discovery rate correction. Differential expression genes were determined if their adjusted P-values were <0.05 and the fold change was >2 or <0.5.

Unsupervised hierarchical clustering

Hierarchical clustering analysis was performed in R (www.r-project.org) using ‘hclust’ function. Expression values (or methylation level) of the different genes were mean-centered, followed by hierarchical clustering using one minus Pearson correlation as distance measure and average linkage.

Gene set enrichment analysis

GSEA (version 2.2.1, Broad Institute, Cambridge, MA, USA) was used to identify significantly enriched functional gene sets.45 DEGs generated above were inputted into GSEA and analyzed using hallmark gene sets or oncogenic signatures database to identify dysregulated cellular pathways in COX-2-induced HCC. The enrichment score was calculated by member genes ranking at the top or bottom of the ranked gene list with 100 permutations. We considered the gene sets that with nominal P-value <5% and false discovery rate <25% as associated with liver tumorigenesis by COX-2 overexpression. The maximum gene set size was fixed at 500 genes, and the minimum size fixed at 15 genes. The leading edge analysis was performed to identify the genes that are shared across the enriched gene sets.

RRBS library preparation

Genomic DNA was isolated from snap-frozen mouse liver tissues and digested with 20 units of MspI (New England Biolabs, Ipswich, MA, USA). After purification, the blunt-end DNA fragments were added with dA, followed by methylated-adapter (Illumina, San Diego, CA, USA) ligation. The libraries were size selected for fractions in 40–120 bp and 120–220 bp range on a 2% agarose gel, and then subjected to bisulfite conversion using EpiTect Bisulfite Kit (Qiagen, Valencia, CA, USA). The final libraries were amplified using HiFi HotStart Uracil+ ReadyMix (Kapa Biosystems, Wilmington, MA, USA) and Illumina Multiplexing PCR Primers. Quality and quantity of RRBS libraries were analyzed by Agilent 2100 Bioanalyzer (Agilent Technologies) and reverse transcriptase–PCR, respectively.

RRBS data processing and mapping

RRBS libraries were sequenced using the Illumina HiSeq2000 platform. All reads were trimmed for adaptor sequences and low-quality sequences using trimgalore (http://www.bioinformatics.babraham.ac.uk/projects/trim_galore/). Bisulfite alignment of the RRBS reads to the mouse reference genome mm9 was performed.46 Features of Repeat Masker annotations, CpG islands and refGene coding loci were downloaded from UCSC.47 The promoter region was defined as the sequence from −1 kb to +1 kb of transcriptional start site.

Differential methylation analysis

DMRs were identified by comparing the methylomes of three HCC and three non-HCC livers from different TG mice using our previous developed software swDMR.48 DMRs were defined as those with at least five adjacent CpG sites, 1.25-fold change in the methylation level and analysis of variance adjusted P-value <0.05.

Enrichment analysis of Kyoto encyclopedia of genes and genomes pathways

Genes dysregulated by epigenetic changes and differential expression levels were included in the Kyoto encyclopedia of genes and genomes pathway enrichment analysis using the Gene Set Analysis Toolkit V2 (http://bioinfo.vanderbilt.edu/webgestalt/). The hypergeometric test statistical method and the Benjamini-Hochberg multiple test adjustment method were used. Mus musculus genome was used as reference. Pathways with at least 10 genes and adjusted P-values of <1e-06 were considered significantly enriched.

Real-time quantitative PCR

RNA was extracted using Trizol reagent (Invitrogen, Carlsbad, CA, USA). Complementary DNA was synthesized using Transcriptor Reverse Transcriptase (Roche Applied Sciences, Indianapolis, IN, USA). Quantitative reverse transcriptase–PCR was performed on a LightCycler 480 real-time PCR system (Roche Applied Science) using the LightCycler 480 SYBR Green I Master (Roche Applied Science). The primers used were listed in Supplementary Table 1. All reactions were done in triplicate. The relative expression level was normalized to β-actin and calculated using 2-ΔΔCt Method.

Bisulfite genomic sequencing

CpG islands were predicted using the CpG island searcher (http://cpgislands.usc.edu/). Bisulfite genomic sequencing primers were designed using MethPrimer (http://www.urogene.org/cgi-bin/methprimer/methprimer.cgi). Briefly, 2 μl of bisulfite-treated DNA by sodium metabisulfite was amplified by BGS primers (Supplementary Table 2). Sequencing was performed using the BigDye Terminator Cycle Sequencing kit version 1.0 (Applied Biosystems, Foster City, CA, USA) and analyzed by using SeqScape software (Applied Biosystems). Methylation level at each CpG site was calculated as the height ratio of cytidine/(cytidine+thymine).

Detection of global 5hmc

Genomic DNA was isolated from snap-frozen mouse liver tissues. Briefly, 100 ng of genomic DNA was used for 5hmC quantification using Quest 5-hmC DNA ELISA Kit (Zymo Research, Irvine, CA, USA).

Dot blot analysis

Serially diluted DNA (200, 100 , 50 and 25 ng) was denatured at 99 °C for 5 min and cool down on ice. DNA samples were spotted on a nitrocellulose membrane, air-dry and ultraviolet cross-linked for 10 min. Membrane was blocked with 5% bovine serum albumin TBST (1x Tris-buffered saline +0.1% Tween-20) for 30 min at room temperature and then incubated with rabbit anti-5hmC polyclonal antibody (Active Motif, Carlsbad, CA, USA) overnight at 4 °C. The membrane was washed for 5 min three times in TBST, and then incubated with horseradish peroxidase-conjugated goat anti-rabbit IgG (Santa Cruz Biotechnology, Dallas, TX, USA) for 1 h at room temperature. The membrane was then washed for 5 min three times in TBST and visualized by chemiluminescence.

Cell culture

Human HCC cell lines used in this study were cultured in Dulbecco's modified Eagle's medium (Invitrogen) supplemented with 10% fetal bovine serum (Thermo Scientific, Waltham, MA, USA) and maintained at a 37 °C in a humidified incubator with 5% CO2.

Construction and transfection of expression vectors

The pcDNA3.1 expression vector containing the full-length open reading frames of human COX-2 (NM_000963.3) was constructed. The sequence of the construct was confirmed. Liver cell line was transfected with pcDNA3.1 carrying COX-2 complementary DNA or empty vectors using Lipofectamine 2000 (Invitrogen). Stably transfected cells were established under selection with neomycin (G418) (Invitrogen).

Small interfering RNA-mediated gene silencing

The sequences of two Ltbp1-small interfering RNAs (Genepharma, Shanghai, China) are as follows: GCAGGGUCCUAUGAUUGUATT (sense: 5′–3′) and UACAAUCAUAGGACCCUGCTT (anti-sense: 5′–3′); GUACUAUGAUCCUGUGAAATT (sense: 5′–3′) and UUUCACAGGAUCAUAGUACTT (anti-sense: 5′–3′). Liver cell lines were transfected using Lipofectamine 2000 (Invitrogen). Knockdown efficiency was determined 4 days post-transfection by quantitative PCR.

Cell viability and colony formation assay

Cells (1 × 103 per well) were seeded in a 96-well plate for four days. MTT (5 mg/ml; Promega, Madison, WI, USA) was added to the culture to equal one-tenth the original culture volume. After 4-h incubation at 37 °C, 100 μl of a solution containing 10% sodium dodecyl sulfate and 0.01 n HCl was added. Absorbance at 570 nm was measured on a microplate reader. For colony formation array, cells (500 per well) were seeded in 12-well plate for 2 weeks. Colonies were fixed with methanol/acetone (1:1), stained with Gentian violet, and counted.

Western blot analysis

Western blot analysis was performed as previously described.44 The primary antibodies used were COX-2, phosphorylated mTOR and β-actin (Santa Cruz Biotechnology); MTOR (Cell Signaling Technology, Danvers, MA, USA); TET1 (GeneTex, Irvine, CA, USA); DNMT3A and DNMT3B (Abcam, Cambridge, UK).

Statistical analysis

All measurements or variables are shown as mean ±s.d. Mann–Whitney U-test, Wilcoxon matched-pairs test and Student's t-test were performed. All statistical tests were performed using Graphpad Prism 5.0 (GraphPad Software Inc., San Diego, CA, USA), and a two-tailed P-value of <0.05 was considered statistically significant.