Patient Samples

Cholangiocytes were obtained from the bile ducts of patients undergoing ERCP. PSC and control patients were recruited from the Section of Gastroenterology at the Department of Transplantation Medicine, Oslo University Hospital Rikshospitalet. Written informed consent was obtained prior to ERCP and ethics approval for the use of cells from human livers and bile ducts was approved by the Regional Ethics Committee (2012-286/2016-1540). Patient inclusion in the study was done independent of gender.

Isolation of Cholangiocytes from Patient Bile Ducts

ERCP brush samples were placed into William’s E + (WE +) medium (Gibco Inc., Waltham, Massachusetts, USA) with 50 ng/ml of epidermal growth factor (EGF) (R&D Systems Inc., Minneapolis, Minnesota, USA) and 10 μM Rho-associated kinase inhibitor Y-27632 (Selleck Chemicals LLC., Houston, Texas, USA) on ice and carefully washed to collect patient material. Isolated ERCP material was centrifuged at 400 g for 4 min, supernatant was discarded, and the cells were seeded onto 24-well cell culture plates in 50 μl droplets of supplemented WE + culture medium [21] and Matrigel (Corning Inc., New York, New York, USA) mixed at a 1:3 ratio for in vitro 3D expansion [17].

Expansion and Maintenance of Cholangiocyte Organoid Lines

Cell culture of COs was performed as previously described [17, 21]. Briefly, cells were cultured with a solubilized basement membrane matrix (Matrigel) and WE + medium, supplemented with a combination of EGF (Bio-Techne Co., Minneapolis, Minnesota, USA), R-spondin 1 (Bio-Techne Co.), and Dickkopf-related protein 1 (DKK-1) (Bio-Techne Co.) to promote the expansion of cholangiocytes in the form of organoids [21]. Organoid media was exchanged every 48 h and the organoid lines were split every 5–7 days depending on organoid quality and quantity. After split 4, all organoid lines were transferred to liquid nitrogen and stored for scRNAseq and further analysis. For scRNAseq, all organoid lines were simultaneously revived from liquid nitrogen storage and synchronized in culture. All experiments were performed using passage 6 organoids unless otherwise stated.

Immunofluorescence Staining of Cholangiocyte Organoids

For staining, organoids were fixed with 4% paraformaldehyde (PFA) in PBS for 20 min. Fixation solution was removed and the organoids were washed twice in PBS for 10 min each. Permeabilization and blocking were performed by adding 0.1% Triton X-100 and 5% donkey serum in PBS for 1 h. Organoids were stained with primary antibodies (Supplementary Table 1) overnight at 4 °C. Samples were washed 3 times in PBS for 45 min each and secondary antibody staining (Supplementary Table 1) was performed for 1 h at room temperature, followed by incubation with Hoechst 33258 in PBS for 20 min. Organoids were washed 3 times in PBS for 45 min each and then imaged using fluorescence microscopy. All immunofluorescence images were acquired using a Zeiss Axiovert 200 M inverted microscope (Zeiss Group, Oberkochen, Germany) or a Zeiss LSM 700 (LSM 710) confocal microscope. The ImageJ 1.48 k software [22] was used for image processing. Immunofluorescence images are representative of at least 3 different CO lines.

Flow Cytometry Analysis

COs were dissociated from Matrigel using Cell Recovery Solution (Corning Inc.) for 30 min at 4 °C and were then centrifuged at 444 g for 4 min and dissociated to single cells using TrypLE Express (Gibco Inc.). Cells were stained for 30 min on ice, washed twice in flow buffer (PBS containing 5% fetal bovine serum and 0.09% sodium azide), and measured and analyzed using a BD FACSVerse Flow Cytometer and the BD FACSuite software (BD Life Sciences comp., Franklin Lakes, New Jersey, USA). A complete list of the antibodies and dilutions used is provided in Supplementary Table 2.

γ-Glutamyl Transferase Activity

γ-glutamyl transferase (GGT) activity was measured in triplicate using the MaxDiscovery GGT Enzymatic Assay Kit (BioScientific, Avondale, Arizona, USA) based on the manufacturer's instructions. Absorbance units refer to light absorbance at a wavelength of 405 nm.

Alkaline Phosphatase Staining

Alkaline phosphatase (ALP) staining was performed using the BCIP/NBT Color Development Substrate (5-bromo-4-chloro-3-indolyl-phosphate/nitro blue tetrazolium) (Promega Co., Madison, Wisconsin, USA) according to the manufacturer’s instructions.

REAL-TIME Quantitative Polymerase Chain Reaction

RNA extraction and quantitative polymerase chain reaction (qPCR) were performed as previously described [23]. A complete list of the primers used is provided in Supplementary Table 3. All qPCR data are presented as the median, interquartile range (IQR), and range (minimum to maximum) of four independent CO lines unless otherwise stated. Values were normalized to the housekeeping gene hydroxymethylbilane synthase (HMBS).

