RNA-Seq transcriptome profiling in three liver regeneration models in rats: comparative analysis of partial hepatectomy, ALLPS, and PVL

The liver is a unique organ that has a phenomenal capacity to regenerate after injury. Different surgical procedures, including partial hepatectomy (PH), intraoperative portal vein ligation (PVL), and associated liver partition and portal vein ligation for staged hepatectomy (ALPPS) show clinically distinct recovery patterns and regeneration. The observable clinical differences likely mirror some underlying variations in the patterns of gene activation and regeneration pathways. In this study, we provided a comprehensive comparative transcriptomic analysis of gene regulation in regenerating rat livers temporally spaced at 24 h and 96 h after PH, PVL, and ALPPS. The time-dependent factors appear to be the most important determinant of post-injury alterations of gene expression in liver regeneration. Gene expression profile after ALPPS showed more similar expression pattern to the PH than the PVL at the early phase of the regeneration. Early transcriptomic changes and predicted upstream regulators that were found in all three procedures included cell cycle associated genes (E2F1, CCND1, FOXM1, TP53, and RB1), transcription factors (Myc, E2F1, TBX2, FOXM1), DNA replication regulators (CDKN1A, EZH2, RRM2), G1/S-transition regulators (CCNB1, CCND1, RABL6), cytokines and growth factors (CSF2, IL-6, TNF, HGF, VEGF, and EGF), ATM and p53 signaling pathways. The functional pathway, upstream, and network analyses revealed both unique and overlapping molecular mechanisms and pathways for each surgical procedure. Identification of molecular signatures and regenerative signaling pathways for each surgical procedure further our understanding of key regulators of liver regeneration as well as patient populations that are likely to benefit from each procedure.

