The Genetic Architecture of Carbon Tetrachloride-Induced Liver Fibrosis in Mice.

BACKGROUND & AIMS: Liver fibrosis is a multifactorial trait that develops in response to chronic liver injury. Our aim was to characterize the genetic architecture of carbon tetrachloride (CCl4)-induced liver fibrosis using the Hybrid Mouse Diversity Panel, a panel of more than 100 genetically distinct mouse strains optimized for genome-wide association studies and systems genetics.
METHODS: Chronic liver injury was induced by CCl4 injections twice weekly for 6 weeks. Four hundred thirty-seven mice received CCl4 and 256 received vehicle, after which animals were euthanized for liver histology and gene expression. Using automated digital image analysis, we quantified fibrosis as the collagen proportionate area of the whole section, excluding normal collagen.
RESULTS: We discovered broad variation in fibrosis among the Hybrid Mouse Diversity Panel strains, demonstrating a significant genetic influence. Genome-wide association analyses revealed significant and suggestive loci underlying susceptibility to fibrosis, some of which overlapped with loci identified in mouse crosses and human population studies. Liver global gene expression was assessed by RNA sequencing across the strains, and candidate genes were identified using differential expression and expression quantitative trait locus analyses. Gene set enrichment analyses identified the underlying pathways, of which stellate cell involvement was prominent, and coexpression network modeling identified modules associated with fibrosis.
CONCLUSIONS: Our results provide a rich resource for the design of experiments to understand mechanisms underlying fibrosis and for rational strain selection when testing antifibrotic drugs.

Genes and pathways underlying susceptibility to liver fibrosis following CCl4 treatment were identified using a systems genetics analysis of a diverse panel of inbred strains of mice. These data provide a rich resource for further mechanistic studies of susceptibility to fibrosis and for testing of antifibrotic drugs.

Liver fibrosis is the wound healing response to various types of chronic injury, including nonalcoholic fatty liver disease (NAFLD) and steatohepatitis (NASH), viral hepatitis (B and C), alcohol abuse, hemochromatosis, and rarer conditions including autoimmune hepatitis and primary biliary cholangitis. Fibrosis is characterized by excess deposition of fibrillar collagen within the hepatic parenchyma. The process is mediated by hepatic stellate cells (HSCs), resident non-parenchymal cells comprising ∼5% of the total liver mass. Liver fibrosis is a multifactorial trait with numerous genetic and environmental factors contributing to disease progression., The heritability of hepatic fibrosis, as assessed by magnetic resonance elastography liver stiffness in a twin study, was estimated to be about 50%, and several other studies also suggest that fibrosis is strongly influenced by genetics.

Because liver biopsy is the gold standard for diagnosing fibrosis and obtaining it percutaneously poses significant risks, only a few genome-wide association studies (GWAS) have analyzed histologically verified fibrosis in humans.7, 8, 9, 10 Huang et al performed genetic association in 420 hepatitis C virus (HCV) patients and developed a risk score that was based on 7 single nucleotide polymorphisms (SNPs), which together with clinical factors predicts cirrhosis development better than clinical factors alone. Chalasani et al performed genetic association in 236 women with NAFLD/NASH and identified a chromosome 7 SNP associated with the degree of fibrosis. Patin et al analyzed SNP associations in 2342 HCV patients and identified 3 loci containing apoptosis-related genes. Patatin-like phospholipase domain-containing protein 3 (PNPLA3), originally discovered to predispose to liver steatosis and hypertriglyceridemia, was later shown to be associated with fibrosis stage.

Studying mice or other animal models can overcome some of the challenges inherent in human studies, including heterogeneous environments and lack of access to tissues. In the present study we used the Hybrid Mouse Diversity Panel (HMDP), a resource consisting of more than 100 inbred strains of mice developed for systems genetics analysis of complex traits., The panel enables high resolution association mapping and genetic replication and has been successfully used in identifying genes and pathways underlying a variety of multifactorial traits.14, 15, 16

The aim of this study was to characterize the genetic factors influencing liver fibrosis in HMDP strains subjected to chemically induced liver damage. To provide a systematically reproducible quantitative measure of liver fibrosis17, 18, 19, 20, 21 for our GWAS, we used an expert digital image analysis system to design an algorithm that fully automates the determination of pathologic collagen in mouse livers. We analyzed genomic association in 693 mice representing 98 HMDP strains that were treated with either carbon tetrachloride (CCl4) or vehicle for 6 weeks. We observed broad variation among the HMDP inbred strains in response to CCl4 treatment. Using a systems genetics approach, we have characterized the genetic architecture of the response to CCl4 and defined high probability candidate genes and molecular pathways underlying the response.

Results

Liver Fibrosis in Response to Liver Injury Is a Genetically Driven Trait

A total of 693 male mice representing 98 HMDP strains were successfully treated with CCl4 (n = 437) or vehicle (n = 256) for 6 weeks, and their livers were harvested for histology and gene expression. For each strain, between 1 and 9 animals were chronically treated with CCl4, and 1–7 animals were treated with vehicle only (average per strain of 4.3 mice and 2.5 mice, respectively).

First, we validated our measurement method by comparing our automated image analysis fibrosis scores with traditional manual scoring for a subset of CCl4-treated mice (93 images from 37 strains). Fibrosis quantification with the automated method strongly correlated with histopathologic scoring of fibrosis performed by a liver pathologist blinded to the strains and collagen proportionate area (CPA%) (Spearman’s r = 0.87, P < 1e-4). We previously validated this method in a mouse model of NASH, showing that CPA% strongly correlated with both hydroxyproline content and a pathologist’s fibrosis score.

