A preliminary model of football-related neural stress that integrates metabolomics with transcriptomics and virtual reality.

Research suggests contact sports affect neurological health. This study used permutation-based mediation statistics to integrate measures of metabolomics, neuroinflammatory miRNAs, and virtual reality (VR)-based motor control to investigate multi-scale relationships across a season of collegiate American football. Fourteen significant mediations (six pre-season, eight across-season) were observed where metabolites always mediated the statistical relationship between miRNAs and VR-based motor control ( p S o b e l p e r m ≤ 0.05; total effect > 50%), suggesting a hypothesis that metabolites sit in the statistical pathway between transcriptome and behavior. Three results further supported a model of chronic neuroinflammation, consistent with mitochondrial dysfunction: (1) Mediating metabolites were consistently medium-to-long chain fatty acids, (2) tricarboxylic acid cycle metabolites decreased across-season, and (3) accumulated head acceleration events statistically moderated pre-season metabolite levels to directionally model post-season metabolite levels. These preliminary findings implicate potential mitochondrial dysfunction and highlight probable peripheral blood biomarkers underlying repetitive head impacts in otherwise healthy collegiate football athletes.

Introduction

Repetitive head acceleration events (HAEs) have been shown to induce neurophysiological and biochemical abnormalities without triggering clinically observable symptoms (Kawata et al., 2017; Mainwaring et al., 2018; Nauman et al., 2020). The biological effects of these “subconcussive injuries” can mirror changes observed following concussion, a condition diagnosed in approximately 1.6–3.8 million people in the US annually (Langlois et al., 2006). Extensive research has begun uncovering the molecular mechanisms of concussion/HAEs; however, their relationship to other functional change, such as neurobehavioral impairment, appears complex (Nauman et al., 2020). This complexity, in the face of minimal clinical symptoms, means many tests are inconsistent or insensitive to subtle behavioral change (Dessy et al., 2017; Martini and Broglio, 2018). Among the more consistent findings are those related to motor control, including disturbances in balance accuracy (accuracy correcting balance), sensory motor reactivity (speed correcting balance), and spatial navigation, together referred to as “motor control” henceforth. These three behavioral features were identified by Alexander Luria (Luria, 1973) as being commonly affected by a broad range of head injury severity. These motor control features have been quantified more recently using virtual reality (VR) as a means to assess behavioral dysfunction potentially missed in standard clinical examinations (Slobounov et al., 2006).

The emergence of concussion-related motor control dysfunction can be interpreted, in part, by the neuroinflammatory model of brain injury, which is supported by transcriptomic (miRNA) (Di Pietro et al., 2021; Papa et al., 2019) and neuroimaging studies in humans (Jang et al., 2019; I K Koerte et al., 2015). The referenced miRNA studies were noteworthy for identifying abnormalities across a panel of biomarkers in collegiate American football athletes (Bhomia et al., 2016; Papa et al., 2019). These miRNAs were not only elevated at the end of the football season, but also prior to the season (Pre), and a subset of these miRNAs were associated with regional cerebral blood (rCBF) flow and postural control at Pre suggesting a chronic neuroinflammatory syndrome (Chen et al., 2020; Papa et al., 2019). These miRNAs increase following TBI (Bhomia et al., 2016; Papa et al., 2019; Sharma et al., 2014) and have established roles in inflammation (Jee et al., 2012b; Paraskevi et al., 2012; Reddycherla et al., 2015; Wu et al., 2011; Zheng et al., 2015), neurodegeneration (Basak et al., 2016; Denk et al., 2015; Khoo et al., 2012; Liu et al., 2017), apoptosis (Li et al., 2016; Lv et al., 2018), and oncogenesis (Che et al., 2019; Kobayashi et al., 2012; Kohram et al., 2018; Li et al., 2011; Ling et al., 2019; Tommasi et al., 2016). Multiple concussion studies support the neuroinflammatory model (Alosco et al., 2020; Aungst et al., 2014; Bailes et al., 2013; Inga K Koerte et al., 2015; McAteer et al., 2016; Shultz et al., 2011; Tagge et al., 2018; Xu et al., 2016), and others suggest that mitochondrial dysfunction may underlie the inflammatory response (Harland et al., 2020; Hernández-Aguilera et al., 2013; Kumar Sahel et al., 2019; Yoo et al., 2020).

Mitochondrial dysfunction has been consistently observed across animal studies, and raised as a hypothesis for human studies, of concussion and HAEs (González-Domínguez, 2016; Hiebert et al., 2015; Vagnozzi et al., 1999; Xiong et al., 1997). Dubbed the powerhouse of the cell, mitochondria are critically important for energy production and cellular respiration (Siekevitz, 1957). Damage to these organelles can result in serious cellular, and potentially systems-level, dysfunction. One major mitochondrial role involves the oxidation of fatty acids (FAs) into functional metabolites that fuel the TCA cycle and downstream energy production (Daley et al., 2016; Schulz, 1991). Preclinical research suggests that metabolites connected to mitochondrial energy metabolism and oxidative stress (Bhat et al., 2015; Bhatti et al., 2017; Ježek et al., 2018) have been linked to neuroinflammation and may be altered in organisms with head injury (Signoretti et al., 2010; Verweij et al., 1997, 2000).

In this study, we sought to determine whether metabolomic measures were altered across a collegiate football season and hypothesized that altered metabolites would be related to changes in neuroinflammatory miRNAs and VR-based motor control in a consistent manner. This hypothesis specifically focused on the relationship of changes in these three types of measures, as well as their baseline relationships prior to the season. We further hypothesized that molecular measures would be statistically related to HAEs and tested, in a post-hoc fashion, if HAEs across a season would mediate/moderate the relationship of Pre-season to Post-season measures. The same subjects were studied at Pre and postseason (Post) given they were under the same nutritional, workout, and lifestyle management program, thus serving as an adequate comparison group (see Methods). For integration of multiple omic measures (i.e., multi-omics (Chu et al., 2019)) with computational behavior (i.e., mathematical psychology and the advanced technology like VR used to study behavior), we focused on a set of 40 metabolites with roles in energetics, inflammation, or exogenous consumption, a candidate panel of nine miRNAs (miR-20a, miR-505, miR-92a, miR-195, miR-93p, miR-30d, miR-486, miR-3623p, and miR-151-5p) connected to neuroinflammation and concussion (Bhomia et al., 2016; Papa et al., 2019; Sharma et al., 2014), and four VR-based metrics of motor control. Twenty metabolites were chosen based on the highest mean decrease accuracy (MDA) (q-value; after correction) from a random forest analysis (see Methods for details) and the other 20 were chosen based on literature survey for hypothesized involvement in neuroinflammatory-related injury (Abdul-Muneer et al., 2015; Lozano et al., 2015; Rodriguez-Rodriguez et al., 2014). Our hypotheses constrained VR performance to be the dependent variable (Y), and we used mediation to move past simple association (Lockhart et al., 2011; Mackinnon, 2008). We hypothesized that metabolomic measures (mediator, M) would consistently carry the relationship between neuroinflammatory miRNAs (independent variable, X) and VR-based motor control (Y). A secondary hypothesis would be used to confirm that miRNAs (M) would not carry the relationship between metabolomic measures (X) and VR-based motor control (Y). If both hypotheses were true, and if a consistent set of relationships were observed across analyses, this would support a statistical model framing the relationship of transcriptome, metabolome, and behavior in athletes sustaining repetitive HAEs. We based our hypothesis on (1) previous work which found miRNAs to be elevated prior to and following the football season (Papa et al., 2019), (2) previous work which showed the same panel of miRNAs consistently acting as the IV, with rCBF as the mediator (Chen et al., 2020), and (3) the assumption that metabolite levels change dynamically (like rCBF) and thus may carry miRNA-behavior relationships. Others have argued that mediation analysis which includes (a) precedence information (i.e., applying theory to select X, M, and Y) and (b) controls to rule out non-spuriousness (i.e., the secondary hypotheses) allow better understanding of potential mechanisms (Cohen et al., 2003; Shrout and Bolger, 2002). When conditions for (a) and (b) are met, mediation is argued to be statistically mechanistic (Shrout and Bolger, 2002) as opposed to simply correlative (Spencer et al., 2005), and the mediator variable lies statistically in the causal pathway between the two other variables (MacKinnon et al., 2007).

Given the large number of two-way relationships tested prior to mediation analysis and our limited human sample with its potential for skewed distributions and outliers, we applied outlier analysis (Cook’s distance) and permutation-based mediation methods (Fisher, 1935). Permutation, unlike standard statistics, produces a true distribution, thus increasing the overall power of the analysis (see Methods for further elaboration). Given the benefits argued for permutation and/or bootstrapping approaches with limited-sized datasets, we implemented a permutation-based mediation (PMe) approach, using the data itself as a framework for addressing the issues of central limit theorem violation and small-to-moderate sample size (Fritz and MacKinnon, 2007; Kroehl et al., 2020).

This preliminary study utilized PMe to determine if any metabolic pathway, and its compounds, was consistently observed in the majority of mediation relationships with neuroinflammatory miRNAs and VR-based motor control, thus pointing to a potential unifying construct for HAE/concussion research and providing probable blood biomarkers of subconcussive injury (Chu et al., 2019). The statistical framework used for these analyses could be applied broadly to integrate diverse molecular biology measures with human behavior and could provide a translational construct for animal studies.

Results

The presented results integrated three variables (miRNA, metabolites, and VR scores of motor control) collected in a cohort of collegiate American football athletes at two time points: 1) prior to the season (Pre) and 2) immediately following the season (Post). Nine miRNAs previously implicated in HAE-related inflammation (Papa et al., 2019) and four VR-based performance metrics sensitive to subtle concussion-related behavioral deficits (Luria, 1973) were included. Metabolites were selected if (1) they were reported in the top 20 of the random forest analysis, and (2) they were involved in biochemical pathways related to energetics, inflammation, or exogenous consumption and they significantly changed from Pre to Post (Wilcoxon signed-rank p-value < 0.05). This led to the selection of 40 metabolites. Permutation-based mediation was used to assess relationships between miRNAs, metabolites, and VR metrics at both Pre and across-season ( = Post–Pre). MiRNAs, metabolites, and VR metrics observed to have significant mediation/moderation relationships were further regressed against HAEs to test for linear relationships. HAEs implicated in these linear relationships were lastly investigated for mediating or moderating effects with metabolites, miRNAs, or VR measures across-season.

