Transcriptomic and proteomic sequencing unveils the role of vitamin D and metabolic flux shifts in the induction of human hepatic organoids

Abstract

Background

Hepatic organoids (HOs), validated through comparative sequencing with human liver tissues, are reliable models for liver research. Comprehensive transcriptomic and proteomic sequencing of HOs throughout their induction period will enhance the platform’s utility, aiding in the elucidation of liver development’s molecular mechanisms.

Methods

We developed hepatic organoids (HOs) from embryonic stem cells (ESCs) through a de novo induction protocol, mimicking the stages of fetal liver development: ESCs to definitive endoderm (DE), then to foregut (FG), hepatoblasts (HB), and finally to HOs stage 1 (HO1), culminating in self-organizing HOs stage 2 (HO2) via dissociation and re-inoculation. The successful establishment of HOs was validated by immunofluorescence staining and RT-qPCR for specific markers. Comprehensive transcriptomic and proteomic sequencing and analysis were conducted on FG, HB, HO1, and HO2.

Results

Our data suggest that several transcription factors (TFs) activated during the HB stage share overlapping target genes with the vitamin D receptor (VDR). Calcitriol, a direct activator of VDR, notably facilitated the FG to HB stage transition by activating VDR and enhancing key TFs, thereby promoting hepatic progenitor cell maturation. Furthermore, our findings revealed a significant transition towards glycolytic energy metabolism at the HO2 stage, characterized by increased glycolytic flux and reduced oxidative phosphorylation. Inhibition of glycolysis using 2-deoxy-D-glucose (2-DG) led to suppressed growth and differentiation at the HO2 stage. Analysis of signaling pathways indicated upregulation of the HIF-1 pathway, which is associated with glycolysis activation, as well as the MAPK and PI3K-AKT pathways, which regulate HIF-1α protein translation.

Conclusions

We elucidated a pivotal role for calcitriol in facilitating the transition from FG to HB by activating VDR and augmenting the expression of critical transcription factors (TFs). Besides, our research underscores a shift in metabolic pathways toward glycolytic energy metabolism in HO2 organoids. Overall, our multiomics approach reveals the intricate molecular regulation during the development of HOs.

Graphical abstract

Background

The liver, essential for glycolipid metabolism, protein synthesis, and bile acid metabolism, undergoes a significant transformation during embryonic development, shifting from an early hematopoietic organ to a metabolic organ in the later stages. Understanding the molecular mechanisms governing hepatocyte fate determination, differentiation, and maturation is crucial for advancing liver regenerative medicine. However, current knowledge on the transcriptomic, proteomic, and metabolic states during liver development is largely based on model organisms. The molecular landscape of the human early embryonic liver remains largely unexplored, constituting a significant “black box” [1].

The successful creation of hepatic organoids (HOs) from pluripotent stem cells (PSCs) provides a robust in vitro for replicating the in vivo liver differentiation and maturation process, guiding PSCs through definitive endoderm (DE) and posterior foregut (PFG) to form hepatoblast (HB) organoids, which are then matured into fully functional HOs through specialized culture media [2,3,4]. Over the past decade, following the successful creation of HOs, hailed as “the optimal in vitro model” [5], they have been effectively employed in a range of studies, including cell drug toxicity testing [4, 6, 7], modeling of genetic and metabolic diseases [2, 8,9,10], viral hepatitis research [11], and investigation into liver regeneration [12, 13].

In pioneering research on HOs construction, several studies have employed transcriptomic bulk-Seq and single-cell sequencing to analyze the simulation of organoids. Guan Y et al. compared the transcriptomic profiles of immature and fully differentiated mature organoids with those of human liver tissue, confirming that the organoids closely resemble liver tissue in genes related to mature liver function [2]. Shinozawa T and colleagues conducted a comparative single-cell sequencing analysis of HOs derived from PSCs and primary human liver cells, revealing that the HOs contain a diverse range of zonal hepatocytic populations, closely mimicking the cellular heterogeneity of primary adult liver cells [4]. However, few studies sequence the entire induction process of HOs in stages, limiting insights into liver fate determination and differentiation. Furthermore, there is a relative scarcity of proteomic sequencing, which could directly shed light on the functional mechanisms at the level of biological macromolecules. Identifying key molecules that regulate stage-specific fate decisions, including transcription factors and small metabolites, is crucial for understanding the fine-tuned control of organ development.

