All animals were handled in strict accordance with good animal practice as defined by the relevant national and/or local animal welfare bodies, and all animal work was approved by the appropriate committee. All experiments in this paper were carried out with UCLA IACUC approval.
Animals and clinical phenotype collection
Male mice from the HMDP panel, approximately 6–10 weeks of age, were purchased from Jackson Labs and were fed Purina Chow (Ralston-Purina Co., St. Louise, MO) at 16 weeks of age. All mice were maintained on a 12 h light/dark cycle. At 16 weeks of age, whole body fat, fluids and lean tissue mass of mice were determined using a Bruker Optics Minispec nuclear magnetic resonance (NMR) analyzer (The Woodlands, TX, USA) according to the manufacturer's recommendations. We also calculated the total mass of the mice, sum of lean mass, free fluid, and fat mass, and body fat percentage, fat mass/total mass. Following a 16-hour fast, mice were weighed and then bled retro-orbitally under isoflurane anaesthesia. Complete blood counts were performed using a Heska CBC-Diff analyzer (Heska Corp, Loveland, CO, USA). An external control sample with known analyte concentration was run in each plate to ensure accuracy. Glucose levels were determined using commercially available kits from Sigma (St Louis, MO, USA). Insulin levels were measured using commercial ELISA kits (ALPCO Diagnostics). All measurements were performed in triplicate according to the manufacturer's instructions. Plasma lipids were determined as previously described 
. Mice were euthanized by cervical dislocation and the mass of individual tissues and fat depots (heart, kidney, retroperitoneal fat pad, epididymal fat pad, subcutaneous fat pad, and omental fat pad) were determined by dissecting and weighing each tissue/pad separately after the mice were euthanized. Following this, liver tissues were dissected out, flash frozen in liquid nitrogen, and kept at −70 degrees until further processing.
RNA isolation, expression profiling, and RNA-Seq experiment
At 16 weeks of age, the liver tissues of the mice were dissected out, flash frozen in liquid nitrogen, and kept at −70 degrees until further processing. For RNA profiling the RNA from 3 mice per strain were hybridized to Affymetrix Mouse Genome HT_MG-430A arrays. Frozen liver samples were weighed and homogenized in Qiazol according to the manufacturer's protocol. Following homogenization, RNA extraction was performed using Qiagen's RNeasy kit (cat# 74104). Ninety two strains of mice had three biological replicates, five strains had two biological replicates and two strains had one biological replicate each. All RNA samples were cleaned using a Biosprint96 (Qiagen, Valencia, CA) with RNA cleanup beads (Agencourt Bioscience, Beverly, MA) following manufacturer's protocol with adaptations for use with the Biosprint. The quality of the total RNA from the samples was monitored by the Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA) and RNA quantity was measured with a NanoDrop (NanoDrop Technologies, Inc. Wilmington, DE) following the manufacturer's instructions. All samples were arrayed into three 96 well microtiter plates following a randomized design format that places samples from the same strain on different plates to better estimate variance across testing strains. All target labeling reagents were purchased from Affymetrix (Santa Clara, CA). Double-stranded cDNAs were synthesized from 1 ug total RNA through reverse transcription with an oligo-dT primer containing the T7 RNA polymerase promoter and double strand conversion using the cDNA Synthesis System. Biotin-labeled cRNA was generated from the cDNA and used to probe Affymetrix Mouse Genome HT_MG-430A arrays. The HT_MG-430A Array plate consists of 96 single MG-430A arrays arranged into standard SBS 96 well plate format. All cDNA and cRNA target preparation steps were processed on a Caliper GeneChip Array Station from Affymetrix. Array hybridization, washing and scanning were performed according to the manufacturer's recommendations. Scanned images were subjected to visual inspection and a chip quality report was generated by the Affymetrix's GeneChipOperating System (GCOS) and Expression console Affymetrix). Two of 288 chips were excluded due to low QC scores. The image data was processed using the Affymetrix GCOS algorithm utilizing quantile normalization or the Robust Multiarray method (RMA) to determine the specific hybridizing signal for each gene. Expression data can be obtained from Geo database (GSE16780). To avoid the effect of SNP on hybridization, we matched the location of ~14 million SNPs from dbSNP database (NCBI) to the location of the individual probes on the genome. If the location of the probe had a matching SNP within it we flagged the probe and exclude it from the cdf file prior to RMA normalization. If a probeset contained SNP in 8 or more 25-mer probes, we excluded the probeset from the analysis. The cleaned datasets were then background corrected and normalized using the affy package (from bioconductor) using rma, pmonly, and median-polish normalization methods.
RNA isolation for Next Generation Sequencing followed the same protocol as the one described above for the microarray data. For the RNA-Seq experiment, two inbred mice (C57BL/6J and DBA/2J) were chosen. Library preparation was carried using illumina's mRNA-Seq 8-Sample Prep Kit protocol (Illumina, cat# RS-100-0801). In brief, 1 to 10 ug of RNA was used for library construction. In the first step the poly-A containing mRNA molecules were purified using poly-T oligo-attached magnetic beads. Next, the purified mRNA was fragmented into small pieces using divalent cations, followed by double stranded cDNA synthesis using random primers, adenylated at the 3′ end and ligated to the sequencing adapters. The ligated products were then separated on 2% agarose gel, 200 bp fragments were selected and PCR amplified using PE 1.0 and PE 2.0, and purified using QIAquick PCR Purification Kit (QIAGEN, part # 28104). The final library concentration was verified by Bioanalyzer. Sequencing reaction was performed by Illumina Genome Analyzer 2.0 at UCLA Human Genetics microarray core. Raw sequences were uploaded onto Galaxy website (at http://main.g2.bx.psu.edu/
) and using the Tophat software 
was aligned against the reference genome (M. musculus
) downloaded from UCSC. Alignment was performed by setting the parameter for misalignment to one. Relative abundance of transcripts (in Fragments Per Kilobase of Exon per Million read sequence units) was estimated using the Cufflink software 
and the Ensembl's Mus_musculus NCBIM37 as the reference annotation file.
Protein isolation and sample preparation
Male mice were euthanized using isoflurane followed by cervical dislocation at 6–10 weeks of age. The liver tissue was immediately frozen in dry ice until further processing. The 97 samples corresponding to different mouse strains plus some extra samples from C57BL/6J mouse were randomized into 10 batches of 10 samples. Each batch was processed separately prior to quantitative LC-MS analysis. Before the LC-MS analysis the batches were put together and the sample list was randomized one more time. The extraction and digestion of the proteins was performed using a commonly used protocol based on denaturation of protein in 8 M urea followed by digestion with trypsin. Briefly, approximately 5 mg of liver tissue was resuspended in 100 ul denaturing solution (8 M urea, 50 mM Tris-HCl pH 8.0 and 1 mM EDTA) and homogenized with a motorized pestle. Upon homogenization, the total protein content was measured by Bicinchoninic Protein Assay (BCA, Pierce, Rockford, IL) and the 500 ug aliquots were taken from each sample for further processing. DTT was added to a concentration of 10 mM in sample, then to solubilize and unfold the proteins the samples were incubated for 30 min at 37oC with shaking. Cysteine residues were alkylated by adding iodoacetamide up to 40 mM concentration and incubating for 1 hour at 37oC, with shaking, in the dark. For protein digestion the samples were diluted 10-fold with 50 mM ammonium bicarbonate (pH 7.8), supplemented with 1 mM CaCl2, 10 ug of trypsin and incubated for 3 hours at 37oC with shaking. The sample digests were purified with solid phase extraction using C18 columns (Discovery DSC-18, SUPELCO, 52601-U), lyophilized and resuspended in 25 mM ammonium bicarbonate pH 7.8. The peptide amounts were estimated with BCA assay. On the average the amount of purified tryptic peptides was 200 ug. To generate the 18O reference sample, 20 ug of each sample was pooled, then boiled 10 minutes, followed by immediate cooling for 10 minutes. The boiling/cooling steps were performed to inactivate trypsin (this step helps to avoid back-exchange of 18O-labeled peptides). The pooled reference was then subjected to solution-phase tryptic 18O exchange, followed by quenching of tryptic activity with formic acid. The pooled sample was then added in equal amounts with each individual sample for quantitation purposes.
Characterization of the mouse liver proteome
Construction of a library of proteins and tryptic peptides present in the liver is an important step for follow-up quantitation. 10 ug aliquots from all 97 strains were pooled together and subjected to LC fractionation by strong cation exchange (SCX) chromatography on a 200 mm×2.1 mm Polysulfoethyl A column (PolyLC, Columbia, MD) preceded by a 10 mm×2.1 mm guard column, using a flow rate of 0.2 mL/min. LC separations were performed using an Agilent 1100 series HPLC system (Agilent, Palo Alto, CA). Mobile phase solvents consisted of (A) 10 mM ammonium formate, 25% acetonitrile, pH 3.0 and (B) 500 mM ammonium formate, 25% acetonitrile, pH 6.8. Once loaded, isocratic conditions at 100% A were maintained for 10 min. Peptides were separated by using a gradient from 0–50% B over 40 min, followed by a gradient of 50–100% B over 10 min. The gradient was then held at 100% solvent B for another 10 min. Following lyophilization, all thirty fractions collected during this gradient were dissolved in 25 mM ammonium bicarbonate and stored at −80OC.
Each SCX fraction was analyzed with an automated custom-built capillary HPLC system coupled online to an LTQ ion trap mass spectrometer (Thermo Fisher, San Jose, CA) by using an electrospray ionization interface. The reversed phase capillary column was prepared by slurry packing 3-µm Jupiter C18 particles (Phenomenex, Torrance, CA) into a 75 µm i.d.×65 cm fused silica capillary (Polymicro Technologies, Phoenix, AZ). The mobile phase solvents consisted of (A) 0.2% acetic acid and 0.05% TFA in water and (B) 0.1% TFA in 90% acetonitrile. An exponential gradient was used for the separation, which started with 100% A, and gradually increased to 60% B over 100 min. The instrument was operated in a data-dependent mode with an m/z range of 400–2000. Ten most abundant ions from each MS scan were selected for further MS/MS analysis by using a normalized collision energy setting of 35%. Dynamic exclusion was applied to avoid repeat analyses of the same abundant precursor ion.
The SEQUEST software (Thermo Fisher) was used to search the MS/MS data against the mouse International Protein Index (IPI) database (version 3.52 http://www.ebi.ac.uk/IPI
). Human keratins and porcine trypsin were added into the database as expected contaminants. Trypsin cleavage specificity was required for all of the considered peptides. The following criteria were used to filter raw SEQUEST results: 1) Xcorr≥1.9 and DeltaCn2≥0.21 for charge state +1; 2) Xcorr≥2.5 and DeltaCn2≥0.26 for charge state +2; 3) Xcorr≥2.8 and DeltaCn2≥0.32 for charge state +3. These criteria provide the maximum number of peptide identifications not exceeding 1% false discovery rate (FDR). To estimate the FDR of peptide identifications we searched against a reversed database as previously described 
Relative protein abundance quantitation
Relative peptide and protein quantitation was based on ratios between intensities of natural 16
O isotope containing peptides and reference peptides labeled with stable 18
O isotope at the carbonyl group at the C-terminus of the peptide. To create a reference sample we pooled together 20 ug aliquots from all strains and labeled the C-termini with 18
O isotopes using trypsin catalyzed exchange in the presence of heavy H218
O water as described above and elsewhere 
. Prior to the LC-MS analysis 3.75 ug aliquots from each individual sample were mixed with the same amount of 18O-labeled reference sample.
The 7 ug aliquots were analyzed on a LTQ-Orbitrap mass spectrometer that was interfaced with a 75 um i.d.×65 cm long LC column packed with 3 um Jupiter C18 particles (Phenomenex). The mobile phase solvents consisted of (A) 0.2% acetic acid and 0.05% TFA in water and (B) 0.1% TFA in 90% acetonitrile. An exponential gradient was used for the separation, which started with 100% A and gradually increased to 60% B over 100 min. LC-MS datasets were analyzed by in-house software VIPER 
that detected features in mass – elution time space and assigned them to peptides in AMT tag database as described elsewhere 
. Typically an LC-MS run identifies ~3,500 16O/18O peptide pairs that co-elute with a 4.0085 Da mass difference.
As we mentioned before, the relative abundances of tryptic peptides were calculated as the ratio between light and heavy isotopes. The relative abundances then were normalized with EigenMS procedure 
to correct systematic biases that may arise for example from unequal sample loading, batch-to-batch differences in sample processing and LC column variability. Briefly, the EigenMS procedure discovers the systematic trends (so-called eigenpeptides) in the data using singular value decomposition and then removes contributions of those eigenpepides from each peptide. For all data analysis purposes the peptide and protein intensities were log2 transformed and zero-centered by subtracting the peptide or protein specific means taken across all the samples.
To determine the protein levels by immunobloting, liver samples were homogenized in RIPA including phosphatase and protease inhibitors (Santa Cruz Biotech sc-24948), and protein determination were done using the Biorad Dc Assay. Protein samples were boiled following addition of Laemmli loading dye, separated on Invitrogen precast gels, and transferred to PVDF membranes. Membranes were rinsed in 1× TBST (Cell signaling #9997) blocked in 5% skim milk-TBST, rinsed in TBST, and incubated with primary antibodies diluted in 3% BSA-TBST for 1 hr at 23degC or overnight at 4degC. Membranes were washed in TBST and incubated with an HRP-conjugated anti rabbit IgG KPL (#474-1516) 1/5000 in 5% skim milk-TBST. Membranes were washed again, incubated in ECL-plus, and signal detected using a Biorad Chemidoc or film. Densitometry was done using the Biorad Quantity One software. The following list of antibodies and working dilutions were used for each protein: Fasn (Cell Signaling cat #3180, 1/2000), Acyl (Cell Signaling cat #4332, 1/2000), Ywhae (Cell Signaling cat #9635, 1/2000), Vim (Cell Signaling cat#3932, 1/1000), Rkip (Cell Signaling cat#5291, 1/2000), Gapdh (Cell Signaling cat#3683, 1/5,000), Glo1 (Sigma Chemical SAB1100242, 1/20,000), GstA4 (Sigma Chemical SAB1100244, 1/20,000), AnxA5 (Sigma Chemical AV36687, 1/2000), Hao1 (Sigma Chemical AV42480, 1/2000), Aldh3A2 (Sigma Chemical HPA014769, 1/20,000), Actin (Sigma Chemical A2066, 1/5,000), Acox1 (Abnova PAB4367, 1/2000).
For transcript data we applied three filtering steps based on 1) genetic heritability, 2) probeset annotation. We have profiled 3 mice per strain which allowed us to estimate the broad sense heritability for each transcript. Broad sense heritability for each transcript was measured using ANOVA where strain information was used as a grouping factor. The broad sense heritability which is defined as the ratio of genetic variance over total variance for the phenotype was estimated by dividing the sum of squares of the strain information factor over total sum of squares in the ANOVA. The significance of heritability was established if the p-value for the strain information term in ANOVA was below the nominal 0.05 threshold. The selection cutoff for including gene in the analysis based on heritability was set to heritability p-value of 0.05. From the 22700 probesets 10186 probesets did not meet this cutoff. For annotation filtering, we acquired the Ensembl Gene ID for each Affymetrix probeset and selected those probesets that were only annotated to only one Ensembl gene. From the initial 22700 probesets, 4401 probesets had ambiguous annotation (either did not map to a gene or mapped to more than 2 Ensembl genes). Overall, 9896 probesets met both filtering criteria (significant heritability and unique Ensembl annotation). For the protein data, the initial filtering steps were based on 1) eliminating peptides with excessive missing values which would otherwise have unreliable mapping information, 2) eliminating peptides with missed internal cleavage sites which cause unreliable measurement. To annotate peptides we utilized the SpliceCenter web-based tool 
to obtain the location of the exon each peptide represents. Peptides which mapped to multiple exons of more than one gene (as determined by SpliceCenter) were excluded from the analysis because of ambiguous annotation. The genomic exon coordinates for each peptide was then used to query the Ensembl database to acquire the Ensembl Gene IDs. In the second step of filtering peptides that had more than the value of 2 for the ratio of total variance over technical variance were chosen. Overall, 1543 peptides met the two-stage filtering described above.
Genotyping and genome-wide association mapping
Inbred strains were previously genotyped by the Broad Institute (http://www.broadinstitute.org/mouse/hapmap
), and they were combined with the genotypes from Wellcome Trust Center for Human Genetics (WTCHG). Genotypes of RI strains at the Broad SNPs were inferred from WTCHG genotypes by interpolating alleles at polymorphic SNPs among parental strains, calling ambiguous genotypes missing. Of the 140,000 SNPs available, 95,854 were informative with an allele frequency greater than 10% and missing values in less than 10% of the strains. These SNPs were used for both protein and transcript genome-wide association analysis.
We applied the following linear mixed model to account for the population structure and genetic relatedness among strains in the genome-wide association mapping 
In the formula, μ represents mean, x represents SNP effect, u represents random effects due to genetic relatedness with Var(u)
K and Var(e)
, where K represents IBS (identity-by-state) matrix across all genotypes in the HMDP panel. A restricted maximum likelihood (REML) estimate of σg2
are computed using EMMA, and the association mapping is performed based on the estimated variance component with a standard F-test to test β≠0. We applied EMMA (Efficient Mixed Model Association) as an R implementation of a linear mixed model. The percent of variance explained for each molecular phenotype was calculated using the SNP effect calculated from EMMA by defining it as 1-(variance of residuals/variance of original phenotypes). It should be noted that since EMMA is orders of magnitude faster than other implementations commonly used, we were able to perform statistical analyses for all pairs of transcripts and genome wide markers in a few hours using a cluster of 50 processors. Both pQTL and eQTL were defined as “local” if the peak association SNP position was within a 4 Mb interval, flanking 2 Mb on either side of the transcription start and end of the gene under regulation. Genome-wide cutoff:
Genome-wide cutoffs were calculated as the false discovery rates using the “qvalue” package for FDR calculation in the R statistical software 
. Due to the computational complexity associated with evaluating q-values for over 400 million p-values, we computed the FDRs by taking the average FDR for 100 samples each containing 5 million randomly selected p-values from the original calculated p-values. FDR calculation was carried out separately for the protein and transcript dataset.
GO analysis and other statistical methods/software
All statistical analyses and data visualizations were carried out using the R statistical software (available at http://cran.r-project.org/
). Classification of proteins and transcripts to various GO categories was accomplished using Mouse Genome Informatics website at Jackson Laboratories and the GO ontology tool at Princeton. For each probe and each peptide, we first obtained the MGI IDs using the MGI batch query tool at http://www.informatics.jax.org/ 
. Using MGI IDs we utilized the GO Term Mapper at http://go.princeton.edu
, which is based on map2slim algorithm 
to obtain the GO annotations and summary statistics. The background geneset used in this analysis was the list of all genes annotated by MGI. To assess the significance of the correlation coefficients observed in the PT-pair GO analysis, for each GO category we created 100,000 bootstrap datasets each equal in size to the number of genes assigned to the GO term. Bootstrapping was carried out randomly and without replacement from the pool of 584 original correlation p-values among the PT-pairs. The significance of the observed average p-values for each GO term is reported as the two-tailed test against the empirical distribution created by the corresponding 100,000 permutation set. All the correlation coefficients and corresponding p-values reported in the paper are calculated using the bicor function in the WGCNA R package 
. The main advantage of using bicor, which performs biweight midcorrelation calculation, over Pearson's correlation is based the robustness of the correlation coefficient measurement to the presence of outliers in the data.