Random forest analysis revealed biochemically important metabolites

Of the 1300+ metabolites analyzed, 20 had the highest mean decrease accuracy (MDA) when distinguishing Pre from Post (Figure S1). Of these metabolites 11/20 (55%) were lipids, 3/20 (15%) were nucleic acids, 3/20 (15%) were xenobiotics, 2/20 (10%) were carbohydrates, and 1/20 (5%) was an amino acid. The metabolite with the highest MDA was 2-hydroxyglutarate (2-HG), followed by three other lipids. All 20 metabolites from the random forest plot were included in subsequent analyses.

Wilcoxon signed-rank analyses revealed across-season metabolite changes

Wilcoxon signed-rank tests were conducted to determine which miRNA, metabolites, and VR scores fluctuated between Pre and Post. There were no significant changes in miRNA levels and VR scores between Pre and Post (p-value > 0.05). Of the measured metabolites, 259 significantly increased or decreased (p-value 0.05) between Pre and Post and 161 survived an FDR correction (q-value ≤ 0.05). Of those 161 metabolites, 40 were selected, including those from the random forest analysis, based on the criteria described in MiRNA quantification and Metabolite quantification (Figure S1 and Table 1). Of the 40 metabolites, 26 increased (negative z-score) and 14 decreased (positive z-score) (Table 1). 48% of the metabolites were lipids and 52% fell in another category (12.5% energy-related, 25% xanthine, 7.5% amino acid, 5% carbohydrate, and 2.5% nucleic acid). Of the lipids, 74% were FAs and 26% were other lipid types. Overall, 35% of the metabolites were FAs and 65% fell in another category. Of the FAs, 79% significantly increased from Pre and Post. Additionally, all TCA-related metabolites (-ketoglutarate, citrate, aconitate, and fumarate) decreased from Pre to Post and all xanthine metabolites (e.g., paraxanthine) increased.

Table 1

Metabolite changes across-season

Metabolite	Super pathway	Sub pathway	p-value	q-value	Z score
corticosteronea	Lipid	corticosteroid	0.0002	0.0028	−3.74
cortisola	Lipid	corticosteroid	0.0001	0.0016	−3.98
Cortisone	Lipid	corticosteroid	0.0447	0.0447	−2.01
Cortoloneglucuronide (cg)	Lipid	corticosteroid	0.0040	0.0194	2.88
Stearidonate	Lipid	long chain polyunsaturated fatty acid	0.0035	0.0194	−2.92
linoleate3n6	Lipid	long chain polyunsaturated fatty acid	0.0013	0.0091	−3.22
sebacatea	Lipid	fatty acid, dicarboxylate	0.0001	0.0016	−3.95
linoleate2n6	Lipid	long chain polyunsaturated fatty acid	0.0097	0.0194	−2.59
Dodecadienoate	Lipid	fatty acid, dicarboxylate	0.0068	0.0194	−2.71
azelatea	Lipid	fatty acid, dicarboxylate	0.0003	0.0036	−3.62
suberatea	Lipid	fatty acid, dicarboxylate	0.0002	0.0028	−3.77
7-hydroxyoctanate	Lipid	fatty acid, monohydroxy	0.0074	0.0194	−2.68
8-hydroxyoctanatea	Lipid	fatty acid, monohydroxy	0.0001	0.0016	−3.86
undecanedioatea	Lipid	fatty acid, dicarboxylate	0.0004	0.0040	−3.53
caproatea	Lipid	medium chain fatty acid	0.0004	0.0040	3.53
heptanoate (7:0)a	Lipid	medium chain fatty acid	0.0008	0.0072	3.35
tridecenedioatea	Lipid	fatty acid, dicarboxylate	0.0003	0.0036	3.65
2-hydroxyglutaratea	Lipid	fatty acid, dicarboxylate	0.0000	0.0000	−4.14
N-palmitoylserine	Lipid	endocannabinoid	0.0009	0.0072	3.34
Alpha-ketoglutarate	Energy	TCA cycle	0.0386	0.0386	2.07
Citrate	Energy	TCA cycle	0.0192	0.0384	2.34
Aconitate	Energy	TCA cycle	0.0068	0.0272	2.71
Fumarate	Energy	TCA cycle	0.0177	0.0384	2.37
Phosphate	Energy	oxidative phosphorylation	0.0051	0.0255	−2.80
Paraxanthine	Xenobiotic	xanthine	0.0220	0.0440	−2.29
1-methylurate	Xenobiotic	xanthine	0.0034	0.0272	−2.93
1,3-dimethylurate	Xenobiotic	xanthine	0.0081	0.0440	−2.65
1,7-dimethylurate	Xenobiotic	xanthine	0.0120	0.0440	−2.51
1-methylxanthine	Xenobiotic	xanthine	0.0009	0.0081	−3.32
5-acetylamino-6-amino-3-methyluracil (aam)	Xenobiotic	xanthine	0.0184	0.0440	−2.36
5-acetylamino-6-formylamino-3-methyluracil (afm)	Xenobiotic	xanthine	0.0123	0.0440	−2.50
O-sulfo-L-tyrosinea	Xenobiotic	chemical	0.0001	0.0011	−3.98
2-isopropylmalatea	Xenobiotic	food component/plant	0.0089	0.0440	2.62
1,2,3-benzenetriol sulfate (2)a	Xenobiotic	chemical	0.0005	0.0050	−3.50
1-carboxyethylphenylalaninea	Amino acid	phenylalanine metabolism	0.0002	0.0002	3.74
prolyl-hydroxy-prolinea	Amino acid	urea cycle; arginine and proline metabolism	0.0000	0.0000	−4.17
1-carboxyethylvalinea	Amino acid	leucine, isoleucine and valine metabolism	0.0002	0.0002	3.74
fructosea	Carbohydrate	fructose, mannose and galactose metabolism	0.0001	0.0001	3.86
mannosea	Carbohydrate	fructose, mannose and galactose metabolism	0.0001	0.0001	−3.95
adenosinea	Nucleotide	polyamine metabolism	0.0333	0.0333	2.13

Metabolite levels between Pre and Post were assessed using the Wilcoxon signed-rank test. Metabolites with significant q-values (<0.05) were reported. Negative z-scores indicate an increase from Pre to Post and positive z-scores indicate a decrease from Pre to Post.

Indicates metabolites that also appeared in the random forest plot. See also Figure S1.

There were significant relationships between miRNAs, metabolites, and VR scores

Linear regressions assessed relations between four VR scores (Balance (Bal), sensory-motor reactivity or Reaction Time (RT), spatial navigation or Spatial Memory (SM), and Comprehensive (Comp, the combination of Bal, RT, and SM)), nine miRNAs (miR-20a, miR-505, miR-3623p, miR-30d, miR-92a, miR-486, miR-195, miR-93p, and miR-151-5p), and 40 metabolites (Table 1). Years of playing experience (YoE) and history of concussion (HoC) were initially considered as covariates but none were significantly correlated with the variables considered for mediation. At Pre, there were nine significant relationships between VR scores and miRNAs, 15 between VR scores and metabolites, and 36 between miRNAs and metabolites (p ≤ 0.05 following Cook’s outlier removal) (Figure S2). Across-season ( = PostPre), there were 12 relationships between VR scores and miRNAs, 20 between VR scores and metabolites, and 31 between miRNAs and metabolites (p 0.05 following Cook’s outlier removal) (Figure S3). In general, both Pre and across-season analyses revealed negative relationships between VR scores and miRNAs, and positive relationships between VR scores and metabolite levels (exceptions: cortoloneglucuronide, corticosterone, adenosine, and tridecenedioate). Linear correlations between miRNAs and metabolites varied directionally, with the majority being negative (79%).

MiRNAs, metabolites, and VR scores formed three-way relationships

As a prerequisite for mediation analysis, only regressions which formed three-way relationships were considered (e.g., relationships with p < 0.05 between miRNAs → metabolites, metabolites → VR scores, and miRNAs → VR scores; see below).

Of the 9, 15, and 36 significant Pre relationships, there were 18 instances where relationships formed a triangular path; of the 12, 20, and 31 across-season relationships, there were 32 triangular paths. In total, there were 50 potential three-way relationships (18 at Pre and 32 across-season). Most three-way relationships included Comp behavior (70%) and 88% included a FA metabolomic measure. Following Cook’s outlier removal, 11 of 18 Pre three-way relationships met significance (p-value < 0.05) for all triangular paths. Across-season, 14 of 32 relationships met this same requirement. Of these 25 significant three-way relationships (Pre + across-season), 92% included a FA metabolomic measure and 88% included Comp. It should be noted that in no case was the SM variable involved in any three-way relationship. In total, there were six miRNAs involved in Pre and across-season three-way relationships: miR-20a, miR-505, miR-92a, miR-151-5p, miR-195, and miR-30d (five at Pre and five across-season). All 25 three-way relationships were tested for mediation.

Mediation analyses revealed significant effects at Pre and across-season

Primary mediation analyses were performed to assess if metabolites (mediator, M) carried the relationship between miRNAs (independent variable, X) and VR scores of motor control (dependent variable, Y). Additionally, secondary mediation analyses were conducted to confirm that miRNAs did not act as mediating variables (all secondary mediation results are presented in the Supplemental Material).

Mediation analysis revealed six significant Pre mediations, out of the 11 tested ( 0.05, Teff 50%; Table 2, bolded and italicized rows). Each significant mediation, including the directionality and significance of each linear regression relationship, is depicted in Figure 1. The plots illustrate the slopes of linear models 1 (X → Y) versus 2 (X + M → Y). In each case, the slope of model 1 was larger than that of model 2, indicating that the metabolite (M) significantly altered the relationship between the miRNA (X) and VR score (Y). The directionality of the linear regression relationships was common across all Pre mediations – negative between miRNA and metabolite, negative between miRNA and VR score, and positive between metabolite and VR score. No secondary mediation analyses met both a significant p-value and a Teff greater than 50% when metabolite was the X, miRNA was the M, and VR score was the Y (i.e., > 0.05, Teff≤ 30%) (Table S1).

