Abstract

Liver regeneration is a tissue growth process after loss or injury of liver tissue, which is a compensatory hyperplasia rather than true regeneration, mainly depending on hepatocyte proliferation. Currently, a large number of studies on hepatocyte proliferation have been conducted. However, studies on the regulation of long noncoding RNA (lncRNA) on hepatocyte proliferation are still limited. To identify specially expressed lncRNA during rat liver regeneration, high-throughput sequencing technology was performed, and a total of 2446 lncRNAs and 4091 mRNAs were identified as significantly differentially expressed. Gene ontology (GO) enrichment analysis was performed to analyze the role of differentially expressed mRNAs, and 695 mRNAs were identified to be related to cell proliferation. Then, an lncRNA-mRNA coexpression network based on the differentially expressed lncRNAs and proliferation-related genes was constructed to analyze the potential function of lncRNAs on hepatocyte proliferation, and ten lncRNAs, NONRATT003557.2, NONRATT005357.2, NONRATT003292.2, NONRATT001466.2, NONRATT003289.2, NONRATT001047.2, NONRATT005180.2, NONRATT004419.2, NONRATT005336.2, and NONRATT005335.2, were selected as key regulatory factors, which may play crucial roles in hepatocyte proliferation during rat liver regeneration. Finally, a protein-protein interaction (PPI) network was established to illuminate the interaction between proliferation-related genes, and ten hub genes (Aurkb, Cdk1, Cdc20, Bub1b, Mad2l1, Kif11, Prc1, Ccna2, Top2a, and Ccnb1) were screened with the MCC method in the PPI network, which may be important biomarkers involved in the hepatocyte proliferation during rat liver regeneration. These results may provide clues for a more comprehensive understanding of the molecular mechanism of hepatocyte proliferation during rat liver regeneration.

1. Introduction

The vast majority of the eukaryotic genomes are transcribed into noncoding RNAs, which can be divided into small noncoding RNAs (<200 bp) and long noncoding RNAs (lncRNAs; ≥200 nt) based on transcript size [1]. lncRNAs can be divided into five categories: sense, antisense, bidirectional, intronic, and intergenic [2]. Initially, lncRNAs were considered to be “dark matter,” as byproduct of transcription of RNA polymerase II, with no biological function. In the past 20 years, genome-wide identification of lncRNAs has become possible with the development of high-throughput technology of RNA-seq, many of which are involved in various biological functions [3]. Increasing lncRNAs have been found to play a critical role in biological processes, like development [4], gene transcriptional regulation [5], chromatin regulation [6], epithelial-to-mesenchymal transition (EMT) [7], and cell proliferation [8].

In rodents and humans, the liver can grow rapidly after partial hepatectomy (PH) or acute chemical injury. This growth process is known as LR, which is a compensatory hyperplasia rather than true regeneration [9]. During LR, quiescent hepatocytes undergo one or two rounds of replication and then return to a nonproliferative state [10]. This process is very complex and regulated by a variety of growth factors, cytokines and noncoding RNAs [11, 12]. Therefore, the study of the molecular mechanism of hepatocyte proliferation is crucially important to understand the process of LR and provide clues for the treatment of liver diseases. Several recent studies have shown that lncRNAs play a critical role in hepatocyte proliferation [1214]. However, the study of hepatocyte proliferation during LR is still largely unknown.

In the present study, high-throughput sequencing technology was used to identify DE lncRNAs and mRNAs during rat LR. Then, functional enrichment analysis of DE mRNAs was performed to screen proliferation-related genes. Finally, the lncRNA-mRNA coexpression network and PPI network were constructed based on DE lncRNAs and proliferation-related genes to elucidate the molecular mechanism of hepatocyte proliferation during LR. These results lay a foundation for understanding the regulatory function of lncRNAs on hepatocyte proliferation and provide an important clue for the study of the LR process.

2. Materials and Methods

2.1. Preparation of 2/3 Hepatectomy Model