While adult hepatocytes are normally quiescent, they show phenomenal replicative potential when parenchymal loss occurs. This has made the liver an excellent model for studying organ regeneration [1][2][3][4] . The regenerative process is typically divided into the priming, proliferating and termination phases, occurring approximately at first, third and seventh days after resection, respectively 5 . Clinically, the regenerative ability of livers provides the cornerstone of many surgical treatments for liver diseases; resection of diseased portions seems to improve survival 6 . However, extensive removal of liver parenchyma, which may be necessary in some cases, may overcome the functional reserve of the remaining hepatocyte population, leading to decompensation, failure and even mortality 7 . To avoid this complication, pre-operative iatrogenic induction of hepatocyte proliferation can be attempted to increase the future remnant tissue's resilience against the anticipated surgical volume loss 7 .
Several options are available to induce the regenerative capacity of the liver. While cellular proliferation may be chemically induced 8 , surgical induction procedures provide additional advantages such as intraoperative Experimental design and animal groups. Rats were divided into groups and anesthetized under isoflurane administration. A middle line incision was performed on the abdomen. (1) Sham group: Abdomen was closed after manipulation of the liver hilum; (2) PVL group: Portal vein ligation was performed on all branches except the branch to the right median lobe using a size 7-0 nylon thread; and (3) ALPPS group: selective portal vein ligation and liver parenchymal partitioning were performed (4) PH group: A 5-0 nylon thread was tied tightly around the median and left lateral lobes (~70% of the total liver mass) and both lobes were then resected, as described previously 9 . There were three rats for each group at each time point (24 h and 96 h). The detail of surgical operations and experimental grouping is provided in our previous study 9 .
RNA sequencing analysis. TRIzol (Invitrogen) was used for isolation of total RNA (15 µg) that was utilized for purifying the poly(A)-containing mRNA molecules, RNA amplification, and synthesis of double-stranded cDNAs according to Illumina's TruSeq RNA Sample Prep guidelines (San Diego, CA). Multiplexed samples were sequenced at 43 bp length on Illumina-based Technology. Triplicate biological replicates were performed for each group. The paired-end reads of each sample were aligned to the rat genome (rn5) using the TopHat 18 . Transcript abundance was estimated by Cufflinks 19 . The quantification and normalization (DESeq method 20 ) and further downstream analyses of identification of differentially expressed genes (DEGs) were done by using Strand NGS 2.7 (Strand Life Sciences, India) and PARTEK Genomics Suite (Partek Inc., St. Louis, MO, USA). Significantly regulated genes across different operation types (ALPPS, PVL, PH, and Sham) and two time points (24 h and 96 h) were determined using two-factor analysis of variance (ANOVA) by taking operation type, time and their interactions into the statistical model. Genes exhibiting false discovery rate (FDR)-adjusted P value<0.05 and the absolute fold changes (FC) > 2 were considered significant. Functional pathway, upstream regulator, and network analyses. Functional, pathway, and gene ontology (GO) enrichment analysis were performed using Database for Annotation, Visualization and Integrated Discovery (DAVID) 21 , Protein Analysis Through Evolutionary Relationships (PANTHER ™ ) classification systems and Ingenuity Pathways Analysis (IPA) (QIAGEN Inc., https://www.qiagenbioinformatics.com/products/ ingenuity-pathway-analysis). We also performed upstream regulator, canonical pathways, and gene network analyses using IPA and Network Analyst 22 . The DEGs lists for each surgical procedure for different time points were mapped to its corresponding gene object in the Ingenuity pathway knowledge base and protein-protein interaction networks. A right-tailed Fisher's exact test was used to calculate a p-value determining the probability that the biological function (or pathway) assigned to that data set is explained by chance alone. The IPA upstream regulator analysis predicts the upstream transcriptional regulators based on the Ingenuity ® Knowledge Base by examining how many known targets of the upstream regulators are present in the differentially expressed gene list. An overlap p-value, based on significant overlap between genes in the list and known targets regulated by the transcriptional regulator, and an activation z-score are computed. The predicted activation state and activation z-score are based on the direction of fold change values that we observed in the gene expression data. The Inter-procedural similarities/differences at 24 h and 96 h post operations. Since the time had the greatest effect on the expression, we next compared the transcriptomes of the procedures at each time point separately. Differentially expressed genes (DEGs) that are specific or commonly dysregulated between the procedures at 24 h and 96 h were obtained using the Venn diagram approach (Fig. 1a). The PCA and hierarchical clustering analyses of the significant genes (FDR < 5% and FC > 2) separated the samples according to the surgical types at each time points (Supplementary Fig. S1).
At 24 h, the PH vs Sham produced the largest number of DEGs, followed by ALPPS vs Sham (Fig. 1a, Supplementary Table S2). There were 1242 DEGs for the PH (703↑ up-regulated, 539↓ down-regulated), 814 for the ALPPS group (536↑, 278↓), and 539 for the PVL (428↑, 111↓), and 411 DEGs (351↑, 60↓) were common to all surgical procedures. The arrow represents an increase (↑) or decrease (↓) in fold compared to Sham. Comparing different procedures, 123, 542, and 71 genes were specifically expressed in ALPPS, PH, and PVL, respectively (Fig. 1a). There were more DEGs commonly shared between the ALPPS and the PH than with the PVL, especially at 24 h. The list of DEGs shared and unique to each procedure is given in Supplementary Table S2. The ALPPS showed more similar expression pattern to the PH than the PVL especially at the early phase of the regeneration at 24 h post operation ( Supplementary Fig. S1). Two-dimensional hierarchical clustering of genes, that are significant in at least one surgical procedure (with respect to sham group) at any of the two time points, and samples. The figure shows the most associated GO biological processes for each cluster of genes. Red and green denote highly and weakly expressed genes, respectively.
We further performed an extensive comparison of functional enrichment and gene networks for the three surgical procedures using different bioinformatics tools to further investigate molecular similarities/differences among the procedures. Gene ontology enrichment analysis revealed that genes related to cell cycle, mitotic cell cycle, M phase, and DNA replication and repair were significantly enriched among the up-regulated genes at 24 h for all surgical procedures ( Supplementary Fig. S1). However, genes related to oxidation-reduction, triglyceride, and steroid metabolic processes were the most critically enriched among the down-regulated genes in ALPPS and PH only (Fig. 2a, Table 1, and Supplementary Fig. S4).
At 96 h post-operatively, there were 657 (493↑, 164↓), 448 (387↑, 61↓), and 211 (179↑, 32↓) genes for ALPPS, PH and PVL, respectively, when compared with the Sham livers (FDR < 5% and FC > 2) (Fig. 1a, Supplementary Table S3). The ALPPS vs Sham produced the largest number of DEGs. Genes responsible for cell cycle, ion transport, mitosis, and organelle fission seemed to be highly enriched in the up-regulated genes in the ALPPS; genes responsible for cell activation, response to wounding, and immune response seemed to be the most significantly enriched in the PH. The down-regulated genes at 96 h, especially for the ALPPS and PH, were enriched for oxidation-reduction, metabolic process and inflammatory response (Fig. 2b). The cell adhesion, immune system process, cellular process and cell proliferation associated genes were significantly enriched after ALPPS and PH, and macrophage activation, system development, and cellular component and morphogenesis were significant after all the procedures at 96 h (Table 1).  Fig. S2A, Tables S2 and S3, respectively). GO enrichment and functional analyses of DEGs at 24 h indicated involvement of genes related to cell cycle, DNA replication and repair, cell division, mitosis, and DNA metabolic process ( Table 1 At 96 h, the DEGs that are commonly regulated in all surgical groups were mostly up-regulated (94%) and mainly involved in cellular movement, system development, extracellular matrix organization, and immune response ( Supplementary Fig. S2B). The gene ontology and network analyses indicated enrichment of genes related to cell death and survival, cellular growth and proliferation, lipid metabolism, and immune cell trafficking, including genes such as LGALS3, FCGR2A, SPARC, Integrin, and Collagens (p-value <0.01) (Fig. 3d-f), at 96 h post-operations.
Functional, pathway and gene network comparison of temporal changes after ALPPS, PH and PVL. We further performed an extensive comparison of functional enrichment and gene networks for the three surgical procedures using different bioinformatics tools to further investigate molecular similarities/differences of the procedures as well as to investigate temporal changes within each surgical procedure at two time points.
Genes related to cell cycle, mitotic cell cycle, M phase, and DNA replication and repair were significantly enriched among the up-regulated genes at 24 h for all surgical procedures (Fig. 2a, Table 2, and Supplementary Bar chart of the most significant GO biological processes that are associated with up-and downregulated genes for each surgical procedure at 24h (a) and 96h (b). X-axis represents the statistical significance of the enrichment (-log10(p-value)). Color-coding represents different surgical types.
www.nature.com/scientificreports www.nature.com/scientificreports/ Fig. S4). However, genes related to oxidation-reduction, triglyceride, and steroid metabolic processes were the most critically enriched among the down-regulated genes in ALPPS and PH only. On the other hand, at 96 h, genes responsible for cell cycle and organelle fission seemed to be highly enriched in the up-regulated genes in the ALPPS; genes responsible for cell activation, response to wounding, and immune response seemed to be the most significantly enriched in the PH. The down-regulated genes at 96 h, especially for the ALPPS and PH, were enriched for oxidation-reduction, metabolic process and inflammatory response (Fig. 2b).
The DEGs for each procedure at each time point were mapped to gene interaction networks in order to obtain deeper insights into the interactions of these genes among various pathways. The functional, pathway and gene network analyses highlighted the potentially critical genes, biological processes, and signaling pathways for the temporal changes within ALLPS, PVL, and PH groups ( Table 2, Supplementary Fig. S4). Cell cycle, mitosis, and DNA replication, process-related genes, such as Cdk1, E2f1, Ezh2, Ccnb1, Cdkn1a, Myc, and p53 pathways were significantly regulated at 24 h in all procedures. The networks for ALLPS and PH were more similar than PVL at early phase of the regeneration, corroborating with the clustering results in Fig. 1 and Supplementary Fig. S1. At 96 h, cell proliferation, cell cycle, mitosis and cell division remained active in the ALPPS when compared to the PH and PVL. Genes related to macrophage activation, anion transport, immune response, lipid metabolic process were also significantly regulated in the ALPPS (Table 1, Supplementary Fig. S4). Significantly altered pathways included the FXR/RXR activation, cell cycle and integrin signaling pathways in the ALPPS group at 96 h. On the other hand, cellular movement, immune cell trafficking, and signal transduction related genes and T-cell activation and Ras pathways were significantly altered in the PH group at 96 h. Top five significant molecular and cellular functions, canonical pathways, and predicted upstream regulators (activated/inhibited) for each surgical procedure at two time points are given in Table 2. Furthermore, the complete list of significant PANTHER and canonical pathways are given in Table 3 and Supplementary Table S4, respectively.

Temporal pattern of upstream transcriptional regulators after ALPPS, PH and PVL.
We identified the transcriptional regulators likely to be involved in the regenerating liver after each surgical procedure at two time points (Fig. 4). The IPA upstream regulator analytic identifies the upstream transcriptional regulators and mechanistic networks that can explain the gene expression changes observed after each procedure. There were 184, 158, and 206 upstream regulators predicted to be activated or inhibited after ALPPS, PH, and PVL, respectively at 24 hr (Fig. 4a). There was a noteworthy overlap of transcriptional regulators between different procedures at 24 h. The inter-procedural variations in ALPPS, PVL and PH are displayed at 24 h and 96 h (in Fig. 4b  www.nature.com/scientificreports www.nature.com/scientificreports/ respectively) for the top 15 activated upstream regulators. The union of top 15 activated upstream regulators in ALPPS, PVL and PH at 24 h and 96 h is used for the heatmap. Similarly, the temporal variations in each procedure are also investigated for ALPPS, PH and PVL, in Fig. 4d-f, respectively. Figure 4 also displays the heatmap of the top 50 activated/inhibited upstream regulators across all samples to investigate the temporal changes for the three procedures (complete list is given in Supplementary Table S5) (Fig. 4g).

Discussion
We performed next-generation RNAseq approach to comprehensively delineate liver regeneration induced by ALPPS, PVL, and PH approaches at two time points (24 h and 96 h). Recently, we have replicated ALPPS in a rat model and presented its quick rejuvenation mode 17 . In this study, we identified pertinent common transcriptomic signatures at critical time points (i.e. early proliferation and late-proliferation phase) of liver regeneration and revealed the inter-procedural and temporal variations in gene expression patterns.
In response to injury, liver regeneration is achieved by the activation of otherwise functional, fully-differentiated hepatocytes as a result of the autocrine and paracrine signaling (e.g., cytokines and growth  www.nature.com/scientificreports www.nature.com/scientificreports/ factors) 12,17 . Thus, major underlying changes in gene networks are expected to be seen at the cellular level. Depending on the initial injurious stimulus (or surgical procedure), such networks may have distinct elements as well as shared factors that are critical for regeneration. Therefore, delineation of the gene networks in liver regeneration is of particular interest. Moreover, understanding of signaling cascades and gene networks of the ALPPS, PVL and PH at 24 hours is vital as these may be critical factors for the initiation of regeneration or resulting in poor recovery, incomplete regeneration, and ultimately, liver failure.
Our data revealed that time-dependent factors appear to be the major source of variation in post-injury alterations of gene expression in liver regeneration. Early transcriptomic changes and the upstream regulators after all the procedures included cell cycle associated genes (E2F1, CCND1, FOXM1, TP53, and RB1), transcription factors (Myc, E2F1, TBX2, FOXM1) [23][24][25][26][27][28] , DNA replication regulators (CDKN1A, EZH2, RRM2) 29 , G1/S-transition regulators (CCNB1, CCND1, RABL6) 30 , cytokines and growth factors (CSF2, IL-6, TNF, HGF, VEGF, and EGF) 9,16,31,32 . At the cellular level, this corresponds to the transition from the quiescent G0 phase to active mitosis. Liver regeneration is governed by numerous growth factors, including HGF and EGF that are responsible priming the parenchymal and nonparenchymal liver cells and boosting their access into the cell cycle to proliferate to restore the original liver size after the surgical procedures 9,16,31,32 . Growth factors and signaling pathways activates cyclin-dependent kinases (CDKs), and upregulates the expression of CCND1 that encodes cyclin D1. CDK1 is essential for DNA replication and downstream formation of replication-initiation complexes in hepatocytes and shown to play a critical role in DNA replication control during rat liver regeneration following PH 33 . It is also required for the activity of CCNB1, a protein from cyclin regulatory proteins family that is essential for cell cycle control during G2/M (mitosis) transition.
Upstream events that can induce the quiescence-to-mitosis transition may potentially be important in liver regeneration, such as activation of E2F1, CCND1, FOXM1, and inhibition of TP53, and RB1. Notably, enhancer of zeste homolog 2 (EZH2) was among the early significantly up-regulated genes. As the functional subunit of Polycomb Repressive Complex 2 (PRC2), which normally regulates development and differentiation in healthy embryonic tissue, it is potentially involved in early de-differentiation that hepatocytes undergo to become highly proliferative 34 . Additionally, its up-regulation is consistent with the decreased lipid metabolism identified by our functional network analysis 35 . As a specialized function of hepatocytes, the decrease in lipid metabolism may well be an indication of dedifferentiation. Besides the critical inducers of dedifferentiation and G0-to-G1 transition, other molecules seem to be critical in sustaining a robust intracellular environment. One example is p400 E1A-associated protein (EP400), which was predicted to be an early upstream regulator. EP400 is essential for the control of reactive oxygen species (ROS) intracellularly, maintaining an oxidative-stress free environment without which DNA damage, senescence, and apoptosis may ensue 36 . Another example is Lipocalin 2 (LCN2), an innate-immunity molecule with iron-sequestering properties. It has a wide expression in various tissues but has been particularly used as an early biomarker for kidney injury 37 . Recently, the increased expression level of LCN2 has been demonstrated in acute liver injury as well [38][39][40] . Indeed, protein and mRNA levels of Lcn2 is significantly increased after partial hepatectomy 41 .
Immune mediators of liver regeneration that showed a significant increase in expression immediately after the operations were IL-6 and VEGF 17,[42][43][44] . These molecules are known mediators of acute and chronic inflammatory response, carrying out essential functions in repairing tissue injury in all parts of the human body. IL-6, secreted by macrophages, is a known pyrogen and a potent stimulator of the acute phase reactants following infection or trauma. VEGF is secreted by endothelial cells to promote vasculogenesis and angiogenesis. In liver injury, VEGF mediates the proliferation and mobilization of liver sinusoidal endothelial cell progenitor cells form the bone marrow to allow for the formation of new sinusoids in the regenerating, highly-vascular liver tissue 45 .
The network analysis of up-regulated genes at 96 h post-operatively revealed an expression profile predominantly in mediating tissue reconstruction, including cellular movement, system development, extracellular matrix organization and restoring specialized hepatocytic functions (e.g., carbohydrate transport and morphogenesis). Of note, this phase seems to be orchestrated by the macrophages (known as Kuppfer cells in the liver). Transforming growth factor-β (TGF-β), a potent inhibitor of inflammation and a stimulator of healing and tissue repair that is secreted by macrophages, was significantly up-regulated at 96 h. Furthermore, the analysis also highlighted genes such as LGALS3, FCGR2A, SPARC, integrin, and collagens, that may play significant roles in tissue regeneration and repair [46][47][48] . These genes closely coordinate their function with TGF-β in tissue organization, matrix structure and cell-to-cell interactions. The knockout mice of secreted protein acidic and rich in cysteine (SPARC) showed a reduction of expression of TGF-β1 and collagen in hepatic tissue 49 .
Among the three examined procedures, as noted in the clustering and upstream regular analyses, ALPPS and PH shared many significantly regulated genes whose expression were not otherwise significantly changed in PVL, especially at the early phase of regeneration. One possible explanation lies in the differences of the type of injury caused by each procedure. Specifically, ALPPS and PH both require the removal of hepatic parenchyma. On the other hand, the main stimulant of regeneration in PVL is ischemia and oxygen deprivation. Clinically, PVL is the  www.nature.com/scientificreports www.nature.com/scientificreports/ least invasive and has the lowest rate of morbidity and mortality but also has a slower rate of regeneration than the other procedures 50 . Compared to PVL, the upstream regulator analysis revealed the activation of TGF-β, KRAS, ERK1/2, IL1, and INS at early phase and Vegf, HGF, Interferan alpha, IL6, IL1B, and NFKB at the later regeneration phase after ALPPS and PH. On the other hand, compared to PH and PVL, the cell cycle associated genes,  www.nature.com/scientificreports www.nature.com/scientificreports/ cytokines, and transcription regulators, such as E2F1, TBX2, FOXM1, and EP400, are still activated after ALPPS at later regeneration phase. Finally, the Notch signaling, which is a complex signaling pathway that is crucial for the development of multiple organs, was also seemed to be activated in ALPPS. In the liver, it controls the hepatic cell differentiation into, and the formation of, the biliary system 51 . In addition, notch proteins have an extracellular Heatmap that displays surgical procedure as well as temporal changes for the top 50 activated/inhibited upstream regulators. The upstream regulators are ordered according to significance in ALPPS procedure. The color-key is for the z-score. Data were analyzed through the use of IPA (QIAGEN Inc., https://www.qiagenbioinformatics.com/products/ingenuitypathwayanalysis). membrane component with epidermal growth-like factor (EGF) repeats, which may be responsible for the high growth rate seen clinically in ALPPS 52 .
The upstream analyses also indicated the predicted inhibition of NUPR1, CDKN2A, Rb, PAX6, and TP53, at initial phase of liver regeneration 53,54 . However, at the later stage of liver regeneration, TP53 is activated, especially after PH. Furthermore, the mechanistic network of upstream regulators and their predicted relationship based on the observed gene expression changes in our data reveals the working mechanism of TP53 at early and late stage of the liver regeneration, as demonstrated in Fig. 5. p53 regulates liver homeostasis, and initiation of cell proliferation through proliferative signaling and disruption of p53 signaling lead to faster recovery 53,54 .
The limitation of the study is that we examined the liver proliferation at two time points, early phase (i.e. 24 h) and late-stage (i.e. 96 h). Our earlier study of ALPPS and PVL indicated higher proliferation index (PI) at 24 h and 48 hr comparing to 96 h, and there was no significant difference between the two time-points (i.e. 24 h vs. 48 h) 17 . The future remnant liver volume (FRLV) ratio was significantly higher in ALPPS comparing to PVL at 24 hr and 96 hr time points, but not significantly different at 48 h or 1 week. Higher FRL ratio is critical factor in improving surgical outcomes and liver regeneration. In addition, the ALPPS model has significantly more inflammatory cells infiltration at 24 hr comparing to 48 hr; higher infiltration of inflammatory might promoted earlier liver regeneration in the ALPPS model. It also indicated higher portal pressure in ALPPS group at 24 hr comparing to other www.nature.com/scientificreports www.nature.com/scientificreports/ time points; having higher portal pressure might act as physical stressor that contribute to ignite the regeneration process. Xu et al. studied the expressed genes in regenerating rat liver after PH, also reported that temporal patterns of gene expression were similar at 48 h and 96 h after PH 55 . Nevertheless, more time points should worth to be assessed to fully investigate the liver proliferation process. While recognizing this limitation, we believe we have largely achieved our aim, as the major objective of the study is to fish out early phase (i.e. 24 h) and late-stage (i.e. 96 h) liver regeneration molecular markers and examining molecular differences between these phases in the ALPSS, PVL and PH models that have provided unique molecular signatures.
In summary, our study presents a comprehensive transcriptomic profiling of three surgical procedures that are commonly used in clinical practice and identified the inter-procedural and temporal variations in gene expression patterns in each surgical procedure. Identification of molecular signatures and signaling pathways specific to each surgical procedure further our understanding of key regulators of liver regeneration as well as patient populations that are likely to benefit from each procedure.

Data availability
The datasets generated and analyzed during the current study is available at this link (https://www.dropbox. com/s/hgewpv07qce6fed/RNASeqTopHataligned_Normdata_sample29.txt?dl=0).