Comparative metabolomics in primates reveals the effects of diet and gene regulatory variation on metabolic divergence

Human diets differ from those of non-human primates. Among few obvious differences, humans consume more meat than most non-human primates and regularly cook their food. It is hypothesized that a dietary shift during human evolution has been accompanied by molecular adaptations in metabolic pathways. Consistent with this notion, comparative studies of gene expression levels in primates have found that the regulation of genes with metabolic functions tend to evolve rapidly in the human lineage. The metabolic consequences of these regulatory differences, however, remained unknown. To address this gap, we performed a comparative study using a combination of gene expression and metabolomic profiling in livers from humans, chimpanzees, and rhesus macaques. We show that dietary differences between species have a strong effect on metabolic concentrations. In addition, we found that differences in metabolic concentration across species are correlated with inter-species differences in the expression of the corresponding enzymes, which control the same metabolic reaction. We identified a number of metabolic compounds with lineage-specific profiles, including examples of human-species metabolic differences that may be directly related to dietary differences.


Sample preparation
For this analysis we used 18 liver samples from three primate species (human, chimpanzee, and rhesus macaque), with 6 individuals from each species. Liver samples from non--human primates were collected at necropsy, within four hours of death, by the Yerkes primate center, the Southwest Foundation for Biomedical Research, and MD Anderson Cancer Center. Additional primate tissues were given to us by Anne Stone. In all cases, we collected liver tissue samples from adult chimpanzees and rhesus macaques that died of natural causes (such as accidents or fights) or were euthanized due to an illness unrelated to liver. The human adult liver samples were collected for us by the National Disease Research Interchange (NDRI), and by the pathology department at Yale University (with IRB approval). Detailed information about all samples is available in Table S2.
Tissue samples were immediately frozen and maintained at --80C. From each liver sample, we excised three small tissue pieces, each 100 mg, from different sections of the original, larger sample. This step was conducted on dry ice to avoid thawing the samples.

Metabolomics methodology
To quantify metabolite levels in liver samples, we have applied gas chromatography/time--of--flight mass spectrometry (GC--TOF MS), as this technique yields the largest overview over metabolites smaller than approximately 500 Da, especially the diversity of carbohydrates (mono--, di-- and trisaccharides), sugar alcohols, hydroxyl acids (including intermediates of the tricarboxylic acid cycle), amino acids, aromatics, free fatty acids, and ranges of miscellaneous compounds such as purines and pyrimidines. While there is overlap in metabolite coverage with complementary techniques such as liquid chromatography/mass spectrometry, GC--TOF MS is superior in separating isomeric compounds such as fructose and glucose and has better command over data processing software such as mass spectral deconvolution, data processing algorithms [1] and mass spectral libraries [2]. As every other metabolomic technique, GC--TOF MS is limited in scope; for example, complex lipids such as phosphatidylcholines or thermodegradable metabolites like ATP cannot be analyzed this way.

