A human multi-lineage hepatic organoid model for liver fibrosis

To investigate the pathogenesis of a congenital form of hepatic fibrosis, human hepatic organoids were engineered to express the most common causative mutation for Autosomal Recessive Polycystic Kidney Disease (ARPKD). Here we show that these hepatic organoids develop the key features of ARPKD liver pathology (abnormal bile ducts and fibrosis) in only 21 days. The ARPKD mutation increases collagen abundance and thick collagen fiber production in hepatic organoids, which mirrors ARPKD liver tissue pathology. Transcriptomic and other analyses indicate that the ARPKD mutation generates cholangiocytes with increased TGFβ pathway activation, which are actively involved stimulating myofibroblasts to form collagen fibers. There is also an expansion of collagen-producing myofibroblasts with markedly increased PDGFRB protein expression and an activated STAT3 signaling pathway. Moreover, the transcriptome of ARPKD organoid myofibroblasts resemble those present in commonly occurring forms of liver fibrosis. PDGFRB pathway involvement was confirmed by the anti-fibrotic effect observed when ARPKD organoids were treated with PDGFRB inhibitors. Besides providing insight into the pathogenesis of congenital (and possibly acquired) forms of liver fibrosis, ARPKD organoids could also be used to test the anti-fibrotic efficacy of potential anti-fibrotic therapies.

