Genetics and Not Shared Environment Explains Familial Resemblance in Adult Metabolomics Data

Abstract

Metabolites are small molecules involved in cellular metabolism where they act as reaction substrates or products. The term ‘metabolomics’ refers to the comprehensive study of these molecules. The concentrations of metabolites in biological tissues are under genetic control, but this is limited by environmental factors such as diet. In adult mono- and dizygotic twin pairs, we estimated the contribution of genetic and shared environmental influences on metabolite levels by structural equation modeling and tested whether the familial resemblance for metabolite levels is mainly explained by genetic or by environmental factors that are shared by family members. Metabolites were measured across three platforms: two based on proton nuclear magnetic resonance techniques and one employing mass spectrometry. These three platforms comprised 237 single metabolic traits of several chemical classes. For the three platforms, metabolites were assessed in 1407, 1037 and 1116 twin pairs, respectively. We carried out power calculations to establish what percentage of shared environmental variance could be detected given these sample sizes. Our study did not find evidence for a systematic contribution of shared environment, defined as the influence of growing up together in the same household, on metabolites assessed in adulthood. Significant heritability was observed for nearly all 237 metabolites; significant contribution of the shared environment was limited to 6 metabolites. The top quartile of the heritability distribution was populated by 5 of the 11 investigated chemical classes. In this quartile, metabolites of the class lipoprotein were significantly overrepresented, whereas metabolites of classes glycerophospholipids and glycerolipids were significantly underrepresented.

BBMRI Metabolomics Consortium

Cohort Collection and Sample Management Group: M. Beekman¹, H. E. D. Suchiman¹, N. Amin², J. W. Beulens^3,4, J. A. van der Bom^5–8, N. Bomer⁹, A. Demirkan², J. A. van Hilten¹⁰, J. M. T. A. Meessen¹¹, R. Pool¹², M. H. Moed¹, J. Fu^13,14, G. L. J. Onderwater¹⁵, F. Rutters³, C. So-Osman¹⁰, W. M. van der Flier^3,16, A. A. W. A. van der Heijden¹⁷, A. van der Spek², F. W. Asselbergs¹⁸, E. Boersma¹⁹, P. M. Elders^20,21, J. M. Geleijnse²², M. A. Ikram^2,23,24, M. Kloppenburg^8,25, I. Meulenbelt¹, S. P. Mooijaart²⁶, R. G. H. H. Nelissen²⁷, M. G. Netea^28,29, B. W. J. H. Penninx^21,30, C.D.A. Stehouwer^31,32, C.E. Teunissen³³, G. M. Terwindt¹⁵, L. M. ‘t Hart^1,3,21,34,35, A. M. J. M. van den Maagdenberg³⁶, P. van der Harst⁸, I. C. C. van der Horst³⁷, C. J. H. van der Kallen^31,32, M. M. J. van Greevenbroek^31,32, W. E. van Spil³⁸, C. Wijmenga¹³, A. H. Zwinderman³⁹, A. Zhernikova¹³, J. W. Jukema⁴⁰

Database & Catalogue: J. J. H. Barkey Wolf¹, M. Beekman¹, D. Cats¹, H. Mei^1,41, M. Slofstra¹³, M. Swertz¹³

Quality Control: E. B. van den Akker^1,42,43, J. J. H. Barkey Wolf¹, J. Deelen^1,44, M. J. T. Reinders^42,43

Steering Committee: D. I. Boomsma^21,45, C. M. van Duijn², P.E. Slagboom¹

¹Department of Molecular Epidemiology, Leiden University Medical Center, Leiden, the Netherlands.

²Department of Epidemiology, Erasmus MC University Medical Center, Rotterdam, the Netherlands.

³Department of Epidemiology and Biostatistics, Amsterdam University Medical Center, Amsterdam, the Netherlands.

⁴Julius Center for Health Sciences and Primary Care, University Medical Center Utrecht, Utrecht, the Netherlands.

⁵Centre for Clinical Transfusion Research, Sanquin Research, Leiden, the Netherlands.

