Identification and MS-assisted interpretation of genetically influenced NMR signals in human plasma

Study population

The KORA study is an independent population-based survey from the general population living in the region of Augsburg, Germany. The KORA S4 study was conducted in 1999-2001 and comprises a total of 4,261 participants [14]. Between 2006 and 2008, a total of 3,080 subjects participated in a follow-up examination, KORA F4, which is the basis for the results presented here. All participants gave their signed informed consent and the local ethics committee approved the studies.

Blood sampling

Blood samples for metabolic analysis and DNA extraction were collected as part of the KORA F4 follow-up study. To avoid variation due to circadian rhythm, blood was drawn in the morning between 08:00 and 10:30 after overnight fasting. A part of the blood was drawn into serum gel tubes, gently inverted twice, and then allowed to rest for 30 min at room temperature (18-25°C) to obtain complete coagulation. The material was then centrifuged for 10 min (2,750 g at 15°C). Serum was divided into aliquots and kept for a maximum of 6 h at 4°C, after which it was frozen at -80°C until final analysis. Another part of the blood was drawn into ethylene diaminetetraacetic acid (EDTA) tubes, gently inverted twice, and left on the Sarstedt Universal mixer <5 min to avoid mechanical hemolysis, followed by centrifugation at 15°C for 10 min at 2,750 g. Thereafter, plasma was separated, divided into 200 mL aliquots and kept at 4°C, after which it was deep-frozen to -80°C. Within 2 weeks, plasma was stored in the gaseous phase of liquid nitrogen at -196°C.

Genotyping

For all individuals profiled from the KORA study, genome-wide single nucleotide polymorphism (SNP) data were already available. These data have been used and described extensively in the past in the context of several GWAS (for example, [4, 5]). Therefore, we summarize only the essential details here. For genotyping, 1,814 randomly selected participants of KORA F4 were included. These samples were genotyped using the Affymetrix Human SNP Array 6.0 (sample call rate >93%). Genotypes were determined using the Birdseed2 clustering algorithm. For quality assurance, the criteria of SNP call rate >95%, minor allele frequency >1%, and P(Hardy-Weinberg) >10^-6 were applied as filters. In total, 655,658 autosomal SNPs satisfied these criteria.

Metabolomics measurements

Metabolic analyses were conducted using clinical biochemistry methods, two distinct MS-based platforms (targeted and non-targeted), and NMR spectrometry (lipid classes and binned spectral data). For the joint analysis, metabolomics and genotype data were available for a total of 1,757 individuals. All metabolite measurements have been reported before. We therefore only summarize here the points that are essential for the present study.

Statistical analysis

To test for associations between genetic polymorphisms and individual NMR signal bins, we created age- and gender-adjusted linear additive models using PLINK (Version 1.07) [16]. Testing all possible ratios between NMR signals all over the spectrum would not be feasible due to the huge amount of computational time that would be needed. Therefore, we performed a simple feature selection. In most cases SNPs do not associate with only one NMR signal bin but with a number of adjacent bins. Thus, we performed a search for local minima on the pseudo-spectra (chemical shift vs. strength of association) resulting from the first GWAS run to pick the positions on the spectra which show the strongest associations. We then used the 500 best-scoring NMR signal positions to compute pair wise ratios (124,750 in total). As before, we used age- and gender-adjusted linear additive models in PLINK to test for associations.

We applied a conservative Bonferroni correction to control for false-positive error rates resulting from multiple testing. We corrected for tests on 655,658 SNPs and 133,350 NMR features at a nominal significance level of 5%, thus obtaining an adjusted P value of 0.05/655,658/133,350 = 5.72 × 10^-13.

Spearman correlations between biochemically and MS determined metabolite concentrations and NMR signal intensities were calculated using the statistical analysis system R (Version 2.15.1). To avoid false-positive associations due to small sample sizes, only metabolic traits with at least 300 non-missing values were included. Furthermore, metabolite concentrations that were more than three times the standard deviation away from the mean were excluded.

Seven genetic associations identified