Rhodamine 123 Assay

Organoids were initially incubated for 30 min in the presence or absence of 50 µM verapamil (Merck and Co Inc., New York, New York, USA). Organoids were then incubated with 100 µM rhodamine 123 (Merck and Co Inc.) and 1 µg/ml Hoechst 33342 in supplemented WE + culture medium for 5 min. Organoids were washed 3 times in WE + medium and then incubated for 40 min in fresh supplemented WE + culture medium. Organoids were imaged using a Zeiss LSM 710 confocal microscope. Images were analyzed using ImageJ [22] with a mean intensity measurement taken for the interior of each organoid. The experiment was run in triplicate with 40 organoids quantified in total per condition.

Cytokine-Stimulation Experiments

COs were stimulated with 50 ng/ml tumor necrosis factor α (TNF-α) (Bio-Techne Co.), transforming growth factor β (TGF-β) (Bio-Techne Co.) or interleukin 17 (IL-17) (Bio-Techne Co.) in supplemented WE + culture medium for 5 days. Medium and stimulation cytokines were exchanged every 48 h. Senescence-associated-β-galactosidase (SA-β-Gal) stain assay (Abcam Limited, Cambridge, UK) of unstimulated and stimulated organoids was performed following the manufacturer’s instructions. Cell culture supernatant was analyzed by Luminex Discovery Assay (R&D Systems Inc.), following the manufacturer’s instructions. The experiment was repeated in n > 3 independent stimulation experiments with individual patients.

Single-Cell RNA Sequencing and Data Analysis

COs were harvested after the sixth passage and processed for scRNAseq using the 10X Genomics scRNAseq platform and Chromium Next GEM Single Cell 3’ Gel Bead Kit v3.1 (10X Genomics Inc., Pleasanton, California, USA). Raw sequencing reads were delivered in FASTQ format.

Initial quality checks, alignment, and transcript quantification were performed using Cell Ranger version 3.1.0. The alignment was performed against the GRCh38-3.0.0 H. sapiens reference genome; a total 41,954 cells were initially called on default acceptance thresholds.

Further quality control checks and filtering were performed, including cells with fewer than 2500 features (genes) were discarded; the threshold was determined based on the nfeature distribution (Supplementary Figs. 1, 2, 3). In addition, cells with > 15% reads mapping to mitochondrial genes, > 40% reads mapping to ribosomal genes, or > 100,000 unique molecular identifiers (UMIs) were removed; the latter corresponded to doublets (Supplementary Figs. 2, 3). Once the cells were filtered, mitochondrial, ribosomal, and Y chromosome genes were excluded from the expression matrix, prior to normalization. Cell-cycle was assessed using standard sets of cell-cycle genes, using dedicated Seurat functions; no significant batch effect was noted; the cell-cycle genes were subsequently excluded from expression matrix. 15,564 cells passed all the quality control filters, across samples; post-filtering, the expression levels were normalized using SCTransform [24].

To evaluate the clustering stability and optimize hyper-parameters, we employed ClustAssess [25]; parameters such as the feature set (most abundant or highly variable genes, and the number of genes not affected by random noise variation), the resolution, and number of clusters were all assessed. The tested number of genes were 500, 1000, 1500, and 2000, in line with the nfeature distributions across samples. The stability is evaluated using Element-Centric similarity; the consistency of outputs on random iterations was summarized using Element-Centric Consistency (ECC) [26], on 100 iterations on distinct random seeds. We identified the transition point from signal to noise at the 1500 most variable features, which we subsequently used for downstream clustering (Supplementary Figs. 4, 5).

Clustering was performed with Seurat v4.3.0. Initially, a Principal Component Analysis (PCA) was performed on the 1500 most variable features, retaining the first 30 principal components. The Louvain clustering algorithm was applied on the PCs, with a resolution of 0.8, revealing 13 distinct clusters. Differential Expression (DE) analysis was performed using the ROC test, considering genes expressed in > 25% of the cells, and a log2FC > 0.25.

scRNA data were visualized using the ShinyCell package from R. Heatmap generation for group comparisons was performed using the webtool ClustVis [27].

Spatial Transcriptomics and Sequencing

Spatial transcriptomics and sequencing have previously been described for the dataset used [28]. Biliary regions were defined by expression of the cholangiocyte markers KRT7 or EPCAM (KRT7 > 2 counts per spot, EPCAM > 1 counts per spot).

Gene expression levels were visualized within pre-defined biliary regions in Loupe Browser software (10X Genomics Inc.).

Statistical Analysis

All statistical analysis for non-scRNAseq data was performed using GraphPad Prism 9 (GraphPad Software, La Jolla, California, USA). Double sided Student’s t test was used to evaluate statistical significance. P values below 0.05 were considered significant. The number of replicates refers to organoid lines unless otherwise stated.