Metabolomics Identifies Novel Blood Biomarkers of Pulmonary Function and COPD in the General Population

Determination of metabolomic signatures of pulmonary function and chronic obstructive pulmonary disease (COPD) in the general population could aid in identification and understanding of early disease processes. Metabolome measurements were performed on serum from 4742 individuals (2354 African-Americans and 1529 European-Americans from the Atherosclerosis Risk in Communities study and 859 Europeans from the Cooperative Health Research in the Region of Augsburg study). We examined 368 metabolites in relation to cross-sectional measures of forced expiratory volume in 1 s (FEV1), forced vital capacity (FVC), their ratio (FEV1/FVC) and COPD using multivariable regression followed by meta-analysis. At a false discovery rate of 0.05, 95 metabolites were associated with FEV1 and 100 with FVC (73 overlapping), including inverse associations with branched-chain amino acids and positive associations with glutamine. Ten metabolites were associated with FEV1/FVC and seventeen with COPD (393 cases). Enriched pathways of amino acid metabolism were identified. Associations with FEV1 and FVC were not driven by individuals with COPD. We identified novel metabolic signatures of pulmonary function and COPD in African and European ancestry populations. These may allow development of biomarkers in the general population of early disease pathogenesis, before pulmonary function has decreased to levels diagnostic for COPD.


Introduction
Pulmonary function, as assessed by forced expiratory volume in 1 s (FEV1), forced vital capacity (FVC), and the ratio of FEV1 to FVC (FEV1/FVC), reflects the physiological state of the lung. These measures are used to diagnose and monitor chronic obstructive pulmonary disease (COPD). Reduced pulmonary function in the general population correlates with systemic inflammation biomarkers [1,2] and has been linked to development of other pathologies including diabetes [3], kidney [4] and cardiovascular diseases [5,6]. Pulmonary function is influenced by genetics [7,8] and environmental factors, most notably, tobacco smoking [9]. Lower pulmonary function is a risk factor for mortality in the general population, independently of smoking and other risk factors and even among individuals with normal spirometry [10][11][12][13]. However, the underlying mechanisms for these diverse impacts of reductions in pulmonary function remain unknown.
Metabolomics systematically identifies and quantifies small molecules in blood and other biologic tissues. It has increasingly been applied to the study of chronic diseases [14][15][16][17]. Various circulating metabolites have been reported as associated with COPD in several clinically selected studies, generally with a modest number of cases [18][19][20][21][22][23][24][25][26][27][28][29][30]. Studies of COPD identified in clinical settings identify more severe disease than population-based studies. Studies of more severe disease may tend to identify markers of disease progression and treatment. The clinical utility of identifying circulating biomarkers of very early stages of COPD has been highlighted [31]. Biomarkers of the earliest stages of the disease process can enable development of new therapies and improve understanding of mechanisms by which reduced pulmonary function contributes to risk of developing chronic diseases. Well-powered studies of circulating metabolites in relation to spirometric measures in general population samples, unselected for disease, can identify reproducible biomarkers of early stages of the disease process before pulmonary function has declined to levels diagnostic of airflow obstruction and COPD. Publications on circulating metabolites in relation to pulmonary function in population cohorts are rare. A study of UK twins identified metabolites, including several amino acids, cofactors and vitamins, related to both FEV1 and FVC, as well as some unique associations for each [32,33]. No study has examined the association of metabolomic profiles with pulmonary function as well as relationships with COPD in a multi-ethnic general population sample.
We conducted meta-analyses of metabolomic analytes in relation to pulmonary function (FEV1, FVC and FEV1/FVC) in 2354 individuals of African ancestry and 2388 individuals of European ancestry from the Atherosclerosis Risk in Communities (ARIC) study and Cooperative Health Research in the Region of Augsburg (KORA). We further examined molecular pathways of specific metabolites and associations of metabolites with COPD.

Study Characteristics
Characteristics of study participants are presented in Table 1. Mean ages were very similar for the 2354 ARIC African-Americans (53.0 years, SD = 5.7), the 1,529 ARIC European-Americans (54.6, SD = 5.8) and the 859 KORA Europeans (53.8, SD = 5.7). Current smoking prevalence were also similar, ranging from 21 to 28%. A total of 393 individuals met our definition of COPD (358 from ARIC and 35 from KORA).