The healthy adult male Sprague Dawley (SD) rats weighing 210∼250 g were provided by the Laboratory Animal Center of Zhengzhou University (Zhengzhou, China). These rats were raised in a controlled temperature room of 19∼23°C with a relative humidity of 50∼70% and an illumination time of 12 h/d (8 : 00 to 20 : 00) and permitted to freely have water and food. A total of 60 rats were taken for the experiment with six rats per group: nine PH groups and one normal group (CG). The rats in PH groups were conducted 2/3 PH according to the method of Xu et al. The rats were anesthetized and condemned to death at 0, 2, 6, 12, 24, 30, 36, 72, 120, and 168 h after operation. The right liver lobes of six rats were mixed at each time point and stored at −80°C. All operations conformed to the Animal Protection Law of China and Animal Ethics.

2.2. RNA Sequencing

RNA sequencing was performed by the Shanghai OE Biotech (Shanghai, China). In brief, the mirVana miRNA Isolation Kit (Ambion) was used to extract the total RNA from liver tissues. The TruSeq Stranded Total RNA with Ribo-Zero Gold (Illumina) was used to construct cDNA libraries. The purified cDNA libraries were sequenced on Illumina HiSeq 2500 following the manufacturer’s instruction. After filtrating the adaptor and low-quality reads, clean reads were obtained for subsequent analysis. The reads were matched to the rat reference genome using the hisat2 (v2.2.1.0) software. The StringTie2 (v1.3.3b) software was used to splice the aligned reads. lncRNA identification included two categories: one is known lncRNA, which completely matches with the known lncRNAs, and the other is candidate lncRNA lacking protein-coding ability, whose length is greater than 200 bp and exon is greater than or equal to 2. The software CPC (v0.9-r2), CNCI (v1.0), Pfam (v30), and PLEK (v1.2) were used to predict the protein-coding ability of transcripts. The expression of transcription was calculated by the fragments per kilobase of exon per million reads mapped (FPKM) method using the bowtie2 (v2.2.9) and eXpress (v1.5.1) software.

2.3. Identification of Differentially Expressed lncRNAs

The counts of lncRNAs in each sample were standardized by the baseMean value using the DESeq (1.18.0) software. Differentially expressed (DE) lncRNAs were identified with fold change ≥2 or ≤0.5 and as the threshold. All DE lncRNAs in nine PH groups underwent hierarchical clustering analysis using the cluster3.0 and treeview software.

2.4. GO Enrichment Analysis

