Integrative analysis of single-cell RNA-seq and ATAC-seq reveals heterogeneity of induced pluripotent stem cell-derived hepatic organoids

Abstract

To gain deeper insights into transcriptomes and epigenomes of organoids, liver organoids from two states (expandable and more differentiated) were subjected to single-cell RNA-seq (scRNA-seq) and single-cell ATAC-seq (scATAC-seq) analyses. Mitochondrial gene expression was higher in differentiated than in non-differentiated hepatocytes, with ATAC-seq peaks increasing near the mitochondrial control region. Differentiation of liver organoids resulted in the expression of transcription factors that act as enhancers and repressors. In addition, epigenetic mechanisms regulating the expression of alpha-fetoprotein (AFP) and albumin (ALB) differed in liver organoids and adult liver. Knockdown of PDX1, an essential transcription factor for pancreas development, led to the hepatic maturation of liver organoids through regulation of AFP and ALB expression. This integrative analysis of the transcriptomes and epigenomes of liver organoids at the single-cell level may contribute to a better understanding of the regulatory networks during liver development and the further development of mature in vitro human liver models.

Keywords: Experimental models in systems biology; Integrative aspects of cell biology; Molecular mechanism of gene regulation; Transcriptomics.

Graphical abstract

Figure 1

scRNA sequencing analysis of liver organoids in HM and DM (A) UMAP derived from scRNA-seq analysis of 17,434 liver organoid cells in HM and 21,866 cells in DM (top). Each point indicates a single cell within each cluster. The profile of each organoid in UMAP (bottom). (B) Experimental design to treat liver organoids with HM and DM (top). Representative images of liver organoids in HM and DM (bottom). (C) Quantification of albumin secretion in each replicate sample. Data are the mean ± SEM (n = 3) compared by Student t tests (∗p < 0.05). (D) Genes expressed in each cluster. (E) Representative immunofluorescence images of organoids cultured in DM and stained with each indicated antibody. (F) Cell trajectory analysis using Slingshot. Each point represents a cluster, with the arrows indicating the inferred direction of differentiation based on trajectory analysis. The pie chart shows the proportion of HM- and DM-associated liver organoids. (G) Gene set enrichment analysis for the inter-DEGs of each cluster using EnrichR.

Figure 2

Integrated analysis of scRNA-seq and scATAC-seq results of liver organoids treated with HM and DM (A) Integration of scATAC-seq with scRNA-seq. tSNE maps of liver organoids in HM and DM: (left) before integration; (center) integration by batch correction; (right) cluster-based clusters predicted after integration. Cells with low-prediction scores (<0.3) were removed using Snaptools. Each point indicates a single cell colored by integrative analysis. (B) Heatmap of normalized motif enrichment scores for each cluster and TF. Because the number of OCRs differed among clusters, the p value of motif enrichment was normalized relative to the number of OCRs of each cluster. (C) Interaction scores were divided into four categories: high (p < 10⁻¹⁰, N = 30,409), middle (10⁻¹⁰ < p < 10⁻³, N = 28,251), low (10⁻³ < p < 0.05, N = 23,182), and unrelated (p > 0.05, N = 33,879). The degree of overlap was compared with the peaks of TFs from ENCODE ChIP-seq data. The arrows indicate enhancer-promoter and TAD loop-associated TFs (i.e., gray: RAD21; red: CTCF; blue: YY1). (D) The percentage of OCRs corresponding to each category. (E) The OCRs of scATAC-seq were divided into cluster-specific and common regions. The inter-DAR selected from a single cluster was set as the specific region, with the remaining inter-DARs set as the common region. Fold changes were calculated for specific and common regions in the center of the peaks of TF ChIP-seq from ENCODE. Arrows indicate three TFs (RAD21, CTCF, and YY1) shown in Figure 2C. (F) The percent overlap between common and specific regions is shown on the boundary of intra-TAD, inter-TAD, and Hi-C-based TAD.

Figure 3

Integrated analysis of inter-DEG and inter-DAR revealing cluster-specific transcriptional activators (A and B) Integrated analysis of inter-DEGs and inter-DARs showing cluster-specific transcriptional activators. (B) Enhancer pairs were defined by three criteria: OCRs with interaction scores >10 (left); positive correlation between inter-DEG and inter-DAR in the total clusters (center); and positive correlation (p < 10⁻²⁰) with TF located in OCR associated with a cluster-specific gene in the scRNA-seq profile (right). (C) Heatmap summarizing the enhancer pairs according to log₂(fold changes) of scRNA-seq and scATAC-seq. (D) Proportions of enhancer pairs of representative TFs for each cluster. (E) Illustration of the EPCAM gene in a genomic browser for scATAC-seq levels for each cluster. (F) Interaction scores (top) and associated TFs (bottom) for each enhancer of the EPCAM gene.

Figure 4