⁶Jon J van Rood Centre for Clinical Transfusion Research, Leiden University Medical Centre, Leiden, the Netherlands.

⁷TIAS, Tilburg University, Tilburg, the Netherlands.

⁸Department of Clinical Epidemiology, Leiden University Medical Centre, Leiden, the Netherlands.

⁹Department of Cardiology, University Medical Center Groningen, University of Groningen, Groningen, the Netherlands.

¹⁰Center for Clinical Transfusion Research, Sanquin Research, Leiden, the Netherlands.

¹¹Department of Orthopedics, Leiden University Medical Centre, Leiden, The Netherlands.

¹²Department of Biological Psychology, Vrije Universiteit, Amsterdam, the Netherlands.

¹³Department of Genetics, University Medical Center Groningen, University of Groningen, Groningen, the Netherlands.

¹⁴Department of Pediatrics, University Medical Center Groningen, University of Groningen, Groningen, the Netherlands.

¹⁵Department of Neurology, Leiden University Medical Center, Leiden, the Netherlands.

¹⁶Department of Neurology and Alzheimer Center, Neuroscience Campus Amsterdam, VU University Medical Center, Amsterdam, the Netherlands.

¹⁷Department of General Practice, The EMGO Institute for Health and Care Research, VU University Medical Center, Amsterdam, the Netherlands.

¹⁸Department of Cardiology, Division Heart and Lungs, University Medical Center Utrecht, Utrecht, The Netherlands Julius Center for Health Sciences and Primary Care, University Medical Center Utrecht, Utrecht, the Netherlands.

¹⁹Thorax Centre, Erasmus Medical Centre, Rotterdam, the Netherlands.

²⁰Department of General Practice and Elderly Care Medicine, VU University Medical Center, Amsterdam, the Netherlands.

²¹Amsterdam Public Health Research Institute, VU University Medical Center, Amsterdam, the Netherlands.

²²Division of Human Nutrition and Health, Wageningen University, Wageningen, the Netherlands.

²³Department of Radiology, Erasmus University Medical Center Rotterdam, Rotterdam, the Netherlands.

²⁴Department of Neurology, Erasmus University Medical Center Rotterdam, Rotterdam, the Netherlands.

²⁵Department of Rheumatology, Leiden University Medical Center, the Netherlands.

²⁶Department of Internal Medicine, Division of Gerontology and Geriatrics, Leiden University Medical Centre, Leiden, the Netherlands.

²⁷Department of Orthopaedics, Leiden University Medical Center, Leiden, the Netherlands.

²⁸Department of Internal Medicine, Radboud Center for Infectious Diseases, Radboud University Medical Center, Nijmegen, the Netherlands.

²⁹Department for Genomics & Immunoregulation, Life and Medical Sciences Institute (LIMES), University of Bonn, Bonn, Germany.

³⁰Department of Psychiatry, VU University Medical Center, Amsterdam, the Netherlands.

³¹Department of Internal Medicine, Maastricht University Medical Center (MUMC+), Maastricht, the Netherlands.

³²School for Cardiovascular Diseases (CARIM), Maastricht University, Maastricht, the Netherlands.

³³Neurochemistry Laboratory, Clinical Chemistry Department, Amsterdam University Medical Center, Amsterdam Neuroscience, the Netherlands.

³⁴Department of Cell and Chemical Biology, Leiden University Medical Center, Leiden, the Netherlands.

³⁵Department of General practice, Amsterdam University Medical Center, Amsterdam, the Netherlands.

³⁶Department of Human Genetics, Leiden University Medical Center, Leiden, the Netherlands.

³⁷Department of Critical Care, University Medical Center Groningen, Groningen, the Netherlands.

³⁸UMC Utrecht, Department of Rheumatology & Clinical Immunology, Utrecht, the Netherlands.

³⁹Department of Clinical Epidemiology, Biostatistics, and Bioinformatics, Academic Medical Centre, University of Amsterdam, Amsterdam, the Netherlands.

⁴⁰Department of Cardiology, Leiden University Medical Center, Leiden, the Netherlands.

