Proteins constitute much of the structure and functional machinery of cells, forming signaling 23 networks, metabolic pathways, and large multi-component complexes. Protein abundance is 24 regulated at multiple levels spanning transcription, translation, recycling, and degradation to 25 maintain proper balance and optimal function. To better understand how protein abundances 26 are maintained across varying genetic backgrounds, we analyzed liver proteomes of three 27 genetically diverse mouse populations. We observe strong concordance of genetic and sex 28 effects across populations. Differences between the populations arise from the contributions of 29 additive, dominance, and epistatic components of heritable variation. We find that the 30 influence of genetic variation on proteins that form complexes relates to their co-abundance. 31 We identify effects on protein abundance from mutations that arose and became fixed during 32 breeding and can lead to unique regulatory responses and disease states. Genetically diverse 33 mouse populations provide powerful tools for understanding proteome regulation and its 34 relationship to whole-organism phenotypes. 35 36 37 38
Key words 39 proteomics, protein complex, systems genetics, pQTL, Collaborative Cross, Diversity Outbred, 40 inbred strains 41
Protein abundance is regulated at multiple levels spanning transcription, translation, recycling, 43 and degradation. It is responsive to genetic variation, as observed in yeast (Picotti et In this work, we obtained multiplexed mass spectrometry (mass-spec) quantification of proteins 79 in liver samples of 116 CC mice representing female/male pairs from 58 strains. We previously 80 collected proteomics data from the livers of 192 DO mice and 32 mice representing the eight 81 founder strains (two animals of each sex per founder strain) (Chick et al., 2016) (Figure 1a) . We 82 map protein quantitative trait loci (pQTL) in both the CC and DO and find consistency of the 83 genetic and sex effects on protein abundance across the populations. We characterize the 84 genetic architecture of protein-complexes in both the CC and DO, and find that complexes with 85 highly co-abundant members are heritable and exhibit non-additive and/or highly polygenic 86 genetic architecture, which is more apparent in the CC. Finally, we identify genetic variants that 87 originated and became fixed during the of the CC and contribute to CC strain-specific 88 phenotypes. Our work demonstrates the complementary strengths of these genetic resources 89
and offers new insights into the genetic regulation of protein abundance. 90 91 RESULTS 92
We quantified the abundance of 6,779 proteins in liver tissue from 58 inbred CC strains, one 94 female and male per strain. We previously reported quantification of proteins from liver tissue 95 of 192 outbred DO mice and 32 mice representing the eight founder strains (two per sex per 96 strain) (Chick et al., 2016) . The data for DO and founder strains were re-analyzed for this study 97
to ensure that all data were processed consistently (Methods), resulting in the quantification of 98 6,588 and 6,950 proteins, respectively. From the 9,226 total proteins detected, 4,434 were 99 observed in all three populations (Figure 1b ; Table S1 ). 100
We estimated protein abundance heritability (h 2 ), which reflects the combined effects of 101 genetic variants, their allele frequencies, and the large-scale genetic background (e.g., inbred or 102 outbred) of each population (Figure 1c) . Heritability estimates were highest in the CC and 103 founder strains, due to capturing variation from non-additive genetic effects (i.e., recessive and 104 epistatic) in inbred strains. Heritability estimates in the DO are limited to additive genetic 105 factors (narrow sense heritability; Lynch and Walsh 1998) that are captured in the genetic 106 relatedness of our sample. Heritability estimates were significantly correlated across 107 populations (r > 0.4, p <2.2e-16), suggesting that much of the controlling genetic variation is 108 conserved across populations. 109
In order to identify the genetic loci that drive variation in protein abundance, we carried out 110 protein quantitative trait locus (pQTL) mapping in the CC and DO (Table S2) . To determine 111 significant pQTL, we first applied a permutation analysis (Doerge and Churchill, 1996) to control 112 genome-wide error rate for each protein and then applied a false discovery rate adjustment 113 (FDR < 0.1) across proteins (Chesler et al., 2005) to establish a stringent detection threshold for 114
pQTL. Using this criterion, we identified 1,199 local and 289 distal pQTL in the CC and 1,652 115 3 116 (Figure 2a & b) , where local is defined as when the pQTL is 132 located within 10 Mbp of the midpoint of the protein-coding gene. We also identified a local 133 pQTL on the mitochondrial genome in the CC for mt-Nd1 ( Figure S1d) . In order to prevent the  134  exclusion of true pQTL due to overly stringent control of false positive and accurately compare  135 pQTL discovery across populations, we carried out a parallel analysis with more lenient FDR 136 control (FDR < 0.5; Figure S1a & b). 137
We next restricted the set of proteins to the 4,556 that were detected and analyzed in both 138 populations in order to compare genetic effects between the CC and DO. (Figure 2c) . Among 139 1,441 local pQTL stringently detected for the shared proteins, 639 were detected in both 140 populations (Figure 2d ). To determine if these local pQTL were driven by the same genetic 141 variants, we compared the estimated haplotype effects at each pQTL (Methods) and found that 142 622 (97.3%) were significantly positively correlated (Figure 2e ; Table S3 ). To assess whether 143
pQTL detected in only one population are population-specific or fell below the detection 144 threshold in the other population but are likely driven by the same genetic variation, we 145 compared the haplotype effects of detected pQTL to effects estimated at the corresponding 146 locus in the other population, finding that 1,239 (86.0%) were significantly positively correlated. 147
This trend for local pQTL holds true for lenient detection (Figure S1f & g; Table S3 ), indicating 148 that local genetic effects on proteins are highly conserved between the CC and DO populations, 149 even when failing to pass the threshold of detection in one of the populations. 150
The founder strains provide additional support for local pQTL in the CC and DO (Figure S1k-n) , 151 particularly for those that are challenging to detect due to rare alleles in the CC or DO (e.g., a 152 founder allele observed in three or fewer CC strains). We selected all genes with rare local 153 founder alleles that did not have a leniently detected pQTL in the CC, representing 2,410 genes. 154 We correlated the haplotype effects estimated at the locus closest to the gene transcription 155 start sites (TSS) with the protein abundance in the founders (Methods) and found significant 156 positive correlation for 196 genes, such as Cyp2j5. The conservation of local genetic effects 157 among the founder strains, CC, and DO can be striking (such as in Gosr2 and Cyp2d22; Figure  158 S2). The three populations together identify a more complete catalog of local pQTL than any 159 one alone. In total we find evidence to support local genetic effects on abundance for 3,004 160 proteins observed across the CC and DO. 161
Of the 183 distal pQTL stringently detected in the CC and 278 in the DO for the shared set of 162 proteins, we found 19 that were stringently detected in both populations. All 19 also had 163 significantly correlated haplotype effects (Figure 2h ; only stringently detected in only one population, we identified an additional 45 shared distal 167 pQTL between the CC and DO populations (Figure 2i ; Table S3 ). they unobserved, such as proteins that failed to be detected by mass-spec or functional non-198 coding RNA. We identified candidate mediators for each of the 19 shared distal pQTL (Figure 199 2j). The same mediator was identified in the CC and DO for 14 of these -we note that BCKDHB, 200 CCDC93, and IKBKAP are each candidate mediators of the distal pQTL for two proteins (CCDC93 201
for Ccdc22 and Fam45a, BCKDHB for Bckdha and Ppm1k for, and IKBKAP for Elp2 and Elp3). For 202 four of the 19 distal pQTL, the best candidate mediator was different for the CC and DO. For 203 example, TUBGCP3 is the strongest mediator of distal pQTL of Tubg1 and Tubgcp2 in the CC but 204 NAXD is the strongest mediator in the DO (Figure 2l) . Given that Tubg1, Tubgcp2, and Tubgcp3 205 are all members of the tubulin superfamily, TUBGCP3 is a stronger candidate based on 206 functional overlap. A true mediator should possess a local pQTL that co-localizes with the distal 207 pQTL. Notably, the local pQTL for Naxd is much stronger in comparison to the local pQTL of 208
Tubgcp3 in the DO (LODNaxd = 40.2 compared to LODTubgcp3 = 11.0) whereas in the CC, they are 209 more comparable (LODNaxd = 9.7 compared to LODTubgcp3 = 7.9). The local pQTL of Naxd in the 210 DO may be acting as a surrogate for the local genotype, and thus outperforming TUBGCP3. 211
Similarly, COPS8 is the strongest mediator for distal pQTL of Cops6 and Cops7a in the CC, 212
whereas in the DO, COPS8 was not detected as a mediator and notably did not have a local 213 pQTL (even leniently detected) in the DO, suggesting that COPS8 may be less accurately 214 measured in our DO sample population. The final case of discordance was for the distal pQTL of 215 Dnajb9, for which UGT3A2 was identified as a mediator in the DO whereas no candidate was 216 found in the CC. Ugt3a2 possesses local pQTL in both the CC and DO, which is notably strong in 217 the DO (LODUgt3a2 = 36.9), suggesting it may be a false positive mediator in the DO and the true 218 mediator was unobserved for both populations. We considered all distal pQTL that were 219 stringently detected (FDR < 0.1) in one of the populations, and evaluated the corresponding 220 pQTL status (stringent, lenient, or not detected) and mediation status (e.g., same or different 221 mediator) in the other population (Figure S1o) , which revealed similar levels of concordance in 222 mediation for pQTL detected in only one of the populations. 223
When we expand our pQTL comparisons to include all leniently detected pQTL, 1,462 proteins 224 have significantly correlated local pQTL (223 more than with stringent detection in one 225 population; Figure S1g ). In contrast, fewer distal pQTL with significantly correlated effects are 226 detected (22 compared to 64) ( Figure S1j ). Some pQTL that fail to meet the stringent threshold 227
in either population still replicate across the CC and DO, as is the case for 12 of the 41 proteins 228
with distal pQTL and correlated haplotype effects (Figure S1i). One interesting example is Ercc3 229 (Figure S3 ), which has distal pQTL near a region on chromosome 7 that contains Gtf2h1 and 230
Ercc2, which all strongly associate with each other based on protein-protein interactions 231 ( (22.3%) proteins in the founder strains. The differences between male and female were 240 overwhelmingly in the same direction for all populations (Figure 1g -i). Gene set enrichment 241 analysis revealed that proteins related to ribosomes, translation, and protein transport gene 242 ontology (GO) terms were more abundant in male livers whereas proteins related to catabolic 243 and metabolic processes, including fatty acid metabolism, were more abundant in female livers 244 in all populations. 245 246
Drivers of variation in the abundance of protein-complexes 247
Members of protein-complexes exhibit varying degrees of co-abundance, with some groups of 248 proteins being tightly co-abundant, i.e., correlated, and others less so (Romanov et al., 2019) . 249
Correlation between members of a complex suggests some degree of co-regulation. We found 250
that Table S6 ). However, protein-265 complexes can be influenced by genetic variation, as we previously reported for the Chaperonin 266
containing T (CCT)-complex that is stoichiometrically regulated in response to the low 267 abundance of a member protein, CCT6A (Chick et al., 2016) . We evaluated the extent to which 268 members of annotated protein-complexes were internally correlated, as well as how genetic 269 factors and sex contribute to variation in their joint abundance. We summarized co-abundance, 270
which we refer to as complex cohesiveness, as the median pairwise correlation between 271 complex members (Romanov et al., 2019) . To assess the contributions from genetic factors and 272 sex, we summarized each protein-complex using the PC1 from member proteins. We first 273 filtered proteins with local pQTL (FDR < 0.5) or strong distal pQTL (FDR < 0.1) to minimize the 274 influence of individual proteins with independent genetic effects in order to focus on the 275 shared genetic effects on a protein-complex. We estimated heritability and the proportion of 276 variation explained by sex for each complex-specific PC1 (Table S7) . is homozygous in seven CC strains, whereas in our DO cohort, there were no mice homozygous 319
for the PWK allele ( Figure S5b ). Heterozygous carriers of the PWK allele in the DO mice 320 exhibited no discernable effect, suggesting that the exosome complex is regulated by a 321
recessive PWK allele at Exosc7 (Figure S5c) , which is consistent with founder PWK mice having 322 low EXOSC7 abundance ( Figure S5d ). Mediation analysis indicates that genetic variation at the 323
Exosc7 locus distally regulates the complex as well as the functionally related genes, Dis3L and 324
Etf1. These two genes were not included in the original complex annotations, and their co-325 regulation suggests new biological interactions discovered via our approach. The complex-326 heritability (with Dis3l and Etf1 included as well as complex members previously filtered out 327 due to possessing pQTL) was 91.3%, and after removing the seven strains with the PWK allele at 328
Exosc7, was reduced but still high at 61.4%, indicating the presence of additional genetic factors 329 that affect the abundance of the exosome. 330 331
Secondary genetic effect on the chaperonin complex 332
Previously we found that the CCT complex was stoichiometrically regulated and driven by low 333 abundance of CCT6A when the NOD haplotype is present ( The DO sample is well-powered to detect the low NOD allele that drives 336 the complex-wide regulation, due to 19 (9.9%) of the mice being homozygous NOD. The CC 337 strains replicate the distal pQTL at the locus of Cct6a through a low NOD effect for complex 338 members (Cct4, Cct5, Cct8, and Tcp1) -but CCT6a itself was not quantified in the CC samples. 339
The effect of the pQTL at Cct6a drives less of the overall variation in the CC sample due to fewer 340 occurrences of NOD homozygotes (12 CC mice in comparison to 19 DO mice). The CC reveals a 341 secondary genetic effect mediated through CCT4, with high NZO and PWK alleles. After 342 including complex members that were previously filtered out for possessing pQTL, the complex-343 heritability was 56.3%, which, after excluding CC strains that are NOD at Cct6a and NZO or PWK 344
at Cct4, was still significant at 46.6% (Figure 5h ), indicating that, as with the exosome complex, 345 additional genetic effects contribute to CCT complex abundance. inducible forms suggests that although individual samples vary in the relative abundance of the 384 constitutive proteasome and immunoproteasome, they predominantly express one of the 385 forms. Across the founder strains, this relationship appears to be genetically regulated, with the 386 WSB, AJ, and NZO strains expressing more immunoproteasome, and the others expressing 387 more of the constitutive form (Figure S6b-g) . In the recombinant CC and DO, we identified a 388 genetic variant that controls the balance between PSMB6 (constitutive) and PSMB9 389 (immunoproteasome). In both the CC and DO, genetic variation near Psmb9 affected PSMB9 390 abundance, as well as PSMB6 abundance, confirmed through mediation analysis (Figure 6i & j) . 391
Consistent with the anti-correlation between inducible and constitutive subunits as well as the 392 balance in the founder strains, the relationship between PSMB9 and PSMB6 is negative, e.g., 393
mice that inherited the WSB haplotype at the Psmb9 locus have high PSMB9 abundance and 394 low PSMB6 abundance (Figure 6k Genetic factors also influence other components of the 26S proteasome. We identified a strong 403 local pQTL that is consistent in both the CC and DO for Psmd9 that does not affect other 404 members of 19S regulator (Figure 6g & h) , thus explaining the lack of cohesiveness of PSMD9 405 within the proteasome, which was previously observed in the DO (Romanov et al., 2019) . A 406 distal pQTL hotspot comprising members of the 20S catalytic core and 19S regulator mapped to 407 a short interval on chromosome 1 in the CC, most strongly observed in Psma1 (Figure 6m) . No 408 single, strong mediator was detected for the hotspot, though PEX19 was the best candidate for 409 a number of the proteins, suggesting there may be multiple drivers of the hotspot or that the 410 true driver was not observed at the protein-level. 411 412
Polygenic regulation of the mitochondrial ribosomal small subunit 413
The mitochondrial ribosomal small subunit was highly cohesive in both populations (Figures 3a,  414  b, 6n, & o) . The complex-heritability is also high, and after including all annotated members, 415
was 74.0% in the CC and 51.9% in the DO. Despite high complex-heritability, we detected few 416 pQTL for individual members of the complex. The one notable exception is Auh, for which we 417 detected a strong local pQTL in the CC and DO (Figure 6r & s) . population (Figures 6p & q) . This distribution contrasts with the bimodal pattern of abundance 451
for the exosome complex, which is driven by a single strong pQTL (Figure 4h) (Figure 1a) . As a proof of concept, we first confirmed the loss of expression of 462 proteins with known deletions, such as the 80 kbp deletion in CC026 that includes C3 and a 15 463 bp deletion for CC042 in Itgal (Figure 7a- protein's coding gene, and this level of enrichment was significant per permutation (p = 3.7e-4). 482
Interestingly, not all of the observed outliers associated with genetic variants were low 483 extremes, as would be expected with a mutation that results in loss of the protein; we also 484 observed increases in protein abundance in strains with private variants, such as Sash1, which 485 harbors a novel SNP allele in CC058 (Figure 7c ). CC004 has a unique SNP allele associated with 486 low abundance of Plek, a gene which also has a weak local pQTL based on genetic variation 487 from the founder strains (Figure S7e & f) . 488
The presence of biologically related proteins with extreme abundances specific to CC strains 489 further supports the impact of strain-specific genetic variants and regulatory patterns. We 490 identified strain-specific dynamics by testing the set of outlying proteins specific to each CC 491 strain for enrichment in GO and KEGG pathway terms (Tables S8 and S9) . In CC013, we 492 observed increased abundance in proteins with GO annotations for the innate immune system 493 (Figure 7d) , leukocytes, and other immune system-related GO terms. CC013 possesses a unique 494 SNP allele in Hcls1 that was associated with increased HCLS1 abundance (Figure 7e) , a gene 495 involved in myeloid leukocyte differentiation that may contribute to the high abundance of 496 immune-related proteins in CC013. CC013 also expressed a unique liver phenotype, 497 characterized by white granules across the tissue (Figure 7f) , which may relate to the excess of 498 immune-related proteins. CC030 has a 4 bp insertion in Brd4 that was associated with high 499 BRD4 abundance and may contribute to increased abundance in proteins related to nuclear 500 chromosome and chromatin (Figure 7g-h) . The strain-specific protein outlier sets were enriched 501 in a wide range of GO biological functions (Figure S7g-j) sex effects on proteins. In comparison to individual proteins, protein-complexes were less 511 consistent across these populations, highlighting the role of non-additive genetic effects in 512
controlling protein-complex interactions. We observed a wide range of genetic effects on 513 protein complexes, ranging from stoichiometric regulation in response to a large-effect locus 514
(exosome and CCT complexes) to multi-locus and highly polygenic regulation (e.g., the 26S 515 proteasome and mitochondrial ribosomal small subunit, respectively). Within the CC, we 516 observed strong effects on protein abundance from new mutations in specific CC strains. Lastly, 517
we highlight individual CC strains with both aberrant protein regulation based on strain-specific 518 mutations that broadly disrupt protein regulatory networks. 519 520
Consistent with expectations, the effects of local genetic variation are highly conserved 522 between the CC and DO (and often with the founder strains). In cases where discordance 523 occurred between populations, we can often explain based on differing founder allele 524 frequencies between the populations, such as for the local pQTL of Ercc2 that is not detected in 525
the DO because only one NOD homozygote was observed in the DO sample. Differences 526 between populations stemming from different allele frequencies has been observed in human 527
populations (Mogil et al., 2018) . These discordant cases highlight how dominance and large-528 scale homozygosity contribute to the genetic effects on proteins, exposed by comparing the CC 529 and DO. Distal (i.e., trans) pQTL are often harder to detect due to weaker effects and are thus 530 more likely to be discordant based on differences in allele frequencies and the presence of non-531 additive effects. Here we showed that the 19 distal pQTL that were detected in both 532 populations at FDR < 0.1 are highly consistent, both in terms of haplotype effects and their 533 mediation candidates. Discordance in effects between the CC and DO was much greater for 534 distal pQTL detected in a single population, which is in striking contrast to the strong 535 concordance of local genetic effects detected in only one population. The extent of discordance 536 in distal pQTL effects compared to local effects suggests that additional factors are contributing 537 beyond differing allele frequencies and may represent distal genetic effects that are not 538 conserved across population. 539
Mediation analysis results mirror the concordance of distal pQTL effects, with largely the same 540 mediators detected for distal pQTL detected with high confidence in both populations. There 541 are some key caveats with mediation analysis, such as its accuracy being dependent on the true 542 mediator being present in the data. If a candidate mediator is correlated with the true but 543 unobserved mediator, possibly due to LD, it will likely be identified as a false positive mediator. 544
This issue is problematic for proteomics studies if the mediator is not captured in the protein 545 19 data, such as non-coding RNAs or lowly abundant proteins that go undetected in the mass-spec 546 analysis. For example, PGD was detected as a candidate mediator of a distal pQTL of Akr1e1 in 547 both the CC and DO; however, previous studies identified zinc finger proteins (Rex2 and Zfp985) 548
as the most likely mediator of the pQTL (Hamilton-Williams et al., 2010; Keele et al., 2020). Zinc 549 finger proteins are lowly expressed and not prevalent in proteomics data, and in their absence, 550 the nearby protein PGD stands out as the best mediator in both populations. Variable 551 measurement error across candidate mediators can also cause preference for a false mediator 552 due to noise reducing the correlation between the distally controlled protein and its true 553
mediator. This is likely occurring in the DO where NAXD is the stronger mediator for the Tubg1 554 distal pQTL than TUBGCP3, the strong biological candidate. Furthermore, the presence of 555 variable measurement error in proteomics is likely, stemming from the overall magnitude of 556 each protein's expression, the size of the protein and the number of peptides used to 557 summarize it, and how specifically those peptides map to the protein. Comparisons of 558 mediation analysis of two independent genetic experiments can provide replication of findings, 559
as well as assess a more complete set of candidate mediators and correct misidentifications. 560 561
Protein-complexes can be viewed as emergent phenotypes that can be driven by independently 563 regulated members or sub-complexes, as well as higher order dynamics like protein-protein 564 stoichiometry that control complex assembly. We observed a spectrum of cohesiveness 565 (Romanov et al., 2019) across the protein-complexes and our three populations. These 566 differences reflected components and sub-complexes that are semi-independent of the greater 567 complex, as well as the underlying genetic architectures of the populations. Because high 568 cohesiveness does not necessarily imply that a protein-complex is genetically regulated (e.g., 569
stoichiometry could produce tight correlation among complex members due to non-genetic 570 factors), we employed genetic analyses (e.g., QTL and heritability) to characterize genetic 571 control of the protein-complexes. New models of aberrant protein functional networks and disease 601
The unique biology of individual mouse strains has led to a variety of discoveries of genetically 602 based disease phenotypes. We leveraged the CC population to identify strains exhibiting 603 aberrant protein dynamics that could not be tied to the founder strains, such as being 604
downstream of mutations unique to specific CC strains. They may also result from unique 605 combinations of alleles of upstream drivers of a shared functional network. Examples include 606 increased abundance of immune-related proteins and unique liver phenotype in CC013, which 607 we are following up, and altered mitochondrial respiratory complex I function in CC007. Due to 608 the replicability of the CC population, these strain-specific protein networks can be followed up, 609
confirmed, and the underlying mechanisms dissected to reveal new biology and develop new 610 models of disease. 611 612
Integrative genetic resource populations 613
Proteins are a more functionally relevant measure of physiology and disease than the more 614 commonly measured gene transcripts. In our examination of protein regulation across three 615 related populations, we found that local genetic and sex effects on protein abundance were 616 consistent, suggesting that molecular phenotype data and their findings can largely be 617 integrated across these populations. It is easy to conceive of investigators querying specific 618 genes of interest to assess whether their proteins possess sex effects, pQTL and their haplotype 619 effects, or mediators of distal pQTL. To enable these queries, we provide interactive QTL 620 analysis tools for both the CC (https://churchilllab.jax.org/qtlviewer/CC/Ferris) and DO 621 (https://churchilllab.jax.org/qtlviewer/DO/Svenson). In contrast, higher order molecular 622 phenotypes or characteristics, such as protein-complexes and regulatory networks, showed 623 greater discordance across these populations. Non-additive genetic effects, either at a single 624 locus or across loci, were apparent in the regulation of protein-complexes, driving the 625 discordance between populations based on the genetic architecture of each. Lastly, we 626
identified CC strain-specific aberrant protein abundances, their phenotypic or system relevance, 627
and their putative consistency with previously described mutations present across these 628 strains. 629
In this work, we used these diverse mouse populations to finely dissect genetic effects on 630 protein-complexes in liver tissue, revealing the presence of non-additive effects and a polygenic 631 spectrum of genetic regulatory patterns. In the future, we envision highly expandable joint 632 population resources, covering a range of molecular phenotypes (e.g., RNA-seq, ATAC-seq, 633 mass-spec proteomics) and tissues relevant to human disease. proteins were digested by Lys-C at a 1:100 protease-to-peptide ratio overnight at room 691 temperature with gentle shaking. Trypsin was used for further digestion for 6 hours at 37°C at 692 the same ratio with Lys-C. After digestion, 50 µL of each sample were combined in a separate 693 tube and used as the 11 th sample in all 12 tandem mass tag (TMT) 11plex. 100 µL of each 694 sample were aliquoted, and 30 µL acetonitrile (ACN) was added into each sample to 30% final 695 volume. 200 µg TMT reagent (126, 127N, 127C, 128N, 128C, 129N, 129C, 130N, 130C, 130N,  696 and 131C) in 10 µL ACN was added to each sample. After 1 hour of labeling, 2 µL of each sample 697 was combined, desalted, and analyzed using mass-spec. Total intensities were determined in 698 each channel to calculate normalization factors. After quenching using 0.3% hydroxylamine, 11 699 samples were combined in 1:1 ratio of peptides based on normalization factors. The mixture 700 was desalted by solid-phase extraction and fractionated with basic pH reversed phase (BPRP) 701 high performance liquid chromatography (HPLC), collected onto a 96 well plate and combined 702 for 24 fractions in total. Twelve fractions were desalted and analyzed by liquid 703 chromatography-tandem mass spectrometry (LC-MS/MS). 704 705
Off-line basic pH reversed-phase (BPRP) fractionation 706
We fractionated the pooled TMT-labeled peptide sample using BPRP HPLC (Wang et al., 2011) . 707 We used an Agilent 1200 pump equipped with a degasser and a photodiode array (PDA) 708
detector. Peptides were subjected to a 50-min linear gradient from 5% to 35% acetonitrile in 10 709 mM ammonium bicarbonate pH 8 at a flow rate of 0.6 mL/min over an Agilent 300Extend C18 710 column (3.5 μm particles, 4.6 mm ID, and 220 mm in length). The peptide mixture was 711 fractionated into a total of 96 fractions, which were consolidated into 24, from which 12 non-712 adjacent samples were analyzed (Paulo et al., 2016a) . Samples were subsequently acidified with 713 1% formic acid and vacuum centrifuged to near dryness. Each consolidated fraction was 714 desalted via StageTip, dried again via vacuum centrifugation, and reconstituted in 5% 715 acetonitrile, 5% formic acid for LC-MS/MS processing. 716 717
Liquid chromatography and tandem mass spectrometry 718
Mass spectrometric data were collected on an Orbitrap Fusion Lumos mass spectrometer 719 coupled to a Proxeon NanoLC-1200 UHPLC. The 100 µm capillary column was packed with 35 720 cm of Accucore 50 resin (2.6 μm, 150Å; ThermoFisher Scientific). The scan sequence began with 721 an MS1 spectrum (Orbitrap analysis, resolution 120,000, 350−1400 Th, automatic gain control Following acquisition of each MS2 spectrum, we collected an MS3 spectrum in which multiple 728 MS2 fragment ions are captured in the MS3 precursor population using isolation waveforms 729 with multiple frequency notches. MS3 precursors were fragmented by HCD and analyzed using 730 the Orbitrap (NCE 65, AGC 3E5, maximum injection time 150 ms, resolution was 50,000 at 400 731 Th). 732 733
Mass spectra data analysis 734 Mass spectra were processed using a Sequest-based pipeline (Huttlin et al., 2010) . Spectra were 735 converted to mzXML using a modified version of ReAdW.exe. Database search included all 736 entries from an indexed Ensembl database version 90 (downloaded:10/09/2017). This database 737 was concatenated with one composed of all protein sequences in the reversed order. Searches 738
were performed using a 50 ppm precursor ion tolerance for total protein level analysis. The 739 product ion tolerance was set to 0.9 Da. TMT tags on lysine residues, peptide N termini 740 (+229.163 Da), and carbamidomethylation of cysteine residues (+57.021 Da) were set as static 741
modifications, while oxidation of methionine residues (+15.995 Da) was set as a variable 742 modification. 743
Peptide-spectrum matches (PSMs) were adjusted to FDR < 0.01 Gygi, 2007, 2010) . 744 PSM filtering was performed using a linear discriminant analysis (LDA), as described previously 745 (Huttlin et al., 2010) , while considering the following parameters: XCorr, ΔCn, missed cleavages, 746 peptide length, charge state, and precursor mass accuracy. For TMT-based reporter ion 747 quantitation, we extracted the summed signal-to-noise (S:N) ratio for each TMT channel and 748
found the closest matching centroid to the expected mass of the TMT reporter ion. For protein-749 level comparisons, PSMs were identified, quantified, and collapsed to a peptide FDR < 0.01 and 750
then collapsed further to a final protein-level FDR < 0.01, which resulted in a final peptide level 751 FDR < 0.001. Moreover, protein assembly was guided by principles of parsimony to produce the 752 smallest set of proteins necessary to account for all observed peptides. PSMs with poor quality, 753
MS3 spectra with TMT reporter summed signal-to-noise of less than 100, or having no MS3 754 spectra were excluded from quantification ( where local is the effect of the local haplotype pair to peptide for sample , 2012). Here we used a leave-one-chromosome-out or "loco" , in which markers from the 818 chromosome the peptide is predicted to be located on are excluded from estimation in order 819
to avoid the kinship term absorbing some of local[ ] (Wei and Xu, 2016) . We then calculated 820 poly = cor( local , q), the Pearson correlation coefficient between  local , the BLUP of  local and 26 q, the incidence vector of the B6 allele among the founder strains (e.g., q = [0 1 0 0 0 0 0 0] 822 for a peptide with a B6 allele that is missing in the other founder strains Selecting which marker to fit the haplotype effects at is complicated by the fact that the CC and 934
DO have different sets of markers and the genomic coordinate of the peak LOD score is also 935 subject to variation. When comparing pQTL detected in both populations, we fit Equation 4 at 936 the markers with the highest LOD score specific to each population, meaning they may not be 937
the markers closest to each other. When comparing pQTL that were detected in only one 938 population, we selected the marker in the population that failed to map the pQTL that was 939 closest to the marker in the population that detected it. 940 29 941
If the genetic effects on a protein are primarily local, the relative abundances for a protein in 943 the founder strains should match the local pQTL effects observed in the CC and DO, which can 944 be used to confirm findings and better support suggestive pQTL in the CC or DO. We evaluated 945 the consistency of local pQTL in the CC with the founder strains, using an approach similar to 946 how we compared pQTL effects between the CC and DO. For the founder strains, rather than 947 fitting pQTL effects ( QTL ), we fit the founder effects as random terms (as described for the 948 local term in Equation 2 for the founder strains) summarized as BLUPs ( strain Founders ). We then 949 calculated the Pearson correlation between local pQTL effects in the CC and founder effects in 950 the founder strains: local = cor( QTL CC ,  strain Founders ). As when comparing QTL effects between the 951 CC and DO, we then tested local > 0, and corrected for multiple testing through the BH 952 procedure. 953 954
We performed mediation analysis on each distal pQTL with LOD > 6 in the CC or DO, which 956 involved a scan analogous to the QTL genome scans. The underlying model is 957 producing a mediation LOD score. The mediation model is fit for all proteins as mediators, 967 excluding protein , resulting in a mediation scan. 968
We assume that the vast majority of candidates are not true mediators of a specific pQTL and 969 thus the distribution of mediation LOD scores approximates a null distribution. Assuming that 970 the null distribution is approximately normal, we calculate the z-scores of the mediation LOD 971 scores, and define mediators of the QTL at marker for protein as proteins with med < -4, 972
where med is the z-score of the mediation LOD score for candidate mediator protein . The 973 rationale being that when testing the QTL term in Equation 5, if the mediator contains some of 974 the information of the QTL, its presence in both the alternative and null models will result in a 975 large drop in the LOD score of the detected pQTL. In order for a mediator to be a clear 976 candidate driver of the distal pQTL, we required that the mediator TSS be within 10 Mbp of the 977 pQTL marker. Strong mediators that were not near the pQTL often represent proteins that are 978 correlated with protein , such as co-regulated members of a protein-complex. 979 980
Proteins that exhibited differential abundance between the sexes, i.e., sex effects, were 982 identified using an LMM similar to the heritability model (Equation 3) for the CC, DO, and  983 founder strains, but instead testing the significance of the sex coefficient: 984 