Long noncoding RNAs in hepatitis B virus-related hepatocellular carcinoma

Abstract

AIM: To study the expression of long noncoding RNAs (lncRNAs) in hepatitis B virus (HBV)-related hepatocellular carcinoma (HCC).

METHODS: The lncRNA profiles between HBV-related HCC tissues and corresponding normal liver tissues were generated using microarray analysis. Datasets were analyzed using multiple algorithms to depict alterations in gene expression on the basis of gene ontology (GO), pathway analysis, and lncRNA levels.

RESULTS: The microarray revealed that 1772 lncRNAs and 2508 mRNAs were differently expressed. The pathway analysis demonstrated that the cell cycle, cytokine-cytokine receptor interaction, chemokine signaling pathway, and phosphoinositide 3-kinase-protein kinase B signaling pathway may play important roles in HCC. Several GO terms, such as cell cycle, DNA replication, immune response, and signal transduction, were enriched in gene lists, suggesting a potential correlation with HBV-related HCC. The upregulated large intergenic noncoding RNA ULK4P2 was physically combined with enhancer of zeste homolog 2. Therefore, the lncRNAs may participate in regulating HBV-related HCC.

CONCLUSION: lncRNAs play important roles in HCC, future studies should verify whether large intergenic noncoding ULK4P2 functions by combining with enhancer of zeste homolog 2 in HCC.

Key Words: Enhancer of zeste homolog 2; Hepatocellular carcinoma; Long noncoding RNAs; Microarray; ULK4P2

Core tip: This manuscript examines the differential expression of long noncoding RNAs in hepatitis B virus-related hepatocellular carcinoma using gene ontology and pathway analyses, and constructing a long noncoding RNA-mRNA network to research the data. Furthermore, RNA immunoprecipitation revealed that the large intergenic noncoding RNA ULK4P2 physically combined with enhancer of zeste homolog 2 (EZH2). EZH2 plays an important role in many cancer types and is critical for cancer cell proliferation, survival, and metastasis, and drug resistance. The combination of ULK4P2 and EZH2 provides a novel avenue for further study of ULK4P2 in hepatitis B virus-related hepatocellular carcinoma.

INTRODUCTION

Hepatocellular carcinoma (HCC) is a common human cancer in many countries, especially China[1]. The mortality rate of HCC is third among cancer-related deaths[2], and hepatitis B virus (HBV) infection is associated with almost half of all HCC cases. A striking similarity exists between the geographic distribution of the rates of chronic HBV infection and HCC[3]. However, the mechanism of HBV in HCC has not yet been clarified, and the five-year survival rates of HCC patients remain poor. Therefore, understanding the underlying mechanism of HBV in HCC is crucial to provide accurate information for the early screening, clinical diagnosis, targeted molecular therapy, and prognosis of patients.

Although a significant portion of the human genome is transcribed, protein-coding genes account for only 2% of the genomic sequence[4]. In recent years, the nonprotein-coding portion of the genome has become important in the basic biology and major pathologies of cancer[5]. Noncoding RNAs include microRNAs and long noncoding (lnc)RNAs; the microRNAs has been well studied and are related to cell differentiation and cancers in recent publications[6]. With the development of lncRNA microarrays, high-throughput sequencing, and bioinformatics, an increasing number of lncRNAs have been discovered and have attracted considerable attention in medical molecular biology. lncRNAs are a class of noncoding RNA transcripts longer than 200 nucleotides, with no or little protein-coding capacity. Recent studies have shown that more genomic sequences have been transcribed into lncRNAs than protein-coding RNAs[7]. Although lncRNAs are among the least well understood of the noncoding RNAs, they cannot be completely dismissed as mere transcriptional “noise’’[8].

