Sex-Specific Relationship between the Cardiorespiratory Fitness and Plasma Metabolite Patterns in Healthy Humans—Results of the KarMeN Study

Cardiorespiratory fitness (CRF) represents a strong predictor of all-cause mortality and is strongly influenced by regular physical activity (PA). However, the biological mechanisms involved in the body’s adaptation to PA remain to be fully elucidated. The aim of this study was to systematically examine the relationship between CRF and plasma metabolite patterns in 252 healthy adults from the cross-sectional Karlsruhe Metabolomics and Nutrition (KarMeN) study. CRF was determined by measuring the peak oxygen uptake during incremental exercise. Fasting plasma samples were analyzed by nuclear magnetic resonance spectroscopy and mass spectrometry coupled to one- or two-dimensional gas chromatography or liquid chromatography. Based on this multi-platform metabolomics approach, 427 plasma analytes were detected. Bi- and multivariate association analyses, adjusted for age and menopausal status, showed that CRF was linked to specific sets of metabolites primarily indicative of lipid metabolism. However, CRF-related metabolite patterns largely differed between sexes. While several phosphatidylcholines were linked to CRF in females, single lyso-phosphatidylcholines and sphingomyelins were associated with CRF in males. When controlling for further assessed clinical and phenotypical parameters, sex-specific CRF tended to be correlated with a smaller number of metabolites linked to lipid, amino acid, or xenobiotics-related metabolism. Interestingly, sex-specific CRF explanation models could be improved when including selected plasma analytes in addition to clinical and phenotypical variables. In summary, this study revealed sex-related differences in CRF-associated plasma metabolite patterns and proved known associations between CRF and risk factors for cardiometabolic diseases such as fat mass, visceral adipose tissue mass, or blood triglycerides in metabolically healthy individuals. Our findings indicate that covariates like sex and, especially, body composition have to be considered when studying blood metabolic markers related to CRF.


Introduction
Cardiorespiratory fitness (CRF) is a health-related component of physical fitness (PF) [1], reflecting the ability of the circulatory, respiratory, and muscular systems to take up, transport, and utilize oxygen during sustained physical exercise (PE) [2]. It is affected by sex [3,4], age [4,5], lean body mass (LBM) [6], heredity [7], and behavioral factors like diet, smoking, or physical activity (PA) [4]. Actually, CRF is not only an objective measure of regular PA [8] but has also emerged as a strong predictor of all-cause and diseasespecific mortality [9]. Representing the main modifiable determinant of CRF, regular PA is known to favorably influence body composition and glucose-insulin homeostasis, as well as the lipoprotein profile [8]. At the muscular level, beneficial adaptations to PA comprise an increased capillarization, a higher mitochondrial density, and enhanced oxidative metabolism, finally leading to an improved endurance capacity [10]. However, despite profound knowledge on the health-promoting benefits of PA, the molecular mechanisms and metabolic pathways involved in the whole-body and skeletal muscle adaptation to PA are still insufficiently understood [11].
The emerging field of metabolomics is a promising approach to systematically investigating exercise-induced changes in human metabolism related to performance and health [12]. By applying nuclear magnetic resonance (NMR) spectroscopy-or mass spectrometry (MS)-based techniques, metabolomics permits the simultaneous analysis of a high number and variety of metabolites, i.e., low-molecular-weight compounds, that represent the end-products of interactions between genes, proteins, and the cellular environment [13]. Thus, metabolomics can help to identify PA-or PF-associated metabolite profiles, possibly hinting at metabolic pathways that are linked to the well-known effects of exercise [12].
Interestingly, recent metabolomic studies have provided the first evidence that higher levels of PA or PF are linked to lower circulating branched-chain amino acid (AA) [14][15][16][17] and higher circulating phosphatidylcholine (PC) concentrations [18][19][20][21]. However, research on the relationship between CRF and the blood metabolome in large populations, including both sexes and with a broad age spectrum, is rather scarce. In fact, only half of the studies that assessed maximal oxygen uptake (VO 2max ) as the gold standard of aerobic fitness conducted correlation or regression analyses [17,19,20,22,23], while the other half examined differences between groups with high or low CRF [15,18,21,24]. Limitations of the former studies are that the results were restricted to young [17,22], middle-aged [19,20], or older [23] individuals and are thus hardly transferable to the general population. Besides, studies either had rather small sample sizes [22,23] or solely included male subjects [17]. Apart from Lustgarten et al., who detected nearly 300 serum analytes [22], the remaining studies focused on a limited number of metabolites. Since several CRF-associated metabolites that were reported in the literature have also been linked to other phenotypical variables such as body composition, adjustments for potential confounders are decisive to determine if correlations can be specifically attributed to CRF [25].
As a way of overcoming those limitations, we applied a multi-platform metabolomics approach and exploratory bi-and multivariate statistical procedures to systematically analyze the relationship between CRF and 427 plasma metabolites in 252 healthy women and men from the cross-sectional Karlsruhe Metabolomics and Nutrition (KarMeN) study. Participants had a wide age range and were thoroughly characterized by anthropometric, functional, and clinical examinations [26]. Therefore, we were able to take a variety of known and potential confounders into account. Since it has already been shown that sex, age, and menopausal status are determinants of CRF [4,27] and are also linked to a discriminatory plasma metabolite profile in the KarMeN participants [28], all analyses were conducted in sex-specific subgroups and adjusted for age and menopausal status. Firstly, both bi-and multivariate associations between CRF and metabolites were calculated, using bivariate correlation or partial least squares (PLS) regression analyses, respectively. Secondly, to identify associations that were independent of other phenotypical or clinical variables, correlation analyses and PLS models were additionally adjusted for parameters related to body composition, clinical blood biochemistry, lung and arterial function, short-term and habitual PA, or diet. Thirdly, cross-validated stepwise regression procedures were conducted, thus selecting sets of phenotypical, clinical, and plasma metabolite variables that contribute to a preferably good explanation of CRF.

Metabolomics Data
427 plasma analytes were included in the final data analysis. Untargeted methods yielded 234 analytes, of which 43 (18.4%) could be identified with sufficient certainty. 193 analytes were derived from targeted analyses and were thus known a priori. Of the 236 identified metabolites, most belonged to lipid metabolism (69.5%) and AA metabolism (17.4%), followed by xenobiotics-related metabolism (5.1%), mammalian-microbial cometabolism (3.0%), carbohydrate or energy metabolism (each 2.1%), and nucleotide or cofactors and vitamins metabolism (each 0.4%). The classification of the identified metabolites to both major and specific metabolic pathways is provided in Table S1.

Basic Characteristics of Study Participants
The study sample consisted of 252 healthy individuals, 150 males (M, 59.5%) and 102 females (F, 40.5%), with a mean age of 45.9 ± 17.1 years and a mean peak oxygen uptake (VO 2peak ) of 38.9 ± 11.7 mL kg −1 min −1 . The characteristics of participants are presented according to sex-specific VO 2peak quarters ( Figure 1). In the radar plots, the respective means of the lowest VO 2peak quarter (1st q) were used as a reference value and the means of the other quarters (2nd q, 3rd q and 4th q) were related to the means of the first quarter. The absolute means, standard deviations, and respective units are provided in File S1.
Metabolites 2021, 11, x FOR PEER REVIEW 3 of 25 related to body composition, clinical blood biochemistry, lung and arterial function, shortterm and habitual PA, or diet. Thirdly, cross-validated stepwise regression procedures were conducted, thus selecting sets of phenotypical, clinical, and plasma metabolite variables that contribute to a preferably good explanation of CRF.