Mice developed varying degrees of fibrosis after CCl4 treatment, with some strains exhibiting little or no increase in picrosirius red staining after treatment and others exhibiting dramatically increased staining (Figure 1). The spectrum of fibrosis across the HMDP strains was broad (Figure 2A, Supplementary Table 1), with approximately 12-fold difference in fibrosis between the most resistant and most susceptible strains, ranging from mean of 0.28% in BXH2/TyJ to 3.33% in DBA/2J. In vehicle-treated mice, fibrosis ranged from mean of 0.045% in BXD13/TyJ to 0.425% in MA/MyJ. Notably, there was no correlation between fibrosis in the CCl4- and vehicle-treated mice (Spearman’s r = 0.11, P = .256). The variation between the strains was larger (standard deviation [SD] = 0.72 in CCl4-treated mice) than the variation among mice within strains (average SD = 0.40), consistent with genetic control. Additive effects of the DNA variants explained 31% of the fibrosis variation among the strains (narrow sense heritability), whereas broad sense heritability was 44%. The increase in broad sense as compared with narrow sense heritability suggests the importance of gene-by-gene interactions.

Figure 1

Liver fibrosis histology in a subset of resistant and susceptible mouse strains. Representative picrosirius red stained liver sections from 5 resistant and 5 susceptible strains. The CPA% value from the entire tissue area, automatically determined using Definiens Tissue Studio software, is indicated. Images are ×4 magnification snapshots from whole-slide images scanned with the Aperio ScanScope AT.

Figure 2

Variation in liver fibrosis and survival across mouse strains after CCl-induced liver injury. (A) Mean CCl4-treated (blue, N = 1–9 mice per strain) and vehicle-treated (red, N = 1–7 mice per strain) fibrosis area, as measured by CPA%, plotted for 98 Hybrid Mouse Diversity Panel strains. Error bars denote standard error of the mean. (B) Survival curves of AXB15/PgnJ (N = 17) and NOD/ShiLtJ (N = 16) mice with vehicle, half-dose (0.25 μL/g), or full-dose (0.5 μL/g) CCl4 treatment. Survival after CCl4 injections was dramatically worse in full-dose treated AXB15/PgnJ and NOD/ShiLtJ than for other strains.

Among the 100 strains, most tolerated chronic CCl4 well, with low overall mortality rates. Two strains (AXB15/PgnJ and NOD/ShiLtJ) were highly susceptible to CCl4 injury and did not survive the full course for 6 weeks (Figure 2B). This effect was dose-dependent because these mice, when treated with a half dose of CCl4, had improved survival (Figure 2B). Because we could not reliably complete 6 weeks of CCl4 injections in these strains, we excluded them from the genetic analysis.

Genome-Wide Association Analyses

Genome-wide association was performed for both vehicle-treated and CCl4-treated CPA%. Before GWAS, data were log-transformed to achieve approximately normal data distribution (Figure 3). Thresholds for significant (P < 4.1e-6) and suggestive (P < 4.1e-5) loci were defined using simulation.,

Figure 3

Distribution of CCl- and vehicle-treated fibrosis data. Histograms showing distribution of fibrosis data used for GWAS. (A) CCl4-treated CPA%, normalized to vehicle-treated CPA%, and then log-transformed. (B) Log-transformed vehicle-treated CPA%.

In the CCl4-treated GWAS (Figure 4A), a genome-wide significant locus on chromosome 13 between ∼17.7 Mb and 20.7 Mb was tagged by peak SNP rs36600914 (P = 3.98e-6) at 19,558,275 base pairs (bp) (Figure 4B). This locus overlaps a locus that was associated with hepatic hydroxyproline levels in a female mouse BALB/cJ and FVB/NJ intercross study of CCl4-induced fibrosis. Table 1 lists the genes residing in the chromosome 13 locus and indicates those with a significant local expression quantitative trait locus (eQTL). In a separate locus, a lone SNP rs13481805 at 47,918,469 bp on chromosome 13, located in an intron of non-coding lincRNA G630093K05Rik and upstream of gene Id4, reached genome-wise significance (P = 3.39e-6). This SNP is located near the peak of a locus spanning 44.2–52.7 Mb that was previously associated with CCl4-induced liver fibrosis in a BXD study. There is also a significant lone SNP in an intron of Maml3 on chromosome 3, rs45898322 (P = 1.71e-6), that has a minor allele frequency <5% in our GWAS. Maml3 was identified as a gene associated with liver fibrosis due to HCV infection in a European human GWAS study primary cohort.

Figure 4

GWAS identifies loci associated with CCl-induced liver fibrosis. (A) Manhattan plot showing GWAS results of CCl4-induced fibrosis normalized to vehicle. Loci with SNPs above the red bar (-log10 (P) > 5.4) are statistically significant, whereas loci above -log10 (P) > 4.4 are suggestive. (B) LocusZoom plot of chromosome 13: 16–21 Mb significant locus. Colors of SNPs indicate r2 measure of linkage disequilibrium with the labeled peak SNP.

Table 1

Genes Expressed in Livers of CCl4-Treated Mice in Chromosome 13 Significant GWAS Fibrosis Locus