lncRNAs regulate gene expression at the epigenetic, transcriptional, and post-transcriptional levels, and are involved in many biologic functions, including genomic imprinting, chromosome dosage-compensation, X-chromosome silencing, chromosome modification, intranuclear transport, transcriptional activation, and interference[9]. lncRNAs are abnormally expressed in cancer cells[10,11], and several lncRNAs play important roles in carcinogenesis. However, only a small fraction of lncRNAs is currently well understood.

The molecular mechanisms underlying HBV-related HCC are not clear, and the contributions of lncRNAs in HCC have only been gradually elucidated, including of HEIH, which is an lncRNA that is highly and specifically expressed in HCC tissues[1]. The downregulated lncRNA GAS5 is associated with HCC prognosis[12]. The lncRNA ATB, a highly expressed gene in HCC, participates in the epithelial-mesenchymal transition induced in HCC metastases and is associated with poor prognosis. However, the important roles of lncRNAs in HCC have yet to be elucidated. In the present study, lncRNA microarrays were used to detect differentially expressed lncRNAs between HBV-related HCC tissues and nontumor tissues. In addition, gene ontology (GO) analysis, pathway analysis, and an lncRNA-mRNA network were used to predict the functions of these abnormally expressed lncRNAs in HCC.

MATERIALS AND METHODS

Preparation of HCC tissues

Tissue samples from five HBV-related HCC patients were used for the microarray analysis, and tissue samples from 14 HCC patients who had hepatectomy in the Renmin Hospital of Wuhan University between May 2013 and May 2014 were used for data validation. After removal from the body, the tissues were snap frozen in liquid nitrogen within 30 min and then stored at -80 °C until use. All patients provided written informed consent. The study was approved by the Human Research Ethics Committee of Renmin Hospital of Wuhan University. Table S1 lists the characteristics of the patients.

RNA extraction

To extract RNA, frozen tissues were ground into powder with TRIzol reagent (Invitrogen of Thermo Fisher Scientific, Waltham, MA, United States). RNA purification was performed with an RNA-containing aqueous phase using the RNeasy minikit (Qiagen, Venlo, Limburg, Netherlands). Quantification and quality evaluation were performed using a Nanodrop and Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, United States), respectively.

Arraystar human lncRNA Microarray V3.0

Arraystar Human LncRNA Microarray V3.0 (Arraystar, Rockville, MD), an updated version of Microarray V2.0, was designed for the global profiling of human lncRNAs and protein-coding transcripts. Approximately 30586 lncRNAs and 26109 coding transcripts can be detected by the third-generation lncRNA microarray. lncRNAs were carefully constructed using well-respected public transcriptome databases (Refseq, UCSC Known Genes, and Genecode) and landmark publications. Each transcript was represented by a specific exon or splice junction probe that can accurately identify individual transcripts. Positive probes for housekeeping genes and negative probes were also printed onto the array for hybridization quality control. Data were extracted and normalized using GeneSpring GX v11.5.1 software (Agilent Technologies). Differentially expressed lncRNAs with statistical significance were identified through Volcano Plot filtering and hierarchical clustering. Differentially expressed genes were identified through the random variance model, and P values were calculated using the paired t-test. The significance thresholds set for the up- and downregulated genes were fold change ≥ 2.0 and P≤ 0.05.

qRT-PCR

Total RNA was reverse transcribed using a Fermentas RT reagent kit (Perfect Real Time; Thermo Fisher Scientific) according to the manufacturer’s instructions. Four distinctively upregulated lncRNAs were randomly selected to validate their expression levels through quantitative real-time reverse-transcription (qRT)-PCR using SYBR Green assays (TaKaRa, Otsu, Shiga, Japan). GAPDH was used as an internal control. Table S2 shows the primer sequences.

lncRNA-mRNA network

An lncRNA-mRNA network was developed to identify the interactions between mRNA and lncRNA[13]. The network was built according to the normalized signal intensities of specific expression levels of genes and lncRNAs (for each gene-lncRNA, gene-gene, or lncRNA-lncRNA pair). Pearson’s correlation was calculated and significantly correlated pairs were used to construct the network.

Pathway and GO analysis

