Identification of Quantitative Trait Loci Underlying Proteome Variation in Human Lymphoblastoid Cells^*

2D DIGE Materials

Sodium chloride, Trizma (Tris base), bromphenol blue, chloroform, N,N-dimethylformamide, phosphatase inhibitor, sodium carbonate, sodium bicarbonate, ammonium bicarbonate, ammonium monobasic phosphate, l-lysine, and α-cyano-4-hydroxycinnamic acid were obtained from Sigma-Aldrich. Urea, thiourea, SDS, DTT, Tris, and iodoacetamide were obtained from GE Healthcare. Protein labeling dyes (CyDye DIGE fluors) with an N-hydroxysuccinimide ester reactive group for reaction with the ε amino group of lysine of proteins were obtained from GE Healthcare. CHAPS was obtained from USB Corp. Complete Mini protease inhibitor mixture tablets were obtained from Roche Applied Science. Methanol, acetonitrile, and HPLC grade water were obtained from Burdick & Jackson. Acetic acid was obtained from Mallinckrodt Baker. CEPH-CEU LCLs were acquired from the Coriell Institute for Medical Research.

Cell Cultures

Twenty-four cell lines from the CEPH-CEU collection were obtained from the Coriell repository. Identification numbers for the 24 cell lines are GM12057, GM07345, GM12145, GM10860, GM11829, GM12056, GM11840, GM11830, GM12004, GM12144, GM10846, GM07357, GM11839, GM12003, GM06994, GM11993, GM10856, GM11992, GM07000, GM07348, GM10854, GM07029, GM10838, and GM10851. Cells were cultured at 37 °C in RPMI 1640 medium supplemented with 15% fetal bovine serum (HyClone), 1 mm l-glutamine (Invitrogen), and 100 units each of penicillin and streptomycin (Invitrogen) in a humidified chamber containing 5% CO_2. Cells were plated at an initial density of 200,000 cells/ml and grown to between 800,000 and 1,100,000 cells/ml at which time they were harvested by washing two times in phosphate-buffered saline. Frozen pellets containing ∼10⁷ cells were stored at −80 °C.

Sample Preparation and Protein Labeling

Upon thawing, cells were lysed with 180 μl of TNE buffer (50 mm Tris-HCl, pH 7.6, 150 mm NaCl, 2 mm EDTA, pH 8.0, 1% Nonidet P-40) with protease and phosphatase inhibitors added. Soluble protein was extracted by centrifugation at 14,000 rpm for 10 min at 4 °C. Protein concentration was determined using the BCA Protein Assay kit (Pierce). Protein disulfide bonds were reduced in 2 mm DTT. For each replicate, 25 μg of each sample was precipitated using 4:1 (v/v) methanol:chloroform. The samples were washed once with methanol and resuspended in L-buffer (7 m urea, 2 m thiourea, 4% CHAPS, 30 mm Tris, pH 8.5). An internal standard was created by pooling 8.5 μg of each sample. The remainder of each sample was then labeled with 3 pmol of either Cy3 or Cy5 dye/ μg of protein. The internal pooled standard was labeled using 3 pmol/μg Cy2. Reactions were stopped by adding 10 mm lysine in equal volume to the volume of dye. For each gel, a Cy3- and Cy5-labeled protein sample was mixed with an equal volume of Cy2-labeled standard. Samples were then reduced and denatured in R-buffer (7 m urea, 2 m thiourea, 4% CHAPS, 13 mm DTT) at 10:1 (v/v) with ampholytes (IPG Buffer) (GE Healthcare). The samples were then applied to a 24-cm pH 4–7 Immobiline DryStrip (GE Healthcare).

2D DIGE and Image Analysis

First dimension isoelectric focusing was performed using an Ettan IPGphor3 (GE Healthcare) for a total focusing time of 64.5 kV-h at 20 °C and 75 μA/strip maximum. After focusing was complete, each strip was equilibrated with 10 ml of equilibration buffer (6 m urea, 75 mm Tris, pH 8.8, 29.3% glycerol, 2% SDS, 0.002% bromphenol blue) with 10 mg/ml DTT added for 15 min. The strips were then equilibrated with 10 ml of equilibration buffer with 25 mg/ml iodoacetamide for 15 min. The strips were then loaded onto 12% homogenous SDS-polyacrylamide gels (NextGen Sciences). Electrophoresis was performed on an Ettan DALTtwelve System Separation Unit (Amersham Biosciences) at 0.5 watt/gel until the bromphenol blue dye front had just run off the bottom of the gel. Once complete, the gels were scanned at 488, 532, and 633 nm using a Typhoon Trio Variable Mode Imager (GE Healthcare) and analyzed using the DeCyder 2D v6.5 software (GE Healthcare).