⁴¹Sequencing Analysis Support Core, Leiden University Medical Center, Leiden, the Netherlands.

⁴²Leiden Computational Biology Center, Leiden University Medical Center, Leiden, the Netherlands.

⁴³Department of Pattern Recognition and Bioinformatics, Delft University of Technology, Delft, the Netherlands.

⁴⁴Max Planck Institute for Biology of Ageing, Cologne, Germany.

⁴⁵Netherlands Twin Register, Department of Biological Psychology, Vrije Universiteit, Amsterdam, the Netherlands.

Metabolites (small molecules involved in biological processes) are important intermediates in understanding how a person’s genotype translates to health or disease. Many different types of metabolites can be distinguished, including amino acids, lipids, and sugars (Adamski & Suhre, 2013). Due to their diversity, metabolites have various functions in the human body, including energy storage, signaling, and forming structures, such as proteins or cell walls (Dunn et al., 2011). Thus, metabolites can be considered the building blocks of the body (Hasirci & Hasirci, 2018). The complete set of metabolites found within a specific biological sample is referred to as the metabolome (Wishart, 2007), the study thereof is termed metabolomics (Fiehn, 2002). The metabolome is downstream of gene transcription, protein translation and protein function; therefore, metabolites are close to observable phenotypes in health and disease (Draisma et al., 2013; Goodacre et al., 2004). Metabolomics has been successful in identifying disease biomarkers, unraveling biological mechanisms, and for drug discovery and development (Pang et al., 2019).

The metabolome differs between people, as metabolite levels are influenced by many exogenous (originating from outside an organism) and endogenous factors (originating from inside an organism). Exogenous factors influencing the human metabolome include lifestyle, diet or medication use. For example, metabolite levels differ between current-, former- and never-smokers, between individuals on a low fat, low glycemic, or very low carbohydrate diet, and among users of various medication classes (Esko et al., 2017; J. Liu et al., 2020; Xu et al., 2013). Endogenous factors influencing the human metabolome include sex, age or body mass index (BMI). For example, metabolite levels differ between males and females, younger and older individuals and obese and nonobese individuals (Chaleckis et al., 2016; Fan et al., 2018; Rangel-Huerta et al., 2019). Endogenous factors also include genetic influences, either directly on metabolite levels or indirectly through the effect on behavior or lifestyle (e.g., smoking; M. Liu et al., 2019). The metabolome can differ between cases and controls — for example, in major depressive disorder (Bot et al., 2019). Here it is observed that associations between metabolites and the case/control status are attenuated by antidepressant use, while the causalities of the associations remain unknown as of yet. Thus, metabolite levels reflect individual differences in genetics, physiology, lifestyle and behavior or responses to environmental changes (Fiehn, 2002).

Genetic factors account for approximately 50% of the individual differences in metabolite levels (Shin et al., 2014; Yet et al., 2016). The average proportion of phenotypic variance in metabolite levels ascribed to genetic factors (i.e., heritability; h ²) differs per type of metabolite. The median heritability for lipids is approximately 37%, with heritability estimates for sphingolipids (e.g., sphingomyelins) and glycerolipids (e.g., triglycerides) often, but not consistently, higher than for phospholipids (e.g., phosphatidylcholines; Bellis et al., 2014; Darst et al., 2019; Frahnow et al., 2017). Similarly, while the median heritability for amino acids is approximately 40% (Darst et al., 2019), amino acids that the body is able (nonessential) or unable (essential) to synthesize de novo differ in mean heritability. Specifically, levels of essential amino acids are less heritable than levels of nonessential amino acids (Rhee et al., 2013). These differences in heritability among metabolite classes also occur for single-nucleotide polymorphism (SNP)-based heritability (i.e., heritability estimates derived from genomewide SNPs; Rhee et al., 2016; Tabassum et al., 2019). In fact, the observed differences in additive heritability estimates among metabolite classes are rarely significant, while differences in heritability estimates based on known genetic variants are frequently significant (Hagenbeek et al., 2020).