Table 2

Permutation-based Pre-season mediation results

Preprimary mediations			Step 1: X predicting Y		Step 2: X predicting M		Step 3: X (with M) predicting Y		Step 3: M (with X) predicting Y		Teff	pSobelperm	N
VR score (Y)	X	M	c	pc	a	pa	c'	pc'	b	pb	%
Comprehensive	miR-20a	2-hydroxyglutarate	−0.531	0.023	−0.529	0.024	−0.256	0.263	0.520	0.032	52	0.002	18
Comprehensive	miR-20a	sebacate	−0.567	0.014	−0.489	0.039	−0.400	0.099	0.343	0.152	30	0.054	18
Comprehensive	miR-20a	8-hydroxyoctanoate	−0.544	0.020	−0.631	0.005	−0.224	0.378	0.506	0.058	59	0.005	18
Comprehensive	miR-20a	undecanedionate	−0.575	0.013	−0.630	0.005	−0.304	0.239	0.287	0.290	32	0.113	18
Comprehensive	miR-505	2-hydroxyglutarate	−0.513	0.021	−0.657	0.002	−0.050	0.819	0.704	0.005	90	0.000	20
Comprehensive	miR-505	8-hydroxyoctanoate	−0.515	0.020	−0.475	0.034	−0.273	0.193	0.509	0.022	47	0.007	20
Comprehensive	miR-92a	8-hydroxyoctanoate	−0.447	0.048	−0.551	0.012	−0.098	0.648	0.634	0.008	78	0.001	20
Comprehensive	miR-151-5p	2-hydroxyglutrate	−0.612	0.005	−0.533	0.019	−0.320	0.110	0.547	0.011	48	0.000	19
Comprehensive	miR-151-5p	8-hydroxyoctanoate	−0.595	0.007	−0.572	0.010	−0.298	0.171	0.547	0.011	50	0.004	18
Balance	miR-151-5p	undecanedionate	−0.578	0.015	−0.575	0.016	−0.213	0.317	0.635	0.008	63	0.001	17
Reaction time	miR-195	8-hydroxyoctanoate	−0.472	0.031	−0.458	0.037	−0.277	0.205	0.426	0.059	41	0.014	21

Cook's outliers from regressions X → Y, X → M, and M → Y were removed to obtain a common set of participants for mediation testing. X = independent variable; M = mediator; Y = dependent variable; Std = standardized beta coefficients for each regression; px = p-value for each regression (a, b, c, and c') at significance level 0.01; permutation-based Sobel p-values = at significance level = 0.05; effect mediated (Teff) is expressed as %; sample size = N. Significant mediations are . See also Figure 1 and Table S1.

Figure 1

Significant Pre-season mediation results

In all analyses, miRNA was X, metabolite was M, and VR score was Y. VR terms were reported as standardized values and adjusted R2 (), p-values (significance level = 0.05), and terms were reported from linear regression analyses (i.e., without common outlier removal that preceded the mediation analyses presented in Table 2). Teff and were reported from each permutation-based mediation analysis.

(A) There was a negative relation between miR-20a and 2-hydroxyglutarate (2-HG), a positive relation between 2-HG and VR composite score (Comp), and a negative relation between miR-20a and Comp. When 2-HG was added to the regression model, the relationship between miR-20a and Comp no longer existed (p-value > 0.05); therefore, 2-HG statistically mediated the relationship ( = 0.002, Teff = 52%). The graph depicts the change in slope between model 1 (X → Y), which plots the slope term for the interaction between miR-20a and Comp, and model 2 (X + M → Y), which plots the slope term for the relationship between miR-20a and Comp when 2-HG was included in the regression model.

(B–F) mimic what has been described in (A).

Across-season mediation analysis revealed eight mediations, out of the 14 tested, meeting our thresholds ( 0.05; Teff 50%; Table 3, bolded and italicized rows). Figures 2A–2H visualizes each across-season mediation, including the directionality of each regression between X and M, M and Y, and X and Y. The plots depict the slopes of linear models 1 (X → Y) versus 2 (X + M → Y). In each case, the slope of model 1 was larger than that of model 2. Mediations depicted in Figures 2A–2C and 2F–2H share the same directionality across regressions: negative between miRNA and metabolite, positive between metabolite and VR score, and negative between miRNA and VR score. In 2D-E, the relations between miRNA (miR-30d) and metabolite (adenosine), as well as metabolite and VR score (Comp), were opposite of those depicted in Figures 2A–2C and 2F–2H. Only one secondary mediation met threshold when metabolite was the X, miRNA was the M, and VR score was the Y (= 0.011, Teff = 68%; Table S2); in this case, the primary mediation analysis also had a significant permutation-based Sobel p-value ( = 0.047), suggesting a more complex relationship.

Table 3

Permutation-based across-season mediation results

Across-season primary mediations			Step 1: X predicting Y		Step 2: X predicting M		Step 3: X (with M) predicting Y		Step 3: M (with X) predicting Y		Teff	pSobelperm	N
VR score (DV)	IV	M	c	pc	a	pa	c'	pc'	b	pb	%
Comprehensive	miR-505	sebacate	−0.579	0.015	−0.605	0.010	−0.289	0.252	0.481	0.066	50	0.008	17
Comprehensive	miR-505	azelate	−0.579	0.015	−0.643	0.005	−0.302	0.265	0.432	0.118	48	0.017	17
Comprehensive	miR-505	suberate	−0.579	0.015	−0.638	0.006	−0.341	0.219	0.375	0.178	41	0.035	17
Comprehensive	miR-505	8-hydroxyoctanoate	−0.524	0.037	−0.526	0.037	−0.217	0.353	0.584	0.023	59	0.007	17
Comprehensive	miR-30d	sebacate	−0.586	0.008	−0.595	0.007	−0.298	0.196	0.485	0.043	49	0.003	19
Comprehensive	miR-30d	suberate	−0.551	0.022	−0.518	0.033	−0.309	0.199	0.467	0.061	44	0.009	17
Comprehensive	miR-30d	heptanoate	−0.586	0.008	−0.618	0.005	−0.328	0.183	0.416	0.097	44	0.009	19
Comprehensive	miR-30d	adenosine	−0.474	0.047	0.550	0.018	−0.191	0.434	−0.515	0.047	60	0.001	18
Comprehensive	miR-92a	suberate	−0.505	0.039	−0.622	0.008	−0.161	0.537	0.553	0.048	68	0.004	17
Comprehensive	miR-92a	heptanoate	−0.479	0.033	−0.523	0.018	−0.240	0.299	0.456	0.058	50	0.006	20
Comprehensive	miR-195	heptanoate	−0.491	0.033	−0.667	0.002	−0.247	0.386	0.366	0.206	50	0.029	19
Comprehensive	miR-151-5p	8-hydroxyoctanoate	−0.616	0.011	−0.653	0.006	−0.279	0.286	0.516	0.060	55	0.007	17
Comprehensive	miR-151-5p	heptanoate	−0.569	0.017	−0.585	0.014	−0.371	0.169	0.338	0.207	35	0.047	17
Reaction time	miR-30d	adenosine	−0.473	0.055	0.529	0.029	−0.200	0.420	−0.516	0.050	58	0.001	17

Cook's outliers from regressions X → Y, X → M, and M → Y were removed to obtain a common set of participants for mediation testing. X = independent variable; M = mediator; Y = dependent variable; Std = standardized beta coefficients for each regression; px = p-value for each regression (a, b, and c) at significance level 0.1; permutation-based Sobel p-values = at significance level = 0.05; effect mediated (Teff) is expressed as %; sample size = N. Significant mediations are . See also Figure 2 and Table S2.

Figure 2

Significant across-season mediation results

In all analyses, miRNA was X, metabolite was M, and VR score was Y. VR terms were reported as standardized values and adjusted R2 (), p-values (significance level = 0.05), and terms were reported from linear regression analyses (i.e., without common outlier removal that preceded the mediation analyses presented in Table 3). Teff and were reported from each permutation-based mediation analysis.

(A) There was a negative relation between miR-505 and sebacate, a positive interaction between sebacate and Comp, and a negative relation between miR-505 and Comp. When sebacate was added to the regression model, the relationship between miR-505 and Comp no longer existed; therefore, sebacate statistically mediated the relationship (= 0.008, Teff = 50%). The graph depicts the change in slope between model 1 (X → Y), which plots the slope term for the interaction between miR-505 and Comp, and model 2 (X + M → Y), which plots the slope term for the relationship between miR-505 and Comp when sebacate was included in the regression model.

(B–H) mimic what has been described in (A).

In total, from Pre- and across-season mediation analysis, there were 14 significant mediations (< 0.05, Teff 50%) comprised of seven metabolites (2-HG, 8-hydroxyoctanoate (8-HOA), undecanedioate (UND), sebacate, suberate, heptanoate, adenosine), six miRNAs (miR-20a, miR-505, miR-92a, miR-151-5p, miR-195, and miR-30d), and three VR scores (Comp, Bal, and RT). Of the seven metabolites, six were FAs. Five of the six FAs (2-HG, 8-HOA, UND, sebacate, and suberate) increased from Pre to Post, while heptanoate, a FA, and adenosine, a nucleoside, decreased (Figure S4). Interestingly, five variables–Comp, 8-HOA, miR-505, miR-92a, and miR-151-5p–were observed in both Pre and across-season mediations.

HAEs showed relationships with miRNA and metabolite levels and moderated metabolite change

Three methods were used to assess HAE relationships: (1) extrapolation of athlete-specific Post HAE metrics to model Pre VR scores, miRNA levels, and metabolite levels; this analysis is purely exploratory and should only be considered for setting up future hypotheses; (2) assessment of how Post HAE metrics affect across-season changes in VR scores, miRNA levels, and metabolite levels; this analysis was not exploratory as it tested expected relationships with HAEs; and (3) mediation and moderation analyses to model directional relationships across variables.

