Regulation of protein abundance in genetically diverse mouse populations

Summary Genetically diverse mouse populations are powerful tools for characterizing the regulation of the proteome and its relationship to whole-organism phenotypes. We used mass spectrometry to profile and quantify the abundance of 6,798 proteins in liver tissue from mice of both sexes across 58 Collaborative Cross (CC) inbred strains. We previously collected liver proteomics data from the related Diversity Outbred (DO) mice and their founder strains. We show concordance across the proteomics datasets despite being generated from separate experiments, allowing comparative analysis. We map protein abundance quantitative trait loci (pQTLs), identifying 1,087 local and 285 distal in the CC mice and 1,706 local and 414 distal in the DO mice. We find that regulatory effects on individual proteins are conserved across the mouse populations, in particular for local genetic variation and sex differences. In comparison, proteins that form complexes are often co-regulated, displaying varying genetic architectures, and overall show lower heritability and map fewer pQTLs. We have made this resource publicly available to enable quantitative analyses of the regulation of the proteome.


INTRODUCTION
Protein abundance in cells is regulated at multiple levels, including transcriptional and various post-transcriptional, translational, and protein degradation mechanisms. 3 Each of these regulatory mechanisms can be influenced by genetic variation, as observed across a range of organisms, including Arabidopsis, 4 yeast, 5,6 mice, 1,7-9 and humans. [10][11][12][13][14] Genetic effects on protein abundance can be broadly divided into two classes: local and distal. Local variation in the vicinity of the coding gene typically influences protein abundance by altering the rate of transcription or stability of the transcript. 15 Distal genetic variation at loci far from the coding gene typically influences later stages of regulation, often acting through a diffusible intermediate such as another protein. Other modes of regulation are possible, and the local versus distal distinction is a useful but imperfect indicator for distinguishing translational from post-translational regulation. The complexity of genetic regulation is compounded for proteins that form multi-unit complexes because stoichiometry can impose varying degrees of constraint. 1,[16][17][18][19] In each genetic context, proteins can be influenced by a single locus (monogenic) or many (multigenic to polygenic), and effects can be additive or dominant and may involve epistatic interactions with other loci.
Resource populations with high levels of genetic diversity can be used to identify and characterize the genetic loci that affect protein abundance. The Collaborative Cross (CC) 20,21 and Diversity Outbred (DO) 22 mouse populations are two genetic resource populations that are descendant from a common set of eight inbred strains (i.e., the founder strains; short names in parentheses): A/J (AJ), C57BL/6J (B6), 129S1/SvImJ (129), NOD/ShiLtJ (NOD), NZO/HlLtJ (NZO), CAST/EiJ (CAST), PWK/PhJ (PWK), and WSB/EiJ (WSB). The founder strains represent three subspecies of the house mouse, Mus musculus, 23,24 and encompass genetic variation from across laboratory and wild mice. Each DO mouse is genetically unique with high levels of heterozygosity and low linkage disequilibrium (LD) that support fine mapping of genetic variants. 25-28 CC mice consist of more than 60 inbred strains that are homozygous at most loci (>99%). 21,29,30 They have larger LD blocks and thus lower mapping resolution because of fewer outbreeding generations in their derivation than DO mice. The genomes of CC strains are inbred and thus replenishable, enabling repeated measurements of genetically identical mice 31,32 within and across experiments as well as characterization of strain-specific phenotypes. 33 severe acute respiratory syndrome (SARS) coronavirus, 39 and peanut allergy. 40 In this study, we quantified protein abundance in liver samples of 116 CC mice representing female/male pairs from 58 strains (Table S1). We previously collected proteomics data from the livers of 192 DO mice and 32 mice representing the eight founder strains (two animals of each sex per founder strain) 1 ( Figure 1A).
Both studies employed tandem mass tag (TMT) multiplexed mass spectrometry (MS) but represent separate experiments with differences that reflect refinements in the protocols (STAR Methods), most notably use of a pooled bridge sample in each TMT plex for the CC. We compare sex differences and quantitative trait loci for protein abundance (pQTLs) between CC and DO mice, finding strong conservation of sex differences and local  Figure S1 and Tables S1, S2, and S3.
pQTLs and, to a lesser degree, distal pQTLs. We examine proteins that form complexes and find fewer local genetic effects. We highlight examples of protein complexes with diverse genetic architectures. We identify proteins showing unusual expression patterns in specific CC strains and associate some of these with de novo mutations in the CC strains and their phenotypic consequences. Our work demonstrates the consistency of MS proteomics data across experiments and conservation of the regulation of protein abundance across these resource populations.

