Liver samples<\/p>\n
All liver samples were collected with written informed consent from Addenbrooke\u2019s Hospital, Cambridge, UK, according to procedures approved by Local Research Ethics Committees (20\/NI\/0109, 16\/NI\/0196). All participants consented to publication of the research results. The liver samples were snap-frozen in liquid nitrogen and stored at \u221280\u2009\u00b0C in the Human Research Tissue Bank of the Cambridge University Hospitals NHS Foundation Trust.<\/p>\n
Background diseased liver tissue was obtained from individuals with A1AT deficiency or hemochromatosis who were undergoing liver transplantation for hepatocellular carcinoma or liver failure (Supplementary Table 1<\/a>). Patients were identified from clinical history, pre-operative investigations and explant histology (Extended Data Fig. 1<\/a> and Supplementary Table 1<\/a>). All patients with hemochromatosis had current or previously elevated serum ferritin with significant hepatocellular siderosis by Perls\u2019 stain. Pre-operative clinical genotyping was only performed for the common C282Y and H63D variants of the HFE gene. Patients with A1AT deficiency were diagnosed based on PiZZ A1AT electrophoretic phenotype and sub-normal serum A1AT levels. The SERPINA1, HFE, HJV, HAMP, TFR2, SLC40A1, BMP2 and RAB6B genotypes3<\/a> were derived from WGS or WES data (Extended Data Fig. 1<\/a> and Supplementary Table 2<\/a>). Clinical and sequencing data from patients with SLD were derived from previous publications5<\/a>,11<\/a>.<\/p>\n The explant liver histology was reviewed by a specialist liver histopathologist (A.D., blinded to the other study results) and scored according to the Kleiner35<\/a> system on formalin-fixed paraffin-embedded samples away from the fresh\u2013frozen block used for the LCM. The Kleiner score, developed for MASLD (to generate a cumulative MASLD activity score), was used in the absence of a validated scoring system for A1AT deficiency or hemochromatosis, allowing comparability between all study samples regardless of disease etiology. The histological findings were scored for siderosis by Perls\u2019 staining using the Scheuer system36<\/a> and dPAS-positive globule deposition (Extended Data Fig. 1<\/a>) as previously described37<\/a>. Fibrosis was assessed using both the Kleiner35<\/a> and the Ishak38<\/a> scoring systems (Supplementary Table 1<\/a>).<\/p>\n Sample preparation<\/p>\n The protocols for preparing liver tissue sections, LCM and subsequent cell lysis, DNA extraction and WGS have been previously described11<\/a>,39<\/a>. In brief, for six biopsies (PD51605b, PD51606b, PD51607b, PD51608b, PD52285b and PD52286b), 20\u2009\u00b5m-thick tissue sections (prepared with a Leica cryotome) were fixed with 70% ethanol. The other four biopsies (PD60232b, PD60233b, PD60802b and PD60803b) were fixed in PAXgene solutions (PreAnalytiX), processed using a Tissue Tek VIP 6 AI tissue processor (Sakura Finetek) and embedded in paraffin, and 16\u2009\u00b5m-thick sections were generated using an Accu-Cut SRM 200 microtome (Sakura Finetek). All LCM sections were mounted on polyethylene naphthalene membrane glass slides (Leica Microsystems) and stained with H&E for subsequent LCM generation using a Leica Microsystems LMD 7000.<\/p>\n LCM<\/p>\n For six biopsies (PD51605b, PD51606b, PD51607b, PD51608b, PD52285b and PD52286b), 48 microdissections were cut with a target area of 86,000\u2009\u00b5m2, with the same x,y-region cut into the same well from two adjacent z-stacks. For the remaining four biopsies (PD60232b, PD60233b, PD60802b and PD60803b), cuts were taken at 166,000\u2009\u00b5m2 without z-stacking. Overall, microdissection sizes corresponded to 800\u20131,000 hepatocytes. Images were taken before and after LCM. The microdissected samples were then lysed using the Arcturus PicoPure DNA Extraction Kit (Thermo Fisher Scientific) following the manufacturer\u2019s instructions. DNA libraries for Illumina sequencing were prepared using a protocol optimized for low input amounts of DNA39<\/a>. The resulting libraries were submitted for paired-end WGS or WES.<\/p>\n Exome sequencing<\/p>\n Exome capture was performed using either SureSelect All Exon v5 bait set (Agilent, S04380110) or Twist Human Core Exome (Twist Biosciences) bait set. Samples were multiplexed and sequenced using 150\u2009bp paired-end reads on an Illumina NovaSeq 6000 (average pool size of 40). Paired-end reads were aligned to human genome assembly GRCh38 using BWA-MEM40<\/a>. Duplicate reads were marked using biobambam41<\/a> and sample contamination estimates were calculated using VerifyBamID42<\/a>. Library complexity and coverage statistics were calculated using Picard (http:\/\/broadinstitute.github.io\/picard<\/a>). The mean on-target coverage across all samples and genes was 48\u00d7. Across donors, median coverage ranged from 39\u00d7 (PD52285b) to 65\u00d7 (PD60233b).<\/p>\n WGS<\/p>\n WGS was performed using 150\u2009bp paired-end reads on an Illumina NovaSeq 6000. DNA sequences were aligned to the GRCh38 reference genome using BWA-MEM. Duplicate reads were marked using biobambam41<\/a>. The calculation of library statistics was performed using the CollectWgsMetrics function of Picard Tools, rather than the CollectHsMetrics function used for exome sequencing data. The mean coverage across all samples was 42\u00d7.<\/p>\n Calling of single nucleotide variants and indels from exome and WGS data<\/p>\n Single nucleotide variants (SNVs) were called using the CaVEMan (cancer variants through expectation maximization) algorithm43<\/a> against an in silico unmatched normal: a BAM file generated from the human reference genome (GRCh38). Indel calling was performed using cgpPindel44<\/a>, with filtering strategies the same as for SNVs. In addition to the default CaVEMan filters, putative SNVs were required to have a median alignment score of reads showing the variant allele of at least 0.87 and fewer than half supporting reads being clipped (median number of soft-clipped bases, CLPM\u2009=\u20090). Duplicate reads and LCM library preparation-specific artifactual variants resulting from the incorrect processing of secondary cruciform DNA structures were removed with bespoke post-processing filtering39<\/a>. Across all samples from the same patient, we force-called the SNVs and indels that were called in any sample, using a cutoff for read mapping quality (30) and base quality (25). To filter shared artifacts, we applied a beta-binomial-based filtering approach as previously described45<\/a>. Mutations with a dispersion estimate of \u22650.1 were considered true somatic variants. Putative germline variants were filtered using a one-sided exact binomial test used on the aggregated counts of reads supporting the variant and the total depth at that site, as described previously45<\/a>. This tests whether the observed variant counts are likely to have come from a germline distribution or from a distribution with a lower true VAF (probably somatic). For sex chromosomes in male individuals, the binomial probability (true VAF) for comparison was set to 0.95 rather than 0.5. The resulting P\u2009values were corrected for multiple testing using a Benjamini\u2013Hochberg correction. Any variant with a q-value of \u226410\u22125 was classified as a putative somatic variant.<\/p>\n SV calling from WGS data<\/p>\n SVs were called with GRIDDS46<\/a> (v.2.9.4), used with default settings. SVs larger than 1\u2009kb in size with QUAL\u2009\u2265\u2009250 were included. For SVs smaller than 30\u2009kb, only SVs with QUAL\u2009\u2265\u2009300 were included. Furthermore, SVs that had assemblies from both sides of the breakpoint were only considered if they were supported by at least four discordant and two split reads. SVs with imprecise breakends (that is, the distance between the start and end positions of >10\u2009bp) were filtered out. We further filtered out SVs for which the standard deviation of the alignment positions at either ends of the discordant read pairs was smaller than five. To remove potential germline SVs and artifacts, we generated the panel of normals by adding in-house normal samples (n\u2009=\u2009350) to the GRIDSS panel of normals. SVs found in at least three different samples in the panel of normals or in matched normal tissues were removed.<\/p>\n Copy number variant calling from WGS data<\/p>\n Somatic copy number variants were called using the Allele-Specific Copy number Analysis of Tumors (ASCAT) algorithm47<\/a> as part of the ascatNGS package48<\/a> (https:\/\/github.com\/Crick-CancerGenomics\/ascat<\/a>). ASCAT was run with default parameters with the exception of a segmentation penalty of 100. A bespoke filtering algorithm, ascatPCA, was used to reduce the number of false-positive calls that can arise when analyzing genome sequences from normal tissue (https:\/\/github.com\/hj6-sanger\/ascatPCA<\/a>). ascatPCA extracts a noise profile by aggregating the logR ratio from across a panel of normal unrelated samples and subtracts this signature from that observed in the sample being analyzed using principal component analysis.<\/p>\n Extraction of mutational signatures<\/p>\n As previously described5<\/a>,11<\/a>, the hierarchical Dirichlet process (HDP) algorithm, as implemented in the HDP R package (https:\/\/github.com\/nicolaroberts\/hdp<\/a>), was used to extract mutational signatures based on a reference catalog of 65 previously identified 192-context-based mutational signatures from the Pan-Cancer Analysis of Whole Genomes (PCAWG) study49<\/a> as well as a novel liver signature identified previously in CLD5<\/a>.<\/p>\n Bayesian Dirichlet process for clustering VAFs across multiple samples<\/p>\n The nonparametric Bayesian HDP was used to group SNVs based on similar VAFs identified across multiple microdissections in each liver sample. This method, called N-dimensional Dirichlet process (NDP) clustering, has been detailed previously11<\/a>. We ran the algorithm with 15,000 burn-in iterations, followed by 25,000 iterations of posterior Gibbs sampling for the clustering process. During each iteration, there is a defined probability that mutations will be allocated to new clusters that did not exist in the previous iteration. Existing clusters can also be eliminated if member mutations are reassigned to another cluster. This adaptive process allows for the dynamic adjustment of the number of clusters throughout the sampling. To prevent the formation of uninformative clusters, we capped the number of SNV clusters at 100 per patient. We also used a multi-threaded version of the ECR algorithm, adapted from the label.switching R package (https:\/\/doi.org\/10.18637\/jss.v069.c01<\/a>), to correct for label switching efficiently. For subsequent reconstruction of phylogenetic trees, we only considered SNV clusters that included at least 50 distinct mutations.<\/p>\n Construction of phylogenetic trees<\/p>\n The statistical pigeonhole principle50<\/a> was used to deduce the phylogenetic clonal relationships among SNV clusters identified in each patient by the NDP algorithm. Here, each cluster forms a branch on a phylogenetic tree. Evidence that a cluster is considered nested within another is considered strong if its mutation-carrying cell fraction is consistently lower than that of another cluster across all sampled microdissections, and if the combined mutant cell fractions exceed 100%. A combined cell fraction of \u2264100% indicates only weak evidence of such nesting. If only certain microdissections show a lower cell fraction for one SNV cluster compared to another, these clusters are deemed independent, not nested within one another. Our analysis was restricted to clusters with a mutant cell fraction of >0.05. We calculated the cell fraction for each SNV cluster by doubling the median VAF in each microdissection, assuming diploidy. SNV clusters with microdissections lacking shared mutations with others in the same cluster were divided into new, independent clusters. These were then reassessed for their phylogenetic relationships to all other clusters from the same patient biopsy, using the pigeonhole principle. Additionally, a naive Bayes algorithm was used to categorize each identified indel into the SNV clusters detected by the NDP algorithm.<\/p>\n Analysis of driver variants<\/p>\n To determine whether any coding variants were under selection in diseased liver tissues, the dN\/dScv method12<\/a> on the gene level was used. This algorithm identifies genes with an excess of non-synonymous mutations relative to the expected number from the synonymous mutation rate. For this analysis, variants called from whole genomes were collapsed to unique events per SNV cluster identified by the NDP algorithm. Any variants identified in exome data that were not already called in whole genomes from the same patient were collapsed to unique mutations per individual and added to the dN\/dScv input. In addition, we conservatively collapsed SNVs separated by \u226410\u2009bp within the same LCM sample into multinucleotide variants. Genes with q\u2009values of <0.1 were considered to be under positive selection.<\/p>\n Extraction of mutational signatures from SNV contexts using HDP<\/p>\n To identify possibly undiscovered mutational signatures in the liver from A1AT deficiency and hemochromatosis patients, the HDP (https:\/\/github.com\/nicolaroberts\/hdp<\/a>) was run on the 96 trinucleotide counts of all microdissected samples, collapsed to unique mutations across samples. To avoid over-fitting, samples with fewer than 50 mutations were not included in the signature extraction. HDP was run with individual patients as the hierarchy, in 20 independent chains, with a burn-in of 20,000 and the collection of 100 posterior samples off each chain with 200 iterations between each. Owing to the lack of novel signatures in this data set, the remainder of the mutational signature analysis was performed by fitting the identified set of signatures from HDP to trinucleotide counts from each microdissection using the R package sigfit (https:\/\/doi.org\/10.1101\/372896<\/a>).<\/p>\n Protein structure prediction modeling<\/p>\n Structural predictions of the Z-K367*, Z-E387* and Z-M398delFS somatic mutant A1AT proteins were performed using Pymol (v.2.5.2, Schrodinger) to map amino acid sequences onto the crystal structure of native human A1AT (PDB 1QLP<\/a>)15<\/a>, identifying secondary structures lost in somatic mutant proteins. To predict the ability of somatic mutant proteins to form mixed polymers with full-length Z-A1AT, their sequences were mapped onto a single protomer within the crystal structure of a human Z-A1AT polymer (PDB 3T1P<\/a>)16<\/a>. Atoms of the 3T1P protomer that are lost through truncation of somatic mutants were removed from the displayed image to highlight any incompatibility with polymerization. The two missense residues in Z-M398delFS (F396Y, L397D) were modeled using the Pymol mutagenesis function, selecting side chain rotamers of Y396 and D397 that best conserved the resolved density in 3T1P. Somatic mutant structural models did not undergo energy minimization or geometry optimization processes and hence are used only to predict the feasibility of polymerization, not changes to the protein fold. Structural prediction of Z-M398delFS A1AT was carried out using AlphaFold 3 (ref. 51<\/a>), in which the top-ranked model is displayed.<\/p>\n Cloning and constructs<\/p>\n Somatic variants of Z-A1AT (Z-K367*, Z-E387*) were generated by site-directed mutagenesis of the Z-A1AT coding sequence (untagged and N-terminally tagged with HaloTag) used previously17<\/a>. The sequences of primers used in this study are provided in Supplementary Table 3<\/a>. Bicistronic vectors for co-expression of A1AT with an ER marker protein were generated from pGL4.2 encoding A1AT tagged with HaloTag (M-A1AT, Z-A1AT, Z-K367*-A1AT and Z-E387*-A1AT variants) and the ER marker protein moxsynGFP-KDEL by Gibson assembly (New England Biolabs) to form pGL4.2_HaloTag-A1AT_ER-moxGFP-KDEL variants. For the expression of two A1AT forms to mimic heterozygosity, the ER-moxGFP-KDEL in pGL4.2_HaloTag-Z-A1AT_ER-moxGFP-KDEL (SJM1134) was replaced with untagged A1AT variant coding sequences using Gibson assembly. Sequencing of the plasmids was carried out at Plasmidsaurus (Oxford Nanopore Technologies).<\/p>\n Mammalian cell culture<\/p>\n CHO cells (Clontech) were cultured in F12 Ham\u2019s nutrient mixture (Merck) supplemented with 10% FBS and 2\u2009mM GlutaMAX (Thermo Fisher Scientific) at 37\u2009\u00b0C and 5% CO2. CHO Tet-ON cells with inducible expression of Z-A1AT were cultured in high-glucose DMEM (Merck) supplemented with 10% tetracycline-free FBS, 2\u2009mM GlutaMAX and non-essential amino acids (Thermo Fisher Scientific). COS7 cells (MERCK) were cultured in high-glucose DMEM supplemented with 10% FBS and 2\u2009mM GlutaMAX.<\/p>\n Live-cell imaging<\/p>\n CHO parental cells were transfected with a bicistronic vector encoding moxGFP-KDEL and HaloTag-A1AT (M, Z, Z-K367* and Z-E387* variants), performed according to the Lipofectamine LTX reagent protocol (Life Technology). Then, 24\u2009h later, cells were sorted for a GFP-positive population using the BD Influx Cell Sorter and then seeded at 1\u2009\u00d7\u2009104 cells per cm2 in an eight-well ibidi slide. Imaging was performed 48\u2009h after seeding. Before imaging, cells were labeled with 0.5\u2009\u00b5M JFX549 (Janelia Fluor) in Opti-MEM for 15\u2009min at 37\u2009\u00b0C. Airyscan images were collected on an LSM880 (Zeiss) confocal microscope using an Airyscan detector and processed using the Zen (v.2.6) software package (Black edition, Zeiss). Quantitation of ER morphology using co-expression of ER-moxGFP-KDEL was performed using standard confocal imaging using the ImageJ (Fiji, v.2.14.0\/1.54f) software package. Cells were counted based on ER morphology and categorized as either reticular ER or as containing ER inclusions. Quantitation was performed using the moxGFP-KDEL channel to avoid bias introduced by variable accumulation (and hence fluorescence intensity) of different A1AT variants. To quantify ER morphology in cells co-expressing variants of untagged antitrypsin with HaloTagged Z-A1AT, as a mimic of heterozygous expression, CHO cells were transfected as described above with the relevant bicistronic vector. Cells were labeled with JFX549 HaloTag ligand after 48\u2009h expression, then imaged using a Zeiss LSM980 in Airyscan mode. Images were processed using Zen blue software before a gamma correction of 0.5 was applied to allow simultaneous visualization of high-intensity and low-intensity features. Cells were then categorized as having either reticular or inclusion-containing ER morphology.<\/p>\n Single-particle tracking<\/p>\n Single-particle tracking of HaloTag-A1AT variants was performed as described previously18<\/a> using Trackmate (v.7) software22<\/a>. In brief, COS7 cells were transfected 4\u2009h after seeding using 1\u2009\u00b5g DNA and 3\u2009\u00b5l Fugene 6 (Promega) as per the manufacturer\u2019s instructions, to drive expression of an ER marker protein (mEmerald-KDEL) and HaloTag-A1AT variants. Cells were imaged 18\u2009h after transfection, following labeling with PA-JF646 photoactivatable ligand52<\/a>.<\/p>\n Sandwich ELISA for A1AT<\/p>\n CHO cells were seeded at a density of 2.08\u2009\u00d7\u2009104 cells per cm2 in six-well cell culture plates (Greiner Bio-One), as described above for CHO parental cells. Cell culture media was substituted with 1\u2009ml of Opti-MEM reduced serum medium (ThermoFisher Scientific) 24\u2009h after transfection. Conditioned Opti-MEM and cell lysates were collected 48\u2009h after transfection and subject to ELISA for A1AT as described previously17<\/a>, using the monoclonal anti-A1AT antibody mAb3C11 (ref. 53<\/a>) (a kind gift from David Lomas).<\/p>\n Native PAGE and SDS\u2013PAGE<\/p>\n