Abstract

   Deficiencies in the development of epithelial structures and delays in
   cellular maturation can increase the susceptibility of neonates to
   disease early in life. To investigate human biliary development and its
   vulnerability to biliary atresia, a severe pediatric cholangiopathy, we
   engineered multi-compartment biliary organoids (MBOs) from co-cultures
   of human liver-derived epithelial organoid cells with human endothelial
   and mesenchymal cells. MBOs derived from normal livers effectively
   replicated the epithelial structure of the bile duct epithelium and
   peribiliary glands (PBGs). Conversely, MBOs from diseased livers
   exhibited defective epithelial layers, a significant
   epithelial-mesenchymal transition (EMT), and an activation of the
   TGF-β/Activin-SMAD2/3 signaling, primarily due to intermediary cell
   sub-populations. Inhibition of TGF-β signaling suppressed EMT and
   promoted biliary epithelial development in human MBOs and suppressed
   the phenotype of experimental biliary atresia in neonatal mice. Thus,
   the modulation of TGF-β-dependent EMT regulates bile duct epithelial
   development and influences the susceptibility of neonates to biliary
   injuries.

   Subject terms: Mechanisms of disease, Transdifferentiation,
   Cholestasis, Liver fibrosis
     __________________________________________________________________

   The study reveals that TGF-β–driven epithelial-mesenchymal transition
   disrupts bile duct development in biliary atresia. Inhibiting this
   pathway restores epithelial structure and reduces disease features,
   highlighting a potential therapeutic target.

Introduction

   The normal development of epithelial compartments and functional
   cellular maturation are important for postnatal adaptation to
   extra-uterine life. There is increasing evidence that delayed
   maturation may render neonates susceptible to diseases, as evidenced by
   the time-restricted onset of biliary injury in neonates with biliary
   atresia (BA), a severe cholangiopathy of childhood. The biliary tract
   is composed of the gallbladder and bile ducts, with an epithelial
   lining that also contains peribiliary glands (PBGs) along the intra-
   (IHBD) and extrahepatic bile ducts (EHBD)^[42]1,[43]2. Although the
   epithelial lining allows for an uninterrupted flow of bile into the
   duodenum, it has other physiological functions, such as mucin
   production and stem cell reservoir^[44]3,[45]4. Thus, biliary
   epithelial integrity is crucial in preventing hepatic and biliary
   diseases^[46]5–[47]9. Although several factors have been linked to the
   etiology of the disease, the mechanisms of epithelial disruption and
   the rapid progression of fibrosis remain undefined.

   Bile duct epithelial injury by the administration of rhesus rotavirus
   (RRV) to neonatal mice recapitulates the tissue injury observed in
   human neonates^[48]10–[49]12. However, the recovery of viruses from
   human BA tissues is inconsistent. The best evidence that the
   developmental defects play a significant role in the prenatal onset of
   disease emerged from the detection of hyperbilirubinemia within days of
   birth in newborn infants who later develop BA^[50]13–[51]15. Although
   harmful mutations in several genes, such as PKD1L1, ADD3, EFEMP1,
   KIF3B, TTC17, and PCNT, have been linked to abnormal bile duct
   development, these mutations are detected in only a small subset of
   patients^[52]16–[53]18.

   Studies using patient-derived cholangiocyte organoids revealed a
   delayed maturation of epithelial cells, with altered tight junctions,
   functional defects, and mal-positioning of cilia^[54]19,[55]20. Based
   on these data, we hypothesized that the delayed maturation of the
   diseased biliary epithelium results from altered crosstalk between
   neighboring epithelial and mesenchymal cells. To test this hypothesis,
   we established human multi-compartment biliary organoids (MBOs) by
   co-culturing cells from normal- and patient-derived cholangiocyte
   organoids with human umbilical vein endothelial cells (HUVECs) and
   mesenchymal stem cells (MSCs). These cells self-organized into
   3-dimensional (3D) organoids, with HUVECs and MSCs localized centrally
   and underneath a layer of epithelial cells and PBGs. To generate a
   human multi-cellular BA model, we generated MBOs from the epithelial
   cells of BA patients (BA-MBOs). The BA-MBOs displayed markers of
   epithelial-mesenchymal transition (EMT) and upregulated transforming
   growth factor-β (TGF-β)/Activin-SMAD2/3 pathway without immune
   microenvironments. The importance of EMT and TGF-β in regulating the
   defects of the diseased epithelium was supported by a decrease in EMT
   and further maturation of epithelial cells upon blocking of TGF-β
   signaling in MBOs, and by a decreased biliary injury and fibrosis in
   neonatal mice treated with a TGF-β inhibitor.

Results

Cellular compartmentalization and formation of PBG-like structure in MBOs

   Normal bile duct development and regeneration requires an interaction
   of liver progenitor epithelial cells with neighboring mesenchymal and
   endothelial cells^[56]21–[57]23. When cholangiocyte organoids are
   derived from human livers, they form spheroids comprised of epithelial
   cells without identifiable PBGs^[58]19. Here, we investigated the
   regulatory role of non-epithelial cells in the development of the
   biliary epithelium and PBGs, and in the delayed epithelial maturation
   associated with BA. First, we dissociated normal cholangiocyte
   organoids and co-cultured them with HUVECs and human MSCs
   (Fig. [59]1a)^[60]24,[61]25. The mixture of three cell types was seeded
   on basement membrane extract matrix and incubated with an expansion
   media (EM)^[62]26 containing R-spondin 1, EGF, FGF10, TGF-β inhibitor
   A8301, HGF, Forskolin, Gastrin and an equal proportion of endothelial
   cell growth medium (EGM)^[63]27,[64]28 for 5 days. This was followed by
   culture in a cholangiocyte differentiation media (DM)^[65]29 containing
   R-spondin 1, EGF, FGF2, DKK1, Forskolin and Dexamethasone with equal
   proportion of EGM for 9 days. A separate group of cholangiocytes was
   subjected to the same culture conditions (type of media and number of
   days in culture) to serve as controls (Supplementary Table [66]1).
   Cholangiocytes cultured in the presence of HUVECs and MSCs
   self-organized into 3D MBOs containing a central cell cluster
   juxtaposed to an outer cell layer. The outer cell layer organized in
   MBOs grew by 4-fold between 5 and 10 days of culture (Fig. [67]1b). In
   contrast, cholangiocytes plated as a single cell type grew as
   mono-lineage organoids (termed “Mono” to reflect the epithelium-only
   compartment) and formed spheroids (Fig. [68]1c, Supplementary
   Movie [69]1).

Fig. 1. Cellular compartmentalization and formation of PBG-like structure in
MBOs.

   [70]Fig. 1
   [71]Open in a new tab

   a Overview of MBO generation. Created in BioRender. Ayabe, H. (2025)
   [72]https://BioRender.com/iy9dzsm. b Change in area of outer cell layer
   in MBOs from days 5-14 in culture (Mean ± SEM, n = 30 organoids, Multi:
   NC-MBO). c Phase contrast images of biliary organoids (Mono: NC-MBO
   without HUVEC and MSC, G: gland-like structure, E: epithelium-like
   structure). Three experiments were repeated with similar results. d
   Whole-mount immunostaining of CK19 (red), CD31 (green), VIM (white) and
   nucleus (blue; Hoechst 33342) in MBOs. Three experiments were repeated
   with similar results. e Immunostaining of CK19 (green), CD31 (white),
   VIM (red) and nucleus (blue; DAPI) in 3-day-old mouse extrahepatic bile
   duct (left) and MBO (right). G: Gland-like structure and E: Epithelium.
   Three experiments were repeated with similar results. f PCA of MBO
   transcriptional phenotype. RNA-sequencing data from MBOs were
   integrated into publicly available data^[73]32. The IHBD organoids were
   cultured in medium according to Huch et al., (1)^[74]26 or extrahepatic
   medium (2)^[75]32. g, h Immunostaining of EPCAM (white) and TROP2 (red)
   in MBOs (n = 3) in (g) and EPCAM (white), TROP2 (red) and nucleus
   (blue; DAPI) in 3 to 4-day-old mouse extrahepatic bile ducts (n = 3) in
   (h). G: Peribiliary gland and E: Luminal epithelium. i FACS of biliary
   organoids. Left: Representative plot for TROP2^+/EPCAM^+ and
   TROP2^low/EPCAM^+ cells. Middle: Frequency of TROP2^low cells in
   EPCAM^+ cells in MBOs (n = 6 organoids). Right: Number of
   TROP2^low/EPCAM^+ cells in each group (Mean ± SD, n = 6 organoids). j
   Immunostaining of SOX9 (red), SOX2 (red), CK19 (green) and nucleus
   (blue; DAPI) in MBOs (n = 3). b, i the p-values were shown in plots and
   were determined using one-way ANOVA with Tukey’s test.

   Immunostaining of MBOs revealed that the central cell cluster consisted
   of CD31 (also known as PECAM1) positive HUVECs and vimentin (VIM)
   positive MSCs. HUVECs maintained CD34 expression, which is a pattern
   previously described in portal vein endothelial cells (Supplementary
   Fig. [76]1a)^[77]30. For MSCs, the reduced expression of CD105 within
   MBOs supports a potential adaptation of MSCs to a hepatobiliary-like
   microenvironment of MBOs (Supplementary Fig. [78]1b)^[79]31. CK19 (also
   known as KRT19) positive cholangiocytes formed an outer cell layer and
   PBG-like structures embedded within the mesenchymal and endothelial
   cells, resembling the anatomic organization of bile duct tissue with
   PBGs (Fig. [80]1d, e)^[81]1. To gather insight into the transcriptional
   phenotype of MBO cells, we applied principal component analysis (PCA)
   to MBOs and published data curated from biliary tissues, and found that
   the MBOs are positioned away from cholangiocyte organoids and in
   proximity to common bile duct, gall bladder, and pancreatic duct
   tissues (Fig. [82]1f)^[83]32. Further immunostaining with antibodies
   against trophoblast surface protein 2 (TROP2 [also known as TACSTD2])
   stained the epithelial layer uniformly, while PBGs lacked TROP2 signal
   or stained faintly^[84]33,[85]34, with both epithelial compartments
   expressing epithelial cell adhesion molecule (EPCAM) (Fig. [86]1g).
   This pattern was similar to the TROP2 immunostaining of mouse
   extrahepatic bile ducts (EHBD), as additional evidence that the
   differential expression of TROP2 in the MBO epithelium supported the
   epithelial compartmentalization and the presence of PBG compartments
   (Fig. [87]1g, h). By fluorescence activated cell sorting (FACS), the
   number of PBG epithelial cells (TROP2^low/EPCAM^+) were ~2.7 and
   ~30-fold higher in a MBO compared to cholangiocytes prior to culture
   (“Chol”) and the mono-lineage (“Mono”) organoid group (Fig. [88]1i).
   Consistent with a property of housing stem cells^[89]1, PBG-like
   structures expressed the biliary tree stem/progenitor cell markers SOX9
   and SOX2 (Fig. [90]1j)^[91]35, and low level of LGR5 (Supplementary
   Fig. [92]1c)^[93]36. Functional biliary differentiation in NC-MBO was
   supported by the expression of CFTR, ASBT [also known as SLC10A2] and
   AE2 [also known as SLC4A2]^[94]4,[95]37 which was detected in a smaller
   frequency in epithelial cells of BA-MBOs (Supplementary Fig. [96]2a).

   To identify the role of each cell type in MBO formation, we examined
   the co-culture outcome of different combinations of the three cell
   types. The co-culture of cholangiocytes with MSCs (without HUVECs)
   self-assembled into a distinct central cluster of MSCs and an adjacent
   but separate epithelial cell population. This formation differed from
   the distribution of MSCs of MBOs, which joined HUVECs to form
   subepithelial cell populations with direct contact with the epithelium.
   Co-culturing cholangiocytes and HUVECs (without MSCs) failed to form
   organoids (Supplementary Fig. [97]2b, c). Altogether, these results
   suggest that the culture of cholangiocytes with HUVECs and MSCs forms a
   suitable microenvironment for the development of epithelial cells and
   PBG-like structures, thus reproducing the two epithelial cellular
   compartments of bile ducts.