Metabolic Associations with Pulmonary Function Measures
Manhattan plots for the three pulmonary function parameters are shown in Figure 1. At FDR < 0.05, 95 metabolites were associated with FEV1 (30 at Bonferroni correction) and 100 (37 at Bonferroni correction) with FVC. Given the large number of findings, in the main text tables, we show only the top 20 associations (based on meta-analysis p value) that meet Bonferroni correction where the metabolite was measured in all three datasets (Table 2 for FEV1 and Table 3 for FVC); all metabolites associated at FDR < 0.05 are in Table S1 for FEV1 and Table S2 for FVC. For metabolites associated at FDR < 0.05, the average effect size per SD difference in metabolite for FEV1 was a 25.8 mL difference (range 17.3 mL-50.7 mL) and for FVC, a 29.9 mL difference (range 19.7 mL-56.9 mL). There were 23 metabolites associated, at Bonferroni significance, with both FEV1 and FVC, all with matching directions of effect. These include eight amino acids and derivatives (glycine, asparagine, tryptophan betaine, 3-(4-hydroxyphenyl)lactate, 3-phenylpropionate (hydrocinnamate), isoleucine, 2-methylbutyrylcarnitine (C5) and glutamine), three carbohydrates (fructose, lactate and mannose), the lipid glycerol, the nucleotide N2,N2-dimethylguanosine, and several gamma-glutamyl amino acids. At FDR < 0.05, 73 metabolites associated with both FEV1 and FVC, all with matching directions (Table S3).
For FEV1/FVC we identified fewer associations (Table 4): 10 at FDR < 0.05, including 5 at Bonferroni significance. On average, per SD change of those 10 FDR significant metabolites is correlated with 0.48% difference (range 0.33% to 0.93%) in the FEV1/FVC, expressed as a percentage. Among the findings are two xenobiotics with inverse associations that reflect behaviors of individuals with low FEV1/FVC. These include the respiratory medication theophylline and its metabolite, 1,3-dimethylurate. Beside theophylline treatment in asthma and COPD patients, caffeine intake has to be considered as a common source for circulating theophylline. In sensitivity analyses of the association between metabolites and pulmonary function after deleting COPD cases, associations between FEV1/FVC and theophylline and its metabolites disappeared. These sensitivity analyses further showed that results for FEV1 and FVC were not driven by individuals with COPD: correlations between beta coefficients from fixed effect meta-analyses before and after exclusion of COPD were 0.96 for FEV1 and 0.99 for FVC. methylbutyrylcarnitine (C5) and glutamine), three carbohydrates (fructose, lactate and mannose), the lipid glycerol, the nucleotide N2,N2-dimethylguanosine, and several gamma-glutamyl amino acids. At FDR < 0.05, 73 metabolites associated with both FEV1 and FVC, all with matching directions (Table S3).

Pathway Analyses
In pathway analyses, 83 out of 95 metabolites related to FEV1 and 84 of 100 related to FVC (FDR < 0.05) were matched to the pathway database and mapped to 49 and 42 pathways, respectively, for testing of enrichment. Among the four enriched pathways identified for either FEV1 or FVC, "aminoacyl-tRNA biosynthesis" and "phenylalanine metabolism" were implicated for both phenotypes (Table 5). We did not perform pathway analysis for FEV1/FVC given the small number of significant associations.

Consistency of Associations
Directions of effect were extremely consistent across the three datasets for metabolites that were statistically significant in the meta-analysis. Across all associations at FDR < 0.05 for FEV1, FVC or FEV1/FVC, only a very small proportion of associations have directions of effect that do not match for all three datasets with the meta-analysis result (14 for FEV1, 10 for FVC and 0 for FEV1/FVC); the majority of the apparently discordant results are in the smallest dataset, KORA (Table S8). For COPD all 3 discordant results are again from KORA which contributed only 9% of COPD cases (Table S8). As expected, for all instances where the direction is flipped, the estimate is very imprecise with P-values ranging between 0.16 and 0.99 and a mean of 0.66 (Table S8). Significant results in the meta-analysis reflect that the fact that larger, more precisely estimated effects in a consistent direction in the other studies are not negated by a low precision estimate in the other direction. Results were also consistent between African and European ancestry populations. We did not identify any metabolites significantly related in the meta-analysis to any of the four traits (FDR < 0.05) where there was even a nominally significant association (uncorrected P < 0.05) in African-Americans where the effect estimate was in the opposite direction.