This study aims to comprehensively explore the molecular mechanisms and metabolic changes during hepatocyte differentiation and maturation in HOs. We use integrated transcriptomic and proteomic sequencing alongside analyses of transcription factor (TF) activity and metabolic flux. Our findings indicate that the transcription factors (TFs) activated during hepatocyte fate transition may be regulated by the vitamin D receptor (VDR) pathway. Additionally, during the HO2 stage, the cells exhibit a glycolysis-dominant energy metabolism state.

Methods

Cell acquisition and HOs culture

Human ESC line H1 [14] was maintained on feeder-free conditions in ncTarget media (shownin, China) in six-well plates (Corning, USA) coated with Matrigel (Corning, 354234). The cells were cultured at 37 °C, 5% CO₂.

The differentiation of ESCs into HOs was conducted following the protocol previously published by the Peltz lab [2]. To generate HO1s, ESCs were dissociated into single cells using Solase (Shownin) and then reaggregated in ultra-low-attachment six-well plates (Corning) with DE differentiation medium from days 0 to 3. This medium included DMEM/F12 supplemented with 0.1 mM nonessential amino acids (Gibco), 1 mM pyruvate (Sigma, USA), 100×ITS (Gibco), 100 ng/mL Activin A (Novoprotein, China), 10 ng/mL BMP4 (Novoprotein), and 2.5 μM Blebbistatin (Shownin). From day 3 to 6, RPMI1640 medium, enriched with B27 (Gibco), 50 ng/mL FGF10 (Novoprotein), and 2% growth factor reduced (GFR) Matrigel (Corning, 354230), was used to induce foregut (FG) spheroids. Between days 6 and 9, the medium was further supplemented with 10 ng/mL FGF10 and 10 ng/mL BMP4 to promote hepatoblast (HB) differentiation. Subsequently, until day 19, HB spheres were induced into primary hepatic organoids (HO1) using a medium composed of HCM (Lonza, Australia), 50 ng/mL HGF (PeproTech, USA), 50 ng/mL OSM (R&D, USA), 10 μM dexamethasone (DEX, Sigma), and 1% GFR Matrigel. On day 20, HO1s were dissociated into single cells using 1 mg/ml collagenase and 0.25% trypsin-EDTA, then reaggregated in GFR Matrigel. The cells were seeded at a density of 1000 cells/35 μl gel per well in the center of a 24-well plate. After Matrigel solidification, 1 mL of growth medium, containing RPMI1640, 50×B27, 250 nM LDN-193,189 (MCE), 3 μM CHIR99021 (Tocris, UK), 10 μM A83-01 (iTGF, MCE), 100 ng/mL EGF (Novoprotein), 10 ng/mL FGF10, and 20 ng/mL HGF, was added and maintained for 5 days. The organoids were then cultured for an additional 6 days in differentiation medium, which consisted of HCM supplemented with 10 μM DAPT (iNOTCH, MCE), 10 ng/mL OSM, 20 ng/mL HGF, 10 μM DEX, and 10 ng/mL BMP4. Mature hepatic organoid 2 (HO2) was harvested using a cell recovery solution (Corning).

RNA extraction and RT-qPCR

Total RNA was isolated using the RNAfast200 RNA Extraction Kit (Fastagen, China), and cDNA synthesis was conducted with the SureScriptTM First-Strand cDNA Synthesis Kit (GeneCopoeia, USA). Real-time quantitative PCR (RT-qPCR) was performed using the BlazeTaq™ SYBR^® Green qPCR Kit (GeneCopoeia) on a Real-Time Thermocycler (Analytik Jena AG, qTOWER3G, Germany), with gene expression levels normalized to GAPDH. The primer sequences utilized in this study are detailed in Table S1.

Immunofluorescence

Cell spheres and organoids at various stages were fixed using 4% paraformaldehyde (PFA) at 4 °C overnight, embedded in OCT compound and sectioned at 8 μm thickness. The sections were permeabilized with 0.3% Triton X-100 for 10 min, and blocked with 5% bovine serum albumin (BSA) for 1 h at room temperature (RT). They were then incubated with primary antibodies overnight at 4 °C, followed by secondary antibodies for 1 h at RT, and finally mounted with DAPI. A comprehensive list of the primary and secondary antibodies utilized is provided in Table S2. All immunofluorescence images were captured using a confocal microscope (LSM980, Germany).

Western blot

The primary and secondary antibodies utilized are detailed in Table S2. Relative densitometric analysis was performed using ImageJ software.

RNA sequencing of cell spheres and organoids

