Animal experimentationAll animal experimentation was done according to local laws, and approved by the veterinary office of the Canton of Vaud, Switzerland. Environmental and biological variability were mitigated by a tightly controlled environment in our animal housing facility and by using inbred strains and regular veterinary checks, respectively. Mice were housed under a 12 h light-dark cycle and temperature of 21 ± 2 °C, with ad libitum access to food (Safe-Lab SAFE 150) and water. The time of a perturbation is denoted with Zeitgeber time (ZT), in accordance with the nomenclature in chronobiology, and ZT0 and ZT12 refers to the time of the switching on or off of lights, respectively. Experiments were performed with 129S6/SvEv (Taconic) mice, abbreviated as 129S in the manuscript. Additional and comparative data from C57BL/6 and CD-1 strains came from previous work, as described in corresponding method and result sections.Handling of human tissueAll human tissue handling, storage and human experimentation was performed according to local laws, and the work was approved by the Cantonal commission of Vaud for ethics in human research. All participants signed a written consent form. The tissue (Supplementary Table 2), assessed by a clinical gastrointestinal pathologist (A.M.), was obtained from the Tissue Biobank Bern (TBB) of the University of Bern, with the criteria to exclude steatotic tissues. The TBB sectioned FFPE at 4 µm and flash frozen samples at 10 µm thickness. Further manipulation of the tissue was performed in house, except H&E staining and scanning, which was done at the University of Bern. Tissue donors were not financially compensated.Analysis of mouse bulk sequencing dataC57BL/6 RNA-seq data was reanalysed from Weger et al.7. Differential expression analysis was previously performed using DESeq275.Analysis of human mRNA data (GTEx RNA-seq)Raw gene read count data was obtained from GTEx V8 (2017-06-05_v8_RNASeQCv1.1.9, dbGaP Accession phs000424.v8.p2). For each tissue, genes with fewer than 10 counts across the samples were excluded. Counts were normalized based on library size and scaled using the edgeR TMM method76, followed by conversion to counts per million (CPM). After the addition of a pseudo-count of 0.25, CPM values were log-transformed (log2). Information on sex, age, and BMI was retrieved from the dbGap metadata (phs000424.GTEx.v8.p2.c1.GRU) for every donor. The history of heart attack and acute coronary syndrome was taken from the metadata category MHHRTATT. Differential expression was assessed between age, sex and BMI groups using a two-sided t-test, and the calculated p-values are reported above corresponding box plots. The p-values that are lower than 0.05 are marked with an asterisk (*), p-values lower than 0.01 with two asterisks (**) and p-values lower than 0.001 with three asterisks (***). The reported p-values were not corrected for multiple testing. We selected an age cutoff of 49 to distinguish between pre- and post-menopausal women to align with the age bins in the publicly available GTEx dataset. Since the health status of individual donors in the postmortem GTEx cohort could confound our analysis, a stratified analysis based on age groups and health conditions was performed to assess whether the effects observed are consistent across different subsets of the population.Analysis of human proteomics dataHuman proteomics data was reanalysed from Jiang et al.61. The normalized protein abundances were used to calculate the mean values and confidence intervals within individual organs.Isolation and cryopreservation of cells for single-cell sequencingFor perfusions, a peristaltic pump was used to withdraw buffered solutions from a water bath, kept at 42 °C, via a system of medical grade silicone tubing. The tubing led to a Graham condenser, which was heated with a circulating water flow from the water bath, utilizing a separate pump. The solutions exited the Graham condenser vial medical grade silicone tubing, which was connected via Luer adapters to a butterfly needle. The system ensured that all solutions had an approximate temperature of 37 °C at the tip of the butterfly needle. To perform the perfusion, mice aged approximately 12 weeks were euthanized with 150 mg per kg pentobarbital. Upon absence of reflexes, the cadavers were opened and the butterfly needle, connected to the peristaltic pump system, was inserted into the inferior vena cava and clamped into place with a surgical hemostat. The perfusion pump was started and the portal vein was immediately cut to prevent an increase of pressure within the liver. First, to flush the blood, the perfusion was performed for 3 min using a buffer without collagenase. Subsequently, the mice were perfused for 3 min with a solution containing collagenase. Intact livers were isolated, placed into approx. 15–20 ml DMEM: F12, and cells released by breaking the Glisson’s capsule. Next, straining through a 70 µm mesh and enrichment for hepatocytes based on centrifugation (50 G for 5 min, with a low breaking speed) was performed. For the centrifugation, each sample was split in two: one falcon was immediately stored as described below. The other falcon was used for resuspension of the pellet, and cell counting and viability assessment (trypan blue), followed by another round of centrifugation. After isolation, we implemented a isopentane cell freezing protocol to preserve cells isolated at the two precise time-points, allowing us to prepare sequencing libraries at a later point. Briefly, after centrifugation, the whole falcon tube containing the pellet enriched for hepatocytes was dipped into isopentane, precooled to −80 °C. The cell viability assessment upon thawing (see below) confirmed the feasibility of the freezing protocol.Preparation of libraries for single-cell sequencingThe cryopreservation reduced hands-on time significantly, reducing technical variability. It also enabled us to thaw the cells immediately before loading onto the microfluidics chip of the scRNA-seq kit, reducing cell sticking and appearance of doublets in the dataset. Briefly, cells stored in 50 ml falcon tubes were placed on ice, and DMEM: F12 (1: 1) medium with 10% FBS, cooled at 4 °C, was immediately added to a final concentration of 2 × 106 cells per ml. Within 5 min, the cells were gently resuspended with a 50 ml serological pipette. The addition of 10% FBS to the medium reduced cell aggregation and the larger serological pipette reduced shear stress on the cells. The cell concentration and viability was assessed and the cells were immediately mixed with the 10x Genomics reagents and loaded on a pre-prepared (glycerol already in empty rows) 10x Genomics Chromium chip. Library preparation was performed using the 10x Genomics Chromium v3 kit according to manufacturer’s instructions, and the quality of the library preparation was monitored with the TapeStation 4200 or Fragment Analyzer instruments at all steps. The assessment of viability post freezing showed virtually no loss of cells.Filtering of mouse scRNA-seq data and clustering of cellsRaw sequencing data were aligned to the mouse genome (mm10) using CellRanger (6.1.2). To minimize the confounding effects of low-quality cells, which can introduce technical noise and obscure true biological signals, we excluded cells that did not meet the following quality criteria. Intron and exon counts were computed using Velocyto77. From all sequenced GEMs (gel beads in emulsion), cells were selected on a per-sample basis by evaluation of the intronic fraction of individual GEMs (refer to78). GEMs with no or low introns were removed for each sample. A cutoff between a low and mid intron distribution was set on a per-sample basis, but at 8.5% for most samples. As hepatocytes are highly metabolically active with a high percent of mitochondrial genes in the transcriptome79, and nonparenchymal cells have lower mitochondrial fractions, a filter of < 40% mitochondrial reads was set. Initial analysis was performed using the Seurat 4.3.0 package. Briefly, samples were normalized (sctransform) and integrated (features = 3000) and the cells clustered based on a t-distributed stochastic neighbor embedding (t-SNE) dimensionality reduction and Louvain algorithm (resolution = 0.5). This approach allowed us to mitigate batch effects and accurately identify individual cell types within the integrated samples. For quantifications of gene expression in identified clusters of individual celly types, non-integrated transcripts (unique molecular identifiers – UMIs) were normalized to the total cytoplasmic gene count (excluding mitochondrial genes) since the mitochondrial fraction changes in function of cell type and lobular zone in hepatocytes. We identified a subcluster of hepatocytes that separated from others by low expression of secreted proteins and which exhibited a low mitochondrial fraction. To split the hepatocytes into those with low and high mitochondrial fractions, a cutoff of 3% mitochondrial genes was chosen. For the high mitochondrial gene subpopulation of hepatocytes (≥ 3%), a lobular position was calculated (spatial reconstruction), as described below. The pseudobulk calculations were done on all samples, with 3 biological replicates per condition. For spatial reconstruction and subsequent plotting of expression profiles, one male ZT22 sample was removed due to overall lower UMI counts.Spatial reconstruction of hepatocyte zonation profiles from scRNA-seqWe developed a probabilistic model to assign a one-dimensional spatial position to every cell based on the expression patterns of established zonation marker genes in individual cells3,25,26. The initial set of markers, taken from ref. 25, included Pck1, Aldh1b1, Ctsc, Sds, Hal, Hsd17b13, Cyp2f2 as portal and Oat, Cyp2e1, Lect2, Cyp2c37, Gulo, Cyp2a5, Glul, Aldh1a1, Cyp1a2, Slc22a1, Slc1a2 as central markers. The spatial latent lobule coordinate x was determined by maximizing the likelihood of the observed raw counts, given x and other model parameters. We used a negative binomial noise model for the read counts to best capture the fluctuations of scRNA-seq data75. The expectation of the negative binomial distribution for a cell \(c\) and gene \(g\) in sample \(s\) was modeled by \({\mu }_{{gcs}}={N}_{{cs}}\exp ({a}_{g}{x}_{{cs}}+{b}_{{gs}})\), where \({a}_{g}\) and \({b}_{{gs}}\) are gene-specific slope and intercept coefficients. Note that the gene slopes \({a}_{g}\) are the same across all samples, while intercepts \({b}_{g,s}\) can be sample dependent. Thus we assumed that there are landmark zonation genes with invariant slopes across all the samples, and these define the lobular coordinates. Finally, \({N}_{{cs}}\) is the number of unique molecular identifiers (UMIs) captured for cell \(c\) in sample \(s\). The dispersion parameter \(\alpha\) was also fitted and it is unique across all samples and cells. The fitted values of the dispersion parameter (0.076 for the initial gene set and 0.096 for the enriched one) are quite low indicating minimal overdispersion, closely approximating a Poisson model.This model displays inherent symmetries without further constraints on the parameter, in fact the model has rank-deficiency \(s+1\). One symmetry is a scaling invariance as the expected values \({\mu }_{{gcs}}\) only depends on the product \({a}_{g}{x}_{{cs}}\). We broke this symmetry by setting the slope \({a}_{g}\) of a selected gene (Cyp2e1) to -1 (which also sets the directionality of the central to portal axis, i.e. portal genes have a positive slope). The remaining \(s\) symmetries concern translational invariance: the shifting of all cells in a sample by the same quantity \({q}_{s}\) redefines the intercept \(b^{\prime}_{{gs}}={b}_{{gs}}+{a}_{g}{q}_{s}\). To set the reference, we fixed \({q}_{s}\) such that cells expressing the highest central marker genes (e.g. Glul) were aligned across all the samples.We initialized cell positions \({x}_{{cs}}\) using principal component analysis (PCA). Specifically, the projection on the first principal component was used. This choice was justified since it was clear that the first principal component correctly separated central and portally expressed genes. The initial gene parameters, \({a}_{g}\) and \({b}_{{gs}}\) were fitted to the PCA positions through negative binomial regression. Once we obtained these complete sets of initial values, we used gradient descent to optimize positions and gene parameters together.In a second step, once the cell coordinates were estimated as described above, we expanded our set of reference zonation markers to enhance the robustness of the spatial assignment. Initially, we chose a subset of 400 genes that showed the most significant deviations from a linear trend in the coefficient of variation versus mean plot, as the initial set of marker genes exhibited substantial divergence from this trend. Next, we fitted the model’s slopes and intercepts for every gene, while keeping the cell positions x previously fitted with the smaller/initial set of genes fixed. Large positive or negative slopes, corresponding to central or portal genes, allowed us to identify potential zonation markers. Here to account for inter-mouse variability, we let the slope coefficient be sample dependent (\({\mu }_{{gcs}}={N}_{{cs}}\exp ({a}_{{gs}}{x}_{{cs}}+{b}_{{gs}})\)) which was not the case in the first model. Finally, we chose genes that had the highest slope’s mean to variance ratio (highest signal to noise) across the different mice. The final enriched set of genes included Ctsc, Aldh1b1, Pigr, Serpina1d, Apoc2, Hal, Gc, Hsd17b13, Itih4, Apoc3, Hpx, Hsd17b6, Apoa4, Sds, Serpina1b, Serpinc1, Alb, Cyp2f2, Pck1 as portal and Rnase4, Glul, Oat, Cyp2c37, Airn, Lect2, Slc22a1, Cyp2e1, Cyp7a1, Gm30117, Aldh1a1, Cyp2c50, Lhpp, Cyp1a2, Slc1a2, Cyp2a5, Pon1, Slco1b2, Cyp2c54, Mgst1, Gulo, Lgr5 as central markers.The plotted expression profiles (Figs. 1–2) were generated by binning the 1D coordinates into 12 bins, representing layers of hepatocytes around the central vein. The expression is represented as counts per million (CPM), normalized by the total gene count excluding mitochondrial genes. Calculating the spatial expression profiles in the pseudobulk space allowed us to mitigate the noise inherent in scRNA-seq data. The expression profiles can be generated using data from S. Data 1.Model selection on mouse scRNA-seq dataAfter spatial inference of hepatocytes, the spatial expression profiles were further characterized by a model selection approach. Genes with > 300 total reads in at least one condition were retained. For these 6235 genes, we employed the ‘glmmTMB’ R package80 to fit a generalized linear mixed-effects (GLMM) model (with a log link function), accounting for negative binomial noise and utilizing the Broyden–Fletcher–Goldfarb–Shanno optimization method. To address numeric instability for genes not expressed in at least one condition (less than 10 total counts), we performed an imputation using random draws from a negative binomial distribution with mean μ = e−4 and dispersion 0.1. Spatial coordinates across all samples were normalized to the range of 0 to 1 for consistent comparisons. Ten models, with incrementally increased complexity, were considered for each gene. These models included fixed effects, with intercepts and slopes being either sex, or time-specific, along with possible interaction terms, and random effects for the intercept and slope in each animal. The library size, computed by summing all gene counts except mitochondrial genes, served as an offset in the GLMM. A Bayesian Information Criterion (BIC) was employed, assigning genes to the model exhibiting the lowest BIC. To optimize computational efficiency, we parallelized the analysis at the gene level using the ‘mclapply’ function from the parallel R package. We kept the estimated fixed effect coefficients from the best-fitting model for each gene. Note that the current simplest model with linear dependency in the positions was not designed to capture genes with peaks in the midlobular region, of which previous analysis identified only very few26. Moreover, the assignment to complex models with interacting effects of sex and time on the shape of zonation profiles did not yield many genes (15 genes in models 9 and 10, BIC, S. Data 2), and for those models the assignments seem sensitive to noise. Therefore we did not consider those genes further.Fenofibrate injectionsFenofibrate was injected intraperitoneally at 80 mg per kg body mass for 5 straight days, at ZT0. Due to its low solubility and instability in water, fenofibrate was dissolved in DMSO in 10x working concentration and frozen at −20 °C. Upon thawing of the solution, medical grade 0.9% NaCl was added to the final volume, resulting in a suspension (rather than solution), which was shaken before loading into the syringe and injecting. Tissue was isolated at two time points that followed the 5th injection, ZT10 and ZT22.Collection and freezing of mouse tissueTissue for histology was isolated from approximately 12-week-old mice, placed into cryo molds, covered with 2% carboxymethylcellulose (CMC) in distilled water (prepared by sonication of CMC in water), and placed into a deep section of dry ice, allowing for even and quick freezing. Alternatively, the tissue was frozen without CMC in the same manner. The tissue was subsequently cryosectioned at 10 µm, placed on glass slides and stored at −80 °C until processing.Immunofluorescence of global zonation markersTriple immunofluorescence was performed manually on mouse liver 10 µm fresh frozen sections. Tissues were fixed with 4% PFA for 15 min at RT before permeabilization with 0.1% Tween in 1x PBS for 10 min. After 30 min of blocking with 1% BSA in 1x PBS, primary antibodies were incubated simultaneously overnight at 4 °C: rabbit anti-CYP2E1 (Abcam ab28146, 1: 300), mouse anti-GS/GLUL (Santa Cruz sc-74430, clone E-4, 1: 100) and goat anti-E-CAD (R&D systems AF478, 1: 200). After washes with 1x PBS, secondary antibodies were applied simultaneously for 40 min at RT: donkey anti-rabbit Alexa568 (Thermo Fisher A10042, 1: 1000), donkey anti-mouse Alexa647 (Thermo Fisher A31571, 1: 1000) and donkey anti-goat Alexa488 (Thermo Fisher A11055, 1: 1000). Sections were counterstained with DAPI and mounted with FluoromountG (Bioconcept 0100-01).VLDLR immunohistochemistry (mouse sections)The immunohistochemistry was performed using the fully automated Ventana Discovery ULTRA (Roche Diagnostics, Rotkreuz, Switzerland). All steps were performed on the equipment with Ventana solutions except if mentioned. Briefly, 10 µm mouse liver fresh frozen sections were fixed with 4% PFA for 15 min at RT, washed with 1x PBS before loading wet on the automated device. The primary antibody, goat anti-VLDLR (R&D Systems AF2258, diluted 1: 50 in 1% BSA in 1x PBS) was incubated 1 h at 37 °C. After incubation with goat Immpress HRP (Ready to use, Vector Laboratories), chromogenic revelation was performed with ChromoMap DAB kit. Sections were counterstained with Harris hematoxylin and permanently mounted. The specificity of the antibody was evaluated in a western blot (Supplementary Fig. 26).VLDLR immunohistochemistry (human sections)Human tissue was processed using the Ventana Discovery ULTRA machine with Ventana solutions except if mentioned. Dewaxed and rehydrated 4 µm paraffin sections were pretreated with heat using the CC1 solution for 40 minutes at 95 °C. For DAB staining, the primary antibody, goat anti-VLDLR (R&D Systems AF2258), was used in a dilution of 1: 25 (1% BSA in 1x PBS), 1 h incubation at 37 °C. After incubation with goat Immpress HRP (Ready to use, Vector Laboratories), chromogenic revelation was performed either with the ChromoMap DAB kit or with the DISCOVERY Purple Kit. Sections were counterstained with Harris hematoxylin and permanently mounted. Given that the VLDLR antibody is not validated for clinical use, we further assessed its specificity in a western blot and in a human control tissue, the pancreas81 (Supplementary Fig. 26). Given a higher background stain in the two needle biopsies (different tissue fixation), these were excluded from the quantification.Staining of neutral lipidsThe staining of neutral lipids was done with BODIPY 493/503 (Cayman Chemical). Flash frozen tissue, cut at 10 µm, was thawed and incubated with BODIPY in 1x PBS for 1 h at room temperature under gentle shaking, protected from light. After incubation, the tissue was washed three times with 1x PBS for 3 min. Sections were counterstained with DAPI, which was added to the first PBS wash. After the final wash, the samples were dried, protected from light, and immediately imaged without mounting.ImagingImaging was done with the Olympus VS120 and Olympus VS200 slide scanners and Leica DM5500 widefield microscope. Images were subsequently analyzed using QuPath82 and FIJI software suites.Pathological assessment of coronary arteriesThe evaluation of atherosclerosis was performed on coronary artery H&E images from GTEx. For each donor, the coronary artery sections were assessed by a clinical vascular pathologist (M.R.) and cross-referenced to the metadata and bulk tissue expression, as described above. For the group of individuals ≤ 49 years of age, 2 samples were excluded based on their medical history (procedures in the area of interest which could affect the results). When assessing the donors of all age groups, no individuals were excluded, as the large sample size prevented individual donors from skewing the data. The prevalence of a condition was calculated cumulatively in function of age (what percentage of individuals below a given age show the condition), with a subsequent linear fit across the data-points for easier visualization. For each age, the odds ratio was calculated with a Chi-squared test.Electron microscopy (EM)Two different preparation methods were used to image VLDL particles in hepatocytes with electron microscopy. A first, optimized for visualizing cell ultrastructure at high resolution with transmission electron microscopy, appears to leave the VLDL particles unstained, therefore light gray. While the second method, optimized for imaging cell ultrastructure in serial images with block face scanning electron microscopy, renders the particles heavily stained, therefore black (Supplementary Fig. 13). This dark staining of these particles is similar to that seen in previous methods where the tissue is stained using osmium tetroxide as the primary fixative54.After pentobarbital euthanasia, and immediately after the heart had stopped beating, approximately 3–5 months old adult mice were perfused with a buffered mix of 1.2% glutaraldehyde and 2.0% paraformaldehyde in 0.1 M phosphate buffer (pH 7.4). The livers were removed 2 h after the perfusion had finished, and were then vibratome sectioned at 100 µm thickness.For transmission electron microscopy, these sections were then washed with cacodylate buffer (0.1 M, pH 7.4), postfixed for 40 min in 1.0% osmium tetroxide with 1.5% potassium ferrocyanide, and then 40 min in 1.0% osmium tetroxide alone. They were then stained for 30 min in 1% uranyl acetate in distilled water before being dehydrated through increasing concentrations of ethanol and then embedded in Durcupan ACM (Fluka, Switzerland) resin. The slices were finally mounted between two glass slides and placed in an oven for 24 hours at 65 °C. Once the resin had hardened, a region of interest from a single slice was cut away from the rest and glued to a blank resin block, with cyanoacrylate glue, and thin (50 nm thick) serial sections were cut with a diamond knife (Diatome, Switzerland). These were collected onto pioloform support films on single slot copper grids and contrasted with lead citrate and uranyl acetate. Images were taken with a transmission electron microscope at 80 kV (Tecnai Spirit, FEI Company with Eagle CCD camera).For serial block face scanning electron microscopy, the vibratome slices of fixed liver tissue were postfixed in potassium ferrocyanide (1.5%) and osmium (2%), then stained with thiocarbohydrazide (1%) followed by osmium tetroxide (2%) alone. They were finally stained overnight in uranyl acetate (1%), washed in distilled water at 50 °C before being stained with lead aspartate at pH 5.5 at the same temperature. The slices were then dehydrated in increasing concentrations of alcohol and embedded in Durcupan resin and hardened at 65 °C for 24 h. Once hardened, the region of interest was cut away from the rest of the slice with a scalpel blade. This piece was glued to an aluminum holder with conductive cement and then placed inside the scanning electron microscope (Merlin, Zeiss NTS) integrated with an in-chamber ultramicrotome device (3View, Gatan). The ultramicrotome was set to remove 50 nm thick slices from the surface of the block and images of the region were taken after each removal. An imaging voltage of 1.7 kV was used and a pixel size of 7 nm. Backscattered electrons were collected with a Gatan backscattered electron detector with a pixel dwell time of 1 µs.Tracer studies129S6/SvEv (Taconic) mice, aged 4-6 months, were placed under an IR lamp 5 min before tail vein injections to ensure vasodilation and visualization of the vein. The mice were then placed into a restrainer, and injected with 200 µl of the buffered solution of VLDL particles containing the DiI (1,1’-dioctadecyl-3,3,3’3’-tetramethylindocarbocyanine perchlorate) fluorophore (Kalen Biomedical). After 15 min, the mice were injected with 150 mg per kg body mass pentobarbital, and dissected 5 min later. At that time-point, all reflexes were absent in all of the mice, and death was confirmed with decapitation. Blood was removed from the liver by perfusion with PBS via the spleen, whereas other imaged organs were cut open and flushed of blood in PBS. Organs were collected, immediately imaged in the IVIS in vivo imaging system (PerkinElmer) in combination with the Living Image 4.5 software.Reporting summaryFurther information on research design is available in the Nature Portfolio Reporting Summary linked to this article.