Modeling human BA using patient-derived MBOs

   We next investigated the mechanisms underlying pathogenesis of BA by
   co-culturing HUVECs and MSCs with cholangiocyte organoid cells derived
   from liver biopsies of BA patients (BA-MBOs). BA-MBOs self-organized
   similarly to MBOs from normal control donors (NC-MBO) but had
   incomplete epithelial lining (Fig. [98]2a, [99]b and Supplementary
   Fig. [100]3a).

Fig. 2. Modeling human BA using patient-derived MBOs.

   [101]Fig. 2
   [102]Open in a new tab

   a Phase contrast and H&E staining of NC- and BA-MBOs. Three experiments
   were repeated with similar results. b Immunostaining of CK19 (green),
   CD31 (white), VIM (red) and Nucleus (blue; DAPI) in MBOs. Abundance of
   CK19^+ cells (n = 12 fields for NC or 14 fields for BA examined over 3
   independent samples) and average size of gland-like structures (n = 11
   fields for NC or 7 fields for BA examined over 3 independent samples)
   were quantified (Mean ± SD). c Immunostaining of ZO1 (green), β-catenin
   (red), CK19 (white) and nucleus (blue; DAPI) in MBOs. Graphs show
   CK19^+ cells exhibiting apical-basal polarity (Mean ± SD; n = 12 fields
   per 3 organoids). d Alcian Blue staining for MBOs. Mucin production in
   epithelium was determined by ratio of Alcian Blue^+ areas in CK19^+
   (from separately examined immunofluorescence data) areas (Mean ± SD;
   n = 15 fields [for NC] or 12 fields [for BA] examined over 3
   independent samples). e Heatmap (left) and Volcano plot (right) of
   differentially expressed genes between NC- and BA-MBOs (up in BA: 350;
   down in BA: 851). f DAB staining of MBOs for MMP-7 and OPN. The
   staining intensities were quantified by imaging analysis (Mean ± SD;
   n = 12 fields per 3 organoids). g Representative plots and
   quantification of FACS for hybrid cells (pan-CK^+/VIM^+) and
   cholangiocytes (pan-CK^+/VIM^−) in Mono and MBO. The number and cell
   ratio were quantified (Mean ± SD; n = 6). For b–d, f, g, the p-values
   were shown in plots and were determined using two-tailed Welch’s
   t-test.

   To investigate if epithelial cells of BA-MBOs displayed cell polarity
   defects^[103]19,[104]20, we first stained MBOs with anti-ZO1 (also
   known as TJP1) and anti-β-catenin antibodies (polarity markers).
   BA-MBOs displayed an aberrant staining pattern when compared to NC-MBOs
   (Fig. [105]2c). In addition, epithelial markers (CFTR, ASBT [also known
   as SLC10A2] and AE2 [also known as SLC4A2])^[106]4,[107]37 and
   acetylated tubulin (cilium protein) were predominantly stained in the
   epithelial layer of NC-MBOs, with decreased expression in BA-MBOs
   (Supplementary Figs. [108]2a and [109]3b). Using the ratio of Alcian
   Blue^+ to CK19^+ epithelium to evaluate the excretory ability of acidic
   mucin, NC-MBOs had strong staining in PBGs of MBO, with nearly absent
   signal in BA-MBOs (Fig. [110]2d). These findings provide evidence that
   the delayed maturation affects PBGs, in addition to the epithelial
   defects previously reported in BA. Despite these differences, the
   number of CK19^+ cells expressing Ki67 or the of CK19 colocalized with
   cleaved caspase-3, surrogate markers of proliferation and apoptosis
   respectively, was similarly low in NC- and BA-MBOs (Supplementary
   Fig. [111]3c, d).

   To examine the effects of long-term culture, we extended the culture of
   MBOs to 28 days. At this juncture, 3D rendering of images captured by
   confocal microscopy revealed that the structural integrity of NC-MBOs
   was preserved, with further growth of epithelial cells and PBGs
   juxtaposed to mesenchymal and endothelial cells. Although BA-MBOs also
   maintained anatomical organization with incomplete epithelial layer and
   rudimentary cells clusters resembling minute PBGs, these organoids had
   an expansion of VIM^+ cells and a ~ 3-fold reduction of epithelial cell
   volume (Supplementary Fig. [112]3e). These observations underscore the
   feasibility of extended MBO culture, highlight the persistent
   epithelial anomalies in BA-MBOs, and reveal an accentuated fibrotic
   phenotype, which typically characterizes disease progression in
   affected children.

EMT in patient-derived organoids

   To investigate the molecular basis of the abnormal epithelial
   development, we analyzed differentially expressed genes (DEGs) from
   RNA-sequencing and identified distinct transcriptional signatures
   showing increased expression of genes encoding matrix metalloproteinase
   7 (MMP-7 [also known as MMP7]) and osteopontin (OPN [also known as
   SPP1]), which play important roles in epithelial injury and the
   accentuated fibrosis of BA (Fig. [113]2e)^[114]38,[115]39. By
   immunohistochemical staining, MMP-7 and OPN expression was limited to
   the epithelial lining in NC-MBOs, while in BA-MBOs the staining was
   more intense in cells immediately adjacent to the thin epithelial
   lining (Fig. [116]2f). Based on these findings and on the role of
   endothelial and mesenchymal cells influencing bile duct development and
   pathogenesis of disease^[117]40–[118]43, we tested the hypothesis that
   epithelial cells of BA-MBOs display EMT. First, to examine if the
   development of MBOs involves EMT, we compared the baseline expression
   level of VIM by pan-cytokeratin (pan-CK) positive cells in NC-MBOs and
   NC-Monos (Supplementary Table [119]1) by FACS. MBOs had a significant
   increase in pan-CK^+/VIM^+ cells (Fig. [120]2g). Next, we investigated
   degree of EMT in patient-derived MBOs.

   To determine if EMT is preferentially increased in BA-MBOs, we first
   queried the RNAseq data from NC- and BA-MBOs using Gene Ontology (GO)
   and Molecular Signatures Database (MSigDB)^[121]44,[122]45. EMT
   signatures were enriched in the GO (TGFβ−3 [also known as TGFB3],
   SNAI1, GREM1) and MSigDB (INHBA, OPN, GREM1) results, in addition to
   TGF-β binding and receptor signaling, extracellular matrix (ECM), and
   collagen functions (Fig. [123]3a, [124]b and Supplementary
   Data [125]1)^[126]46, with the simultaneous increase of TGFβ−3, SNAI1
   and GREM1 in BA-MBO pointing to TGF-β pathway
   activation^[127]47,[128]48. In NC-MBOs, TGF-β3 expression was detected
   primarily in the immediate subepithelial compartment, which increased
   by 4-fold in BA-MBO (Fig. [129]3c). Additional immunostaining for
   DESMIN (DES), a MSC marker^[130]49, detected a similar pattern of
   DES^+/CK19^− cells in both MBO groups, while percentage of CK19
   overlapping with DES was increased in BA-MBOs (Supplementary
   Fig. [131]4a). Furthermore, immunostaining to evaluate for EMT showed
   increased expression of SNAI1, and ZEB2 in CK19^+ cells in BA-MBOs,
   with a decrease in epithelial co-staining for E-cadherin and CK19
   (ECAD^+/CK19^+ cells) (Supplementary Fig. [132]4b)^[133]47,[134]50.