Metabolomics Data
427 plasma analytes were included in the final data analysis. Untargeted methods yielded 234 analytes, of which 43 (18.4%) could be identified with sufficient certainty. 193 analytes were derived from targeted analyses and were thus known a priori. Of the 236 identified metabolites, most belonged to lipid metabolism (69.5%) and AA metabolism (17.4%), followed by xenobiotics-related metabolism (5.1%), mammalian-microbial cometabolism (3.0%), carbohydrate or energy metabolism (each 2.1%), and nucleotide or cofactors and vitamins metabolism (each 0.4%). The classification of the identified metabolites to both major and specific metabolic pathways is provided in Table S1.

Basic Characteristics of Study Participants
The study sample consisted of 252 healthy individuals, 150 males (M, 59.5%) and 102 females (F, 40.5%), with a mean age of 45.9 ± 17.1 years and a mean peak oxygen uptake (VO2peak) of 38.9 ± 11.7 mL kg −1 min −1 . The characteristics of participants are presented according to sex-specific VO2peak quarters ( Figure 1). In the radar plots, the respective means of the lowest VO2peak quarter (1st q) were used as a reference value and the means of the other quarters (2nd q, 3rd q and 4th q) were related to the means of the first quarter. The absolute means, standard deviations, and respective units are provided in File S1. Radar plots visualizing the basic characteristics of KarMeN participants according to sex-specific VO2peak quarters. The means of the 1st q were used as reference values to the means of the 2nd, 3rd and 4th q. * : significant differences between quarters according to the Welch ANOVA. °: n = 23 (1st q), n = 25 (2nd q); ⱽ: n = 36 (1st q); ⱽⱽ: n = 35 (1st q), n = 37 (2nd q). AEE: activity energy expenditure; BMC: bone mineral content; BMI: body mass index; BP: blood pressure; dia: diastolic; FEV1: forced expiratory pressure in one second; FM: fat mass; Hb: hemoglobin; HDL: high-density lipoprotein; HEI-NVS: Healthy Eating Index (modified version); HRrest: resting heart rate; LBM: lean body mass; LDL: lowdensity lipoprotein; MET: metabolic equivalent of task; PWV: pulse wave velocity; q: quarter; sys: systolic; TGs: triglycerides; VATM: visceral adipose tissue mass; VCmax: maximal vital capacity; VO2peak: peak oxygen uptake.
In both females and males, there were statistically significant differences between VO 2peak -related quarters with regard to age, weight, body mass index (BMI), fat mass (FM (%)), and visceral adipose tissue mass (VATM (in kg)). Furthermore, differences for clinical parameters like fasting blood glucose, HbA1c, triglycerides (TGs), high-density lipoprotein (HDL) and low-density lipoprotein (LDL) cholesterol, systolic and diastolic blood pressure (BP), pulse wave velocity (PWV), maximal vital capacity (VC max ), and forced expiratory pressure in one second (FEV1) were observed across the quarters of the sex-specific VO 2peak . In females, additional differences between VO 2peak -related quarters could be observed for height and the total metabolic equivalent of task (MET), while in males, differences were detected for LBM, resting heart rate (HR rest ), and the activity energy expenditure (AEE). The menopausal status in women also significantly differed between VO 2peak quarters, with increasing ratios of pre-to post-menopausal women from the 1st q to the 4th q (see File S1 (Table 1)). As the non-modifiable factors age and menopausal status were associated with the sex-specific VO 2peak , and since they have already been shown to determine plasma metabolite patterns in KarMeN subjects [28], these variables were treated as confounders in all subsequent analyses.

Sex-Specific Relationship between CRF and Phenotypical/Clinical Variables
To examine sex-specific relations between the VO 2peak and selected phenotypical as well as clinical variables, correlations adjusted for age (and menopausal status in females) were calculated. A visual comparison of the pairwise correlations in women and men is provided in Figure 2. After adjustments for the above-mentioned confounding factors, the VO 2peak in females showed positive correlations with HDL cholesterol (r = 0.43) and negative correlations with the FM (%) (r = −0.61), VATM (r = −0.44), PWV (r = −0.38), and TGs (r = −0.33), that were all significantly different from zero. In males, not only HDL cholesterol (r = 0.29) but also the AEE (r = 0.20) showed positive correlations with the VO 2peak . Compared to females, negative correlations to the VO 2peak that were significantly different from zero were not only detected for the FM (%) (r = −0.62), VATM (r = −0.57), PWV (r = −0.32), and TGs (r = −0.25) but additionally for HR rest (r = −0.30), diastolic BP (r = −0.27), LDL cholesterol (r = −0.22) and insulin (r = −0.18).
Metabolites 2021, 11, x FOR PEER REVIEW 4 of 25 (%)), and visceral adipose tissue mass (VATM (in kg)). Furthermore, differences for clinical parameters like fasting blood glucose, HbA1c, triglycerides (TGs), high-density lipoprotein (HDL) and low-density lipoprotein (LDL) cholesterol, systolic and diastolic blood pressure (BP), pulse wave velocity (PWV), maximal vital capacity (VCmax), and forced expiratory pressure in one second (FEV1) were observed across the quarters of the sex-specific VO2peak. In females, additional differences between VO2peak-related quarters could be observed for height and the total metabolic equivalent of task (MET), while in males, differences were detected for LBM, resting heart rate (HRrest), and the activity energy expenditure (AEE). The menopausal status in women also significantly differed between VO2peak quarters, with increasing ratios of pre-to post-menopausal women from the 1st q to the 4th q (see File S1 (Table S1)). As the non-modifiable factors age and menopausal status were associated with the sex-specific VO2peak, and since they have already been shown to determine plasma metabolite patterns in KarMeN subjects [28], these variables were treated as confounders in all subsequent analyses.

Sex-Specific Relationship between CRF and Phenotypical/Clinical Variables
To examine sex-specific relations between the VO2peak and selected phenotypical as well as clinical variables, correlations adjusted for age (and menopausal status in females) were calculated. A visual comparison of the pairwise correlations in women and men is provided in Figure 2   Confounder-adjusted sex-specific correlations between the VO 2peak and phenotypical/clinical variables. * Pearson correlations were performed on Van der Waerden (VdW)-transformed data adjusted for age (and menopausal status in females) and respective correlation coefficients (r; dots) and 95% confidence intervals (CIs; bars) are illustrated. AEE: activity energy expenditure; BMC: bone mineral content; BP: blood pressure; FEV1: forced expiratory pressure in one second; FM: fat mass; Hb: hemoglobin; HDL: high-density lipoprotein; HEI-NVS: Healthy Eating Index (modified version); HR rest : resting heart rate; LBM: lean body mass; LDL: low-density lipoprotein; MET: metabolic equivalent of task; PWV: pulse wave velocity; TGs: triglycerides; VATM: visceral adipose tissue mass; VC max : maximal vital capacity.
After additionally adjusting for phenotypical and clinical variables, 59 metabolites (0.20 ≤ |r| ≤ 0.39) in females and 24 metabolites (0.16 ≤ |r| ≤ 0.25) in males still exhibited weak to moderate correlations with the VO 2peak that were significantly different from zero. The majority of the VO 2peak -related plasma metabolites in both females and males belonged to lipid metabolism, followed by AA metabolism and xenobiotics and related metabolism. However, only a few correlations with the same directions in both sexes were detected (e.g., U1.156). The top 10 of sex-specific positive and negative partial correlations between the VO 2peak and plasma metabolites are summarized in Table 2. Table 2. Top 10 of sex-specific partial correlations between the VO 2peak and plasma metabolites.