In this exploratory framework, analysis of HAEs against Pre measures revealed six significant relationships between HAEs and miRNAs (5/6 related to HAEs above 80G) and two relationships between HAEs and metabolites (2/2 related to HAEs above 25G) (Table S3). In testing non-exploratory hypotheses of HAE relationships with across-season molecular measures, we found six significant relationships between HAEs and omic measures (5/6 related to HAEs above 25G). All relationships were positive (+) indicating that as HAEs increased miRNA and metabolite levels were increasingly higher at Post as compared to Pre (Table S3, Figure S5). Relationships with miRNAs were related to normalized or averaged HAEs (i.e., aHAE25G and aHAE80G) while relationships with metabolites (sebacate and suberate) were specifically related to HAEs exceeding 25G (cHAE25G and aHAE25G). Significant HAE regressions are plotted in Figure S5 and their regression terms with across-season mediations are depicted in Figures 2A–2I.

For the third analysis involving mediation and moderation, mediation showed no significant effects. Moderation was conducted on metabolites implicated in across-season mediation to assess whether the interaction of HAEs and Pre metabolite levels could directionally model Post metabolite levels. In total, there were five significant moderations ( 0.05; 0.05); 4/5 (80%) implicated HAEs at, or exceeding 25G and 1/5 implicated sessions (Table S4). 3/5 moderation findings included the metabolite suberate, a dicarboxylic FA. Following a Bonferroni correction, 3/5 of the moderations remained statistically significant ( 0.01). It is important to note that moderation moves past simple association, indicating that the interaction between the X and moderator (Mo) predict, or directionally model, the Y.

Post-hoc miRNA-gene network analysis

To improve the interpretability of our findings, network analyses were performed to establish relationships between the nine miRNAs and genes involved in mitochondrial functions (i.e., fatty acid oxidation, the TCA cycle, the electron transport chain (ETC), oxidative phosphorylation, general aerobic respiration) and neuroinflammation. When brain tissue was selected, 43 genes involved with mitochondrial processes appeared in the nodal network with miRNAs. It should be noted that only three (miR-20a, miR-9-3p, miR-30d) of the nine miRNAs of interest existed in the miRNet database for human brain tissue. This could, in part, be due to insufficient research of these miRNAs in human brain tissue. Of the 43 genes, 33 had direct nodal relationships with a given miRNA (Table S5, Figure S6). When tissue was not specified, there were 87 genes in the miRNA-gene network. Of the 87, 61 had direct nodal relationships with a given miRNA (Table S6, Figure S7).

Out of the 37 neuroinflammatory-related genes investigated, 11 showed a nodal relationship with one or more of the three miRNAs when human brain tissue was selected as the target tissue (Table S7, Figure S8). For the described network analyses, it should be noted that although not every gene had a direct, one-nodal, relationship with a given miRNA, the miRNAs still existed in the overall miRNA-gene network.

Discussion

This preliminary study tested the hypothesis that collegiate American football athletes, exposed to repetitive HAEs, would exhibit significant PMe effects between neuroinflammatory-related miRNAs, energy-related metabolites, and VR-based motor control metrics, where one measure was consistently the mediator (M) and carried the relationship between the independent variable (X) and dependent variable (Y). This hypothesis was confirmed with five general findings. (1) Metabolomic analysis identified 40 compounds that significantly changed across-season. (2) There were 14 significant mediations where miRNA was X, metabolite was M, and VR score was Y (i.e., metabolites always sat in the statistically causal pathway between transcriptome and behavior). (3) Three miRNAs (miR-505, miR-92a, and miR-151-5p) and one metabolite (8-HOA) appeared in both Pre and across-season mediations. (4) All seven metabolomic compounds implicated in mediation relationships are involved with some aspect of energy metabolism. Specifically, five of the seven are involved in FA oxidation, a sixth compound is involved with diminishing mitochondrial respiration and respiration-coupled ATP production in cancer cells (i.e., 2-HG), and a seventh is a core structural unit of ATP (i.e., adenosine). (5) miRNAs and metabolite levels were related to accumulated HAEs, and the interaction between HAEs and Pre metabolite levels predicted, or directionally modeled, Post metabolite levels. Altogether, these findings point to a shift in mitochondrial metabolism, away from mitochondrial function, and consistent with known mitochondrial disorders (Longo et al., 2016; Wajner and Amaral, 2016).

Metabolite changes indicative of mitochondrial dysfunction and altered metabolism

Metabolites of the TCA cycle consistently decreased across-season while monohydroxy (7-HOA, 8-HOA) and dicarboxylic (sebacate, suberate, UND) FAs increased and were consistently implicated in PMe analyses. This observation suggests a dysfunction in energy metabolism, specifically related to FA oxidation. Typically, long- and medium-chain FAs are processed in peroxisomes, small organelles that oxidize FAs via -, -, and -oxidation (Lodhi and Semenkovich, 2014). -oxidation replaces the methyl terminus with a hydroxy group (producing a monohydroxy FA) which is then further oxidized into a carboxy group (dicarboxylic FA). Once in dicarboxylic form, the FA can be further processed via -oxidation in both the peroxisome and mitochondria (Bartlett and Eaton, 2004). Depending on FA length and oxidation location(s), different sets of enzymes metabolize it into a smaller, functional metabolite, namely acetyl-CoA. Acetyl-CoA enters the TCA cycle to produce energy (GTP) and energy-storing (NADH and FADH2) metabolites (Fernie et al., 2004). NADH and FADH2 are used in the electron transport chain to produce ATP – the most fundamental energy source. Here, increased levels of medium-chain monohydroxy and dicarboxylic FAs, and decreased levels of TCA metabolites, suggest that these FAs cannot be further oxidized into acetyl-CoA for normal TCA cycle respiration (Figure 3). This is relevant in that increased serum levels of 7-HOA, 8-HOA, suberate, and sebacate are observed in patients with medium-chain acyl-CoA dehydrogenase deficiency (MCADD) – a genetic mitochondrial disorder (Bodman et al., 2001; Lee et al., 2005; Merritt and Chang, 1993). Similar to our findings, another feature of this disorder is increased serum levels of acyl-carnitine derivatives with 6, 8, and 10 carbons. These metabolites were increased in our cohort, however with FDR-corrected q-values > 0.05, so a larger cohort is needed to determine the validity of this observation (data not shown). Lastly, heptanoate, a monohydroxy FA, was observed to decrease, rather than increase, across-season. Increased heptanoate, via metabolism of triheptanoin, has been reported as beneficial in the treatment of -oxidation related metabolic disorders (Marin-Valencia et al., 2013), suggesting that the observed decrease in heptanoate could be detrimental.

Figure 3

Summary of the observed metabolic disturbances

There were significant increases in medium-chain monohydroxy and dicarboxylic fatty acids (FAs) (suberate, sebacate, UND, and 8-HOA) from Pre to Post. Increases in these FAs are associated with genetic disorders related to impaired -oxidation, which can result in an accumulation of medium-chain FAs that cannot be further oxidized into smaller, functional, metabolites, such as acetyl-CoA. Acetyl-CoA is a critical input for the TCA cycle – a major source of energy-rich molecules that are fed into further energy-producing processes (e.g., electron chain transport system). Here, TCA-related metabolites (citrate, aconitate,-KG, fumarate, and malate) all decreased, suggesting a problem with the initial step of the cycle (i.e., lack of acetyl-CoA). Additionally, there were alterations in energy-rich molecules such as adenine, adenosine, nicotinamide, and phosphate, suggesting a state of energy imbalance. Lastly, 2-HG increased. Regardless of its role as a known oncometabolite, its increase suggests a state of oxidative stress. Together, it is suggested that there are dysfunctional -oxidative mitochondrial processes in this cohort of collegiate football athletes leading to subsequent issues with energy production.

Taken together, the observed increase in medium-chain monohydroxy and dicarboxylic FAs, decrease in heptanoate, and concurrent decrease in TCA metabolites suggest an impairment in -oxidation, potentially stemming from HAE-induced mitochondrial dysfunction, which has been reported in several TBI-related studies (Fischer et al., 2016; González-Domínguez, 2016; Hiebert et al., 2015; Signoretti et al., 2010; Vagnozzi et al., 1999).

Other across-season metabolic changes suggest energy-related dysfunction

2-Hydroxyglutarate (2-HG) increased from Pre to Post. 2-HG is dicarboxylic FA generated from -ketoglutarate (i.e., oxoglutaric acid), a component of the TCA cycle, via the mutated IDH gene (Ye et al., 2018). 2-HG has been implicated in tumor suppressor inactivation and oncogenic activation (Dang et al., 2009; Intlekofer et al., 2015), and has been associated with oxidative stress and poorer outcomes following brain injury (Latini et al., 2003; Magalhães et al., 2019). L-2-HG structurally resembles glutamate and thus may affect neurotransmission (Farias et al., 2012). The observed increase in 2-HG is likely linked to oxidative stress resulting from repetitive HAEs and may even alter neurotransmission and subsequent network connectivity (Manning et al., 2017; Zhu et al., 2015).

A decrease in adenosine was also observed at Post. Adenosine is a potent vasodilator that increases rCBF, thereby decreasing energy demands. It is also implicated as a robust neuroprotector and has been reported to suppress ATP release (Laketa et al., 2015). Adenosine dysregulation has been reported following TBI and is linked with numerous TBI-related pathologies and comorbidities (Lusardi, 2009). In particular, severe TBI is associated with increased adenosine levels, which is hypothesized to play a neuroprotective role in humans (Headrick et al., 1994) but is associated with a loss of neuroprotection in rats (Cui et al., 2013). These and other findings suggest an interplay between protective effects of adenosine receptors, A1 and A2A, and further research is required to elucidate time course effects of adenosine across the full spectrum of TBI. Although we observed an overall decrease in adenosine when using nonparametric statistical tests (i.e., median comparison), adenosine did increase overall when parametric tests were applied (e.g., t-test which compares mean). This point can be observed in Figure S3 which depicts high levels of adenosine at Post. This observation emphasizes the need for a larger sample size to truly understand how adenosine levels are affected in these athletes.

Another energy-related metabolite, phosphate, increased and may be associated with reduced ATP synthesis, ATP depletion, and/or ATP hydrolysis. In cases of ischemia, mitochondria hydrolyzed ATP more frequently (Classen et al., 1989) and ATP synthesis was impaired in animal models of brain injury (Pandya et al., 2019). Impaired ATP synthesis may be related to depleted energy stores and a subsequent lack of energy sources. Contrary to these observations, serum phosphate levels decreased in patients with severe TBI (Porter et al., 2017), yet football athletes often experience repetitive subconcussive HAEs. More research is required to elucidate the role of phosphate in the full TBI spectrum.

