Muscle contractile activity is supported by a cascade of biochemical and metabolic events fueled by the conversion of energy substrates to ATP. To study the effects of controlled exercise on muscle energy metabolism, we developed an in vivo stimulation protocol that mimics moderate exercise in fed anesthetized mice (Figure 1). The protocol consists of controlled bursts of muscle activity alternating with periods of rest, repeating over 22 min. When applied to control and Il15ra^–/– mice, the stimulation evoked the simultaneous contraction of all the muscles of the anterior and posterior compartments of the targeted limb. For each animal, one limb underwent stimulation while the contralateral was used as unstimulated control. The protocol did not induce complete muscle exhaustion, as indicated by the only partial decrement of EMG amplitude and root mean square (RMS) (Figure 1B). Consistent with moderate exercise, peripheral capillary oxygen saturation (SpO₂) remained stable, while heart rate (HR) increased moderately as a consequence of the repeated stimulation (Figure 1B).

Previous studies performed on a B6;129×1-based Il15ra^–/– model had raised the concern of potential differences in background compared to the control strain used (B6129SF2/J). To address this concern, we generated a new C57BL/6-based Il15ra^–/– line by crossing IL15RA-flox mice to mice expressing Cre driven by the ubiquitous promoter EIIA (Figure 2A) and performed every experiment using control and Il15ra^–/– littermates. Genomic PCR detected the successful recombination and full excision of exons 2/3 (Figure 2B), the region coding for the IL15 binding site in IL15RA. We confirmed the successful excision in all tissues, including muscle, brain, liver, spleen and adipose tissue (Figure 2C). To test whether the C57BL/6-Il15ra^–/– line (from now on called Il15ra^–/–) retained the previously described muscle phenotype, we compared EDL muscles from Il15ra^–/– and littermate controls in terms of isometric strength (Supplementary Table S1), resistance to fatigue (Figure 2D) and recovery capacity (Figure 2E). Measurements of muscle strength and weight were similar between genotypes (Supplementary Table S1). Consistently with the previous reports, Il15ra^–/– EDL retained higher fatigue resistance (Figure 2D) and recovery capacity (Figure 2E), suggesting differences in energy metabolism. In vivo, Il15ra^–/– retained the overall normal circadian rhythm of locomotor activities, although they were more active than controls when tested for cumulative wheel running activity (Figure 2F) and beam breaks (Figure 2G, – top panels), with higher oxygen consumption and energy expenditure especially during the dark phase (Figure 2G, – bottom panels).

Because changes in transcriptomics or proteomics are unlikely to drive the acute response to a single bout of muscle stimulation by themselves, our approach focused on the analysis of muscle metabolome. We performed targeted metabolomics on unstimulated and stimulated EDL muscles (Figure 3A). A total of 321 metabolites were detected and identified by GC-TOF MS analysis (Supplementary Table S2). To model the effects of stimulation and genotype, we applied a supervised partial least squares discriminant analysis (PLS-DA). The significance of class discrimination, assessed using a permutation test based on prediction accuracy, was calculated to be P = 0.007. The PLS-DA model distinguished between the contribution of two main components, accounting respectively, for 54 and 12% of the total variance. Component 1 best separated the effect of stimulation, while component 2 separated the effect of genotype (Figure 3B). A loading plot representation (Figure 3C and Supplementary Table S3) highlighted the contribution of carnitine species and TCA cycle intermediates to the variance associated with component 1 (stimulation), and the contribution of gut microbiota metabolites to the variance of component 2 (genotype). A ranked representation of the 25 highest variable importance projection (VIP) scores for each component confirmed the contribution of carnitines, fatty acids, and TCA intermediates to the variation associated with component 1 (Figure 3D and Supplementary Table S4). Carnitines and fatty acids were also among the 25 highest VIPs for component 2 (Figure 3E and Supplementary Table S4), suggesting an interaction between genotype and stimulation.