confirms that theT36M mutation was introduced into 3 different donor iPS cells. A silent mutation was introduced to fulfill sgRNA recognition PAM sequence. b, The pipeline for the multiplex analysis of the scRNA-Seq data generated from control HOs prepared from 3 unrelated individuals (C1-C3), and from ARPKD organoids that were prepared from isogenic iPSCs (M1, M2, M3) with an engineered PKHD1 M36 mutation. Pooled samples of control and ARPKD organoids were separately analyzed, and the multiplexed data was deconvoluted based upon analysis of allelic differences using 'demuxlet' software 2 . c, Dot plots show a summary the number of cells identified for each subject by 'demuxlet' in control (left) and ARPKD (right) pooled samples. Dots that are close to base line represent doublet or ambiguous droplets that were not analyzed. d, A heat map showing marker genes for the 15 clusters. Each cluster is represented by the top 5 marker genes, which are color-coded for each cluster; and representative genes for each cluster are shown on the left. Each column is an individual cell within the indicated cluster at the top; and each gene is in a row. A Wilcoxon rank sum test was used to test for differential expression of the mRNAs; all markers were expressed in >25% of the cells; and the threshold for selecting a marker was set at a minimum of log2 (fold-change) > 0. 25 cholangiocytes (KRT10, SOX9), myofibroblasts (PDGFRB, COL1A1) and primary cilium (PKDH1, PKD1) during HO development. scRNA-Seq data was generated from iPSC, day 9 hepatoblast (HB), and in primary (HO1) and secondary (HO2) organoid cultures. For comparison purposes, scRNA-Seq from primary human hepatocyte (PHH), and published scRNA-Seq data for cholangiocytes (Cho) and mesenchymal cells (MSC) in human liver tissues was used for this analysis. The secondary organoid (HO2) cultures are formed by dissociating HO1 organoids into single cells, which will then reform organoids in 12 days. Each dot shows the level of expression of the indicated mRNA at the indicated stage as determined by analysis of the scRNA-Seq data. PKDH1 and PKD1 mRNAs are expressed at a lower level than other mRNAs in the organoid cultures; and because of this, their expression levels in human hepatocytes and in cholangiocytes in human liver are shown as separate graphs. Of note, the primary cilium mRNAs are expressed at the hepatoblast and HO stages in the organoid cultures, which means they are expressed during the period when hepatoblasts differentiate into hepatocytes and cholangiocytes. i, Cells in day 9 control and ARPKD organoid cultures have equivalent levels of expression of mRNAs encoding myofibroblast markers. RNA was prepared from day 9 (hepatoblast stage) differentiating control and ARPKD organoid cultures (n=3 donors, each measured in duplicate). The box plots show the following: center line, median; box limits, upper and lower quartiles; whiskers, 1.5interquartile range. RT-PCR was used to measure the level of expression of mRNAs encoding 5 genes that are expressed in myofibroblasts, and the level of PKDH1 mRNA expression is also shown. Each dot shows the measurement obtained from one culture; the vertical line is the average expression level; and the box plot covers the 25 to 75 percentiles. Statistical differences between the groups were determined using the unpaired two-tailed t-test. 1.5interquartile range. Statistical differences between the groups were determined using the unpaired two-tailed t-test. N.S. indicates that the p > 0.05. Figure 9. A proposed model for the pathogenesis of ARPKD liver disease. Within a normal cholangiocyte (dotted box), the FPC-mediated interaction between the TGF receptor and the E3-ubiquitin complex will lead to TGF receptor (TGFR) degradation. The ARPKD mutation (FPC*) interferes with this interaction, which inhibits TGFR degradation in ARPKD cholangiocytes. ARPKD cholangiocytes also have an increased level of MMP-2 and MMP-9 expression, which converts latent TGF to its active form. Increased TGF and TGFR expression activates TGF-associated signaling pathways in ARPKD cholangiocytes. TGF acts in concert with LIF and PDGF, which are constitutively produced by multiple different cells types, to promote mesenchymal cell differentiation into collagen-producing myofibroblasts. The myofibroblasts, which have increased levels of LIF receptor and PDGFR receptor expression and STAT3 pathway activation, mediate the pathogenesis of ARPKD liver fibrosis.  4 . The comparisons were performed using the Seurat label transfer function 7 . Based upon the highest level of concordance between the clusters and the sequences present in the previously defined cell types identified in human liver: Clusters 0, 1, 3, 4, 6 and 7 are identified as scar-associated mesenchymal cells (SAMe); cluster 2 as mesothelia; cluster 10 as cholangiocytes; and cluster 14 as endothelia. Since the reference sequences only include a very limited number of hepatocytes, cell clusters containing hepatocytes could not be identified by this analysis. However, because mRNAs (ALB, TTR, HNF4A and SERPINA1) encoding known hepatocyte markers were expressed in cluster 5, they were identified as hepatocyte precursor cells (see Fig. S4)). Although a high proportion of cluster 8 and 9 cells were labeled as endothelial cells by the Seurat analysis, because the Pearson analysis indicates that their correlation with reference endothelial cells is ~0.5, they were identified as early endothelial cells. The identities of clusters 11, 12 and 13 could not be assigned with certainty. Numbers within the parentheses located next to cell type indicators correspond with the subpopulations annotated in 4 .  Figs. 1b, c). In addition to the cells with canonical hepatocyte and cholangiocyte markers, cells expressing endothelial cell (PECAM1, ICAM1 and KDR) and

Gene
fibroblast (PDGFRB, COL1A1, DES and ACTA2) markers were present. (Supplementary Fig.   1e-f). To plot the progression of iPSC to hepatic lineage differentiation, the "Monocle" package 8 was also used to construct single sell trajectories. It drew a branched single-cell trajectory beginning with iPSC, with a transition at the hepatoblast stage, and it terminated at the PHH and other cell types present in HO (Supplementary Fig. 1h). A CyTOF analysis was performed to characterize the expression of multiple protein antigens on the cells at different differentiation stages within organoid cultures. The resulting bhSNE map clustered the cells into distinct populations, which included Tra1-81 + iPS cells, CK19 + hepatoblast cells, and CD26 + hepatocytes and cholangiocytes (Supplementary Fig. 1i). The protein expression signature also revealed that there was a continuous path of cellular differentiation within the HO cultures ( Supplementary Fig. 1j). We previously demonstrated that hepatic organoids could be dissociated, and that that individual cells could then reform organoids, which are referred to as secondary organoids (HO2) 1 . These HO2 have more epithelial cells and are more similar to an epithelial organoid 9 , the polarized epithelial cells reproduce some aspects of the native epithelium. To do this, the cells from a dissociated organoid are cultured in growth media (GM) for 6 days, and then in differentiation media for 6 more days to produce the HO2. After enzymatic dissociation of an epithelial organoid into single cells and secondary culture in growth media; the dissociated cells reorganize, proliferate and then reform organoids. In some supplemental figures, the scRNA-Seq and CyTOF data analyzes HO2 cells that were cultured in GM or DM.