Lastly, all xanthine metabolites increased at Post suggesting either (a) an increase in caffeine consumption (Bonati et al., 1982), (b) activation of xanthine oxidases (Solaroglu et al., 2005), or (c) a combination of both (Signoretti et al., 2010). Caffeine has various anti-inflammatory and neuroprotective effects in both animals and humans, especially in relation to neurodegenerative disease and TBI (Kolahdouzan and Hamadeh, 2017). The timing of caffeine administration may be critical; however, as was seen in a study where caffeine administered prior to TBI was beneficial while caffeine administered after TBI was detrimental (Boison and Lusardi, 2012). Increased caffeine consumption in this cohort could be the result of school-related stress and lack of sleep and/or compensation for HAE-induced brain damage and metabolic dysfunction. Given the current data, it is not possible to identify the exact reason(s) for the increased levels of xanthines and further research is required. However, it should be noted that previous studies have reported caffeine consumption in contact athletes with the main motive being regaining energy lost during exercise (Badaam, 2013).

Metabolites statistically mediate the relationship between neuroinflammatory miRNAs and behavior

This study used what some call an “exact inference” framework (Ernst, 2004), based on integration of permutation statistics with mediation analyses (PMe), allowing for a moderately sized cohort to be studied (N < 25). Using PMe, it was found that metabolomic measures were a required variable for the effect of miRNA on behavior in all significant cases. This analysis confirmed the absence of confounds, such as miRNA levels mediating the relationship between metabolite and behavior; in such contexts, others have argued that mediation is statistically mechanistic (Shrout and Bolger, 2002) as opposed to purely associative (Spencer et al., 2005). Per MacKinnon (MacKinnon et al., 2007), “mediating variables…transmit the effect of one variable to another variable”, indicating that the mediator lies in a statistically causal pathway between two other variables (i.e., miRNA → metabolite → behavior). The 14 mediation effects observed were all greater than 50%, indicating that metabolomic measures filled an important role carrying the effect of miRNA levels to motor control. It is likely relevant that this relationship of neuroinflammatory miRNAs and behavioral performance was carried by a set of metabolites involved in some aspect of energy metabolism, in particular FA metabolism.

Previous studies have reported increased levels of neuroinflammatory miRNAs prior to contact play with minimal change in miRNA levels across season (Chen et al., 2020; Papa et al., 2019). MiRNAs in this panel have known roles in cancer progression (Che et al., 2019; Hsu et al., 2016; Huang et al., 2018; Kobayashi et al., 2012; Kohram et al., 2018; Li et al., 2011, 2016; Ling et al., 2019; Lv et al., 2018; Mogilyansky and Rigoutsos, 2013), inflammation and inflammatory diseases (Paraskevi et al., 2012; Reddycherla et al., 2015; Wu et al., 2011; Zheng et al., 2015), and neurodegeneration (Basak et al., 2016; Denk et al., 2015; Jee et al., 2012b; Khoo et al., 2012; Liu et al., 2017). These miRNAs have also been shown to target genes that have neuroinflammatory functions and roles in diseases with neuroinflammatory pathophysiology (Cao et al., 2017; Lisha Chang et al., 2020; Gao et al., 2019; Gayen et al., 2020; Hu et al., 2019; Jee et al., 2012a; Jiang et al., 2018; Ma et al., 2014; Mao et al., 2019; Martinez and Peplow, 2017; Naqvi et al., 2016; Navaderi et al., 2019; Peng and Ying, 2014; Shi et al., 2018; Sinha et al., 2012; Sujith et al., 2018; Thangaraj et al., 2021; Wang et al., 2021; Xie et al., 2020; Yang et al., 2018; Yao et al., 2014; Zhao et al., 2017; Zhong et al., 2021; Zhou et al., 2021). Post-hoc miRNA-gene network analyses further solidified the involvement of these nine miRNAs in mitochondrial processes and neuroinflammation (see Tables S6–S8, Figures S6–S8).

It should also be noted that three of the miRNAs identified herein (miR-20a, miR-92a, and miR-30d) were involved in mediation relationships with regional brain perfusion changes in the basal ganglia of the same subjects, grounding the current findings in clear neurophysiological abnormalities in brain regions that are important for motor control and spatial behavior (Chen et al., 2020). Behavioral changes alone are rarely observed and difficult to replicate in studies of contact athletes without diagnosed concussion (Bailes et al., 2013; Nauman et al., 2020).

At Pre, there were six significant mediations where metabolites (2-HG, 8-HOA, and UND) statistically mediated the relationship between chronically elevated miRNAs and behavior. The fact that FAs consistently mediated miRNA-motor control relationships suggests that dysfunctional -oxidation may underlie subtle alterations in behavior. Moreover, the fact that these relationships were observed prior to contact practice suggests chronic metabolic dysregulation and neuroinflammation that may persist long after season commencement (Figure 4).

Figure 4

Shift in metabolism and its relationship to neuroinflammatory miRNAs and complex behavior

We observed changes indicative of mitochondrial distress as evidenced from both the accumulation of medium-chain FAs and decreases in TCA-related metabolites. Mitochondrial dysfunction can lead to numerous physiological disturbances, some of which were observed in the present study: 1) oxidative stress (i.e., increased levels of 2-HG), 2) impairment in beta-oxidative processes (i.e., increased levels of medium-chain FAs), and 3) increased metabolic demands (i.e., decreased TCA and energy-related metabolites). Together, these processes may be related to the observed elevation in neuroinflammatory-related miRNA molecules (specifically miR-20a, miR-505, miR-151-5p, miR-30d, miR-92a, and miR-195). In fact, these miRNAs were significantly correlated with the metabolites shown in Figure 3. In addition, the metabolites were shown to mediate the relationship between elevated miRNA levels and VR-based motor control. This complex relationship may explain why obvious behavioral changes in subconcussed athletes are not routinely observed, but how repetitive, long-term exposure to HAEs, chronic elevation of neuroinflammatory miRNAs, and acute, but deleterious changes in energy metabolites could result in behavioral disturbances later in life.

Across-season, there were eight significant mediations. Again, the majority (6/8) of mediating metabolites were FAs which have been implicated in dysfunctional mitochondrial -oxidation by others (González-Domínguez, 2016; Merritt et al., 2018; Wajner and Amaral, 2016; Wanders et al., 2011). Additionally, adenosine was observed in two mediation relationships. The observation of significant mediations using across-season () measures indicated that small alterations in miRNAs, from an elevated Pre-season level (Papa et al., 2019), affect motor control in a negative fashion via changes in metabolite levels.

The consistency of (1) metabolomic measures being mediators between miRNA and behavior and (2) the observed metabolites relating to altered FA oxidation, suggests this pathway may be important for explaining how HAEs may diminish motor control, an indicator of neurological function. It should be noted that the majority of mediations involved the Comp score; this is likely because Comp is the aggerated form of the other three scores (Bal, RT, and SM), thereby increasing its reliability as compared to Bal, RT, or SM alone. The metabolites identified in this study, and their relative changes, support the hypothesis of metabolic dysregulation in these athletes, potentially related to dysfunctional mitochondria and subsequent dysfunction in cellular respiration (i.e., TCA cycle and downstream energy-producing processes). It should be noted that this type of multi-omic approach might generally work well in translational human and animal studies, particularly as there is no face validity issue of the studied human clinical model.

HAEs were related to changes in metabolites and miRNA

Two sets of findings were observed: (1) using standard linear regression, miRNAs, metabolites, and VR-based motor control measures were related to accumulated HAE exposures over the football season and (2) moderation analyses revealed that the interaction between HAEs and Pre-season metabolite levels successfully ( 0.05, ) modeled Post-season metabolite levels; this may indicate that HAEs played a critical role in metabolite change across-season. Results from standard regression analyses revealed that all across-season relationships were positive, indicating that increased HAE exposure was related to larger increases in neuroinflammatory miRNA levels and increased metabolite measures. Pre results indicate an opposite trend where higher extrapolated HAEs were related to relatively lower levels of miRNA, raising the possibility that lower levels of neuroinflammatory-related miRNAs may be related to more HAE exposures from the previous season of play. However, these results are purely explorative because HAEs from the prior season were not recorded, and should only be considered for future hypothesis testing (however, see (Alosco et al., 2017) for a study that extrapolates preseason measures to other measures taken during the season).

Previous research, specifically in the neuroimaging field, has demonstrated the role head impacts play in neurophysiological changes, such as altered connectivity (rs-fMRI and fMRI), decreased brain volumes (T1-and T2-weighted MRI), neurochemical alterations (MRS), altered CBF (perfusion imaging), and axonal injury (DTI) (Bernick et al., 2020; Chen et al., 2020; Nauman et al., 2020; Vagnozzi et al., 2010; Wang et al., 2008). In this study, multi-omic blood biomarkers (metabolites and miRNAs) corresponded with HAE exposure. Interestingly, no computational behavior (i.e., VR measures) was related to HAEs. This absence of relationships may be due to inherent limitations of the sensors, the moderate sample size, and/or the fact that HAEs were not recorded during competition – an unavoidable limitation as to not disturb game preparations.

Limitations of the study

There are several limitations to this study, one being the absence of intermediate sample collection time points between Pre and Post. Blood collections and VR tests were only conducted twice – once prior to contact practice (Pre) and once directly following the end of the regular competitive season (Post). To observe more transient changes in these metrics, or more complex dynamics in the structure of change, it would be beneficial to collect data at more time points during and after the season. This study also lacked age-matched, non-contact athlete controls with similar nutritional plans, body types, athletic regimens, and lifestyle management programs to the studied cohort. However, individual baseline measures taken at Pre served as appropriate comparison group for Post measures given these players were under strict nutritional, workout, and lifestyle management programs at each collection time point. Given that participants follow a strict dietary regimen before and during the season, data surrounding diet and energy demands were not collected. However, this is an important domain for future research as past studies on athlete diets have noted changes in fat composition, energy availability, and metabolic profiles (e.g. (Al-Khelaifi et al., 2018; Binkley et al., 2015; Cole et al., 2005; Hinton et al., 2004; Kelly et al., 2020; Logue et al., 2017)), making it important to tease apart metabolomic changes in response to diet, exercise, and head injury. A significant caveat to this work is that it only involved male participants; future studies need to also incorporate females and include more sophisticated metrics of sex steroid cycles.