Discussion
By combining data from over 4700 individuals across two population-based studies, encompassing both African-American and European ancestry, we identified numerous metabolites related to quantitative pulmonary function measures. Individual metabolites related to FEV1 or FVC highlighted four enriched pathways that collectively reflect associations with amino acids and their gamma-glutamyl derivatives.
There have been various studies of circulating metabolomic profiles in relation to pulmonary function or COPD, predominantly in individuals selected in clinical setting for COPD phenotypes, including exacerbations [18,28,[34][35][36]. Most included fewer than 100 individuals with COPD and several included less than 20. In clinically selected samples, it can be difficult to enroll a suitable control group truly representative of the population that gave rise to those cases. Given the large number of metabolites examined and the limited sample size of most studies, correction for multiple testing for individual metabolites was often not possible, making replication difficult. Circulating biomarkers of early COPD have been lacking and could be useful in development of interventions to limit disease progression. Examination of variations in circulating metabolites related to pulmonary function in the general population hold the potential to develop circulating markers of early stages of COPD pathogenesis that predate diagnostic reductions in pulmonary function [31]. However, studies of metabolomics profiles and quantitative pulmonary function traits across the general population are scarce. Our findings generally reflect associations across the normal range of these parameters, as they were largely similar when deleting individuals with COPD, despite the resulting smaller sample size and attenuation of the range of pulmonary function. A previous study by Menni et al. [32] with wide coverage (280 metabolites) in a population-based cohort of twins in the UK, identified at Bonferroni significance, 18 metabolites related to FEV1 and 10 related to FVC. From the combined list of 21 metabolites identified for either trait, we replicate 11 in our meta-analysis at FDR < 0.05 (3-phenylpropionate, asparagine, gamma-glutamylvaline, glycerate, glycine, indolepropionate, N-acetylglycine, pseudouridine, scyllo-inositol, succinylcarnitine, and threonate) and an additional three at nominal significance (serine, proline and glutamate). Four of their 21 metabolites were not included in our meta-analysis (butylcarnitine, C-glycosyltryptophan, CMPF-3-carboxy-4-methyl-5-propyl-2-furanpropanoate and pyroxidate). Three metabolites identified by Menni et al. did not replicate in our data (alpha-tocopherol, benzoate and bilirubin(ZZ)). However, gamma-tocopherol, which is tightly inversely correlated with alpha-tocopherol [37] and also strongly influenced by supplement use, was associated with FVC at FDR < 0.05 in the corresponding opposite direction from their alpha-tocopherol association. In addition, biliverdin, which, like bilirubin, is a product of heme breakdown, was related to both FEV1 and FVC at FDR < 0.05 in the same direction as their bilirubin association. Thus, our data provide very substantial replication of the population-based cohort findings of Menni et al. [32]. In addition to replication of previous findings, our meta-analysis identified many novel associations not reported by Menni et al. [32]. We identified at Bonferroni significance, 30 novel associations for FEV1 and 37 for FVC, including 22 that overlap (Table S3). At FDR < 0.05, we identified 111 novel metabolites related to either FEV1 or FVC (included in Tables S2  and S3). The effect sizes we observed for metabolites are not trivial in the context of the natural history of pulmonary function traits. Once maximum lung growth has been achieved, between 20-30 years of age, FEV1 and FVC inexorably decline with age. The mean effect size across significant associations, per SD difference in metabolite, of 25.8 mL, is greater than the estimated average decline in lung function among health nonsmoking adults of 18-20 mL per year [38]. It is also larger than the 10 year impact on FEV1 of exposure to secondhand tobacco smoke, estimated at 15 mL [39]. Modest differences can be meaningful because lower FEV1 and FVC at a given age increases the risk of declining below threshold levels diagnostic of COPD and lower lung function, even within the normal range is related to increased mortality, independently of other risk factors.
For FEV1/FVC, we identified far fewer associations than for either FEV1 or FVC. Identification of many more metabolites associated with FEV1 or FVC than FEV1/FVC suggests that our findings reflect coordinate effects on either preservation or deterioration in lung size as opposed to airflow obstruction. Metabolites significantly associated with FEV1 substantially overlapped with those for FVC with matching direction and similar magnitudes of effects. Accordingly, for the ratio of FEV1 to FVC, these associations generally canceled out. Although we identified many metabolites associated with FEV1 and FVC, speculation regarding mechanisms for most of these associations is difficult because the determinants of variation between individuals in the general population are not known for most metabolites. Experimental studies with manipulation of diet, physical activity, and other lifestyle factors in humans and animal models will be needed to gain a full understanding of normal inter-individual variation.
Among the interesting novel associations with FEV1 and/or FVC in our study is a cluster of metabolites previously associated with metabolic health. This cluster-containing the branched-chain amino acids (BCAA: isoleucine, leucine and valine), BCAA metabolic byproducts (2-methylbutyrylcarnitine, 3-hydroxyisobutyrate), aromatic amino acids (tyrosine, phenylalanine), glycine and glutamate/glutamine-was initially described as a metabolic signature in obese individuals that correlates with insulin resistance [40]. Subsequent studies demonstrated that this metabolic signature is prognostic for incident type 2 diabetes and predicts response to weight loss interventions [41][42][43][44]. Notably, the overrepresented pathways associated with FEV1 and FVC are enriched with these metabolites, as well. In general, higher BCAA/metabolic byproduct and aromatic amino acid concentrations, and lower glycine concentrations, reflect worse metabolic health. Thus, the inverse associations of BCAA/metabolic byproducts and aromatic amino acids, and positive association of glycine, with FEV1 and FVC, suggest that these metabolites might be reporting on the metabolic health status of individuals which, in turn, may impact pulmonary function. The glutamyl amino acids gamma-glutamylglutamine, gamma-glutamylglutamate, gamma-glutamylleucine and gamma-glutamylvaline were also associated with FEV1 and/or FVC, which likely reflects abundance of their cognate amino acids. Interestingly, a prior study of serum metabolomics in the ECLIPSE COPD cohort showed a decrease, rather than increase, in BCAA/metabolic byproducts in individuals with GOLD stage IV COPD compared with control individuals [36]. No association was seen with less severe COPD and similar findings have been seen in smaller clinically selected COPD studies [30,45]. In ECLIPSE, circulating glutamine, which largely reflects muscle turnover, was higher in COPD patients than controls [36], whereas in our general population study glutamine was positively associated with pulmonary function. Glutamine, the most abundant amino acid in the circulation, is a precursor for synthesis of glutathione, a major endogenous antioxidant, and plays important roles in response to injury and immune function [46,47]. Because patients with advanced COPD often exhibit reduced caloric intake, some degree of muscle wasting or frank cachexia, findings for BCAA and glutamine in these studies can be ascribed to the generalized state of protein deficiency and/or increased protein turnover [36].
A recent paper that examined urinary metabolites identified associations with histidine in relation to FEV1 and FEV1/FVC [48]. Although histidine was measured in our serum samples, we did not find associations with any of our phenotypes. Low lung function is associated with asthma and forms the basis for diagnosis of COPD. Thus, one may speculate, that biomarkers associated with low lung function may also be associated with lung disease. Three of the four major pathways enriched for both FEV1 and FVC in our analysis, "aminoacyl-tRNA biosynthesis" and "nitrogen metabolism" and "alanine-aspartate-glutamate metabolism", were found to be associated with asthma phenotypes in 6-14 year old children [49]. Different carnitine compounds and amino acids, e.g., tyrosine, glutamine, glycine, asparagine, phenylalanine, leucine, and valine, as well as succinate have been associated with COPD in clinical studies [25,36] and were also detected in our analysis.
A recent study of 136 individuals from the COPDGene study also found differences in amino acids and to be related to exacerbations severity and frequency [50]. We are not able to examine exacerbations. However, the fact that we identified many of these associations in relation to quantitative traits of FEV1 and FVC, and these associations were not driven by individuals with COPD in our study, and that COPD in our study was skewed toward less severe disease than COPDGene, suggests that our findings are not impacted by failure to account for exacerbations.
Identifying subtle metabolic differences related to spirometry in the general population may shed light on why variation in pulmonary function, even within the normal range, is related to development of metabolic and cardiovascular disease outcomes [3][4][5][6]. Of interest, a recent study identified alterations in metabolic activity in cardiac tissue related to FEV1/FVC [51]. With respect to early metabolic changes that may presage COPD, in an experimental model of emphysema in the mouse, circulating isoleucine and valine were increased very early in the disease process [52] consistent with the inverse associations we observed for these branched chain amino acids and pulmonary function parameters in our general population study. Some recent studies have suggested the importance of sphingolipids and sphinomyelins in subphenotypes of COPD [18,29,45]. As noted, to be included in the meta-analysis, metabolites had to have measured in at least two of the three datasets. Because the three datasets were not measured at the same calendar time, the number of metabolites reported is greater in the newer datasets. No sphingolipids or sphingomyelins were reported for the KORA, thus our meta-analysis included only two sphingomyelins, both measured in ARIC: stearoyl sphingomyelin which was inversely associated with FVC in our data (FDR < 0.05) and palmitoyl sphingomyelin which was not associated with any phenotypes. In a previous study [29], sphingomyelins were related to emphysema but not to the same phenotypes that we examined in our study.
Some previous metabolomics studies used data reduction methods such as Principal Component Analysis (PCA) and Partial Least Squares Discriminant Analysis (PLS-DA) to summarize results. These are generally smaller studies with limited power to reliably estimate associations with individual metabolites. These methods can be used when reporting results from single studies, but are not suitable for meta-analysis results. Because these methods are dataset driven, results are not transportable from one dataset to another. For example, the first principle component (PC1) from a given study is not necessarily similar to PC1 from another. These methods must be run within a single dataset and do not accommodate meta-analyzed results and thus could not be used here.
Our study has limitations and strengths. The many metabolites analyzed represent only a fraction of the whole metabolome. Further, as is commonly done in these analyses, we only included metabolites with values reported for at least 25% of each population. Measurements were not done in the three datasets at the same time points. The number of metabolites reported by the laboratory increased over time. Thus, there are some metabolites that are not available (denoted as NA in tables) in one dataset or another. Further, it is possible that a given metabolite might meet threshold for inclusion in one dataset (nonmissing in <25%) but in another dataset perhaps 26% are missing, so it is coded as NA in the table. Our meta-analyses further excluded those with measurements only in one population. Thus, we may have missed many additional metabolites related to our phenotypes. We are limited in our ability to comprehensively assess how associations compared by ethnicity because analyses were mostly done in different analytic batches. The cross-sectional measures of lung function we analyzed reflect both the maximal growth by early adulthood and the rate of decline with age; we cannot disentangle these two processes. The three population-based samples had, as expected, small proportions of COPD cases and thus we have much more limited power to find associations with this disease phenotype compared to the quantitative spirometry traits. However, the number of individuals with COPD is larger than in most previous studies. We acknowledge that because our studies, like most general population cohorts, did not have post-bronchodilator measures of pulmonary function, our COPD definition does not meet GOLD criteria. However, bronchodilator testing for COPD diagnosis, has been shown to be poorly reproducible, especially in milder disease [53], the majority of COPD in our population. We also note that major recent studies of trajectories of lung function in COPD development used pre-bronchodilator measures [54] and large scale genetic studies of COPD have also defined disease based on pre-bronchodilator measures [55]. COPD in this general population study was less severe than in previous clinically-selected metabolomic studies. This enhances our ability to identify potential biomarkers of early disease processes, as opposed to effects of severe disease, including reduced protein and caloric intake, increased muscle turnover, and cachexia. Our relatively large population-based cohort samples enabled us to examine associations with metabolomic profiles across the range seen in the general population, unselected for any diseases. Further, classifying disease status within larger population-based studies enables comparison to the optimal representative control population that gave rise to the cases. We also studied populations from two different ancestral populations: European and African ancestries. Few data are available in non-European populations, and are able to find some common metabolic signatures for lung function.
In summary, in these two population-based cohorts, that span two continents and two ancestral populations, we both confirmed a number of previously reported associations with quantitative lung function traits and identified many novel ones. Many of these are also related to COPD. Patterns of the identified metabolites overlap with some previously reported for metabolic disorders, traits that are also associated with reduced lung function. These findings may shed light on mechanisms of reduction in lung function and its links to metabolic and cardiovascular disease. In addition, identifying subtle variation in circulating metabolites related to pulmonary function in the general population may identify circulating biomarkers of the early stages of COPD pathogenesis before spirometry has reduced to levels diagnostic of disease. Such biomarkers can help discover potential therapeutic interventions.

Study Populations
The Atherosclerosis Risk in Communities (ARIC) study is an epidemiological study designed to investigate the etiology of cardiovascular disease in European and African Americans. ARIC enrolled 15,792 individuals aged 45-64 years from four U.S. communities (Forsyth County, NC; Jackson, MS; suburbs of Minneapolis, MN; and Washington County, MD) at baseline (1987)(1988)(1989) [56]. Metabolite levels were measured in fasting serum samples collected at the baseline examination among a random sample of African and European Americans. Pulmonary function was measured at this same visit.
The KORA F4 study is a population-based cohort from Augsburg city and the two adjacent counties in southern Germany (KORA-Cooperative Health Research in the Region of Augsburg). Between 2006 and 2008, 3080 individuals aged 32 to 81 years participated in KORA F4, the seven-year follow-up of the KORA S4 survey. Details of the KORA platform and the study design have been described previously [57,58]. Pulmonary function was measured on a random sample of 1321 subjects aged 41 to 62 years and metabolite levels were measured in fasting blood samples collected at this visit.

Assessment of Pulmonary Function and Covariates
In ARIC, spirometry was conducted with a Collins Survey II water-seal spirometer (Collins Medical, Inc., Braintree, MA, USA) [62]. In KORA F4, spirometry was measured using a pneumotachograph-type spirometer (Masterscope PC, CareFusion/VIASYS Healthcare/Erich Jaeger, Höchberg, Germany). Height and weight were measured at the time of spirometry. Information on smoking was assessed by questionnaires. For both ARIC and KORA, predicted values were calculated using reference equations from the Global Lung Function Initiative [63]. COPD was defined using prebronchodilator spirometry applied to the Global Initiative for Chronic Obstructive Lung Disease (GOLD) criteria for diagnosis of COPD of at least moderate severity (FEV1 < 80% predicted and FEV1/FVC < 0.7) [55]. Estimated glomerular filtration rate (eGFR CKD-EPI ) was calculated based on measured creatinine levels [64].

Human Subjects
ARIC was approved by the Institutional Review Board at the University of Minnesota, Johns Hopkins University, Wake Forest University, University of North Carolina, University of Texas Health Sciences Center at Houston, and University of Mississippi Medical Center. KORA was approved by the Ethics Committee of the Bavarian Medical Association (Bayerische Landesärztekammer) and the Bavarian commissioner for data protection and privacy (Bayerischer Datenschutzbeauftragter). Written consent was obtained from all ARIC and KORA study participants.

Statistical Analyses
In ARIC and KORA, metabolites with values reported for at least 25% of the study participants were analyzed. Metabolite missing values were imputed with half of the lowest value of the corresponding metabolite in each study by wave, and were standardized with mean at 0 and standard deviation at one (Methods in Supplementary Materials). The total number of metabolites included in the analyses was 721 for ARIC European Americans (measured in 2014), 355 for ARIC African Americans (measured in 2010) and 266 for KORA (measured 2009-2010). The reason for the differing number of metabolites in the three datasets is improvement in both instrumentation and the size of the authentic standard library over time.
We analyzed data from each of the 4742 study participants (1529 ARIC EA, 2354 ARIC AA and 859 KORA) with metabolite measurements and complete data on lung function and covariates. Linear regression models were fitted for each standardized metabolite and lung function parameter separately in each cohort, adjusting for the factors traditionally included in pulmonary function prediction equations (age, age 2 , sex, height, height 2 ) [65], weight, smoking status (never, former, current as two categorical variables), lifetime cigarette pack years, cigarettes currently smoked per day, and the eGFR CKD-EPI [64]. For the 368 metabolites detected in at least two of the three populations, study-specific associations with pulmonary outcomes were meta-analyzed using fixed effect inverse variance weighting [66]. Consequently, metabolites only available in one of the cohorts were excluded from the analyses. For correction for multiple testing we applied false discovery rate (FDR) of 0.05 using the Benjamini-Hochberg procedure and the stricter Bonferroni corrected p < 1.4 × 10 −4 based on 368 tests, which assumes tests are independent. All analyses were performed using R, version 3.4.1 [67].
For metabolites associated with pulmonary function parameters at FDR < 0.05, we performed metabolic pathway analysis using MetaboAnalyst 3.0 [68] to examine if a particular pathway was represented more than expected compared to metabolic pathways based on the Human Metabolome Database 3.0 [69]. Statistical significance of enriched pathways was assessed by FDR < 0.05 using the hypergeometric test.