Multivariate Association Analyses
The multivariate association between the VO 2peak and all 427 plasma analytes was assessed based on rank products obtained in cross-validated PLS models. Rank products were calculated by the geometric mean of the ranks of the regression coefficients of each metabolite in PLS models, across 20 random splits. For each metabolite variable, the importance of its contribution to the multivariate association was determined by permutation tests (see Section 4.7 (ii) for more details). A graphical overview of all sex-specific associations is provided in File S2. In Table S2, metabolites are listed together with their rank products of contribution to multivariate associations. Regarding females, PLS analysis showed that metabolites highly contributing to the confounder-adjusted multivariate association with CRF included several diacyl-and acyl-alkyl-PCs, the LCFA C16:1 9cis, as well as a number of unknowns. In males, metabolites with high contributions to the multivariate association comprised specific lysoPCs (C18:2, C18:1, C17:0), SMs (C18:0, C18:1), and diacyl-PC C40:6 next to unidentified analytes, when adjusting for age. As opposed to confounder-adjusted multivariate association analyses, the mean of the root mean square errors (RMSEs) based on the test samples was not always higher in the permutations than by using the original data, if additionally controlling for phenotypical and clinical variables (File S4). Consequently, the multivariate relationship between CRF and all 427 plasma analytes lost importance after applying additional adjustments.
To visualize the sex-specific association patterns of plasma metabolites, the results of both bi-and multivariate association analyses were combined in a volcano plot (see Figure 3 for confounder-adjusted findings (*) and File S5 for confounder-and phenotypical/clinical variables-adjusted findings (**)). In the upper right and left corners, metabolites with moderate (|r| ≥ 0.25) bivariate correlations and significant contributions to multivariate associations are detectable. Metabolites in the upper middle region showed weak (|r| ≤ 0.25) bivariate correlations but significant contributions to multivariate associations, i.e., their relationship with the VO 2peak depended on all other considered analytes. In contrast, the lower right and left corners include metabolites with moderate bivariate correlations but no relevant contributions to multivariate associations, i.e., their relationship with the VO 2peak lost relevance if other metabolite variables with possibly redundant information were taken into account. For a subsequent metabolic interpretation of CRF-related metabolite patterns, plasma metabolites with either relevant bivariate correlations (|r| ≥ 0.25) or significant contributions to multivariate associations were considered ( Figure 4). . Volcano plots illustrating sex-specific plasma metabolite patterns associated with the VO2peak. F: females; M: males. * confounder (age/menopausal status)-adjusted. The y-axis represents the significance of the contribution of each metabolite variable to the multivariate association with the VO2peak, expressed as the negative logarithm of the relative frequencies of permutationobtained rank products below measured rank products. The x-axis illustrates the direction and strength of partial correlations between the VO2peak and metabolite variables, expressed as Pearson correlation coefficients (r) of Van der Waerden (VdW)transformed variables. The classification of plasma analytes to major metabolic pathways is color-coded as follows: amino acid metabolism (dark blue); carbohydrate metabolism (yellow); cofactors and vitamins metabolism (dark green); energy metabolism (light blue); lipid metabolism (brown); mammalian-microbial cometabolism (orange); nucleotide metabolism (purple); xenobiotics and related metabolism (light green); unknown (black). . Volcano plots illustrating sex-specific plasma metabolite patterns associated with the VO 2peak . F: females; M: males. * confounder (age/menopausal status)-adjusted. The y-axis represents the significance of the contribution of each metabolite variable to the multivariate association with the VO 2peak , expressed as the negative logarithm of the relative frequencies of permutation-obtained rank products below measured rank products. The x-axis illustrates the direction and strength of partial correlations between the VO 2peak and metabolite variables, expressed as Pearson correlation coefficients (r) of Van der Waerden (VdW)-transformed variables. The classification of plasma analytes to major metabolic pathways is color-coded as follows: amino acid metabolism (dark blue); carbohydrate metabolism (yellow); cofactors and vitamins metabolism (dark green); energy metabolism (light blue); lipid metabolism (brown); mammalian-microbial cometabolism (orange); nucleotide metabolism (purple); xenobiotics and related metabolism (light green); unknown (black).

Figure 4.
Classification of relevant CRF-associated plasma metabolites to major and specific metabolic pathways. Top: females; bottom: males. Major metabolic pathways are color-coded as follows: amino acid metabolism (dark blue); carbohydrate metabolism (yellow); energy metabolism (light blue); lipid metabolism (brown). * findings from confounder (age/menopausal status)-adjusted bi-/multivariate association analyses; ** findings from bivariate association analyses additionally adjusted for phenotypical/clinical variables. CRF: cardiorespiratory fitness; PC: phosphatidylcholine; SM: sphingomyelin. Metabolites with bivariate correlations to CRF that were significantly different from zero in both females and males are indicated in bold.