The moderate sample size (N = 23) was also a considerable limitation, owing to the preliminary nature of these findings. We did, however, apply multiple methods to account for this, i.e., multiple comparisons correction and permutation. Although we applied a multiple comparisons correction, there is always a possibility that Type I error still exists given our moderate sample size. Future work with larger datasets is necessary. Further, we incorporated a rigorous permutation-based statistical approach to avoid the required assumptions regarding data distributions in smaller samples. Permutation, unlike standard statistics, produces a true distribution, thus increasing the overall power of the analysis. Permutation statistics only require one assumption as opposed to the many that are required for standard parametric and non-parametric statistics; the violation of any of these assumptions leads to higher false positives. Collection of many types of detailed data in humans can lead to smaller samples, which is less likely to violate the central limit theorem if inference testing follows a permutation framework, specifically of the data itself (Fisher, 1935). Permutation-based statistics, that permute the actual data used in mediation, produce a meaningful distribution for exact inference (Ernst, 2004). The small sample size and nonparametric nature of this data can also impact our observations–namely, that there was an overall decrease in adenosine when comparing median–but given a larger dataset, it is probable that a larger number of subjects would exhibit high levels of adenosine at Post.

For this study, and others using mediation, it is important to acknowledge that there is a divergence in the literature about using cross-sectional data, change measures, and temporal data with mediation (Cole and Maxwell, 2003; Gollob and Reichardt, 1987; Kraemer et al., 2002), and both viewpoints have been published extensively in peer-reviewed neuroscience and physiology journals. We followed a Pearl-MacKinnon viewpoint (Lockhart et al., 2011; Mackinnon, 2008; MacKinnon et al., 2007, 2012; Pearl, 2009, 2012, 2014), although the MacArthur viewpoint (Cole and Maxwell, 2003; Gollob and Reichardt, 1987; Kraemer et al., 2002, 2008; Maxwell and Cole, 2007) is also important, and future work will incorporate the latter. Lastly, this study used peripheral blood for metabolomic and miRNA quantification, and it could be argued that they are not direct brain measures. This is a recognized limitation; however, many studies have shown breakdown of the blood-brain-barrier in cases of repetitive mild TBI which could lead to leakage of larger metabolites (namely, lipids) into peripheral circulation (Marchi et al., 2013; Sulhan et al., 2020; Weissberg et al., 2014). Further work is needed to translate these findings with animal models where metabolites/miRNA can be collected from specific tissues.

Conclusions

This preliminary study describes significant relationships between energy-related metabolites, neuroinflammatory miRNAs, and VR-based motor control in collegiate athletes Pre and across-season. The majority of mediation findings involved fatty acids (2-HG, 8-HOA, UND, sebacate, suberate, and heptanoate), which increased across-season. In parallel, TCA metabolites significantly decreased across-season and HAEs showed both correlative relationships with metabolite and miRNA levels across-season, as well as moderating relationships where the interaction between HAEs and Pre metabolite levels directionally modeled Post metabolite levels – indicating HAEs were critical in the observed relationships. Given the context of elevated neuroinflammatory miRNAs and mitochondrial dysfunction, these findings suggest a state of chronic HAE-induced neuroinflammation. The consistent mediation findings suggest a model where metabolites sit in the statistically causal pathway between transcriptome and behavior. Integrative permutation-based mediation approaches used in this study can be applied to other integrative studies of human disease and compliments translational in vitro and in vivo research (Nestler and Hyman, 2010). Together, these findings (1) suggest that metabolic (e.g., 2-HG) and miRNA (e.g., miR-30d, miR-505, miR-151-5p, etc.) disturbances discovered in cancer studies could be a primary focus for head injury research, (2) highlight potential biomarkers for longitudinal assessment of head injury, or biomarkers that implicate head injury in the absence of clinically diagnosable symptoms, and (3) point to mitochondrial dysfunction as a central organizing principle for repetitive head injury in humans.

STAR★Methods

Key resource table

Resource availability

Lead contact

Further information and requests for resources should be directed to and will be fulfilled by the lead contact, Hans Breiter (h-breiter@northwestern.edu).

Materials availability

This study did not generate new unique reagents.

Experimental model and subject details

Human

Twenty-three male collegiate American football athletes (mean age = 21 2 years) participated in this study. The selected athletes were all seasoned starters who experience, and have experienced, a high frequency of HAEs throughout the athletic season (i.e., 16/23 athletes were non-speed linemen who sustain a high volume of impacts during practices and games) (Lee et al., 2020; Lehman, 2013). Written informed consent was obtained from each participant in accordance with the Pennsylvania State University Institutional Review Board and in alignment with the Declaration of Helsinki. Athletes were excluded if they had any known neurological disorder, any lower body injury that affected balance, or if they were unable to participate in the football season.

Method details

Data collection

All data were collected between July 2015 and December 2015 at Pennsylvania State University. Demographic information was obtained from each participant and confirmed by a team physician (Table S4). Blood samples were taken prior to any contact practices (Pre) and within one week of the last regular season game (Post). Our scientific construct specifically used the same participants before and after a season of play who were under the same nutritional program, workout program, and lifestyle management program at both experimental time points; selecting matched controls was infeasible given players’ strict regimens. Pre measures thus served as an appropriate comparison group for Post measurements. None of the athletes had a diagnosed concussion in the nine months preceding Pre data collection. Only one participant was taken out of seasonal play due to a significant peripheral injury that occurred early on in the season; they were retained in the study cohort. Blood samples were prepared and sent out for metabolite and miRNA quantification. Concurrent with blood collection timepoints, athletes also completed virtual reality (VR) testing.

Virtual reality (VR) testing of motor control

VR testing was first described and validated by (Slobounov et al., 2006) and Teel et al. (2016), Teel and Slobounov, (2015). The software used to display the virtual reality animations was developed and provided by HeadRehab, LLC (Chicago, IL). The HeadRehab Performance Test Software allows for the use of a range of interactive devices and depending on the system used, a variety of options are viable. These different systems have all been validated against one another. Overall, the subject will interact with the HeadRehab SideLine Performance Test Software modules via interactive devices or motion tracking devices. This technology was used to create “moving room” experiments with two conditions: (1) the virtual room could move as a whole structure or (2) separate components of the room could move in isolation (e.g., only the front wall moves).

In this study, a 3DTV system was used. A laptop was connected to the TV with an HDMI cord and the Sideline v8.1 test TV program was used. Display resolution was set at 1920 x 1080 (high definition, 1080p) with a 16:9 aspect ratio for use with the 3DTV system. Graphics without Stereo Effect were also used. A VCube head-tracking device and MoBar interactive device (in this case an Xbox controller) were used for interaction with the 3DTV. The VCube uses acceleration and translation to track the subject’s head movement. The MoBar allows the subject to navigate and interact with the software. The program has a subprogram built in (Start Device Manager) that is required to be run to ensure that the devices are connected and functioning properly.

Once the athlete was properly set up, they underwent three virtual reality tasks: 1) balance (Bal), 2) reaction time (RT), and 3) spatial memory (SM). The scores from each test were normalized and combined to produce a comprehensive score (Comp). These tests were based off initial findings from Dr. Alexander Luria and have been validated to detect functional abnormalities following mild traumatic brain injury (Luria, 1973; Slobounov et al., 2006).

The Bal module tested an athlete’s ability to maintain posture with a changing virtual environment (reported and validated by Teel and Slobounov, (2015)). Each visual scene was 30 s in duration and athletes were instructed to hold a tandem Romberg position for each trial (Romberg, 1853). For the first baseline trial, athletes were instructed to remain as still as possible while there were no changes to the visual field. The next six trials consisted of changes to the virtual environment in one of three directions (yaw, pitch, or roll) via two methods: (1) whole-room forward-backward oscillations with 18 cm displacement at 0.3 Hz or (2) whole-room lateral rolls at 10–30 degrees and 0.3 Hz. Deviances from baseline in the yaw, pitch, and roll directions, were recorded and interpolated (scale 0–10) such that a higher score indicated better task performance (10 being the best).

The RT module tested the time it took for each athlete to adjust their posture in the direction of an altered virtual environment. Athletes stood shoulder-width apart with their hands on their hips and were instructed to move their body with the direction of the changing virtual environment (along the anterior-posterior axis at 0.2 Hz). Randomly, the room would shift to shifts in the medial-lateral axis. When this occurred, athletes were asked to bend at the waist in the same direction as the virtual environment (i.e., bend left or right). Both the time and direction of the room shifts were randomized. Response times, in milliseconds, to the abrupt changes in the virtual environment were retrieved using accelerometers and results were reported as latency to shift the body. Five total trials were completed, with the first trial being a practice trial and the remaining four trials being used in the score calculation. The measured RT and errors in anticipation were calculated, interpolated, and concerted into a whole-body RT score ranging from 0 to 10. Higher scores were indicative of better performance.

The SM (virtual corridor) module tested an athlete’s ability to recall a virtual pathway. Athletes were shown a randomized virtual pathway (to avoid practice effect, which in fact was documented and published in Slobounov & Sebastianelli (Slobounov et al., 2014) with multiple turns leading to a door, followed by a return trip. Then, athletes were asked to repeat the pathway from memory using the remote. It should be noted that athletes were allowed to practice using the remote prior to any trials. If the athlete navigated the path correctly, testing was finished. If the athlete navigated the path incorrectly, the computer would replay the desired route and the athlete was given three total chances to correctly navigate. Outcomes were reported as (1) the average time of task completion and (2) the number of errors made. A score between 0 and 10 (10 being the best) was calculated based on the average time a participant stayed on the correct path. For each error made, a 3.33-point reduction was applied. If the athlete was unsuccessful after three attempts, their score was zero.

The Comp score was calculated by combining the three module scores (Bal, RT, and SM) into a ten-point scale (0 = worst, 10 = best).

Serum extraction

Five mL of venous blood were collected from each athlete at Pre and Post. Samples were placed in a serum separator tube, allowed to clot at room temperature, and then centrifuged. Serum was extracted from each tube and pipetted into bar-coded aliquot tubes. Serum samples were stored at −70°C until they were transported to 1) a central laboratory for blinded miRNA batch analysis (Papa et al., 2019) and 2) Metabolon (Morrisville, NC, USA) for blinded metabolite quantification.