Gene symbol	Gene name	Best cis-eQTL rsID (P value)a	eQTL for rs36600914 (P value)b	Correlation of cis component of expression with fibrosis: bicor (P value)c	Correlation with fibrosis: bicor (P value)d	Deleterious mutation types?e
Inhba	inhibin beta-A	—	.612312	—	0.040 (.69)	no
Gm48493	predicted gene, 48493	—	.128485	—	0.009 (.93)	no
Gm48496	predicted gene, 48496	—	.173772	—	0.078 (.44)	no
Sugct	succinyl-CoA glutarate-CoA transferase	—	.0876457	—	–0.115 (.26)	no
Gm48499	predicted gene, 48499	—	.12715	—	0.044 (.67)	no
Gm48501	predicted gene, 48501	—	.319134	—	0.052 (.61)	no
Mplkip	M-phase specific PLK1 interacting protein	rs29872460 (4.1e-08)	.00628365	–0.104 (.31)	–0.050 (.63)	no
Cdk13	cyclin-dependent kinase 13	—	.0100412	—	0.305 (.0022)	yes
Rala	v-ral simian leukemia viral oncogene A (ras related)	—	.0295267	—	0.037 (.72)	no
Gm5628	predicted gene 5628	—	.0158592	—	–0.015 (.88)	no
Yae1d1	Yae1 domain containing 1	rs47401911 (8.3e-07)	.133346	.178 (.08)	0.194 (.056)	yes
Gm48800	predicted gene, 48800	—	.00278421	—	0.035 (.73)	no
Vps41	VPS41 HOPS complex subunit	rs30064860 (5.9e-08)	.00700974	.010 (.92)	0.028 (.79)	no
Gm47659	predicted gene, 47659	—	.164951	—	–0.068 (.5)	no
Tcrg-C2	T-cell receptor gamma, constant 2	—	.577282	—	–0.022 (.83)	yes
Tcrg-C4	T cell receptor gamma, constant 4	—	.111916	—	0.211 (.037)	no
Stard3nl	STARD3 N-terminal like	rs29631766 (6.2e-06)	.000105595	–0.251 (.013)	0.015 (.88)	yes
Epdr1	ependymin related protein 1 (zebrafish)	rs37221784 (4.5e-13)	1.48205E-05	–0.009 (.93)	0.035 (0.73)	yes
Gpr141	G protein-coupled receptor 141	—	.415455	—	0.080 (.43)	yes
Gm19154	predicted gene, 19154	—	.741509	—	0.229 (.023)	no
Elmo1	engulfment and cell motility 1	—	.253941	—	0.170 (.094)	yes
Aoah	acyloxyacyl hydrolase	rs29599428 (1.5e-10)	.000333089	0.024 (.82)	0.030 (.77)	yes

NOTE. Parameters to prioritize genes for further study include significant: acis-eQTL; beQTL at peak GWAS SNP; ccorrelation between cis component of expression and CCl4-induced fibrosis; dcorrelation of overall expression with CCl4-induced fibrosis, and epresence of a potentially deleterious mutation. Here, missense, nonsense, and splicing site variants reported in 1 or more strains (129S1/SvlmJ, A/J, AKR/J, BALB/cJ, C3H/HeJ, CBA/J, DBA/2J, FVB/NJ, LP/J, C57BL/6J) at https://www.sanger.ac.uk/sanger/Mouse_SnpViewer/. Significant cis-eQTL (P < 1e-4); bicor (P < .05).

Additional suggestive loci were identified on other chromosomes (Figure 4A). One of 2 loci on chromosome 3 is at 128–132 Mb, with peak SNP rs37226384 (P = 1.41e-5) near genes Enpep and Pitx2 (Figure 5A). The other chromosome 3 locus at 139–145 Mb has a peak SNP at rs6390575 (P = 3.30e-5) (Figure 5B). On chromosome 4, lone SNP rs27853444 (P = 1.49e-5) is in an intron of Svep1. A broad locus overlapping this SNP was also associated with CCl4-induced liver fibrosis in the BXD study. Finally, a broad locus spanning 55–70 Mb on chromosome 7 falls just below the significance threshold (Figure 6). It has 2 peak SNPs, rs32424606 (P = 2.16e-5 at 56,597,169 bp) nearest to the Oca2 gene and rs31808490 (P = 3.40e-5 at 65,679,712 bp) in an intron of Tarsl2. Although these SNPs are in high linkage disequilibrium, they may tag 2 overlapping loci. This same locus (48.2–74.2 Mb) was also associated with CCl4-induced liver fibrosis in the BXD study.

Figure 5

Suggestive loci on chromosome 3 associated with CCl-induced liver fibrosis. LocusZoom plot of the (A) 128–132 Mb region locus and (B) 139–145 Mb region locus. Colors of SNPs indicate r2 measure of linkage disequilibrium with the labeled peak SNP. SNPs with -log10 (P) > 5.4 are statistically significant and with -log10 (P) > 4.4 are suggestive.

Figure 6

Suggestive broad locus on chromosome 7 associated with CCl-induced liver fibrosis. LocusZoom plot of the 55–70 Mb region locus, showing 2 peak SNPs that may tag 2 overlapping loci. Colors and shading of SNPs indicate r2 measure of linkage disequilibrium with the 2 labeled peak SNPs. SNPs in the locus are labeled with the color of the peak SNP with highest shared r2. SNPs with -log10 (P) > 5.4 are statistically significant and with -log10(P) > 4.4 are suggestive.