OCRs in the mitochondrial genome (A) Illustration of OCRs of mitochondrial genomes in liver organoids, showing the control region (light green, inner circle) and the HNF4a motif (yellow, inner circle) located in MR1. (B) Genome browser showing the OCRs of each cluster in the scATAC-seq and bulk liver ATAC-seq data from ENCODE in the mitochondrial genome. (C and D) Log₂(fold change) in inter-DEG and inter-DAR for each cluster. (E) Motifs in the MR1 and MR2 regions. Motifs were analyzed after dividing MR1 and MR2 by 400 bp. (F) Representative FACS analysis and morphology of liver organoids in DM stained with MitoTracker (left). mRNA expression levels of mitochondrial-encoded genes and mature hepatocyte markers in low- and high-MitoTracker populations and PHHs (right). β-actin was used as an internal control. Data were reported as the mean ± SEM (n = 3) and compared by Student t tests. ∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001.

Figure 5

Enhancers and repressors of each organoid in HM and DM, as shown by the integrated analysis of DEGs and DARs at the bulk level (A and B) Strategy to identify enhancers and repressors differing in liver organoids treated with HM and DM at the bulk level. (B) Gene expression and scATAC-seq peaks were divided into DEGs and DARs of organoids treated with HM and DM at the bulk level, with OCRs associated with a specific gene defined as having an intersection score of p < 0.05 (left). The four patterns of regulation of gene expression in liver organoids treated with HM and DM (right). DM enhancer: both bulk-DEG and bulk-DAR increased in DM; DM repressor: bulk-DEG increased in DM, but bulk-DAR decreased in HM; HM repressor: bulk-DEG increased in HM, but bulk-DAR increased in DM; HM enhancer: both bulk-DEG and bulk-DAR increased in HM. The pair type counts were based on the number of pairs between OCRs related to the genes. (C) Gene-OCR pairs determined by log₂(fold change) of scRNA-seq and scATAC-seq at the bulk level. (D) Representative TFs summarized for each pair type. (E) Genomic browser showing the bulk-level scATAC-seq for OCRs corresponding to enhancers and repressors of the AFP gene (chr4:74,273,152-74,328,858). (F) Log₂(fold change) of the bulk-level DAR of OCRs corresponding to Figure 5E. (G) Interaction score (top) and enhancer and repressor TF (bottom) in the OCRs corresponding to Figure 5E.

Figure 6

Differences in gene expression and OCR of AFP and ALB genes in liver organoids and tissues (A) scRNA-seq data showing the correlation of AFP and ALB expression. The R value is the Pearson correlation coefficient. (B) Selection of adult liver-specific OCRs. OCRs common to the adult liver and liver organoids were excluded. (C) A portion of the overlap of adult liver-specific OCRs with liver TF ChIP-seq data from ENCODE. (D) Genome browser showing liver organoid scATAC-seq, adult liver ATAC-seq, and ENCODE TF ChIP-seq results at chr4: 74,200,000–74,330,000 (top). The scATAC-seq peaks for liver organoids treated with HM and DM are displayed at the bulk level. The overlapping OCRs for liver organoids and adult liver are shown in green boxes and adult-specific OCRs in yellow boxes. The adult liver-specific OCRs near the AFP and ALB genes were named A1 to A4. Genome browser showing publicly available data from the GSE131336 dataset and ENCODE of ATAC-seq data from embryo to adult and CTCF ChIP-seq data in mouse adult liver at chr5: 90,398,000–90,520,000 (bottom). Additionally, the regions found above were marked in the mouse genome by liftover (dashed lines). (E) Total RNA-seq data of ENCODE, showing the expression of the Afp and Alb genes from embryo to adult mice. (F) Rad21, Smc3, and CTCF ChIP-seq data of R6 under control (CTRL) and Nipbl (NIPBL) deletion conditions in adult mouse livers.

Figure 7

PDX1 suppression induced further maturation of liver organoids (A) Experimental scheme of PDX1 knockdown by lentiviral shRNA transduction. (B) Representative morphology and mCherry expression of shNT- and shPDX1-transduced organoids 3 days after transduction. (C) mRNA expression levels of each indicated gene in shNT- and shPDX1-transduced organoids on day 10. β-actin was used as an internal control. (D) Quantification of albumin and alpha-fetoprotein secretion by shNT- and shPDX1-transduced organoids. (E and F) (E) ChIP-PCR and (F) ChIP-qPCR analyses of H3K27ac and PDX1 at the R6 region in shNT- and shPDX1-transduced organoids. (G) ChIP-PCR and (H) ChIP-qPCR analyses of H3K27ac at A1 and A4 regions in shNT- and shPDX1-transduced organoids. The relative enrichment of H3K27ac and PDX1 was normalized to the input DNA. All experimental data are the mean ± SEM (n = 3) and were analyzed using Student’s t test. ∗p < 0.05, ∗∗p < 0.01, and ∗∗∗p < 0.001. (I) Illustration of genetic and epigenetic modulation of the R6 region by PDX1 knockdown in organoids.