Supplementary note 2: Second Harmonic Generation (SHG) microscopy. The spatial
distribution and morphology of collagen fibers within organoids were monitored by SHG microscopy 10 , 11 . SHG has been used to analyze human liver diseases that include fatty liver 12 and liver fibrosis 13,14 . With the use of a high-numerical aperture objective, early signs of fibrosis (i.e. when submicron-sized collagen fibers start to aggregate) can be detected. Of importance, the information provided by SHG microscopy cannot be deduced by fluorescence microscopy of immunostained collagen, since immunostaining also detects non-fibrous collagen. Picrosirus red staining (imaged using polarization microscopy) is not capable of resolving thin collagen fibers; it requires thin-sectioning of samples so that no information can be provided on the 3D distribution of fibers; and it has been shown to produce data that is significantly different from that provided by the SHG signal 15 . SHG analysis revealed that control HOs had a network of thin cross-linked collagen fibers (average diameter< 1.5 m) surrounding the cells in isolated regions. Coherent anti-Stokes Raman scattering (CARS) microscopy also showed that HOs had micron-sized intracellular lipid droplets (Fig. 1j). These results indicate that, besides synthesizing pro-collagen, cells within HOs have the enzymatic machinery required for crosslinking collagen to form fibers, as well as a functional lipid storing mechanism.

Supplementary note 3: developmental effect of ARPKD mutation.
To identify the developmental stage affected by the ARPKD mutation, we analyzed scRNA-Seq data generated from ARPKD and isogenic control iPSCs, hepatoblasts, and organoid cultures prepared from the 3 unrelated individuals. In total, the transcriptomes of 10,000 ARPKD and 10,000 isogenic control cells were analyzed. The developmental trajectories of ARPKD and isogenic control cells are quite similar at the iPSC (day 0) and hepatoblast (day 9) stages, but significantly differed at the organoid stage (Fig. 4c, Supplementary Figs.3a-e). We also analyzed previously obtained scRNA-Seq data 16 to determine when the primary cilium genes, which are mutated in ARPKD or ADPKD (PKD1), are expressed during organoid development. PKHD1 mRNA is expressed in iPSCs, but was significantly decreased at the hepatoblast stage, and then increased at the organoid stage. PKD1 mRNA is expressed at the hepatoblast and HO stages ( Supplementary   Fig. 3f). RT-PCR analyses indicated that equivalent levels of the mRNAs encoding four mesenchymal cell markers (COL1A1, PDGFRB, Vim, ACA2) were expressed in ARPKD and control hepatoblast cultures (Supplementary Fig. 3g). Taken together, these results indicate that the ARPKD mutation-induced effect on mesenchymal populations probably occurs when hepatoblasts differentiate into the cells that are present the mature liver organoid.

