Hepatocytes undergo punctuated expansion dynamics from a periportal stem cell niche in normal human liver

,


Introduction
The liver possesses an extensive playbook for cell renewal and recovery from injury.Decades of partial hepatectomy studies have established that mature hepatocytes, regardless of their anatomic location, are capable of restoring the vast majority of the liver mass. 1,2Hepatocytes are seemingly inexhaustible in their replicative capacity, without signs of malignant transformation 3,4 and are capable of dedifferentiation, expansion and redifferentiation, both in vitro 5 and in rodent damage models, 6,7 and thus, may themselves be capable of regenerating the biliary epithelium. 8,9Self-renewal and multipotentiality, properties of stem cells, are unusual for a differentiated cell.As such, the necessity for a dedicated hepatic stem cell niche, as is present in many other organs has been questioned.4][15] This has been the source of some controversy, with several lineage-tracing studies failing to detect significant proportions of progenitor-derived hepatocytes upon regeneration from injury. 16,17Evidence is now mounting that widespread impairment of hepatocyte replication is required for activation of the progenitor cell regenerative compartment. 13,14,18This finding is particularly relevant for human disease in which widespread hepatocyte replicative impairment is commonplace.
Whilst there is a reasonable degree of clarity regarding cellular responses upon liver injury, there is still conflicting opinion regarding homeostatic dynamics.Early studies utilised tritiated thymidine to track hepatocytes in normal rat livers over time. 19,20In these studies, hepatocyte migration along the portal tract-central vein axis was observed; a phenomenon coined 'the streaming liver'.Subsequent studies have proposed homeostatic neo-hepatocyte generation from all three zones in the hepatic lobule; periportal, 21 mid-zonal 17,22 and pericentral. 235][26][27] Beyond rodent models, there is evidence supporting the existence of a periportal progenitor niche in the human liver.Huch et al. were able to generate bipotential organoids from EPCAM + ductal cells isolated from normal human livers. 28Similarly, a periportal, bipotent, organoid-forming cell population was identified from a single-cell RNAsequencing analysis of normal human liver cells. 29ur approach to studying homeostatic dynamics in normal human liver centres on demarcation of clonal boundaries within these tissues.In contrast to diseased liver, clones within the typically quiescent normal human liver are often modest in size 30 however, intriguingly, clonal expansions are indeed detectable. 31,32Lineage tracing studies using a visually recognisable deficiency of the mitochondrial enzyme cytochrome c oxidase (CCO), caused by acquisition of somatic mutations, have demonstrated the presence of clonally related cells presumably derived from long-lived stem cells. 32,33tudies utilising this technique in the liver suggested a spatial association between the clonal patches and the PT; however, further evidence that the periportal region is a stem/ progenitor niche in the normal human liver is lacking.Whether potential progenitor cells exist as SOX9 + hybrid hepatocytes, 34 or as a subpopulation of biliary epithelial cells, is currently unknown.
Mitochondrial DNA (mtDNA) mutations clonally expand within cells by random intracellular drift and mitochondria do not require cell division for the generation of somatic mutations. 35Cell division facilitates the spread of mutant mitochondria through tissues and the rate of cell division determines the speed and the physical direction this takes. 36hould cell division be constant then large clones comprised of clonally related cells will be detectable over time; however, if cells within a clonal patch are quiescent, each will generate its own independent mitochondrial mutation repertoire not shared with its neighbours.In this study, we exploit this unique feature of mutation spread, using it to determine the behaviour of clonal expansions, as a proxy for cellular dynamics within the normal human liver.
Herein, we demonstrate, for the first time, the spatial dynamics and the origins of human hepatocytes in vivo.We utilise a range of techniques including CpG methylation sequencing, mitochondrial DNA next-generation sequencing (NGS) and mathematical modelling to establish their origins and the nature of their expansion.Collectively, our data support the existence of a periportal stem cell niche that feeds the clonal expansion of hepatocyte patches in the normal human liver.These patches most likely share a common cell of origin with the biliary epithelium and expand in a slow and punctuated fashion, with a portal to central trajectory, and are largely quiescent, potentially for decades thereafter.