Metabolite quantification

Serum was sent to Metabolon (Morrisville, NC, USA) for metabolite quantification, following well-validated and established procedures. Upon arrival, samples were assigned a unique identifier via an automated laboratory system and stored at −80°C. Samples were prepared for subsequent analyses using an automated MicroLab STAR® system (Hamilton Company, Reno, NV, USA). Proteins were precipitated out of each sample using methanol and a shaker (Glen Mills GenoGrinder2000), and centrifuged. The resulting extract was then divided into five fractions for various analyses: 1) Two fractions for analysis by two separate reverse phase (RP)/UPLC-MS/MS methods with positive ion mode electrospray ionization (ESI), 2) one for analysis by RP/UPLC-MS/MS with negative ion mode ESI, 3) one for analysis by HILIC/UPLC-MS/MS with negative ion mode ESI, and 4) one reserved for backup. To remove organic solvent, samples were briefly placed on a TurboVap® (Zymark); samples were then stored under nitrogen prior to analyses. Serum metabolites were quantified using Ultrahigh Performance Liquid Chromatography-Tandem Mass Spectroscopy (UPLC-MS/MS). All methods utilized a Waters ACQUITY UPLC and a Thermo Scientific Q-Exactive high resolution/accurate mass spectrometer interfaced with a heated electrospray ionization (HESI-II) source and Orbitrap mass analyzer operated at 35,000 mass resolution. The sample extract was dried and reconstituted in solvents compatible to each of the listed analyses. Each reconstitution solvent contained a series of standards at fixed concentrations to ensure injection and chromatographic consistency. One aliquot was analyzed using acidic positive ion conditions, chromatographically optimized for more hydrophilic compounds. In this method, the extract was gradient eluted from a C18 column (Waters UPLC BEH C18-2.1 × 100 mm, 1.7 μm) using water and methanol, containing 0.05% perfluoropentanoic acid (PFPA) and 0.1% formic acid (FA). Another aliquot was also analyzed using acidic positive ion conditions; however, it was chromatographically optimized for more hydrophobic compounds. In this method, the extract was gradient eluted from the same afore mentioned C18 column using methanol, acetonitrile, water, 0.05% PFPA and 0.01% FA and was operated at an overall higher organic content. Another aliquot was analyzed using basic negative ion optimized conditions using a separate dedicated C18 column. The basic extracts were gradient eluted from the column using methanol and water, however with 6.5mM Ammonium Bicarbonate at pH 8. The fourth aliquot was analyzed via negative ionization following elution from a HILIC column (Waters UPLC BEH Amide 2.1 × 150 mm, 1.7 μm) using a gradient consisting of water and acetonitrile with 10mM Ammonium Formate, pH 10.8. The MS analysis alternated between MS and data-dependent MSn scans using dynamic exclusion. The scan range varied slighted between methods but covered 70–1000 m/z.

Peak analysis was conducted using a bioinformatics system which consisted of four major components: 1) the Laboratory Information Management System (LIMS - a system used to automate sample accession and preparation, instrumental analysis and reporting, and data analysis), 2) the data extraction and peak-identification software, 3) data processing tools for quality control and compound identification, and 4) a collection of information interpretation and visualization tools.

Raw data were extracted, peak-identified, and QC processed using Metabolon’s hardware and software. Compounds were identified by comparison to library entries of purified standards. Biochemical identifications were based on three criteria: 1) retention index (RI) within a narrow RI window of the proposed identification, 2) accurate mass match to the library (±10 ppm), and 3) the MS/MS forward and reverse scores between the experimental data and authentic standards. The MS/MS scores were based on a comparison of the ions present in the experimental spectrum to the ions present in the library spectrum. While there may have be similarities between these molecules based on one of these factors, the use of all three data points was utilized to distinguish and differentiate more than 3,300 registered biochemicals.

Peaks were quantified using area-under-the-curve. A data normalization step was performed to correct variation resulting from instrument inter-day tuning differences (i.e., variation between pre and postseason analyses). Specifically, each compound was corrected in run-day blocks by registering the medians to equal one (1.00) and normalizing each data point proportionately. Data were then log-transformed.

Metabolites were included in statistical analyses if 1) they significantly changed from Pre to Post following an FDR correction (q 0.05), 2) they are involved in energetics, inflammation, or exogenous consumption (see Introduction) (Abdul-Muneer et al., 2015; Lozano et al., 2015; Rodriguez-Rodriguez et al., 2014), and 3) they appeared in the random forest importance plot (as described below in Random forest analysis).

MiRNA quantification

Serum was sent out for miRNA quantification. Nine miRNAs were selected based on previous findings by Papa and colleagues: miR-20a, miR-505, miR-195, miR-30d, miR-92a, miR-93p, miR-3623p, miR-486, and miR-151-5p (Papa et al., 2019). The miRNA data described in this study were first analyzed by Papa et al. (2019) (Papa et al., 2019). Serum samples collected at Pre and Post were used to isolate and quantify levels of RNA. All samples were given a unique identifier and experimenters were blinded to the conditions (Pre vs. Post). 100 L of serum was aliquoted, and RNA was isolated using a serum/plasma isolation kit (Qiagen Inc., Venlo, Netherlands) as per the manufacturer’s protocol. RNA was eluted in 20L of DNAse/RNAse-free water and stored at −80°C until further use.

Droplet digital PCR (ddPCR; Bio-Rad Inc., Hercules, CA, USA) was then used to quantify absolute levels of the nine miRNA (Papa et al., 2019). Prior to ddPCR analysis, RNA was checked for quality using a bioanalyzer assay with a small RNA assay. After quality confirmation, 10 ng of RNA was reverse transcribed using specific miRNA TaqMan assays as per the manufacturer’s protocol (Thermo Fisher Scientific Inc., Waltham, MA, USA) in a 15 L total reaction volume; 5 L of reverse transcribed product was used to set up the real-time PCR reaction using TaqMan assays, and 20 L of the final real-time PCR reaction was mixed with 70 L of droplet oil in a droplet generator (Bio-Rad Inc., Hercules, CA, USA). After droplet formation, the PCR reaction was performed as per the recommended thermal cycling conditions. The final PCR product was analyzed using a droplet reader (Bio-Rad Inc., Hercules, CA, USA). Total positive and negative droplets were quantified, and from this, the concentration of miRNA/L of the PCR reaction was reported. All reactions were performed in duplicate.

Head acceleration event (HAE) monitoring

HAEs were monitored at all contact practice sessions (max = 53) using the BodiTrak sensor system from The Head Health Network (Slobounov et al., 2017). Sensors were individually mounted on the inner surface, between the shell and padding, of each active player’s helmet prior to contact practice beginning. To avoid game disruption, no games were monitored; however, previous studies have revealed that most HAEs accumulate from practice sessions, not games (Lee et al., 2020). Sensors were monitored throughout the season for integrity and functionality and outputs included peak translational acceleration (PTA; G-units) and impact location. HAEs were quantified as 1) cumulative hits >25G to <80G (25G) and >80G (80G) (cHAE25G and cHAE80G; Equations 1 and 2) cumulative hits exceeding 25G and 80G normalized to the total number of sessions per player (aHAE25G and aHAE80G; Equation 2). The G-unit thresholds (Th) were selected based on previous reports of impacts related to brain health and injury (McCuen et al., 2015).

Quantification and statistical analysis

All statistical analyses were performed in STATA (StataCorp, 2021) apart from the random forest analysis which was performed by Metabolon (Durham, NC, USA) and the permutation-based mediation analyses and moderation analyses which were performed in R (R Core Team, 2021). Analyses were performed with data from all 23 collegiate football athletes. The sample size was based on 1) the limited availability of participants and 2) a power analysis using expected mean metabolite values at Pre and Post which yielded a power of 86% (significance level = 0.05, sample size = 23). For analyses involving miRNA, some of the serum samples did not have a sufficient quantity for miRNA quantification, thus decreasing the number of samples from 23 to between 20 and 22 for Pre analyses and 18–20 for across-season analyses, depending on the serum availability of a given miRNA. For Pre-season analyses: n = 20 for miR-20a; n = 21 for miR-505; n = 23 for miR-363-3p, miR-30d, miR-92a, miR-486, miR-195, miR-9-3p, miR-151-5p. For across-season analyses: n = 19 for miR-20a; n = 20 for miR-505; n = 20 for miR-363-3p, miR-30d, miR-92a, miR-486, miR-195, miR-9-3p, miR-151-5p. There were no other missing data.

Random forest analysis

Random forest was used to assess the importance ranking of biochemicals (i.e., how well a given metabolite can distinguish Pre from Post) (Breiman, 2001). For a given decision tree, a random subset of the data, with identifying true class information, was selected to build the tree (“bootstrap sample” or “training set”). The remaining data, the “out-of-bag” (OOB) variables, were then passed down the tree to obtain a class prediction for each sample. This process was repeated thousands of times to produce the forest. The final classification of each sample was determined by computing the class prediction frequency for the OOB variables over the whole forest. When the full forest was grown, the class predictions were compared to the true classes, generating the “OOB error rate” as a measure of prediction accuracy. Thus, the prediction accuracy was an unbiased estimate of how well a sample class was predicted in a new dataset. To determine which metabolites made the largest contribution to the classification, a variable importance measure, Mean Decrease Accuracy (MDA), was computed. MDA was determined by randomly permuting a variable, running the observed values through the trees, and then reassessing the prediction accuracy. If a variable was not important, the procedure produced little change in the accuracy of the class prediction (permuting random noise gave random noise). By contrast, if a variable was important to the classification, the prediction accuracy dropped; this was recorded as the MDA. The top 20 metabolites were reported.

Wilcoxon signed-rank tests