Cell spheres and organoids were collected in duplicate for RNA extraction using TRIzol (Sigma). cDNA libraries were constructed through PCR amplification and subsequently sequenced on the Illumina HiSeq X10 platform (Annoroad Genomics, China). Adapter-contaminated reads and low-quality reads from the initial sequencing data were filtered out, and the resulting high-quality clean reads were aligned to the human reference genome (GRCh38) using HISAT2 (version 2.2.1), a tool designed for hierarchical indexing for spliced alignment of transcripts. All RNA sequencing samples were processed and analyzed by Xiuyue Biol (China).

Mass spectrometry analysis of cell spheres and organoids

Cell spheres and organoids samples were dispatched to the Experimental Research Center, China Academy of Chinese Medical Sciences, where they underwent protein hydrolysis and subsequent peptide separation using FASP method. These peptides were then analyzed using liquid chromatography coupled with tandem mass spectrometry (LC-MS/MS) on a nanoflow HPLC instrument (EASY-nLC 1000 system, Thermo, USA) linked on-line to a Q-Exactive™ mass spectrometer equipped with a nano electrospray ion source (Thermo). The raw data were processed through the DIA-NN 1.8.1 proteomics platform for Data-independent acquisition (DIA) data searching. The fragmentation spectra were searched with the precursor and fragment mass tolerances set at 10 ppm and 20 mmu for Q-Exactive (QE) data, respectively. Peptide ions were filtered with a cut-off score of 20 and the Percolator algorithm, which employed a P-value threshold of < 0.01. The false discovery rate (FDR) was maintained at 1% for peptide identifications.

Bioinformatics

Gene expression levels were quantified using the fragments per kilobase per million mapped fragments (FPKM) method. Differentially expressed genes (DEGs) were identified with DESeq2, employing a log2 fold change threshold of > 1 and an adjusted P-value cutoff of < 0.05.

Label-free quantitative analyses of the proteome were conducted using the peak area method (Peak Intensity). Precursor ion areas were extracted via DIA-NN1.8.1, with a mass precision of 10 ppm, ensuring accurate recording of experimental m/z and retention times for precursor area quantification. Differentially expressed proteins (DEPs) were selected using the “DEP” package (v3.8), with criteria set at a log2 fold change > 1 and an adjusted P-value < 0.05.

Enrichment analyses, including GO, KEGG, GSEA, and GSVA, were performed to identify significantly enriched terms associated with the DEGs and DEPs. Visualization of heatmaps and volcano maps was achieved through the online tool hiplot (https://hiplot.com.cn/). The Protein-to-Transcript Ratio (PTR) was utilized to assess post-transcriptional regulation of protein abundance, calculated as follows: PTR = Protein(log2LFQ) / mRNA(log2TPM) [15]. The “decoupleR” package (v2.1) was employed to infer transcription factor activity from transcriptomic and proteomic data [16].

Organoid viability assay

To assess the role of calcitriol in hepatocyte differentiation, a concentration gradient was first established. Following this, 5μM calcitriol (MCE, USA) was incorporated into the induction medium from FG to HB cells. The viability of the organoids was evaluated using the Organoid Viability Detection Kit (AIMINGMED, China), adhering to the manufacturer’s instructions. The absorbance of each group was quantified at 560 nm to determine the variations in organoid viability.

Periodic acid-Schiff (PAS) staining of HBs

To evaluate the glycogen storage function of HBs, we conducted periodic acid-Schiff (PAS) staining on paraffin sections of the HBs, following the guidelines provided with the PAS (Beyotime, China) staining kit. The stained sections were then captured using a microscope, and the areas of positive staining were analyzed and compared using Image J software.

Quantification of lactate levels

The LA Assay Kit (Solarbio, China) was used to determine the content of lactate in organoid lysates according to the manufacturer’s instructions.

Seahorse extracellular acidification rate (ECAR) and oxygen consumption rate (OCR) analysis

ECAR and OCR were measured using a Seahorse XFe24 analyzer (Agilent, USA) following protocols to profile the bioenergetics of organoids using Seahorse [17]. Glycolysis rates were measured in response to sequential injections of glucose (10mM), oligomycin (5𝜇M) and 2-DG (100mM) (Agilent) and respiratory rates were measured in response to sequential injections of oligomycin (5𝜇M), FCCP (2𝜇M) and Rotenone/Antimycin-A (2𝜇M) (Agilent). All measurements were normalized to cell number.

Statistical analysis

Statistical analysis was conducted using GraphPad Prism 8.0.2 (GraphPad Software, USA). Quantitative data were first evaluated for normal distribution and reported as the mean ± standard error, followed by a variance homogeneity test. The t-test or one-way ANOVA was employed to ascertain statistically significant differences between two or among multiple homogeneous groups. For significant ANOVA outcomes, Tukey’s multiple comparison test and Dunnett’s multiple comparison tests were applied to analyze the differences between these groups. Results with P values less than 0.05 were considered statistically significant.

Results

Protein and mRNA expression patterns in differentiating HOs

Figure 1A shows the differentiation stages and sequencing patterns of HOs, along with bright-field images. Samples were taken on day 6 (FG), day 9 (HB), day 19 (HO1), and day 29 (HO2) for transcriptomic and proteomic analysis (n = 3). Immunofluorescence indicated liver-specific proteins (ALB) and progenitor markers (KRT19, EPCAM) at various stages (Fig. 1B). RT-qPCR showed a decrease in pluripotency markers (NANOG, POU5F1) and peak expression of foregut-specific CDX2 at FG stage. KRT19 peaked at HB stage and declined, while ALB and AFP increased in mature HOs (Fig. 1C). RNA-seq and mass-spec data confirmed these trends, with increasing liver cell markers (ALB, AFP, APOA4, APOB, HNF4A, CYP3A5) and biliary/progenitor markers (KRT19, TBX3, EPCAM) and decreasing pluripotency markers (NANOG, POU5F1) (Fig. 1D).

Fig. 1

Changes in marker expression during hepatic organoid induction. A.)Experimental design schematic (created by Figdraw). Gene expression at the four stages was analyzed by RNA-seq and Mass-spec. B. Expression of marker proteins (green) during differentiation, Scale bar = 50 μm. C. mRNA expression of marker genes during differentiation. D. Abundance of RNA (orange) and protein (blue) for liver markers, biliary/progenitor markers, and pluripotency markers at different developmental stages. RNA abundance is represented as Fragments Per Kilobase per Million (FPKM), and protein abundance is represented as log2 intensity