Samples
Human liver samples were obtained by informed consent from patients undergoing hepatic resection.Samples were obtained through the Queen Mary University of London Cancer Tissue Bank (REC reference 14/LO/2031 renewed 19/LO/1700, more information on https://www.cancertissuebank.org, and project approval 2015/01/QM/SM/E/FrozenTissue) and University of Cambridge (16/NI/0196) and University of Padova, Padova, Italy (IRB CESC 3312/AO/14).Additional clinical details can be found in Table S1.

Preparation of tissue
Normal liver, as identified by a pathologist, was taken from each resection, distal to the metastatic or diseased tissue.Resections were cut to size, submerged in isopentane and frozen in liquid nitrogen.Frozen samples were either sectioned onto UV-treated P.A.L.M. membrane slides (Zeiss, Cat # 415190-9041-000) or charged glass slides.Tissues and slides were stored at -80 C until required.

Immunohistochemistry
Dual immunohistochemistry was performed on frozen tissue sections using anti-SOX9 (ABCAM, Cat # ab185966) and anti-KRT19 (DAKO Cat # M0888) according to standard procedures.Anti-Ki67 (ABCAM, Cat # ab16667) immunostaining was performed using the Discovery Ultra Ventana apparatus.A detailed methodology can be found in Supplementary Methods S3.

Enzyme histochemistry
Enzyme histochemistry to identify clonal CCO deficient patches/cells was performed as previously described. 55Enzymatic chromogen development of DAB and NBT for CCO and succinate dehydrogenase (SDH) activity were typically halted after 50 and 25 min, respectively.

CCO patch-structure association and digital cell counting
Serial sections stained with H&E or for CCO/SDH, SOX9/ KRT19 or Ki67, were scanned with a NanoZoomer S210 (Hamamatsu).Slides were opened and analysed using detection tools and scripts within QuPath.See Supplementary Methods S4 for a full description of this methodology.

Mitochondrial copy number analysis
Tissue and DNA were extracted as per the mtDNA NGS methodology below, using the same tissue cut sizes and quantities.The Human Mitochondrial DNA Monitoring Primer Set (Takara Bio, Cat # 7246) was used to determine mtDNA copy number.Analysis was performed as per the manufacturers' instructions.

Mitochondrial sanger sequencing
Single hepatocytes or bulk stroma were subjected to LCM and tissue was lysed according to the Arcturus ® PicoPure ® DNA Extration Kit protocol (Thermo Fisher Scientific, Cat # KIT0103).Amplification of the mitochondrial genome was performed as previously reported. 56PCR products were Sanger sequenced by Eurofins Genomics.Sequence files were aligned against the revised Cambridge reference sequence to identify mutations.Visualisations and comparisons were performed using the "A plasmid Editor" software available at https://jorgensen.biology.utah.edu/wayned/ape/.

DNA methylation
DNA methylation data for the CpG islands in the promoter region of CSX (cardiac-specific homeobox) were generated as described previously. 40

3' RNA seq
Periportal, pericentral and centrilobular regions of the liver were laser-microdissected from CCO/SDH-stained frozen liver sections.Libraries were prepared using the 'fresh-frozen' variant of the Smart-3SEQ protocol, 57 with slight modifications to the protocol including: the use of 7xN unique molecular identifiers (UMI), an 8 nucleotide index for increased plexity, as well as a 30 min lysis time and a 2 lM 1S primer.The libraries were pooled and single-end sequenced on the Illumina NextSeq 2000 platform using a P2 kit for 100 cycles.
Adapter and polyA trimming were performed using Trimmomatic.Sequences were trimmed to 75 nucleotides then, using umi_homopolymer, UMIs were extracted and appended to the read metadata and G-overhangs were discarded.Alignment was performed using STAR and deduplication, using UMI-tools.Differential expression was determined using DESeq2 and pathway enrichment was performed using the fgsea package.

mtDNA NGS
Four to five areas in a portal-central line were laser-capture microdissected from large CCO-deficient patches that spanned from the portal tract to the central vein.DNA was extracted, then PCR amplification of mitochondrial genomes and subsequent library preparation was conducted according to the Human mtDNA Genome Guide (Illumina, 15037958).Libraries were sequenced on the Illumina MiSeq v2 or NextSeq 500 mid output platforms.A full description of the mtDNA NGS library generation and data processing can be found in Supplementary Methods S5.1.

mtDNA NGS data processing
Data quality was assessed by FastQC and reads were aligned to the revised Cambridge Reference Genome (GenBank accession number NC_012920.1).Single-nucleotide variants (SNVs) were identified using the deepSNV R package. 58Variants identified by DeepSNV with <10 supporting reads in the test sample were filtered.Germ-line mutations, defined as those present in both test duplicate samples and in the control at a frequency of >1% were removed.Variants not called independently in both technical replicates were discarded.Public variants were defined as those present in all samples along the portal-central axis of a CCO-deficient patch.Conversely, private variants were defined as those which were called in a subset of samples along the portal-central axis.A full description of the DNA NGS data processing can be found in Supplementary Methods S5.2.

Mathematical modelling and Bayesian inference to integrate model and data
We investigate spatial patterns of mtDNA mutation in a simulated 2-dimensional liver lobule, comprised of periportal "PT" cells and hepatocyte-like "non-PT" cells.We assume three alternative scenarios, where PT cells divide slower, faster and at the same rate compared to non-PT cells.The rates of PT and non-PT cell division are related by r (PT) = b $ r (non−PT) , where b > − 0 is a model parameter.When b = 1, both PT and non-PT cells divide at the same rate and there are no special niches for PT cells.The regime 0 < − b < 1 leads to "non-localised" expansion with the mutation accumulation over space dominated by the non-PT cells.Values of b >1 lead to "localised" dynamics, where the mutation spatial pattern is driven by the expansion of PT cells in the lattice centre, which corresponds to the existence of periportal progenitor niches.
By modelling these alternative hypotheses under a range of b values and simulating the mtDNA mutations in cells located in space, we can predict different possible spatial patterns of mutation accumulation under different hypotheses.Finally, we compare our spatial patterns estimated by multiple simulations with our NGS data and use Bayesian statistics to infer the most likely range of b for each sample.Detailed methodology can be found in Supplementary Methods S1.

Statistical analyses
Two-group individual gene expression statistics were performed with the non-parametric Wilcoxon test using the R package 'ggpubr'.All other two-group analyses were performed within Prism (Graphpad) using the unpaired, nonparametric Mann-Whitney U test.Correlation analyses were statistics calculated using a Pearson correlation test.Gene set enrichment p value statistics were Benjamini-Hochberg adjusted.Statistics relating to computational modelling are described within the supplementary information.

Characterisation of CCO-deficient clonal expansions
We previously established the use of CCO deficiency to investigate dynamics of clonal cell populations in human tissue, including the liver. 31,32,36This method enables histochemical identification of clonal patches of cells that have inherited a lack of CCO activity from a long-lived cell of origin.Herein, we have utilised this technique to visualise CCO-deficient clonal hepatocyte expansions (blue), henceforth referred to as 'patches'.Observations from large-area samples reveal clones extending from the portal tree (Fig. 1A).To characterise the patches in our cohort, we investigated their abundance with age.Detection of CCO deficiency by enzyme histochemistry requires mitochondrial homoplasmy, or at least a high degree of heteroplasmy that takes decades to develop stochastically. 35Indeed, in our cohort, the number of patches observed per mm 2 of tissue increased significantly with patient age (Fig. 1B; p = 0.0017).Despite CCO deficiency, blue patches were previously determined to be metabolically indistinguishable from their adjacent, brown, CCO-proficient counterparts. 31We digitally counted the number of nuclei in blue, and equally sized neighbouring, brown areas and found the former had significantly more nuclei/area (Fig. 1C; p = 0.0067).These nuclei were of similar size indicating a similar ploidy profile (Fig. 1D).Hepatocytes collected by LCM from throughout the lobule had a 2.2-fold greater mtDNA copy number in blue vs. brown regions (Fig. 1E; p = 0.0009).This is a previously reported compensatory mechanism that attempts to restore residual wild-type mtDNA to optimum levels. 37No mtDNA copy number differences were observed when comparing periportal vs. pericentral hepatocytes within brown (Fig. S1 A,B) or blue (Fig. S1 C,D) regions.We also evaluated proliferation by Ki67 immunostaining, calculated as a percentage of total cell nuclei rather than by area due to the increased cell density in blue patches.As expected in normal liver, overall Ki67 positivity was low, with no significant difference in the proportion of proliferative cells in CCO-deficient patches vs. adjacent CCO-proficient regions (Fig. 1F).
To further determine if CCO deficiency impacts hepatocyte biology, we performed Smart-3SEQ (3' RNAseq) on lasercapture microdissected areas taken from CCO-deficient and proficient hepatocytes in periportal, centrilobular and pericentral regions (Fig. S2A,B).After correcting for patient-patient variability, we investigated the impact of patient therapy, the presence of patient metastases (i.e.not within our samples) and the effects of CCO deficiency on gene expression (Fig. S3A-D).It is important to note that our liver samples were obtained distal to any cause of resection (metastases or otherwise).Using principal component and unsupervised heatmap analyses, we found no informative clustering of samples based on therapy or the presence of patient metastases.We did observe a degree of clustering based on CCO deficiency (Fig. S3C,D).
Hepatocyte zonation remained intact (Fig. S4A,B) as we observed a host of previously reported genes with periportal zonation such as HAL, C7, SDS, RAPGEF5, NAV2 and LEPR and pericentral zonation such as SLCO1B3, CYP1A2, CYP3A4, GLUL and SLC1A2. 29,38We compared CCO-deficient with CCO-proficient microdissections using 'fast gene set enrichment analysis' and observed significant negative enrichment of the oxidative phosphorylation pathway in pericentral CCOdeficient samples as expected.A negative periportal enrichment was observed, however, this was not significant.The apoptosis pathway was significantly negatively enriched within CCO-deficient samples of both comparisons, suggesting CCO deficiency is not pro-apoptotic.The reactive oxygen species pathway, a crude indicator for localised tissue damage, was not significantly altered, irrespective of location within the lobule.We additionally compared all (periportal, centrilobular and pericentral) CCO-deficient and proficient microdissections and found both the oxidative phosphorylation and apoptosis pathways to be significantly negatively enriched.The reactive oxygen species pathway was not significantly enriched (Fig. S4C-E).In conclusion, our data broadly reflects previous observations that CCO deficiency does not have major deleterious effects on hepatocyte biology.
We next sought to localise the source of CCO-deficient clonal expansions, particularly as stem cells or committed progenitors are thought to be sufficiently long-lived to acquire homoplasmy or the high level of heteroplasmy required for visualisation as CCO-deficient and thus be the cells of origin of CCO-deficient patches. 31Observationally, CCO-deficient patches appear to emanate from terminal PTs (Fig. 1A).This propensity for periportal patch localisation was proposed in our previous investigations of histologically normal human livers, 31,32 but never quantified.Here, we have determined the closest structure, either PT or central vein (CV), for 377 CCOdeficient patches of 10 cells or greater in 17 normal human livers (Fig. 1G).In these livers, an average of 76.9% of patches were located closer to PTs than CVs, whilst only 21.5% were closer to CVs.A much smaller proportion (1.6%) were sufficiently large to associate with both structures (Fig. 1G).Patch localisation thus favours the presence of a periportal stem/ progenitor niche in phenotypically normal human livers.

Hepatocytes and cholangiocytes share a common somatic ancestor
If CCO-deficient patches originate periportally, the biliary epithelium or nearby progenitor-like hybrid hepatocytes 34 may well be a source.To investigate this possibility, we studied 54 patient samples and observed only three with CCO-deficient biliary epithelium and nearby deficient hepatocyte patches.In one such case, several CCO-deficient patches were found adjacent to terminal PTs that branched from a conducting PT containing the CCO-deficient biliary epithelia (Fig. 2A,B).CCOproficient (Fig. 2B; Patch 1) and deficient hepatocyte patches (Fig. 2B; Patch 2&3) as well as the nearby CCO-deficient biliary epithelium (Fig. 2B) were laser-microdissected and Sangersequenced.Two CCO-deficient biliary areas (Fig. 2C-E) shared a m.8251G>A mutation (in MT-CO2 coding for subunit COII of CCO) that was not present within in the stroma of other PTs from the same tissue (Fig. 2I,J).This mutation was also undetected within CCO-proficient areas of the biliary epithelium and hepatocytes within Patch 1 (Fig. S5).Two CCO-deficient areas tested within Patch 2 shared a m.1214A>G mutation (in MT-RNR1) but did not have the m.8251G>A mutation present within the CCO-deficient biliary epithelium (data not shown).However, the m.8251G>A mutation was present within CCOdeficient hepatocytes within Patch 3 and absent from CCOproficient hepatocytes nearby (Fig. 2F-J).The chances of a common mutation independently arising in these two cell types are extremely low; 2.48×10 9 :1. 32,36We were unable to determine the clonal relationship within other samples with CCOdeficient biliary epithelium (Fig. S6A-D) due to insufficient CCO-deficient cells in serial sections cut for LCM.Regardless, based on these data, the existence of a somatic, periportal, bipotential, common ancestor of biliary epithelial cells and hepatocytes seems a strong likelihood.To our knowledge, this is the first such clonal demonstration in normal human liver.SOX9 + and Ki67 + hepatocytes are present in normal human liver but SOX9 + hepatocytes are predominantly periportal Normal liver is largely quiescent.However, if hepatocyte turnover in this tissue is actively fed from periportal cells, this region may be more proliferative than the parenchyma at the CV.We annotated all PTs and CVs within normal livers from 13 patients using the pathology software QuPath, and we detected Ki67 + cells within 50 lm expansions from the boundaries of these structures (Fig. 3A,B).We observed no significant difference in the proportion of Ki67 + cells at the PT vs. CV (Fig. 3C).We also investigated the possibility of a SOX9 + progenitor as a source of neo-hepatocyte generation in normal liver, irrespective of CCO status. 34As SOX9 is also expressed in KRT19 + ductal cells, we performed dual immunohistochemistry for SOX9 (Brown; DAB) and KRT19 (Blue; Vector Blue AP) to enable their exclusion (Fig. 3D,E).QuPath was trained to differentiate between SOX9 + KRT19 + ductal cells (combined blue & brown immunostaining) and SOX9 + KRT19 -hepatocytes (stained brown only) (Fig. 3E).Blue and brown digital detection masks overlay and highlight these detections within 50 lm expansions surrounding each CV and PT (Fig. 3E).SOX9 + hepatocytes were detected in significantly greater abundance within 50 lm of PTs relative to CVs; a mean of 14.85 vs. 5.89 cells/mm 2 (Fig. 3F; p <0.0001).Collectively, these results are in line with the observations of Font-Burgada. 34

Methylation diversity supports the proposal of a periportal stem cell niche
We next sought to investigate the expansion dynamics of clonal patches in normal human liver.For this, we used methylation diversity as a dynamic clonal markerone that changes with time and size of clonal expansion. 39,40Methylation patterns record cell ancestry as they are inherited with high fidelity by daughter cells.The greater the methylation diversity between cells, the more distant their ancestral relationship. 39We used two markers of diversity: intra-patch distance (calculated as the minimum number of changes required for each sample to have the same sequence) and the number of unique methylation sequences (tags).First, we measured methylation diversity within large periportal and pericentral regions without regard to clonality and found no difference in the epigenetic distance between groups (Fig. S7A-C).Next, CCO-deficient periportal, centrilobular and pericentral cuts were analysed from each patch and no differences were observed (Fig. S7D-F).If the streaming liver hypothesis 19 is correct, and hepatocytes migrate along the portal-central axis, the ancestral and epigenetic distance should increase with anatomical distance from the stem cell niche.Cuts were made within the patches, at varying distances from the PT along the portal-central axis (Fig. S7G-I), but no correlation was observed between epigenetic and anatomical distance.These data indicate the epigenetic diversity measurable using this strategy was maximised.Smaller extractions, remaining close to the putative stem cell niche are more likely to remain within the limits of this diversity measure.To this end, we investigated methylation diversity in patches of increasing size that abutted the portal tract.Measuring diversity is heavily influenced by the number of cells analysed thus micro-dissections were size-controlled.Extractions were made within small, medium and large CCOdeficient patches directly abutting PTs (Fig. 4A,B) and CVs.Larger clones had multiple dissections, each the same size but analysed independently.Methylation tag diversity increased with patch-abuttal size at the PT (Fig. 4C).A significant, positive correlation was observed between epigenetic distance and patch PT abuttal length (Fig. 4D).No such correlation was observed for patches abutting the CV (Fig. 4E).These data demonstrate an ancestral relationship between hepatocytes at the PT but not the CV, where it dissipates due to the number of cell divisions taken for clones to expand to the CV.
We hypothesise there are no differences in expansion dynamics between CCO-deficient and proficient patches.Indeed, we found no significant difference in the % methylation of our CpG sites between CCO-deficient and proficient cells (Fig. S7J) indicating patch dynamics are similar.Given the relationship determined within known clonal boundaries, we also measured methylation diversity in size and spatial location-matched CCO-proficient PT and CV abutting areas.Intra-patch distance vs. increasing abuttal size were non-significant and more weakly correlated than their CCO-deficient counterparts (Fig. S7M,N).This observation demonstrates the benefit of examination within known clonal boundaries, which are otherwise difficult to determine in normal liver.Similar overall results were achieved using another measure of diversity; the simple number of unique methylation patterns or 'tags' observed within a given extraction (Fig. S7C,F,I,J,K,L,O,P).

Clonal hepatocytes display large numbers of unique mtDNA variants across PT-CV axes
Our data thus far indicated the presence of a periportal progenitor niche being the origin of CCO-deficient patches; however, the spatio-temporal dynamics of these expansions could not be determined with the limited resolution of our methylation assay.To address this, we performed deep LCM NGS of mtDNA genomes, a technique with far greater resolution, enabling analysis across the whole patch.In a similar fashion to methylation diversity, mtDNA mutations can be exploited to determine the relatedness of cells across a patch.As mitochondria will continue to replicate in the absence of cell division, mutations can accumulate even in largely quiescent tissues such as the normal liver.
If hepatocytes continuously stream from their cell of origin, we hypothesise that we would detect a distinct stepwise pattern of inheritance of mtDNA variants across each patch.Older cells, hypothetically at the CV, would have had more time to develop independent mitochondrial variants compared to cells at their PT origin.Faster streaming would see an increase in public (clonal) variants detected across the whole patch as there is less time for large proportions of independent variants to develop.Alternatively, should patches grow in a discontinuous or punctuated fashion, with intervening quiescence, we would observe a large degree of independent variants across the whole patch.Four to five 'cuts' of size-matched hepatocytes were laser microdissected along the PT-CV axis in CCO-deficient patches that bridged PTs and CVs (Fig. 5A).After DNA extraction, mtDNA genomes were sequenced by NGS and the variants within each cut were called.In total, 68 cuts across 14 PT-CV axes were extracted from 10 patches across five patient samples.The loci of all PT and CV variants were also plotted for each patch and no discernible differences were observed (Fig. 5B).Individual patch data are shown in Figs.S8-10.The mutation types for cuts adjacent to PTs and CVs were compared to investigate whether differences between these locations could be observed (Fig. 5C).Overall, the base substitutions were very similar, with the T>C substitution predominant in cuts from both locations.This substitution has been associated with aging, 41 a phenomenon coinciding with the development of the large patches used for this analysis (Fig. 1A).
Most PT-CV axes contained a large percentage of unique private variants; (those not shared by neighbouring pairs of cuts) (Fig. 5A).By virtue of their shared CCO deficiency, all cells within a given patch share a common ancestor, yet these data indicate cells in each cut are largely genetically distant from their physical neighbours.We propose that this pattern is suggestive of punctuated expansion of clonal patches, with extensive periods of quiescence.This would allow sufficient time for the cells in each cut to develop variants that are unique when compared to their neighbours (Fig. 5D).
Additionally, we sequenced CCO-proficient areas, some directly adjacent to the CCO-deficient areas previously examined, to determine whether CCO deficiency impacted on the variant composition of CCO-deficient clones.We observed no notable differences in the overall distribution of variants (Fig. S11A).We did, however, observe a change in the variant signature; while T>C variants remained dominant, T>A variants showed an increase in CCO-proficient areas relative to their deficient counterparts (Fig. S11B).Furthermore, we also observed an increase in the overall number of unique variants in CCO-proficient areas (Fig. S11C).This can be explained by multiclonal sampling due to an inability to discern clonal boundaries.Crossing such boundaries also drastically reduced the number of clonal variants detected (Fig. S11D).

Spatial modelling supports neo-hepatocyte generation from a periportal niche
We built a spatial simulation model of liver tissue dynamics to explore scenarios that further test our hypothesis of a periportal niche.We compared all possible mtDNA spatial patterns predicted under different model assumptions with the patterns observed in NGS data (Fig. 5A).We explored two modes of clonal expansion: homeostatic cell turnover dominated by proliferation of randomly distributed 'non-PT' parenchymal cells or, alternatively, faster replication by a small number of stationary 'PT cell' progenitors at the centre of the model (Fig. 6A).These two regimes are described by 'b', with the value of b <1 and b >1, respectively.Cycles of homeostasis are characterised by the death of non-PT cells with subsequent replacement and mitochondrial mutation (Fig. 6B).Our model explored randomly distributed non-PT death as well as a CVcentric cell killing model.In both, cell death leaves open space, allowing for occupation during cell replacement.
Dividing cells displace neighbours into nearby empty spaces, made available by prior cell death, to create space to divide (Fig. 6C).Cell movement in our model is thus passive, generated via our displacement mechanics, with no active migration.A traditional straight-line cell-pushing algorithm ultimately produced very similar dynamics to this adapted algorithm (Fig. S12A).After each round of homeostasis, non-PT cells in discrete zones along a random radial axis are sampled to record their mutational burden (Fig. 6D).As mutational burden depends on the number of cells sampled (N s ) per zone, we used N s = 60 in accordance with the average cell number in each cut of our NGS experiment.Our model also captures the phenomenon of mitochondrial mutations without cell division.We used a literature-based mitochondrial mutation rate of 0.1 per cell division. 42,43Further details are in provided in the materials and methods and Supplementary Methods S1.
Each simulation was run over 100 cycles of homeostasis under various combinations of b and number of PT cells (N PT ).To compare our simulated zonal mtDNA mutation patterns with those along each PT-CV axis obtained in our experimental NGS data (Fig. 5A,B), Approximate Bayesian computation rejection sampling was implemented.This was applied to estimate the joint posterior distribution of the b and N PT parameters that best simulate our experimental NGS data.Approximately 2,000 b and N PT parameter pairs were sampled from the posterior distribution for each experimental sample, and geometric medians were used as a point estimator for the best-suited combination of parameters (Fig. S9-10).
With randomly distributed non-PT death, our experimental samples were best fit by b >1, signifying the presence of small numbers of static PT progenitors that replicate faster than the randomly distributed parenchymal cells (Fig. 6E).A low sensitivity to the N PT parameter in our model meant we did not infer the number of PT cells.With N PT = 10, we compared these results to our CV-centric killing model.All b values investigated resulted in a greater total variant burden compared with the prior, homogeneous non-PT death model.Importantly, both models converge at below b z 5. Given we predict values of 1 < b < 5 for the majority of our experimental NGS samples, this implies a 'PT'-dominant expansion niche is present irrespective of the distribution of cell death (Fig. S12B).Collectively, this modelling supports our hypothesis for the generation of CCOdeficient patches from a confined niche and not solely from randomly distributed hepatocytes.

Clonal patches expand slowly and quiesce
We further tested our expansion dynamics hypotheses by simulating the distribution of mtDNA variants across a simpler, 1-dimensional PT-CV axis (Supplementary Methods S2).The simulation takes place over two phases; first, a clonal expansion occurs either slowly or rapidly to populate the axis.As was determined in our prior model, generation of new hepatocytes by the PT-cell positioned at the simulated PT occurs more rapidly than by randomly distributed hepatocytes along the axis.At the conclusion of clonal expansion, a second phase consisting of either hepatocyte quiescence or streaming occurs (Fig. 7A,B).As the duration between clonal expansion and observation of the patch within the biological samples is unknown, the second phase is simulated across various time periods (1 month, 5 years, 10 years and 20 years).During simulated quiescence, all cells stop dividing; however, mtDNA mutation persists as mitochondria are capable of replication in the absence of cell division.Under streaming dynamics, PT-cell replication persists, and as new hepatocytes are generated, those at the simulated CV are pushed out and lost from the system, akin to an active conveyor belt of hepatocytes.Extensive methodological details and parameters are provided in Supplementary Methods S2.At the conclusion of both phases, the axis is divided into five sections of equal cell number, mimicking the patch cuts in the experimental data (Fig. 5A).Unique private variants are detected within each section and compared in neighbouring pairs, as in our biological data.In our simulations, regardless of phase 2 dynamics, the proportion of unique variants in each pair of sections approaches 1 (100%) as the time since phase 1 (clonal expansion) increases (Fig. 7A,B).Interestingly, quiescent dynamics (Fig. 7A) serve to decrease the proportion of public variants (those shared across all sections) over time whilst a conveyor-like streaming dynamic (Fig. 7B) achieves the opposite in this respect.In all instances, an initial rapid expansion appears to accelerate the trend formed under a slow expansion system.Whilst there is some variability in the unique variant plots obtained from our NGS data (Fig. 7C), most approach 100% unique variants across the patch.This  average NGS data (Fig. 7C; black line & boxplot), we found the simulated scenario which minimised the Euclidean distance to the average NGS data with respect to the unique private variants curve and public variant proportion.A slow expansion with 20 years quiescence was determined to be the best fit for both the unique private variants and public variants data.As such, for the average patch, this modelling lends support to a slow/ punctuated clonal expansion with long periods of quiescence.

Discussion
There is a broad understanding that injury-related expansion dynamics are different to those within the normal homeostatic liver; however, these dynamics are poorly understood.The current evidence is contentious and contradictory, with few studies investigating human livers.Under disease conditions, large clones are known to dominate, but these are smaller in the normal liver. 30As such, investigation of normal liver dynamics requires a means of accurately revealing clonal boundaries.We achieve this by visualising clones that share an inherited CCO deficiency.Here we uncover an informative spatial bias of these clonal expansions and within their confines, we measure gene expression, methylation diversity and mtDNA variants, gaining insights into their ancestry.We pair these findings with mathematical models and Bayesian inference to identify the most likely hepatocyte origins and expansion dynamics within our ostensibly normal human liver samples.Some commentators have been dismissive of the possibility of facultative stem cells, particularly within the biliary epithelial population of the liver 16 .Models of hereditary tyrosinaemia type1 have recorded the massive capacity for expansion of transplanted hepatocytes in such severely damaged livers, 4 while similar large scale hepatocytic repopulation abilities have been noted in livers severely poisoned by a hepatotoxic transgene targeting the liver in albumin-urokinase-type plasminogen activator mice. 44,45More recently, studies employing dietary models of liver damage combined with cell lineage tracing, found no evidence of a significant input of biliaryderived cells to the resultant re-population.These investigations led to the conclusion that hepatocyte selfduplication was the sole mechanism responsible for liver regeneration after injury. 16However, all these models failed to recapitulate the seemingly key role of hepatocyte senescence in determining the nature of the regenerative response. 46At least in liver regeneration after injury, a non-competitive environment is required for biliary cells to contribute to hepatocyte replenishment as seen when Mdm2 (negative regulator of p53) is deleted in hepatocytes, 15 or when p21 is over-expressed in hepatocytes 13 or when b-catenin is deleted specifically in hepatocytes. 47Such findings suggest that under certain circumstances the liver can be considered as a maturational lineage system emanating from the periportal area, similar to the lineage systems prevailing in other epithelia, notably the intestine 48 and epidermis.
These observations beg the question as to whether a similar compartmentalisation operates under more homeostatic conditions, where the requirement for hepatocyte replacement is perhaps less urgent?Indeed, such an organisation was first proposed 30 years ago based upon gradients in cell size, ploidy, growth potential, gene expression and matrix deposition. 49,50Further studies revealed that small cells with stem cell attributes such as clonogenicity resided in the proximal biliary tree (canal of Hering), and might be the origin of biphenotypic transit amplifying hepatoblasts and subsequently hepatocytes. 51,52here is evidence for homeostatic hepatocyte turnover from specific liver zones 17,21,23 as well as from non-zonal, randomly distributed hepatocytes. 24,25In our study, CCO-deficient patch location, SOX9 + hepatocyte proximity to the PT, methylation diversity analysis and spatial modelling all support the existence of a periportal progenitor niche as the origin of clonal hepatocyte expansion.We demonstrate clonal patches have a periportal spatial bias, suggesting this as their origin.As CCO deficiency likely initially develops in long-lived stem/progenitors, 32 this localisation likely designates a periportal progenitor niche.In support of this niche, we demonstrate a clonal, common somatic ancestor of the biliary epithelium and hepatocytes.This finding raises the likelihood of a periportal bipotential progenitor as a source of clonal hepatocyte expansions.We also identify the presence of putative SOX9 + KRT19 - progenitor-like hepatocytes 34,53 in greater abundance immediately surrounding PTs.Using methylation as a tool to determine ancestry, we were able to establish a relatedness of hepatocytes at the PT that was absent at the CV, again indicative of the PT as a site of origin for clonal expansions.Furthermore, the pattern of mtDNA SNVs detected across PT-CV axes within clonal patches best recapitulate our computational modelling scenario of portal progenitor-based expansion, as opposed to a system dominated by proliferation of randomly distributed hepatocytes.
In addition to determining the origins of hepatocyte clonal expansions, our experiments offer insights on the temporal nature of clonal expansions.Changes in methylation are realised only through cell replication over time. 40This technique, in theory, can then also provide information on expansion dynamics.Should a patch expand, and stop, there is insufficient cell replication for the methylation pattern to appreciably alter, thus the diversity of patterns observed across the patch should remain low.Subsequent quiescence, as is generally observed in the normal liver, would ensure that the methylation diversity remains largely unaltered.A continuous streaming scenario would involve a greater amount of replication, in theory allowing differences to be observed between the cell origin and terminus.Unfortunately, our methylation methodology did not have the resolution to distinguish these dynamics.This could be improved using methylation systems that enable a greater number of CpG sites to be sampled within a given target promoter.Regardless, our methylation technique was beneficial for its support of periportal clonal origins as mentioned.
By contrast, mitochondria will replicate in the absence of cell division and develop SNV diversity over time, even within quiescent cells.Patch expansion followed by lengthy quiescence allows even neighbouring clonal cells to largely individualise their mtDNA variant profiles, greatly reducing the proportion of public variants (those present across the whole PT-CV axis).Alternatively, a conveyor-like streaming liver model 19,34,48 serves to greatly increase the proportion of public variants over time.Both these hypotheses are confirmed by our simulations.Within our experimental data, the average patch has a low level of public variants with largely individualised variant pools, which is best matched by our long-term (20 years) quiescence and slow patch expansion simulations.The slow patch expansion scenario translates to a periportal cell replication rate approximately 3 times greater than that of randomly distributed hepatocytes in the 1-dimensional model, a value broadly consistent with the Bayesian estimates from our 2-dimensional modelling.An untested but plausible scenario that may still explain our data would be punctuated streaming with long periods of intervening quiescence.Furthermore, limiting the simulations to 1 dimension will likely accentuate the effects of streaming dynamics.If hepatocytes do exhibit streaming in vivo, they would not always strictly travel along the 1-dimensional PT-CV axis.Last, whilst lattice and off-lattice models of diseased livers have been previously compared, how well a lattice approximation recapitulates the dynamics of a highly ordered normal human liver, remains unexplored.Supplementary Methods S1.3 contains a full discussion of such model limitations.These factors may explain some of the discordance between experimental measurements and the simulated data generated under streaming dynamics, as well as variance observed at an individual patch level.
Questions remain as to the cause of these patch expansions.Whilst we have observed differences in expression, mtDNA copy number and nuclear density in CCO-deficient clones compared to proficient counterparts, these are not expected to cause a selective advantage. 35As such, it is reasonable to assume that these clonal expansions are a liverwide phenomenon; not restricted to regions made visible by CCO/SDH staining.Our methylation data indeed indicate that CCO-proficient and deficient hepatocytes are equally replicative.Additionally, CCO-proficient and deficient hepatocytes are indistinguishable in terms of Albumin, cytochrome P450 1A2 immunostaining and Ki67 positivity. 31Given we are analysing histologically normal tissue, these clonal expansions could arise from homeostatic processes and our data support longterm quiescence as is expected in normal tissue.Whilst histologically normal, our tissues could have been affected by the conditions that necessitated their resection.The majority, but not all of our tissues, were obtained from patients undergoing resection for metastases and not all patients received systemic therapy.Crucially, this allowed us to examine whether these factors had a significant influence on the biology of our hepatic tissues taken distally from the cause of resection.Our gene expression data indicated no appreciable contribution of these potentially confounding factors on the biology of our samples.It does, however, remain possible that clonal patches result from repeated bouts of acute injury from environmental or lifestyle exposures/stimuli; however, chronic injury is usually required for new hepatocytes to extend far beyond their niche.Examining clonal patches from younger donors who have yet to experience repetitive hepatotoxic injury could in theory shed further light on this matter; however, normal tissue from this age cohort is exceptionally rare.Additionally, our lineagetracing methodology, combined with our selection criteria of patches bridging the PT and CV, by its nature, predisposes our analyses to older patients who have more abundant patches, but are unlikely to have had a life free from damage.Indeed, the youngest donor used in our NGS analyses was 62 years of age at the time of resection.
Hepatocyte clonal dynamics in humans is understudied, especially within the normal liver, yet it is essential as a baseline for comparison to disease dynamics.Furthermore, identification and characterisation of a progenitor niche could ultimately be therapeutically beneficial, particularly since hepatocyte transplantation in humans has met with only modest success. 54We have demonstrated the presence of a periportal progenitor niche in phenotypically normal human livers and have provided evidence that neo-hepatocyte expansions spread from this niche.Further work is required to understand what initiates these expansions and how their dynamics may be altered.
Punctuated hepatocyte expansion from a common cell of origin with biliary epitheliumHepatocytesPortal vein f = Number of private mutations t = Time = Mutant cells Hepatic artery

Fig. 1 .
Fig. 1.Hepatocyte clonal expansion characteristics.(A) A low-power image of a CCO/SDH-stained liver section with outlined examples of CCO-deficient cell 'patches' (black border) seemingly emanating from portal tracts (green) that branch from the CPT (green dashed).Central veins are highlighted in red.(B) The number of patches normalised by tissue area vs. patient age at resection; n = 22 patients and 688 patches, ***p < 0.001; testing by Pearson correlation coefficient.(C,D) The number and size of nuclei within CCO-deficient and proficient patches.Data in (C) are from 4 patients, n = 130 patches, **p < 0.01 and from 11 patients for (D).(E) The relative number of mtDNA copies in CCO-deficient and proficient hepatocytes; n = 20 microdissections, ***p < 0.001.(F) The percentage of CCO-deficient and proficient cells that are Ki67 + ; n = 13 patients.(G) The % of patches in a given patient that were closer to PTs, CVs or both structures; n = 17 patients.Statistical tests were performed using Mann-Whitney unless otherwise stated.CCO, cytochrome c oxidase; CPT, conducting portal tract; CV, central vein; mtDNA, mitochondrial DNA; PT, portal tract.

Fig. 2 .
Fig. 2.An ancestral relationship between biliary epithelium and hepatocytes.(A) An H&E and (B) serial CCO/SDH-stained section highlighting three microdissected regions (Patch 1-3; yellow outlines) and a partially CCO-deficient bile duct (blue outline).(C) The partially deficient bile duct with corresponding high-power images pre (D) and post (E) laser microdissection.(F) Patch 3 with corresponding high-power images pre (G) and post (H) laser microdissection.(I,J) Sanger sequencing electropherograms from distal stroma and laser microdissected regions outlined in D&E and G&H.* designates the location of the m.8251G>A mutant.CCO, cytochrome c oxidase; SDH, succinate dehydrogenase.

Fig. 4 .
Fig. 4. Methylation diversity demonstrates a periportal hepatocyte origin.(A,B) Left to right; representative CCO-deficient patches of increasing abuttal length, before (A) and after (B) laser-capture microdissection.Red and yellow outlines highlight the equally sized CCO-deficient and CCO-proficient regions captured.(C) Representative lollipop diagrams for small, medium and large patch-abuttal length.Columns represents a CpG site and row represents a sequenced cell.Solid circles denote CpG methylation.(D,E) Methylation diversity vs. patch abuttal length correlations; n = 15 PTs and 10 CVs.Statistics calculated using Pearson correlation test.CCO, cytochrome c oxidase; CV, central vein; PT, portal tract.

Fig. 5 .
Fig. 5. Next-generation sequencing of CCO-deficient hepatocytes.(A) CCO/SDH immunostained sections are laser microdissected across the PT to CV axis.Samples are sequenced and analysed for SNVs.(B) Distribution of SNVs across mitochondrial genome.Blue and red data points represent PT and CV-abutting samples respectively.Outer histogram displays combined SNV data for all PT and CV samples (100 bp bin width).(C) SNV distributions in PT and CV samples.(D) Proposed model of hepatocyte clonal expansion in normal human liver.A patch of clonal hepatocytes expands, then quiesces for extended periods.Over time, cells acquire unique mtDNA variants (different cell colours) independently of cell division.CV, central vein; mtDNA, mitochondrial DNA; PT, portal tract; SNV, single nucleotide variant.

Fig. 6 .
Fig. 6.Spatial modelling of homeostatic liver dynamics indicates the presence of a stem cell niche.(A) Low (< − 1) b regime dominated by non-localised non-PT dynamics.Regimes of b >1 denote PT-cell dynamics localised at the centre of system.(B) Simulated lobule dynamics.After initialising the lobule, all cells carry zero mutations (single colour).Following time-steps of random killing and repopulation, cells develop mutations (different colours).(C) Cell pushing algorithm.Dividing cells (star) push neighbours towards a nearby empty lattice point to create an adjacent empty space, then divide and occupy the space with its progeny.(D) Number of unique variants measured along random radii outwards from system's centre.(E) Estimated b values for each portal-central axis assayed by mitochondrial NGS.Error bars are SEM.PT, portal tract.
is most closely represented by longer phase 2 (5-20 years) simulations.Our biological data has public variant proportions largely <10%, which is most consistent with quiescent phase 2 dynamics.To determine the simulation that most closely resembled our

Fig. 7 .
Fig. 7. Simulated and experimental public and unique private variant dynamics in spatially aligned microdissected areas of CCO clones.(A,B) Unique private variants across sections spanning 1-dimensional simulated portal-central axes.(A) Quiescent and (B) streaming liver dynamics assuming either slow or rapid initial patch growth.Each panel displays data from 1,000 simulated portal-central axes (transparent curves).Solid curves represent mean values across all simulations.Clonal hepatocytes are first generated at either a slow or rapid rate, followed by 1 month, 5, 10 and 20 years of quiescent or streaming dynamics.Proportion of public variants (those shared across all cuts) is displayed as a bar alongside each panel (mean±SD).(C) Data for non-simulated experimental samples and a table of public variant percentages for each PT-CV axis.The solid curve and boxplot (mean ± SD) represent the averages of all pateches.CCO, cytochrome c oxidase; CV, central vein; PT, portal tract.
Affiliations 1 Centre for Cancer Genomics and Computational Biology, Barts Cancer Institute, Queen Mary University of London, London, UK; 2 School of Mathematical Sciences, Queen Mary University of London, London, UK; 3 Centre for Tumour Biology, Barts Cancer Institute, Queen Mary University of London, London, UK; 4 Department of Surgery, Oncology and Gastroenterology, University of Padova, Padova, Italy; 5 Cancer Research UK Cambridge Institute, University of Cambridge, Cambridge, UK; 6 Histopathology, Barts Health NHS Trust, London, UK; 7 Department of Cellular Pathology, University College London, London, UK; 8 UCL Cancer Centre, University College London, London, UK; 9 Cancer Tissue Bank, Barts Cancer Institute, Queen Mary University of London, London, UK; 10 Barts and the London HPB Centre, The Royal London Hospital, Barts Health NHS Trust, London, UK; 11 Group of Theoretical Biology, The State Key Laboratory of Biocontrol, School of Life Science, Sun Yat-sen University, Guangzhou, China