The main functions of the differentially expressed genes were analyzed using GO analysis, a key functional classification of NCBI, GO can organize genes into hierarchical categories and uncover the gene regulatory network on the basis of biologic processes and molecular functions[14].

Pathway analysis was used to determine the main pathway of the differentially expressed genes according to KEGG, Biocarta, and Reatome. Fisher’s exact test and χ² tests were used to select the main pathway, and the significance threshold was defined with P value and FDR[15].

RIP assay

RNA immunoprecipitation experiments were performed using a Magna RIP RNA-Binding Protein Immunoprecipitation Kit (Millipore Corp, Billerica, MA, United States) according to the manufacturer’s instructions using enhancer of zeste homolog 2 (EZH2) antibodies (Cell Signaling Technology Inc., Danvers, MA, United States). The coprecipitated RNAs were detected via qRT-PCR. RIP assays were performed in biologic triplicates.

Statistical analysis

All statistical data were analyzed using SPSS 17.0 software (SPSS Inc., Chicago, IL, United States). Differences in lncRNA expression between the tumor and corresponding non-tumor tissues were analyzed using Student’s t-tests. Statistical significance was considered at P < 0.05.

RESULTS

Different expression profiles of lncRNAs and mRNAs in HBV-related HCC and corresponding non-tumor tissues

Five HBV-related patients (HBV surface antigen-positive) were selected to visualize the differentially expressed lncRNAs and mRNAs. Volcano plot analysis was used to directly detect lncRNAs and mRNAs abnormally expressed between the HCC tissues and corresponding non-tumor tissues (Figure 1). Hierarchical clustering is a simple and commonly used clustering technique to analyze gene expression data. Figure 2 shows the differentially expressed lncRNAs and mRNAs. Data analysis revealed 1772 differently expressed lncRNAs (Table S3) and 2508 differentially expressed mRNAs (Table S4) in HCC tissues compared with corresponding normal tissues. Of the 1772 differentially expressed lncRNAs, 637 were upregulated and 1135 were downregulated. Of the 2508 differentially expressed mRNAs, 1194 were upregulated and 1314 mRNAs were downregulated. The distinctively expressed lncRNAs in HCC are listed in Table 1, several of which, particularly MEG3 and GAS5, reportedly play roles in cancer development.

Table 1 Partly distinctive expressed long noncoding RNAs in hepatocellular carcinoma.

lncRNA	Fold change	Regulation	RNA length	Chromosome
ENST00000577848	12.43885	Up	1455	18
ENST00000426413	101.25230	Up	1577	7
NR_024478	14.71172	Up	2282	6
uc010kwq.2	12.91793	Down	4744	7
TCONS_00005259	10.69086	Down	3742	2
uc002nbr.3	29.97769	Down	1413	19
TCONS_00019684	19.73914	Down	1147	11
NR_027133	12.69278	Down	2597	15
NR_033957	17.82503	Down	2877	10
MEG3	10.96779	Down	1351	14

lncRNA: Long noncoding RNA.

Figure 1 Long noncoding RNAs and mRNA profile comparisons between the hepatocellular carcinoma and adjacent non-tumor tissues. Volcano plots used to distinguish the different expressed A: Long noncoding RNAs; and B: mRNAs. The vertical green lines correspond to 2.0-fold up- and downregulation, respectively, and the horizontal green lines represent P = 0.05. The red points represent the differentially expressed genes.

Figure 2 Hierarchical clustering of long noncoding RNAs and mRNAs in hepatocellular carcinoma. Hierarchical clustering analysis of A: 1772 Long noncoding RNAs; and B: 2508 mRNAs that were differentially expressed between hepatocellular carcinoma samples (A1-A5) and the adjacent non-tumor samples (B1-B5). The red and the green shades indicate the expression above and below the relative expression, respectively, across all samples.

qRT-PCR validation