Data were first tested for normality and equal variances using the Shapiro-Wilks test (Shapiro and Wilks, 1965) and Bartlett’s test (Bartlett, 1937), respectively. Because (1) data were not normally distributed (Shapiro-Wilks test; p 0.05) and/or data did not have equal variances (Bartlett’s test; p 0.05) and (2) data were paired by athlete across-season, a Wilcoxon signed-rank test (Wilcoxon, 1945) was utilized to assess across-season changes. Changes were analyzed at a significance level of 0.05. In total, there were four VR scores, nine miRNA, and 1300+ metabolites. Metabolites were further selected based on previously described criteria (see Metabolite quantification). If p-value < 0.05, metabolites were then grouped based on super pathway and FDR-corrected using the Benjamini-Hochberg (Benjamini and Hochberg, 1995) method for multiple comparisons (Table 1). Changes were considered significant if q-value 0.05.

Linear regression analysis

Linear regressions were conducted to assess significant relations between VR scores and miRNAs, miRNAs and metabolites, and VR scores and metabolites, where VR score (i.e., motor control) was always the dependent variable. When regressing miRNA and metabolite, metabolite was designated as the dependent variable; this decision was based on previous findings where miRNA levels were elevated in this cohort at Pre when compared to controls (Papa et al., 2019). After each initial regression analysis, Cook’s distance was calculated to reveal outliers which drastically influenced the regression (i.e., large shift in ) (Cook, 1977). Outliers were removed if Cook’s distance >4/N and regressions were re-run (Altman and Krzywinski, 2016). For Cook’s distance assessments, the decision of which variable was dependent variable and which one was independent variable was critical, as this would influence the resultant outliers; we thus list the sample size (N) for each mediation (see following section) and the number of outliers removed in the supplemental tables. Regressions were run for both Pre and across-season (=Post-Pre) measures and were considered significant if p-value< 0.05, and standardized beta coefficients (Std. ) and adjusted R2 values () were also reported.

Permutation-based mediation analysis

Based on results from the linear regression analyses, data were prepared for mediation analysis by first assessing three-way relationships. To be included in mediation analysis, all individual paths must have met significance (p-value< 0.05) to form a three-way relationship. Paths were as follows: A) miRNA → metabolite, B) metabolite → VR score, and C) miRNA → VR score. It should be noted that subjects had been in a static regimen for workouts, absence of impact, and diet for the three months before Pre-season data collection. We thus considered the assumption of stationarity to be plausible for Pre-season mediation testing (i.e., the causal structure among variables is not changing over time), and we followed the guidance of MacKinnon (2008) when using difference measures (Mackinnon, 2008). Cook’s outliers were removed across all regressions to achieve the same set of athletes, and regressions were re-run with common athletes. Mediation was utilized to clarify the causal relationship between the independent variable (X) and dependent variable (Y) with the inclusion of a third mediator variable (M). The mediation model proposes that instead of a direct causal relationship between X and Y, the X influences M, which then influences Y. Beta coefficients and their standard error () terms from the following linear regression equations, following the four step process of Baron and Kenny (1986) (Baron and Kenny, 1986; Kenny, 2021), were used to calculate Sobel p-values and mediation effect percentages (Teff):

Step 4: Sobel’s test was then used to test if was significantly lower than using the following equation:

Using a standard 2-tail z-score table, the Sobel p-value was determined from the Sobel z-score and the mediation effect percentage (Teff) was calculated using the following equation:

Directed mediation analysis was performed with the VR score always acting as the dependent variable (Y); this a priori hypothesis has been observed in previous studies (Chen et al., 2020). Secondly, for the primary mediation analyses, miRNA was always the independent variable (X) and metabolite was always the mediator (M):

The selection for X was, in part, based off work by Papa and colleagues where miRNA levels were observed to be elevated in collegiate football athletes, compared to controls, at both Pre and Post (Papa et al., 2019). Therefore, metabolites, which can fluctuate dynamically and are synthesized downstream of miRNA processes, were hypothesized to mediate the relationship between elevated miRNAs and VR performance. For the secondary mediation analyses, VR score was Y, miRNA was M, and metabolite was X:

The measurement matrices defined below schematize the variables used: miRNAs, VR scores, and metabolites respectively.where P is the total number of variables in matrix (i.e., the total number of miRNAs) and Q and R are the total number of variables for matrices i.e., the total number of VR scores) andi.e., (the total number of metabolites) respectively. N denotes the number of participants and matrices were measured at two time points with representing Pre and representing Post measurements.

Across-season measures for miRNAs, VR scores, and metabolites were calculated as

Permutation approaches do not require (i) random selection of samples, (ii) sample independence, (iii) Gaussian distributions, (iv) sample sizes large enough for the central limit theorem to work, or (v) similar variance in samples being compared. Permutation methods only require the assumption that the two samples being compared are interchangeable when setting up the null hypothesis (Fisher, 1935). Permutation tests re-sample observations from the original data multiple times to build empirical estimates of the null distribution for the test statistic being studied (Camargo et al., 2008). Permutation-based tests are especially well-suited for studies with small-to-moderate sample sizes as they estimate the statistical significance directly from the data being analyzed rather than making assumptions about the underlying distribution. First, the test statistic was obtained from the original dataset, then the data were randomly permuted multiple (S) times and the test statistic was computed on each permutated dataset. The statistical significance was computed by counting the number of times the statistic value obtained in the original dataset was more extreme than the statistic value obtained from the permuted datasets (K), and then dividing that value by the number of random permutations (K/S) (Camargo et al., 2008).

For this study, permutation-based mediation analysis was performed for only variables that formed three-way relationships following the steps listed below:

The non-permuted data variables were assigned as independent, dependent, and mediator variables, respectively. Linear regressions →, →, and → were run to obtain Cook’s outliers from each regression. All outliers indicated by the three regressions were removed to obtain a common set of participants.

A four-step process for mediation analysis (Baron and Kenny, 1986; Kenny, 2021) was run on non-permutated data variables with all outliers removed to obtain the reference Sobel z-score (. (For Pre analysis, were used to obtain the reference Sobel z-score ()).

Data permutation: values were randomly sampled without replacement from and to assign to and . (For Pre analysis, values from were shuffled to get .)

Across-season measures were computed from the permuted dataset . (Note that for Pre analysis, the difference was not computed.)

Similarly, and were computed. (For Pre analysis, , were computed by randomly shuffling values in .)

Mediation analysis following steps 1 through 2, as above, was performed on the permuted dataset , , and the test statistic () was obtained. (For Pre analysis, mediation was performed on the permuted dataset , , and was obtained.)

The counter variable was incremented by one if the absolute value of was greater than the absolute value of

Steps 2–6 were repeated: times. Here, .

Permutation-based Sobel p-value () was calculated as the proportion of values that were as extreme or more extreme than i.e.,

Primary mediation results were considered significant if p-values associated with terms a, b, and c (terms derived from unpermuted data) were <0.1; and < 0.05, Teff was >50%, and Teff for the secondary mediation was <30%.

Mediation was also utilized to assess if HAEs carried/mediated the relationship between Pre and Post metabolite levels.

The codes used to conduct permutation-based mediation is available as a text file in Data S1.

Moderation analysis for metabolites and HAE

Statistical moderation was used to assess whether HAEs were and important factor to address the significant change in metabolite level across-season. Moderation proposes that the strength and direction of the relationship between X and Y is controlled by the moderator variable (Mo). Moderation is characterized by the interaction term between X and Mo using a linear regression equation:

Moderation is significant if and , where indicates whether is significantly different than zero using a t-test and is the -value associated with the overall F-test for the regression equation and indicates a significant linear relationship. indicates a significant interaction between X and Mo and indicates that the interaction term of X∗Mo directionally models Y. Here, metabolite level at Pre was X, HAE was the Mo, and metabolite level at Post was Y. A Bonferroni correction was applied by dividing 0.05 () by five (the number of HAE measures tested); this correction was applied to to protect for the overall effect.

Post-hoc miRNA-gene network analysis

To enhance the interpretability of the results, post-hoc miRNA-gene network analyses were conducted, using miRNet (Le Chang et al., 2020), to investigate relationships between the miRNAs of interest and genes implicated in mitochondrial processes (fatty acid oxidation, the TCA cycle, the electron transport chain (ETC), oxidative phosphorylation, general aerobic respiration) and neuroinflammation (Alston et al., 2017; Bult et al., 2019; Houten et al., 2016; Stelzer et al., 2016; Wanders et al., 2010). Genes of interest were selected using the Mouse Genome Database (Bult et al., 2019). Analyses were conducted selecting H. Sapiens as the species of interest and selecting either brain tissue as the tissue of interest, or by not specifying a tissue of interest. Network analyses were conducted separately for genes involved in mitochondrial processes and those involved in neuroinflammation.

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Biological samples
Human venous blood	Peripheral human blood obtained by Slobounov and maintained by Papa	N/A
Critical commercial assays
TurboVap	Zymark	N/A
ACUITY UPLC	Waters	N/A
Q-Exactive high-resolution/accurate mass spectrometer	Thermo Scientific	N/A
Electrospray ionization (HESI-II) source	Thermo Scientific	N/A
Orbitrap mass analyzer	Thermo Scientific	N/A
C18 column; UPLC BEH C18-2.1 × 100mm, 1.7 μm	Waters	N/A
HILIC column, UPLX BEH Amide 2.1 × 100m, 1.7 μm	Waters	N/A
Serum/plasma isolation kit	Qiagen Inc.	N/A
Droplet digital PCR (ddPCR)	Bio-Rad Inc.	N/A
Small RNA assay
miRNA TaqMan assay	Thermo Scientific	N/A
Droplet generator	Bio-Rad Inc.	N/A
Droplet reader	Bio-Rad Inc.	N/A
Bioanalyzer assay
BodiTrak sensor system	The Head Health Network	N/A
MicroLab STAR system	Hamilton Company	N/A
3D system, head-mounted accelerometer	Head Rehab LLC	N/A
GenoGrinder2000	Glen Mills	N/A
Deposited data
Vike_iScience_multiomics_dataset	Mendeley Data	Mendeley Data: https://doi.org/10.17632/5tygm8sn24.1
Software and algorithms
Permutation-based mediation analysis code	Available as Supplemental Material	N/A
STATA	StataCorp	https://download.stata.com/download/
R	R Core Team	https://www.R-project.org/
VR software	HeadRehab LLC	N/A