We adopted a two-stage approach to detect associations between genetic variants and NMR signals. In the first step, we performed a GWAS using the NMR intensity readouts of the ¹H-NMR spectra from 1,757 plasma samples. To this end, we binned the spectral region ranging from 0 to 9 ppm at a high resolution of 0.001 ppm. The region surrounding the water peak (δ = 4.6-5 ppm) was discarded. NMR intensities were log₁₀-scaled and extreme outliers were removed (see Methods). Usually, NMR-based GWAS either use a selection of chemical shifts or perform spectral annotation to reduce the NMR data to the underlying metabolite concentrations (for example, [6, 7, 17]). In contrast, our approach uses the signal intensities at almost all chemical shift positions of the NMR spectra.

In the second step, we selected the 500 bins in the NMR spectrum that exhibited the strongest signal of association. In cases where neighboring bins all displayed associations of comparable strength to the same single nucleotide polymorphism (SNP), the bin with the lowest P value was selected. In a second GWAS, we then tested all possible ratios between the intensities at these 500 chemical shift positions for associations with all genetic variants. To our knowledge, this is the first GWAS that uses a ratio-based, hypothesis-free approach with raw NMR spectral data as a phenotype. For both GWAS, age and gender were used as covariates. In the subsequent analysis, we considered only robust associations of frequent SNPs with minor allele frequencies (MAF) >10%.

Overlap with NMR-derived lipid subclasses

Strong correlations of NMR signals to MS-based and further metabolite concentrations

A ¹H-NMR spectrum of a biological sample is the superposition of the resonance spectra of all individual metabolites measured in that sample. The signals of different metabolites thus overlap in the NMR spectrum. Identification and annotation of individual metabolites is a challenging task. Here we compute 'correlation spectra' to visualize the degree of correlation between NMR bins and MS-derived or biochemically determined metabolite concentrations. Note that all measurements (NMR, MS, clinical biochemistry) were performed on blood samples from the same draw and subject. Correlation spectra thus link the chemical identity of the metabolites detected by other platforms to chemical shifts in the NMR spectra, where the chemical shifts may correspond to the same metabolite, but also to a metabolite on a related biochemical pathway. As an example, Figure 1 shows the correlation spectrum for MS determined glycine. We also observed a number of correlations with phosphatidylcholines (PCs) (Figure 2), and with parameters that were determined independently using standard clinical biochemistry methods, namely triglycerides (TG), high-density lipoproteins (HDL), low-density lipoproteins (LDL), total cholesterol (TC), and glucose (Figure 3), and a number of other metabolites. Most of the correlations we observed are positive, meaning that a stronger NMR signal at the given chemical shift goes with a higher concentration of the correlated metabolite, as one would expect in cases where MS and NMR are targeting identical or closely-related metabolites. In a few cases, we also observe negative correlations between NMR bins and MS-determined metabolites. This may indicate cases where the NMR signal is indirectly related to a MS-measured metabolite. It could, for instance, correspond to the product or substrate of a metabolite quantified on the MS platform.

In the case of correlation with lipid-related parameters, the correlation spectra generally show wider areas of the NMR spectrum corresponding to these metabolites (for example, PCs in Figure 2, triglyceride levels in Figure 3). In contrast, correlations with non-lipid metabolites are visible as individual, sharp peaks (for example, glycine in Figure 1 and glucose in Additional file 2). These peaks correlate, among others, with the concentrations of glucose, lactate, proline, and glycine with squared Spearman correlation coefficients (r_s²) up to 0.61. In most cases, we can verify the validity of the correlations using experimental spectra of pure compounds from HMDB [18]. For glycine, we conducted an additional spiking experiment to confirm that the NMR signal at δ = 3.599 ppm indeed is driven by glycine. However, since many metabolites are by nature interlinked in metabolic networks, some signal correlations do not reflect just a single compound but a mixture of co-regulated (and therefore auto-correlated) metabolites. Also, due to overlapping signals in NMR spectra, there may be multiple, non-related metabolites showing correlations at the same chemical shift position.