To validate the microarray data, four upregulated lncRNAs were randomly selected and analyzed for their expression levels in tissue samples from 14 HBV-related HCC patients. The corresponding non-tumor tissues were analyzed via qRT-PCR (Figure 3), in which the obtained results are consistent with the microarray data.

Figure 3 Quantitative real-time reverse-transcription PCR. Long noncoding (lnc)RNAs GAS5, uc003ycp.3, uc004bdv.3, and ULK4P2 were differentially expressed between hepatitis B virus-related hepatocellular carcinoma tissues (ca) and the paired non-tumor samples (p).

lncRNA classification and subgroup analysis

According to the relationship between lncRNAs and their associated protein-coding genes, lncRNAs can be classified into natural antisense, intronic antisense, bidirectional, exon-sense overlapping, intron-sense overlapping, and intergenic. In order to systematically predict the function of lncRNAs, lncRNA subgroup analyses were performed, including antisense lncRNA, large intergenic noncoding (linc)RNA, Hox Loci lncRNA, T-UCR, and enhancer-like lncRNAs analyses based on those classifications. The profiles concerning enhancer-like lncRNAs (Table S5), HOX Loci lncRNAs (Table S6), and lincRNAs were examined (Table S7).

GO and signaling pathway analyses of differentially expressed mRNAs

Due to the lack of a comprehensive annotation database for lncRNAs, a profile of mRNAs in HCC was constructed to indirectly predict the function of the lncRNAs. A total of 1270 filtered mRNAs (threefold change) and 392 lncRNAs (fivefold change) were included in GO and Signaling pathway analyses (Figure 4). The GO analysis showed the distinctive functions of up- and downregulated involve the cell cycle, cell division, immune response, and DNA replication. The Pathway analysis showed that aberrantly genes are involved in the cell cycle, mitogen-activated protein kinase pathway, phosphoinositide 3-kinase-protein kinase B signaling pathway, and cytokine-cytokine receptor interaction.

Figure 4 Gene ontology and pathway analyses. A total of 1270 differentially expressed mRNAs were chosen with a fold change > 3; the column graphs represent the enrichment of these mRNAs. The (-lg^P) value is a positive correlation with gene ontology (GO) (A); Signaling pathway enrichment for upregulated mRNAs (B); GO (C); Signaling pathway enrichment for downregulated mRNAs (D). The (-lg^P) values above five are presented.

lncRNA-mRNA coexpression network

Coexpression network analysis was performed between the filtered 392 lncRNAs and the filtered 348 mRNAs with significant enrichment in GO and Signaling pathway analyses. The network structure of the HCC tissues and corresponding non-tumor tissues samples was markedly different (Figure 5).

Figure 5 Long noncoding RNA-mRNA coexpression network. A: The long noncoding (lnc)RNA-mRNA network containing the 349 filtered mRNAs and 392 filtered aberrant expressed lncRNAs in HCC tissues; B: The lncRNA-mRNA network containing the 349 filtered mRNAs and 392 filtered aberrant expressed lncRNAs in paired non-tumor tissues. Upregulated RNAs are shown in red, and downregulated RNAs are presented in blue. Nodes without a ring represent mRNAs, nodes with a ring represent lncRNAs, and node size represents the degree of centrality.

RIP assay

EZH2 is a component of the polycomb repressive complex 2 (PRC2). EZH2 has been reported in numerous cancers, including HCC[16]. It is reported that 20% lincRNAs can combine with EZH2, including HOTAIR, HEIH, and H19[17,18] which may regulate gene expression in cancers. Thus, we hypothesized that EZH2 may also have some relationship with lincRNAs that are differentially expressed in HBV-related HCC. Significantly expressed lincRNAs were selected to validate their combination with EZH2 using a RIP assay. The data showed that ULK4P2 is physically combined with EZH2 (Figure 5). Next, a ULK4P3-mRNA subnetwork was constructed to characterize the role of ULK4P2 (Figure 6). In the coexpression network, ULK4P2 is connected to six lncRNAs and nine mRNAs that are enriched for gene products involved in tumor cell proliferation and metastasis.

