RNA-Seq Reveals Different Gene Expression in Liver-Specific Prohibitin 1 Knock-Out Mice

Prohibitin 1 (PHB1) is an evolutionarily conserved and ubiquitously expressed protein that stabilizes mitochondrial chaperone. Our previous studies showed that liver-specific Phb1 deficiency induced liver injuries and aggravated lipopolysaccharide (LPS)-induced innate immune responses. In this study, we performed RNA-sequencing (RNA-seq) analysis with liver tissues to investigate global gene expression among liver-specific Phb1−/−, Phb1+/−, and WT mice, focusing on the differentially expressed (DE) genes between Phb1+/− and WT. When 78 DE genes were analyzed for biological functions, using ingenuity pathway analysis (IPA) tool, lipid metabolism-related genes, including insulin receptor (Insr), sterol regulatory element-binding transcription factor 1 (Srebf1), Srebf2, and SREBP cleavage-activating protein (Scap) appeared to be downregulated in liver-specific Phb1+/− compared with WT. Diseases and biofunctions analyses conducted by IPA verified that hepatic system diseases, including liver fibrosis, liver hyperplasia/hyperproliferation, and liver necrosis/cell death, which may be caused by hepatotoxicity, were highly associated with liver-specific Phb1 deficiency in mice. Interestingly, of liver disease-related 5 DE genes between Phb1+/− and WT, the mRNA expressions of forkhead box M1 (Foxm1) and TIMP inhibitor of metalloproteinase (Timp1) were matched with validation for RNA-seq in liver tissues and AML12 cells transfected with Phb1 siRNA. The results in this study provide additional insights into molecular mechanisms responsible for increasing susceptibility of liver injuries associated with hepatic Phb1.


INTRODUCTION
Liver disease is a comprehensive term, including diverse stages of the disorder related to liver injury, and is one of the main causes of illness and death worldwide. According to Globocan 2018 by the WHO, the incidence and mortality of liver cancer are the sixth and third highest in both sexes and all ages, respectively (Bray et al., 2018). Although the risk factors include hepatitis B and C virus, alcohol, obesity, etc., the cellular regulatory mechanisms for liver disease and accountable evidence for physiological responses are yet to be explored (Smalling et al., 2013). It may be crucial to identify specific genes and pathway changes by liver injuries through whole transcriptome analyses of a relevant animal model. RNA-sequencing (RNA-seq) can proceed with a wholegenome survey of gene expression and provide a digital measure of the presence and prevalence of transcripts from known and previously unknown genes (Wang et al., 2010;Benjamin et al., 2014). The objective of the present study was to conduct global gene expression analyses of liverspecific Prohibitin 1 (Phb1) knockout (KO) mice as a liver disease model. PHB1 is an evolutionarily conserved and ubiquitously expressed protein that stabilizes mitochondrial chaperone (Yang et al., 2018). The protein is located to inner mitochondrial membrane and nucleus, so it regulates several important cellular processes, including apoptosis, cell proliferation, and transcriptional regulation by interacting with retinoblastoma protein (Rb), p53, and E2F transcription factors (Nijtmans et al., 2000;Ramani et al., 2016).
Liver-specific Phb1 homozygous KO (referred to as Phb1 −/− ) mouse showed positive staining for OV-6, an oval cell marker, and glutathione S-transferase Pi (GST Pi), a preneoplastic marker, at 3 weeks old and exhibited spontaneous liver injury, fibrosis, and hepatocellular carcinoma (HCC) from 8 months of age (Ko et al., 2010). A follow-up research study demonstrated that lipopolysaccharide (LPS)-innate immune responses were exacerbated by Phb1 deficiency, which mimics Phb1 −/− , in murine macrophages (Jung et al., 2020). Furthermore, another research showed that there were no significant differences in body weight between liver-specific Phb1 heterozygote (referred to as Phb1 +/− ) and wild-type (WT, referred to as Phb1 +/+ ) mice, but Phb1 +/− showed more severe liver damage when feeding a methionine or choline-deficient diet, compared with WT (Heo and Ko, 2019).
Collectively, these results suggested that liver-specific Phb1 deficient mice may be a suitable model to investigate molecular mechanisms related to hepatotoxicity. In an earlier study, we observed phenotypic differences and differential gene expression, using microarray assay between Phb1 −/− and WT (Ko et al., 2010). On the basis of the result, we have also investigated the vulnerability of hepatotoxicity in deletion of Phb1 (Heo and Ko, 2019;Jung et al., 2020). As part of research elucidating the susceptibility, it is necessary to observe the genetic change in liver-specific Phb1 +/− , compared with WT to determine which genetic factors increases the susceptibility of liver injuries. In the current study, we aimed to investigate genome-wide gene expression among liver-specific Phb1 −/− , Phb1 +/− , and WT mice by using the RNA-seq method, focusing on the differentially expressed (DE) genes between Phb1 +/− and WT. Also, we verified the expression level of hepatic disease-related genes in normal murine hepatocytes transfected with Phb1 siRNA to compare RNA-seq data. Through this comparison, DE genes showing greater differences between Phb1 +/− and WT may provide sharper perceptions of alteration of physiological responses in a liver-specific Phb1-deficient mouse as a liver disease model.