RNA-protein multiomics analysis of differentiating HOs

We performed whole transcriptomic and proteomic analyses to study gene regulation during liver differentiation, identifying 7,364 proteins via mass spectrometry (Fig. S1A). PCA showed distinct clustering of RNA and protein sequencing results for each group (Fig. 2A). The correlation matrix and hierarchical clustering heatmap indicated that FG and HB stages have similar gene expression patterns, while HO1 and HO2 are more closely aligned (Fig. 2B and C). Correlation analysis showed that starting from the HB stage, HOs have a strong correlation with primary liver tissue (Fig. S1B). Liver metabolic markers suggest that HOs closely resemble the fetal liver phenotype (Fig. S2). A volcano plot highlighted differential genes and proteins at each stage (Fig. 2D). The distribution of differential genes (DEGs) and proteins (DEPs) varied: FG and HB had similar patterns with mostly upregulation, while HO1 vs. HB and HO2 vs. HO1 showed significant disparities, indicating varying transcriptome and proteome correlations at different stages.

Fig. 2

RNA-Protein multiomics analysis of differentiating HOs. A. Principal component analysis (PCA) of transcriptome and proteome of differentiating HOs. B. The correlation matrix for each stage. C. The hierarchal clustering heatmap of the results of transcriptome and proteome. D. Volcano plots of differential genes (DEGs) and proteins (DEPs) derived from comparison matrices constructed for adjacent stages. Blue dots show significantly upregulated genes or proteins and red dots show significantly upregulated genes or proteins. The colored dots have P < 0.05 and fold change > 2 (|log2 fc|>1). Black dots show detected non-significant (ns) differently expressed genes or proteins

Enrichment analyses for DEGs and DEPs

Gene Set Enrichment Analysis (GSEA) showed that DEGs and DEPs in HO samples from adjacent stages were enriched in specific biological pathways (Fig. 3A). Between HB and FG, up-regulated genes were linked to epithelial differentiation, while down-regulated genes were associated with DNA transcription. KEGG analysis revealed that up-regulated DEGs in the HB stage were enriched in liver-specific pathways like ‘Bile secretion’ and ‘Drug metabolism-cytochrome P450’ (Fig. S3A).

Fig. 3

Enrichment and correlation analyses for differential genes (DEGs) and proteins (DEPs). A. GSEA enrichment analysis results of DEGs and DEPs in HOs at adjacent developmental stages. B-D. Correlation analysis of DEGs and DEPs in each comparison matrix

In the transition from HB to HO1, up-regulated genes were enriched in TNF signaling, apoptosis, and p53 pathways, with a downregulation in cell cycle and lipid metabolism pathways. Conversely, HO2 compared to HO1 showed increased activity in lipid metabolism and PPAR pathways and enhanced liver-specific functions. HO2 organoids also exhibited pathways related to cell morphogenesis, highlighting the importance of the extracellular matrix in organoid development and maturation.