Figure 6 Large intergenic noncoding RNA ULK4P2 combines with enhancer of zeste homolog 2. A, B: ULK4P2 combined with enhancer of zeste homolog 2 (EZH2) was measured in hepG2 cells via an RNA immunoprecipitation assay using an anti-EZH2 antibody (input = total RNA, IgG = control); C: ULK4P2 subnetwork in the hepatitis B virus-related hepatocellular carcinoma coexpression network. Upregulated RNAs are shown in red, downregulated RNAs are presented in blue, and nodes without a ring represent mRNAs, whereas nodes with a ring represent long noncoding RNAs. Node size represents the degree centrality.

DISCUSSION

Recently, studies have shown that more genomic sequences are transcribed into lncRNAs than protein-coding RNAs, and lncRNAs are being characterized at a rapid pace. With an increase in the number of well-characterized cancer-associated lncRNAs, the study of lncRNAs in cancer is now generating new hypotheses about the biology of cancer cells. For HCC, deregulated expression of both protein-coding genes and microRNAs has been suggested to have considerable potential for predicting the prognosis of HCC patients[19]. However, the research concerning lncRNAs in HBV-related HCC is still in a preliminary stage.

In this study, the aberrant lncRNAs in five HBV-related HCC and corresponding non-tumor tissues were screened using lncRNA microarray V3.0, which is designed for the global profiling of human lncRNAs and protein-coding transcripts, and it is updated from the previous Microarray V2.0. The lncRNAs are carefully constructed using the most highly respected public transcriptome databases (Refseq, UCSC knowngenes, Genecode, etc.), as well as landmark publications. Each transcript is represented by a specific exon or splice junction probe that can accurately identify an individual transcript. Positive probes for housekeeping genes and negative probes are also printed onto the array for hybridization quality control. Data was extracted and normalized using the GeneSpring GX v11.5.1 software package (Agilent Technologies) to guarantee the accuracy, which showed that there are 637 upregulated lncRNAs and 1135 downregulated. Four distinctively expressed lncRNAs were randomly selected from the microarray data to validate their reliability in 14 HBV-related HCC samples; the result is consistent with microarray data. Meanwhile, the differentially expressed mRNAs were also examined.

Recently, some classes of lncRNAs, such as lincRNAs and HOX lncRNAs, have been identified with specific functions in human cells[20]. Thus the data are displayed as three clusters of lncRNAs.

lncRNAs with enhancer-like function were identified using GENCODE annotation of the human genes. Profiling data of all probes for lncRNAs with enhancer-like function revealed that many enhancer-like RNAs were bidirectional, lacked a polyA tail, and had very low copy numbers[21]. Depletion of this type of lncRNA leads to decreased expression of their neighboring protein-coding genes, as enhancer-like lncRNAs can combine with enhancer-binding proteins to regulate gene activation. Furthermore, detailed functional analysis of noncoding RNA adjacent to the Snai1 locus, the master regulator of hematopoiesis, using reporter assays demonstrated a role for this noncoding RNA in an RNA-dependent potentiation of gene expression[22]. These studies suggest a role of enhancer-like lncRNAs in positive regulation of protein-coding genes.

lincRNAs are another class of newly discovered lncRNAs with dysregulated expression in many tumors. lincRNAs have been involved in diverse biologic processes, including cell-cycle regulation, immune surveillance, and embryonic stem cells[23], but their mechanism is elusive. The results of this study show aberrant expression of 374 lincRNAs, which are nearby (distance < 300 bp) coding gene pairs. Recently, it was reported that 20% of lincRNAs are bound by PRC2, and that additional lincRNAs are bound by other chromatin-modifying complexes to regulate gene expression[24]. The mechanism concerning lincRNAs is intriguing, with recently reported novel functions that indicate they work in a more conventional way[25].