By applying a two-way mixed ANOVA analysis (Genotype × Stimulation) with statistical significance at α ≤ 0.05 after FDR adjustment, 195 metabolites were significantly affected by stimulation (123 more than 1.5×), while only 7 were significantly different between genotypes (5 more than 1.5×) (Supplementary Table S2). Because of the magnitude of the changes associated with stimulation, we initially dropped the genotype factor and focused on the comparison between stimulated and unstimulated muscles. The majority of stimulation-induced changes were toward higher metabolite levels (107 up and 16 down – Figure 4A). The most represented metabolite superclasses were fatty acids (34.1%), amino acids (23.6%), carbohydrates (14.6%), and nucleotides (8.9%) (Figure 4B). Consistent with the PLS-DA model, carnitines were predominant within the group of upregulated metabolites (Figure 4C). Pathway analysis suggested that the most represented pathways impacted by stimulation were amino acids metabolism (Asp – P = 4.20E-12; Phe, Tyr – P = 1.79E-09; Glu – P = 1.03E-06; Arg, Pro – P = 8.57E-07), purine metabolism (P = 1.13E-09), citric acid cycle (P = 5.00E-07), urea cycle (P = 2.30E-06), glycolysis (P = 1.36E-04) and oxidation of branched-chain fatty acids (P = 1.10E-06) (Figure 4D). Due to the lack of Kegg identifiers for the fatty acids in our dataset, both measures of enrichment and impact on the entries associated with fatty acids metabolism were likely underpowered. A ranked visualization of the most significantly regulated metabolites highlighted the important role of TCA and purine metabolism in the response to stimulation in fast muscles (Figure 4E). The most differentially regulated metabolites in the dataset were inosine monophosphate (IMP, 54×, P = 2.34E-11) and deoxyadenosine monophosphate (dAMP, −5.7×, P = 3.10E-05). IMP and dAMP are intermediates in the reactions catalyzed by adenylate kinase (which converts 2 ADP – >ATP and AMP) and myoadenylate deaminase (AMP or dAMP – >NH₃ and IMP), which are activated during medium to high-intensity exercise in response to ATP utilization (−1.8×, P = 5.85E-05) and changes in pH (Sahlin et al., 1989; Tullson and Terjung, 1991, 1999). Metabolism of amino acids such as glutamine, glycine, and aspartate plays an essential anaplerotic role in muscle recovery by providing intermediates for the TCA, e.g., by converting aspartate (−4×, P = 2.69E-07) to fumarate (7×, P = 1.48E-09) (Gibala et al., 1997). The moderate exercise-mimicking stimulation also decreased the creatine-phosphate/creatine ratio by ∼50% and caused a 1.6× accumulation of intramuscular lactate in both genotypes. Together, these results show that the activation of oxidative metabolism following the application of our in vivo stimulation protocol is consistent with a bout of moderate-intensity exercise.

A pathway diagram of the main muscle energy metabolites affected by stimulation (Figure 5) suggested that, upon moderate exercise, glucose derived from increased uptake and glycogenolysis was shunted from glycolysis to the pentose phosphate pathway, generating reducing equivalents and ribose-5-phosphate that entered the purine nucleotide pathway as ZMP (also known as AICAR; 3.78×, P = 2.77E-04), a potent activator of AMPK and a known regulator of energy metabolism (Ezagouri et al., 2019). Pyruvate levels were not affected by stimulation, suggesting that the increased acetyl-CoA and TCA cycle intermediates were mainly derived from fatty acid degradation and anaplerotic reactions. Interestingly, glycolytic intermediates accumulated upstream of GAPDH. GAPDH enzymatic activity may be allosterically inhibited by multiple mechanisms, including nitrosylation in presence of high levels of NADH (Mohr et al., 1996) (25.95×, P = 1.05E-03), or by increased concentrations of fumarate (Blatnik et al., 2008) (7.48×, P = 1.48E-09). The application of our in vivo stimulation protocol allowed us to map the metabolic consequences of a single bout of moderate intensity exercise, highlighting the important contribution of fatty acid metabolism and anaplerotic reactions to muscle energetics.