Suggestive loci were also found in the vehicle-treated fibrosis GWAS, but none of these loci overlapped with those from the CCl4-treated GWAS (Figure 7A). A signal on chromosome 3 spanning 81.6–85.6 Mb with peak SNP rs30271424 (P = 6.07e-6) is close to a cluster of fibrinogen-encoding genes including Fgg, Fgc, and Fgb (Figure 7B). A suggestive locus on chromosome 16 (peak SNP rs4202102, P = 2.32e-5) is centered just upstream of Robo1 (Figure 8A). A final suggestive locus is on chromosome 4 (peak SNP rs27515792, P = 4.2e-6) at 111–113 Mb (Figure 8B). The top SNPs in this locus are in or near the genes Spata6, Slc5a9, Bend5, and Agbl4.

Figure 7

GWAS identifies suggestive loci associated with baseline liver fibrosis in vehicle-treated mice. (A) Manhattan plot showing GWAS results of fibrosis in vehicle-treated mice. The red bar (-log10 (P) = 5.4) indicates the significance threshold; loci with -log10 (P) > 4.4 are suggestive. (B) LocusZoom plot of the chromosome 3: 81–86 Mb suggestive locus with peak SNP rs30271424 (P = 6.07e-6). Colors of SNPs indicate the r2 measure of linkage disequilibrium with the labeled peak SNP.

Figure 8

Suggestive loci associated with baseline liver fibrosis in vehicle-treated mice. (A) LocusZoom plot of the chromosome 16: 69–75 Mb region locus, with peak SNP rs4202102 (P = 2.32e-5). (B) LocusZoom plot of the chromosome 4: 111–113 Mb locus, with peak SNP rs27515792 (P = 4.2e-6). Colors of SNPs indicate the r2 measure of linkage disequilibrium with the labeled peak SNP. SNPs with -log10 (P) > 5.4 are statistically significant, and with -log10 (P) > 4.4 are suggestive.

Differentially Expressed Genes and Pathways Associated With Liver Fibrosis

Total RNA from the livers of 1 mouse per strain in both the vehicle- and CCl4-treated groups was analyzed by RNAseq. Figure 9A shows a volcano plot for genes differentially expressed in response to CCl4 across all strains. Prominently down-regulated genes were members of the mouse urinary protein (MUP) family encoded by a cluster of genes on chromosome 4. These genes are regulated by testosterone and insulin resistance and have been associated with energy expenditure. MUPs were also recently reported to be some of the most down-regulated genes in response to CCl4 treatment of male C57BL/6N mice. Various solute carrier family members including glutamate transporter 1 (Slc1a2) and organic cation transporter (Slc22a3) were also markedly down-regulated. A number of highly up-regulated genes are suggestive of HSC and macrophage involvement: macrophage metalloproteinase Mmp2, a receptor for the innate immune system (Trem2), type 1 transmembrane glycoprotein (Gpnmb), and a type 1 collagen (Col1a1), the major component of fibrosis produced only by activated HSCs.

Figure 9

The top differentially expressed genes between CCl- and vehicle-treated mice across all 98 strains. (A) Volcano plot of statistical significance (-log10Q value) versus log2 (fold change) for differentially expressed genes in CCl4- versus vehicle-treated mice. Red (lower expression versus vehicle) and green (higher expression versus vehicle) genes have absolute log2 (fold change) greater than 1 and Q values < .05. (B and C) Pathway enrichment analysis annotations and their associated -log10 (P) values for the top 300 genes (ranked by Q value, excluding genes with absolute fold change < 1.5) with (B) higher and (C) lower expression in the CCl4-treated mice.

For a global view of pathways affected by CCl4 across all strains, we performed gene-enrichment analysis on the top differentially expressed genes between vehicle- and CCl4-treated mice using Metascape. The top 300 genes with increased expression (fold increase > 1.5) in response to CCl4 treatment, as ranked by Q value (Supplementary Table 2A), were enriched in pathways involving extracellular matrix (ECM), wound healing, and immune response (Figure 9B). These are pathways most associated with Kupffer cells and activated HSCs. In addition, 48 of these 300 genes were identified as top HSC state-defining genes during a CCl4 treatment time course study in C57BL/6J mice.

The top 300 genes with decreased expression (fold decrease > 1.5) in response to CCl4 treatment, as ranked by Q value (Supplementary Table 2B), were enriched in pathways including biological oxidations (involving mostly cytochrome P450 genes and including Cyp2e1, the primary gene that metabolizes CCl4 to a toxic free radical), steroid and lipid metabolism, cellular response to xenobiotic stimulus, drug metabolism, and solute-coupled transport, pathways most associated with hepatocytes (Figure 9C).

We next identified genes whose expression levels across the HMDP strains correlated with fibrosis. The top 500 genes whose expression was positively correlated with fibrosis in the CCl4 dataset (P < .05, Supplementary Table 3A) showed significant enrichment in several pathways including ECM organization and proteoglycans, blood vessel development, ossification, skeletal system development, elastic fiber, and collagen formation (Figure 10A). These pathways have been previously associated with liver injury and confirm involvement of HSCs. In a cell-type specific proteomics study of the liver, pathway terms including ossification, collagen, and ECM were uniquely found to be associated with HSCs and not hepatocytes, Kupffer cells, or liver sinusoidal endothelial cells. There was also significant enrichment in several pathways for genes whose expression was negatively correlated with fibrosis (P < .05, Supplementary Table 3B). These pathway annotations include metabolism of lipids, monocarboxylic acid and isoprenoid metabolic process, small molecule catabolic process, and peroxisome (Figure 10B). These terms were uniquely associated with hepatocytes in the cell-type specific liver proteomics study and thus may represent a loss of hepatocyte function in the injured livers.