Fig. 3. Presence and cellular lineage of EMT in patient-derived organoids.

   [135]Fig. 3
   [136]Open in a new tab

   a Top enriched terms of MSigDB Hallmark 2020^[137]45,[138]79 and GO
   sets^[139]44,[140]78 in genes upregulated in BA-MBOs vs NC-MBOs (MF:
   Molecular Function, BP: Biological Process, n = 3 organoids). b Boxplot
   visualization of EMT-related gene expressions in MBOs (n = 3
   organoids). The box shows median, 25th and 75th percentile values and
   the whiskers show 1.5 IQR from the edges of the box. c DAB staining of
   MBOs for TGF-β3. The staining intensities were quantified by imaging
   analysis (Mean ± SD; n = 12 fields per 3 organoids). d Representative
   plots and quantification of FACS of day 5 and 14 MBOs. Number of
   cholangiocytes (pan-CK^+/VIM^-), hybrid cells (pan-CK^+/VIM^+) and
   mesenchymal cells (pan-CK^−/VIM^+) in MBOs were quantified (Mean ± SD,
   n = 6 organoids). e UMAP plots of scRNAseq for MBOs colored by cell
   type (n = 4 organoids). f Top 5 enriched terms (plus ties) of MSigDB
   Hallmark 2020 sets in genes upregulated in Intermediate clusters (n = 4
   organoids). g Percent cell proportions of each cluster in MBOs (n = 4
   organoids). The box shows median, 25th and 75th percentile values and
   the whiskers show the furthest point that is within 1.5 IQR from the
   edges of the box. The outliers are plotted outside the whiskers. h EMT
   score overlaid in the UMAP plot (n = 4 organoids). The p-values were
   shown in plots and were determined using DESeq2 model (b) and
   two-tailed Welch’s t-test (c, d, g).

   To determine the EMT temporal dynamics during organoid formation, we
   quantified VIM^+ and pan-CK^+ cells using FACS. After 5 days of
   culture, no changes were observed between NC- and BA-MBOs. At 14 days,
   pan-CK^+/VIM^+ hybrid cells (transitional state or partial EMT^[141]43)
   were ~2-fold higher, pan-CK^−/VIM^+ mesenchymal cells were ~5-fold
   higher and pan-CK^+/VIM^– epithelial cells (cholangiocytes) were
   ~3-fold lower in BA-MBO compared to NC-MBO (Fig. [142]3d). Despite
   these changes in epithelial and mesenchymal cells, the population of
   endothelial cells was similar in NC- and BA-MBOs, as evidenced by
   whole-mount immunostaining using anti-CD31 and CK19 antibodies.
   Examining confocal images of organoids with the Imaris Filament Tracer
   Tool, both MBO groups contained cells forming short endothelial
   networks (Supplementary Fig. [143]4c).

Cellular lineage of EMT in patient-derived organoids

   To investigate the cell lineage of EMT, we performed single-cell RNA
   sequencing (scRNAseq) for MBOs. Transcriptome analysis of 27,645 cells
   identified 3 distinct primary cell clusters that predominantly
   expressed features of epithelial (KRT7, KRT8, KRT19, EPCAM),
   endothelial (PTX3, CDCA2, ANKRD37, PGF) and mesenchymal cells (FN1,
   VIM, TFPI2), 5 smaller clusters representing intermediate populations,
   and 1 other cluster showing indeterminate identity (S100A4, CD74,
   ARHGDIB, HLA-DRA, IFI27) (Fig. [144]3e, Supplementary Data [145]2). To
   initially visualize these intermediate cell populations, we used
   intensity-based co-localization analysis for percentage of CK19 signal
   overlapping with VIM signal. Using these markers, we localized these
   cells predominantly in BA-MBOs, where they populated the sites near the
   defective epithelium and in the subepithelial space (interspaced among
   mesenchymal and endothelial cells) (Supplementary Fig. [146]5a).

   Among these Intermediate clusters (Fig. [147]3f and Supplementary
   Fig. [148]5b), intermediate cell groups 1, 2 and 4 were enriched for
   EMT genes, and Intermediate 1 included partial EMT or transitory
   cholangiocyte genes (e.g., COL7A1, VEGFA, PDGFRB, COL6A3 and ZEB2),
   representing a subpopulation of epithelial cells with high mesenchymal
   potential. Intermediate 2 contained cells expressing mediator genes
   implicated in TGF-β-induced EMT (e.g., LGALS1, IGFBP7, TAGLN, SPARC and
   ZEB1), and cells in Intermediate 4 expressed ECM remodeling genes
   involved in TGF-β-induced EMT (e.g., LOX, FN1, TGFBI, MMP2, SERPINE1,
   COL1A2, COL3A1, SPARC, INHBA, TGFB1) (Supplementary Data [149]2 and
   [150]3). Cells in Intermediate clusters 2 and 4 were significantly
   over-represented in BA-MBOs relative to NC-MBOs (Fig. [151]3g). BA-MBOs
   also demonstrated an increased prevalence of cells scoring highly for
   top EMT-associated genes when compared to their normal control
   counterparts (Fig. [152]3h, and Supplementary Fig. [153]5c)^[154]51.
   Moreover, cell-cell interaction prediction identified that both NC and
   BA-MBOs had TGF-β ligand-receptor interactions (e.g., TGF-β1, INHBA and
   INHBB) in Intermediate 1 and 4, but only BA-MBOs had this same
   interaction in the Intermediate 2 population (Supplementary
   Fig. [155]6a, [156]b).

   Since MBOs showed an upregulation of genes related to TGF-β-induced EMT
   and ECM remodeling, we assessed the activation of fibrosis in different
   cell types and calculated ECM scores using an extracellular matrix gene
   list including ECM glycoproteins, collagens, and proteoglycans from
   GSEA^[157]52. This approach identified populations of mesenchymal and
   intermediate cells with higher ECM scores and high correlation with EMT
   scores (r = 0.934, p < 0.0001), especially in BA-MBOs (Fig. [158]4a–c).
   To further investigate the temporal dynamics of the mesenchymal cell
   population in MBOs, we quantified MBO cells using FACS at days 5 and
   14. Pan-CK^+/VIM^+ hybrid cells and pan-CK^−/VIM^+ mesenchymal cells
   increased in BA-MBOs (Fig. [159]4d, [160]e). Measuring the surface area
   of organoids at day 14, BA-MBOs had a modest decrease in size compared
   to NC-MBO (Fig. [161]4f), potentially indicating a contraction by the
   expanded population of Intermediate and Mesenchymal cells
   (Fig. [162]4f). These findings suggest that there is activation of
   EMT-related ECM substrates in the expanding mesenchyme and intermediate
   populations of BA-MBOs (Fig. [163]4c). To examine a potential role of
   endothelial cells in organoid structure and its potential influence on
   the abnormal phenotype of BA-MBOs, we performed ligand-receptor
   analysis using scRNAseq. The endothelial cell cluster had 25
   interactions involving VEGF ligands^[164]53–[165]55 in NC-MBOs and 27
   interactions in BA-MBOs. Visualization of these interactions inferred
   strong crosstalk implicated in VEGF signals in Mesenchymal, Epithelial,
   Endothelial, and Intermediate cell clusters 1, 3, and 4, with similar
   patterns observed in both NC and BA-MBOs. Conversely, no relevant VEGF
   interaction was detected in Intermediate clusters 2, 5, and Other
   (Supplementary Table [166]2).

Fig. 4. Analysis of MSC and EMT-expressing populations in MBOs.

   [167]Fig. 4
   [168]Open in a new tab

   a ECM score overlaid in the UMAP plot (n = 4 organoids). b Violin plots
   showing EMT (upper) and ECM (lower) scores in clusters of NC- and
   BA-MBOs (n = 4 organoids). Source Data file contains exact p-values for
   comparison between cell types. c Scatter plot showing correlation
   between EMT and ECM scores (n = 8 organoids). The circles were filled
   by test groups and colored by clusters. d MBO cell number kinetics by
   FACS at days 5 and 14, with graphs showing numbers of cholangiocytes
   (pan-CK^+/VIM^-), hybrid cells (pan-CK^+/VIM^+) and mesenchymal cells
   (pan-CK^−/VIM^+) in each MBO (mean ± SD, n = 6 organoids). e Ratios of
   VIM^+ cells to pan-CK^+ cells in MBOs on day 5 and 14. VIM^+ cells were
   divided by pan-CK^+ cells (mean ± SD, n = 6 organoids). f Surface area
   of MBOs at day 14 in culture (Mean ± SD, n = 18 organoids). The
   p-values were shown in plots and were determined using two-tailed
   Wilcoxon rank sum test with Benjamini–Hochberg correction (b), cor.test
   function in R, which utilizes an asymptotic t approximation (c),
   one-way ANOVA with Tukey’s test (d), and two-tailed Welch’s t-test (e,
   f). Two-tailed Spearman’s correlation test was used to assess the
   relationship between EMT and ECM scores (c).

   Altogether, these data identify phenotypic and transcriptional evidence
   of EMT in cell sub-populations. Additionally, the prominent signature
   of TGF-β/Activin-SMAD2/3 genes in intermediate clusters is consistent
   with its function in biliary development, remodeling, and
   injury^[169]56–[170]59, and may regulate crosstalk between neighboring
   cells in the diseased tissue.