Whereas the contribution of genetic variants to metabolite levels is fairly well established, also through genome-wide association and (exome-) sequencing studies (Hagenbeek et al., 2020; Kastenmüller et al., 2015; Yazdani et al., 2019), the contribution of the shared or common environment (c ²) to metabolite levels is not as well characterized. Not all studies investigating the genetic contribution to metabolite levels were based on the classic twin design (e.g., Draisma et al., 2013; Tremblay et al., 2019). Instead, these types of methods estimate the familial resemblance (or generalized heritability) of metabolite levels, which comprises both additive genetic effects and common environment effects shared by family members (Rice, 2008). Studies using the classic twin design to investigate the contributions to metabolite levels vary widely with respect to how metabolites are influenced by common environment as well as in the estimate of the effects of the common environment. Overall, it would appear that studies in smaller samples more frequently report larger contributions of the common environment to metabolite levels (e.g., Frahnow et al., 2017; Kettunen et al., 2012). Large-scale twin studies that estimate the contribution of shared environment tend to be scarce.

In the current study, we aim to expand our understanding of the contribution of common environment shared by family members to variation in fasting blood metabolites and analyzed data from multiple metabolomics platforms from a large cohort of twins (between 1037 and 1407 twin pairs per platform), representing a general population. First, a series of power analyses were performed, estimating the statistical power to detect shared environment in the classical twin design, given the number of monozygotic (MZ) and dizygotic (DZ) twin pairs available in our study. The power to detect shared environment in quantitative genetic studies, employing the classical twin design, is influenced by effect size, the heritability of the trait, the sample size, the probability level that is chosen, and the homogeneity of means and variances in the MZ and DZ groups of the sample (Martin et al., 1978; Posthuma & Boomsma, 2000). Sample size and probability were given, and we investigated different values for the proportion of variation explained by shared environment against a background of different heritability values that were chosen based on what is typically reported in the literature for metabolomics. Next, we determined the heritability of all metabolic traits by structural equation modeling where contributions to additive genetic effects (A), shared environmental effects between siblings (C) and unique environmental effects (E) were estimated. By computing the significance of the C variance component when comparing ACE model outcomes to AE model outcomes, we were able to assess whether C contributes to the total variance observed in all metabolic traits. Finally, to obtain insight into the distribution of chemical classes our metabolites belong to over the range of calculated heritabilities and contributions of shared environment, we performed enrichment analyses. The metabolic traits were grouped in heritability and shared environment estimate quartiles. By counting the number of metabolites of a given class per quartile and comparing these to the counts of the entire range of heritabilities or common environment contributions, we determined class enrichment factors per quartile and assessed their statistical significance by Fisher’s exact tests.

Methods

Participants

At the Netherlands Twin Register (NTR; Ligthart et al., 2019) metabolomics data for twins were available for 886 complete MZ pairs and 601 complete DZ pairs (c.f. Table 1). All measurements were performed in blood samples that were collected from participants of the NTR biobank projects (Willemsen et al., 2010; Willemsen et al., 2013). Blood samples were collected after overnight fasting. Fertile women were bled in their pill-free week or on day 2−4 of their menstrual cycle. For the Nightingale Health metabolomics platform (see below), data were acquired in several shipments (subsets). After completing the preprocessing of the metabolomics data, each platform subset (if applicable) was merged into a single per platform dataset, randomly retaining a single observation per platform whenever multiple observations were available. Characteristics for the sample of individuals included in the analyses can be found in Table 1. Informed consent was obtained from all participants. Projects were approved by the Central Ethics Committee on Research Involving Human Subjects of the VU University Medical Centre, Amsterdam, an Institutional Review Board certified by the US Office of Human Research Protections (IRB number IRB00002991 under Federal-wide Assurance-FWA00017598; IRB/institute codes and NTR 03-180).

Table 1. Participant characteristics by the Nightingale, lipidomics and NMR-LUMC platforms