RNA-seq and Data Analysis
The RNA quality based on the RNA integrity number (RIN) was assessed using the RNA R6K assay for the Agilent 2200 TapeStation (Agilent Technology, Santa Clara, CA, USA). RIN scores were ranged between 8.5 and 9.5. For the RNA-seq analysis, library preparation and sequencing analyses were carried out at the Research Technology Support Facility of Michigan State University (East Lansing, MI, USA). Transcriptome (n = 6 for each WT, Phb1 +/− , and Phb1 −/− mice, including three males and three females per group) was analyzed, using a 1 × 50 bp single-end-read method of Illumina HiSeq system (Illumina Inc., San Diego, CA, USA) as described in earlier reports (Kong et al., 2017;Lee et al., 2020). The raw reads were aligned with the mouse reference genome (GRCm38.p6) downloaded from NCBI (https://www.ncbi.nlm. nih.gov/assembly/GCF_000001635.26) using ArrayStar program in Lagergene software package (DNAStar Inc., Madison, WI, USA). Total mapped counts were transformed into log 2 values of the number of reads per million (RPM) to stabilize the variance, and then further quantile normalization was performed. Normalized values were subjected to further statistical analyses performed by JMP Genomics 9 (SAS Institute Inc., Cary, NC). The one-way ANOVA statistics was used to compare among WT, Phb1 +/− , and Phb1 −/− mice. Genes showing >2-fold changes (1 in log 2 ) differences and <0.05 false discovery rate (FDR) in the comparisons were considered DE genes. The raw data were submitted to the National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA) database under accession numbers: SRR14876823, SRR14876822, SRR14876815, SRR14876812, SRR14876818, SRR14876817, SRR14876813, SRR14876821, SRR14876820, SRR14876819, SRR14876816, SRR14876814, SRR14876807, SRR14876811, SRR14876810, SRR14876809, SRR14876808, and SRR14876806.

Bioinformatics Pathway Analyses
For pathway analyses, ingenuity pathway analysis (IPA; Qiagen, Valencia, CA; http://www.ingenuity.com) software was used for functional annotation, canonical pathway analysis, upstream analysis, and network discovery. The IPA core analyses are based on previous knowledge of the associations of upstream regulators and their downstream target genes archived in the ingenuity knowledge base. The p-values were calculated by Fisher's exact test for the upstream regulator analysis.

Validation for RNA-seq Using RT-PCR in Liver Tissues
For cDNA synthesis and PCR, a Verso cDNA synthesis kit (ThermoFisher, Waltham, MA, USA) and Promega PCR master mix (Promega Corp., Madison, WI, USA) were used. PCR reactions were carried out at 95 • C for 2 min, followed by 22 or 32 cycles with denaturation at 95 • C for 30 s and annealing and elongation at 55 • C for 45 s and 72 • C for 1 min. Band density was determined, using Image J software (National Institutes of Health, Bethesda, MD, USA). The gene expression levels were normalized to the expression levels of the GAPDH housekeeping gene. Relative quantification was reported as fold changes compared to control samples. The primer sequence of RT-PCR used in the experiment is shown in Supplementary Material 1.

qRT-PCR for Expression of Target Genes
Total RNA was isolated from AML12 cells transfected with siRNAs by using TRIzol reagent (ThermoFisher, Waltham, MA, USA) according to the guideline of the manufacturer. First-strand complementary DNA (cDNA) was prepared from 2 µg of total RNA, using RevertAid First Strand cDNA Synthesis Kit (ThermoFisher, Waltham, MA, USA) according to manufacturer's protocol. Quantitative PCR (qPCR) was performed with the Maxima SYBR Green/ROX qPCR Master Mix (2X) (ThermoFisher, Waltham, MA, USA). The reaction was performed in a volume of 10 µl, containing 50 ng of cDNA, 0.3 µM of each primer (Supplementary Material 2) and 5µl maxima SYBR Green/ROX qPCR Master Mix (2X), which were run in duplicate, using a QuantStudio3 thermocycler (ThermoFisher, Waltham, MA, USA). Amplification was carried out with a template denaturation at 95 • C for 10 min, followed at 95 • C for 15 s and 60 • C for 1 min for 40 cycles. The gene expression levels were normalized to the expression levels of the β-actin housekeeping gene and calculated using the 2 − Ct method.

Statistical Analysis
Data on validation using RT-PCR and qRT-PCR were presented as the mean ± standard deviation. Differences among groups were determined by one-way ANOVA analysis with Duncan posthoc test, using SAS 9.4 (SAS Institute Inc., Cary, NC, USA). Values of p < 0.05 were considered statistically significant.

RNA-seq Results
A total of 18 RNA-seq libraries were constructed using RNA samples extracted from the liver obtained from WT, Phb1 +/− , and Phb1 −/− . After filtering low read counts and normalization, a total of 8,982 transcripts remained. Of those, only mRNAs showing FDR <0.05 and log 2 fold change >1 or <−1 (linear fold change >2 or <−2) were considered differentially expressed (DE) genes for further bioinformatics analyses and PCR validation. Using these criteria, this study was able to select 78 (15 up-and 63 downregulation) DE genes in Phb1 +/− mice, compared with WT. The complete gene list for 78DE genes is provided in Supplementary Material 3. Several genes, including optic atrophy 1 (Opa1), mitofusin 1 (Mfn1), and mitofusin 2 (Mfn2), which are involved in mitochondrial fusion and fission, contributing to mitochondrial morphology, were not differentially expressed between Phb1 +/− and WT. Expression of mitochondrial cytochrome c oxidase (Cox) responsible for aerobic energy generation was also no difference in Phb1 +/− , compared with WT (data not shown). But these mitochondrialrelated genes are downregulated in Phb1 −/− (data not shown). This result can indicate that half of the reduction in Phb1 expression may not affect mitochondrial dysfunction.

The 10 Most DE Genes in Hetero Compared With WT
The full name of the 10 most DE genes in the Phb1 +/− group compared WT (plus FC values and FDR in KO compared with WT) and their functional characteristics are shown in Tables 1, 2. Of those, upregulated genes were involved in nucleolar functions (Nol3, Hist1h2b, and Snora17), xenobiotic enzyme (Sult1e1), ubiquitin peptidase (Usp2), and circadian functions (Ciart, Per2). The 10 most downregulated genes appeared to be involved "Hetero" and "KO" in tables represent "Phb1 +/− " and "Phb1 −/− ", respectively.
The mRNA expression levels of all 5 DE genes were well-matched between Phb1 +/− and WT in terms of the direction (Table 3; Figures 1A,B). But, in in vitro validation, the mRNA expression levels of Foxm1 and Timp1 were only matched with RNA-seq and RT-PCR validation in liver tissues (Table 3; Figure 1C), which may be due to differences between animal tissues and cell culture. These results confirmed both validity of RNA-seq data for further pathway analyses and the reliability of our liver disease model, indicating that half deficiency of Phb1 in the liver can elucidate the altered molecular mechanisms related to liver injuries.

Bioinformatic Pathway Analyses: Upstream Regulators and Molecular Networks
The 78 DE genes in Phb1 +/− compared with WT were subjected to core pathway analyses, using ingenuity pathway analysis (IPA) software. Upstream regulator analysis was conducted with enriched pathways and p-value computed by IPA literature review algorithm. The resulting causal networks were scored for known physiological relationships with diseases, functions, genes, or chemicals in the liver tissue and liver cell lines. Causal networks derived by participating the regulator molecule, which controls the expression effector molecules in the given dataset. Thus, causal networks are small hierarchical networks of regulators that control the expression of the given data set. Upstream regulator analyses containing networks provide positive and negative z-scores (≥2 and ≤−2, respectively), which are considered significant predictions of activation and inhibition, respectively. Results showed that sterol regulatory element-binding transcription factor 2 (Srebf2) (zscore = −3.12), Srebf1 (z-score = −3.08), SREBP cleavageactivating protein (Scap) (z-score = −2.98), insulin receptor (Insr) (z-score = −2.62), ATPase copper transporting beta (Atp7b) (z-score = −2.65), and nuclear receptor coactivator 2 (Ncoa2) (z-score = −2.00) were predicted to be inhibited, while peroxisome proliferator-activated receptor alpha (Pparα) (z-score = 2.16) was predicted to be activated in liver of Phb1 +/− mice ( Table 4). Figure 2 showed these three potential interactions were centered with: (1) Scap, Srebf1, Srebf2; (2) Atp7b, Ncoa2; (3) Pparα, Ncoa2. Similarly, causal network

Gm40498
• Predicted gene Hist1h2bm • Is a protein which is intronless and encodes a member of the histone H2B family Sult1e1 • Is an estrogen sulfotransferase (EST), which is a major sulfotransferase isoform • Expresses in multiple human tissues • Regulates estrogen homeostasis by sulfonating and deactivating estrogen (Guo et al., 2015) Ciart • Forms a complex with other clock proteins and modulates the circadian machinery (Hatanaka and Takumi, 2017) • Associates mutant human CIART gene with grade 4 astrocytoma in human brain  Gm38832 • Predicted gene • Encodes the coagulation factor XIII A subunit which is a proenzyme of plasma transglutaminase consisting of catalytic A (FXIII-A) and non-catalytic B subunits (FXIII-B) • Congenital or acquired FXIII-B deficiency may result in increased bleeding tendency through impaired fibrin stabilization (Souri et al., 2015) Timp1 • A natural inhibitor of the matrix metalloproteinases (MMPs), a group of peptidases involved in degradation of the extracellular matrix, but able to exert multiple effects on cell growth, proliferation, differentiation and apoptosis, in an MMP-independent manner • TIMP1 deficiency exacerbates carbon tetrachloride-induced liver injury and fibrosis in mice  • Displays complete pro-tumor activities through its upregulation in liver tissue and serum from patients and mouse model with liver disease (Yoshiji et al., 2000) • Shows MMP-independent role of TIMP1 at the blood brain barrier during viral encephalomyelitis (Savarin et al., 2013), or regulation of adipogenesis of adipose-derived stem cells via Wnt singnaling pathway (Wang et al., 2020) Shc4 • Involves in coupling receptor tyrosine kinases to the Ras-mitogen activated protein kinase signaling pathway, and to have a predominant cytoplasmic distribution (Ahmed and Prigent, 2014) • Acts non-canonically to promote phosphorylation of select epidermal growth factor receptor (EGFR) residues (Wills et al., 2014) Abcc12 • Is a member of the superfamily of ATP-binding cassette (ABC) transporters and the MRP subfamily which is involved in multi-drug resistance • In all the breast cancer tissues, the expression of ABCC12 gene 3.74 times higher than in controls (Esmaeili et al., 2018) Tubb2b • Microtubules, which is a key participant in processes such as mitosis and intracellular transport, are composed of heterodimers of alpha-and beta-tubulins • Among taxane-based chemotherapy group, cases with higher β-tubulinIII expression were associated with a significantly more favorable prognosis compared with those having lower β-tubulinIII expression (Aoki et al., 2009) regulator analysis showed that Srebf2 (z-score = −3.16), Srebf1 (z-score = −3.16), Scap (z-score = −3), Insr (z-score = −2.65), Atp7b (z-score = −2.65), insulin-like growth factor 1 (Igf1) (zscore = −2.53) were predicted to be inhibited, while arachidonic acid (z-score = 2.31) was predicted to be activated (Table 5).
Srebf1, Srebf2, and Scap, which are involved in regulating cellular lipid metabolism, were listed as the highest ranks of causal network regulators. The predicted downstream target molecules of arachidonic acid were mitogen-activated protein kinase 14 (Mapk14) and Srebf1, followed by downregulating molecules involved in cholesterol and triglycerides (TG) synthesis such as farnesyl diphosphate synthase (Fdps), cytochrome P450 family 51 subfamily A member 1 (Cyp51a1), and Scd1 ( Figure 3A). Igf1 was predicted to be inhibited and downregulated similar genes associated with cholesterol and triglycerides via Insr ( Figure 3B). Also, Timp1 is one of the target molecules of Igf1 via signal transducer and activator of transcription 3 (Stat3) (Figure 3B). The most critical diseases and biofunctions (Table 6), molecular and cellular functions ( Table 7), and physiological system development and functions ( Table 8) were identified with influences of 78 DE genes, and specific genes involved in all diseases and functions were listed in Supplementary Material 6. Of diseases and biofunctions listed in Table 6, hepatic system diseases were ranked as the second. The molecules, including Scd1, Acta2, collagen type I alpha 1 chain (Col1a1), Timp1, Foxm1, Sult1e1, and plasminogen activator (Plat), are also involved in hepatotoxicity pathways, such as liver fibrosis, liver hyperplasia/hyperproliferation, liver proliferation, and liver necrosis/cell death (Table 9). Molecules associated with hepatotoxicity were shown in Figure 4. Hepatic DE genes in Phb1 +/− mice showed increased risks of hyperplasia of the liver, apoptosis of liver cells, hepatic injury, and hepatoma and decreased functions in the proliferation of liver cells and fibrosis. The hepatotoxicity analysis indicates that Phb1 +/− mice may become more susceptible with liver damages/injuries. Molecular network analysis produced five potential molecular interactions among 78 DE genes ( Table 10). The most interesting network #1 centered with regulatory molecules, including Pparα and Srebf1, which were discussed in upstream regulators above, in addition to tumor necrosis factor (TNF), which is predicted to be inhibited. Molecules interacting network #1 are connected to functions of lipid metabolism, small molecule biochemistry, and vitamin and mineral metabolism (Table 10; Figure 5).

DISCUSSION
In this study, we examined the hepatic transcriptome of liverspecific Phb1 deficient mice and identified biological functions and genes associated with hepatotoxicity, lipid metabolism, and metabolic disorders. IPA suggested a network related to hepatotoxicity, which is composed of Foxm1, Timp1, Usp2, Scd1, and Sult1e1. These data exhibited that decreased Phb1 in hepatocytes contributed to abnormal proliferation of various liver cell types and cancer cell transformation by upregulated or downregulated genes. Also, IPA analysis with DE genes suggested increasing arachidonic acid and suppressing Igf1-mediated signaling as causal networks in liver-specific Phb1 +/− , compared with WT. They had in common with downstream targets as Timp1, Scd1, Fdps, and Cyp51a1 (Figure 3), which are mainly associated with the extracellular matrix, fatty acid metabolism, and steroid biosynthesis (Ge et al., 2020;Yan et al., 2020) and downregulated via Mapk14, Srebf1, Stat3, and Insr. This finding is in line with the blocking of activation of phosphoinositide 3kinase by insulin via p38 Mapk in hepatocytes with arachidonic acid (Talukdar et al., 2005) and evidence that arachidonic acid has the potential to decrease insulin-mediated activation of Srebp-1c by inhibiting liver X receptor (Lxr) activation in rat hepatocytes (Chen et al., 2004). Also, Igf1 and Insr-mediated Cyp51a1 were predicted as a downregulated enzyme, and related research reported that hepatocyte Cyp51KO mice showed inflammation and mild-to-moderate portal fibrosis and abnormal hepatic sterol metabolism (Urlep et al., 2017). Collectively, these results may elucidate physiological changes by Phb1 and provide a comprehensive understanding of prognosis and prevention of liver injuries, containing liver cancer.

Liver Disease-Related Genes
Of 78 DE genes, which were filtered out from liver-specific Phb1 +/− vs. WT, Foxm1, Timp1, Usp2, Scd1, Sult1e1, are known to be associated with hepatotoxicity (Figure 4). Foxm1, a transcription regulator, is located to the nucleus and highly expressed in proliferating normal cells and various cancer cells (Teh, 2012). The knockdown of FOXM1 in human hepatocellular carcinoma (HCC) cell lines significantly alleviated cell cycle arrest and cell growth suppression (Hu et al., 2014). In addition, overexpression of FOXM1 is associated with an aggressive tumor feature and poor prognosis of HCC (Sun et al., 2011). Unlike earlier results of positive correlation between Foxm1 and HCC, our transcriptomic profile demonstrated that fold change of Foxm1 was rather reduced in Phb1 +/− compared with WT (Phb1 +/− vs. WT, log 2 FC = −1.6) (Supplementary Material 3) and then increased in Phb1 −/− compared with Phb1 +/− (Phb1 −/− vs. Phb1 +/− , log 2 FC = 4.6) (data not shown). This finding was consistent with validation using the PCR method in liver tissues and AML12 cells transfected with siPhb1 (Figure 1). Deletion of Foxm1 is correlated with Ras-induced HCC with stem cell features by accumulating reactive oxygen species (ROS) (Kopanja et al., 2015). Another study revealed a negative association with Foxm1/nuclear factor kappa B (NF-kB) and methionine adenosyltransferase 1A (Mat1a) that is responsible for catalyzes in the conversion of methionine to homocysteine (Li et al., 2020). It was reported that Foxm1 directly interacts with Mat1a, which is a tumor suppressor in the liver. Collectively, it can be interpreted that Phb1 leads to a unique change of the expression pattern of Foxm1 and may be a key regulator, modulating the transcriptional activity of Foxm1 in terms of tumorigenesis. The mRNA expression of selected five DE genes in each group of three males and three females (n = 6/ genotype) was detected by RT-PCR. (B) Band density from RT-PCR was quantified with Image J software and represented as the fold change of control and normalized Gapdh. (C) The mRNA expression levels of selected five DE genes in AML12 cells transfected siPhb1 (n = 9) were detected by qRT-PCR and normalized b-actin. "LE" means "low efficiency", represented Phb1 +/− , and "HE" means "high efficiency", represented Phb1 −/− . All data were expressed as means ± standard deviation. One-way ANOVA followed by Duncan's post hoc test was performed, and differences were considered statistically significant. *p < 0.05 vs. wild type (WT) and control, respectively. # p < 0.05 vs. Phb1 +/− and low efficiency (LE). $ p < 0.05 vs. scramble.
Likewise, Timp1 expression in the transcriptomic profile and PCR validation was similar to the results of Foxm1 expression. Previous studies explained increasing Timp1 expression as a marker of extrahepatic and intrahepatic tumors from liver tissues and serum (Ylisirniö et al., 2000;Yoshiji et al., 2000;Yukawa et al., 2007). Also, Thiele et al. (2017) suggested that hepatic Timp1 mRNA expression from WT mice was upregulated in HCC tissue compared with adjacent paired normal tissue when co-treated with diethylnitrosamine and CCL 4. In agreement with these studies, our transcriptomic profile showed increased Timp1 expression in Phb1 −/− , compared with Phb1 +/− (log 2 FC = 10.8) ( Table 1). On the contrary, Timp1 expression was rather decreased in Phb1 +/− , compared with WT, which was similar to Foxm1 (log 2 FC = −4.4) ( Table 1). Wang et al. (2011) demonstrated the newly hepatoprotective role of Timp1 during acute and chronic liver injuries, which is positively regulated by the interleukin-6 (IL-6)/STAT3 signaling pathway. Our liver-specific Phb1-deficient mouse showed different Timp1 expression patterns according to Phb1 depletion, which encompassed both protumor-and antitumor activity of Timp1. It was previously reported that Timp1 was increased in the liver from adipocyte-specific nuclear form of SREBP-1c (nSREBP-1c) transgenic mice, showing human nonalcoholic steatohepatitis (NASH) (Kakino et al., 2018). However, Upregulation of Timp1 was reduced in the liver from the Tnf −/− nSREBP-1c transgenic mice (Kakino et al., 2018). Tomita et al. explained that control mice-fed-MCD diet rapidly increased mRNA expression of Timp1 in the whole liver, compared with both TNF receptors 1 (Tnfr1) and 2 (Tnfr2) knockout mice, representing Tnfr-double KO mice    (Tomita et al., 2006). We can interpret that decreased Timp1 may act as an initial signal of liver injuries by Phb1 +/− , and then chronic liver injuries by Phb1 −/− may represent a rapid increase of Timp1. Since we recognize the importance to uncover a direct correlation between Phb1 and Timp1, further study is underway to assess the association. Nevertheless, both RNA-seq and PCR validation in siPhb1-transfected hepatocytes showed indirect inverse relevance between Phb1 and Timp1 in this research. Especially, Ko et al. already explained the relation using microarray analyses (Ko et al., 2010), showing an increase in the expression of Timp1 in liver-specific Phb1 −/− mice. Therefore, Timp1 may be an important enzyme that predicts the extent of liver injuries by lack of the Phb1 gene. Further investigations are needed to characterize functional connections between Phb1 and Timp1 in the liver disease model. Despite increased expression of Foxm1 and Timp1 in an environment with liver injuries, another possibility to interpret the differential gene expression in Phb1 +/− vs. Phb1 −/− may be deciphered by genetic compensation by the gene knockout (El-Brolosy and Stainier, 2017). Genetic compensation means that another gene takes over the role of the knockedout gene (De Souza et al., 2006). This concept is used for explaining the discrepancy between knockdown and knockout phenotypes. For example, it was reported that upregulated Emilin 3, which shares the functional domain with egfl7 represented for endothelial extracellular matrix genes, was observed only in the zebrafish egfl7 knockout, not in the knockdown, demonstrating that Emilin 3 may compensate for the egfl7 knockout (Rossi et al., 2015). Besides, both Prnp and Sprn proteins share biological functions in early embryogenesis according to the need for either Prnp or Sprn expression, showing that one can compensate for the absence of the other in PrP-knockout mammals (Young et al., 2009). According to these cases in which gene expression is reversely regulated in the knockout and knockdown comparison of the unique expression pattern of Foxm1 and Timp1 in Phb1 +/− and Phb1 −/− may be explained by genetic compensation. To specify the concept in our liver disease model, further intensive research is necessary in the environment with regulating the expression of Phb1. Ubiquitin carboxyl-terminal hydrolase 2 (Usp2), which is upregulated in Phb1 +/− compared with WT (log 2 FC = 1.3) ( Table 1), is a deubiquitinating enzyme and involved in controlling activity and levels of protein under certain physiological conditions. Usp2 has been reported to stabilize mouse double minute 2 (Mdm2) and Mdm4, which are protooncogenes with both p53-dependent or independent activities, to degrade p53 (Benassi et al., 2012;Young et al., 2019), maintaining the expression of p53 at a low level. Specifically, protein-protein interaction assays using the bacterial two-hybrid system showed that USP2 can deubiquitinate MDM2 and promotes p53 degradation, showing the association between USP2 and MDM2 (Stevenson et al., 2007). Especially, it was reported that transfection of tumor-derived cell lines with siUsp2 resulted in increasing p53 protein expression and its target gene, p21 (Stevenson et al., 2007). PHB1 has an ability to physically interact with p53 and enhance p53-mediated transcriptional activation by promoting its recruitment to promoters in two breast cancer cell lines where it co-localizes with p53 (Fusaro et al., 2003). Our transcriptome analysis showed that the expression level of Usp2 was upregulated in Phb1 +/− compared with WT (log 2 FC = 1.3) and then decreased in Phb1 −/− compared with Phb1 +/− (log 2 FC = −1.8) ( Table 1). We can infer that Phb1 may participate in enhancing Usp2-mediated proteosome activity, and Phb1 deficiency may result in abnormal p53 stabilization and activation (Haupt et al., 1997;Honda and Yasuda, 2000). Thus, we can explain that there is no defense response against a high-stress environment, such as an increasing susceptibility of liver injuries by the drastic shortage of Phb1, represented to Phb1 −/− . Stearoyl-CoA desaturase-1 (SCD1) plays an important role in the conversion of saturated fatty acids to monosaturated fatty acids and further synthesis of TG (Lounis et al., 2016). We expected that Scd1 expression in Phb1 +/− , compared with WT, also increased in transcriptome analysis, since our previous research on the effect of Phb1 deficiency with siPhb1 knockdown on impaired lipid metabolism demonstrated that palmitic acid promoted Scd1 mRNA expression levels (unpublished data, under review). In the current study, Scd1 expression was decreased in Phb1 +/− compared with WT (log 2 FC = −1.4) (Supplementary Material 3), and there was no statistically significant difference between Phb1 +/− and Phb1 −/− ( Table 3). Wang et al. (2002a) reported that PHB1-mediated transcriptional repression required histone deacetylase (HDAC), and additional corepressors like N-CoR are involved. Also, E2F transcription factors, which play a key role in regulating mammalian cell cycle progression, were controlled by PHB1, interacting with retinoblastoma protein (Rb) (Wang et al., 2002b;Mishra et al., 2006). Recently, combining research using lipidomics and transcriptome analyses in Rb depletion in mouse embryonic fibroblast has shown that Rb deficiency increases the concentration of fatty acid, acylcarnitine, phosphatidylcholine, and arachidonoyl ethanolamine (Muranaka et al., 2017). Also, Scd1 is most strongly controlled by Rb, possibly through E2F and SREBP transcription factors (Muranaka et al., 2017). Collectively, a decrease in Scd1 expression in Phb1 +/− compared with WT may be due to the interaction between PHB1 and SCD1 at the protein level, maintaining a balance between lipid catabolism and anabolism.
Sulfotransferase family 1E member 1 (SULT1E1) is a cytosolic enzyme, which is called as an estrogen sulfotransferase and is responsible for the inactivation of β-estradiol (E2) involved in changing estrogen metabolism and liver function (Li and Falany, 2007;Li et al., 2009;Xu et al., 2018). According to a study on cystic fibrosis, which is an inherited disorder and affects the lung, intestine, pancreas, and livers, hepatic Sult1e1 activity was increased in cystic fibrosis transmembrane conductance regulator (CFTR) −/− mice compared with CFTR +/+ (Li and Falany, 2007). Also, the estrogen receptor α protein level was reduced in CFTR −/− mice with high Sult1e1 activity, showing high affinity between E2 and substrate of SULT1E1 and thereby altering the level of estrogen-regulated proteins (Song, 2001;Li and Falany, 2007). Wang et al. (2004) demonstrated the association of PHB1 with E2F transcription factor 1 (E2F1) and the repressive function of estrogen antagonists in human breast cancer cells (Wang et al., 2004). Moreover, the PHB1 gene and protein expression were both markedly enhanced by estrogen antagonists (Wang et al., 2004). In agreement with these studies, Sult1e1 expression in the current transcriptomic profile was increased by depletion of Phb1 in between Phb1 +/− and Phb1 −/− compared with WT (log 2 FC = 2.7 and 2.4, respectively) (Supplementary Material 3). Both transcriptome analyses and validation, using AML12 cells transfected with siPhb1, showed a negative relationship between Phb1 and Sult1e1 ( Figure 1C). We can infer that Phb1 deficiency is likely to have relevance to increasing expression of Sult1e1. Also, Phb1 depletion may be a crucial status to elevate hepatic Sult1e1 activity responsible for repressing estrogen-related genes in terms of its transcriptional activities. Thus, the association of Phb1 depletion with its effect on estrogen metabolism in livers needs to be investigated further.

Lipid Metabolism and Metabolic Diseases
In addition to the hepatotoxicity of 5 DE genes between liver-specific Phb1 +/− and WT, there are some associations with lipotoxicity. The role of fatty acid metabolism in cancer initiation, progression, and drug resistance through the FOXO3-FOXM1 axis had been mentioned (Saavedra-García et al., 2018). Besides, compared with adipose tissue from controls, toll-like receptor 2 (Tlr2) KO mice fed a high-fat diet that decreased levels of Timp1, collagen 1, and transforming growth factor-β1 (TGFβ1) (Song et al., 2018). Usp2 enhances the stability of    fatty acid synthase (FASN) by impeding proteasome-dependent degradation in human HCC (Calvisi et al., 2011;Kitamura and Hashimoto, 2021). The high-fat diet-induced NAFLD rat model showed increased protein levels of SULT1E1, compared with mice-fed normal standard diet through proteomics analysis and western blots analysis (Cong et al., 2021). From upstream regulator analysis and causal network regulators, Insr, Scap, Srebf1, and Srebf2 were predicted to be inhibited in liver-specific Phb1 +/− compared with WT (Tables 4, 5). SREBPs regulate the biosynthesis of triglyceride (TG), fatty acids (FAs), and cholesterol, which can increase the expression of genes related to lipid synthesis and lipid uptake (Brown and Goldstein, 1997). SREBF1 and SREBF2 encode three SREBP isoforms: SREBP-1a, SREBP-1c, and SREBP-2, respectively. SREBP-1a and−1c are usually expressed in the liver, adipose tissue, and adrenal gland of mice and humans, whereas SREBP-2 is ubiquitously expressed in cell lines, spleen, and intestinal tissues (Jeon and Osborne, 2012;Moslehi and Hamidi-Zad, 2018). When cellular sterol levels are high, SCAP, which is called ER membrane protein, interacts with SREBPs in the endoplasmic reticulum (ER) membrane, and the SREBP/SCAP complex moves to Golgi to translocate to the nucleus and bind to the target gene promoters (Xiao and Song, 2013;Moslehi and Hamidi-Zad, 2018). SREBP-1a and−1c mostly activate fatty acid and TG synthesis, and SREPB-2 increases transcription of genes related to cholesterol synthesis and uptake. Remarkably, SREBP-1c is known to be controlled by insulin treatment. Hegarty et al. (2005) showed that full induction of mature and transcriptionally active form of SREBP-1 in rat hepatocyte resulted from insulin treatment (Hegarty et al., 2005). Many studies reported that insulin is one of the most potent activators of SREBP-1c (Shimomura et al., 1999;Azzout-Marniche et al., 2000;Zeng et al., 2017). We found that Insr, FIGURE 4 | A network of DE genes affecting functions related to "hepatotoxicity" between Phb1 +/− and WT. Red, green, orange, and blue represent upregulation and downregulation, predicted upregulation, and predicted downregulation, respectively.
SREBPs, and Scap were predicted to be inhibited in liver-specific Phb1 +/− , compared with WT (Tables 4, 5). Potential decrease of Insr may interpret Insr-related insulin resistance in non-adipose tissues like the liver, contributing to metabolic diseases. The liver regulates various metabolism by insulin, which is triggered by INSR (Samuel and Shulman, 2012;Knebel et al., 2018). Liver-specific Insr-KO mice showed dramatic insulin resistance, abnormal glucose homeostasis, and hyperinsulinemia (Michael et al., 2000;Miao et al., 2014). In addition, according to a study on patients with simple steatosis and NASH, downregulation of SREBP-1c may be associated with the development of burned-out NASH through decreasing quantification of SREBP-1c positive hepatocyte nuclei and increasing mature SREBP-1c levels by immunoblot analysis (Nagaya et al., 2010). The study demonstrated that hepatic expression of SREBP-1c is increased in simple steatosis but gradually decreases with fibrosis progression. Based on these findings, we can speculate that Phb1 is relevant to the downregulation of hepatic expression levels of SREBPs and Insr, which control lipid metabolism and metabolic homeostasis by insulin. However, SREBPs were not predicted as upstream regulators between Phb1 −/− and Phb1 +/− (data not shown); further study is necessary to investigate the hepatic expression level of Insr and SREBPs in our liver disease model for elucidating contribution by gradual depletion of Phb1 on Insr-related insulin resistance and pathogenesis of liver injuries.
Peroxisome proliferator-activated receptor alpha (Pparα) was predicted as an increased upstream regulator in liver-specific Phb1 +/− compared with WT (z-score = 2.16) ( Table 4). PPARα is a key regulator of lipid oxidation, which regulates genes involved in lipid and glucose metabolism and inflammation in the liver (Li and Palinski, 2006). PPARα agonists, such as fenofibrate and bezafibrate, are widely used to treat dyslipidemia, preventing hepatic steatosis and improving insulin sensitivity (Chou et al., 2002;Harano et al., 2006). We expected that if Pparα is an upstream regulator, it will decrease its function by Phb1 deficiency. Because our previous study showed that the Pparγ mRNA expression level, one of the markers of hepatic steatosis, was increased in normal murine hepatocytes transfected with siPhb1, which mimics liver-specific Phb1 −/− (unpublished data, under review). On the other hand, our current study depicted that a slight decrease in Phb1 may attribute to increase Pparα. This trend was consistent with SREBPs as other upstream regulators. In a previous study on a patient with liver diseases, hepatic expression of SREBP-1c was increased in simple steatosis (SS) but gradually decreased in mild NASH and advanced NASH. Likewise, the hepatic mRNA level of PPARα was significantly increased in SS and decreased in the mild NASH and advanced NASH (Nagaya et al., 2010). Taken together, we can infer that Phb1 deficiency at the early stage reflects a defective mechanism via Pparα against abnormal lipogenesis and/or lipid Drug metabolism, glutathione depletion in liver, endocrine system development and function FIGURE 5 | A network relevant to "lipid metabolism", "small molecule biochemistry", "vitamin and mineral represented in a gene network (network ID 1). Red, green, orange, and blue represent upregulation and downregulation, and predicted upregulation and downregulation, respectively.
catabolism. But dramatic depletion of Phb1 is no capable of preventing impaired lipid metabolism and rapidly exacerbates liver injuries.

CONCLUSIONS
In summary, our study demonstrated that the degree of depletion in Phb1 reflects different physiological responses.
We characterized that Foxm1 and Timp1 can provide a basis for investigating a molecular mechanism on the increased susceptibility of liver injuries by Phb1 deficiency. Additionally, SREBPs and Pparα may be key genes to elucidate the pathogenesis of liver diseases relevant to Phb1 deficiency in terms of lipid metabolism. Taken together, these insights may lead to the establishment of novel therapeutic strategies against liver diseases containing liver cancer.

DATA AVAILABILITY STATEMENT
The data generated in this study can be accessed from NCBI using the accession numbers of the project and samples PRJNA738620, SAMN19735929.

ETHICS STATEMENT
The animal study was reviewed and approved by Cedars-Sinai Medical Center IACUC (The IACUC protocol number is 9370).