Multiple Regression Analyses
To additionally investigate the sex-specific relationship between the VO2peak and selected metabolite variables in the presence of phenotypical and clinical data, multiple linear regression analyses were performed. By adjusting all included variables for age and menopausal status, the examined associations were independent of these non-modifiable variables. As described in the methods section, sex-specific models were calculated based on three different sets of phenotypical, clinical, and metabolite variables. The results of the cross-validated stepwise regression analyses are presented in Table S3.
Finally, the most suitable combination of phenotypical, clinical, and metabolite variables for a preferably good explanation of the VO2peak was selected. The included variables ("approach 1-3 selection") and the adjusted coefficients of determination (R 2 (adjusted)) for the evaluation of the sex-specific final models are summarized in Table 3. With regard to approach 1 selection, seven phenotypical or clinical variables were selected for the final model of females, resulting in an R² (adjusted) of 0.40. For the final model of males, six phenotypical or clinical variables were selected, leading to an R 2 (adjusted) of 0.43. When including all 21 phenotypical and clinical variables, an R 2 (adjusted) of 0.36 for women and an R 2 (adjusted) of 0.39 for men was obtained (see Table S3). With approach 2, the added value of plasma metabolites for explaining the confounder-adjusted VO2peak in the presence of all 21 phenotypical and clinical variables was considered.
Nine or ten metabolites, respectively, were additionally selected for the final models, leading to a comparatively higher performance in both sexes (F: R 2 (adjusted) = 0.72; M: . Classification of relevant CRF-associated plasma metabolites to major and specific metabolic pathways. Top: females; bottom: males. Major metabolic pathways are color-coded as follows: amino acid metabolism (dark blue); carbohydrate metabolism (yellow); energy metabolism (light blue); lipid metabolism (brown). * findings from confounder (age/menopausal status)-adjusted bi-/multivariate association analyses; ** findings from bivariate association analyses additionally adjusted for phenotypical/clinical variables. CRF: cardiorespiratory fitness; PC: phosphatidylcholine; SM: sphingomyelin. Metabolites with bivariate correlations to CRF that were significantly different from zero in both females and males are indicated in bold.

Multiple Regression Analyses
To additionally investigate the sex-specific relationship between the VO 2peak and selected metabolite variables in the presence of phenotypical and clinical data, multiple linear regression analyses were performed. By adjusting all included variables for age and menopausal status, the examined associations were independent of these non-modifiable variables. As described in the methods section, sex-specific models were calculated based on three different sets of phenotypical, clinical, and metabolite variables. The results of the cross-validated stepwise regression analyses are presented in Table S3.
Finally, the most suitable combination of phenotypical, clinical, and metabolite variables for a preferably good explanation of the VO 2peak was selected. The included variables ("approach 1-3 selection") and the adjusted coefficients of determination (R 2 (adjusted)) for the evaluation of the sex-specific final models are summarized in Table 3. With regard to approach 1 selection, seven phenotypical or clinical variables were selected for the final model of females, resulting in an R 2 (adjusted) of 0.40. For the final model of males, six phenotypical or clinical variables were selected, leading to an R 2 (adjusted) of 0.43. When including all 21 phenotypical and clinical variables, an R 2 (adjusted) of 0.36 for women and an R 2 (adjusted) of 0.39 for men was obtained (see Table S3). With approach 2, the added value of plasma metabolites for explaining the confounder-adjusted VO 2peak in the presence of all 21 phenotypical and clinical variables was considered. Nine or ten metabolites, respectively, were additionally selected for the final models, leading to a comparatively higher performance in both sexes (F: R 2 (adjusted) = 0.72; M: R 2 (adjusted) = 0.62). Contrary to approach 2, both phenotypical and clinical parameters, as well as metabolite variables, entered the model in a competing stepwise manner in approach 3. Actually, FM (%) was the only phenotypical/clinical variable being included in the final models of both sexes. In females, nine plasma analytes were additionally selected, resulting in an R 2 (adjusted) of 0.68. In males, eight plasma analytes completed the model which finally showed an R 2 (adjusted) of 0.59. In summary, approach 2, as well as approach 3, selection demonstrated an improved performance in both sexes, compared to the initial models solely based on phenotypical/clinical variables. While the acyl-alkyl-PC C40:3 was present in the females' final models for both approaches 2 and 3, diacyl-PC C36:3, malic acid, as well as glutamate, were selected for the final models of males.

Discussion
The major finding of our systematic association analyses is that the VO 2peak was related to sex-specific sets of plasma metabolites that primarily belong to lipid metabolism. However, the observed correlations were rather moderate, and, independently of other clinical or phenotypical variables considered, only a small number of metabolites were significantly correlated with CRF. Multiple regression analyses revealed that models explaining the sex-specific VO 2peak could be improved when including selected plasma metabolites in addition to clinical and phenotypical parameters. For a metabolic interpretation of CRF-related metabolite patterns, a graphical overview of identified metabolites with relevant bivariate correlations (|r| ≥ 0.25) or significant contributions to multivariate associations with CRF and their pathway classification is provided in Figure 4.
Apart from detecting bi-and multivariate associations between CRF and plasma metabolites, which are discussed in the following sub-sections, we were also able to prove well-known relationships between CRF and several health-related clinical or phenotypical variables [4] in the KarMeN population. After correcting for age and menopausal status, we confirmed previous studies showing that CRF correlated negatively with the FM (%) [29], VATM [30], PWV [31], and TGs [32,33] and positively with HDL cholesterol [32,33] in both females and males. Equally consistent with the literature, but only present in men, were negative correlations of the age-adjusted VO 2peak with the HR rest [5], diastolic BP [29,32], LDL cholesterol [34], and insulin [5]. In summary, our data provide evidence that even in a study sample consisting of metabolically healthy individuals and independent of age, individuals with higher CRF generally demonstrated lower values of clinical parameters, some of which are recognized as traditional risk factors for chronic metabolic or cardiovascular diseases.
3.1. Sex-Specific Plasma Metabolite Patterns Related to CRF 3.1.1. Age-and Menopausal Status-Adjusted Findings As shown by the results of age-and menopausal status-adjusted (*) correlation and PLS regression analyses, the CRF-related plasma metabolite pattern in females mainly comprised PCs, all of which were individually positively correlated with the VO 2peak . PCs represent the most abundant phospholipid component in cellular membranes and plasma lipoprotein classes, being important for cell integrity or the assembly and stability of lipoproteins [35]. Besides, acyl-alkyl-PCs were proposed as antioxidants preventing lipoprotein oxidation [36]. Previous studies have reported lower levels of acyl-alkyl-PCs in obese or insulin-resistant individuals [36,37] and showed that specific acyl-alkyl-PCs were related to higher HDL cholesterol and a lower risk for type 2 diabetes [38]. In line with our results, Wientzek et al. demonstrated age-and sex-adjusted positive associations between CRF and several serum PCs in middle-aged adults [20]. Likewise, higher plasma levels of four acyl-alkyl-PCs were observed in adults with high CRF compared to less fit adults, when controlling for age and BMI [21]. While those relationships seemed to be independent of sex, associations between CRF and specific PCs in the KarMeN population were more pronounced in women than in men.
With regard to men, age-adjusted association analyses revealed that the CRF-related plasma metabolite pattern was dominated by three lysoPCs (C18:2, C18:1, C17:0), all of which were positively linked to the VO 2peak , and two SMs (C18:0, C18:1), which were negatively linked to the VO 2peak . In fact, lysoPC C18:2 also showed a relevant positive correlation with CRF in females. LysoPCs represent hydrolysis products from PCs, with relevant roles for cell signaling. As a major component of oxidized LDL, lysoPCs are also supposed to regulate the pathophysiological processes underlying atherosclerosis [39]. It is noteworthy that saturated lysoPCs are assumed to exert pro-inflammatory effects, whereas polyunsaturated lysoPCs such as C18:2 do not seem to possess inflammatory properties [40]. Actually, circulating lysoPCs were found to be reduced in obese individuals [41], and in particular, lysoPC C18:2 has been linked to a lower risk for type 2 diabetes [38] or cardiovascular disease [42]. As is consistent with our findings, Wientzek et al. showed positive correlations between CRF and serum lysoPCs C18:2 and C18:1, after controlling for sex and age. Thus, it is possible that lysoPC C18:2, in particular, might provide a potential link between CRF and its protective effects on chronic diseases. Our data moreover suggest a sex-dependent regulation of phospholipid metabolism in relation to CRF status. Although PC hydrolysis and lysoPC formation have been shown to be generally higher in men than in women (possibly due to differences in enzymatic activities, body composition, and/or hormonal or lifestyle factors [43]), the exact mechanisms underlying sex-related differences in CRF-associated phospholipids are largely speculative. However, we assume sex to be an important factor when studying how PC metabolism is linked to the health-beneficial effects of high CRF. In addition to lysoPCs, SMs are present in cell membranes or linked to lipoproteins [44]. Higher plasma SMs have been proposed as independent risk factors for cardiovascular diseases [45]. In particular, SMs with saturated acyl chains (C18:0 to C24:0) were closely correlated with parameters of obesity or insulin resistance [46]. As is consistent with our results, previous studies reported negative correlations between CRF and blood SM C18:0 [22,47] or SM C18:1 [20,47] in young [22] and middle-aged [20] adults or patients with coronary artery disease [47]. Despite sex-related differences, our findings indicate that even in healthy individuals and across a broad age range, higher CRF tends to be associated with lower values of potential novel blood biomarkers of the pathophysiological processes underlying cardiometabolic diseases.
Further relevant CRF-related plasma metabolites in females were two acylcarnitines (C14:2, C10:2), the LCFA C24:0, glyceric acid, acetate, and citric acid, all of which were positively linked to the VO 2peak , and two LCFAs (C16:1 9cis, C18:1 11cis), which showed negative correlations with the VO 2peak . Of those metabolites, palmitoleic acid C16:1 9cis also significantly contributed to the multivariate association with CRF. C16:1 9cis is an abundant fatty acid in human blood and adipose tissue. It can be ingested through diet or endogenously produced, and is assumed to act as a beneficial lipokine that prevents the negative effects of adiposity on insulin sensitivity [48]. Since circulating C16:1 9cis has been shown to be proportional to FM [48], the negative correlation between CRF and C16:1 9cis might be explained by the generally lower FM (%) in fitter females. Acylcarnitines are intermediates in the transport of LCFAs to mitochondria or byproducts of β-oxidation [49]. Although blood acylcarnitines have been identified as markers of insulin resistance, mitochondrial overload, and incomplete fat oxidation [50], they are also physiologically elevated in conditions with high lipolytic rates [51]. Recently, it has been shown that a training-induced rise in fasting levels of muscular long-and medium-chain acylcarnitines (e.g., C14:2 and C10:2) were related to improved CRF and potentially reflective of a more robust carnitine buffering system [52]. Glyceric acid is a sugar acid that connects several pathways, e.g., glycerolipid metabolism and glycolysis/gluconeogenesis. Even if its biological relevance with regard to PE has, to our knowledge, not yet been described, Lustgarten et al. also demonstrated a positive association between CRF and circulating glyceric acid in healthy females [22]. Equally in line with our results, blood acetate was shown to be higher in physically active adults [14]. Acetate is either directly formed from pyruvate, providing a source for acetyl-coenzyme A [53], or can be produced by the intestinal microbiota, finally entering circulation [54]. It has been suggested that a higher PF is linked to a greater abundance of gut bacteria with positive health effects, being reflected in the release of fermentation metabolites like acetate [55,56]. In addition to its anti-inflammatory and vasodilatory properties [54], acetate has been proposed as an important energy substrate during endurance exercise in mice [57]. However, the extent to which the microbiome indeed affects the PE capacity in humans, and whether resting blood acetate might mirror this association, requires further investigation. While the tricarboxylic acid (TCA) cycle intermediate citric acid only showed a relevant positive correlation with the VO 2peak in females, malic acid and succinic acid tended to be positively linked to CRF in both sexes. Previous studies revealed a rise in fasting plasma malic acid after weight loss and PE intervention in obese women [58], a training-induced increase in muscular succinic acid in subjects at risk of metabolic disease [52], or a slightly positive correlation between CRF and serum succinic acid in healthy young men [59]. Despite weak bivariate correlations, our results support the suggestion that specific TCA cycle intermediates might be potentially interesting blood markers of the beneficial effects of chronic PA occurring at a muscular level, such as an increased mitochondrial density or TCA cycle capacity.
To conclude, we provided evidence of sex-specific CRF-associated plasma metabolite patterns after adjusting for age and menopausal status. Sex-related differences were especially observed for lipid metabolism-related PCs, which generally showed relevant associations with the VO 2peak in females but not in males. However, weak to moderate correlations between CRF and lysoPC C18:2, LCFA C16:1 9cis, glyceric acid, as well as the energy or carbohydrate metabolism-related succinic acid, malic acid, or acetate, were observed in both sexes. Accordingly, our results suggest that CRF-related adaptations in glycerophospholipid metabolism might vary between sexes, whereas the consequences of a high CRF on, e.g., the TCA cycle seem to be present in females and males.

Age-, Menopausal Status-and Phenotypical/Clinical Variables-Adjusted Findings
A higher CRF is generally concomitant with a healthier body composition, adaptations in heart and arterial function, or cardiometabolic risk parameters. Hence, we next aimed to identify metabolite patterns that were independent of further assessed phenotypical and clinical variables associated with CRF. When adjusting for potential covariates (**), a comparatively smaller number of metabolites were correlated with the VO 2peak in both sexes. Thus, it can be suggested that many of the observed associations cannot be primarily and specifically attributed to CRF. As no relevant CRF-related information seemed to remain in the overall plasma metabolite profile, findings from PLS regression analyses are not further discussed.
After applying additional adjustments, seven acyl-alkyl-PCs, C5-carnitine, choline, glyceric acid, acetylornithine, and diacyl PC C28:1 still showed moderate positive correlations with the VO 2peak in females. Two of these acyl-alkyl-PCs (C40:3, C44:3) were also present in a cluster of 19 PCs that were related to CRF in the study of Wientzek et al. when controlling for sex, age, BMI, and waist circumference, among others [20]. The fact that we adjusted for more precise body composition measures and parameters of lipoprotein metabolism, which are known to be linked to both CRF [60] and PCs [35,38], could explain the study-specific results. Choline participates in multiple pathways of lipid or AA metabolism, serving as a precursor for PC species, the neurotransmitter acetylcholine, or betaine. Circulating choline can result from diet or PC breakdown [61] and was shown to be either lower [18] or higher [59,62] in fitter individuals. Even though a more efficient conversion of PCs to choline in trained individuals was assumed [62], it cannot be excluded that female KarMeN subjects with high CRF were also characterized by a higher dietary choline intake. In summary, sex-specific differences in CRF-related plasma metabolites were still detectable, e.g., for some PCs that remained positively correlated with female CRF but were not or even slightly negatively linked to male CRF. Moreover, the xenobiotic tartaric acid showed sex-specific correlations with CRF. This might be due to differences in the dietary intake of wine, vinegar, or grapes [63,64] in more or less fit women or men. The AA alanine showed a slightly negative correlation to the female VO 2peak but a slightly positive correlation to the male VO 2peak . In the literature, circulating alanine was mostly negatively linked to CRF status [15,23].

Sex-Specific CRF Explanation Models
Contrary to a previous attempt to explain the variability of CRF in the KarMeN population based on urinary metabolites [65], we were now able to identify sets of plasma metabolites that, together with clinical and phenotypical variables, contributed to a relatively good explanation of the VO 2peak in both sexes, after adjusting for age and menopausal status. In all approaches, the FM (%) entered the models as the first variable, already explaining 33.5% (females) or 42.3% (males) of the confounder-adjusted VO 2peak . Six (females) or five (males) further phenotypical or clinical parameters slightly improved the models (approach 1 selection). However, when phenotypical and clinical parameters, as well as metabolite variables, entered the models in a competing manner, plasma analytes led to fairly improved performance in both sexes, as demonstrated by an R 2 (adjusted) of 0.68 for females and 0.59 for males (approach 3 selection). Regarding females, acyl-alkyl-PC C40:3 was the second most important determinant of the VO 2peak . Interestingly, this PC sum parameter also showed relevant bi-and multivariate associations with CRF. As the relationship between acyl-alkyl-PC C40:3 and the VO 2peak did not appear to be influenced by other assessed plasma analytes, and the bivariate correlation persisted independently of adjustments, the ability of acyl-alkyl-PC C40:3 to predict CRF in females should be given special consideration in further studies. Regarding males, the TCA cycle intermediate malic acid was the second variable being included in the model after the FM (%) (approach 3 selection). In a study by Lustgarten et al., blood metabolites-based CRF explanation models were also sex-dependent. Seven serum metabolites in females and five serum metabolites in males explained 58 or 80% of CRF variability, respectively. However, as regression models were not adjusted for age, body composition, diet, or PA [22], it is uncertain if those metabolites are specifically indicative of CRF. As demonstrated by our stepwise regression models, the VO 2peak was largely determined by the FM (%) in the KarMeN study, support-ing the assumption that some associations between CRF and plasma metabolites might be mechanistically linked to clinical or phenotypical traits that were also influenced by chronic PA. In fact, Kujala et al. have shown that comparatively few blood metabolites remained significantly associated with CRF after adjusting for the body fat content in healthy young men [17]. Nevertheless, our findings could emphasize the additional value of metabolomics data for explaining the variability inherent in the VO 2peak , hinting at some possibly relevant plasma metabolites that could help to infer an individual's CRF status. Further studies will be needed to verify if those metabolites are indeed specific for CRF, and to what extent they are influenced by other functional or morphological characteristics of the human organism.

Strengths and Limitations
The major strength of this study is that it provides a systematic overview of the relationship between CRF and the plasma metabolome in a relatively large population consisting of both women and men, with a wide age range. The KarMeN study was characterized by highly standardized anthropometric and clinical examinations, as well as a strictly controlled procedure for blood collection, CRF assessment, and metabolomics analyses. To minimize the variability in metabolomics measurements, plasma samples were collected in the fasting state, and pre-menopausal women were examined within the luteal phase of their menstrual cycle. As the KarMeN study focused on healthy, non-smoking subjects with a normal to moderately high weight, excluding individuals with supplement use or hormonal treatment, the metabolic variation related to diseases, medication, or metabolic disorders was also markedly reduced. To control for known confounders from the very beginning, we conducted sex-specific analyses adjusted for age and menopausal status. Owing to the comprehensive characterization of study participants, further potential confounding factors related to body composition, clinical blood biochemistry, lung and arterial function, short-term and habitual PA, as well as diet could be considered. Another strength of this study is the applied multi-platform metabolomics approach, allowing the detection of a large number of plasma analytes from a broad range of biochemical classes and pathways. Limitations of the study include the cross-sectional design, as it does not allow the deriving of causal relationships. Furthermore, some of the plasma analytes showing relevant associations with CRF could not be identified with sufficient certainty and thus, regrettably, remained unknown.

Subjects and Study Design
The cross-sectional KarMeN study was conducted between March 2012 and July 2013 at the Division of Human Studies of the Max Rubner-Institut in Karlsruhe, Germany. Details on inclusion and exclusion criteria, as well as a comprehensive description of the study design and examination procedures, have already been published [26]. Briefly, 301 healthy, non-smoking individuals (172 men, 129 women) between 18 and 80 years of age were included. All subjects visited the study center three times and were thoroughly characterized by anthropometric, clinical, and functional examinations. Moreover, data on PA, diet and the menopausal status of female subjects were collected. Since the menstrual cycle is known to affect metabolite profiles [66], all premenopausal women were scheduled for examinations within their luteal phase. On the morning of the second study day, fasting plasma samples were collected using 9 mL EDTA plasma tubes (S-Monovette, Sarstedt, Nümbrecht, Germany). The plasma samples were immediately centrifuged at 1850× g at 4 • C, aliquoted into small portions, and cryopreserved at −196 • C until metabolomics analyses. Serum samples (S-Monovette Z-gel, Sarstedt, Nümbrecht, Germany) were collected for standard clinical biochemistry analyses.

Anthropometry and Body Composition Assessment
Bodyweight and height were measured in underwear and without shoes (Seca 285, Hamburg, Germany), and the BMI was calculated by dividing the body weight in kilograms by the height in meters squared. The body composition was assessed by dual-energy Xray absorptiometry (Lunar iDXA, GE Healthcare, München, Germany) and the LBM, FM, and VATM, as well as the bone mineral content (BMC), were determined with the supplementary software enCOREv16. The FM (%) was calculated by dividing the total FM by the total body weight. Approval for dual-energy X-ray absorptiometry measurements was received from the Federal Office for Radiation Protection (Z5-22462/2-2011-043).

PF and PA Assessment
As a measure of PF, CRF was assessed by a standardized incremental exercise test on a bicycle ergometer (Ergobike medical, Daum, Fürth, Germany). All participants started pedaling at 25 Watts and the workload was then augmented by 25 Watt every 2 min until individual exhaustion, as previously described [65,67]. The respiratory gas exchange was measured breath-by-breath by using a cardiopulmonary exercise testing system (MetaMax 3B, Cortex, Leipzig, Germany). Since the VO 2max could not be determined with certainty, as it requires the presence of a plateau in oxygen uptake, the VO 2peak, as the highest attained oxygen uptake during the test, was assessed and expressed relative to the bodyweight in mL kg −1 min −1 . During the entire procedure, the heart rate was recorded (T31 coded, Polar Electro GmbH Deutschland, Büttelborn, Germany). In addition, continuous hemodynamic monitoring was conducted by running a 12-channel electrocardiogram (CardioDirect 12, DelMar Reynolds GmbH, Feucht, Germany) and by measuring the BP every 2-3 min on the right upper arm (Boso-Carat Professional, Bosch + Sohn, Jungingen, Germany). As measures of PA, both short-term and habitual PA were determined [67]. Briefly, the level of short-term PA was assessed during a period of seven consecutive days by combined accelerometry and heart rate measurements (Actiheart, CamNtech, Cambridge, UK). The average AEE during the study week was finally calculated by the supplied software (Version 4.0.103) and given in kcal/day. To obtain the habitual PA for the last three months, participants filled out the standardized international physical activity questionnaire. The average weekly PA was calculated and expressed in MET-min/week.

Dietary Assessment
Food consumption for the day prior to blood sampling was assessed by conducting a 24-h recall in a personal interview, using the software EPIC-Soft, as described in detail elsewhere [26,68]. In order to evaluate the diet quality of participants, a modified version of the Healthy Eating Index was calculated, which was initially applied in the second German National Nutrition Survey ("Nationale Verzehrsstudie (NVS) II") [69] and adapted with minor modifications in the KarMeN study. The so-called HEI-NVS evaluates the overall diet quality, with scores ranging from 0 (low quality) to 110 (high quality).

Clinical Examinations
Clinical parameters like the HR rest , as well as systolic and diastolic BP, were measured after a resting period of at least five minutes in a sitting position (Boso-Carat Professional, Bosch & Sohn, Jungingen, Germany). The pulmonary function was assessed by spirometry (FlowScreen, CareFusion, Hoechberg, Germany) and the VC max , as well as the FEV1, were recorded. Moreover, arterial stiffness was determined (ARTERIOGraph, Medexpert, Budapest, Hungary) and the PWV was calculated. Standard clinical biochemistry analyses in fasting serum samples (e.g., hemoglobin (Hb), glucose, HbA1c, TGs, HDL-and LDL-cholesterol) were carried out by the certified medical laboratory MVZ Labor PD Dr. Volkmann (Karlsruhe, Germany) and insulin concentrations were determined with an enzyme-linked immunosorbent assay (ME E-0900, LDN, Nordhorn, Germany).

Metabolomics Analyses
To obtain a broad coverage of the plasma metabolome, a number of complementary (non-)targeted analytical techniques were applied. Quality control (QC) samples prepared by pooling plasma samples from all KarMeN participants were analyzed, along with plasma samples in all applied methods. The following section provides a brief summary of the different methods; further details are available in the supplement of Rist et al. [28].
Non-targeted two-dimensional gas chromatography (GC × GC)-MS analysis. Plasma samples were analyzed by a non-targeted GC × GC-MS-based approach using a Shimadzu GCMS QP2010 Ultra instrument equipped with a ZOEX ZX2 modulator [70]. The GC × GC-MS raw data files were then processed by non-targeted alignment using in-house developed R-modules [71]. By means of regularly injected QC samples, signal intensity drift, i.e., intra-and interbatch effects occurring during the measurement period, were corrected. With this method, a wide range of metabolites (including AAs, amines, organic acids, sugars, sugar alcohols, or polyols) could be detected.
Targeted gas chromatography (GC)-MS analysis of fatty acids. The chromatographic separation of plasma fatty acids usually requires the application of specialized polar columns and can thus not be conducted adequately by using a standard apolar × mediumpolar GC × GC column setup. Therefore, a previously described method [72] was applied to detect plasma fatty acids such as methyl esters, with minor modifications. By using a GC single quadrupole instrument (Shimadzu GCMS QP2010 Ultra) and a BPX90 column (Trajan Scientific), 48 fatty acids were finally determined in a quantitative manner.
Liquid chromatography (LC)-MS metabolite profiling using the Absolute IDQ™ p180 kit. Plasma samples were also utilized for targeted metabolite profiling using the Absolute IDQ™ p180 kit developed by Biocrates AG (Innsbruck, Austria). The general preparation and quantification procedure has already been described [73]. For the chromatographic separation of AA and biogenic amines, a Zorbax Eclipse XDB-C18 column (3 × 100 mm, 3.5 µm; Agilent, Waldbronn, Germany) equipped with a SecurityGuard™ column (C18, 4.0 × 3.0 mm; Phenomenex, Aschaffenburg, Germany) was used. PCs and SMs were analyzed by flow injection analysis into the analytical system, comprising a Nexera UHPLC system (Shimadzu) coupled to an API QTRAP15500 mass spectrometer (AB Sciex, Darmstadt, Germany). With this method, a variety of acylcarnitines, AAs, biogenic amines, SMs, and PCs were detected. While lysoPC species always have one acyl-bound fatty acid, other included PC species are characterized by two acyl-bounds (PC aa) or one acyl-and one alkyl-bound (PC ae) fatty acid, respectively. In general, each analyzed PC represents a sum parameter of different PC species with identical residue sums (e.g., PC ae C32:1 may consist of PC ae C16:0/C16:1, PC ae C18:0/C14:1, etc.).
Targeted LC-MS analysis of methylated amino compounds. The quantification of seven amino compounds in plasma was conducted by ultra-performance liquid chromatography-tandem MS, using an Acquity H-Class UPLC coupled to a Xevo TQD triple quadrupole MS (both from Waters, Eschborn, Germany), as previously established [74]. Plasma samples were diluted with acetonitrile after protein precipitation and separated by an inverse acetonitrile gradient on a polar hydrophilic interaction liquid chromatography column (Acquity BEH Amide, Waters, Eschborn, Germany). Target analytes, as well as deuterated internal standards, were monitored by positive electrospray ionization in multiple reaction monitoring mode.
Targeted LC-MS analysis of bile acids. 14 bile acids in plasma samples were quantified using a 1200 series HPLC system (Agilent, Waldbronn, Germany) coupled to a Q-Trap 3200 MS (AB Sciex, Darmstadt, Germany) [75].
Non-targeted NMR analysis. Plasma samples were analyzed by 1D-1 H-NMR spectroscopy. Briefly, they were measured at 310 K on an AVANCE II 600 MHz NMR spectrometer equipped with a 1H-BBI probe head and a BACS sample changer (Bruker BioSpin GmbH, Rheinstetten, Germany). All obtained plasma spectra were automatically phased with the Bruker AU program apk0.noe. Using the program AMIX 3.9.14 (Bruker BioSpin GmbH, Rheinstetten, Germany), they were then referenced to the EDTA signal and buck-eted graphically, so that buckets contained only one signal or group of signals and no peaks were split between buckets, whenever possible. Buckets were either annotated to a previously known and identified analyte, or registered as unknown. The identification of metabolites was carried out with Chenomx NMR Suite 8.1 (Chenomx, Edmonton, Canada). The detected analytes included organic acids, AAs, amines, and sugar alcohols.

Data Handling and Statistical Analysis
The data from the different analytical platforms were integrated into a common data matrix, consisting of 301 samples and 657 plasma analytes. With respect to the study participants, we excluded 49 individuals due to missing spiroergometry data (n = 40), technical errors during analyses (n = 7), implausibly low HR rest data (n = 1), and a missing plasma sample (n = 1). If identified metabolites were measured by more than one of the analytical platforms, those metabolites that were detected by the less quantitative method were excluded (n = 149). Further analytes were deleted if they had a detected frequency lower than 20% in either the female or male subgroup (n = 81). Thus, the final data matrix contained 252 individuals and 427 plasma analytes. The dataset, including metabolite data and metadata, is provided in Table S1. Prior to statistical analyses, metabolite data were transformed into Van der Waerden (VdW) scores. By using this rank-based inverse normal transformation, the data were converted into ranks, transformed to a scale between 0 and 1, and then, the corresponding standard normal quantiles were calculated. This transformation took the issue of values below the limit of detection into account and led to a uniform scale for all analytes, i.e., they were finally comparable between analytical platforms.
Based on sex-specific VO 2peak quartiles, the sex-specific VO 2peak data were divided into four quarters (q), and differences in basic characteristics between the subgroups of the corresponding participants of the quarters were examined by Welch ANOVA (chisquared test) for numeric (categorical) variables. All subsequent statistical analyses were conducted separately for sexes and the non-modifiable factors age (and menopausal status in females) were treated as confounders. Random forest regression algorithms considering age, phenotypical, and clinical variables were applied to impute missing values for AEE and PWV. Similar to the metabolite data, both the VO 2peak and considered phenotypical and clinical parameters were transformed into VdW scores prior to statistical analyses. To examine the sex-specific relationship between the VO 2peak and selected phenotypical and clinical parameters, Pearson correlations adjusted for age (and menopausal status in females) were calculated. Correlations were considered statistically significantly different from zero when the 95% confidence intervals (CIs) did not include zero.
Regarding metabolomics data analysis in sex-specific subgroups, three major aims were pursued: (i) to investigate the relationship between the VO 2peak and single plasma metabolites (bivariate association analyses), (ii) to determine the relationship between the VO 2peak and all plasma metabolites simultaneously (multivariate association analyses) and (iii) to identify a set of plasma metabolites possibly improving the explanation of the VO 2peak in the presence of phenotypical and clinical data (multiple regression analyses).

(i) Bivariate association analyses
To examine the sex-specific relationship between the VO 2peak and single plasma metabolites, adjusted for age and menopausal status, partial Pearson correlation coefficients (r) with 95% CIs were calculated. In a second step, correlations independent of further phenotypical and clinical variables were assessed. More specifically, we performed sexspecific correlation analyses by adjusting not only for the above-mentioned confounders but also for the following phenotypical and clinical parameters: LBM, FM (%), VATM, BMC, height, Hb, glucose, insulin, HbA1c, TGs, HDL and LDL cholesterol, HR rest , systolic and diastolic BP, PWV, VC max , FEV1, AEE, total MET, and HEI-NVS.

(ii) Multivariate association analyses
To analyze the relationship between the VO 2peak and all 427 plasma metabolite variables simultaneously, PLS regression analyses using nested cross-validation were conducted separately for women and men. The PLS analyses were either applied on confounder-adjusted metabolite variables or on metabolite variables additionally adjusted for the above-listed phenotypical and clinical parameters. The outer loop contained 20 random splits in a calibration dataset (containing 80% of all samples) and a test dataset (containing the remaining 20% of all samples). The data were preprocessed, including the formation of VdW-scores, respective adjustments, and unit variance scaling based on the calibration data. As the inner loop, a single random eight-fold cross-validation was used to tune the PLS model, based on the RMSE. Thereby, the number of predictive components was restricted to being at most ten. A rank for the obtained PLS regression models was assigned to each metabolite variable according to the negative absolute value of its regression coefficient. By calculating the geometric means of the ranks across the 20 random splits, a final rank product for each metabolite variable was obtained. The model performance was evaluated by the mean of RMSEs on the test samples across the 20 random splits. Moreover, 2500 permutations of the VO 2peak values were run, and the relative frequency of permutation-obtained rank products below the previously calculated rank products was assessed. If the relative frequency was ≤0.05, the contribution of a metabolite variable to the VO 2peak in a multivariate association was considered significant.

(iii) Multiple linear regression analyses
To assess the relationship between the VO 2peak and sets of phenotypical, clinical, and metabolite variables, three different exploratory multiple linear regression models were calculated for each sex, with the confounder-adjusted VO 2peak as the dependent variable, and with stepwise forward-selected confounder-adjusted clinical, phenotype, and metabolite variables as independent variables. In detail, the following approaches were applied: While in approach 2, all phenotypical/clinical variables entered the model as fixed variables before considering the plasma metabolites, in approach 3 all variables entered the model in a competing manner. To obtain a ranking of confounder-adjusted phenotypical, clinical, or metabolite variables according to their contribution for explaining the adjusted VO 2peak , the models were built by maximizing the coefficients of determination (R 2 ).
In addition to the ranking of variables, a single linear multiple regression model was calculated for each of the approaches 1 to 3 in order to obtain a manageable number of variables that, in combination, explained CRF. For variable selection, the previously described stepwise multiple linear regression analyses were performed on a calibration dataset (containing 80% of all samples) and the predictive accuracy of each step was assessed on the test dataset (containing the remaining 20% of the samples). The selection was stopped if the predictive accuracy decreased for the first time. In total, the analysis was repeated 1000 times with random assignments of samples into calibration and test datasets. Finally, the number of times each variable was present in those cross-validated stepwise regression models was counted (a relative frequency of 1 means that the variable was always considered in stepwise variable selection). All variables with a relative frequency ≥ 0.05 were then included in a final model with respect to approaches 1 to 3. The obtained final models were described by the adjusted R 2 and compared within each subpopulation. Statistical analysis was performed by using SAS JMP 11.0.0 (SAS Institute Inc. 2013, Cary, NC, USA) and the software R Version 4.0.0 [76], using the packages named caret [77], openxlsx [78], and missForest [79]. Figures were generated in Excel, 2016 or R, using the packages named ggplot2 [80], ggrepel [81], and ggpubr [82]. Supplementary PDF files were generated with LaTeX, using the package knitr [83]. R-scripts for the calculation of results and generation of figures and PDF files are provided in Files S6-S9.

Metabolite Classification
For the biological interpretation of CRF-related metabolite profiles, identified and putatively annotated metabolites, i.e., compounds with Metabolomics Standards Initiative (MSI)-level 1 or 2 [84], respectively, were manually assigned to 8 major and 32 specific pathways of human metabolism based on the information provided by the human metabolome database 4.0 [85] and the Kyoto Encyclopedia of Genes and Genomes PATH-WAY database [86]. Annotation of analytes detected using untargeted approaches, showing relevant associations with CRF, was performed as follows: spectra of GC × GC-MS analytes were matched against an in-house spectral library as well as against the FiehnLib and the NIST17 libraries. If matching was unsuccessful, structural hypotheses were derived depending on the presence of known diagnostic fragments (see File S9 in Ulaszewska et al. [87]) and additionally, especially in the case of sugars and sugar-like compounds, based on the compound's position in the two-dimensional chromatogram. Unfortunately, the spectra of CRF-related unknown analytes were mostly unspecific, which hampered structural elucidation. In the case of NMR, the identification of CRF-associated unknown analytes was not possible because the particular buckets contained either unspecific signals or overlapped peaks.

Conclusions
In summary, our findings demonstrated sex-dependent relationships between CRF and specific plasma metabolites in the KarMeN population. Apart from proving wellknown associations between CRF and, further, partly health-related phenotypical or clinical variables in both sexes, we could identify a number of PCs, lysoPCs, and SMs as being associated with the VO 2peak in either females or males, when controlling for age and menopausal status. However, independently of selected clinical or phenotypical variables, sex-specific CRF tended to be correlated with a rather small number of plasma metabolites primarily related to lipid-, AA-, or xenobiotics-related metabolism. Hence, many of the observed associations between CRF and metabolites were likely to be mediated by the considered clinical or phenotypical parameters. Although the variability of CRF was largely determined by the FM (%) in both sexes, our stepwise regression analyses revealed certain sets of plasma metabolites able to improve sex-specific VO 2peak explanation models. In particular, acyl-alkyl-PC C40:3 could be identified as a possibly interesting metabolite parameter for conclusions on CRF status in healthy females. Remarkably, CRF-associated metabolites have already been discussed as being reflective of exercise-induced adaptations in muscular energy metabolism (e.g., malic acid, succinic acid, acylcarnitines) or inversely linked to the development of cardiometabolic diseases (e.g., PCs, lysoPCs). Therefore, those metabolites might represent potential mediators of the performance-or health-enhancing effects of chronic PA. However, more research is needed to further clarify the mechanisms and metabolic pathways underlying sex-specific differences in CRF-associated metabolite profiles. Finally, we recommend future studies examining blood metabolic markers related to CRF to conduct sex-separated analyses and to consider age, menopausal status, body composition, and other health-related variables as covariates.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/metabo11070463/s1, File S1: Basic characteristics of the KarMeN participants according to VO 2peak quarters, File S2: Graphical overview of associations between the VO 2peak and plasma metabolites, File S3: Classification of metabolites with significant bivariate correlations to metabolic pathways, File S4: Evaluation of PLS approaches, File S5: Volcano plots of VO 2peak -associated plasma metabolite patterns, File S6: Calculation, File S7: Results, File S8: Output, File S9: Supplement, Table S1: Dataset including metabolite data and metadata, Table S2: Results of bi-and multivariate association analyses, Table S3: Results of multiple linear regression analyses.  peak oxygen uptake