Note: Age (SD) denotes the mean age and standard deviation, BMI (SD) the mean BMI and standard deviation, F_Fml the fraction of female subjects, F_Smkng the fraction of current smokers, F_LLMd the fraction of subjects using lipid lowering medication, N _MZPrs the number of complete monozygotic twin pairs and N _DZPrs the number of complete dizygotic twin pairs.

Metabolite Profiling

Below, we briefly describe the methods used for metabolite profiling. For more detailed information, see Hagenbeek et al. (2020).

Nightingale Health ¹H-NMR platform

Metabolic biomarkers were quantified from plasma samples using high-throughput proton nuclear magnetic resonance spectroscopy (¹H-NMR) metabolomics (Nightingale Health Ltd, Helsinki, Finland). This method provides simultaneous quantification of routine lipids, lipoprotein subclass profiling with lipid concentrations within 14 subclasses, fatty acid composition, and various low-molecular weight metabolites, including amino acids, ketone bodies and glycolysis-related metabolites in molar concentration units. Details of the experimentation and epidemiological applications of the NMR metabolomics platform have been reviewed previously (Soininen et al., 2015; Würtz et al., 2017).

UPLC-MS lipidomics platform

Plasma lipid profiling was performed at the division of Analytical Biosciences at the Leiden Academic Center for Drug Research at Leiden University/Netherlands Metabolomics Centre. The lipids were analyzed with an Ultra-High-Performance Liquid Chromatograph directly coupled to an Electrospray Ionization Quadruple Time-of-Flight high-resolution mass spectrometer (UPLC-ESI-Q-TOF; Agilent 6530, San Jose, CA, USA) that uses reference mass correction. For liquid chromatographic separation, a ACQUITY UPLC HSS T3 column (1.8 μm, 2.1 × 100 mm) was used with a flow of 0.4 ml/min over a 16-min gradient. Lipid detection was performed using a full scan in the positive ion mode. The raw MS data were pre-processed using Agilent MassHunter Quantitative Analysis software (Agilent, Version B.04.00). Detailed descriptions of lipid profiling and quantification have been described previously (Dane et al., 2014; Gonzalez-Covarrubias et al., 2013).

Leiden ¹H-NMR platform (for small metabolites)

The Leiden ¹H-NMR spectroscopy experiment of ethylenediaminetetraacetic acid plasma samples used a 600-MHz Bruker Advance II spectrometer (Bruker BioSpin, Karlsruhe, Germany). The peak deconvolution method used for this platform has been previously described (Demirkan et al., 2015; Verhoeven et al., 2017).

Metabolomics Data Preprocessing

To ensure our data were consistent with Hagenbeek, Pool, van Dongen, Draisma, Boomsma et al. (2020), we excluded participants if they were on lipid-lowering medication at the time of blood draw or if they had not adhered to the fasting protocol (~4 % of the sample of each platform). Preprocessing of the metabolomics data was executed for each of the platforms and measurement/shipment batches per platform separately. Metabolites were excluded from analysis when the mean coefficient of variation exceeded 25% or the missing rate exceeded 5%. Metabolite measurements were set to missing if they were below the lower limit of detection or quantification or could be classified as an outlier (five standard deviations greater or smaller than the mean). Metabolite measurements that were set to missing because they fell below the limit of detection/quantification were imputed with half of the value of this limit, or when this limit was unknown with half of the lowest observed level for this metabolite. All remaining missing values were imputed using multivariate imputation by chained equations (‘mice’; van Buuren & Groothuis-Oudshoorn, 2011). On average, nine values had to be imputed for each metabolite (SD = 12; range: 1−151). Data for each metabolite on the lipidomics platform and both ¹H-NMR platforms were normalized by inverse normal rank transformation (Demirkan et al., 2015; Kettunen et al., 2016).

We computed heritability for 237 single metabolic traits (i.e., no ratios or composite variables, see Supplementary Table S1). These traits are members of 11 different chemical classes, as listed by the human metabolome database (Hagenbeek et al., 2020; Wishart et al., 2018). As is shown in Table 2, most metabolites are lipid species.

Table 2. Chemical class counts for each platform used

Note: The numbers in the platform columns represent the number of compounds of a chemical class, where ‘−’ means that there are none.