RESULTS
Heritability and sex effects on proteins are shared across the CC, DO, and founder strains We quantified the abundance of 6,798 proteins (Table S2) in liver tissue from 58 inbred CC strains, one female and one male per strain. We previously reported quantification of proteins from liver tissue of 192 outbred DO mice and 32 mice representing the eight founder strains (two per sex per strain). 1 The data for DO and founder strains were re-analyzed for this study to ensure that all data were processed consistently (STAR Methods), resulting in quantification of 6,745 and 6,897 proteins (Table S2), respectively. Among the 9,235 proteins detected in total, 4,503 were seen in all three populations ( Figure 1B; Table S3).
We estimated protein abundance heritability (h 2 ), which reflects the combined effects of genetic factors relative to the precision of protein abundance estimation ( Figure 1C). Heritability was higher on average in the CC and founder strains compared with the DO strain, likely because of the combined effects of their inbred genetic architecture and improved precision of the MS measurement for the CC samples. Despite the differences in average heritability, heritability of individual proteins was correlated significantly across populations (Pearson correlation coefficient [r] > 0.43, p < 2.2eÀ16), suggesting that the underlying genetic factors are conserved.
Protein abundance can differ between sexes. 1,18 We characterized sex effects in the CC, DO, and founder strains (STAR Methods; Table S3) and detected significant sex effects (false discovery rate [FDR] < 0.1) for 3,750 (55.2%) proteins in the CC strains, 4,520 (67.0%) proteins in the DO mice, and 1,583 (23.0%) proteins in the founder strains. Sex-specific differences in protein abundance were overwhelmingly in the same direction for all populations (Figures 1D-1F). Gene set enrichment analysis revealed that proteins related to ribosome, translation, and protein transport Gene Ontology (GO) terms were more abundant in males, 18 whereas proteins related to catabolic and metabolic processes, including fatty acid metabolism, were more abundant in females in all populations. 41 Genetic regulation of proteins is shared between the CC and DO To identify the genetic loci that regulate variation in protein abundance, we carried out pQTL mapping in CC and DO mice (Table  S4). To determine significant pQTLs, we first applied a permutation analysis 42 to control the genome-wide error rate for each protein and then applied an FDR adjustment (FDR < 0.1) across proteins 43 to establish a stringent detection threshold for pQTLs. Using this stringent criterion, we identified 1,087 local and 285 distal pQTLs in CC mice and 1,706 local and 414 distal pQTLs in DO mice (Figures 2A and 2B). We defined local pQTLs as being located within 10 Mbp of the midpoint of the protein-coding gene. Although this wide local window may result in some nearby but distal pQTLs being misclassified as local, it accounts for the large LD blocks in the CC strains and yields more consistent classification of local pQTLs between CC and DO mice. We also identified a local pQTL on the mitochondrial genome in CC mice for mt-Nd1 (Figure S1D). Stringent control of false positive rates can result in a high rate of false negative results. Therefore, to compare pQTL discovery across populations, we carried out a parallel analysis with more lenient FDR control (FDR < 0.5; Figures S1A and S1B).
We compared genetic effects between CC and DO mice by focusing on the 4,654 proteins that were detected in both populations ( Figure S1C). Among 1,427 local pQTLs detected in either population, 636 were detected in both ( Figure 2C). To determine whether the shared local pQTLs were driven by the same genetic variants, we compared the estimated haplotype effects at each (I) The 20 distal pQTLs detected in CC and DO mice. Arrows connect candidate drivers identified through mediation analysis to their targets (proteins with distal pQTLs). Gene names in black represent the top candidates identified in both CC and DO mice. Red and blue gene names indicate top candidates specific to CC or DO mice, respectively. TUBGCP3, a top candidate from the CC and the stronger biological candidate based on shared membership in protein families, is underlined. The red asterisk denotes PGD as a likely false positive mediator because of the true mediator being unobserved, seen in both CC and DO mice. (J) Mediation analysis of the Snx4 distal pQTL, an example of agreement between CC and DO mice. Panels showing pQTL LOD scores for CC (pink) and DO (blue) mice are overlayed with mediation conditional LOD scores (gray dots). Mediation scores were evaluated for all proteins genome-wide. Candidate mediators of interest with low mediation scores are labeled. Horizontal lines at a LOD score of 6 are included as a reference point across genome scans. (K) Causal diagram consistent with the relationships revealed by QTL and mediation analysis for SNX4 and SNX7. (L) Mediation analysis of the Tubg1 distal pQTL, an example of disagreement between CC and DO mice, details as described above. (M) Causal diagram consistent with the relationships revealed by QTL and mediation analysis for TUBG1, TUBGCP3, TUBGCP2, and NAXD. See also Figures S1-S3 and Tables S4 and S5. pQTL (STAR Methods) and found that 628 (98.7%) were significantly positively correlated (FDR < 0.1; Figure 2D; Table S4). To assess whether pQTLs detected in only one population are population specific, we compared the haplotype effects of detected pQTLs with effects estimated at the corresponding locus in the other population regardless of significance and found that 1,314 (92.1%) were significantly positively correlated (FDR < 0.1; Figure 2E). The concordance of local pQTLs also holds for lenient detection (Figures S1F-S1G; Table S4). Based on these analyses, we find that local genetic effects on proteins are highly conserved between the CC and DO populations.
The founder strains can provide additional support for local pQTLs in CC and DO mice ( Figures S1K and S2), particularly for those that are hard to detect because of rare alleles in CC or DO mice (e.g., observed in three or fewer CC strains). We selected all genes with a rare local founder haplotype that did not have a leniently detected pQTL in CC mice, representing 2,439 genes. We correlated the haplotype effects estimated at the locus closest to the gene transcription start sites (TSSs) with the protein abundance in the founders (STAR Methods) and found significant positive correlation for 194 genes (FDR < 0.1). The three populations together provide evidence of local genetic effects at 2,905 proteins.
Of the 186 distal pQTLs detected in CC mice and 294 in DO mice, 20 were detected in both populations, all with significantly correlated haplotype effects (FDR < 0.1; Figure 2G; Table S4). Overall, the distal pQTLs are weaker than the local pQTLs, 1,44 which may contribute to an increased rate of false negatives. By comparing haplotype effects for pQTLs that were detected in only one population, we identified an additional 55 shared distal pQTLs ( Figure 2H). Based on the stringent criteria, a total of 75 distal pQTLs had consistent effects of the total (16.3%) detected in CC or DO mice compared with 19 (0.5%) for lenient criteria ( Figure S1J). The reduction in concordance for the lenient criteria is likely due to increased numbers of false positives and weak distal pQTLs, although it did uniquely identify a shared distal pQTL for Ercc3 ( Figure S3).
Mediators of strong distal genetic effects detected in CC and DO mice are concordant The genetic variants that drive distal pQTL are generally thought to act through diffusible intermediates that are under local genetic control at the pQTL. We used mediation analysis 1,45-47 to identify candidate mediators of distal pQTLs (STAR Methods; Table S5). Mediation analysis is best used for prioritizing candidate mediators and, in this study, is limited to evaluating candidates that exert their effects through changes in protein abundance. Therefore, we cannot exclude the possibility of mediation by proteins that were not detected, by non-coding RNAs or by protein variants that affect function without altering abundance. For example, in this study, PGD was found to mediate a strong Akr1e1 distal pQTL in CC and DO mice. However, Zfp985 has been identified previously as the mediator based on gene expression in CC mice 47 but was not detected at the protein level in this study.
We identified candidate mediators for each of the 20 shared distal pQTLs ( Figure 2I), of which only four had best candidate mediators that differed between CC and DO mice. For example, TUBGCP3 is the strongest candidate mediator of the distal pQTL for Tubg1 in CC mice, but NAXD is the strongest candidate mediator in DO mice ( Figures 2L and 2M). Given that Tubg1 and Tubgcp3 (as well as Tubgcp2, which also has a co-mapping distal pQTL) are members of the tubulin superfamily, TUBGCP3 is a strong functional candidate, suggesting that NAXD is likely a false positive mediator in DO mice.
We also examined candidate mediators for all distal pQTLs that were detected (FDR < 0.1) in only one of the populations and evaluated the corresponding pQTL status (stringent, lenient, or not detected) and mediation status (e.g., same or different) in the other population. We found the same candidate mediator between CC and DO mice for 21 of 460 distal pQTLs mapped across both populations, suggesting that mediation is more accurate for strong distal pQTLs that are detected in both populations ( Figure S1L).