Supplementary note 4: similarities with commonly occurring forms of human liver fibrosis.
To investigate whether the ARPKD organoid fibrosis mechanistically resembled that in the commonly occurring forms of human liver fibrosis, 254 genes whose expression was increased in the cluster 0 cells in the ARPKD organoid were used to form a myofibroblast-specific expression signature (Supplementary Table 3). Gene Set Enrichment Analysis (GSEA) 17 was used to assess whether this myofibroblast expression signature was present in other types of fibrotic liver tissue. GSEA has been used to identify genes/pathways associated with treatment response or disease prognosis [18][19][20] , and to identify stem cell signatures in human cancer tissues 21,22 . GSEA calculates a normalized expression score (NES), which indicates whether myofibroblast signature genes are enriched in fibrotic liver tissue. GSEA analysis was performed using expression data obtained from 10 normal and 10 hepatitis C virus infectioninduced cirrhotic liver tissues (GSE6764 23 ). The myofibroblast expression signature was very strongly associated with cirrhotic liver (NSE 2.56, false discovery rate (FDR) 0), but not normal liver (NES -2.55, FDR 0) (Fig. 5k). We next investigated whether the myofibroblast signature was associated with non-alcoholic steatohepatitis (NASH), which is now the most common cause of chronic liver disease 24,25 . Although NASH is triggered by an abnormal triglyceride accumulation; fibrosis develops and progresses as NASH liver disease advances. Myofibroblast activation is key to its pathogenesis [26][27][28] , and the extent of liver fibrosis is the major determinant of NASH outcome 29,30 . Therefore, a gene expression dataset (GSE83452) containing 98 normal and 126 NASH liver tissues was analyzed. The myofibroblast expression signature was strongly associated with NASH liver (NSE 1.65, FDR 0), but not with normal liver tissue (NES -1.64, FDR 0). Of importance, in the absence of liver fibrosis, the myofibroblast expression in the cholangiocytes in ARPKD organoids is consistent with prior observations in ARPKD rodent models 31,32 . In cultured cells, the ARPKD mutation disrupts the interaction between FPC and the NEDD4 family member ubiquitin E3 ligase complex; and this enhances TGF- signaling by impairing TGF- receptor degradation 33 . Moreover, an increase in MMP-2 and MMP-9 mRNAs in ARPKD cholangiocytes could increase the conversion of latent TGF- into its active form 34,35 , which could further amplify the effect of the ARPKD mutation on TGF- signaling. STAT3 pathway activation in ARPKD myofibroblasts is consistent with evidence suggesting that this pathway plays a role in ADPKD, which is caused by a mutated membrane protein (polycystin-1, PC1) that forms a complex with FPC 36 . Increased expression of STAT3 pathway-associated mRNAs (Myc, PDGFR, PIM-1) in ARPKD myofibroblasts is also of interest. The marked increase in PDGFR protein expression is consistent with the well-known role of the PDGFR/STAT pathway in promoting hepatic fibrogenesis 37 , and PDGFR cross-linking leads to STAT3 phosphorylation and activation 38 . Myc, which is induced by PDGF in a STAT3-dependent manner 39 , promotes the proliferation of hepatic stellate cells and their conversion into myofibroblasts 40 . ARPKD myofibroblasts also had increased SOCS3 mRNA levels, which, under normal conditions downregulates the STAT3 signaling pathway 41,42,43,44,45 .
ARPKD myofibroblasts also have an increased level of LIF receptor mRNA expression, and they (along with cholangiocytes and possibly other cell types) produce LIF. The LIF receptor forms a cell membrane-localized complex with gp130 46 , whose expression is down-regulated by SOCS3 45 . Although ARPKD fibroblasts have increased SOCS3 mRNA levels, this "brake" on the system, which would normally reduce both gp130 (as an inducer of the STAT3stimulating pathway) and downstream elements thereof [41][42][43][44][45][46] could be overwhelmed by the combined activation of STAT3 via PDGFR and LIF receptor signals. Whereas TGF-1 alone was not able to induce cultured rodent ARPKD cholangiocytes to differentiate into mesenchymal cells 31 , it acts in concert with LIF to induce (in a STAT3-dependent manner) fibroblasts to develop into cells with activated and invasive properties 47 . This information along with our organoid data can be assembled into a potential model for the pathogenesis of ARPKD liver disease (Supplementary Fig. 9). In brief, the ARPKD mutation in PKHD1 generates cholangiocytes that produce an increased amount of TGF-1, as well as the mesenchymal cellderived enzymes and proteins involved in thick collagen fiber generation. The TGF-1 produced by ARPKD cholangiocytes acts in concert with LIF and downstream phospho-STAT3 to jointly stimulate mesenchymal cells to become activated myofibroblasts. There is also evidence that ARPKD cholangiocytes produce other factors that promote liver fibrosis 48 . As depicted in this model, TGF and STAT3 signaling pathways are known to interact 49 50 . Ligand binding by TGF receptors activates the signal transducing SMAD proteins 51