EMT suppresses epithelial maturation via TGF-β signaling

   To directly determine if TGF-β signaling regulates EMT in BA, we
   cultured MBOs in A8301, a TGF-β inhibitor (Fig. [171]5a). In this
   experiment, the control groups were cultured for 5 days in EM/EGM
   (containing 2.5 μM A8301)^[172]60, followed by 9 days in DM/EGM (with
   DMSO and no A8301), while the experimental group of BA-MBO was under
   the same culture conditions containing a higher concentration of 10 μM
   A8301 throughout the 0- and 9-day periods. A8301 decreased the number
   of pan-CK^+/VIM^+ hybrid cells in BA-MBOs and increased pan-CK^+/VIM^-
   cholangiocytes without changing the number of pan-CK^−/VIM^+ cells
   (Fig. [173]5b). The use of a higher concentration of A8301 (100 μM
   A8301) in separate experiments showed that the treatment reduced
   pan-CK^+/VIM^+ hybrid cells and pan-CK^−/VIM^+ mesenchymal cells, but
   did not increase pan-CK^+/VIM^- cholangiocytes (Supplementary
   Fig. [174]6c–e). 3D rendering of whole-mount stained MBOs showed that
   A8301 restored the epithelial composition in BA-MBOs and formed a
   surface area of similar dimension of NC-MBOs. In addition, PBGs were
   easily identified and had similar size NC-MBOs (Fig. [175]5c). TGF-β
   inhibition also corrected epithelial defects of BA-MBOs, shown by
   a ~ 2-9-fold increase in the percentage of ECAD^+ cells compared to
   vehicle-treated BA-MBOs, and significantly improved the percent of
   ZO1^+ cells and mucin production (Fig. [176]5d–f). We subsequently
   investigated the impact of including A8301 in the EM/EGM expansion
   medium during the initial five days at a concentration of 2.5 μM on the
   culture of BA-MBOs. Omitting A8301 during this period did not
   significantly affect the volume of cholangiocytes and mesenchyme in
   BA-MBOs, nor the colocalization ratio of CK19 and VIM (Supplementary
   Fig. [177]7a, [178]7b). Altogether, these data show that the inhibition
   of TGF-β/Activin signaling reduced EMT in BA-MBOs and allowed for
   epithelial differentiation and expansion.

Fig. 5. TGF-β-induced EMT is a candidate regulator of BA phenotypes.

   [179]Fig. 5
   [180]Open in a new tab

   a Overview of in vitro A8301-mediated TGF-β signal inhibition in
   BA-MBOs. Created in BioRender. Ayabe, H. (2025)
   [181]https://BioRender.com/q1pk62z. b Representative plots and
   quantification of FACS of day 14 MBOs treated with A8301. Number of
   hybrid (pan-CK^+/VIM^+), cholangiocyte (pan-CK^+/VIM^-) and mesenchymal
   (pan-CK^−/VIM^+) cells were quantified (Mean ± SD; n = 12 organoids). c
   Images of MBOs treated with A8301 captured before collection and after
   whole-mount immunostaining (CK19: red, CD31: green, VIM: white).
   Abundance of CK19^+ structures (Mean ± SD; n = 9 organoids examined
   over 3 independent experiments) and size distribution of PBG-like
   structures were quantified (n = 39-104 clusters examined over 3
   independent experiments). d Immunostaining of ECAD (green), CK19 (red)
   and Nucleus (blue; DAPI) in MBOs treated with A8301. ECAD^+ cells in
   CK19^+ cells in the field were quantified (Mean ± SD; n = 12 fields
   from 3 organoids). e Immunostaining of ZO1 (green), CK19 (red) and
   Nucleus (blue; DAPI) in MBO treated with A8301. ZO1^+ cells in CK19^+
   cells per field were quantified (Mean ± SD; n = 12 fields from 3
   organoids). f Alcian blue staining of MBOs treated with A8301. Mucin
   production was determined by ratio of Alcian Blue^+ areas in CK19^+
   (from separately examined immunofluorescence data) areas (Mean ± SD;
   n = 12 fields from 3 organoids). g Top enriched terms of MSigDB
   Hallmark 2020^[182]45,[183]79 and GO^[184]44,[185]78 sets in serum
   proteins (reanalyzed from Lertudomphonwanit et al.^[186]38) upregulated
   in BA subjects (n = 175 patients) vs age-matched normal controls
   (n = 9). MF: Molecular Function, BP: Biological Process. h Heatmap of
   enriched proteins in the EMT pathway in MSigDB Hallmark 2020 sets.
   (n = 9 NC and 175 BA subjects). i Serum TGF-β/Activin ligand levels
   represented by the relative fluorescent unit (RFU; Mean ± SD). (n = 9
   NC and 175 BA subjects). The p-values were shown in plots and were
   determined using one-way ANOVA with Tukey’s test (b–f) and two-tailed
   Welch’s t-test (i).

   Next, we analyzed serum samples from BA patients at the time of
   diagnosis using large-scale proteomics data reported by our group
   previously^[187]38; serum from age-matched healthy subjects served as
   control samples (normal controls or NC). We found a prominent EMT
   signature in the serum of BA patients using MSigDB hallmark 2020 gene
   set (including INHBA, TGFBI, TIMP1), in addition to inflammation
   (CXCL6, CXCL8 and CCL5) and ECM (COL23A1, LAMA1 and BGN) using GO terms
   (Fig. [188]5g and Supplementary Data [189]4). The heatmap generated by
   the EMT hit genes distinctly divided the serum samples into two
   clusters of NC and BA (Fig. [190]5h). Analyzing the expression of each
   TGF-β/Activin ligand on EMT in BA patients, we found that TGF-β1,
   TGF-β2, TGF-β3 and INHBA expression levels were significantly increased
   above normal controls (Fig. [191]5i). Next, we directly investigated
   the relevance of this pathway by inhibiting TGF-β signaling in a mouse
   model of experimental biliary atresia.

Inhibition of TGF-β signaling decreases biliary injury and liver fibrosis in
experimental BA

   To determine if TGF-β-dependent EMT regulates epithelial injury in
   experimental BA, we infected BALB/c mice with Rhesus rotavirus (RRV)
   within 24 hours of birth; phosphate buffered saline was injected in
   control mice (Fig. [192]6a)^[193]10–[194]12. Analysis of extrahepatic
   bile ducts (EHBDs) by immunostaining showed that the percent of
   TGF-β1^+/CK19^+ cells and TGF-β3^+/CK19^+ cells increased ~3.8 and
   ~1.8-fold, respectively, above controls 7 days after RRV infection
   (time of biliary epithelial injury and duct obstruction) (Supplementary
   Fig. [195]7c, [196]d). Then, we repeated the experiments and
   administered A8301 (1 mg/kg body weight) intraperitoneally daily
   beginning 1 day after infection with RRV; phosphate-buffered saline was
   injected in a separate group of control mice. A8301 injections reduced
   the percentage of RRV-infected mice with jaundice, promoted long-term
   survival, and decreased serum bilirubin and alanine aminotransferase
   (Fig. [197]6b–e). Histologically, the biliary epithelium of EHBDs was
   largely intact, with a mild subepithelial inflammation. In liver
   sections, portal spaces had mild expansion, which contrasted to the
   full expression of the BA phenotype in untreated mice (Fig. [198]6f).
   Tissue staining for EMT showed that the coexistences of ECAD or CK19
   and VIM induced by RRV infection was suppressed by A8301
   (Fig. [199]6g). Applying a similar strategy to a neonatal mouse model
   of chronic hepatitis induced by a delayed infection of RRV^[200]61, we
   found that A8301 suppressed the development of cholestasis, liver and
   biliary injury, and hepatic fibrosis (Fig. [201]6h–m). Altogether, the
   rescue of the BA phenotypes by TGF-β inhibition point to a key role of
   the TGF-β pathway in pathogenesis of experimental BA.