Hox genes, a highly conserved subgroup of the homeobox superfamily, have crucial roles in development, regulating numerous processes including apoptosis, differentiation, motility, and angiogenesis. The abnormal expression of HOX genes plays an important role in tumorigenesis, oncogenesis, or tumor suppression[26]. In mammals, 39 HOX transcription factors are clustered on four chromosomal loci, termed HOXA through HOXD[27]. Transcription of many lncRNAs has been observed in human HOX loci, described herein as HOX lncRNAs, which mainly participate in epigenetic regulation. The HOX cluster profiling reported in the present study includes 95 probes designed for the four HOX loci to target lncRNAs and coding transcripts. The results show that HOTAIR, the best-studied HOX lncRNA, is upregulated in cancer tissues compared to corresponding non-cancer tissues, which is consistent with a previous report[28].

To confirm the differential gene expression between HCC and corresponding non-tumor tissues, 348 mRNAs involved in function functions and signal pathways and deregulated lncRNAs were used to construct the lncRNA-mRNA coexpression network. The results show that the network structure of the HCC tissues and corresponding non-tumor tissues samples was markedly different, which indicates that the coexpression patterns of mRNA and lncRNAs between HCC and corresponding non-tumor tissues are different. GO analysis of the filtered differentially expressed mRNAs showed that they are enriched in cell cycle, DNA replication, immune response, and signal transduction pathways. The Signaling pathway analysis showed that cell cycle, cytokine-cytokine receptor interaction and chemokine and phosphoinositide 3-kinase-protein kinase B signaling pathways were correlated with HBV-related HCC.

EZH2 is the core enzymatic subunit of PRC2, and evidence shows that EZH2 is overexpressed in many cancer types and is critical for cancer cell proliferation, metastasis, and drug resistance[29]. To detect a possible combination between the candidate lincRNAs and EZH2, RIP was performed using an EZH2 antibody and nuclear extracts of hepG2 cells. The results demonstrate that the lincRNA ULK4P2 is physically combined with the EZH2 antibody, indicating that it may have function by combining with EZH2 in HCC. ULK4P2 is an lincRNA with little-known function located on chromosome 15. Because the coexpression network may correspond to biologic pathways, and many functions of protein-coding RNAs could be found in NCBI RefSeq[30], an ULK4P2-mRNA subnetwork was constructed, which indicated its relationship with some cell proliferation and metastasis genes. These data suggest that ULK4P2 may have an important role in HBV-related HCC initiation, growth, or metastasis.

Recent advances in microarray technologies and the rapid development have transformed the study of tumors into the genome era. As a result, the unstudied part of noncoding RNAs have been implicated as oncogenes and suppressor genes in tumors. The result in this paper gives us more information to further study lncRNAs in HBV-related HCC.

COMMENTS

Background

Hepatocellular carcinoma (HCC) is a common and deadly cancer, but the mechanism regarding hepatitis B virus (HBV)-related HCC is still not entirely understood. Long noncoding RNAs (lncRNAs) are a group of noncoding RNAs that are longer than 200 nucleotides. Recent study of lncRNAs demonstrates their involvement in human disease, including cancers.

Research frontiers

lncRNA is a hotspot in the field of RNA, and shows involvement in many human diseases.

Applications

The result of the lncRNA microarray provide us with many more lncRNAs study in greater depth; the result may be useful for HBV-related patients. For example, ULK4P2 may participate in HCC by combining with EZH2.

Terminology

Enhancer of zeste homolog 2 (EZH2) is the core enzymatic subunit of polycomb repressive complex 2, and evidence shows that EZH2 is overexpressed in many cancer types and is critical for cancer cell proliferation, metastasis, and drug resistance. ULK4P2 is a large intergenic noncoding RNA.

Peer-review

The manuscript is well presented and of interest. The design of study is appropriate, the study was done well, and its results can contribute to knowledge of this topic. All parts of the manuscript are well organized and valuable conclusions are provided. References are also appropriate, relevant and updated.