Metabolomics data acquisition
2 mg liquid--nitrogen frozen liver tissue samples were homogenized at 25 Hz with 3 mm i.d. steel balls for 30 s and immediately placed back into liquid nitrogen, afterwards.
Samples were taken out one by one and 1 ml of a carefully degassed --20°C cold isopropanol/acetonitrile/water mixture (3:3:2, v/v/v) were added and shaken for 5 min at 4°C to extract metabolites and simultaneously precipitate proteins. After centrifugation at 12,800 x g for 2 min, 90% of the supernatant was removed, separated into two equal aliquots, transferred to a 1.5 ml Eppendorf tube and concentrated to dryness in a vacuum concentrator. Samples were cleaned up by adding 500 ul of a degassed acetonitrile/water mixture (1/1, v/v), vortexing, centrifugation, decanting the supernatant and drying in a speed vac concentrator. This step removes triglycerides and most of the complex lipids, but not phytosterols or free fatty acids and is needed because otherwise, the involatile matrix lipids would interfere with the derivatization reaction of primary amines and amino acids. C08--C30 fatty acid methyl esters in chloroform were added as internal retention index (RI) markers. Subsequently, metabolites were derivatized by adding 10 ul of a solution of 20 mg/ml of 98% pure methoxyamine hydrochloride (CAS number 593-56--6, Sigma--Aldrich) in pyridine for 90 min at 28°C. Afterwards 90 ul of N--methyl--N-trimethylsilyltrifluoroacetamide (MSTFA, Sigma--Aldrich) was added for trimethylsilylation of acidic protons and shaken at 37°C for 30 min. The reaction mixture was transferred to a 2 ml clear glass auto--sampler vial with micro--insert (Agilent) and closed using a 11 mm T/S/T crimp cap (MicroLiter). A Gerstel automatic liner exchange system with a MPS2 dual rail multi--purpose sampler was used in conjunction with a Gerstel CIS cold injection system. For every 10 samples, a fresh multi--baffled liner was inserted (Gerstel #011711--010--00). Before and after each injection, the 10 ul injection syringe was washed three times with 10 ul ethyl acetate. 1 ul sample was filled using 39 mm vial penetration at 1ul sec--1 filling speed, injecting 0.5 ul at a 10 ul s--1 injection speed at an initial temperature of 50°C which was ramped by 12°C sec--1 to a final temperature of 250°C and held for 3 min. The injector was operated in split--less mode, opening the split vent after 25 sec. Samples were injected between 2-24 h after derivatization using randomized sequences controlled by the laboratory information management and The instrument performed auto--tuning for mass calibration using FC43 (perfluorotributyl-amine) before starting analysis sequences.

Metabolomics data processing
Chromatogram acquisition, data handling, automated peak deconvolution, and export of spectra was automatically performed by the Leco ChromaTOF software (v2.32).
Peak picking was achieved in ChromaTOF (v2.32) at signal/noise levels of 10:1 throughout the chromatogram with baseline subtraction just above the noise level, no smoothing, 3 s default peak widths, automatic mass spectral deconvolution and peak detection and export of result spectra as *.csv files in addition to export of raw data in open--access *.cdf formats. Data were further processed using the algorithms implemented in the open--source BinBase metabolome database [3]. This algorithm used the settings: validity of chromatogram (<10 peaks with intensity >10^7 counts sec--1), unbiased retention index marker detection (MS similarity >800 and exceeding thresholds for ion ratio abundances for high m/z marker ions), retention index calculation by 5th order polynomial regression. Spectra were cut to 5% base peak abundance, and matched to database entries from most--to least--abundant spectra using the following matching filters: retention index window ±2000 units (equivalent to about ±2 sec retention time), validation of unique ions and apex masses (unique ion must be included in apex masses and present at >3% of base peak abundance), mass spectrum similarity that must fit criteria dependent on peak purity and signal/noise ratios, optional ion ratio settings to distinguish peaks with high similarity, and a final isomer filter (annotating the isomer spectrum with the closest RI fit). Novel spectra not yet included in BinBase were automatically entered as new database entries if their signal--to--noise ratio >25, purity <1.0 and presence in the biological study design class was >80%. This filter ensured that (i) signals were reported that had never been detected previously in any other sample, but (ii) only signals were reported that can be assumed to be biologically relevant using relatively abundant and pure signals and ensuring that these are positively detected in most of the biological replicates. Signal intensities were reported as peak heights using the unique ion as default, unless an alternative quantification ion was manually set in the BinBase administration software Bellerophon. All known artifact peaks such as internal standards, column bleed, plasticizers or reagent peaks were assigned by BinBase but not exported for further statistical calculations. Metabolites were identified using the Fiehnlib libraries consisting of over 1,200 authentic compounds and referenced using PubChem identifiers [2]. Daily quality controls were used comprising method blanks and five--point calibration curve samples of 31 pure reference compounds which were repeatedly analyzed over the full analytical sequence in addition to injection of one QC sample for every 10 biological samples. A quantification report table was produced for all database entries that were positively detected in more than 50% of the samples of a study design class. This procedure results in 10-30% missing values, which could be caused by true negatives (compounds that were below the detection limit in a specific sample) or false negatives (compounds that were present in a specific sample but that did not match quality criteria in the BinBase algorithm). A subsequent post--processing module was employed to automatically replace missing values from the *.cdf files with the following parameters: for each positively detected metabolite, the average retention time was calculated for the day of analysis. Subsequently, for each chromatogram and each missing value, the intensity of the quantification ion at this retention time was extracted by seeking its maximum value in a retention time region of ± 1 s and subtracting the minimum (local background) intensity in a retention time region of ±5 s around the peak maximum. The resulting report table therefore did not have any missing values.

Metabolomics data normalization and statistics
Result files were normalized by calculating the sum intensities of all structurally identified compounds for each sample (i.e. those signals that had been positively identified in the data pre--processing schema outlined above), and subsequently dividing all data associated with a sample by the corresponding metabolite sum. The resulting data were multiplied by the average sum of all identified metabolites detected in the study (total average metabolome transformation), disregarding unknown metabolites as these might potentially also comprise artifacts. Intensities of identified metabolites with more than one peak (e.g. for the syn--and anti--forms of methoximated reducing sugars, or amino acids with different derivatization status of amine groups) were summed to only one value in the transformed data set. The original non--transformed data set was retained for retrospective analysis. When comparing classes of samples with biologically different sum concentrations of identified metabolites (p<0.05), these class averages were used for mean transformations. The final metabolomics data consists of three replicates for each sample, for a total of 54 (3*6*3) measurements for each metabolite, and is available as a separate data file (Dataset S1). The dataset included information for 399 metabolites detected in all samples, of which 177 had an associated name, and 153 were also found in the KEGG database.

Expression data
We used expression data that we previously collected and analyzed [4]. Briefly, the data was collected using a multi--species microarray, containing orthologous probes from three primate species: human, chimpanzee, and rhesus macaque. The array contains probes for 18,109 genes (368,678 probes in total). The data includes gene expression estimates from six individuals from three tissues (liver, kidney cortex and heart muscle), from each of the three species. Complete information on sample collection, study design, array hybridizations, low--level analysis, and quality control is available in Blekhman et al. (2008).

Low--level analysis and quality control of the metabolite data
First, we normalized the metabolite data using the quantile normalization approach, with the normalize.quantiles function in the R library preprocessCore version 1.8. To assess to quality of the data, we performed principal component (PC) analysis on the normalized data, using the princomp function in R. and plotted the first vs second principal components (figure S1C). We note that the PCs separate the samples based on the batch (replicate) in which it was run, with PC1 splitting the first batch from the rest, and PC2 separating the second and third batches. Since this variation is technical, and could lead to spurious results, we decided to remove this effect and correct for it before moving on with the analysis. We also note that PC3 does not correspond to any expected biological division of data, while PC4 separates the samples based on their species (see figure S1D). Since the species effect is expected to be the main factor explaining the variance between samples, we concluded that PC3 is a technical effect, and decided to correct for this unexplained effect as well.
To do so, we analyzed the normalized data with a metabolite--wise linear model including the replicate and third PC as factors. We generated the corrected data by summing the intercept and residuals for each metabolite, thus regressing out the model factors, namely the batch and the third PC.
In addition to correcting for non--biological variance, we also used the PC analysis, together with a heatmap plot of sample pairwise correlations, to identify outliers in the metabolite data. We found three outlying samples that did not cluster with their other two replicates: (1) human sample 56720 replicate 2 (can be clearly seen in figure S1C), (2) human sample 56655 replicate 3, and (3) rhesus macaque sample YN05--349 replicate 3.
After excluding these samples from the data, we repeated the PC analysis, and found no visible confounding effect and outliers (figure S1A and S1C). Moreover, we see the species as the first and second PCs. In addition, a heatmap plot of the normalized, corrected, outlier--excluded data shows a clear separation ( figure S2).

Metabolic pathway data
All metabolic pathway information was downloaded from KEGG [5]

Identifying differences in metabolite concentrations between species
To identify metabolites with differences in concentrations between species, we applied the following linear model on the levels of each metabolite:

ysij = μs + αsi + εsij
where ysi is the normalized log concentration level from species s in individual i and replicate j, μs is the expected log concentration level in species s, αi is a random effect to capture the variance between individuals within each species, and εsi is the error term.
To identify metabolites differentially concentrated (DC) between species we used a set of three likelihood ratio tests, each comparing the likelihood of the full model above with that of a reduced model, which assumes there is a similar concentration level across the species. We calculated a likelihood ratio (LR) statistic for each pairwise comparison, and estimated a P--value under the assumption that the LR has a chi--square distribution with one degree of freedom. We then calculated the false discovery rate associated with these P--values using the approach of [6]. Finally, we used a threshold of q--value < 0.05 and identified 122, 96, and 29 metabolites with different concentrations between human and chimpanzee, human and rhesus macaque, and chimpanzee and rhesus macaque, respectively (see Figure  S3). Given the small number of differences between chimpanzee and rhesus macaque, we focused most of the following analyses on the human-chimpanzee and human--rhesus macaque comparisons, which were more informative.
Next, we aimed to find metabolites with a human--specific concentration level, i.e., metabolites that show a significantly higher (or lower) level in human, but a similar level in the two non--human species. To do so, we combined the results of the three pairwise tests described above, and picked genes with a P<0.01 for human--chimpanzee and human-rhesus macaque comparison, and P>0.1 for chimpanzee--rhesus macaque. We identified 35 metabolites with the expected pattern, of which 21 have a known name (shown in Figure   1).

Human--specific metabolic concentrations in pathways
We then wanted to identify KEGG pathways that show evidence for human--specific metabolic concentrations. To do so, for each pathway in KEGG we combined the pairwise P--values for differential concentration over all the metabolites in the pathway using Fisher's method, and considered pathways with a combined, Bonferroni--corrected P<0.05 for the human--chimpanzee and human--rhesus macaque differences, and combined P>0.1 for the chimpanzee--rhesus difference. Table S3 lists the most enriched pathways, including the metabolites included in the analysis in each pathway, and the combined pairwise P--values.

Identifying differentially expressed genes
For all subsequent analyses we excluded probes that did not have corresponding orthologs in all three species (i.e., we only consider probes that have the human, chimpanzee, and rhesus macaque species--specific versions on the array -we refer to these as the "corresponding orthologous probes"). Following this step, we excluded genes that were represented by fewer than three corresponding orthologous probes across all species. Thus, the total number of genes included in all subsequent analyses was 17,231 (95% of genes originally included on the array). Expression estimates were obtained from [4].
To identify genes that are differentially expressed (DE) between human and chimpanzee within a tissue, we used likelihood ratio (LR) tests within the frame work of nested mixed linear models, as previously described [4]. Briefly, we estimated the maximum likelihood of the full model as well as that of a reduced (null) model, in which we assume that the expression level in human and chimpanzee is similar. We then calculated −2•(log--likelihood ratio) between the fits of the reduced and full models. We expect genes that deviate from the null (i.e., genes that are truly differentially expressed between human and chimpanzee) to have higher values of this statistic.
Correlation between the number of reactions and differential concentration In order to test correlations among the number of metabolic reactions a metabolite is involved in (commonly called "connectivity") and the level of DC between species, we first examined the distribution of the number of reactions per metabolite extracted from KEGG. As can be seen in Figure S4, the distribution is enriched with low--reaction metabolites, but has a long tail of high--connectivity metabolites. We show results using a cutoff of 20 reactions to distinguish the two groups, but results remain unchanged for multiple other cutoff choices. After splitting the metabolites to low-- and high--connectivity compounds, we plotted the distribution of likelihood ratios for differential concentration between species pairs in Figure S5.

Permutation test for difference in medians
We used a permutation test to estimate the significance of observed differences between the medians of values for different groupings of metabolites (e.g. likelihood ratio for metabolites involved in reactions controlled by differentially expressed vs. non-differentially expressed enzymes). For this purpose, we define the difference between the medians (or means) of the two groups as D. Then, we randomly divided all the values into two groups of same sizes as observed, and calculated the medians of the random groups.
This permutation was repeated 10,000 times, and each time the difference between the medians of the two randomly selected groups (Di) was recorded. The test P--value was defined as the number of times where Di ≥D, divided by 10,000.

Biosynthesis of phenylpropanoids
phenylalanine, fumaric acid, 3-phosphoglycerate, citric acid, ... Table S4. Cellular and molecular functions in the Ingenuity database enriched among the metabolites shown in Figure 2 (top 5 categories are shown).

Functional Category P--value Metabolites
Cell  Table S6. Metabolites found in Fu et al. [7] and also included in our analysis.  2 Whether the metabolite has a human--specific pattern in the CBC (Fu et al.). 3 Whether the metabolite has a human--specific pattern in the PFC. (Fu et al.) 4 Whether the metabolite was identified and included in our analysis. 5 Whether we found the metabolite to have a human--specific pattern in the liver. Figure S1. Principal component (PC) analysis of the normalized metabolite data, showing the first vs. second PC (A,C) and third vs. fourth PC (B,D) on the x and y axes, respectively. (A) Data after removal of outlier sample (HSA_56720_f_2); (B) data after removal of outlier and following correction for batch effect; (C) full data, before removal of outlier, PC1 vs. PC2; (D) full data, before removal of outlier, PC3 vs. PC4. Numbers correspond to the three batches, and colors correspond to the species (blue--human, red--chimpanzee, black--rhesus macaque). Figure S2. Heatmap of pairwise (Spearman) correlations of normalized and corrected metabolite levels between all samples. Colors represent different R 2 values as indicated in the color key at the top, with lighter colors indicate higher correlation. The dendrogram, which was used to order the rows and columns of the heatmap, which is shown at the top and left sides, was generated by clustering the pairwise correlation matrix using a Euclidean distance metric and complete agglomeration. Figure S3. Proportion of DC metabolites between each pair of species (H: human, C: chimpanzee, R: rhesus macaque).  Figure S5. Distribution of likelihood ratios for differential concentration between pairs of species, broken down into high--reaction metabolites (orange) and low--reaction metabolites (grey).  Figure S7. Similar to Figure 2D in the main text, but considering metabolites involved in the same reactions where at least one is DC between the species (and not necessarily both metabolites, as in Figure 2D).

Human-Chimpanzee
Metabolite change in same direction  Figure S10. Similar to Figure S9, showing on the y--axis the LR for enzyme differential expression between human and rhesus macaque.