Fig. 6. Clinical relevance and therapeutic potential of TGF-β-induced EMT in
BA.

   [202]Fig. 6
   [203]Open in a new tab

   a RRV-induced BA mouse model with biliary injury (biliary-injury
   model). Newborn were injected intraperitoneally (IP) with PBS or RRV
   followed by daily PBS or A8301. Created in BioRender. Ayabe, H. (2025)
   [204]https://BioRender.com/b5m1rv9. b–e Jaundice and survival rates
   (n = 9 control [PBS PBS]/23 BA [RRV PBS]/12 BA-treated [RRV A8301]
   mice), and total bilirubin and plasma ALT (mean ± SD; n = 6 control/6
   BA/5 BA-treated mice) in biliary injury model. f H&E-stained EHBD and
   liver from day 12-14 biliary-injury model. g Immunostaining of VIM
   (green), ECAD or CK19 (red) and Nucleus (blue; DAPI) in day-4 and day-7
   biliary-injury model. VIM and ECAD or CK19 colocalization in
   control/BA/BA-treated mice were quantified (ECAD: n = 32/32/41 fields
   examined over 3 mice [day-4], n = 48/69/60 fields examined over 4 mice
   [day-7], CK19: n = 24/23/31 fields examined over 3 mice [day-4],
   n = 49/66/60 fields examined over 4 mice [day-7]) (mean ± SD). h
   RRV-induced BA mouse model with liver fibrosis (fibrosis model). 3-day
   old mice were injected IP with PBS or RRV followed by daily PBS or
   A8301. Created in BioRender. Ayabe, H. (2025)
   [205]https://BioRender.com/6fpm0xa. i–l Jaundice and survival rates
   (n = 24 BA/13 BA-treated mice), total bilirubin and plasma ALT
   (mean ± SD; n = 6 control/6 BA/5 BA-treated mice) in fibrosis model.
   i ^‡p < 0.01. Source Data file contains exact p-values. m
   Representative images of H&E-stained liver and EHBD (n = 6 control/8
   BA/6 BA-treated mice) and Sirius Red-stained liver (n = 6 mice in each
   group). Experiments were repeated with similar results. n Tissue
   expression of TGF-β in human BA. Immunostaining of CK19 (C; red),
   TGF-β1 (T; white), VIM (V; green) and DAPI (N; blue; nucleus) in a
   liver biopsy from BA patient. Graphs depict colocalization of CK19 and
   TGF-β1 or VIM (mean ± SD; n = 24 fields from 3 subjects). The p-values
   were shown in plots and were determined using two-tailed Fisher’s exact
   test (i) with Bonferroni correction (b), Log-rank test (j) with
   Bonferroni correction (c), one-way ANOVA with Tukey’s test (d, e, g, k,
   l) and two-tailed Welch’s t-test (n).

   To further examine the expression of TGF-β signaling in human liver
   biopsies, we immunostaining histological sections for TGF-β1 and TGF-β3
   from BA patients. The percent of TGF-β1^+/CK19^+ cells and
   TGF-β3^+/CK19^+ increased ~5.3 and ~11-fold in BA patients
   (Fig. [206]6n and Supplementary Fig. [207]8a). Analysis of liver gene
   expression data ([208]GSE46995; reported by our group
   previously^[209]62) showed increased expression of SMAD2/3 pathway
   ligands of TGFB1, TGFB2, TGFB3, and the receptor of
   TGFBR1^[210]63–[211]67 in BA patients when compared to liver biopsies
   from normal subjects (NC) or from patients with intrahepatic
   cholestasis (IHC; diseased control) (Supplementary Fig. [212]8b,
   [213]c, and Supplementary Data [214]5). Using gene set enrichment
   analysis, EMT emerged as the most enriched term from the MSigDB
   hallmark set in genes upregulated in BA patients compared to NC and IHC
   (Supplementary Fig. [215]8d). Searching for evidence of EMT in BA
   livers, we immunostained histological sections for CK19 and VIM. The
   percent of VIM^+/CK19^+ cells was ~1.8-fold higher in human BA
   (Fig. [216]6n). Altogether, these data indicate that the EMT signature
   and TGF-β expression identified in diseased biliary epithelium of MBOs
   are also detected in tissues of BA patients.

Discussion

   Our study reveals that the co-culture of biliary epithelial cells with
   endothelial and mesenchymal cells leads to the spontaneous
   self-organization into an organoid exhibiting an epithelial layer with
   peribiliary glands embedded within a subepithelial compartment with
   endothelial and mesenchymal cells. This arrangement replicates the dual
   compartment of the biliary epithelium, thus forming a more
   comprehensive multi-lineage system to enable studies that take into
   account molecular processes and signals that more closely approaches
   the complexity of the human biliary system. In patient-derived BA-MBOs,
   the epithelium was fragmented and PBGs were small, had abnormal
   intercellular junction (ZO1 expression), and expressed markers of EMT,
   with a prominent expression of TGF-β3. Upon TGF-β inhibition, EMT
   decreased in the diseased epithelium, the population of PBGs increased,
   and ZO1 expression improved. In vivo, TGF-β inhibition suppressed the
   BA phenotype and substantially decreased liver fibrosis in neonatal
   mice. These data point to a regulatory role of mesenchymal and
   endothelial cells in the development of PBGs and in the modulation of
   TGF-β-driven EMT of the diseased biliary epithelium in BA.

   The emergence of PBGs in MBOs underscores the critical role of
   mesenchymal and endothelial cells in regulating the development of
   epithelium and specialized epithelial glandular structures such as
   PBGs, which contain a progenitor cell niche to maintain epithelial and
   gland population^[217]1,[218]3,[219]36. While the mechanisms by which
   non-epithelial cells promote PBG development remain elusive, the
   activation of Wnt signaling has been identified as crucial for
   regulating the positional identity of bile duct
   epithelium^[220]36,[221]68. In our experimental system, data from
   BA-MBOs show that TGF-β signaling is functionally linked to abnormal
   PBG development and to EMT of the epithelium proper. This is consistent
   with detailed phenotyping of diseased tissues from affected patients,
   which showed abnormal remodeling of hilar ducts, active fibroplasia,
   and attempts to renew the epithelium^[222]69–[223]71. The emergence of
   PBGs and epithelial expansion after EMT suppression via inhibition of
   TGF-β signaling highlights the significance of coordinated cellular
   diversity in epithelial development and the potential reversibility of
   EMT in disease states.

   EMT plays a role in physiologic tissue remodeling and in regulating the
   fibrosis of multiple organs^[224]43,[225]72. Staining of EMT markers
   and FACS in NC-MBOs showed evidence of EMT in epithelial cells with
   anatomical and population features of cholangiocytes, suggesting
   continued normal remodeling at day 14. A potential cellular source
   emerged from scRNAseq of NC- and BA-MBOs, which identified the
   Intermediate 1, 2 and 4 clusters containing transcriptional enrichment
   of ECM-related genes on a variety of levels. ECM and EMT create
   feedback loops which can be regulated by TGF-β
   signaling^[226]46,[227]73–[228]75. Future analysis of the
   transcriptional program uniquely upregulated in these intermediate
   clusters, such as VIM, SNAI1, SNAI2, ZEB2 and FN1, may give insight
   into the molecular circuitry driving EMT and the developmental delay in
   epithelial and glandular in diseased states.

   The engineering of multilineage human biliary organoids to study
   epithelial development elucidates an interplay between epithelial cells
   and their non-epithelial neighbors, underscoring the pivotal role of
   endothelial and mesenchymal cells in the formation and maintenance of
   epithelial and glandular structures within the bile duct. This MBO
   formation occurs despite the fact that MSCs and HUVECs are not derived
   from livers. Given the specificity of the niche cells in the liver,
   with different mesenchymal and endothelial subtypes specializing in
   their functions in vivo, MBOs therefore may not fully recapitulate the
   full phenotype of cell populations of the liver. Using patient-derived
   MBOs and virus infection of neonatal mice to model human BA, we
   demonstrate how disruptions in this interplay, particularly through
   TGF-β-dependent EMT, contribute to delayed development and lead to
   epithelium pathology and tissue fibrosis. In rotavirus-infected mice,
   residual tissue inflammation likely stems from the persistent
   activation of TGF-β-independent circuits triggered by rotavirus
   infection, which are mediated by tissue macrophages, dendritic cells,
   lymphocytes, and NK cells^[229]7. This notwithstanding, the
   effectiveness of TGF-β inhibition in suppressing EMT markers and
   improving the epithelial and glandular architecture points to a
   critical role of TGF-β in the pathogenesis of BA and as a potential
   therapeutic target. Our study also identifies areas of future research
   to fully elucidate the molecular mechanisms governing epithelial
   crosstalk with neighboring mesenchymal and endothelial cells during
   tissue development and disease states.

Methods

Human biopsy samples

   Patients were recruited sequentially as they visited at a center in the
   consortium, which comprises of 15 centers. Human samples were obtained
   after informed consent from patients’ parents/legal guardians. All
   procedures involving the use of human specimens were approved by the
   Research Ethics Committees at Cincinnati Children’s Hospital Medical
   Center (CCHMC). Liver biopsy fragments used to generate cholangiocyte
   organoids were obtained from male and female subjects younger than 18
   years of age, as follows: 9 from BA patients (n = 3 at diagnosis, ages:
   5 weeks to 2 months; and n = 6 at the time of transplantation, ages: 6
   months to 2 years of age), and 6 from normal deceased donors (one of
   which had urea cycle defect; ages: 6 months to 17 years of age). Liver
   tissues used in histological analysis were from male and female
   pediatric subjects with BA (n = 3, ages: 8 months to 7 years of age)
   and normal deceased donors (n = 3, ages: 13-17 years of age).
   Demographic information was unavailable for one normal liver transplant
   donor (Supplementary Table [230]3).

Animals

   BALB/c mice used for the studies were obtained from Charles River
   Laboratories. Mice were housed in the animal vivarium with a 12-hour
   dark-light cycle and ad libitum access to normal diet and water at
   CCHMC in accordance with the guidelines recommended by the
   Institutional Animal Care and Use Committee (IACUC; Protocol
   IACUC2020-0006).

Mesenchymal cells and endothelial cells

   Human mesenchymal stem cells (MSCs) and human umbilical vein
   endothelial cells (HUVECs) were kindly gifted by Dr. Akihiro Asai
   (CCHMC) and cultured in mesenchymal stem cell growth medium (MSCGM;
   Lonza, PT-3001) and endothelial cell growth medium (EGM; Lonza,
   CC-3124).