Protein Identifications

A pool of 10.5 μg of protein from each sample was used to generate a preparative 2D gel for spot excision and protein identification. Visualization of proteins was performed by poststaining using a 1:200 dilution of Deep Purple Total Protein Stain (GE Healthcare) for 1 h. The image was then imported into the DIGE analysis using the DeCyder software and matched to the prior analytical gels. Gel plugs corresponding to 492 proteins were isolated using an Ettan Spot Picker (GE Healthcare). Dehydrated gel spots were digested using Trypsin Gold (Promega) by overnight incubation. The digests were then desalted and concentrated using C₁₈ P10 ZipTips (Millipore). One-third of each sample was loaded onto the ZipTip by pipetting up and down five times into a fresh 96-well plate. The sample was then washed four times with 10 μl of 0.1% TFA and eluted with 2 μl of 50% acetonitrile, 0.1% trifluoroacetic acid. A total of 1 μl of eluate was pipetted onto a clean MALDI plate and covered with 1 μl of α-cyano-4-hydroxycinnamic acid MALDI matrix. Mass spectra of each spot were acquired using an Applied Biosystems 4700 Proteomics Analyzer MALDI-TOF-TOF instrument. Data were acquired in reflectron positive mode for the mass range 900–3500 Da. Shots per spectrum for both MS and MS/MS totaled 2000. Data-dependent MS/MS analysis was performed on the top 10 peaks from each MS spectrum. Peaks excluded from analysis correlated with trypsin autolysis peaks at 844.0, 1047.5, and 2257.5 Da. Peak lists were generated using the 4000 Series Explorer v3.0. The data were searched using Mascot v1.9.0 (Matrix Sciences) against the human NCBI protein RefSeq database build 36.3 consisting of 37,742 sequences. Search parameters included trypsin enzyme specificity with two permitted miscleavages. Variable modifications included carbamidomethyl-Cys and oxidized Met. Mass tolerances for both precursor and fragment ions were 1 Da. Identification criteria included a Mascot score >65 (selected based on a corrected p value <0.01), at least 5 orders of magnitude difference in the expect score between non-isoformic first and second hits, at least two unique MS/MS peptides scoring at least 20 (p < 0.01), and the protein having a molecular weight similar to experimentally derived values. Where the second hit was an isoform of the first and had identical scores, the first is reported. Where the second hit was an isoform with a lesser score, the higher scoring isoform is reported. Where a peptide matched to multiple proteins in the database, the peptide was assigned to the higher scoring protein. All peptide mass matches for an accepted protein are reported in supplemental Table 2 regardless of whether the matching peptide from the MS/MS search achieved the score cutoff.

Association Analysis

Protein expression levels from 544 spots derived from the average between two technical replicates were used as protein expression phenotypes for genetic association analysis. Genotype information from 3.8 million SNPs was downloaded from the International HapMap Project web site. Because of our small sample size, we removed SNPs with minor allele frequency <0.1 and genotype missingness >0.1 prior to association. This filtering resulted in a set of 1.7 million SNPs that were used in the final association analysis. For the cis-only analysis, SNPs were limited to the 200-kb region surrounding the transcript-coding allele. For each protein, expression levels were log-transformed and treated as the response variable. Based on the allele composition of the genotypes, each SNP (coded as 0, 1, and 2) constituted the single dependent variable, and a linear regression was performed using the genetic toolkit PLINK (44). Statistical significance of the strength of the relationship between the response variable and the dependent (genotype) variable was measured using the Wald test. Multiple testing corrections were made using the Benjamini and Hochberg false discovery rate (BH-FDR) procedure (45). Results were visualized and annotated using WGAViewer (46). Power calculations were performed using the genetic power calculator (47). We used the following settings: (a) locus-specific heritability ranging from 0.03 to 0.8, (b) sibling correlation = 0.5, (c) linkage disequilibrium (LD) (D′) = 1, (d) QTL increaser allele frequency = 0.3, (e) marker allele frequency = 0.3, (f) sample size = 16, sibship size = singletons (both parents genotyped), and (g) type I error rate (α) = 0.0001 and 0.0000001. Higher sibling correlation, a measure of total heritability of the traits, increases the power. Hence, we considered a sibling correlation of 0.5 as a reasonable setting: not too high, not too low, and typical for many traits. Locus-specific heritability values were considered from 0.03 (very small) to 0.8 (large). Because we considered 1.7 million SNPs, there was a high chance that the markers and the pQTL were near to each other; therefore, we set LD (D′) to 1. LD is the non-random association of alleles at two or more loci. D′ is a standardized measure of the degree of LD with a maximal value of 1. Supplemental Table 5 summarizes the power for cis and genome-wide analyses. We considered an α of 1.00e−07 for genome-wide analysis because it involved 1.7 million SNPs per trait. The cis-only analysis involved far fewer SNPs per trait; therefore, power calculations considered an α of 0.0001 for cis-only analysis.