Drivers of variation in the co-abundance of protein complexes
Members of protein complexes exhibit varying degrees of coabundance, to which we refer as cohesiveness of the complex. We quantify cohesiveness as the median Pearson correlation between complex members. 18 A high level of cohesiveness suggests that co-regulation of protein abundances across a complex is maintaining stoichiometry. We found that individual proteins that are part of a complex 48-50 (Table S6) are less heritable and have fewer pQTLs than proteins that are not part of a complex ( Figures S4A-S4D). We evaluated the extent to which members of protein complexes were inter-correlated as well as how genetic factors and sex contribute to variation in their joint abundance. To assess the contributions from genetic factors and sex, we performed principal-component analysis (PCA) on the abundances of proteins for each protein complex 51 and took the first principal component (PC1) as a summary (STAR Methods). We then estimated heritability and the proportion of variation explained by sex for each complex PC1.
Protein complex cohesiveness was correlated between CC and DO mice (r = 0.68, p < 2.2eÀ16) ( Figures 3A and 3B), and within each population, it is correlated with complex heritability (r = 0.28, p = 1.2eÀ4 in CC mice and r = 0.12, p = 0.11 in DO mice) ( Figure S4E), suggesting that cohesiveness reflects some degree of shared genetic regulation in CC mice. Complex heritability is consistently higher in CC mice than in DO mice (115 of 155 complexes, 74.2%) and is uncorrelated with complex heritability in DO mice (r = 0.09, p = 0.26) ( Figures 3C and 3D). This lack of correlation between heritability of complexes for the CC and DO mice contrasts with the highly correlated heritability of individual proteins (r = 0.44, p < 2.2eÀ16). The proportion of variation in complex abundance (as summarized by PC1) explained by sex is correlated between CC and DO mice (r = 0.70, p < 2.2eÀ16; Figures 3E and 3F). Protein complexes that have been shown previously to have sex-specific abundance in DO mice, 18 such as eIF2B, were confirmed in CC mice. See also Figure S4 and  Figure S5D). Among the DO mice with one copy of the PWK haplotype, there is no reduction in the abundance of the exosome complex, which suggests that the PWK haplotype effect on Exosc7 is recessive ( Figure S5E). We also note that inbred founder PWK mice have low EXOSC7 abundance ( Figure S5F). Mediation analysis identifies EXOSC7 as a candidate distal regulator of the complex as well as the functionally related genes Dis3l and Etf1. These two genes were not included in the complex annotations, but our findings suggest that they are maintained in stoichiometric balance with the annotated complex members. The complex heritability (with Dis3l and Etf1 now included) was 91.9%, and after removing the seven CC strains with the PWK haplotype at Exosc7, it was reduced but still high at 65.1% ( Figure 4E), indicating the presence of additional genetic factors that affect the abundance of the exosome.
Secondary genetic effects on the chaperonincontaining T (CCT) complex Previously we reported that the CCT complex was stoichiometrically regulated by low abundance of CCT6A when the NOD haplotype is present. 1 Figures 3C and 3D). The DO sample includes 19 (9.9%) mice homozygous for the NOD haplotype ( Figure S6C). The CC strains, six of which are homozygous for NOD at the Cct6a locus ( Figure S6D), replicate this distal pQTL for the complex members Cct4, Cct5, Cct8, and Tcp1, although CCT6a itself was not detected in the CC samples. The effect of the pQTL at Cct6a in CC mice drives less of the overall variation in the CCT   Figure S6H), indicating that, as with the exosome, additional genetic effects contribute to CCT complex abundance.
Independent genetic effects on the subcomplexes of the 26S proteasome The 26S proteasome is composed of a 20S proteasome catalytic core (PSMA and PSMB proteins) that, in the constitutive form, incorporates subunits PSMB5, PSMB6, and PSMB7 and is capped by two 19S regulators (composed of the PSMC and PSMD proteins). The constitutive form can be modified by replacing the PSMB subunits with the three immunoproteasomeinducible subunits (PSMB8, PSMB9, and PSMB10) and the 19S regulators with the 11S regulators, composed of PSME proteins 52 ( Figure 5A). The immunoproteasome is a highly efficient form of the proteasome that is predominantly, but not exclusively, expressed in immune cells. 53 This alternation between two different forms of the proteasome is apparent in the correlations among the inducible and immune components in the CC, DO, and founder strains ( Figures 5B, 5C, and 5I). Individual mice appear to predominantly express one of the proteasome forms, as suggested by the anti-correlation between the constitutive and inducible components. Across the founder strains, this dynamic appears to be regulated genetically, with the WSB, AJ, B6, and NZO strains expressing more immunoproteasome components and the other founder strains expressing more of the constitutive components (Figures S7). In CC and DO mice, we identified genetic variation that controls the balance between PSMB6 (constitutive) and PSMB9 (immunoproteasome) (Figures 5F and 5G). The WSB haplotype at Psmb9 drives higher PSMB9 abundance as well as lower abundance of its constitutive analog PSMB6, confirmed through mediation analysis (Figures 5H and 5I). The Psmb9 local pQTL only appears to drive the balance between PSMB9 and PSMB6 and does not directly affect the other interchangeable members of the proteasome (PSMB5, PSMB7, PSMB8, and PSMB10), which do not map strong pQTLs. The distinct correlation patterns among the interchangeable components suggests that they are still co-regulated across the three populations.
Some of the other non-interchangeable components of the 26S proteasome are regulated by genetic variation independent of the Psmb9 locus. We identified a strong local pQTL for Psmd9 that is present in CC and DO mice and does not affect other members of the 19S regulator ( Figures 5D and 5E), which explains the lack of cohesiveness of PSMD9 with the rest of the proteasome, as noted previously in these DO mice. 18 Polygenic regulation of the mitochondrial ribosomal small subunit (MRSS) The MRSS is highly cohesive in CC and DO mice ( Figures 3A, 3B, 6A, and 6B). Complex heritability is also high in CC (73.4% [67.3%-79.1%]) and DO (44.8% [22.2%-70.6%]) mice. Despite its high complex heritability, we detect few pQTLs for individual members of the complex. One exception is Auh, which has a local pQTL in CC and DO mice ( Figures 6E and 6F). AUH is not a core member of the MRSS but has been associated with it 54 and is included in the annotations. AUH's local pQTL and lack of cohesion with the core MRSS proteins indicate that it is regulated separately from the core MRSS. Similarly, RPS15 and PPME1 were annotated with the complex but not correlated with the core proteins, suggesting that they are also not co-regulated with the MRSS. MRPS27 and MRPS28, on the other hand, were missing from the annotations and are thus not included in Figure 6, but we found them to be highly cohesive with the core MRSS. For the core MRSS proteins, CC strains have an overall abundance that is highly strain specific, as represented in the complex PC1 and even for individual proteins (r = 0.82 between males and females for the complex PC1 and r = 0.76 for MRPS7). Furthermore, the variation across CC strains is highly continuous ( Figures 6C  and 6D). This distribution contrasts with the bimodal abundance pattern for the exosome complex, which is driven mostly by a single strong pQTL ( Figure 4E), suggesting that many loci with small effects influence the overall abundance of the MRSS.
CC strain-specific variation affects protein abundance Inbred mouse strains accumulate mutations that can have phenotypic consequences. [55][56][57][58] New mutations have arisen and become fixed in the CC strains 29,30 ( Figure 1A), which may affect protein abundance. We confirmed functional effects on protein abundance for genes with CC strain-specific deletions, including the 80-kbp deletion in CC026 that includes the C3 coding gene and a 15-bp deletion in the Itgal gene (Figures 7A and 7B) that occurred in CC042 and has been shown to increase susceptibility to tuberculosis 59 and Salmonella. 60 We estimated strain-specific abundance levels for every protein detected in CC mice and identified CC strains where the male and female had a distinctly low or high abundance of a given protein (|Z score | > 2.5; STAR Methods), which we refer to as strain-specific protein outliers. In total, we identified 5,907 strain-specific protein outliers representing 4,267 proteins across all 58 CC strains. Of these, 67 strain-specific protein outliers occur in strains with a unique genetic variant in or near the coding gene. 30 Furthermore, not all genes with a strain-specific protein outlier and matching strain-specific genetic variants were associated with low protein abundance, as is the case for the deletions; we observed high abundance associated with strain-specific variants, such as Sash1, which harbors a novel SNP allele in CC058 ( Figure 7C).
Outlying protein abundance patterns specific to a CC strain may represent larger biological pathway dynamics that result from the strain's genetic background. We defined sets of proteins for each strain based on having low or high strain-specific protein outliers and observed strain-specific enrichments for biological functions based on GO and KEGG pathway terms (Table S7). In CC013, we observed increased abundance in proteins related to the innate immune system ( Figure 7D), leukocytes, and other immune system-related GO terms. CC013 possesses a unique SNP allele in Hcls1 that was associated with increased HCLS1 abundance (Figure 7E), a gene involved in myeloid leukocyte differentiation that may contribute to the high abundance of these immune-related proteins. During the process of tissue collection, it was noted that CC013 had a unique liver phenotype characterized by white granules throughout the tissue. We examined additional mice to confirm this ( Figure 7F) and hypothesize that the liver granules are related to an excess of immune-related proteins. Additional CC strains (Figures S8C-S8E) with multiple outlier proteins that are functionally related include CC007 ( Figure 7G), which has low-and high-abundance proteins in mitochondrial respiratory complex I. The replenishable inbred CC strains capture these dynamics and allow deeper interrogation of unique strain-specific networks of functionally related proteins with perturbed abundance and a better understanding of their phenotypic consequences.