Cholangiocyte organoids

   Cholangiocyte organoids were generated from human liver
   tissues^[231]19,[232]26. Minced human liver tissues were washed in
   ice-cold basal medium containing Advanced Dulbecco’s Modified Eagle
   Medium (DMEM)/F12 (Thermo Fisher Scientific, 12634028) supplemented
   with 100 U/mL Penicillin/Streptomycin (Thermo Fisher Scientific,
   15140122), 1% GlutaMAX (Thermo Fisher Scientific, 35050061), and 10 mM
   HEPES (Thermo Fisher Scientific, 15630080). Tissue fragments were
   suspended in Cultrex Reduced Growth Factor Basement Membrane Extract
   Type 2 Pathclear (BME; Bio-Techne, 3533-005-02) and plated into 24 well
   plates as a 50 μL droplet/well. After incubating at 37 °C to solidify,
   the tissue domes were initially cultured using 500 μL/well of Isolation
   Medium (IM) containing basal medium supplemented with 1× B27 Supplement
   minus vitamin A (Thermo Fisher Scientific, 12587010), 1× N-2 Supplement
   (Thermo Fisher Scientific, 17502048), 1 mM N-acetylcysteine
   (Sigma-Aldrich, A7250), 10% vol/vol homemade R-spondin 1 conditioned
   media^[233]19, 10 mM Nicotinamide (Sigma-Aldrich, N0636), 10 nM
   Recombinant human Gastrin (Sigma-Aldrich, G9145), 50 ng/ml recombinant
   human Epidermal Growth Factor (EGF; PeproTech, AF100-15), 100 ng/ml
   recombinant human Fibroblast Growth Factor 10 (FGF10; PeproTech,
   100-26), 25 ng/ml recombinant human Hepatocyte Growth Factor (HGF;
   PeproTech, 100-39H), 10 μM Forskolin (Bio-Techne, 1099), 5 μM A 83-01
   (Sigma-Aldrich, SML0788), 25 ng/ml recombinant human Noggin
   (Bio-Techne, 6057-NG), 30% vol/vol homemade Wnt Family Member 3 A
   (WNT3A) conditioned medium^[234]19, and 10 μM Y-27632 (Bio-Techne,
   1254) for 5 days. IM was then replaced with Expansion Medium (EM)
   consisting of IM without Noggin, WNT3A and Y-27632. After reaching high
   density, the generated organoids were mechanically dissociated into
   small fragments and transferred to new plates.

Generation of human multi-compartment biliary organoids (MBOs)

   To establish the MBOs, cholangiocyte organoids, MSCs and HUVECs were
   separately expanded. Then, cholangiocyte organoids were dissociated
   using Dissociation Buffer (Accutase [Innovative Cell Technologies,
   AT-104] and 0.1 mg/mL Deoxyribonuclease I [Sigma-Aldrich, DN25] and
   10 µM Y-27632) for 10 minutes at 37 °C. Cell suspensions were washed
   and passed through 40 μm cell strainer (Fisher Scientific, 22-363-547).
   Suspension medium was replaced with EM containing 10 µM Y-27632. MSCs
   and HUVECs were trypsinized with 0.05% Trypsin-EDTA (Thermo Fisher
   Scientific, 25300054). Viable 4.7 × 10^5 cholangiocyte organoid cells,
   3.3 × 10^5 HUVECs, and 1.0 × 10^5 MSCs were mixed, and cell pellets
   were gently plated on a 48-well plate coated with 200 µL of BME/EGM
   mixture (BME:EGM = 1:1). Cells were co-cultured for 5 days in EM/EGM
   (EM:EGM = 1:1), followed by 9 days in Differentiation Medium (DM)/EGM
   (DM:EGM = 1:1) to establish the MBO. The DM contained 10 mM
   Nicotinamide, 17 mM Sodium bicarbonate (Sigma-Aldrich, S5761), 0.2 mM
   2-phospho-L-ascorbic acid trisodium salt (Sigma-Aldrich, 49752),
   0.63 mM Sodium pyruvate (Thermo Fisher Scientific, 11360070), 14 mM
   Glucose (Sigma-Aldrich, G8644), 20 mM HEPES, 1% GlutaMax, 100 U/mL
   Penicillin/Streptomycin, 1% ITS+ Premix (Corning, 354352), 10 μM
   Forskolin, 10 nM Dexamethasone (Bio-Techne, 1126), 10% vol/vol in-house
   R-spondin1 conditioned media, 20 ng/mL recombinant human EGF, 5 ng/mL
   recombinant human FGF basic (Bio-Techne, 233-FB), and recombinant human
   Dkk-1 protein (Bio-Techne, 5439-DK) in William’s E medium (Thermo
   Fisher Scientific, 12551032). The cholangiocyte organoid cells in MBOs
   were collected from multiple human subjects, and the same MSCs and
   HUVECs were used throughout this study.

   Cholangiocyte organoids were dissociated and phenotyped by FACS prior
   to being subjected to further culture (forming a group termed “Chol”).
   Dissociated 4.7 ×10^5 cholangiocytes were subjected to the same culture
   conditions (type of media and number of days in culture) to serve as a
   control group termed “Mono-lineage biliary organoid” (or “Mono”)
   (Supplementary Table [235]1). Unless otherwise noted, all MBO
   experiments were performed on day 14 of culture. Images were captured
   using either a compound (Olympus, IX71) or stereo (Leica, M205A)
   microscope. CellSens (RRID: SCR_014551, Olympus) or Leica Application
   Suite (RRID: SCR_016555, Leica) were used for visualization. Image
   J/Fiji (RRID: SCR_002285)^[236]76 was used for image quantification.
   Live imaging was performed using Celldiscover 7 (ZEISS) and visualized
   using ZEN software (RRID:SCR_013672, ZEISS).

Whole-mount immunostaining

   MBOs were fixed in 4% paraformaldehyde (PFA; Electron Microscopy
   Sciences, 15710)/DPBS, followed by permeabilization with 4% PFA/0.5%
   Triton X-100 (Sigma-Aldrich, T9284)/DPBS and 0.1% Tween 20 (Fisher
   scientific, BP337)/DPBS. MBOs were, then, cleared using PEGASOS
   protocol^[237]77. In brief, the MBOs were decolorized overnight with
   25% Quadrol (Sigma-Aldrich, 122262)/5% Ammonium hydroxide (Fisher
   scientific, 423305000) and washed with DPBS. After washing, the MBOs
   were treated with blocking solution (5% normal donkey serum [Jackson
   ImmunoResearch, 017-000-121]/5% normal goat serum [Jackson
   ImmunoResearch, 005-000-121]/0.3% Triton X-100) at 4 °C overnight and
   incubated at 4 °C for 2 days with primary antibodies diluted in 1%
   Bovine Serum Albumin (Jackson ImmunoResearch, 001-000-162)/0.3% Triton
   X-100/DPBS. Samples were then incubated at 4 °C overnight with
   fluorescent dye-conjugated secondary antibodies followed by nuclear
   staining with Hoechst 33342 (Thermo Fisher Scientific, H1399). Specific
   antibody details are listed in the Supplementary Table [238]4 and
   Reporting summary. MBOs were sequentially delipidated using 30%, 50%
   and 70% tert-Butanol (tB; Sigma-Aldrich, 360538) solution containing 3%
   Quadrol and dehydrated with tB-PEG medium (75% tB/22% poly [ethylene
   glycol] methacrylate [PEGMMA; Sigma-Aldrich, 409529] and 3% Quadrol)
   for 1 day. The MBOs were finally treated with benzyl benzoate (BB;
   Sigma-Aldrich, B6630)-PEG clearing medium (75% BB/22% PEGMMA/3%
   Quadrol) until the tissue turned transparent. For cilia staining, MBOs
   were fixed in 4% PFA/DPBS and washed with DPBS, followed by
   permeabilization and blocking with 1% BSA/0.5% Triton X-100/DPBS (P/B
   buffer) at room temperature for 1 hour. MBOs were then incubated at
   4 °C for 1 day with primary antibodies in P/B buffer and washed with
   0.2% Triton X-100/DPBS. After washing, MBOs were incubated with
   fluorescent dye-conjugated secondary antibodies/Hoechst 33342 in P/B
   buffer at 4 °C overnight. Whole-mount immunofluorescence images were
   captured using confocal microscope (Nikon, ECLIPSE Ti; Zeiss, LSM880)
   and NIS Elements (RRID: SCR_014329, Nikon), Black Zen software
   (RRID:SCR_018163, Zeiss) and Imaris (RRID: SCR_007370, Oxford
   Instruments) were used for visualization and quantification of images.