To account for confounding by age and sex, we used the residuals of the linear fit of model M _i~Age + Sex for each metabolite i (M_i) as input for the statistical analyses. The above data processing steps were performed in the Jupyterlab environment (v0.35.4), running an IPython kernel (v5.1.1: python v3.7.3) and utilizing modules pandas (v0.24.2), scipy (v1.2.1) and statsmodels (v0.9.0). The complete lists with full names of all detected metabolites that survived QC and preprocessing for all platforms can be found in Supplementary Table S1.

Statistical Analyses

Power analyses

We calculated covariance matrices for multiple combinations of heritability and the proportion of variation explained by common environment (C). Background heritability differed between 0.2 and 0.7, which are typical values estimated for metabolites (Hagenbeek et al., 2020). Within heritability class, the effect of C was increased from 0.1 to 0.3. Power analyses were carried out in the statistical software package Mx (Maes et al., 2009; Neale, 1997), with estimation of parameters by normal theory maximum likelihood. Goodness-of-fit testing was based on likelihood ratio tests. First, an ACE model was considered, and next the influence of C was constrained at 0 and power was obtained for a 1 degree of freedom test with p = .05. Note that we do not take into account the number of tests performed, which would lower the probability level (and hence lower the power to detect C).

Genetic analyses

For each metabolite M _i, we estimated the additive genetic contribution (A) and the contributions of common and unique environment (C and E). We applied genetic structural equation modeling using maximum likelihood estimation on both ACE and AE models. By comparing the outcomes of both models, for each metabolite, we applied a threshold of 3.84 on the χ² statistic, above which we considered the contribution of C significant. The analyses were performed in the R software package (v3.5.2) using the OpenMx (v2.13.2.161) library (applying the NPSOL optimizer) for running the ACE and AE models (Boker et al., 2011).

Analyses of chemical class enrichment over heritability percentile groups

We subdivided the outcomes of the additive genetic components of the heritability or common environment variance component into three groups: (1) the 0−25 percentile group, (2) the 25−75 percentile group and (3) the 75−100 percentile group. In each group, we counted the number of class memberships of the metabolites. Within each group, we performed Fisher’s exact tests (Fisher, 1922) of the group chemical class counts with respect to their counts in the entire sample.

Results

Power Analyses

The outcomes of the power analyses are summarized in Table 3. As can be observed from these results, for all three platforms, when A is in the range of 0.2–0.6, a larger sample size is needed to enable detection of C ≤ 0.2. This implies that when we observe a statistically significant C variance component, its minimal value needs to be approximately 0.2.

Table 3. Statistical power estimates as a function of values for A, C, r_MZ and r_DZ for the metabolomics platforms in this work

Note: A, additive genetic effects; C, shared environmental effects between siblings; E, unique environmental effects; MZ, monozygotic; DZ, dizygotic. For determining the estimates, we applied N(MZ pairs) = 848 and N(DZ pairs) = 559 for the Nightingale platform, N(MZ pairs) = 612 and N(DZ pairs) = 504 for the lipidomics platform, and N(MZ pairs) = 648 and N(DZ pairs) = 389 for the NMR-LUMC platform.

Genetic Analyses

The results of the ACE models for each platform and across the platforms are listed in Table 4. The results of the ACE and AE models are summarized graphically in Figure 1. On average, the analyses included 694 (range 612–848) and 503 (range 389–559) MZ and DZ pairs, respectively. For 6 of the 237 metabolites, we observe a significant contribution of c ² ranging from 0.175 to 0.423 (mean = 0.254). Over 6 ACE models and 231 AE models, the mean additive genetic variance component a ² is 0.456 (range 0.136–0.735) and the mean of the residual error and unique environment variance component e ² is 0.527 (range 0.265–0.864). Supplementary Figures S1 and S2 depict the model outcomes for the AE and ACE models, respectively. Supplementary Figure S2 shows the MZ and DZ twin correlations for all metabolites.