Lymphoblast Proteome

An overview of the experimental approach is shown in Fig. 1. A 2D DIGE analysis was performed to determine base-line proteome expression levels in LCLs derived from 24 individual cell lines derived from 16 unrelated individuals and their eight children. Gel spots were numerically identified and matched between gels using the internal standard. Gel images were manually examined to determine the quality of the matching. Incorrectly matched spots were eliminated or corrected if possible. Spots were then eliminated from analysis if they appeared to be artifacts from streaking or were affected by streaking of another spot. Spots were also eliminated if not detected in more than one gel. As a result of these quality control measures, a total of 544 protein spots were subsequently analyzed (Fig. 2).

A separate preparative gel of the pool of all 24 samples was run and poststained with Deep Purple. On this gel, 492 of the 544 proteins appeared with high spot resolution and abundance. Gel plugs for these 492 proteins were isolated and subjected to an in-gel trypsin digestion followed by MALDI-MS/MS analysis. Identification of proteins was determined by combined peptide mass fingerprinting and MS/MS fragmentation pattern using Mascot. Using a Mascot cutoff score of 65 (false positive rate <1%), proteins within 311 spots were identified. These 311 spots corresponded to 243 different proteins. Predicted and actual molecular weights and pI values for proteins were compared for verification. All protein identifications are presented in supplemental Table 1. Protein and individual peptide IDs are presented in supplemental Table 2.

Proteome Transcription Correlation Analysis

The correlation of basal expression levels of the proteins with corresponding transcripts was analyzed using the publicly available raw transcript data as reported by Morley et al. (20) and acquired from the Gene Expression Omnibus accession number GSE2552. Pearson's correlation coefficient was calculated to compare protein expression values for the 16 individual parents with the gene expression values for the same 16 individuals. Each expression value was derived from the average of two replicates. Of the 311 identified spots, 251 could be directly compared with transcript values from the published microarray analysis. Of these, 41 showed a significant and highly positive correlation with transcript data (r > 0.5 and p < 0.05), and only two showed a significant negative correlation (r < −0.5 and p < 0.05). The remaining proteins did not correlate significantly with transcript levels (supplemental Table 3).

Sources of Variation

To address a concern that lower intensity spots would exhibit greater technical variability and therefore reduce the power to determine sources of biological variation among lower abundance proteins, a linear regression was performed to determine the role of spot volume in population variance. Averaged volume of the same gel spot detected across the population of 24 individuals was used as the predictor variable. Population variance for each spot calculated from the standardized and normalized log expression abundance values across the 24 individuals was used as the dependent variable. As shown in Fig. 3A, following normalization and standardization, spot volume accounts for only 1% of the experimental differences (r² = 0.01).

Another variable that could contribute to protein expression variation is the cell population doubling time (PDT). As reported by Choy et al. (48), many transcript phenotypes in these same lymphoblast cell lines were linked to the growth rate and ATP levels of the individual cell lines. Such an observation confounds the ability to find true genetic effects. To determine the PDT, cell lines were plated at 2⁵ cells/ml and were counted daily until cells reached 10⁶ cells/ml. PDTs ranged from 0.21 to 0.82 doublings/day across the 24 cell lines. A regression analysis was performed using expression values for each protein as the response variable and the PDT of each cell line as the predictor variable. The expression levels of 53 proteins were significantly related to PDT (p < 0.05, r² > 0.5). The 53 proteins represent ∼10% of the total number of proteins analyzed, slightly more than that reported for transcripts by Choy et al. (48). Protein expression was evenly split between positive PDT associations and negative associations (supplemental Table 4). Two examples of PDT-associated protein expression patterns are shown in Fig. 3, B and C. Among these 53 spots, we were able to identify 26 proteins. Functional annotation of these proteins revealed enrichment for nucleotide-interacting proteins and proteins with ATPase activity (p = 0.003; not shown).