Gene ontology (GO) enrichment analysis of the DE mRNAs was conducted using David Bioinformatics Resources 6.8 (https://david.ncifcrf.gov/). The enrichment analysis consisted of three parts: biological process (BP), molecular function (MF), and cellular component (CC). is considered statistically significant, which was calculated by the EASE score.

2.5. lncRNA-mRNA Coexpression Analysis

To explore the relationship between DE lncRNAs and proliferation-related genes, a coexpression analysis was performed. Pearson’s correlation coefficients (PCCs) were calculated between the DE lncRNAs and the proliferation-related genes, and only lncRNA-mRNA pairs with PCC ≥0.8 and were selected and considered as coexpression. Then, these lncRNA-mRNA pairs were used to construct a coexpression network, which was visualized by the Cytoscape v3.6.1 software. The node degree was determined by the number of directly connected neighbors to the topological property of the network.

2.6. Construction of PPI Network and Screening of the Key Gene

To illustrate interactions between proliferation-related DE mRNAs, the string database (https://string-db.org/) was used to construct a protein-protein interaction (PPI) network. Only the interacting pairs with combined score ≥0.4 were selected and considered to be significant. The PPI network was visualized by using the Cytoscape v3.6.1 software. Then, a Cytoscape plugin cytoHubba was used to identify the key genes adopting the MCC method.

2.7. RT-PCR Validation

RNA-seq results were validated by RT-PCR, and the primers are listed in Table 1. cDNA was synthesized using the cDNA Reverse Transcription Kit (Takara, Tokyo, Japan). The qRT-PCR was performed using Q-SYBR Green Supermix (Bio-Rad). Primers were also designed to amplify β-actin as an endogenous control. The expression of each lncRNA was represented as fold change using 2−ΔΔCt methods.

3. Results

3.1. Sequencing and Identification of lncRNAs during Rat LR

To identify the expression of lncRNAs during rat LR, 10 cDNA libraries were constructed from the regeneration rat liver at different time points after surgery. The Illumina HiSeq X Ten platform was used to sequence these cDNA libraries, and a total of 1116M raw reads were produced. After filtering adaptor sequences and low-quality reads, 1084.09M clean reads were obtained. The percentage of clean reads varied from 90.73 to 95.29%, and the percentage of GC content varied from 48.77 to 51.08% in each library. Approximately 95.91–97.83% of clean reads were selected for further research after mapping the clean reads to the rat reference genome (Table 2).

3.2. Identification of DE lncRNAs and DE mRNAs

The expression abundance of the lncRNAs was evaluated by FPKM (fragments per kB per million reads) using DESeq. A total of 2446 lncRNAs were determined to be differentially expressed during rat LR, with 1120 upregulated, 731 downregulated, and 595 up/downregulated lncRNAs (Figure 1(a) and Table S1). To explore the similarity of gene expression, hierarchical clustering was adopted to analyze the expression of DE lncRNAs (Figure 2(a)). To further explore the interactions of DE lncRNAs at different stages: initial stage (2–6 h), proliferation stage (12–72 h), and termination stage (120–168 h), a Venn diagram was constructed using these DE lncRNAs (Figure 2(c)). Among them, 272 DE lncRNAs were common to all three stages.

Through high-throughput RNA-seq, the expression profile of 28635 mRNAs was measured. Among them, 4091 mRNAs were found to be differentially expressed, of which 2,256 were upregulated, 1,686 were downregulated, and 149 were up/downregulated (Figure 1(b) and Table S2). Hierarchical clustering was employed to analyze the expression similarity of DE mRNAs (Figure 2(b)). Venn analysis was conducted to explore the differences of DE mRNAs at different stages (Figure 2(d)).

3.3. GO and KEGG Enrichment Analysis

GO enrichment analysis was employed to determine the function of deregulated genes during rat LR, and 689 significant GO terms (451 under BP, 116 under CC, and 122 under MF) were enriched. The top 30 GO terms of the three groups are listed in Figure 3(a) (Table S3). The most enriched BP terms were response to drug, cell division, and chromosome segregation. As for CC, the most enriched terms were cytoplasm, nucleus, and nucleoplasm. The most enriched MF terms were related to binding activity. Of these BP terms, 41 were associated with cell proliferation involving 695 mRNAs, corresponding to 585 genes (Figure 3(b) and Table S3).

3.4. Coexpression Network Construction Based on DE lncRNAs and Proliferation-Related DE mRNAs

The coexpression analysis of the screened proliferation-related DE mRNAs and DE lncRNAs found 18343 significant coexpression pairs, of which 17547 pairs were positively correlated (PCC ≥ 0.8) and 796 pairs were negatively correlated (PCC ≤ −0.8) (Table S4). Subsequently, a coexpressed network was constructed using screened lncRNA-mRNA pairs, and it was found that some lncRNAs could interact with multiple mRNAs. According to nodes and connectivity, the top 10 lncRNAs were selected, which may exert important regulatory roles in hepatocyte proliferation, lncRNAs NONRATT003557.2 (degree = 178), NONRATT005357.2 (degree = 178), NONRATT003292.2 (degree = 177), NONRATT001466.2 (degree = 176), NONRATT003289.2 (degree = 176), NONRATT001047.2 (degree = 175), NONRATT005180.2 (degree = 175), NONRATT004419.2 (degree = 173), NONRATT005336.2 (degree = 173), and NONRATT005335.2 (degree = 169) (Figure 4). The expression levels of the top 10 lncRNAs during rat LR are shown in Table 3. Among them, NONRATT005357.2 had the highest difference at 24 h in regenerating the liver, and the fold change was 119.

3.5. Construction of PPI Network Based on Proliferation-Related Genes

To better elucidate the interaction network of proliferation-related genes, a PPI network containing 551 nodes with scores greater than or equal to 0.4 was constructed using the string database (Figure 5(a)). Moreover, a cytoHubba plugin was used to select the top 10 key genes from the PPI network using the MCC method. The 10 key genes were Aurkb, Cdk1, Cdc20, Bub1b, Mad2l1, Kif11, Prc1, Ccna2, Top2a, and Ccnb1 (Figure 5(b)). The expression levels of the 10 key genes during rat LR are shown in Table 4. The median of multiple mRNAs corresponding to one gene is taken as its expression level. Among them, Cdk1 had the highest difference at 24 h in regenerating the liver, and the fold change was 109.

3.6. Quantitative Real-Time PCR Validation

We validated the high-throughput RNA-seq results by performing qRT-PCR analysis of differentially expressed lncRNAs; the expression patterns revealed similar conclusions (Figure 6). Our RNA-seq results showed lncRNAs NONRATT003289.2, NONRATT001466.2, NONRATT004419.2, and NONRATT005336.2 were upregulated during rat LR.

4. Discussion

In rodents and humans, the liver can grow rapidly after partial hepatectomy or acute chemical injury. This growth process is known as LR, which is a compensatory hyperplasia rather than true regeneration, mainly depending on the proliferation of hepatocytes [9]. To explore the regulatory roles of lncRNAs in hepatocyte proliferation during rat LR, high-throughput RNA-seq was performed to identify lncRNAs and mRNAs. In the present research, 2446 DE lncRNAs and 4091 DE mRNAs were identified. To investigate the function of DE mRNAs during rat LR, GO enrichment analysis was performed. The result indicated that a large number of GO terms were associated with response to stress, cell proliferation, oxidation-reduction, regulation of transcription, metabolism, and apoptosis, which were considered to be important activities in LR [1517]. Of these GO terms, 41 were associated with cell proliferation, involving 695 mRNAs.

Through coexpression network analysis, a total of 18343 coexpression pairs were obtained based on DE lncRNAs and proliferation-related DE mRNAs. Among these coexpression pairs, 95.7% were positively correlated and 4.3% were negatively correlated. These results indicate that lncRNAs may regulate the expression of genes mainly in a positive way. It was also found that some lncRNAs could be coexpressed with multiple mRNAs, suggesting that lncRNAs could regulate multiple mRNAs. According to nodes and connections, the top 10 lncRNAs (NONRATT003557.2, NONRATT005357.2, NONRATT003292.2, NONRATT001466.2, NONRATT003289.2, NONRATT001047.2, NONRATT005180.2, NONRATT004419.2, NONRATT005336.2, and NONRATT005335.2) were selected. These lncRNAs were connected to 158 common DE mRNAs, respectively, such as Slc25a16, Cdk1, Rfc3, Brca2, Mcm7, Cdca7, Hat1, Chaf1b, Rfc4, and Mcm4. These findings indicated that these lncRNAs and their interaction genes may play an important role in hepatocyte proliferation during rat LR.

cytoHubba, a Cytoscape plugin, could provide 11 topological analysis methods to explore important nodes in biological networks including maximal clique centrality (MCC). MCC, a new method, has a better precision in identifying hub proteins. Through the MCC method, 10 hub mRNAs were selected including Aurkb, Cdk1, Cdc20, Bub1b, Mad2l1, Kif11, Prc1, Ccna2, Top2a, and Ccnb1. Aurkb (aurora kinase B) was involved in regulation of chromatin, glycolysis metabolism, regulation of telomerase activity, regulation of zygote development, and turnover of kinetochore microtubules [1822]. It was also associated with AKT signaling pathway, MAPK signaling pathway, inflammatory pathway, DNA damage response pathway, Wnt signaling pathway, and PI3K/Akt/NF-κB signaling pathway [21, 2327]. Cdk1 (cyclin-dependent kinase 1) was involved in cell cycle, mitosis, DNA end resection, homologous recombination, DNA damage checkpoint, anaphase spindle dynamics, cytokinesis, chromosome biorientation, apoptosis, cell proliferation and survival, and DNA synthesis [2835]. It also played an important role in AMPK signaling pathway, AKT signaling pathway, Golgi checkpoint signaling, PDK1-PI3K/Akt signaling pathway, type I interferon signaling, and p53 signaling pathway [3640]. Cdc20 (cell division cycle 20), an APC activator protein, could regulate mitosis, cell cycle, and presynaptic differentiation [4144]. It was also involved in Wnt/β-catenin signaling pathway and p38 MAPK signaling pathway. Bub1b (BUB1 mitotic checkpoint serine/threonine kinase B) and Mad2l1 (mitotic arrest deficient 2 like 1) were involved in the spindle checkpoint during mitosis [4548]. Kif11 (kinesin family member 11), a microtubule motor, played a vital role in regulating the transport of β-actin mRNA and cell motility through physically interacting with ZBP1, which could govern the direction of migration by responding to directional cues in chemotaxis [49, 50]. Prc1 (protein regulator of cytokinesis 1), a cell cycle protein, played an important role in ensuring proper cell division through directly acting on protein (FIP), and P53 could regulate the transcription of Prc1 [51, 52]. Ccna2 (cyclin A2) was a key factor in cell cycle, and Ccna2 repression regulated by miR-22 could inhibit HCC cell proliferation and tumorigenesis [53, 54]. It was also a prognostic biomarker for several cancers, such as ER+ breast cancer, pancreatic ductal adenocarcinoma, and colorectal cancer [5557]. Top2a (DNA topoisomerase II alpha), a key enzyme in DNA replication, was recruited to ultrafine anaphase bridges (UFBs) by TopBP1 to ensure faithful separation of sister chromatids, and these proteins played an important role in maintaining genome stability [58]. Top2a was overexpressed in hepatocellular carcinoma and associated with early age onset, shorter patient survival, and chemoresistance [59]. Ccnb1 (cyclin B1) was a key factor in cell cycle and regulated by STAT3 via the E2F modulating G2-M phase checkpoint [60]. High expression of Ccnb1 was associated with poor prognosis in HCC patients, and knockdown of Ccnb1 could significantly inhibit cell proliferation, migration, and invasion in HCC [61]. Therefore, these genes may play a vital role in hepatocyte proliferation during LR.

5. Conclusions

In this study, the comprehensive expression abundance of lncRNAs and mRNAs was identified by RNA-seq analysis during rat LR. The lncRNA-mRNA coexpression network and PPI network based on lncRNAs and proliferation-related genes were constructed, and 10 key lncRNAs and 10 key mRNAs were determined that may play crucial roles in hepatocyte proliferation. Our study provides a new idea to better understand the mechanism of LR.

Data Availability

The data used to support the findings of this study are included within the supplementary information file(s).

Conflicts of Interest

The authors declare that there are no conflicts of interest.

Acknowledgments

This study was financially supported by the National Natural Science Foundation of China (no. 31572270) and the National Natural Science Foundation of China Youth Science Fund Project (no. 31601038).

Supplementary Materials

Table S1: results of analyzing differentially expressed lncRNAs at 2, 6, 12, 24, 30, 36, 72,120, and 168 h after PH was compared with CG. Table S2: results of analyzing differentially expressed mRNAs at 2, 6, 12, 24, 30, 36, 72,120, and 168 h after PH was compared with CG. Table S3: results of GO enrichment analysis of differentially expressed mRNAs. Table S4: screened coexpression pairs between DE lncRNAs and proliferation-related DE mRNAs based on the Pearson correlation coefficient. (Supplementary Materials)