FACS analysis

   To obtain single cells, MBOs were washed several times with DPBS and
   dissociated as mentioned above. Dissociated cells were washed and
   passed through 40 µm cell strainer and stained in DPBS using the Zombie
   Violet Fixable Viability Kit (BioLegend, 423113). Flow cytometry (FACS)
   buffer (DPBS/2% Fetal Bovine Serum [FBS; Thermo Fisher Scientific,
   26140079]) was used for subsequent wash steps. For pan-CK^+/VIM^+
   staining, cells were incubated with BUV395-conjugated CD31 antibody (BD
   Bioscience, 565290) and PerCP-Cy5.5-conjugated CD73 antibody (BD
   Bioscience, 561260), washed with FACS buffer and fixed in IC Fixation
   Buffer (Thermo Fisher Scientific, 00-8222). Fixed cells were
   subsequently washed with 1× permeabilization buffer (Thermo Fisher
   Scientific, 00-8333), followed by intracellular-staining with Alexa
   fluor 488-conjugated pan-CK (Thermo Fisher Scientific, 53-9003-82) and
   VIM antibodies (Abcam, ab92547) in 1× permeabilization buffer.
   Following washing, VIM-specific signals were obtained by incubation
   with Alexa fluor 647-conjugated anti-rabbit IgG antibody (Cell
   Signaling Technology, 4414). For TROP2^+/EPCAM^+ staining, cells were
   incubated with PE-conjugated TROP2 (BioLegend, 363804) and EPCAM
   antibodies (R and D Systems, AF960) in FACS buffer after the live/dead
   staining described above. After washing, EPCAM-specific signals were
   obtained using Alexa fluor 647-conjugated anti-goat IgG antibody
   (Jackson ImmunoResearch Labs, 705-605-147). Washed cells were subjected
   to FACS analysis using LSRFortessa Cell Analyzer (BD Bioscience).
   Appropriate single stains and fluorescence minus one (FMO) controls
   were also performed. Cells were sequentially gated based on size and
   granularity (forward scatter [FSC] versus side scatter [SSC]) and
   singlets (FSC-Area versus FSC-Height). Dead cells were excluded using
   Zombie Violet. Cholangiocytes (pan-CK^+/VIM^–), mesenchymal cells
   (pan-CK^–/VIM^+) and hybrid cells (pan-CK^+/VIM^+) were selected from
   CD31^– population. TROP2^+ epithelial cells (TROP2^+/EPCAM^+) and
   TROP2^– epithelial cells (TROP2^–/EPCAM^+) were selected from EPCAM^+
   population (Supplementary Fig. [239]9, Supplementary Table [240]4 and
   Reporting summary). Data were analyzed using FlowJo (RRID: SCR_008520,
   BD Bioscience).

Experimental mouse BA models

   For testing in vivo effects of TGF-β signal inhibition, we used two
   murine models. For modeling biliary injury and obstruction in BA,
   newborn BALB/c mice were intraperitoneally injected with 1.5×10^6
   fluorescent-forming units (ffu) of RRV within 24 hours after birth
   (Fig. [241]6a)^[242]10,[243]12. To model liver fibrosis in BA, 3-day
   old BALB/c mice were intraperitoneally injected with 1.875 × 10^6 ffu
   of RRV (Fig. [244]6h)^[245]61. Mice injected with DPBS were used as
   normal controls. To inhibit TGF-β signaling, 1 mg/kg of A8301 sodium
   (MedChemExpress, HY-10432A) in 20 µL of DPBS was injected
   intraperitoneally and daily beginning 6 hours after RRV inoculation.
   Mice were monitored daily for jaundice, body weight, and survival.
   EHBDs and livers were harvested, and blood was collected with
   Microvette (Sarstedt, 16.444.100) to obtain plasma following
   IACUC-approved protocols. Plasma bilirubin and ALT were analyzed using
   Total Bilirubin Reagent Set (Pointe Scientific, B7538-120) and
   DiscretPak ALT Reagent Kit (Catachem, V164-0A) and levels were
   quantified using a Synergy H1 Hybrid Reader (BioTek).

Imaging histological staining and fluorescence

   MBOs were collected in DPBS by manual pipetting, fixed in 4% PFA/DPBS
   at 4 °C overnight and placed in 70% ethanol at 4 °C after DPBS
   washings. Mouse EHBDs and livers were fixed in 10% formalin (Fisher
   scientific, SF100-4) at room temperature. Routine paraffin embedding
   and serial 5 µm sectioning were performed at the Research Pathology
   core at CCHMC. Sections were baked at 60 °C and deparaffinized using
   Histo-Clear (National Diagnostics, HS-200), followed by re-hydrating
   with graded ethanol-water solutions. Heat-induced epitope-retrieval was
   performed using 1 × Target Retrieval Solution pH 6.1 (Agilent, S1699)
   or 1 × Target Retrieval Solution pH 9.0 (Agilent, S2367). For
   immunohistochemistry, sections were incubated in 3.0% hydrogen peroxide
   (Sigma-Aldrich, 216763) and stained with primary antibodies and an
   avidin-biotin complex kit (Vector Laboratories, PK-6101 or PK-4002)
   with peroxidase substrate (Vector Laboratories, SK-4100). For
   immunofluorescence, cells were washed in DPBS then fixed in 4% PFA and
   washed with DPBS. MSCs were permeabilized by 0.1% Triton X-100/DPBS and
   HUVECs by 0.2% Tween 20/DPBS. The permeabilized cells or
   antigen-retrieved sections were blocked using serum from the same host
   species as the secondary antibodies and incubated with primary
   antibodies diluted in Antibody diluent (Agilent, S0809), followed by
   secondary antibodies. Antifade mounting medium with DAPI (Thermo Fisher
   Scientific, [246]P36931) was used for covering the slides. To detect
   human TROP2 and human EPCAM, mouse anti-human TROP2 antibody (DSHB,
   CPTC-TACSTD2-2) and goat anti-human EPCAM antibody (R and D Systems,
   AF960) were used. For mouse, goat anti-mouse TROP2 antibody (R and D
   systems, AF1122) and rat anti-mouse EPCAM antibody (BD Bioscience,
   552370) were used. Specific antibody details are listed in the
   Supplementary Table [247]4 and Reporting Summary. Conventional
   Hematoxylin (Poly Scientific R&D, s212) & Eosin (Sigma-Aldrich,
   HT1101128) staining to characterize morphology, Alcain blue
   (Richard-Allan Scientific, 88043) to identify acidic mucin production,
   and Sirius red (Dudley Corporation, 74189) staining for evaluation of
   liver fibrosis were performed. Bright field images were captured using
   a light microscope (Olympus, BX43) and immunofluorescence images were
   captured using either fluorescence (Olympus, BX51; Nikon, ECLIPSE Ni-E)
   or confocal (Nikon, ECLIPSE Ti; Zeiss, LSM880) microscopes. Imaris and
   Image J/Fiji were used for visualization and quantification of images.

Image analysis

   Whole-mount immunostaining images of MBOs were analyzed using Imaris
   (Oxford Instruments) software. Surface rendering objects of the red
   channel (CK19) were created by the Surface Creation Wizard to obtain
   epithelial cell volumes and PBG sizes. Total number and average surface
   area, length, and branching number of CD31^+ green channel surfaces in
   MBOs were quantified using Imaris Filament Tracer Tool. Staining
   intensity level in the positive area of DAB staining was measured using
   ImageJ/Fiji software. Color deconvolution for DAB staining and
   hematoxylin were applied and mean grey values in the positive area were
   measured using the deconvoluted DAB images. Areas of images were
   measured using ImageJ/Fiji software and ratio of positive cell was
   quantified by counting the number of positive cells per CK19^+ cells.
   Colocalization of immunostaining images were measured using the Coloc
   modules in Imaris software.

Library preparation and sequencing for bulk RNA-seq

   MBOs were collected in TRIzol Reagent (Thermo Fisher Scientific,
   15596026), subjected to total RNA extraction using miRNeasy Micro Kit
   (Qiagen, 217084) and verified for integrity by 2100 Bioanalyzer
   (Agilent, G2939BA). Total RNA was poly A selected and
   reverse-transcribed using TruSeq RNA library preparation kit V2
   (Illumina, RS-122). Each sample was fitted with one of the 24 adapters
   containing a different 6 base molecular barcode for high level
   multiplexing. After 12 cycles of PCR amplification, completed libraries
   were sequenced on an Illumina NovaSeq 6000 (Illumina), generating
   around 20 million high quality 125 base-long reads per sample.

Data Processing for bulk RNA-seq

   The quality assessment and pre-processing of MBO RNA-seq reads were
   performed using FASTQC (RRID:SCR_014583, version 0.11.2), Trim Galore
   (RRID:SCR_011847, version 0.4.2) with Cutadapt (RRID:SCR_011841,
   version 1.14), and SAMtools (RRID:SCR_002105, version 1.10.0). Reads
   were then aligned to hg38 genome using Bowtie2 (RRID:SCR_016368,
   version 2.3.3.1). Low-quality alignments and PCR duplicates were
   removed using SAMtools (version 1.10.0) and Picard MarkDuplicates tool
   (RRID:SCR_006525, version 2.18.22). Gene expression was quantified
   using htseq-count (RRID:SCR_011867, version 0.6.1).

Bulk RNA-seq analysis

   Differential expression analysis was performed using Bioconductor
   DESeq2 package (RRID:SCR_015687, version 1.34.0) in R 4.1.2
   (RRID:SCR_001905). Volcano plot visualization of the differential
   expression results was generated using Bioconductor EnhancedVolcano
   package (RRID:SCR_018931, version 1.12.0). Heatmap visualization of
   differentially expressed genes was generated using R pheatmap package
   (RRID:SCR_016418, version 1.0.12). Gene expression was normalized using
   DESeq2 variance stabilizing transformation and z-score scaled. For
   pathway enrichment analyses, overrepresentation of MSigDB Hallmark 2020
   Genesets and GO terms in the upregulated genes in RNA-seq of BA vs. NC
   was tested using Enrichr (RRID:SCR_001575) and ToppGene
   (RRID:SCR_005726)^[248]44,[249]45,[250]78,[251]79. Genes with adjusted
   p value < 0.05 and log2 fold change ≥ 1.00 were selected for the
   pathway enrichment analysis. Expression of EMT related genes in RNA-seq
   of BA-MBO vs NC-MBO was visualized with Boxplot. Gene expression values
   are represented as log2-transformed DESeq2 normalized counts.