Correlation analysis between DEGs and DEPs

We analyzed the correlation between DEGs and DEPs for 7,541 proteins with corresponding mRNA data to identify genes with matching or differing mRNA and protein expression patterns. The correlation was strongest between HB and FG stages (Fig. 3B) and weakest between HO2 and HO1 stages (Fig. 3D). Functional enrichment analysis showed that in the HB vs. FG comparison, upregulated KEGG pathways included ECM-receptor interaction, indicating that Matrigel addition affected both ECM protein levels and transcription. Downregulated GO pathways were related to foregut-specific gene expression, such as digestive tract development (Fig. 3B). In the HO1 vs. HB comparison, certain GO pathways were downregulated at the mRNA level but upregulated at the protein level, suggesting post-transcriptional or post-translational regulation (Fig. 3C). For the HO2 vs. HO1 comparison, downregulated GO pathways included “extracellular matrix organization,” indicating that Matrigel suppresses endogenous ECM organization in HO2 (Fig. 3D).

Protein-to-mRNA ratio (PTR) analysis at each stage

In the central dogma, the steady-state abundance of proteins is governed by four key rates: transcription and translation, which are integral to the synthesis pathway, and dilution and degradation, which are part of the decay pathway. This process is inherently dynamic. The Protein-to-mRNA ratio (PTR) is utilized to gauge the efficiency of transcription and translation, as well as the stability of mRNA and protein [18, 19] (Fig. 4A). PTR values were classified into high (above the 75th percentile) and low (below the 25th percentile) categories (Fig. 4B, Fig. S3B). Low PTR genes showed a conserved profile across stages, enriched in pathways like “Spliceosome,” “Ribosome,” and “Proteasome,” indicating active transcription, splicing, and protein decay processes (Fig. 4G, Fig. S4). This indicates active transcription and spliceosome processing pathways for these genes, a significant reserve of RNA in ribosomal proteins, or the existence of a conserved protein decay pathway contributing to their degradation.

Fig. 4

Protein-to-mRNA Ratio (PTR) analysis. A. Definition and calculation method of protein-to-mRNA ratio (PTR). B. Scatter plots of PTR at the stage of FG. C-F. GO enrichment analysis for high/low PTR. G. Enrichment analysis of protein interaction networks of low PTR proteins at the stage of FG

In contrast, genes with high PTR exhibited stage-dependent developmental characteristics, with relatively higher protein stability ensuring stage-specific cellular functions: During the FG stage, enrichment was observed in various substance metabolism pathways (Fig. 4C, Fig. S3C); the HB stage featured enrichment in both substance metabolism and RNA processing pathways, highlighting the roles of both protein stability and RNA decay in high PTR outcomes (Fig. 4D); the HO1 stage showed enrichment in several cell cycle progression pathways (Fig. 4E,, Fig. S3C); and in the HO2 stage, high PTR genes were enriched in various RNA processing pathways, suggesting that RNA decay might predominantly drive high PTR in this stage (Fig. 4F, Fig. S3C).

Calcitriol promotes the fate transition from FG To HB stage

To identify key signals regulating the transition from the foregut (FG) to hepatic progenitor (HB) stage, we analyzed transcription factor (TF) activity using the decoupleR package. TFs activated or repressed in HB compared to FG are shown in Fig. 5A. We performed KEGG pathway enrichment analysis on the activated transcription factors (TFs). Interestingly, the enrichment results from both the transcriptome and proteome data indicated the “Parathyroid hormone synthesis, secretion, and action” pathway (Fig. 5B). Analysis of DEGs and DEPs at various stages suggested that the activation of this pathway is specific to the transition from the FG to HB stage (Fig. S5A). The most critical exogenous regulatory signals identified was 1,25-dihydroxyvitamin D, also known as calcitriol. Calcitriol is the most active metabolite of vitamin D and an agonist for the VDR, which is a powerful TF. The expression of VDR transitions from absent to present during the FG to HB stages (Fig. S5B). Our data indicates that several TF activated during the HB stage share target genes with VDR (Fig. 5C).

Fig. 5