DISCUSSION
We carried out proteomics profiling of mice from the CC and DO strains and their founder strains. Despite the challenges imposed by separate experiments and the relative nature of MS proteomics, 61 the data provided consistent results that supported comparative analysis. Genetic regulation of the proteome in the liver is highly conserved across these distinct but related genetic reference populations. The concordance is exceptionally strong for local genetic regulation and sex differences but also for distal genetic effects strongly detected in both populations. Discordance between CC and DO mice can often be attributed to chance differences in allele frequencies or to dominance effects that manifest differently in the inbred versus outbred populations, similar to what has been observed in Drosophila for gene expression. 62 Mediation analysis of distal pQTLs identified many of the same candidate mediators in CC and DO mice. Proteins that form complexes are generally less affected by local ge-netic regulation compared with other proteins. Complexes displayed a wide range of cohesiveness that was more highly conserved across populations than complex heritability; nevertheless, genetic effects on protein complexes can manifest in remarkably different ways. The CC strains enable discovery of extreme strain-specific abundance of individual proteins and of functionally related groups of proteins, which can be recaptured and studied further because inbred strains are replenishable.

Implications
Protein complex members are often synthesized in proportions that are consistent with stoichiometric balance. 17 Genetic perturbations of one or more member proteins can introduce imbalances that need to be compensated, typically  Article ll OPEN ACCESS through degradation pathways that recycle unassembled units of protein complexes. In this way, a genetic variant that severely reduces transcription of a gene and, consequently, the protein abundance of a single complex member can have effects that propagate through the entire complex and be detected as a shared distal pQTL, as was the case for the exosome and CCT complexes. At the other extreme, the abundance of a highly cohesive complex such as the MRSS can be heritable with few or no detectable pQTLs, which is consistent with polygenic regulation by multiple small-effect loci. Even the exosome and CCT complex display significant residual heritability after accounting for their largeeffect pQTL, indicating that polygenic effects on complex abundance are pervasive. We observed some inconsistent results between CC and DO mice using the PC1 to estimate complex heritability and sex effects. This is likely due to differences in how the main axes of variation differ, suggesting that a single summary measure is insufficient to capture the behavior of many complexes. Examination of the entire correlation matrix of protein complexes can reveal a more detailed picture of the regulatory structure, such as internal heterogeneity in the relative balance of components. For example, our analysis of the 26S proteosome reveal known subcomplexes and individual proteins that are regulated independently. In addition, our analysis of the MRSS highlights some shortcomings of current protein complex definitions that can potentially be corrected based on the correlation structure in shotgun proteomics data.
De novo mutations specific to individual CC strains are clearly responsible for outlying abundance patterns for some proteins, but we identified a large number of such outliers, and it seems implausible that the majority of these would be due to mutations. Based on the functional similarity of protein outliers within specific strains, we propose that these represent perturbations of interacting networks of proteins, whether they are due to de novo mutations or to multi-locus allelic combinations that are fixed in specific CC strains. Epistasis, particularly among interacting proteins, could contribute to these CC strain-specific networks. Regardless of their underlying origin, CC strains with single or functional groups of protein outliers can serve as models for further investigation of biological mechanisms and disease.

Limitations of study
Our study has insufficient power to reliably detect and characterize small distal genetic effects, which likely contributes to the reduced concordance of distal pQTLs between CC and DO mice. This study also highlights some of the caveats of mediation analysis. Successful mediation analysis requires that the true mediator is present in the data and that its effects are mediated through variation in abundance and not through other functional changes in the protein. A protein that is correlated with the true but unobserved mediator may be identified incorrectly as a candidate mediator. In addition to unobserved proteins, other factors, such as non-coding RNAs that could mediate distal pQTLs, may not be measured for a given experiment. Proteins with strong local pQTLs can appear to be mediators, as was likely the case for mediation of the Tubg1 distal pQTL through NAXD in DO mice. Comparison of mediation analysis across independent genetic experiments can correct and refine candidate mediators.