Figure 10

Pathways associated with genes whose expression in CCl-treated mice is most correlated with CCl-induced fibrosis. Top pathway enrichment analysis annotations and their associated -log10 (P) values for genes in CCl4-treated mouse livers with expression significantly (bicor P <.05) (A) positively (top 500 genes) and (B) negatively correlated with fibrosis.

Finally, we looked at the genes in the vehicle-treated group that were either positively (Supplementary Table 4A) or negatively (Supplementary Table 4B) correlated with fibrosis in the CCl4-treated group, which could give clues into initial (pretreated) states that predispose or protect from fibrosis. Top pathways positively associated with CCl4-induced fibrosis included transfer RNA modification and glycan degradation (Figure 11A), whereas pathways negatively associated involved RNA and DNA metabolism and regulation of the cellular response to stress (Figure 11B). Modification of transfer RNAs is known to influence cellular responses to toxins by affecting translation fidelity and overall translation levels.

Figure 11

Pathways associated with genes whose expression in vehicle-treated mice is most correlated with CCl-induced fibrosis. Top pathway enrichment analysis annotations and their associated -log10 (P) values for genes in vehicle-treated mouse livers with expression significantly (bicor P <.05) (A) positively and (B) negatively correlated with fibrosis in CCl4-treated mice.

Genetics of Hepatic Gene Expression and Fibrosis

To prioritize genes at GWAS loci and to test for genes whose regulation is associated with fibrosis, we performed eQTL analysis. The eQTLs were classified as local (within 1 Mb of the gene) and likely cis-regulated (cis-eQTL) or trans-regulated (trans-eQTL). We also identified genes throughout the genome whose local, or cis, component of expression variance correlates with the degree of CCl4-induced fibrosis (Supplementary Table 5). Such genes are particularly strong causal candidates because the association cannot be due to the effect of the trait on gene expression (that is, a reactive relationship). Table 1 lists genes expressed in the CCl4-treated mice that are located in the chromosome 13 fibrosis GWAS locus, with significant cis-eQTL (P < 1e-4) noted. In this locus, 6 genes had significant cis-eQTL, and the cis component of the variance in expression of 1 of these genes, Stard3nl, significantly correlated with fibrosis, making this a potential candidate for further study. Supplementary Tables 6A–D provide details of expressed genes, eQTLs, and their correlations with fibrosis for the significant and suggestive CCl4-induced fibrosis GWAS loci.

To visualize liver eQTL across the genome, the position of each associated SNP was plotted versus the position of the regulated gene (Figure 12A and B). The dots on the diagonal correspond to the positions of locally regulated genes, and the dots off the diagonal correspond to trans regulated genes. Overall, the most significant eQTLs were robust and present in both the vehicle- and CCl4-treated mice. Linear vertical patterns of eQTL, denoting SNPs that influence the expression of many genes in trans, were also visible. On chromosome 1, two such “eQTL hotspots” in the vehicle-treated samples disappear in CCl4-treated mice, whereas a hotspot in the middle of chromosome 14 is present in both the vehicle- and CCl4-treated mice. Additional micro hotspots that are affected by treatment not readily visible here can be identified in the bulk eQTL data.

Figure 12