Il15ra^–/– mice perform significantly better than controls when tested for muscle fatigability and recovery capacity. We, therefore, focused on the differences associated with genotype in unstimulated and stimulated muscles. With this approach, the number of significantly different metabolites between control and Il15ra^–/– muscles were 8 before and 11 after stimulation (Figure 6). Consistent with the PLS-DA model, the metabolites showing the most significant variations between genotypes were the microbiota metabolites 3-phenylpropionate, hippuric acid, and indoxyl sulfate, suggesting an impact of IL15RA deficiency on gut immune system and microbiota. Some fatty acids, as well as intermediates of the glucuronic acid pathway (part of the pentose phosphate pathway), were significantly higher in Il15ra^–/– than controls. In line with previous reports showing that lack of IL15RA has a role in bone homeostasis (Djaafar et al., 2010; Loro et al., 2017), glucosamine, a metabolite important for joint health, was induced 2.2-fold by stimulation and was 1.8 times higher in the stimulated Il15ra^–/– than stimulated controls (Figure 6A). The AMPK activator ZMP trended to be higher in Il15ra^–/– than in control muscles (1.88× in unstimulated and 2.36× in stimulated muscles). Of 48 detected fatty acids, 24 were significantly higher after stimulation in control muscles (Figure 6B). In Il15ra^–/– unstimulated muscles, the basal pool size of long-chain fatty acids (LC-FA) was higher than in controls, but did not change after stimulation (Figures 6B,C). Given the previously described enhancement of FA uptake into Il15ra^–/– isolated muscle mitochondria (Loro et al., 2015), we considered whether the pool of acyl-carnitines was significantly elevated in Il15ra^–/– compared to control muscles. The levels of 19 detected acyl-carnitines (C2–C20) increased significantly upon stimulation and were higher in Il15ra^–/– compared to control muscles (Figure 6D).

To obtain a visual indication of whether IL15RA ablation can function as exercise mimetic, we then calculated the Euclidean distance between the levels of each metabolite in Il15ra^–/– muscles and in the unstimulated and stimulated control conditions. We found that approximately 30% of the metabolites detected in Il15ra^–/– unstimulated muscle, including fatty acids, amino acids, and nucleotides, were closer to the levels found in stimulated control muscles (Figure 7A, distance closer to 1).

Taken together, these results demonstrate that lack of IL15RA primes basal muscle energy metabolism to a state consistent with a moderate exercise condition. In such a state, the previously described activation of lipases and the AMPK pathway (Loro et al., 2015) cooperates with the increased availability of fatty acids to promote reliance on fatty acids and support the higher energy expenditure and exercise capacity (Figure 7B).

The predominantly fast EDL muscle, composed of approximately 80% type 2B-2X glycolytic fibers, responded to stimulation with a strong elevation of intermediates of oxidative metabolism. Upon stimulation, 40% of the detected metabolites were significantly changed. 13% of them were down-regulated compared to unstimulated condition. Among the downregulated metabolites were energy substrates such as ATP, GTP and amino acids. Amino acids such as aspartate and glutamine are substrates for the anaplerotic generation of TCA intermediates (e.g., succinate). Consistent with a moderate-intensity exercise protocol, the majority (87%) of the up-regulated metabolites belonged to FAO, TCA cycle, purine nucleotide metabolism, and pentose phosphate pathways. Similarly to our study, exercise in the form of spontaneous wheel running or treadmill increase lipolysis during the resting phase (Monleon et al., 2014; Ezagouri et al., 2019; Sato et al., 2019). Despite higher stimulated glucose levels (2.83×, P = 8.2E-04) from uptake and/or glycogenolysis), glycolytic intermediates accumulated upstream of GAPDH, suggesting that its enzymatic activity could be reduced by the accumulation of NADH (25.95×, P = 1.05E-03) (Mohr et al., 1996) or fumarate (7.48×, P = 1.48E-09) (Blatnik et al., 2008).

The most up-regulated metabolite upon stimulation was the nucleoside IMP (inosine monophosphate, 54.1×, P = 2.34E-11). The deamination of AMP to IMP occurs during moderately intense muscle activity and is catalyzed by the enzyme adenosine deaminase (AMPD), the first and rate-limiting step of the purine nucleotide cycle. The conversion of AMP to IMP is advantageous for the muscle, as it maintains the ATP:AMP ratio and spares adenine nucleotides from degradation to purine bases, allowing their quick reamination during muscle recovery (Sabina et al., 1980). Supporting the role of IMP and purine metabolism during muscle activity, AMPD activity was recently found to be elevated in EDL muscles of mice lacking HDAC3, accounting for their increased fatigue resistance (Hong et al., 2017). Another source of IMP derives from the conversion of glucose-6-phosphate to ribose-5-phosphate in the pentose phosphate pathway, which is then converted to ZMP (3.78×, P = 2.77E-4 induction by stimulation) and enters the purine metabolism cycle. During this process, the concomitant generation of NADPH may support the exercise-induced activation of antioxidant enzymes (Frasier et al., 2013). ZMP is also a potent activator of AMPK, which in turn regulates energy metabolism and exercise capacity by promoting glucose uptake and FAO (Narkar et al., 2008). ZMP-mediated AMPK activation could therefore participate in the regulation of the network of kinase reactions recently identified by phosphoproteomic investigations of exercised muscles (Hoffman et al., 2015; Potts et al., 2017).