Additional file 3, Table S1 lists all metabolites that correlate with NMR signals at r_s² = 0.20 (or above) and the chemical shift where the correlation coefficient is highest. The corresponding correlation spectra are provided in Additional file 2.

Singular NMR intensities (peaks) correlate with MS-based glycine measurements and provide comparable genetic associations (CPS1 locus)

The carbamoyl-phosphate synthase 1 (CPS1) controls the first step in the urea cycle. Klaus et al. report two mutations of the CPS1 gene that contribute to the onset of CPS1 deficiency, an inborn error of metabolism that causes hyperammonemia [19]. In 2010, Illig et al. reported an association of SNP rs2216405 with glycine concentrations at a nominal P value of 2.59×10^-26 [4]. Glycine is metabolically related to carbamoyl phosphate, which in turn is the product of CPS1. Thus, Illig et al. presented a genetically influenced metabotype that might represent a mild form of CPS1 deficiency.

Here, we find an association of the same SNP with the NMR signals at δ = 3.599 ppm (P = 4.46×10^-14). The correlation between MS determined glycine and the NMR signals at that chemical shift is r_s² = 0.22. Although this correlation is modest, the correlation spectrum for glycine shows a distinct peak at δ = 3.599 ppm (Figure 1). A reference spectrum for glycine taken from HMDB [18] shows a single peak at δ = 3.54 ppm. While the spectra used in this study were referenced to lactate and the sample pH was 7.4, the HMDB spectrum was referenced to 4,4-dimethyl-4-silapentane-1-sulfonic acid (DSS) and the sample pH was 7.0. Thus, the 0.06 ppm distance between the putative glycine peak in this study and the HMDB glycine peak can be explained by different experimental conditions. In a separate spiking experiment comparing two spectra with and without addition of glycine, we confirmed that the signal at δ = 3.599 ppm indeed corresponds with glycine.

Thus, in this case we find that information obtained from a GWAS with NMR data provides comparable results to what was obtained with MS, albeit the strength of the association observed with NMR data is weaker in this case. The NMR chemical shift δ = 3.599 ppm could potentially be used as a marker for mild forms of perturbations in the ammonia metabolism caused by genetic variations in the CPS1 locus.

Associations with ratios between NMR intensities indicate spectral regions which are representative of triple and quadruple fatty acid desaturation (FADS1 locus)

The fatty-acid desaturase 1 (FADS1) gene product catalyzes the desaturation reaction of 8,11,14-eicosatrienoyl-CoA to arachidonoyl-CoA (C20:3 → C20:4). Genetic variants in this locus have been linked to Crohn's disease [8] and risk factors for cardiovascular disorders, namely cholesterol and triglyceride levels [9]. Recently, two independent studies showed that dietary intake of long-chain polyunsaturated fatty acids (for example, C20:3 and C20:4) modulates the association between genetic variation in FADS1 and serum lipid levels, and thereby potentially also modifies the risk of cardiovascular disease [20, 21]. In their previous GWAS, Gieger et al. [3] and Illig et al. [4] identified strong associations of SNP rs174547 (intronic region of FADS1) with a number of glycerophospholipids, many containing lipid side chains with C20:3 and C20:4 polyunsaturated omega-3 and omega-6 fatty acids (PUFAs). The authors observed an exceptionally large increase in the strength of the association when ratios between phospholipids containing PUFAs with <4 double bonds and ≥4 double bonds were testedin pairs [4]. This observation can be explained by the fact that ratios between substrate: product pairs approximate the underlying enzymatic reaction rate and also because ratios between related metabolites reduce the overall variance that is observed between individuals with different levels of overall blood PUFA concentrations [12].

