Natural Genetic Variation Differentially Affects the Proteome and Transcriptome in Caenorhabditis elegans *

Natural genetic variation is the raw material of evolution and influences disease development and progression. An important question is how this genetic variation translates into variation in protein abundance. To analyze the effects of the genetic background on gene and protein expression in the nematode Caenorhabditis elegans, we quantitatively compared the two genetically highly divergent wild-type strains N2 and CB4856. Gene expression was analyzed by microarray assays, and proteins were quantified using stable isotope labeling by amino acids in cell culture. Among all transcribed genes, we found 1,532 genes to be differentially transcribed between the two wild types. Of the total 3,238 quantified proteins, 129 proteins were significantly differentially expressed between N2 and CB4856. The differentially expressed proteins were enriched for genes that function in insulin-signaling and stress-response pathways, underlining strong divergence of these pathways in nematodes. The protein abundance of the two wild-type strains correlates more strongly than protein abundance versus transcript abundance within each wild type. Our findings indicate that in C. elegans only a fraction of the changes in protein abundance can be explained by the changes in mRNA abundance. These findings corroborate with the observations made across species.

Natural genetic variation in gene expression shapes the diversity in phenotypic traits and is the raw material for evolutionary processes (1). Variation in gene expression can be very extensive across individuals with different genotypes. The additive effects (narrow-sense heritability) of independent loci on gene expression variation can reach 35% in humans (2). The broad-sense heritable variation in gene expression has been estimated to be up to 70% in the nematode Caenorhabditis elegans (3,4) and up to 80% in yeast (5). This high heritability and the ability to construct genetically segregating populations facilitate mapping of gene expression regulation and subsequent detection of expression quantitative trait loci (eQTL) 1 (5)(6)(7)(8)(9)(10)(11). eQTLs are genomic regions containing a polymorphism associated with variation in transcript abundance between genotypes (12). eQTL analysis provides insight into the underlying genetic architecture of complex traits and is valuable for the identification of pathways and gene networks (10,(13)(14)(15). A key question is whether gene expression variation is translated into variation at the proteome level and whether it affects functionally relevant proteins. Genetic model species provide an ideal platform to explore the relationship between gene expression variation and variation at the proteome level due to their tractability. Although it is well established that there is a correlation between transcript and protein abundances, the relationship between natural variation in gene expression and variation in protein abundance is less well understood. The proteome provides information required to understand the functioning of cells (16), and combined with natural genetic variation, it allows for mapping of protein expression regulators (17). In yeast, gene expression variation in 354 genes between genetically different strains was found to lead to variation in protein abundance for a limited set of proteins reflecting mainly transcription-independent mechanisms (18). In mice, the natural variation of transcript abundance and protein levels among inbred strains correlates slightly (0.27) (19). Over half of the identified eQTLs in yeast contributed to changes in protein levels of regulated genes, but several protein-QTLs did not correlate with their cognate transcript levels (20).
C. elegans is an ideal model organism to analyze the effect of the genetic background on its transcriptome (8,9,(21)(22)(23) and proteome because comprehensive proteome catalogues have been generated (24,25). The quantitative proteome comparison in C. elegans and Drosophila melanogaster showed that the interspecies' protein abundance correlation was higher than the intraspecies' correlation between protein and mRNA abundance. This suggests that protein levels are under evolutionary selection acting post-transcriptionally (26). In a recent interspecies comparison of the two nematodes C. elegans and Caenorhabditis briggsae, the protein and transcript changes were found to be conserved throughout development (27).
To complete these comparative interspecies expression analyses with an intraspecies analysis and to investigate the effect of natural genetic variation from the transcriptome to the proteome, the genetically highly divergent C. elegans wild-type strains, N2 (isolated in Bristol, UK) and CB4856 (from Hawaii), were compared quantitatively. Transcriptome data were acquired by microarray analysis, and for quantitative proteome analysis, the SILAC labeling method was used, which was recently established for C. elegans (28,29). We analyzed animals at the developmental larval stage 4 (L4), because at this stage the largest differences in gene expression levels were detected between N2 and CB4856 (3, 4, 30 -32). The genomes of these two wild-type strains, as well as polymorphisms between them, have been characterized extensively (33)(34)(35)(36)(37)(38)(39)(40). These genetic polymorphisms often lead to gene expression differences, which are playing a role in a whole range of phenotypic differences between N2 and CB4856, such as life history traits (4,(41)(42)(43)(44), behavior (40,(45)(46)(47)(48)(49), and life span (9,41,50). Here, we asked whether natural genetic variation causes variation in protein levels between N2 and CB4856 and whether it is related to variation in transcript abundance.
Using SILAC on three independent biological replicates, we quantified a total of 3,238 proteins, 2,485 proteins were identified in at least two biological replicates. We found 129 proteins to be significantly differentially expressed between N2 and CB4856 with at least a 1.3-fold change in abundance and a p value Ͻ 0.000457. They were enriched for genes that function in insulin-signaling and stress-response pathways, underlining strong divergence of these pathways in nematodes. The protein abundance of the two wild-type strains correlated more strongly than the protein abundance versus transcript abundance within each wild type.
Transcriptome Comparison-For RNA isolation we used a Maxwellா 16 AS2000 instrument with a Maxwellா 16 LEV simplyRNA tissue kit (both from Promega Corp., Madison, WI). The mRNA isolation was preceded by a modified lysis step. In short, 200 l of homogenization buffer, 200 l of lysis buffer, and 10 l of a 20 mg/ml stock solution of proteinase K were added to each sample. The samples were then incubated for 10 min at 65°C, 1000 rpm in a Thermomixer (Eppendorf, Hamburg, Germany). After cooling on ice for 1 min, the samples were pipetted into the cartridges, and the protocol as recommended by Promega was continued. After mRNA isolation, the "Two-color Microarray-based Gene Expression Analysis, Low Input Quick Amp Labeling" protocol, version 6.0, was followed, starting from step 5 (4,30,31).
The microarrays used were C. elegans (V2) Gene Expression Microarray 4ϫ44K slides, manufactured by Agilent Technologies, Santa Clara, CA. mRNA isolation, labeling with cyanine-3 and cyanine-5, and hybridization were performed as recommended by Agilent. The microarrays were scanned using an Agilent High Resolution C Scanner, using the settings as recommended. Data were extracted with the Agilent Feature Extraction Software version 10.5, following the manufacturer's guidelines.
For processing the data of the RNA microarrays, the "Limma" package for the "R" environment was used. No background correction of the RNA-array data was performed as recommended by Ref. 51. For the "within-array normalization of the RNA-array data" the Loess method was used, and for the "between-array normalization" the Quantile method was used. The obtained normalized intensities were used for further analysis.
SILAC Labeling-E. coli AT713 strain (lysine/arginine auxotrophic; E. coli Genetic Stock Center, CGSC number 4529, Yale) was grown in M9 Minimal Salts Medium (30.0 g of Na 2 HPO 4 , 15.0 g of KH 2 PO 4 , 2.5 g of NaCl, 5.0 g of NH 4 Cl, H 2 O to 1 liter) supplemented with 150 mg/liter of either light (Arg-0, Lys-0) or heavy amino acids (Arg-10, Lys-8; Cambridge Isotope Laboratories). The cultures were kept for 2 days at 37°C on a Lab-Therm Kü hner orbital shaker at 230 rpm with a shaking diameter of 5 cm and harvested at an A 600 between 1.5 and 3.0. Bacteria were pelleted at 1,800 ϫ g for 15 min at 4°C; the supernatant was aspirated, and aliquots of 50 ml were frozen at Ϫ20°C.
Adult N2 and CB4856 worms were bleached, and ϳ20,000 larval stage 1 animals were transferred to NGM plates freshly seeded with light or heavy labeled bacteria at 20°C. For the first biological replicate, CB4856 animals were fed with heavy labeled AT713. For the second and third replicates, N2 animals were fed with heavy labeled AT713. Nematode populations were grown for two generations, and proteins were isolated from animals at L4. To determine the labeling efficiency, N2 heavy protein extracts were analyzed on an LTQ Orbitrap XL mass spectrometer (Thermo Scientific). The mgf files were searched against the C. elegans 6,239 database (07/03/2010, database (07/03/2010, Functional Genomics Center Zurich), using the Mascot software. The library was downloaded from the UniProt database, contains 24,362 entries, and was supplemented inhouse with 259 common MS contaminants. In total, 3,908 assigned peptide spectrum matches were analyzed with a score higher than 30, which yielded 1,562 peptides. Of these, 1,552 (99.6%) peptides were heavy labeled.
Protein Isolation, SDS-PAGE, Protein Digestion-To extract proteins, worm samples were homogenized with glass beads (G1277 acid-washed beads, diameter of 212-300 m, Sigma-Aldrich) in freshly prepared cell lysis buffer (8 M urea, 5% 1 M Tris-HCl, pH 8.3) in a ratio of 1:1:2 (worms/beads/buffer) at 4°C for 30 s at 5 m/s four times (FastPrepா-24, MP Biomedicals). The lysates were centrifuged three times at 20,000 ϫ g for 10 min at room temperature to remove debris. Protein concentration in the supernatant was determined using the Bradford reagent (Sigma-Aldrich), and 1 g of each protein sample was checked for equimolarity by SDS-PAGE.
The L4 animal lysates of N2 and CB4856 were combined in a 1:1 ratio. Protein disulfide bridges were reduced with 5 mM dithiothreitol (DTT) at 60°C for 30 min and alkylated with 15 mM iodoacetamide in the dark at 37°C for 1 h. 600 g of proteins were digested with trypsin (modified sequencing grade porcine, Promega) in a ratio of 1:50 w/w overnight at 37°C. Samples were kept frozen at Ϫ20°C until further processing.

LC-MS/MS Analysis-LC-MS/MS
was performed on a reversedphase nano-LC system (Eksigent) at pH 3. Peptides were separated on a self-packed reverse-phase column (75 m ϫ 10 cm) packed with C18 beads (Magic C18, AQ, 3 m, 200 Å, Bischoff GmbH, Leonberg, Germany) at a flow rate of 200 nl/min. The column was equilibrated with 95% solvent A (0.1% formic acid in water) and 5% solvent B (0.1% formic acid in ACN). Peptides were eluted using the following gradient: 0 -1 min; 5-9% B, 1-56 min; 9 -40% B, 56 -60 min; 40 -50% B and 60 -64 min; 50 -95% B. Peptides were analyzed on an LTQ Orbitrap XL mass spectrometer in the data-dependent acquisition mode. High accuracy mass spectra were acquired in the mass range of 300 -1,800 m/z. In parallel, up to six data-dependent MS/MS were recorded in the linear ion trap of the most intense ions with charge state 2ϩ, 3ϩ, and 4ϩ using collisioninduced dissociation. Target ions already selected for MS/MS were dynamically excluded for 60 s. Each sample was measured twice. The second analysis was performed with an exclusion list containing all precursor values of the first analysis with an elution time window of Ϯ2.5 min.
MaxQuant Analysis-Data acquired on the LTQ Orbitrap XL MS were analyzed with MaxQuant version 1.3.0.5 (Max Planck Institute of Biochemistry Munich (52), searching the C. elegans 6,239 database at the Functional Genomics Center Zurich, Zurich. Search parameters were as follows: cysteine carbamidomethylation as fixed modification, protein N-terminal acetylation and methionine oxidation as variable modifications; SILAC labeling (Arg-10, Lys-8) as heavy labels; enzyme trypsin; two missed cleavages were allowed; and a minimum of six amino acids per identified peptide were required. The precursor ion mass tolerance was set to 20 ppm, and the fragment mass tolerance was set to 0.5 Da. The peptide FDR was set to 1%; protein FDR was set to 5%. In total, 40 raw data files were generated for every biological replicate and combined for database searching with the match between runs set to 2 min. Proteins with all shared peptides were combined into one protein group. The "Majority Protein IDs" column was used for protein assignments. The annotated MS/MS spectra with the lowest posterior error probability for single peptide identifications were extracted from 2,485 protein dataset using R package protViz (53) and are shown in supplemental Figs. S3-S5. The lists of corresponding single peptides are included in supplemental Tables S6 -S8 for all three replicates.
Correction of Arginine-to-Proline Conversion-To estimate the effect of the arginine-to-proline conversion, all raw files were analyzed with the Progenesis QI for Proteomics software (Nonlinear Dynamics). The MS1 m/z feature maps filtered for the charge states from 2 to 7 were generated using the automatic method for peak picking with a sensitivity value of 5. The number of fragment ion counts was limited to 200, and deisotoping and charge deconvolution were applied. The generated mgf files were searched against the C. elegans database with Mascot for the identification of heavy proline peaks. The precursor mass tolerance was set to 10 ppm, and the fragment ion mass tolerance was set to 0.6 Da. Only fully tryptic termini with two missed cleavages were considered. The isotopic labeling of arginine (Arg-10) and lysine (Lys-8), heavy proline (Pro-6), heavy glutamate (Glu-6) labels and methionine oxidation were set as variable modifications. Cysteine carbamidomethylation was set as fixed modification. The analysis yielded 301 peptides with heavy proline and heavy arginine and/or heavy lysine. Based on this, the contribution of the heavy proline signal to the heavy arginine-lysine signal was estimated to 20% on average for peptides carrying one proline (supplemental Fig.  1A). This contribution to the signal intensity was proportional to the number of prolines per peptide (supplemental Fig. 1, A and B). Using a Perl script kindly provided by Jacob D. Jaffe (Broad Institute, Cambridge, MA), new protein group log2 heavy/light (H/L) ratios were calculated. In a first step, the intensity values of heavy peptides with prolines were corrected by adding to the measured signal intensity a value, which corresponds to 20% ϫ n of heavy peptides abundance, where n is the number of prolines in the peptide. In a second step, adjusted peptide log2 H/L ratios were built and combined into the protein groups based on the MaxQuant protein group identifications. For each protein group median log2 H/L ratios were calculated. The protein group log2 H/L median ratios were calculated separately for each biological replicate; median was normalized to 1 and used for further analysis (supplemental Table 2). The effect of the arginine-toglutamate conversion was negligible, because only 15 peptides identified had a heavy glutamate with heavy arginine and/or heavy lysine. Corrected peptide intensity values were used to calculate the corrected iBAQ values (intensity-based absolute quantification). The C. elegans database was digested in silico with trypsin using the Skyline software (version 1305, MacCoss Lab (54)), and the number of tryptic peptides per protein was extracted. The length of the peptides ranged from 6 to 30 amino acids; the N-terminal amino acids and potential ragged ends were not excluded. Using an in-house R function, corrected intensities of peptides per protein were summed up and divided by the corresponding number of tryptic peptides. Heavy and light iBAQ values were calculated separately and normalized using an intersect of 0 and slopes of 1.01896, 1.00456, and 1.00487 for the first, second, and third replicates, respectively. The proline corrected log2 iBAQ and protein abundance ratios correlated well with the uncorrected ratios (0.66 and 0.71, respectively), and the correlations for the differentially expressed proteins were even stronger (ϳ0.9, supplemental Fig. 1, C and D). Proline corrected heavy and light peptide intensities correlated well between the replicates (r Ն0.74, see supplemental Fig. 2).
Experimental Design and Statistical Rationale-mRNAs were quantified in three biological replicates. In total, six worm samples were analyzed, three for each strain (N2 and CB4856). Gene expression differences between the genotypes were determined by a linear model using "log2 hybridization intensities ϳ genotype ϩ error." Correction for multiple testing was done using a permutation determined threshold (3,4,31,55). For each probe the values were permuted once and analyzed by the same model. The FDR rate was set to a ratio of false to true positives of 0.013, which occurred at p Ͻ ϭ 0.0032 (Ϫlog10(p) ϭ 2.5). False positive was the number of genes with a significant p value from the permuted set, and true positives were the number of genes with significant values from the original set. Proteins were quantified in three biological replicates, including label switching. For every biological replicate, two biochemical and two technical replicates were analyzed. In total, six combined protein samples were measured. To identify proteins that are differentially expressed between N2 and CB4856, SILAC data were filtered for proteins that were quantified at least twice with at least a 1.3-fold difference in abundance (the biological significance of the 1.3-fold change cutoff has been shown by Ref. 56, and a p value Ͻ 0.000457 (z-test, permutation determined threshold with the FDR cutoff of 0.05)).
Enrichment Analysis-The enrichment analyses were undertaken using the following databases with annotations. Gene ontologies were downloaded from Ensemble (version 83) using the BioMart package in R. The enrichment of GO terms was analyzed using the topGO package from the bioconductor suite in R (58). GO terms with less than two annotated genes were excluded from the analysis. The enrichments were computed using the Fisher test. The elim algorithm was applied to account for the underlying GO graph topology (57). The anatomy terms, protein domains, and gene classes were obtained via WormMart of the WS220 WormBase release. Genes from WormBook chapters were obtained from the 2012 version of Worm-Book. eQTLs (8,9) were obtained from WormQTL (21)(22)(23). KEGG pathways were obtained from release 65.0 of the Kyoto Encyclopedia of Genes and Genomes. These were tested for the enrichments by a hypergeometric test using R. Specific gene sets were obtained from supplementary data of the following papers: the aging set (59); DAF-16 sets (60, 61); and the DAF-16/PQM-1 set (62). The enrichments for specific gene sets were performed using Benjamini-Hochberg multiple testing corrected two-sided Fisher exact (aging and DAF-16 datasets) and hypergeometric tests (DAF-16/PQM-1 dataset). Note that the hypergeometric test and one-sided Fisher exact test are basically the same tests (63).
Data Storage-The transcript profiles and protein levels were stored in WormQTL (21)(22)(23) and can be accessed on line.
The MS proteomics data have also been deposited to the Pro-teomeXchange Consortium (64) via the PRIDE partner repository with the dataset identifier PXD002010.

RESULTS
Transcriptome Comparison-As the strongest mRNA expression differences between N2 and CB4856 were observed at stage L4 (3, 4, 30 -32), we measured the genome-wide transcription levels of three biological replicates of L4 stage synchronized larvae using microarrays. We found 1,532 genes to be differentially expressed (7.4%, Ϫlog10(p) Ͼ 2.5; FDR ϭ 0.013) of which 712 showed higher expression in CB4856 and 820 higher expression levels in N2 (Fig. 1). Differentially expressed genes were enriched for groups of genes involved in binding of proteins, sugars, and DNA, like f-box, math/bath/ btb, clec, and nhr genes (supplemental Table 4) as observed previously (3,4,32). These groups of genes are highly polymorphic between wild isolates (4,33,34). No major differences in genes linked to L4 development (30) could be detected between CB4856 and N2. We identified specific enrichments between the genes expressed higher in N2 compared with CB4856 (supplemental Table 4). For example, the GO-term "monooxygenase activity" was specifically enriched for genes expressed higher in CB4856. Because many of the gene expression differences between CB4856 and N2 are caused by their genetic differences, most of these will also be present in recombinant offspring. The loci causal for gene expression differences can be observed as eQTLs, which yield information about the genetic and regulatory architecture of gene expression. As expected, the group of 1,532 genes differentially expressed between CB4856 and N2 was highly enriched for genes with eQTLs that were previously identified in these recombinant inbred line populations (6,8,9) (928/1,532; p Ͻ Ͻ 10 Ϫ100 (8), 431/1,532; p Ͻ Ͻ 10 Ϫ100 (9)). To investigate whether these mRNA expression level differences are also found on the protein level, we measured protein abundance in both strains.
Protein Abundances Are Very Similar in N2 and CB4856-We quantitatively compared protein abundances between the two C. elegans wild-type strains N2 and CB4856 in the L4 stage using SILAC and mass spectrometry. Protein extracts were analyzed in three biological replicates, including label switching. Almost complete labeling (99.6%) of the second generation L4 worms was observed. In total, we quantified 3,238 distinct proteins, ranging from 2,536 to 2,622 proteins per SILAC experiment. 2,002 proteins (61.8%) were quantified in all three biological replicates and 2,485 proteins in at least two (protein FDR 5%, Fig. 2A, and supplemental Table 1). We compared our SILAC data to the integrated PaxDB dataset (65), which includes protein abundance data generated us-ing various fractionation methods. We quantified 24.5% of the C. elegans proteins presented in PaxDB (13,224) in N2 and CB4856, distributed over 6 orders of magnitude, with a clear bias for high abundance proteins (Fig. 2B). Protein abundance between the two strains correlated very strongly  0.99, Fig. 3C). To identify proteins that are differentially expressed between N2 and CB4856, the combined SILAC data (3,238 proteins) were filtered for proteins that were quantified at least twice (2,485 proteins) with at least a 1.3-fold difference in abundance and a p value Ͻ 0.000457 (z-test, FDR ϭ 0.05, Fig. 2C). Using these criteria, 129 proteins (5.2%) were found to be differentially expressed between the strains distributed over 5 orders of magnitude in expression. 70 proteins were up-regulated, and 59 proteins were down-regulated (Fig. 5).
Local and Distant Regulation of the Conserved mRNA and Protein Expression-In total, 2,481 quantified proteins (14 proteins with two isoforms each included) overlapped with 2,467 corresponding unique mRNAs in both strains (supplemental Table 3). The mRNA and protein expression patterns of N2 and CB4856 were very similar, whereas the correlation between mRNA and protein abundances within the strains was quite variable (r ϭ 0.54 for both, Fig. 3A). By contrast, mRNA abundance, protein abundance, as well as protein-to-mRNA ratios correlated strongly between the strains (r Ͼ 0.97, Fig. 3, B-D). This suggests that for the detected genes both transcriptional and translational regulation of gene expression are well conserved between N2 and CB4856. The conservation of protein expression was also observed between distantly related species (26).
We also compared the protein abundance differences versus mRNA abundance differences between the two strains ( Fig. 4). Of the 2,481 quantified genes, only 24 genes were found to be changed significantly at both the mRNA and protein levels between N2 and CB4856, and 105 proteins were different only at the protein level and 89 only at the mRNA level. To investigate the genetic loci possibly responsible for the mRNA and protein differences, we compared our data to eQTLs found in Vinuela et al. (9) and Rockman et al. (8). The vast majority of genes (16/23 and 23/25) found to be differentially expressed at both mRNA and protein levels had a QTL (cis-or trans) associated with them, with cis-eQTLs outnumbering trans-eQTLs with a ratio of about 2 to 1. By contrast, for the genes with no expression difference, eQTLs were relatively infrequent (ϳ7-21%), with a slight bias toward trans-eQTLs. Genes varying at the protein or mRNA level only had intermediate eQTL frequencies, with the number of trans-eQTLs being higher for the genes differentially expressed only on a protein level (Table I).
Aging and Stress Response-related Proteins Are Enriched among the 129 Differentially Expressed Proteins-To characterize the function of the 129 differentially expressed proteins, we first analyzed their GO annotations. In our dataset, the categories carbohydrate metabolisms, innate immune response, defense response to bacterium, and cysteine biosynthesis from serine were over-represented. We have also found the NPL4 domain and domains involved in sugar metabolism over-represented among other enriched multiple protein domains (supplemental Table 4). Interestingly, npl4 genes were shown to be involved in the life span of C. elegans (8,66). As altered sugar metabolism, oxidative stress, and pathogen resistance pathways were shown to be associated with aging (67)(68)(69)(70), and because the average life span of CB4856 is shorter than that of N2 (41,50,71), we compared our list with recently published aging and stress response transcriptome datasets. Transcripts differentially regulated during aging (p ϭ 0.007 (59) and transcripts responsive to the DAF-16 transcription factor (p ϭ 0.0002 (60,61)) were over-represented (Fig. 5  and supplemental Table 5). 2,485 proteins were used as the background list, and Benjamini-Hochberg multiple testing correction was applied. We further compared our data to a recent study that showed that the DAF-16-associated element-binding transcription factor PQM-1 complements DAF-16 aging-related functions (62). Genes transcriptionally regulated by PQM-1 alone were enriched in the 2,485 quantified proteins (p Ͻ 2.15e Ϫ13 , data not shown), whereas DAF-16 and PQM-1 targets were not (p ϭ 0.76, data not shown). By contrast, among the 129 differentially expressed proteins, enrichment was found only for the DAF-16 and PQM-1 targets (p ϭ 0.0008), whereas the PQM-1-specific targets were not enriched (p ϭ 0.24, Fig. 5 and supplemental Table 5). This shows that highly abundant proteins have a higher chance to be PQM-1 targets, and the variation in protein abundance acts downstream of DAF-16 and PQM-1. In total, 51 proteins were enriched in at least one of the mentioned stud- ies, and 18 proteins were identified in at least two of them ( Fig. 5 and supplemental Table 5). Together, these observations suggest that the 129 differentially expressed proteins might contribute to the differences in life span of CB4856 and N2.
The 129 differentially expressed proteins were also enriched for eQTLs identified in aging worms (9) (47 proteins, p Ͻ Ͻ 1e Ϫ16 ) and for stage-specific eQTLs (8) (62 proteins, p Ͻ Ͻ 1e Ϫ16 ) ( Table I). Eleven proteins were found to have an eQTL under mildly stressful conditions of 24°C (Y37A1B.5, Y48G1A.4, SSP-16, GRSP-3, MCE-1, MPPB-1, TATN-1, T06F4.1a, NDX-4, SPDS-1, and LEC-8) (6). Combined, these eQTL enrichments provide further evidence suggesting that at least part of the protein expression level differences can be explained by mRNA expression level differences between CB4856 and N2, and they might play a role in aging.  Genes with eQTLs on the npr-1 Locus Are Enriched among the 129 Proteins-CB4856 favors lower oxygen concentrations, mainly caused by an allelic difference in NPR-1 (39,40), and this is reflected in its clumping behavior (72,73). We found an enrichment for genes with an eQTL on the npr-1 locus in the 129 candidates, suggesting that at least part of the gene expression differences lead to differences in protein levels. Furthermore, two cysteine synthases, CYSL-2 and CYSL-3, proposed sensors of cellular H 2 S levels upon hypoxia (74), showed up-regulated protein levels in CB4856 (Fig.  5). The oxygen-binding protein GLB-1, the transcript of which is up-regulated under hypoxic conditions and in hif-1 mutants (75,76), was also up-regulated in CB4856 (Fig. 5), as well as the oxygen-sensing globin protein GLB-5 (72,73). These differences on protein level might contribute to the different oxygen responses of N2 and CB4856.

DISCUSSION
In this paper, we analyzed the influence of genetic background and natural variation on the proteome and the transcriptome in two genetically highly divergent C. elegans wildtype strains, N2 (Bristol) and CB4856 (Hawaii). Using SILAC, we quantified in total 3,238 distinct proteins in three biological replicates; 2,485 proteins were quantified in at least two replicates. These numbers are comparable with other SILAC studies in C. elegans. More than 1,400 proteins were reliably quantified in worm samples using size exclusion chromatography (28), and 3,470 proteins were quantified using LysC digestion followed by hydrophilic interaction liquid chromatography separation of the peptides (29). Between 1,461 and 4,072 proteins were identified in different developmental stages of C. elegans and C. briggsae after separation on SDS gels (27). To quantify less abundant proteins, additional separation methods would be required, as was done for the quantitative comparison of distantly related species (26). In our study, primarily highly expressed proteins were quantified. These show relatively little variation between N2 and CB4856, also at the transcript level.
Among the 2,485 quantified proteins, 129 (5.2%) proteins distributed over 5 orders of magnitude in expression were differentially expressed between N2 and CB4856. By contrast, of all the genes, 1,532 genes (7.4%) were differentially expressed at the mRNA level. Of the 2,481 genes for which proteins were detected and that overlapped with the mRNA dataset, only 113 (4.6%) were different at the mRNA level. The overall correlations for both protein and mRNA abundances were close to one between the strains. This is consistent with the quantitative comparison of C. elegans and C. briggsae (27), but it is different for C. elegans and D. melanogaster, where the conservation in expression was found only on protein level (26). Thus, in this intraspecies comparison, we observed a strong selective pressure to maintain protein levels, an effect that was also shown for closely (27) and even distantly related species (26).
The 129 differentially expressed proteins were strongly enriched for downstream targets of the insulin-signaling pathway (DAF-16 and PQM-1 targets) (60,61,62). An increased expression divergence of DAF-16 targets was also observed in a quantitative comparison of gene expression between the two nematode species C. briggsae and C. elegans (27). These results imply that differential regulation of the insulin-signaling downstream targets can be observed already between strains of the same species, underlining the strong divergence of this signaling pathway in nematodes. Importantly, the genes differentially expressed between N2 and CB4856 were also strongly enriched in eQTLs associated with aging. These results suggest that at least a subset of our differentially expressed genes might contribute to the longevity differences between N2 and CB4856.