Table 4. Summary of the AE (a) and ACE models (b and c) per platform and combined across platforms. Table 4b lists the ACE component means and variances only for the traits that exhibited a signifacant contribution of C. Table 4c lists the ACE component means and variances irrespective of the significance of the contribution of C

Note: A, additive genetic effects; C, shared environmental effects between siblings; E, unique environmental effects; MZ, monozygotic; DZ, dizygotic. For each metabolomics platform, the additive genetic variance component a ² and the variance component of the common environment c ² are listed. These outcomes are based on N _M metabolites of which N _M(ACE) have a significant C component. The analyses are based on N _MZPairs and N _DZPairs of MZ and DZ twin pairs. The bottom row lists the statistics over all platforms. Values for the individual metabolites are listed in Supplementary Table S1.

Note: A, additive genetic effects; C, shared environmental effects between siblings; E, unique environmental effects.

Fig. 1. Heritabilities for twin AE or ACE models. For the six metabolites showing a red bar, a significant contribution of C was observed (χ² > 3.84). For all other metabolites, only estimates for the additive genetic component are shown, as determined by AE models. On the top of bar plot, two color bars are depicted that indicate the chemical class of the metabolite (top bars) and the metabolomics platform of the metabolite was reported from (bottom bars). The x-axes denote the indices of the metabolites, listed in Supplementary Table S1.

Fig. 2. Parallel coordinates plots of the top and bottom heritability quartlies showing the enrichment of the lipoprotein chemical class in the top percentile and its underrepresentation in the bottom percentile groups. Note that classes glycerophospholipids and glycerolipids are underrepresented in the top percentile group.

Analyses of Chemical Class Enrichment Over Percentile Groups of Heritabilities and Shared Environment Variance Components

Only six metabolic traits exhibit a significant shared environment variance component, rendering it futile to perform an enrichment analysis over C percentile groups. As can be observed in Figure 1, it seems that the lipoprotein chemical class is overrepresented in the top heritability percentile group and underrepresented in percentiles 25–50 and 0–25. This is confirmed by enrichment analyses of the chemical groups across the different heritability percentile groups (Table 5). Moreover, chemical classes glycerophospholipids and glycerolipids are significantly underrepresented in the top percentile group 75–100. The parallel coordinates plots in Figure 2 summarize these findings graphically.

Table 5. Representation of chemical classes across the percentiles of the additive genetic component (A) values as determined by the AE models

Note: The Count columns list the metabolite counts per chemical class. The enrichment statistics are coming from Fisher’s exact test performed over all metabolite class counts per subset with respect to the entire sample. OR columns list the odds ratios per chemical class, p columns the associated p values and p (fdrb) the Benjamini–Hochberg false discovery rate corrected ones.

Discussion

Metabolites have an important role in the relationship between the genotype and health and disease. Therefore, characterizing the factors influencing metabolite levels is a vital first step toward elucidating the mechanisms underlying health and disease. It is well established that metabolite levels are influenced by a complex interplay of genetic and environmental influences; however, the role of the common, or shared, environment in metabolite levels remains unclear. Here, we investigated the contribution of the common environment to variation in fasting blood metabolite levels in a large twin cohort. Although our study had sufficient power to detect a contribution of the common environment to the metabolite levels of 20% or higher, we found little evidence for a contribution of this size. In contrast, but in line with expectations, we found that additive genetics contributed significantly to metabolite levels. We found that the top 25% most heritable metabolites included mostly lipoproteins, while lipoproteins were underrepresented in the other.

Our conclusions need to be placed in the context of the design, platforms and biofluid that were used. First, while we established that with a fixed sample size, a probability level of 5% and depending on the background heritability, a contribution of common environment of 20% can be detected with reasonable power. A contribution of the common environment of 30% had very high statistical power for all platforms. However, when the contribution of the common environment is small, even for the largest sample sizes in our study (Nightingale platform with 848 MZ and 559 DZ twin pairs), the power to detect common environment shared by twins is low. Second, our study was cross-sectional, not longitudinal. While the fasting metabolite levels of individuals are generally stable over time (Lacruz et al., 2018), age is a well-established factor influencing metabolite levels, with higher levels observed in adults as compared to children (Ellul et al., 2019). Third, this work involves metabolic traits measured in blood. Results reported here do not necessarily represent those expected for other tissue types. Fourth, we investigated ‘fasting state’ metabolite levels. As such, this can be considered a relatively stable homeostatic state, in contrast to states of high energy consumption or after food and/or drink intake. Finally, this study, as most others, includes platforms favoring lipids, limiting the conclusions we may draw with regard to nonlipid species (Hagenbeek et al., 2020).