In our data, the correlation spectra for glycerophospholipid levels show distinct differences depending on whether the degree of saturation of the side chains is below four double bonds (precursor of substrates of FADS1 enzymatic reaction) or equal or above four double bonds (products downstream of FADS1) (Figure 2). When testing binned NMR intensities for association, SNP rs174547 displays associations at a wider range of chemical shifts, with the strongest signal at δ = 2.801 ppm (P = 3.96×10^-35). Interestingly, when using ratios between the NMR intensities at δ = 2.801 ppm and 2.017 ppm, the strength of association increases by nearly 60 orders of magnitude (P = 1.1×10^-94), similar to the MS-based case. The NMR ratio δ = 2.801 ppm/δ = 2.017 ppm is therefore a likely proxy for the ratio between three- and four-fold desaturated long chain fatty acids. As a metabolic marker, this ratio might be a readout for the efficacy of the FADS1 enzymatic reaction in the context of dietary intake of PUFAs.

A ratio between NMR intensities in the triglyceride and in the glucose range of the spectrum constitutes an integrated pleiotropic diabetes risk marker (GCKR locus)

The product of the glucokinase regulator (GCKR) gene both transports and regulates glucokinase, a key enzyme of glucose metabolism [22]. GCKR is a genetic risk locus for diabetes-related traits [10].

In our study, SNP rs780094 (intronic region of GCKR) associates with NMR intensities at δ = 1.370 ppm. At this spectral position, we observe a high correlation between NMR signals and TG levels (r_s² = 0.61) (Figure 3). However, this association is not of genome-wide significance (P = 1.2×10^-10). When testing ratios for association, the strength of association increases by five orders of magnitude, with a P value of 2.8×10^-15 for the association of SNP rs780094 to the NMR ratio δ = 1.370 ppm/δ = 3.286 ppm. The correlation spectrum with glucose shows a distinct signal at δ = 3.286 ppm (r_s² = 0.28). This chemical shift is indeed located within a known glucose-related spectral region. Furthermore, this NMR ratio correlates with the triglyceride/glucose ratio (r_s² = 0.42). Interestingly, GCKR is a diabetes risk locus that inversely modulates triglyceride and fasting glucose levels [4, 23–26].

Thus, the combination of these two chemical shifts likely provides a combined measure of two main diabetes-risk readouts that have been independently associated with the GCKR locus before. We suggest that this ratio therefore may constitute an integrated biomarker for the pleiotropic biological processes related to perturbations in the GCKR pathway.

Ratios between neighboring NMR intensities act as a local baseline correction and strengthen the genetic association (PYROXD2 locus)

PYROXD2 is a probable pyridine nucleotide-disulphide oxidoreductase gene. Nicholson et al. recently found this locus to be associated with dimethylamine concentrations in plasma, also using NMR spectrometry [6]. At the time when our study was conducted, metabolomics information was only accessible using NMR methods. Meanwhile, an association with an MS determined, biochemically non-identified metabolite has been found, thus replicating this GIM on a different metabolomics platform [27]. However, a link of this metabotype to a phenotype of clinical relevance has not been reported so far.

We find an association between rs4488133 in PYROXD2 and the signal intensities at δ = 2.757 ppm. With a P value of 7.3×10^-12, this association only slightly misses the threshold for genome-wide significance. At the given chemical shift, we do not find any noteworthy correlation to one of the known MS determined metabolites. Most interestingly, similar to the FADS1 case, this locus displays an exceptionally increase in strength of associations when using ratios (>80 orders of magnitude). However, as opposed to the FADS1 locus, the ratio here is between two chemical shift positions in direct neighborhood (δ = 2.757 ppm and 2.755 ppm). An inspection of the detailed spectrum differentiated by genotype reveals that only the signal at δ = 2.757 ppm actually shows clear genotype dependence, with a peak only detected in the major allele homozygote group (Figure 4). This peak is located on the shoulder of a much larger underlying NMR signal that correlates with cholesterol levels. We argue that in this case the ratio between signals at two neighboring chemical shifts, one of them in a localized peak, is equivalent to applying a local baseline correction, thus elevating the peak out of the background noise. Wei et al. have described this effect in a recent publication where they present a ratio-based approach that effectively raises the signal-to-noise ratio in NMR spectra [28].

Since PYROXD2 is a relatively uncharacterized genetic locus, combining the NMR ratio presented here with findings from other studies and platforms could be used to further investigate its biological function.