Outlook
Unbiased profiling of the proteome provides a unique window into the molecular processes that are active in cells and tissues, a view that is complementary to and often more directly relevant to function than transcriptome profiling. Although many proteins are responsive to transcriptional regulation, they can also be regulated by a variety of post-translational mechanisms. Analysis of proteomics data in genetically diverse populations provides causal perturbations in the form of genetic variation that introduce variability in protein abundance across all levels of regulation. Genetic mapping and correlation analyses can identify co-regulated proteins and key drivers that regulate other proteins and protein complexes. Technologies that measure protein abundance are developing rapidly but are already capable of delivering accurate and reproducible data. Together with the demonstrated consistency of genetic effects on proteins across distinct but related mouse resource populations, this suggests that we can extrapolate findings across these genetic reference populations with some confidence. Resource data such as described here can be co-analyzed with future data through meta-analyses; for example, comparing genetic effects across different tissues. These findings suggest that imputation of locally regulated proteins could be an option when direct profiling of proteins is not available.

STAR+METHODS
Detailed methods are provided in the online version of this paper and include the following: (F) CC013 has a unique liver phenotype characterized by white granules, indicated with a red arrow. (G) Abundance of proteins related to the mitochondrial respiratory chain complex I are shown for female and male CC mice with the outlier strain CC007 indicated. See also Figure S8 and

RESOURCE AVAILABILITY
Lead contact Further information and requests for resources should be directed to and will be fulfilled by the lead contact Gary Churchill (gary. churchill@jax.org).