In summary, in vivo muscle stimulation allowed us to map the metabolic reaction (oxidative and anaplerotic) activated to ensure the replenishment of energy substrates and the recovery from fatigue. We then compared such map with the metabolic reactions activated in Il15ra^–/– muscles.

The enhanced fatigue resistance and energy expenditure (Figure 2) of Il15ra^–/– mice is accompanied by a significant reprogramming of muscle metabolism, with increased basal IMTG and higher AMPK activation. Such changes are independent of the exercise state and lead to higher basal mitochondrial uptake and oxidation of fatty acids in isolated mitochondria (Loro et al., 2015).

Here, we hypothesized that if IL15RA blockade has exercise mimetic effects, it would cause a shift in basal energy metabolism toward that of stimulated control muscles. Indeed, the levels of several metabolites in the unstimulated Il15ra^–/– muscles were intermediate between unstimulated and stimulated controls (Figure 7A). A significant effect of genotype was found in microbiota metabolites, amino acid metabolism intermediates, and in the pool sizes of fatty acids and carnitines (Figure 6). Glycolytic intermediates were similar between Il15ra^–/– and controls, suggesting normal muscle glucose utilization during exercise.

The higher energy output required for sustaining Il15ra^–/– exercise capacity is therefore due to enhanced IMTG utilization, fatty acids mobilization, higher mitochondria content (Pistilli et al., 2011, 2013) and cristae complexity (Loro et al., 2018), similarly to what has been described in the muscles of highly trained athletes (Nielsen et al., 2017). In contrast to other experimental models (e.g., overexpression of PGC1-alpha), ablation of IL15RA exerts such effects without changing muscle composition (Loro et al., 2018). It is noteworthy that such small but significant metabolic enhancements translate into clear benefits in terms of muscle contractility, fatigability, and energy expenditure during the mostly inactive light phase (ZT4-8). The magnitude of these differences will be likely bigger during the active dark phase.

The mechanisms responsible for the increased FAO in Il15ra^–/– could be mediated by AMPK and CPT1B activity. AMPK activation in Il15ra^–/– muscle was not accompanied by substantial changes in the AMP/ATP ratio, suggesting that other factors such as ZMP (elevated by 1.88× and 2.36× in unstimulated and stimulated Il15ra^–/–) or calcium-dependent signaling could be causing AMPK phosphorylation (Funai et al., 2013). Together with changes in metabolites, protein phosphorylation is an important mediator of the circadian control of metabolism and physiology (Robles et al., 2017) and its kinetics are consistent with the timescale of events happening during a short bout of exercise (Hoffman et al., 2015; Potts et al., 2017). Hence, this study can serve as a starting point for future investigations correlating multiple parameters such as phosphoproteomic and metabolomic responses to different stimulation intensities, different phases of the day, or multiple rodent disease models.

In summary, our data show that IL15RA ablation elevates basal metabolic rates to levels consistent with moderate exercise, in turn favoring FAO reactions as the primary fuel. Together with our previously published work on exercise capacity and resistance to diet-induced obesity (Pistilli et al., 2011; Loro et al., 2015), these findings strengthen the rationale for the future development of new pharmacological exercise mimetic strategies by blocking IL15RA.

All animal experiments were conducted in accordance with protocols approved by the Institutional Animal Care and Use Committee of the University of Pennsylvania. Mice were housed at 22°C under a 12:12-h light-dark cycle (7:00 AM – 7:00 PM) with food and water provided ad libitum. All founder mouse lines were purchased from Jackson Laboratories. C57BL/6-Il15ratm2.1Ama/J (IL15RA-flox, stock number 022365) and B6.FVB-Tg(EIIa-cre)C5379Lmgd/J (EIIA-cre, stock number 003724) were crossed to obtain a progeny with ubiquitous deletion of the floxed alleles in all tissues. The correct genotype was confirmed by PCR of floxed alleles (according to Jackson protocol for strain 022365), Cre-recombinase (according to Jackson protocol for strain 003966), Il15ra exon 2 (included in the region deleted by EIIa-Cre, Forward 5′-ATTGAGCATGCTGACATCCG-3′, Reverse 5′- TCCAGTGGGCAACATTTGTG-3′, product size 145 bp) and Il15ra exon 4 (outside the region deleted by EIIa-Cre, Forward 5′-TGCAGCAACAATGACCTTGG-3′, Reverse 5′- TATGTGGGAAACTGGCCTGTC-3′, product size 96 bp). These mice were then backcrossed to C57BL/6J (stock number 000664) mice to obtain heterozygous Il15ra^± on a C57BL/6J background, used for the experimental breeding. Ten-week-old male control (Control) and Il15ra^–/– littermates were used for the experiments. Spontaneous wheel running (Columbus Instruments) was monitored for approximately 7 days during which mice were single-housed in cages with ad libitum access to food and water. Energy expenditure measurements were performed on a CLAMS setup (Columbus Instruments) at the Mouse Phenotyping, Physiology and Metabolism Core of the University of Pennsylvania. In line with previous studies, all in vivo experiments were performed between ZT4 and ZT8.