We observed a limited contribution of common environment to the total variation in metabolite levels: we estimated a significant contribution of the common environment only for 6 out of 237 metabolites. For these 6 traits, the average contribution to the common environment was 0.25 (range: 0.17−0.34). Our findings are in line with previous studies reporting an influence of the common environment for a minority of the investigated metabolites, with average contributions of the common environment frequently lower even than the 0.25 we observed (Kettunen et al., 2012; Menni et al., 2013; Yet et al., 2016). The limited contribution of the common environment on metabolite levels is also consistent with observations for other molecular traits, such as expression (Ouwens et al., 2020; Wright et al., 2014) or methylation levels (van Dongen et al., 2016), which were measured in the same group of twin pairs as in this work. This is in stark contrast to a recent family-based study that reported substantial familial resemblance in metabolite levels to which common environment had the strongest contribution (Tremblay et al., 2019). A strong contribution of the common environment has also been observed for specific lipid species, while other lipid species had no or only small contributions of the common environment (Frahnow et al., 2017). It must be noted that those studies reporting substantial influence of the common environment frequently have small sample sizes and are generally younger than our own samples.

The observation that an effect of the common environment is more frequently reported in younger participants is in line with observations for other complex traits. For several complex traits with a strong influence of the common environment during childhood, this decreases greatly with age, often completely disappearing in adulthood (Haworth et al., 2010; Lamb et al., 2010). Currently, longitudinal genetic investigations of metabolite levels starting in childhood are lacking. Therefore, we may not definitely conclude that a similar pattern holds for metabolite levels. Moreover, in general for studies including adult samples, the question arises what the common environment comprises (Finnicum et al., 2019; Zaitlen et al., 2013). In adult twin pairs living apart the common environment may represent lasting influences of their time spend living together. Alternatively, the common environment could reflect similar living environments, including the same (or similar) neighborhood, exposure to similar pollutants or shared lifestyles. Thus, while many environmental factors influencing metabolite levels have been identified, the distinction between common and unique environmental factors remains unclear.

We found an overrepresentation of lipoproteins with heritability estimates in the highest quartile. Interestingly, we observed no differences in total heritability estimates between lipoproteins and other lipid classes in the same sample using genome-wide SNP data; however, the portion of heritability due to known metabolite loci did differ between lipid classes (Hagenbeek et al., 2020). This could indicate that the genetic relationship estimates from genome-wide SNP data do not comprise enough information to assess this enrichment. Combining the enrichment results with our previous findings that class-specific metabolite loci have a larger contribution to metabolite heritability than nonclass loci suggests the merit of investigating the genetics of biochemical pathways rather than individual metabolites. The strength of leveraging knowledge of biochemical pathways in genetic investigations has been previously established. Genetic investigations of metabolite ratios reflecting enzymatic conversions in biochemical pathways have led to stronger associations than studying single metabolites (Suhre et al., 2016). Similarly, multivariate genomewide association studies of correlated metabolic networks have proven effective (Inouye et al., 2012).

Conclusion

Based on structural equation modeling in twins, a significant heritability was found for the blood concentrations of 231 metabolites, explaining 46% of their variance on average. We find that in the top 25% most heritable metabolites, those of class lipoprotein are overrepresented, while metabolites of classes glycerophospholipids and glycerolipids are underrepresented. In contrast, we find little evidence of shared environment influences on individual differences in metabolite levels. These findings show that familial resemblance in the concentrations of metabolites is due to genetic factors with minimal contribution of the shared environment.