Calcitriol promotes the fate transition from foregut (FG) to hepatoblast (HB) stage. A. Transcription factor activity analysis at the HB stage compared to the FG stage. B. KEGG enrichment analysis of activated transcription factors at the HB stage compared to the FG stage. C. Activated transcription factors in the transcriptome and proteome co-regulate several target genes with the vitamin D receptor VDR. D. PAS staining of HB with calcitriol intervention. E. Calcitriol intervention alters mRNA expression in HB. F, G. Calcitriol intervention changes protein expression in HB. Full-length blots/gels are presented in Supplemental Materials_3. Data are presented as the mean ± SD of three replicates. *P<0.05, **P<0.01, ***P<0.001

We hypothesized that calcitriol enhances HB transition by activating VDR. Adding calcitriol to the culture medium during FG to HB induction confirmed this. The addition of calcitriol did not affect the viability of the organoids (Fig. S5C). It increased glycogen storage (Fig. 5D) and upregulated liver markers ALB and HNF4A, while decreasing AFP and KRT19, indicating maturation of HB cells (Fig. 5E). CYP2E1, which promotes hepatocyte apoptosis via oxidative stress, was significantly downregulated [20]. Calcitriol boosted VDR expression. Transcriptomic and proteomic analyses revealed upregulation of HB-stage TFs HNF1A, CEBPA, RXRA, RXRB, CEBPG, and FOXA2. Western blot analysis (Fig. 5F, G) supported these findings, showing elevated ALB and progenitor cell marker EPCAM [21], while KRT19 decreased. Transcription factors CEBPA and HNF1A, which regulate hepatocyte differentiation, were also upregulated.

HO2 is characterized by a predominant glycolytic energy metabolism

We used the scFEA package to infer metabolic flux from transcriptomic data, identifying significant metabolic changes in the HO2 stage (Fig. S6A). The glycolytic flux increased markedly in HO2, with elevated levels of glycolytic intermediates and lactate production (Fig. 6A). GSVA and GSEA analyses confirmed the upregulation of glycolysis and a decline in oxidative phosphorylation (OXPHOS) during HO2 (Fig. 6B, C; Fig. S6B-E). This was validated by performing Seahorse bioenergetic analysis in HO1 and HO2 (Fig. 6D, E; Fig. S6F, G) and measuring lactate levels (a byproduct of glycolysis) (Fig. 6F). The glycolytic pathway activation suggests rapid energy production, with enrichment in the pentose phosphate pathway for nucleotide biosynthesis in HO2 (Fig. S6H). The glycolysis inhibitor 2-DG can significantly inhibit the growth and differentiation of HO2 (Fig. 6G).

Fig. 6

HO2 is characterized by a predominant glycolytic energy metabolism. A. Changes in metabolic flux between the HO2 and HO1 phases in glycolysis and the tricarboxylic acid cycle pathways. Red arrows indicate upregulation, while blue arrows indicate downregulation. B, C. GSEA enrichment results for glycolysis/gluconeogenesis and oxidative phosphorylation pathways in HO2 compared to HO1. D, E. ECAR measurement of a glycolysis stress test of HO1 and HO2. F, Lactate levels in the HO1 and HO2 phases, respectively. G. Light microscopy images of HO2 after glycolysis inhibition with 2μM 2-DG, and mRNA expression of liver markers. H. GSEA enrichment results for the HIF-1 pathway in HO2 compared to HO1. I. Cluster heatmap of transcription factor activity in HO2 compared to HO1. J, K. Downstream genes activated by HIF1 in the transcriptome and proteome of HO2 compared to HO1. *P<0.05, **P<0.01, ***P<0.001, ****P<0.0001

To identify pathways governing glycolysis activation in HO2, we performed GSEA analysis, revealing significant upregulation of the HIF-1 signaling pathway, which regulates glycolysis (Fig. 6H, Fig. S6I). We also focused on the MAPK and PI3K-AKT pathways that regulate HIF-1α translation. Using the decoupleR package, we observed activation of these pathways at both transcriptome (Fig. S6J) and proteome (Fig. S6K) levels, though with some differences. HIF1A, a key transcription factor, was significantly upregulated in HO2, leading to increased expression of glycolytic enzymes such as ENO1, ENO2, HK2, and PFK2 (Fig. 6I-K).

Discussion