Principal component analysis (PCA) of bulk RNA-seq

   Raw sequencing files of RNA-seq of in vivo tissues were downloaded from
   Rimland et al.^[252]32. Reads were pre-processed, aligned, and
   quantified using the same pipelines and reference genome as described
   for the MBO RNA-seq data. Gene expression in NC-MBO and in vivo tissue
   samples were normalized using DESeq2 variance stabilizing
   transformation. Batch differences between in vitro and in vivo samples
   were corrected using the ComBat algorithm in the Bioconductor sva
   package (RRID:SCR_012836, version 3.42.0)^[253]80. Principal component
   analysis (PCA) was performed using the PCA function in the FactoMineR
   package (RRID:SCR_014602, version 2.4) using the top 2000 most highly
   variable genes and visualized using the ggplot2 package (RRID:SCR_
   014601, version 3.3.6).

Library preparation and sequencing for single-cell RNA sequencing (scRNAseq)

   To obtain single cells, MBOs were washed several times with DPBS and
   dissociated using Dissociation Buffer. The dissociated cells were
   passed through a 30 µm cell strainer (pluriSelect, 43-50030-01) and
   resuspended in William’s E medium containing 0.01% Bovine Serum Albumin
   (Sigma-Aldrich, A3059). ScRNAseq was performed according to the
   Chromium Single Cell 3’ Reagent Kit v3.1 protocol (10× Genomics) with
   the Chromium platform (10× Genomics), aiming for 10,000 cells. The
   barcoded cDNA was cleaned with Silane DynaBeads and amplified using 14
   PCR cycles. Full-length, barcoded cDNA was then enzymatically
   fragmented, sized-selected using SPRIselect Reagent, adapter-ligated,
   and amplified for library construction. During the library
   construction, the Illumina R2 primer sequence, paired-end constructs
   with P5 and P7 sequences, and sample indexes were added. The samples
   were pooled and run on the NovaSeq 6000 system sequencer (Illumina)
   using S4 flow cells and paired-end NovaSeq 6000 S4 Reagent Kit v1.5
   (200 cycles) (Illumina, 20028313), aiming for a minimum of 325 M reads
   per sample.

Dataset Processing for scRNAseq

   The raw sequencing data were processed using CellRanger
   (RRID:SCR_023221, version 6.0.0) to generate FASTQ files and counts
   matrices following standard practices. For quality control, each sample
   was subjected to a filter of 500 minimum UMIs detected, 250 minimum
   features detected, and less than 50% of features mapping to the
   mitochondrial genome. This resulted in 27,645 cells remaining for the 8
   samples, of which 4 were NC-MBOs and 4 were BA-MBOs. The gene
   expression data were log normalized and scaled with Seurat (RRID:SCR_
   016341, v4.0.3) according to the authors recommendations. Following
   normalization and scaling, integration was done with Seurat to merge
   the samples into a unified dataset using 2000 integration features.
   Clustering was performed at a resolution of 0.5, which resulted in 15
   original clusters. Dimensionality reduction TSNE and UMAP plots were
   generated for visualization purposes.

Cell classification

   Differential expression was performed to identify genes unique to each
   of the clusters using default parameters for minimum fraction
   expression and log fold change threshold, which indicated the need to
   merge similar clusters due to few unique features. Expression of known
   marker genes from the literature and expert recommendations, together
   with Enrichr and ToppGene-based pathway analysis and scoring cells from
   reference signatures, were used to make the final
   classifications^[254]78,[255]79. Genes with adjusted p value < 0.05 and
   log2 fold change > 1.00 were selected for the pathway enrichment
   analysis. Redundant, mature epithelial and mesenchymal clusters were
   each merged into single clusters, while the small immune-like
   population, the endothelial-like population, and each of the proposed
   epithelial-mesenchymal transition (EMT) populations were left
   separated. Due to the low cell number and low DEG number, one cluster
   was referred as “Other”. This cluster expresses genes related to
   mesenchyme^[256]81, histocompatibility complex^[257]82,[258]83 and
   EMT^[259]81,[260]84. However, CD45 (PTPRC)^[261]85 is not included in
   the DEGs indicating that the cell in this cluster are not immune cells
   but cells having immune-like gene expression. Bar plots displaying
   differential proportions between normal control and biliary atresia
   organoids were generated using the ggplot2 package (version 3.4.4).

EMT and ECM scores

   Identification of clusters along the EMT and ECM production was done by
   examining each cluster for high expression of the top 100
   EMT-associated genes from the EMTome resource and the 275 ECM-related
   genes from core matrisome genes^[262]51,[263]52. To score cells for EMT
   and ECM activity, we adapted a previously described method for combined
   gene expression scoring^[264]86. Briefly, each cell in the data set is
   scored based on the combined normalized expression of all genes in the
   provided signature. To do this, for each gene in the signature, the
   scoring function first defines a set of 100 control genes whose
   expression is similar to the provided gene. The average expression of
   the control genes is then subtracted from the expression of the
   signature gene per cell for normalization. The normalized values for
   all gene are then averaged to get an overall score for each cell in the
   dataset. To generate the UMAP plot, an array of colors is generated,
   ranging from the minimum score per dataset to the global maximum score,
   which allows for accurate visual comparison between datasets of
   different conditions. Code for reproducing this scoring can be found in
   the project GitHub repository in the code availability section.

Cell-cell interaction prediction

   Receptor-ligand interactions were predicted using CellPhoneDB
   (RRID:SCR_017054) with default parameters^[265]87. Chord diagrams were
   used to visualize the interactions between cell types in the NC-MBO and
   BA-MBO organoids separately with the circlize (RRID:SCR_002141, version
   0.4.15) and chorddiag (version 0.1.3) packages in R^[266]88.

Proteomics data

   The data were downloaded from our previously published article
   (Lertudomphonwanit et al.^[267]38). Briefly, serum proteins from a
   total of 175 BA subjects and age-matched normal controls (NC, n = 9)
   were measured by high-throughput proteomic assay (SOMAscan).
   Statistical analyses between the BA and NC groups were conducted using
   Student’s t test with FDR correction (Benjamini-Hochberg procedure)
   using stats (version 4.1.0). Proteins with adjusted p value < 0.1 and
   log2 fold change ≥ 0.585 were selected for pathway enrichment analysis
   using Enrichr and ToppGene. Heatmap were generated using pheatmap
   (version 1.0.12) and employing log2 transformed serum protein levels
   that were enriched in the EMT pathway in MSigDB Hallmark 2020 Gene
   sets. Serum TGF-β/Activin ligand levels were indicated using relative
   fluorescent unit (RFU).

Quantification and statistical analysis

   For organoid culture experiments, quantifications were carried out at
   least by measuring 3 subjects per condition to increase
   reproducibility. For animal experiments, a minimum of three mice were
   included. Unless otherwise noted, GraphPad Prism 9 (RRID: SCR_002798,
   GraphPad Software) was used for all statistical analysis with
   two-tailed, an unpaired t-test for two group comparisons. In cases
   where more than two groups were being compared, one-way ANOVA was
   performed followed by Tukey’s test. p-values of <0.05 were considered
   statistically significant.

   For scRNAseq, the results were analyzed using two-tailed, unpaired
   Welch’s t-test for two group comparisons in parametric data and
   non-parametric data were analyzed using Wilcoxon Rank Sum test. To
   analyze the relationship between EMT score and ECM score, Spearman’s
   Rank Correlation test was used. Differential expression with at least a
   2-fold change (apeglm shrinkage) and false discovery rate <0.05 was
   considered significant in bulk RNA-seq analysis. Data are represented
   as mean with standard deviation unless otherwise noted. In Enrichr and
   ToppGene analyses, the p-values were defined by Fisher’s exact test and
   adjusted as q-values by Benjamini-Hochberg method. In the enrichment
   analysis of bulk-RNAseq, scRNAseq, and microarray data, genes with log2
   fold change greater than or equal to 1 and adjusted p-value less than
   0.05 were considered differentially expressed. In the enrichment
   analysis of proteomics data, to identify 76 DEGs, genes with log2 fold
   change greater than or equal to 0.585 and adjusted p-value less than
   0.1 were considered differentially expressed. Gene sets with less than
   15 or more than 200 genes were not used for ToppGene analysis and the
   gene sets were selected using Enrichr analysis^[268]45,[269]89. The
   combined score is the z-score multiplied by the natural log of the
   p-value^[270]79. The bubble plots were generated using Prism software.
   Details referring to the specific numbers of replicates for each
   experiment can be found in the figure legends.

Reporting summary

   Further information on research design is available in the [271]Nature
   Portfolio Reporting Summary linked to this article.

Supplementary information

   [272]Supplementary Information^ (106MB, pdf)
   [273]41467_2025_61442_MOESM2_ESM.pdf^ (196.8KB, pdf)

   Description of Additional Supplementary Files
   [274]Supplementary Movie 1^ (7.2MB, mp4)
   [275]Supplementary Movie 2^ (43.8MB, mp4)
   [276]Supplementary Movie 3^ (29.7MB, mp4)
   [277]Supplementary Movie 4^ (47.3MB, mp4)
   [278]Supplementary Data 1^ (119.2KB, xlsx)
   [279]Supplementary Data 2^ (354.5KB, xlsx)
   [280]Supplementary Data 3^ (147KB, xlsx)
   [281]Supplementary Data 4^ (139.9KB, xlsx)
   [282]Supplementary Data 5^ (33.1MB, xlsx)
   [283]Reporting Summary^ (204.2KB, pdf)
   [284]Transparent Peer Review file^ (44.1MB, pdf)

Source data

   [285]Source Data^ (255.2KB, xlsx)

Acknowledgements