cis-Acting Association Tests for peQTL

Traditionally, an eQTL is defined as a mappable locus that affects the expression level of a particular transcript. However, a protein appearing on a 2D gel at a location other than expected due to a splice variation or non-synonymous coding polymorphism that alters the overall charge of the protein or because it is differentially modified by phosphorylation or glycosylation represents important functional consequences of polymorphisms not detected if only considering protein expression levels alone. Therefore, we define a peQTL as a mappable locus that affects the level, molecular weight, charge, or modification of a protein. We first performed an association study of each of the 243 identified proteins and SNPs that resided within 200 kb of the genetic coding region for each protein (cis-acting). For association analysis, only data for the 16 unrelated parents were used. Where multiple spots corresponded to a single protein, each spot was considered independently. Each spot trait (i.e. expression level) was regressed against the cis genotypes using PLINK (44). From this analysis, expression levels of 19 spots corresponding to 11 proteins were significantly associated with polymorphisms in cis with the protein-coding gene region (BH-FDR < 0.01). To compare the regression models developed from the association of expression values with individual SNPs within the parent population with that of parents and children considered together, we computed the Pearson's correlation coefficient (r²) for each model with the expectation that there should be little change in r² when the children are added to the model. This r² value also represents the locus specific heritability measure for each phenotype (49). A full list of the 19 proteins or protein forms, the expression of which is associated with cis-acting polymorphisms, is found in Table I.

Among the proteins associating with cis-acting SNPs was a six-spot cluster in the 2D DIGE gel identified as L-plastin (LCP1) (spots 675, 677, 687, 688, 691, and 693). Several polymorphisms in the region of the LCP1 gene were observed to be significantly associated with expression of each of these six spots (for an example of spot 691, see Fig. 4A). One of these (rs4941543) is a non-synonymous single nucleotide polymorphism (NS-SNP) that causes a Lys to Glu alteration at position 533 in the LCP1 protein. To verify genetic association with this SNP, we determined r² for the 16 parental cell lines in addition to the eight children. Relative concordance was observed between the results of the parents alone and parents plus children (Table I, last two columns). In Fig. 4B, the LCP1 protein abundance from all 24 individuals in the population for all six spots clearly shows differential allele-specific expression patterns based on the presence of the Lys-coding allele. This substitution is predicted to shift the pI of LCP1 protein from 5.3 to 5.2. Such a shift is evident on the 2D DIGE gels and segregates with Lys-coding versus Glu-coding genotypes selected from the population (Fig. 4C).

Another of the proteins associating with a cis-acting QTL was Copine-1 (CPNE1). Expression levels of the CPNE1 protein (spot 780) associated with a cis-acting peQTL on chromosome 20 (rs3746410, p value = 4.94e−09). This region continues for about 600 kb upstream through intron 1 of the CPNE1 transcript coding region (Fig. 4A) and coincides with the transcript eQTL identified for the CPNE1 transcript when microarray analysis of these same cell lines was performed by two independent groups (15, 20). To determine the exact overlap between the peQTL and the eQTL, we obtained the raw transcript expression data for CPNE1 available from the GeneVar database (15, 17) and repeated the cis association analysis using data from only the 16 parental lines used in our peQTL analysis. Fig. 5A shows the genomic region that is significantly associated with both protein (bottom panel) and transcript expression (top panel) of CPNE1 and indicates that the regions are identical. Correlation of both the protein expression and transcript expression with a representative SNP (rs3746410) is shown in Fig. 5, B and C.