We carried out exhaustive sequencing and methodical analysis of both the transcriptome and proteome throughout the entire induction timeline of HOs. This strategy revealed DEGs and DEPs at successive stages of HOs development, along with their corresponding enrichment profiles. Correlational studies between DEGs and DEPs showed a decreasing connection as the induction progressed and as HOs gained functional complexity. PTR assessments indicated that proteins with elevated PTRs demonstrated stage-specific enrichments, while those with lower PTRs exhibited evolutionary conservation. Our research highlighted the pivotal role of calcitriol, a vitamin D metabolite, in facilitating the developmental transition from the FG to HB stage. The inclusion of calcitriol in the culture medium resulted in enhanced glycogen accumulation within HB stage organoids and regulated the expression of critical liver-specific genes, signaling a progression towards liver maturation. Metabolic flux analysis underscored a significant shift towards glycolytic energy metabolism during the HO2 phase, evidenced by increased glycolytic activity and reduced levels of oxidative phosphorylation intermediates. It is important to note that the transcriptional, translational, and metabolic changes occurring during organ fate determination and differentiation have unique characteristics. Although induction protocols vary, they commonly use pathways activated or inhibited at critical stages with cytokines or inhibitors. Our findings, based on a specific protocol, align with the universal liver development timeline and can serve as references for other HOs induction protocols or in vivo liver studies.

Although transcript abundance is often seen as an indicator of protein levels, significant discrepancies frequently exist between mRNA and protein expressions. This divergence is mainly due to various RNA and protein modification mechanisms active throughout transcription and translation, from initiation to termination. For many genes, the correlation between mRNA and the resulting protein is relatively low, typically ranging from 0.21 to 0.37 [18, 22]. Our correlation analysis between DEGs and DEPs yielded similar findings, with coefficients ranging from 0.2 to 0.33 across different groups. These coefficients tend to decrease as differentiation progresses, supporting Liu et al.‘s perspective that steady-state protein levels are mainly determined by mRNA levels but with a delay in synthesis during state transitions [23]. High PTR proteins exhibit stage-specific functions during cellular differentiation, tailored to each developmental stage. Conversely, low PTR proteins show conservation across various stages of differentiation. Enrichment of pathways such as spliceosomes, proteasomes, and ribosomes in low PTR proteins suggests that translation processes and post-translational modifications either inhibit their translation or lead to rapid degradation. This aligns with findings by D. Inverso and colleagues [15]. Notably, proteins essential for the oxidative phosphorylation pathway consistently show low PTR across all stages, indicating tight regulation of energy metabolism through stringent transcriptional and translational controls. These modifications within the oxidative phosphorylation pathway are highly responsive to the cell’s varying physiological states. For instance, during viral infection, specific mRNAs involved in oxidative phosphorylation are translationally upregulated to facilitate viral replication [24].

In the central dogma, each step from DNA to RNA to protein involves modifications that cause variations in RNA and protein expression levels. These modifications add complexity to key genes and proteins during organ development. For organoids, changes in cellular functions should be tracked at both the RNA and protein levels. When critical genes show discrepancies between RNA and protein expressions, post-transcriptional and post-translational modifications must be examined. Two key differences between organoids and in vivo conditions should be considered. First, organoids lack interaction between different organs. Second, the absence of a circulatory system leads to the accumulation of metabolic products and creates an oxygen gradient from the exterior to the interior of spherical organoids. This modified environment is a source of post-transcriptional and post-translational modifications. To develop more accurate organoid models, efforts should be made to simulate the in vivo environment as closely as possible, such as by co-culturing with multiple organs or establishing circulation on microfluidic chips. This approach will help reduce discrepancies with the in vivo environment and enhance model fidelity.

HOs derived from PSCs manifest fetal liver phenotypes and encounter challenges related to their immaturity, which presents obstacles in more faithfully replicating disease models [25]. The endeavor to facilitate the progression of HOs towards adult liver phenotypes and enhance their functional maturation remains a focal point of current research. Epidemiological research has indicated that maternal vitamin D deficiency elevates the risk of childhood obesity in later life [26, 27]. Experimental data from C57Bl/6J mice support these findings, showing that prenatal vitamin D deficiency alters the liver’s cellular composition and downregulates several metabolic pathway genes, including those involved in lipid catabolism, xenobiotic metabolism, and linoleic acid metabolism [28]. Building upon these foundational insights, Fang Y et al. have shown in 2D-cultured iPSC-derived hepatocytes that calcitriol supplementation enhances the functional maturation of these cells [29]. Our research extended this approach by incorporating calcitriol into the induction medium for FG to HB, confirming that calcitriol not only augments the hepatocyte markers within organoids but also diminishes the bile duct markers, concurrently advancing the functional maturation of hepatocytes. A more adult-like HOs model is fundamental for successfully simulating adult liver diseases and conducting drug screening experiments. We propose a maturation protocol that includes the addition of calcitriol, which yields functionally more mature HOs. We hypothesize that these advanced HOs are better suited for modeling adult liver diseases and could serve as regenerative medicine replacements for liver dysfunction caused by alcohol, drugs, or viruses.

