Abstract The use of differentiating human induced pluripotent stem cells (hiPSCs) in mini-tissue organoids provides an invaluable resource for regenerative medicine applications, particularly in the field of disease modeling. However, most studies using a kidney organoid model, focused solely on the transcriptomics and did not explore mechanisms of regulating kidney organoids related to metabolic effects and maturational phenotype. Here, we applied metabolomics coupled with transcriptomics to investigate the metabolic dynamics and function during kidney organoid differentiation. Not only did we validate the dominant metabolic alteration from glycolysis to oxidative phosphorylation in the iPSC differentiation process but we also showed that glycine, serine, and threonine metabolism had a regulatory role during kidney organoid formation and lineage maturation. Notably, serine had a role in regulating S-adenosylmethionine (SAM) to facilitate kidney organoid formation by altering DNA methylation. Our data revealed that analysis of metabolic characterization broadens our ability to understand phenotype regulation. The utilization of this comparative omics approach, in studying kidney organoid formation, can aid in deciphering unique knowledge about the biological and physiological processes involved in organoid-based disease modeling or drug screening. Keywords: induced pluripotent stem cells, kidney organoids, transcriptomics, metabolomics, serine metabolism Introduction Personalized induced pluripotent stem cells (iPSCs) is an emerging technology that increases the number of options available for regenerative medicine applications ([47]Robinton and Daley, 2012; [48]Cossu et al., 2018). A major advance in recent years has been the use of stem cell-derived organoids as in vitro models to study development and disease ([49]Lancaster and Knoblich, 2014). The ability of iPSCs to differentiate into various cell types has been utilized recently to generate organoids ([50]Wiegand and Banerjee, 2019). Organoid culture has been established for various organs and tissues, including the intestine, heart, liver, brain, and kidney ([51]Kretzschmar and Clevers, 2016). Unlike zebrafish, which have a single unit of the nephron that can be repaired or regenerated to a certain extent ([52]Diep et al., 2011; [53]Kroeger and Wingert, 2014), the human kidney contains about 1–1.5 million nephrons serving as the functional units within each kidney. Hence, generating complex kidney tissues, such as kidney organoids, in a large quantity from personalized human induced pluripotent stem cells (hiPSCs) have broadened our ability to study human kidney development, disease, and even perform tissue engineering ([54]Xia et al., 2013; [55]Taguchi et al., 2014; [56]Takasato et al., 2014, [57]2015; [58]Freedman et al., 2015; [59]Morizane et al., 2015; [60]Taguchi and Nishinakamura, 2017; [61]Przepiorski et al., 2018; [62]Low et al., 2019). These methods have established a concrete foundation for generating kidney organoids through the modulation of wingless-related integration site (WNT), fibroblast growth factor (FGF), and transforming growth factor β (TGF-β) signaling. Although the mechanisms involved in the different protocols have been well-studied, iPSCs derived from different donors perform differently between similar protocols. Additionally, batch variation was not accounted for, and during the culture process, quality feedback from the derived kidney organoids was not performed adequately. Emerging studies of stem cell metabolism have demonstrated that cellular metabolic pathways influence the proliferation, reprogramming, or differentiation process ([63]Zhang et al., 2012, [64]2018; [65]Xu et al., 2013; [66]Ryall et al., 2015). Monitoring the metabolic profile is an effective way to describe the iPSC status because there are dynamic changes in mitochondrial morphology and metabolic shifts under different culture conditions ([67]Shyh-Chang et al., 2013; [68]Son et al., 2013; [69]Wanet et al., 2015). Moreover, iPSC differentiation to other cell types is generally accompanied by cellular metabolism alterations. For example, a high glycolytic level stimulates embryonic stem cell (ESC) proliferation ([70]Kim et al., 2006), whereas differentiating ESCs switch to oxidative phosphorylation (OXPHOS) ([71]Folmes et al., 2012). The metabolic switch from OXPHOS to glycolysis occurs during the reprogramming process and cooperates with the epithelial-mesenchymal transition to regulate pluripotency ([72]Sun et al., 2020a). In contrast, the bioenergetic profile changes from glycolysis to OXPHOS during cardiomyogenic differentiation ([73]Lopaschuk and Jaswal, 2010), neuronal differentiation ([74]Zheng et al., 2016), and induction of hepatocyte maturation ([75]Hopkinson et al., 2017). Consistently, during the maturation process, nephron progenitor cells (NPCs) have significantly higher glycolysis compared with matured NPCs, and the inhibition of glycolysis promotes nephrogenesis in embryonic kidneys ([76]Liu et al., 2017). Therefore, we hypothesized that characterizing the metabolic profile of the kidney organoids derived from iPSCs might facilitate the discovery of the essential metabolic element. Metabolic profiling can be a more direct method to monitor the quality of organoid differentiation and maturation, and it can also describe the dynamic changes in cellular status during drug screening. Here, in this study, we used comparison studies between transcriptomics and untargeted metabolomics. The identification and correlation of time-dependent metabolic changes in iPSC-derived kidney organoids and the metabolic regulatory mechanism during differentiation were established. Furthermore, we identified that glycine, serine, and threonine metabolism might play an essential role during kidney organoid differentiation. Finally, we verified that during the kidney organoid differentiation, observation of DNA methylation enhancement was mediated through the regulation of serine metabolism-related bioactive metabolite (SAM). Results Generation of Human iPSC-Derived Kidney Organoids To obtain kidney organoids in a scalable manner, we modified a protocol from [77]Przepiorski et al. (2018). The detailed procedure was explained in the methods, and the three stepwise differentiation procedure was illustrated ([78]Figure 1A). Briefly, the assay’s initial setup involved growing hiPSC colonies in a monolayer to about 80% confluence. On day 0 (D0), the colonies were detached with dispase, then scraped and transferred into a 6-well ultra-low attachment plate for suspension culture to form aggregated embryoid bodies (EBs) ([79]Figures 1B,C). Next the formed EBs were cultured in “Stage II” medium from day 3 (D3) until tubule formation which was visible under bright-field microscopy after day 8 (D8), and tubule structuralized organoids at around day 14 (D14) ([80]Figures 1D–F). Subsequently, the quantification of gene expression for lineage markers such as NPHS1 (podocytes), LRP2 (proximal tubule), SLC12A1 (thick ascending limb), CALB1 (collecting duct), GATA3 (collecting duct), and CD31 (endothelial cells) was carried out at different time points using quantitative RT-PCR (qRT-PCR) ([81]Figure 1G). The quantification showed that all these lineage or cell type-specific markers were upregulated from D0 to D14, which recapitulated the developmental program of nephrogenesis. Tubule formation might take a different lineage approach depending on differentiation condition, hence, apart from GATA3 (collecting duct) and LRP2 (proximal tubule) which had biphasic responses, all other markers had continual increases during the assays with the highest expression levels detected at D14. FIGURE 1. [82]FIGURE 1 [83]Open in a new tab Schematic representation and characterization of the kidney organoid differentiation process. (A) Roadmap of kidney organoids derived from hiPSCs. (B,C) Starting iPSC colonies and iPSC colonies digested by dispase. (D) Embryoid bodies. (E,F) The phenotype of kidney organoids at day 8 (D8) and day 14 (D14); formation of tubules is visible on the surface (arrows). (G) qRT-PCR analysis for selected markers during the differentiation at day 0, 3, 8, and 14. Data are presented as mean ± SEM from biological triplicates. **P < 0.01, ***P < 0.001. (H–L) Immunostaining of frozen sections of day 14 organoids showing NPHS1+ podocytes and CD31+ endothelial cells, LTL+ proximal tubules, and CDH1+ distal tubules, LRP2+ proximal tubules, and UMOD+ thick ascending limb segments, GATA3+ collecting duct structures, and MEIS1/2/3+ interstitial cells. Nuclear counterstain: DAPI. (M) Schematic representation of kidney organoids (day 14). (N,O) Transmission electron micrograph (TEM) images of day 14 organoid sections: podocytes (p, podocytes; fp, primitive foot processes); and tubular (te, tubular epithelium; bb, brush border). Furthermore, immunostaining analysis of kidney organoids (D14) with a range of physiology-related markers also confirmed the results observed during nephron differentiation. The identification of NPHS1+ podocytes and CD31+ endothelial cells showed that kidney organoids contained glomerular structures and were surrounded by immature vascular cells ([84]Figure 1H). Proximal segments that stained for Lotus tetragonolobus lectin (LTL) and distal portions labeled with cadherin-1 (CDH1) demonstrated the formation of the proximal and distal tubules ([85]Figure 1I). The UMOD+ thick ascending limb segments were identified with a short segment existed between LRP2+ proximal tubules, which could be a part of the loop of Henle ([86]Figure 1J); GATA3 had positive staining for collecting duct structures ([87]Figure 1K); MEIS1/2/3+ interstitial cells were also identified that were scattered throughout the tissue ([88]Figure 1L). These data collectively displayed that organoids (D14) consist of podocytes in the glomerular structures, proximal tubular, distal tubular, collecting duct structures, and interstitial cells ([89]Figure 1M). The maturation status was further determined by examining the glomerular and tubular structures using transmission electron micrograph (TEM) analysis. TEM results demonstrated characteristic structures, such as podocytes with primitive foot processes in the glomerulus ([90]Figure 1N) and tubular structures with multi-layered epithelial structures and brush-borders in organoids at D14 ([91]Figure 1O). Together, these results confirmed that we established a differentiated kidney organoid culture system consistent with the maturation levels published in the literature. Transcriptional Profiling Reveals the Metabolic Alteration That Occurs During Kidney Organoid Differentiation To understand the gene expression dynamics during the different stages of organoid differentiation, RNA-seq analysis of organoids was conducted at different time-lapse sampling points. A total of 11 samples from four different time points were sequenced for RNA expression (D8 group had only n = 2 biological repeats because of sequence library construction, the other groups had n = 3 repeats). According to Principal Component Analysis (PCA) based on counts per million (CPM) values of genes, replicate samples from similar time points were clustered closely ([92]Figure 2A). These data illustrated that the correlation coefficients of samples were achieved ([93]Supplementary Figure 1). FIGURE 2. [94]FIGURE 2 [95]Open in a new tab Transcription profiling of the kidney organoid differentiation process. (A) Scores plot for transcriptome samples in the PCA analysis. (B) k-Means clustering of expressed genes during different differentiation phases (the top 12,000 variable genes were used). Known markers of kidney development are highlighted. (C) Top 10 overrepresented GO terms in each cluster with functions in biological processes. (D) Top 10 enriched pathways in each cluster based on KEGG pathway enrichment analysis. Pathways with adjusted P-value (adjPval) below 0.01 (2C and 2D). Subsequently, a k-Means clustering analysis was carried out, which divided the top 12,000 of the most variable genes ranked by the standard deviation into four clusters. Interestingly, each cluster had specific marker genes referring to different stages of nephrogenesis. In detail, POU5F1, NANOG, and SOX2 served as pluripotency markers in cluster A and were highly expressed in D0 organoids as expected. Simultaneously, mesoderm markers TBXT, HAND1, and GSC in cluster B were enriched in D3 organoids. Additionally, cluster C had both intermediate mesoderm markers (WT1, LHX1, and SIX2) and mature markers (NPHS1, LRP2, and CD31) enriched in D8 to D14 organoids, while cluster D had additional mature markers (CALB1, AQP2, and MEIS3) enriched in D14 organoids ([96]Figure 2B and [97]Supplementary Figure 2). These findings indicated that kidney organoids progressed through the mathematically-grouped clusters. As indicated, this grouping matched up with the development-related stages, which indicated maturation of renal cell types and correlated with the varied sampling times. Furthermore, Gene Ontology (GO) biological progress annotation was exploited for each cluster gene ([98]Figure 2B), which demonstrated that each cluster gene had a functional role in different biological processes. In particular, genes within cluster A were enriched for DNA replication and cell development progress. In contrast, genes within cluster B were significantly enriched in protein translation-related progression, while cluster C genes were related to the processes involved in differentiation, such as tissue development, tube development, and animal organ morphogenesis ([99]Figure 2C). Additionally, we performed Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis for each cluster gene ([100]Figure 2B) and found that 190 genes in cluster B were strikingly enriched in the metabolic pathway ([101]Figure 2D). These data implied that metabolic gene expression dynamics appeared to have distinct profile stages during the differentiation of kidney organoids, which was consistent with our hypothesis that there might be a correlation between metabolic changes and progression of differentiation. Classification of Metabolic Genes and Changes in Metabolite Dynamics As reported in the literature, cellular metabolism can play crucial roles in maintaining pluripotency, differentiation, and reprogramming of iPSCs ([102]Zhang et al., 2018). A metabolic switch from the glycolysis to OXPHOS state was described consistently in cells that differentiated from pluripotent states initially. Hence, we further investigated the perturbation of the intracellular metabolism stimulated by kidney organoid differentiation. The expression pattern of 190 genes enriched in the metabolic pathway from cluster B ([103]Figure 2D) was divided into three clusters (MGCluster, R Mfuzz package), which included 46, 59, or 85 upregulated genes from D0 to D3, respectively, corresponding to three diverse expression patterns ([104]Figure 3A). MGCluster 1 contained downregulated genes at D8 but upregulated at D14, and genes in MGCluster 2 declined continuously from D3 to D14, while the expression of genes overrepresented in MGCluster 3 remained at a relatively high level from D3 to D14. FIGURE 3. [105]FIGURE 3 [106]Open in a new tab Dynamic changes in metabolite-related genes. (A) Time-course expression profiles of 190 genes enriched in the KEGG term metabolic pathways. (B) KEGG functional analysis of three overrepresented clusters. Pathways with a P-value below 0.05 are shown. (C) Heatmap of relative gene expression related to OXPHOS from three overrepresented clusters. (D) Heatmap of relative gene expression in the TCA cycle related to OXPHOS and in the glycolysis pathway. Genes highlighted in the red frame encoding the rate-limiting enzymes. Heatmap shows averaged values from biological replicates. Subsequently, we performed KEGG pathway enrichment analysis to characterize the function of these three MGClusters ([107]Figure 3B). Consistently, we found that every MGCluster was significantly enriched in OXPHOS, indicating the dynamic transcriptional response of OXPHOS-related genes involved in kidney organoid differentiation, which were associated with the accelerated generation of intracellular energy. This finding confirmed the relative expression level of genes enriched in OXPHOS during the progression of kidney organoid differentiation ([108]Figure 3C). MGCluster 3 contained the major complexes that encoded genes in the OXPHOS. Incorporating genes were upregulated consistently during differentiation and this was the dominant dynamic pattern among the three clusters. We paid particular attention to the expression levels of SDH, IDH2 encoding isocitrate dehydrogenase (which is regarded as a rate-limiting enzyme in the tricarboxylic acid (TCA) cycle), and SUCLA2 encoding succinyl-CoA synthetase beta subunit, because of their essential role in the TCA cycle or OXPHOS. The gene expression level was verified using qRT-PCR via the selected marker IDH2 ([109]Supplementary Figure 3). The verified RNA-seq data were plotted in the heatmap and this demonstrated that the transcriptional level of SUCLA2 and IDH2 consistently improved. This improvement possibly benefited from the provision of more reducing equivalents and an energy substrate for OXPHOS by elevating the carbon flux through the TCA cycle ([110]Figure 3D). Moreover, the altered expression level of SDH could provide alternative evidence for the activation of OXPHOS during the later phases of differentiation. To systematically characterize and evaluate the metabolic changes during transcriptional regulation, we simultaneously performed untargeted metabolomic analysis to identify key metabolic profiles. PCA analysis demonstrated that samples from all four-time points were separated into four distinct groups, as shown in the scores plot ([111]Supplementary Figure 4A), suggesting a recognizable alteration in metabolism during the different stages of differentiation. The heatmap of 96 annotated metabolites among different samples is shown in [112]Supplementary Figure 4B. The change of pool sizes of these metabolites could also be divided into three clusters ([113]Supplementary Figure 4C), and the result of KEGG functional analysis was shown ([114]Supplementary Figure 4D). Accordingly, several metabolic pathways held a consistent alteration pattern on both a transcriptional and metabolic level, these include purine, pyrimidine, alanine, aspartate, and glutamate metabolism ([115]Supplementary Figure 5), which might associate with a high requirement of nucleotides in the early stage of differentiation (from D0 to D3). Collectively, to a certain extent, these data demonstrated that alterations in metabolism measured directly in this study were consistent with the transcriptional levels of enriched genes, which served as the cross-references to verify these findings.