Expression quantitative trait locus (eQTL) hotspots regulating liver gene expression in vehicle- (- ( Cis-eQTL (within 1 Mb of gene) with P < .0021, and trans-eQTL with P < 1e-9, plotted by gene position against SNP position. The -log10 (P) values of eQTL are colored according to the legend at right. Cis-eQTLs are clustered on the diagonal. The visual appearance of vertical lines of eQTLs represents eQTL hotspots.

Network Modeling of Fibrosis

We analyzed the expression data using weighted gene coexpression network analysis (WGCNA) to model coexpression networks in liver and to understand the association of gene networks with fibrosis. WGCNA is a global analysis aimed at identifying, in an unbiased manner, genetic pathways associated with clinical traits and is used to aggregate gene expression into groups of highly coexpressed genes, called modules. The first principal component of each module was then tested for correlation to liver fibrosis to identify gene clusters associated with this trait. We identified 23 coexpressed gene modules in the vehicle-treated mice (ranging in size from 31 to 5677 genes) and 24 gene modules in the CCl4-treated mice (ranging in size from 31 to 6049 genes), inclusive of the grey (unassigned genes) module (Figure 13A, Supplementary Tables 7A and B).

Figure 13

WGCNA identifies 24 coexpression modules correlated with CCl-induced liver fibrosis. (A) List of 24 modules identified by WGCNA, with annotations, correlation with CCl4-induced fibrosis (bicor), correlation P value, and number of genes in each module. Annotations were based on gene descriptions and pathway enrichment terms (Reactome and Metascape analysis). Modules with bicor P < .05 are considered significant. (B) Genes in the green module; edges indicate coexpression correlations with P < 1e-5.

Only 1 module in the vehicle-treated mice, the dark green module enriched for genes encoding MUPs, was significantly (P < .05) correlated with CCl4-induced fibrosis (bicor 0.20; P = .047). Five modules in the CCl4-treated mice were significantly correlated with fibrosis (Figure 13A). The most significantly correlated module, the green module, was positively correlated with fibrosis and enriched in collagens and other ECM genes (Figure 13B). The other 3 significant and positively correlated modules were enriched in interferon signaling (dark green), inflammation and ECM (tan), and signaling (turquoise) genes. The purple module was significantly negatively correlated with fibrosis and was enriched in genes involved in cholesterol and fatty acid metabolism.

To identify genetic loci that influence module expression, QTL mapping was performed on the first principal components of the detected modules in the CCl4-treated mice. Significant loci (5% false discovery rate, P < 3.95e-5) were identified for several modules (dark green: rs31012842, P = 2.33e-6; red: rs51233520, P = 5.20e-7; dark red: rs26921889, P = 1.10e-6; royal blue: rs49350632, P = 1.0e-7; cyan, dark turquoise, grey60, light yellow, magenta, midnight blue: multiple loci), which could be explored further in future studies (Supplementary Tables 8A–X).

Discussion

We report an integrative genetics analysis of liver fibrosis in response to CCl4 across 100 diverse inbred mouse strains of mice. We observed about 12-fold range of variation in response and a significant genetic component. Loci contributing to the response were identified and integrated with global transcriptomic data. The global expression data were also used to identify pathways associated with the degree of fibrosis. Our results provide a rich resource for future experimental studies of hepatic fibrosis. Below, we discuss these points in turn.

We validated and then performed CPA% analyses using a novel fully automated digital pathology system that excludes vascular and capsular collagen, thereby quantifying only pathologic collagen from the entire tissue section. A clear advantage of this automated image analysis approach is that no post hoc analysis by a human pathologist is needed to correct for possible overestimation of liver fibrosis. The CPA% among the strains exhibited continuous variation from less than 0.3% to about 3.3% (Figure 2A). Mice treated with the vehicle also exhibited variation in fibrosis levels, averaging about 0.3%. There was no relationship between fibrosis in the vehicle- versus the CCl4-treated mice. On the basis of heritability estimates, about 44% of the variance in fibrosis levels was explained by genetic differences between strains.

Our study confirms the diversity in the fibrotic response and highlights that multiple mouse strains should be considered when studying fibrosis. The most commonly used laboratory strain, C57BL/6J, lies more on the fibrosis resistant side of the spectrum when challenged with CCl4 (Figure 2A). This raises the question of how many antifibrotic agents, trialed only in C57BL/6J, failed to meet endpoints because that strain is not particularly susceptible to fibrosis. CRISPR/Cas9 now makes it feasible to work with any or multiple genetic backgrounds, and data from this study can be used to guide strain selection. For example, to study the effect of knockdown or overexpression of a candidate fibrosis gene, an investigator could choose strains on the high or low end of basal or delta expression, respectively, to maximize effect. In addition, an investigator could test whether an antifibrotic is effective across a variety of strains on the spectrum of fibrosis; broad actors may be more likely to play an important role in humans.

We performed GWAS of fibrosis using a mixed model to correct for population structure. One significant locus (on chromosome 13) and 3 suggestive loci were identified. Previous studies have mapped CCl4-fibrosis in mice to genetic loci, some of which overlap with those in this study. Because of differences in study design, analysis methods, allele frequency, and sex effects, differences in the loci detected in each study are expected. However, overlapping loci are particularly strong candidates for further study. In an intercross F2 study with BALB/cJ and FVB/NJ parental strains, loci on multiple chromosomes were found to be significantly or suggestively linked with hydroxyproline and fibrosis. The chromosome 13 locus from that study (14.5–24.5 Mb) overlapped with the most significant locus in our study (17.7–20.7 Mb). There was no overlap in loci between our study and another intercross F2 study with BALB/cJ and A/J parental strains, likely because of multiple factors mentioned above. Finally, in an association study with 35 strains of BXD mice, 9 significant and many suggestive loci were identified. In our study, which includes 29 BXD strains, lone SNPs on chromosomes 4 and 13 overlap with loci in the previous BXD study. In addition, the broad suggestive chromosome 7 locus in our study (55–70 Mb) overlaps with 48–74 Mb locus in that study.

We prioritized candidate genes at loci on the basis of co-mapping of eQTL as well as correlation with fibrosis response. In the CCl4-treated mice, the cis-component of the variance in expression of Stard3nl at the chromosome 13 locus was significantly correlated with fibrosis, making it a potential gene of interest. In the chromosome 3: 139–145 Mb suggestive locus, Clca3a2 was the only gene whose cis-component of the variance in expression was correlated with CCl4-induced fibrosis, making it a top candidate gene. Clca3a2 encodes a calcium-activated chloride channel accessory subunit that has been reported to protect keratinocytes from apoptosis in response to hyperosmotic stress. Finally, 3 top candidate genes in the chromosome 7 locus (Mphosph10, Chsy1, and Lins1) had significant cis-eQTL in the CCl4-treated mice, expression significantly correlated with fibrosis, and the cis-component of their variance in expression was significantly associated with fibrosis.

Our global gene expression data also provided insight into the genes and pathways associated with CCl4-induced fibrosis across many genetic backgrounds. We identified thousands of transcripts with significantly different expression between vehicle- and CCl4-treated mice as well as genes correlated with CCl4-induced fibrosis. By WGCNA, we identified and mapped gene coexpression modules associated with fibrosis. Overall, these analyses confirmed the involvement of HSCs and loss of hepatocytes as being associated with fibrosis. In addition, they highlighted associations between fibrosis and the ECM and immune response. In vehicle-treated mice, genes involved in DNA and RNA metabolism and the cellular response to stress were most associated with fibrosis in the CCl4-treated mice, with lower expression of these genes correlated with higher fibrosis. These pathways are involved in the initial response to liver damage, including repair and regeneration., Recently, Ghallab et al reported a loss of hepatic pericentral gene expression in response to CCl4 treatment in mice. One reflection of such cell composition changes is “eQTL hotspots” in which a SNP is associated with expression changes in many genes. Several such trans-eQTL hotspots were present in both the vehicle- and CCl4-treated samples that could be further explored.

The chromosome 13 locus is highly conserved in the syntenic human region on chromosome 7 and has been associated in several human GWAS studies with fibrotic conditions. This locus in humans was the most significant signal in 2 GWAS of Dupuytren’s contractures, a benign fibrosing disease of the palmar fascia, commonly found in patients with alcoholic cirrhosis., Although the causative gene in this locus is not yet known, EPDR1 is a top candidate that has been functionally implicated in myofibroblasts., This locus (at ELMO1) was also a human GWAS hit for primary biliary cholangitis. ELMO1 (engulfment and cell motility 1) has been reported to confer susceptibility to human diabetic nephropathy, which is a disease of excess ECM deposition, similar to fibrosis. Elmo1 also induced collagen and fibronectin expression when overexpressed in a monkey kidney cell line. Finally, this syntenic locus (at AOAH) was also associated with staining of alpha-smooth muscle actin, a marker of activated HSCs that is functionally involved in hepatic fibrosis, in a small GWAS of mostly non-fibrotic donor liver tissues. AOAH encodes acyloxyacyl hydrolase, a lipase that inactivates lipopolysaccharide. AOAH was recently shown to influence lipopolysaccharide-induced lung and primary hepatocyte injury in mice. Finally, genes specific for gamma-delta T cells are present at this locus near the peak SNP. These cells were recently shown to contribute to steatohepatitis in mice, and their activation was required for the development of liver fibrosis and inflammation in Mdr2 knockout mice, a model of cholestatic liver disease. Although none of these genes were top candidates in our analysis based on liver gene expression, some have potentially deleterious mutations in the HMDP mouse strains that may affect protein function.

A limitation of this study is that we do not know the relevance of our results in mice to human exposure to liver toxins, other than some apparent overlap of GWAS loci. Also, although our study has identified high confidence candidate genes for fibrosis, follow-up with knockdown/overexpression in specific cell types in the liver will be required for full validation.

In conclusion, we present a resource for the molecular genetic dissection of liver toxicity. All data are available to interested researchers through our website systems.genetics.ucla.edu or by direct request.

Methods

Mice

The HMDP male mice were generated and maintained at UCLA or obtained from the Jackson Laboratory. Mice were fed a standard rodent chow diet ad libitum and housed with woodchip bedding. CCl4 (Sigma-Aldrich, St Louis, MO) was administered intraperitoneally at a dose of 5 μL of 10% solution in olive oil (Sigma-Aldrich) per gram of body weight. CCl4 injections were started at 10–12 weeks of age and given twice weekly for 6 weeks. Olive oil was used for vehicle-injected mice. Mice were euthanized 72 hours after the final injection and after an 8- to 12-hour fast, and their livers were harvested and divided into aliquots for RNA isolation (snap-frozen) and histology (formalin-fixed and paraffin-embedded).

Liver Fibrosis Histology and Image Analysis

Paraffin-embedded livers were sectioned at 4-μm thickness and stained with picrosirius red to specifically highlight fibrillar collagen. Sections were deparaffinized by incubation in xylenes and rehydrated by a graded series of decreasing amounts of ethanol in water. Slides were then stained with Weigert's hematoxylin (Sigma-Aldrich) for 8 minutes and washed for 10 minutes under running tap water. They were stained with 0.1% picrosirius red for 1 hour (Direct Red 80 and 1.3% picric acid from Sigma-Aldrich) and washed in 2 changes of acidified water (10 seconds in the first and the remaining 10 minutes in the second, 5 mL of acetic acid per 1 L of water). Slides were then dehydrated in 3 changes of absolute ethanol, cleared in 2 changes of xylene, and mounted with Permount (Fisher Scientific, Pittsburgh, PA). Slides were scanned at ×20 magnification with Aperio ScanScope AT (Aperio, Vista, CA). Using Definiens Tissue Studio (Definiens AG, Munchen, Germany), we designed an image analysis algorithm to quantify fibrosis as CPA% of the whole section, excluding the normal vascular wall and liver capsular collagen. This was done by excluding 20 μm adjacent to any white space and quantifying red pixels from the remaining tissue area. The algorithm was automatically applied to all liver sections. To validate the algorithm, we had an expert pathologist (S.F.) blindly assess 93 slides from the CCl4-injected mice and grade fibrosis by using a scoring system based on the pattern of injury and accounting for the degree of pericentral, periportal, and bridging fibrosis.

Genome-Wide Association Analysis and Estimation of Heritability

Genotyping of the mouse strains was performed by using the Mouse Diversity Array. SNPs were filtered for 5% minor allele frequency and 10% missingness rate, resulting in 208,514 informative markers. A linear mixed model as implemented in the FaST-LMM program was used to perform the association with fibrosis area. The mouse strains in the HMDP are genetically different and can be treated as “individuals” in the analysis. However, the same individuals appear in the case and control samples, making these 2 samples not independent and requiring pairing of the mouse strains across the case and control samples. Pairing cases and controls by mouse strain is complex because of the difference in the number of cases and controls within each strain, with the number of controls usually being smaller. We selected the following approach to pair a control with each case. Within each strain, for each case, we selected a control randomly and with replacement when necessary. Using this approach, the cases are all preserved, and the variability among the controls is retained. Thus, to normalize CCl4-treated fibrosis values to vehicle-treated values, mice from each group were randomly sampled without replacement to create pairs. In cases where there were unequal numbers of CCl4- or vehicle-treated mice for a given strain, mice from the smaller group were sampled without replacement and assigned to members of the larger group, and then random sampling with replacement from the smaller group was used to assign partners to remaining members of the larger group. The vehicle-treated fibrosis area in each pair was then subtracted from the CCl4-treated fibrosis area, and a single integer value was added to all results to make all values positive. Traits were then log-transformed to normalize the distribution for GWAS.

Log-transformed vehicle-normalized CCl4-treated fibrosis values were also used for all heritability and correlation analyses. Broad sense heritability was calculated by using the R package “heritability”. Narrow sense heritability was calculated by using the Genome-wide Complex Trait Analysis Tool (GCTA).

Global Gene Expression Analysis

Total RNA from mouse liver samples (1 mouse per strain) was purified by using the Qiagen miRNeasy kit (Qiagen cat#217004) per the manufacturer’s instructions. RNA sequencing libraries were prepared by using the Illumina (San Diego, CA) TruSeq Stranded mRNA Sample Preparation protocol and sequenced in paired-end mode to a length of 2 × 76 bp by using the Illumina HiSeq2500 platform. Reads were quantified against the GRCm38.p6 mouse reference transcriptome (Ensembl release 97) using kallisto version 0.46.0 with 100 bootstrap replicates. For all analyses except differential expression (where program default low abundance filtering parameters were used), only genes with expression values >0.1 transcripts per million (TPM) in at least 20% of samples and ≥6 reads in at least 20% of samples were included.

Pathway and Network Modeling

Correlations between liver fibrosis and gene expression were analyzed by using the bicor package in R. Pathway analysis of top correlated genes was performed with human pathway settings by using Metascape. Differential expression was performed with the sleuth R package (version 0.30) using a likelihood ratio test comparing a full model accounting for strain and treatment with a reduced model accounting only for strain.

Expression quantitative trait loci were mapped from TPM values by using FaST-LMM after log2-transformation. To calculate the cis-component of expression correlated with CCl4-induced fibrosis, the expression values for each relevant gene were partitioned into groups based on the genotype of the most significant cis-eQTL with P <1e-4. Medians were calculated for each resulting group, and the medians (replicated once for each individual within each group) as a whole were correlated against fibrosis.

Network analysis was performed by using the WGCNA R package. WGCNA was used to group highly coexpressed genes into modules. To generate a coexpression network, an adjacency matrix was created by calculating pairwise gene-gene correlations and then raising the Pearson correlation to the 8th power, selected using the scale-free topology criterion. A topological overlap matrix-based dissimilarity measure was used for hierarchical clustering of the genes. Gene modules corresponded to the branches of the resulting dendrogram and were defined by using the “Dynamic Hybrid” branch cutting algorithm. For module generation, “cut height” was set to 0.995 and “minimum module size” to 30. Using the bicor package in R, the first principal component of each module was analyzed for correlation with fibrosis to identify associated gene clusters and then mapped by using FaST-LMM.

Statistical Analysis

A balanced study design with an appropriate sample size was proposed to achieve adequate statistical power for all analyses. However, the random effects of mouse fertility in our breeding program, as well as some losses, resulted in some differences in sample sizes among the hybrid panel breeds as well as a lack of balance between the cases and controls. The reduced sample sizes and lack of balance can impact the power to detect genetic effects, but there is no anticipation of an increase in the detection of false genetic effects beyond what is expected using our conservative levels of significance.

The nonparametric Mann-Whitney test was used to analyze differences between 2 groups. The Spearman method was used to analyze correlation between CPA% and histopathologic scores and between strain average CPA%. Analyses were performed by using GraphPad Prism (GraphPad Software, La Jolla, CA) and in R. P <.05 was considered significant for these tests and for bicor analyses. P <1e-4 was considered significant for cis-eQTL. For GWAS, the threshold value for genome-wide significance (P < 4.1e-6) was based on simulation. All reported P values are based on a two-sided hypothesis.

Ethics Statement

All animal experiments were approved by UCLA Chancellor's Animal Research Committee of the IACUC, under ARC#2011-142-11. Mice were expected to lose up to 10% of body weight in the first 10 days after initiation of CCl4 (prior literature/our observations). A monitoring protocol was designed to minimize pain/distress to the animals and was based on observation and weights. After CCl4 injections, mice were reexamined within 6–10 hours to euthanize any moribund animals. Mice were weighed every other day during the first 10 days after initiation of the CCl4 protocol and at least 2× per week during the remainder of the 6-week period. Mice that lost 10%–20% body weight were monitored and weighed daily. If they lost >20% body weight, they were euthanized. Slow, shallow, or labored respirations with immobility were also criteria for premature euthanasia, performed by isoflurane overdose followed by cervical dislocation.

In vivo experiments are reported in accordance with the ARRIVE guidelines. All authors had access to the study data and reviewed and approved the final manuscript.