The modulation of hepatocyte fate transformation by calcitriol is likely mediated through the VDR, which regulates downstream target genes and interacts with other TF. Nurminen V, et al. have demonstrated that calcitriol modulates the partial activity of the pioneer transcription factor CEBPA [30], which is known to regulate liver self-renewal during liver regeneration [31]. We have additionally noted that calcitriol exerts a similar regulatory effect on the expression of the pivotal transcription factor HNF1A in HOs, which is crucial for hepatocyte differentiation. The precise molecular interplay between the VDR and HNF1A remains to be fully understood, highlighting an area for further investigation into the mechanisms governing hepatocyte maturation within organoid cultures. It is important to note that the VDR does not directly determine liver fate. Instead, it regulates this process at a higher hierarchical level by activating transcription factors associated with liver fate. Consequently, calcitriol is likely to affect non-hepatic cells as well, which is consistent with our experimental results showing decreased KRT19 expression. A more precise cell-specific mapping following calcitriol intervention requires single-cell level analysis, and we look forward to achieving spatial and single-cell resolution in future studies.

Mature HOs predominantly display a glycolytic energy metabolism, a state influenced by processes such as cell replication, differentiation, and self-assembly, which collectively require rapid energy metabolism and material synthesis. This metabolic adaptation is likely reflective of the physiological state of the fetal liver during development. Research on human fetal liver tissues indicates that after 20 weeks of gestation, there is a significant increase in the activity of key enzymes linked to both glycolysis and gluconeogenesis, with glycolysis becoming more dominant as birth approaches [32]. Metabolomic studies of mouse livers at various stages of gestation show that the fetal liver maintains a high glycolytic activity until about embryonic day 10–15 (E10-15), which then notably decreases by E15-18 [33]. We hypothesize that the stage-specific predominance of glycolysis arises from two main factors: firstly, the energy demands of rapidly dividing cells; and secondly, the accumulation of metabolic byproducts (such as lactate), which serve as sources for protein modifications. These byproducts create a feedback loop with gene expression, regulating the transcription, translation, and function of key genes involved in developmental processes. For instance, Nozomu Takata and colleagues found that lactate produced through aerobic glycolysis mediates the upregulation of key eye-field transcription factors by epigenetically modifying histone H3 acetylation, thereby regulating mammalian eye morphogenesis [34]. The difficulty in identifying anatomical structures in early embryos restricts our comprehension of the individual organ’s energy metabolism state during the initial phases of embryonic development. Furthermore, the liver’s adaptation to the energy metabolism needs of other organs throughout embryogenesis introduces an additional level of intricacy to the examination of liver metabolism. HOs possess the potential to serve as a crucial model for investigating energy metabolism at the early stages of organogenesis, either in isolation or in synergy with other embryonic-like organoids through co-culture. Additionally, by integrating protocols that promote the maturation of HOs, it is possible to construct HOs that more closely resemble adult liver tissue. When combined with organoids from other metabolic organs, complex multi-organoid systems can be established to simulate various metabolic phenotypes and intricate in vivo metabolic disease states, such as metabolic syndrome. In regenerative medicine, maintaining the glycolytic activity of HOs might offer a more functionally optimal transplantation strategy for liver replacement.

Although our research provides significant insights, it is important to acknowledge its limitations. Notably, we did not conduct single-cell resolution analyses of the transcriptomic and proteomic changes occurring during the HOs induction process. Furthermore, future research could explore the effects of calcitriol in diverse HOs disease models and optimize organoid culture conditions to more accurately replicate in vivo environments, offering potential avenues for advancing our understanding and applications in this field.

Conclusions

In summary, our study has elucidated the pivotal role of glycolytic metabolism during the HO2 stage and the regulatory influence of calcitriol in facilitating the transition from the FG to HB stage. These insights emphasize the critical involvement of metabolic and signaling pathways in the development of HOs, with profound implications for both fundamental research and clinical applications. The prospect of leveraging these pathways for therapeutic interventions underscores the transformative potential of our research in the realms of stem cell biology and regenerative medicine. As we delve deeper into the intricate interplay between metabolism, signaling, and cell fate, our findings lay the groundwork for pioneering approaches to manage liver diseases and enhance organoid-based therapeutic strategies.

Data availability

The mass spectrometry proteomics data have been deposited to the iProX [35, 36] with the dataset identified PXD054162. The gene expression data of RNA-seq data are deposited to the GEO with the dataset identified GSE264422.