Two spots (495 and 502) identified as hematopoietic cell-specific lyn substrate 1 (HCLS1) proteins are associated with SNPs located in cis to the HCLS1 gene. These spots were part of a four-spot cluster that included spots 517 and 490 (Fig. 6A). Spot 517 was also identified as HCLS1, and spot 490 remains technically unidentified. When expression values derived from spot 490 were used as a phenotype trait in association with cis-acting SNPs in the HCLS1 gene, the same set of SNPs was highly significant (e.g. rs1128163, p value = 9.39e−09), suggesting that the protein represented by spot 490 is likely to be HCLS1. Additional SNPs in the 200-kb region surrounding the HCLS1 gene were significantly associated with expression (Fig. 6B) and ranged from the 5′-upstream region to the 3′-UTR. The correlation of SNP rs1128163 to protein expression level across the full 24 individuals in our population for spots 490, 495, and 502 (Fig. 6, C–E) supported the model. The exception was spot 517 (Fig. 6F). We also analyzed raw transcript data for HCLS1 transcript from analyses with these same 16 parental cell lines from the GeneVar database and did not find SNPs within the HCLS1 locus that were associated with transcript variation (rs1128163, p value = 0.28).

GWA Tests for peQTL

A full GWA analysis was then performed for each of the 544 analyzed spots regardless of whether the corresponding protein was identified. Each individual protein trait (i.e. spot) in the 16 parental cell lines was tested by linear regression against a total of 1.7 million SNPs available from the International HapMap Project. We applied the BH-FDR correction for multiple testing, setting a cutoff of 0.01. Because of the small sample size and potential for false positive results due to random association, we also considered only peQTL that were identified with three or more SNPs within a 100-kb region. To locate the QTL associated with each trait, the −log(p value) for association of each SNP with each protein trait was plotted according to the position of the SNP in each chromosome. As with the cis-only analysis, to verify the model, a representative SNP identified through association analysis with the 16 parents was plotted against protein expression from all 24 individuals, and a new Pearson's correlation coefficient was calculated. This resulted in the identification of eight proteins with associating QTL meeting GWA significance (Table II).

Several of the cis-acting QTL for LCP1, CPNE1, and HCLS1 that were identified through the cis-targeted analysis also met statistical significance in the GWA study. The standardized abundance of three spots, 688, 691, and 693, that were identified as the LCP1 protein showed a statistically significant association with rs4941543 (as described above) in addition to a second trans-acting QTL. This second strongly associating QTL region for all three of these spots was an intragenic region of chromosome 3 (p values <10e−08). However, when protein abundance from all 24 individuals in the population was plotted against genotype for the most significant SNPs in the chromosome 3 QTL region, the correlation was less apparent than when only the 16 parents were used (supplemental Fig. 1). The cis-acting QTL identified for spots 780 (CPNE1) and 490 (HCLS1 inferred) were also significant in the GWA analysis with no additional trans-acting QTL identified.

Two putative trans-acting peQTL were identified only within the GWA analysis (Table II). The expression level of vinculin (spot 175) was significantly associated with a trans-acting QTL located within the exocyst complex component 4 (EXOC4) gene and directly upstream of LRGUK, a leucine-rich repeats and guanylate kinase-containing gene on chromosome 7 (rs12539993, p = 3.75e−10). The EXOC4 gene has been previously identified in a GWA analysis of arterial stiffness phenotypes within the Framingham Heart Study (51). However, the correlation between the genotype of a representative SNP (rs12539993) and vinculin protein expression with children included was weaker (parents r² = 0.94 versus all r² = 0.40; Table II).

The expression level of LCP1 protein in spot 675 was associated more strongly with a trans-acting allele on chromosome 14 than it was for the cis-acting allele identified in our cis-only results (compare the top spot 675 results in Tables I and II). At the chromosome 14 location, five SNPs located around or within the NCP2 gene were significantly associated (p = 1.7e−08) with expression levels of spot 675. However, as was the case with vinculin, the correlation between genotype and expression values for spot 675 with children included was weaker than when considering the parents alone (parents r² = 0.88 versus all r² = 0.37; Table II and supplemental Fig. 2).

Spots 1415, 729, and 741 were associated with QTL on chromosomes 10, 11, and 2, respectively, but the proteins in these spots remain unidentified (Table II). The QTL associated with the expression level of spot 1415 lies entirely within the coding region of the MPP7 gene encoding a membrane-associated guanylate kinase, and the QTL associated with the expression level of spot 729 lies entirely within the gene coding for GUCY1A2, a guanylate cyclase protein. The region on chromosome 2 that is associated with expression levels of spot 741 is intergenic with no nearby gene features. Correlations between a representative genotype and protein expression levels among these unknowns are shown in supplemental Fig. 3.