Materials availability
The CC strains used in this study (key resources table) are available from the UNC Systems Genetics Core (https://csbio.unc.edu/ CCstatus/index.py?run=availableLines). Many of the strains are also available from the Jackson Laboratory.

Data and code availability
The mass-spec proteomics data for the CC liver samples reported here have been deposited in ProteomeXchange (http://www. proteomexchange.org/) via the PRIDE partner repository (ProteomeXchange: PXD018886). This study also makes use of existing, publicly available liver proteomics data from the DO and founder strains (ProteomeXchange: PXD002801). 1 All analyses were performed using the R statistical programming language (v3.6.1). 2 The analysis pipeline used to generate the results, starting from the raw data, scripts to process the raw data, the processed data, and scripts to analyze the processed data and generate the figures, has been made publicly available (figshare: https://doi.org/10.6084/m9.figshare.12818717).

Mice
We received pairs of young mice from 58 CC strains from the UNC Systems Genetics Core Facility between the summer of 2018 and early 2019. Mice were singly housed upon receipt until eight weeks of age. More information regarding the CC strains can be found at https://csbio.unc.edu/CCstatus/index.py?run=availableLines. Mouse studies were performed in strict accordance with the recommendations in the Guide for the Care and Use of Laboratory Animals of the National Institutes of Health. All mouse studies at University of North Carolina at Chapel Hill (UNC) (Animal Welfare Assurance #A3410-01) were performed using protocols approved by the UNC Institutional Animal Care and Use Committee (IACUC) in a manner designed to minimize pain and suffering in infected animals. Any animals that exhibited severe disease signs was euthanized immediately in accordance with IACUC approved endpoints. Mice were kept on a 12h:12h light:dark cycle, with temperature maintained between 68 F and 74 F and 30% and 70% humidity. Mice were provided unrestricted access to food (LabDiet(R) Select Verified Rodent Diets 5V0F -Select Rodent 50 IF/6F auto) and water. Mice were fasted for 6 h before euthanasia and tissue collection.

METHOD DETAILS
Mouse genotyping, founder haplotype reconstruction, and gene annotations The 116 CC mice were genotyped on the Mini Mouse Universal Genotyping Array (MiniMUGA), which includes 11,125 markers. 66 Founder haplotypes were reconstructed using a Hidden Markov Model (HMM), implemented in the qtl2 R package, 65 using the ''risib8'' option for an eight founder recombinant inbred panel. Heterozygous genotypes were omitted, and haplotype reconstructions are limited to homozygous states, smoothing over a small number of residual heterozygous sites that remain in the CC mice. The genotyping and haplotype reconstruction for the DO mice were previously described; 1 briefly, genotyping was performed on MegaMUGA (57,973 markers), 67 and founder haplotypes were reconstructed using the DOQTL R package. 68 Ensembl version 91 gene and protein annotations were used in the CC, whereas version 75 was previously used in the DO and founder strains data. If the gene symbol or gene ID differed for a protein ID between versions 75 and 91, we updated them to version 91 in the DO and founder strains. When comparing results (e.g., heritability, sex effects, and pQTLs) between the CC, DO, or founder strains, we merged based on protein ID. For comparing the more complicated mediation analysis, we allowed matches based on mediator gene symbol rather than mediator protein ID if the target protein IDs matched.

Sample preparation for proteomics analysis
We analyzed liver tissue in the CC to match the previously collected liver data in the DO and founder strains. Sample preparation and mass spectr (MS) analysis for the DO and founder strains were previously described. 1 Singly housed CC mice had their food removed six hours prior to euthanasia and tissue harvest. Tissues were dissected, weighed, and snap frozen in liquid nitrogen. Pulverized CC liver tissue were syringe-lysed in 8 M urea and 200 mM EPPS pH 8.5 with protease inhibitor and phosphatase inhibitor. BCA assay was performed to determine protein concentration of each sample. Samples were reduced in 5 mM TCEP, alkylated with 10 mM iodoacetamide, and quenched with 15 mM DTT. 200 mg protein was chloroform-methanol precipitated and re-suspended in 200 mL 200 mM EPPS pH 8.5. The proteins were digested by Lys-C at a 1:100 protease-to-peptide ratio overnight at room temperature with gentle shaking. Trypsin was used for further digestion for 6 hours at 37 C at the same ratio with Lys-C. After digestion, 50 mL of each sample were combined in a separate tube and used as the 11 th sample in all 12 tandem mass tag (TMT) 11plex. 100 mL of each sample were aliquoted, and 30 mL acetonitrile (ACN) was added into each sample to 30% final volume. 200 mg TMT reagent (126, 127N, 127C, 128N, 128C, 129N, 129C, 130N, 130C, 130N, and 131C) in 10 mL ACN was added to each sample. After 1 hour of labeling, 2 mL of each sample was combined, desalted, and analyzed using MS. Total intensities were determined in each channel to calculate normalization factors. After quenching using 0.3% hydroxylamine, 11 samples were combined in 1:1 ratio of peptides based on normalization factors. The mixture was desalted by solid-phase extraction and fractionated with basic pH reversed phase (BPRP) high performance liquid chromatography (HPLC), collected onto a 96 well plate and combined for 24 fractions in total. Twelve fractions were desalted and analyzed by liquid chromatography-tandem mass spectrometry (LC-MS/MS).

Offline basic pH reversed-phase (BPRP) fractionation
We fractionated the pooled TMT-labeled peptide sample using BPRP HPLC. 69 We used an Agilent 1200 pump equipped with a degasser and a photodiode array (PDA) detector. Peptides were subjected to a 50-min linear gradient from 5% to 35% acetonitrile in 10 mM ammonium bicarbonate pH 8 at a flow rate of 0.6 mL/min over an Agilent 300Extend C18 column (3.5 mm particles, 4.6 mm ID, and 220 mm in length). The peptide mixture was fractionated into a total of 96 fractions, which were consolidated into 24, from which 12 non-adjacent samples were analyzed. 70 Samples were subsequently acidified with 1% formic acid and vacuum centrifuged to near dryness. Each consolidated fraction was desalted via StageTip, dried again via vacuum centrifugation, and reconstituted in 5% acetonitrile, 5% formic acid for LC-MS/MS processing.

Liquid chromatography and tandem mass spectrometry
Mass spectrometric data were collected on an Orbitrap Fusion Lumos mass spectrometer coupled to a Proxeon NanoLC-1200 UHPLC. The 100 mm capillary column was packed with 35 cm of Accucore 50 resin (2.6 mm, 150Å ; ThermoFisher Scientific). Peptides were separated using a 2.5 h gradient of 9$35% acetonitrile gradient in 0.125% formic acid with a flow rate of $400nl minÀ1. The scan sequence began with an MS1 spectrum (Orbitrap analysis, resolution 120,000, 350À1400 Th, automatic gain control (AGC) target 5E5, maximum injection time 50 ms). SPS-MS3 analysis was used to reduce ion interference. 71 Article ll OPEN ACCESS 60 ms), and isolation window at 0.5. Following acquisition of each MS2 spectrum, we collected an MS3 spectrum in which multiple MS2 fragment ions are captured in the MS3 precursor population using isolation waveforms with multiple frequency notches. MS3 precursors were fragmented by HCD and analyzed using the Orbitrap (NCE 65, AGC 3E5, maximum injection time 150 ms, resolution was 50,000 at 400 Th).
Mass spectra data analysis Mass spectra were processed using a Sequest-based pipeline. 73 Spectra were converted to mzXML using a modified version of ReAdW.exe. Database search included all entries from an indexed Ensembl database version 90 (downloaded:10/09/2017). This database was concatenated with one composed of all protein sequences in the reversed order. Searches were performed using a 50 ppm precursor ion tolerance for total protein level analysis. The product ion tolerance was set to 0.9 Da. TMT tags on lysine residues, peptide N termini (+229.163 Da), and carbamidomethylation of cysteine residues (+57.021 Da) were set as static modifications, while oxidation of methionine residues (+15.995 Da) was set as a variable modification.
Peptide-spectrum matches (PSMs) were adjusted to FDR < 0.01. 74,75 PSM filtering was performed using a linear discriminant analysis (LDA), as described previously, 73 while considering the following parameters: XCorr, DCn, missed cleavages, peptide length, charge state, and precursor mass accuracy. For TMT-based reporter ion quantitation, we extracted the summed signal-to-noise (S:N) ratio for each TMT channel and found the closest matching centroid to the expected mass of the TMT reporter ion. For protein-level comparisons, PSMs were identified, quantified, and collapsed to a peptide FDR < 0.01 and then collapsed further to a final protein-level FDR < 0.01, which resulted in a final peptide level FDR < 0.001. Moreover, protein assembly was guided by principles of parsimony to produce the smallest set of proteins necessary to account for all observed peptides. PSMs with poor quality, MS3 spectra with TMT reporter summed signal-to-noise of less than 100, or having no MS3 spectra were excluded from quantification. 76

QUANTIFICATION AND STATISTICAL ANALYSIS
Filtration of peptides that contain polymorphisms Peptides that contain polymorphisms are problematic for protein quantification in genetically diverse samples because the variant peptides cannot be quantified simultaneously. Polymorphisms can result in reduced intensity or non-detection events for peptide isoforms that do not match the reference mouse genome. This in turn can affect protein abundance estimation from peptides and can either obscure the signal of a true pQTL or create a false local pQTL. Therefore, we filtered out polymorphic peptides based on the genome sequences of the founder strains. We further confirmed the presence of the expected polymorphisms by examining the distribution of peptide intensities across samples from the founder strains.
To determine whether peptides with polymorphisms matched their expected allele distribution pattern, the peptide data was standardized within batches and adjusted for batch effects. Each peptide was scaled by a sample-specific within-batch scaling factor: Batch effects were removed from the processed peptide data using a linear mixed effect model (LMM) fit with the lme4 R package. 64 Peptides unobserved for all samples within a batch were recorded as missing (NA). If greater than 80% of samples were missing for a polymorphic peptide, it was removed from the batch correction step and the subsequent evaluation. The following model was fit to peptide intensity data for the CC mice: where m is the intercept, b covar are the fixed effects of covariates, x T i; covar is the i th row of the covariate design matrix, u strain½i is the effect of the strain of sample i, u b½i is the effect the batch of sample i, and ε i is the error for sample i with ε i $ Nð0; s 2 Þ. The strain and batch effects were estimated as random effects: u strain $ Nð0; It 2 Strain Þ and u b $ Nð0; It 2 b Þ. For the CC and founder strains, sex was included as a covariate. A similar model was fit for the DO mice but with no strain effect and diet was included as a covariate along with sex. The batch effects, estimated as best linear unbiased predictors (BLUPs) using restricted maximum likelihood estimates (REML), 77 were subtracted from each peptide measurement:ỹ where local i is the effect of the local haplotype on peptide k for sample i, u kinship i represents a random kinship effect to account for overall genetic relatedness, and all other terms as previously defined. For the CC mice, local½i = p T i b local where p T i is the founder haplotype probability vector at the marker closest to the gene TSS (e.g., ordering the founder strains as AJ, B6, 129, NOD, NZO, CAST, PWK, and WSB, p T i = ½0 1 0 0 0 0 0 0 for a CC mouse i that is B6/B6 at the locus). For the DO, is the founder haplotype dosage vector, scaled to sum to zero, at the marker closest to the gene TSS (e.g., d T i = ½0:5 0:5 0 0 0 0 0 0 for a DO mouse i that is AJ/B6 at the locus). For the founder strains, local i = x T i; strain b local where x T i; strain is the founder strain incidence vector for mouse i (e.g., x T i;strain = ½0 1 0 0 0 0 0 0 for a B6 mouse). b local is an eight-element vector of founder haplotype effects, fit as a random effect: b local $ Nð0; It 2 local Þwhere I is an 838 identity matrix and t 2 local is the variance component underlying the local effects. The kinship effect is included for the CC and DO mice and modeled as u kinship $ Nð0; Gt 2 G Þwhere G is a realized genomic relationship matrix and t 2 G is the variance component underlying the kinship effect, accounting for population structure. [78][79][80][81] Here we used a leave-one-chromosome-out (LOCO) G, in which markers from the chromosome the peptide is predicted to be located on are excluded from G estimation in order to avoid the kinship term absorbing some of local½i. 82 We then calculated r poly = corð b b local ; qÞ, the Pearson correlation coefficient between b b local , the BLUP of b local and q, the incidence vector of the B6 haplotype among the founder strains (e.g., q = ½01000000 for a peptide that contains a B6-specific allele that is missing in the other founder strains). Sets of peptides with polymorphisms were defined based on having r poly >0:5 for each of the CC, DO, and founder strains, to be excluded from further analysis because they would bias protein abundance estimation. where M is the set of peptides that map to protein j, 1 i;m is the indicator function that peptide m was observed in mouse i, and q i is the scaling factor previously defined. 73  For the DO and founder strains, there was no bridge sample, and proteins were instead normalized as:ỹ prot j i = log 2 ðy prot j i + 1Þ. Batch effects were removed from the protein data using the LMM described for the peptide data (Equation 1). If more than 50% of samples were missing a protein, it was removed from further analysis in order to avoid false downstream findings. Batch effects, estimated as BLUPs, were then removed:ỹ

Heritability estimation
We estimated heritability for all proteins in the CC, DO, and founder strains. The heritability model is similar to Equation 2, but for proteins instead of peptides and without the local½i term: where terms are as previously defined. The genomic relationship matrix G -corresponding to the kinship term u kinship $ Nð0; Gt 2 G Þ for the CC and DO -is estimated from all markers, i.e., non-LOCO G -because there are no other genetic factors in the model. In the founder strains, G = X strain X T strain where X strain is the founder strain incidence matrix. Sex was modeled as a covariate for all three populations, and diet as well in the DO. Heritability is then calculated as h 2 = t 2 G t 2 G + s 2 . The estimate in the DO is for the narrow sense heritability, representing the contributions of additive genetic effects. For the CC and founder strain mice, the estimate represents broad sense heritability, incorporating non-additive genetic effects, due to the presence of replicates.

QTL analysis
In the CC and DO mice, we performed a genome-wide pQTL scan for each protein, testing a QTL effect at positions across the genome, using a model similar to Equation 2: wherez prot j i is the standard normal quantile returned by the inverse cumulative distribution function of the normal distribution on the uniform percentiles defined by the ranks ofỹ prot j , i.e., the rank-based inverse normal transformation (RINT) 83  For the CC mice, we mapped pQTLs based on strain averages wherez prot j i is the average ofỹ prot j male; strain i andỹ prot j female; strain i followed by RINT across the strains. Founder haplotype probabilities were reconstructed at the level of individual mice and averaged for strainlevel mapping. No covariates were included when mapping on strain averages. We tried mapping pQTLs in the CC mice on individuallevel data, which returned largely consistent results, but notably fewer and weaker pQTLs. In the CC, we also mapped pQTLs to the mitochondrial genome and Y chromosome by testing whether the founder origin of the mitochondria or Y chromosome was associated with protein abundance. We fit Equation 4, treating the mitochondrial genome or Y chromosome as a single locus QTL Y ½i and QTL MT ½i, respectively, using the non-LOCO G for the kinship effect. The founder strain of origin for the Y chromosome was determined for all CC strains. For the mitochondrial genome, six strains (CC031, CC032, CC041, CC051, CC059, CC072) possessed ambiguity between AJ and NOD, which we encoded as equal probabilities ðp T i;MT = ½0:5 0 0 0:5 0 0 0 0Þ.

QTL significance thresholds
We estimated significance thresholds for pQTLs using permutations. 42 We accounted for missing data by performing 10,000 permutations of the normal quantiles for each level of observed missingness in the CC and DO mice (ranging from 0 to 50%). Genome scans of the permuted data used the model in Equation 4, excluding covariates and the kinship term. We first applied a genome-wide error rate correction across marker loci and then applied an FDR correction to account for testing multiple proteins. 43 We modeled the maximum (genome-wide) LOD scores from the permutation scans using a generalized extreme value distributions (GEV) 84,85 specific to each level of missingness, to compute genome-wide permutation p-value for each protein: where F GEV; nNA½prot j is the cumulative density function for the GEV fit from the permutations of quantiles with n NA number missing values, corresponding to the number missing for protein j, and max LOD½prot j is the maximum LOD score from the genome scan of protein j. We then used the Benjamini-Hochberg (BH) procedure 86 to calculate FDR q-values across the permutation p-values, and applied interpolation to find the permutation p-value that corresponds to FDR < a: p interp perm; a where a˛[0.1, 0.5]. Significance thresholds on the LOD scale, specific to FDR < a and n NA missing data points, were calculated: l nNA FDR < a = F À1 GEV; nNA ð1 À p interp perm; FDR < a Þwhere F À1 GEV; nNA is the inverse cumulative density function for the GEV with n NA missing data points. As a final step to reduce random variation between sets of permutations, we regressed the estimated thresholds for a population and FDR level on the number of missing data points n NA , and created a table of fitted thresholds: b l nNA FDR < a for a˛[0.1, 0.5] for both the CC and DO mice. Whether a pQTL met FDR < a significance, the threshold corresponding to a with the n NA for protein j was used. Defining local/distal status of QTL Detected pQTLs were classified as local if their position was within 10 Mbp upstream or downstream of the middle of the coding gene. If they did not fall within this local window, they were classified as distal. The broad local window was used because the CC have larger LD blocks than the DO due to fewer outbreeding generations. With a narrower definition, it would be more likely to have ''distal'' pQTL in the CC that align and have consistent effects with ''local'' pQTL in the DO. On the other hand, this lenient definition of local may absorb some distally acting pQTLs that happen to be within 10 Mbp the gene on which they act.

Consistency of QTL between the CC and DO
We evaluated the consistency of local and distal pQTLs between the CC and DO by comparing their haplotype effects. We first had to define pQTLs that were detected in both the CC and DO and thus pair them for effect comparison. Local pQTLs were paired based on simply having matching protein IDs. For distal pQTLs, we also required the pQTL positions to be within 10 Mbp of each other.
Haplotype effects were estimated at the pQTL marker using the model in Equation 4. To stabilize the effects, they were modeled as a random effect: b QTL $ Nð0; Haplotype effects for a pQTL are fit at a specific marker. Selecting which marker for effect comparison is complicated by the fact that the CC and DO have different sets of markers and the genomic coordinates of the peak LOD scores also vary. When comparing pQTLs detected in both populations, we fit the Equation 4 model at the markers with the highest LOD score specific to each population. When comparing pQTLs that were detected in only one population, we selected the marker in the population that failed to map the pQTL that was closest to the marker in the population that detected it.
Consistency of local QTL in the CC with the founder strains If the genetic effects on a protein are primarily local, the relative abundances for a protein in the founder strains should match the local pQTL effects observed in the CC and DO. We evaluated the consistency of local pQTLs in the CC with the founder strains, using an approach similar to how we compared pQTL effects between the CC and DO. For the founder strains, rather than fitting pQTL effects ðb QTL Þ, we fit the founder effects as random terms (as described for the local term in Equation 2 for the founder strains) summarized as BLUPs ðb Þ. As when comparing QTL effects between the CC and DO, we then tested r local > 0, and corrected for multiple testing through the BH procedure.

Mediation analysis
For each distal pQTL (lenient threshold) in the CC or DO populations, we performed a mediation analysis which involved a scan analogous to the QTL genome scans. Instead of scanning through genetic markers as putative QTLs, we scan through proteins as putative mediators of a given distal pQTL. The model is where QTL½iis as defined for QTL m ½iin Equation 4 but fixed at the peak marker m of the distal pQTL for target protein t and conditioned on protein q, with mediator q ½i representing its effect on protein t for individual i, and all other terms as previously defined. The effect of the mediator is modeled as mediator q ½i = b prot qz prot q i , where b prot q is the regression coefficient for the mediator protein q andz prot q i is the RINT quantity of protein q for individual i. The likelihood of Equation 6 model is compared to a null QTL model that excludes the QTL i term, producing a mediation conditional LOD score. The mediation model is fit for all proteins as individual mediators, excluding protein t, resulting in a mediation scan.
We assume that most of the proteins evaluated as candidates are not true mediators of the pQTL and thus the distribution of mediation conditional LOD scores approximates a null distribution, roughly centered around the LOD score of the distal pQTL that was first detected. We calculate the z-scores of the mediation conditional LOD scores and then define strong candidate mediators of the pQTL for protein t as proteins with z med q < À4, where z med q is the z-score of the mediation LOD score for candidate mediator protein q. The rationale being that when testing the QTL term in Equation 6, if the mediator contains much of the information from the pQTL, its presence in both the alternative and null models will result in a large drop in the LOD score of the detected pQTL. For a protein to be declared as a candidate mediator of the distal pQTL, we required that the mediator TSS be within 10 Mbp of the pQTL marker. Strong mediators that were not near the pQTL often represent proteins that are correlated with the target protein t, which are often co-regulated members of a protein complex or pathway.
Sex effects on protein abundance analysis Proteins that exhibited differential abundance between the sexes, i.e., sex effects, were identified using an LMM similar to the heritability model (Equation 3) for the CC, DO, and founder strains, but instead testing the significance of the sex coefficient: where b Male is the effect on protein j of being male, x i; Male is an indicator variable of being male, and all other terms as defined previously. Other covariates and the specification of u kinship i for the different populations are the same as described for heritability.
A p-value for the sex effect was calculated by comparing the model in Equation 7 to a null model without the sex effect through the likelihood ratio test (LRT): p prot j sex = PrðX > b c 2 prot j Þ where Prð:Þ denotes the c 2 ð1Þ probability density function and b c 2 prot j is the observed LRT statistic for protein j. The LMM was fit with the qtl2 R package, 65 using maximum likelihood estimates (MLE) for parameters rather than REML, which are more appropriate for asymptotic-based significance testing of fixed effects. Proteins with significant sex effects were selected based on FDR < 0.1 using the BH procedure. 86 We performed gene set enrichment analysis using the clusterProfiler R package. 63 We defined gene sets based on q sex < 0.01 and split them further into subsets based on having higher abundance in males or higher abundance in females. We used the quantified proteins in each population as the background gene set. Hypergeometric tests for enrichment of GO and KEGG terms were performed with FDR multiple testing control. 87 Enriched GO and KEGG terms were selected based on having q set < 0:1.

Protein complex analysis
We assigned proteins to protein complexes using annotations. 49 For each protein complex, we quantified how tightly co-abundant, i.e., cohesive, the members are, by calculating the median pairwise Pearson correlation for each protein with the other members of the complex. We summarized cohesiveness within a complex by recording the median and interquartile range across the median correlations for the individual proteins.
To assess whether genetic factors or sex regulated protein complexes, we estimated the complex heritability and complex sex effect size based on the PC1 from PCA 51 of the abundances of the proteins annotated to the complex. We first filtered out proteins with local pQTLs (FDR < 0.5) or strong distal pQTLs (FDR < 0.1) to minimize the influence of proteins with independent genetic effects in order to focus on the shared effects on a protein complex. We also regressed out effects of covariates from the individual proteins prior to PCA in order to keep the PC1 summary from reflecting their effects. To estimate complex heritability, we removed the effect of sex in the CC, and both sex and diet in the DO. For complex sex effect size, we removed the effect of diet from the DO. We estimated complex heritability using the model in Equation 3, with no covariates and the complex PC1 as the response variable.
To estimate the complex sex effect size: f 2 sex = 1 À terval estimates for complex heritability and complex sex effects represent 95% subsample intervals. We randomly sampled without replacement 80% of the CC and DO data 1,000 times and estimated the complex heritability and complex sex effects for each subsample as well as the 2.5 th and 97.5 th quantiles across the subsamples. We estimated summaries for protein complexes that had four or more proteins observed in the CC or DO, after removing proteins with local pQTLs (FDR < 0.5) or distal pQTLs (FDR < 0.1), thus limiting the potential that the PC1 reflected a strong pQTL not shared by other members of the complex.

Strain-specific outlier proteins
To identify proteins with low or high abundance characteristic to individual CC strains, we fit the following LMM: with all terms as previously defined. Effects for all CC strains for each protein j ð b u prot j strain Þ were estimated as BLUPs, which were then transformed to z-scores per protein ðz prot j strain Þ. We defined a strain-specific protein outlier to be a protein j in CC strain i for which z prot j strain i >2:5. This represents a lenient threshold because we aim to cast a wide net and identify interesting characteristics of CC strains, potentially due to subtle effects across many proteins. We intersected the strain outliers with known CC strain-specific genetic variants based on CC strain identity and the annotated coding gene, 30 identifying variants that likely have local effects on protein abundance. For each CC strain i, we defined sets of proteins that had consistently low, high, and extreme (low or high) abundance based on their strain effects: U high strain i = fprot j : z prot j strain i >2:5g c j, U low strain i = fprot j : z prot j strain i < À 2:5g c j, and U extreme strain i = fprot j : z prot j strain i