For in vivo stimulation, mice (4 controls and 4 Il15ra^–/–) were anesthetized with 1% isoflurane. Depth of anesthesia was confirmed by toe pinch and by monitoring HR and SpO₂ with a pulse oximeter (PhysioSuite, Kent Scientific). The sciatic nerve was exposed at both sides, crushed distally to the muscle terminations to prevent retrograde propagation of the stimuli, and connected to a custom-made bipolar stimulating electrode driven by a Grass S48 stimulator. Two recording platinum-iridium (bare 0.13 mm, coated 0.02 mm) needle-electrodes were applied at fixed distance (5 mm), one in the gastrocnemius and the other near the Achilles tendon. The ground electrode was connected to the skin of the back of the mouse. The recording electrodes were then connected to a Warner Instruments DP-311 differential amplifier and finally to an A/D converter (ADInstruments). Traces were acquired with LabChart 8 (ADInstruments). One leg of each mouse was electrically stimulated, the other was used as sham-operated unstimulated control and dissected prior to the start of the stimulation, paying attention to minimize bleeding. The EMG signal was used to monitor the efficiency and consistency of the high-intensity stimulation protocol. The sequence of stimuli was based on the protocol originally described by Baar and Esser (1999) and was previously used to investigate phosphoproteomic changes in exercised muscle of mice (O’Neil et al., 2009; Potts et al., 2017). Throughout the experiment, vital parameters were monitored with a pulse oximeter (Figure 1A). The stimulation protocol used for the study consisted of 3-s trains of 0.1 ms pulses repeating at a frequency of 100 Hz (Figure 1B, – block A). Each 3-s train was repeated 6 times, alternating with 10 s of rest (Figure 1B – block B). This stimulation pattern was repeated 10 times with 1 min of rest between each repetition, for a total duration of 22 min (Figure 1B, – block C).

Stimulated muscles were collected immediately at the end of the protocol. EDL muscles were quickly flash frozen after dissection and processed for metabolite extraction. An independent set of two control and two Il15ra^–/– mice (not included in the current study) was used to validate the efficacy of the stimulation protocol and the capacity of our metabolomics approach to detect changes associated with muscle activity.

Muscle physiological analysis was performed on isolated EDL muscles using an Aurora Mouse 1200A System equipped with Dynamic Muscle Control v.5.415 software. EDL muscles were dissected and analyzed in constantly oxygenated Ringer’s solution (100 mM NaCl, 4.7 mM KCl, 3.4 mM CaCl₂, 1.2 mM KH₂PO₄, 1.2 mM MgSO₄, 25 mM HEPES, 5.5 mM D-glucose) at 24°C. The twitch stimulation protocol applied was a single stimulus with a duration of 0.2 ms. Muscle length was adjusted to obtain the maximal twitch response and this length was measured and recorded as optimal length (L₀). For measuring tetanic maximal force generation, the stimulus was repeated at a frequency of 120 Hz for 500 ms. 5 min were allowed between two tetanic contractions to ensure muscle recovery. For induction of fatigue, 5 min after the last maximal tetanic contraction, muscles were stimulated every second for 8 min using 40-Hz pulses lasting 330 ms. Following the fatigue protocol, a burst of 50 maximal tetanic contractions (120 Hz for 500 ms) was applied to ensure complete fatigue across all samples. The recovery protocol started 1 s after the last burst contraction. A maximal tetanic stimulation (120 Hz for 500 ms) was given every 5 min for 30 min, and the force recovery was expressed as the percentage of the maximal isometric tetanic force. Muscle cross-sectional area (CSA) of EDL muscles was calculated dividing the muscle mass by the product of the muscle density coefficient (1.06 g/cm³), muscle L₀, and the fiber length coefficient (0.45 for EDL). Specific force was determined by normalizing maximum isometric tetanic force to CSA.

Metabolomics was performed at Penn DRC Metabolomics Core. To extract metabolites from tissue samples, frozen tissue samples were ground at liquid nitrogen temperature with a Cryomill (Retsch, Newtown, PA, United States). The resulting tissue powder (∼20 mg) was weighed and then extracted by adding −20°C extraction solvent (40:40:20 methanol:acetonitrile:water), vortexed, and immediately centrifuged at 16,000 × g for 10 min at 4°C. The volume of the extraction solution (μL) was 40× the weight of tissue (mg) to make an extract of 25 mg tissue per mL solvent. The supernatant was collected for LC-MS analysis.

A quadrupole-orbitrap mass spectrometer (Q Exactive, Thermo Fisher Scientific, San Jose, CA, United States) operating in negative or positive ion mode was coupled to hydrophilic interaction chromatography via electrospray ionization and used to scan from m/z 70 to 1000 at 1 Hz and 75,000 resolution. LC separation was on a XBridge BEH Amide column (2.1 mm × 150 mm, 2.5 μm particle size, 130 Å pore size; Waters, Milford, MA, United States) using a gradient of solvent A (20 mM ammonium acetate, 20 mM ammonium hydroxide in 95:5 water: acetonitrile, pH 9.45) and solvent B (acetonitrile). Flow rate was 150 μL/min. The LC gradient was: 0 min, 85% B; 2 min, 85% B; 3 min, 80% B; 5 min, 80% B; 6 min, 75% B; 7 min, 75% B; 8 min, 70% B; 9 min, 70% B; 10 min, 50% B; 12 min, 50% B; 13 min, 25% B; 16 min, 25% B; 18 min, 0% B; 23 min, 0% B; 24 min, 85% B; 30 min, 85% B. Autosampler temperature was 5°C, and injection volume was 3 μL. Data were analyzed using the MAVEN software.

Missing or equal to zero values (11/5136) were replaced with a small value, calculated as the half of the minimum detection limit (minimum positive values in the original data). The assumption of this approach is that most missing values are caused by low abundance metabolites (i.e., below the detection limit). The values were for the following metabolites: 6-phosphogluconate (1), C24:3 (2), NADH (2), acetyl-CoA (1), deoxyadenosine (4 = 2 control unstimulated and 2 control stimulated), Tetradecadienoyl carnitine C14:2 (1).

For all the analyses, metabolomics data were log₂-transformed. Heatmap of log₂-transformed, row-wise normalized (mean of 0 and variance of 1) values was generated with the ComplexHeatmap R package (Gu et al., 2016) considering euclidean distance and average linkage. PLS-DA and pathway analysis were performed using the Metaboanalyst online platform (Chong et al., 2018) with standard settings. PLS-DA is a supervised method that uses multivariate regression techniques to maximize the covariance between exploratory and categorical codependent variables. The significance of class discrimination was assessed using a permutation test (2000 permutations, P = 0.007). The variance on each component in relation to metabolites can be graphically represented as loading plots. The PLS-DA analysis also returned a VIP score. The VIP score describes the weighted contribution of each metabolite to the variation of each component. The 25 most important features in each component, as identified by PLS-DA, were plotted. Pathway representation in Figure 5, was obtained with Cytoscape 3.7.1. Euclidean distance for each detected metabolite in Figure 7A, was calculated as the absolute value of the difference between the log₂-transformed mean values at the different conditions. The obtained distances were then rescaled to a range from 0 (unstimulated control) and 1 (stimulated control).

Data are presented as mean ± standard error of the mean. All statistical analyses were performed in R 3.5.1 (R Foundation for Statistical Computing, Vienna, Austria¹) using the packages lme4 (Bates et al., 2015) and lmerTest (Kuznetsova et al., 2017). Data visualization was performed with the ggplot2 R package (Gómez-Rubio, 2017) or GraphPad Prism 8 (GraphPad Software, La Jolla, CA, United States). All statistical tests were considered significant at α ≤ 0.05, unless stated otherwise. Univariate analysis of the metabolomics data was performed by 2-way ANOVA (Genotype × Stimulation) followed by the Benjamini-Hochberg FDR adjustment. To compare the effects of stimulation in the two genotypes, pairwise multiple comparisons were performed using the R package emmeans and adjusting the raw P-values using the Sidak method. ^∗P < 0.05, ^∗∗P < 0.01, and ^∗∗∗P < 0.001.